旋转圆盘上分数阶粘弹性流体非稳态流动

赵金虎，王志刚

（阜阳师范大学 数学与统计学院，安徽 阜阳 236037）

近年来，分数阶的本构关系模型在非牛顿流体的表征中应用非常广泛，不仅能利用较少的参数准确模拟粘弹性行为，也为描述不同物质的记忆和继承性提供了强有力的工具[1-4]。一般而言，分数阶的本构关系模型都是从经典的流体力学模型（例如Maxwell 模型等）演化而来，即通过将应力与应变之间的整数阶导数替换成为分数阶导数。Bai 等[5]研究了加速平板上分数阶Maxwell磁流体的速度滑移流动，获得了幂律函数边界条件下的数值解。Liu 等解决了分数阶粘弹性流体拉伸板上边界层流动的变厚度问题[6,7]以及分布阶时间导数模型的传热问题[8]。Chen 等[9]在拉伸板上粘弹性磁流体边界层流动中考虑双参数分数阶Maxwell 模型，求出了特解的数值精度。Jiang等[10]研究了多孔介质中具有霍尔效应分数阶流体的非稳态磁流动与传热传质，利用差分方法和快速算法得到了数值解。Shen 等[11,12]研究了含有Cattaneo 热流量模型的分数阶粘弹性基纳米流体，通过有限差分方法得到了速度场和温度场的数值解。Yang 等[13]使用分布阶时间分数阶导数模型研究了非稳态自然对流传热，利用中心求积方法获得了离散数值解。Long 等[14]引用分布阶时间导数模型讨论了马兰戈尼边界层传热问题，采用有限差分方法讨论了方程组的可解性。更多分数阶粘弹性流体的详细研究可参见文献[15-20]。

本文在圆盘旋转流动的本构关系中引入分数阶Maxwell 模型，旨在获得粘弹性流体三维非稳态流动的数值解。控制方程组不仅含有多项时间分数阶导数，而且在对流项中包括混合时间-空间算子。该方程组关于速度分量是强耦合的，在原点位置具有奇异性。本文为获得较高准确度的数值解，推出了隐式有限差分格式，分析了分数阶导数参数对速度场的影响。

1 数学模型

考虑半径为R 的旋转圆盘上粘弹性流体的三维非稳态对称流动，绕z轴的角速度为Ω。圆柱坐标系记作(r,φ,z)，其中z轴通过圆盘中心垂直向上，如图1 所示。粘弹性流体充满了圆盘上侧的半无限区域。初始时刻，圆盘和流体处于静止状态。

图1 坐标示意图

假设：（1）圆盘无渗透且无滑移；（2）∂p/∂r=0且∂p/∂z=0。控制方程组为：

其中u，v和w分别表示径向、角向和轴向的速度分量，τzr和τzφ分别表示剪切应力分量，ρ表示密度且为常数。

分数阶Maxwell 模型下关于τzr和τzφ的本构关系式按照如下形式给出[21]：

其中λ是松弛时间，μ是动力粘度，α是分数阶导数参数，表征了粘弹性模型中高分子材料依赖于频率的复模量[22]，是Caputo 分数阶导数算子，定义为[23，24]：

其中Γ(·)是Gamma 函数。结合方程（1）、（2）和（3），可以得到分数阶的动量方程组：

其中vf=μ/ρ表示运动粘度。

初始条件和边界条件为：

为了简化方程组，引入下列无量纲变量：

其中Re 表示雷诺数。变量z 和轴向速度w使用Re 进行了无量纲化，主要是为了保证各坐标量和速度分量在数值计算中保持在同阶量级，这样消除了无量纲后方程组中的Re，以便于集中分析分数阶导数参数的作用。

最后，获得了无量纲的控制方程组（为了简便，这里往后省略了无量纲记号“*”）：

无量纲的初始条件和边界条件分别为：

2 数值格式

本部分给出边界层控制方程组的有限差分格式，其中径向上的网格点采用半移位形式。定义ri=(i-1/2)Δr，i=1,2,…,M；zj=jΔz，j=0,1,2,…,N；tk=kΔt，k=0,1,2,…,K；其中Δr=2/(2M-1)和Δz=Zmax/N分别为空间步长，Δt是时间步长。在网格点(ri,zj,tk)的径向速度数值解记为，角向速度记为，轴向速度记为。

方程（9）中角向速度的时间导数采用向后差分：

速度对流项经过线性化处理后，空间导数采用隐式差分格式：

方程右端的扩散项采用CN 格式：

对于方程（8）中径向速度的各项导数采用相似的处理过程。方程组中br和bφ分别代表了离心力和自转偏向力，这两项是由于轴向转动速度引起的，离散格式为[25]：

对于分数阶导数，本文引用L1 近似格式进行离散[26]：

其中 0＜α＜1，αs=(s+1)1-α-s1-a，s=0,1,2,…,K.把整数阶导数的离散格式（9-11）、（13）和（14）代入到方程（15）中，可以得到各项分数阶导数的离散格式，例如：

其中初始速度为0。

连续性方程（7）的离散格式为：

最后，得到了迭代差分方程组。方程（9）中，每个内节点处的径向速度的数值解记为，该离散方程的系数矩阵是(N-2)×(N-2)的三对角稀疏矩阵。在所有网格的迭代截断误差控制在10-4以内，计算区域的两个边界长分别为Rmax=1 和Zmax=12，其中Zmax对应于z→∞且位于速度边界层之外。为了平衡数值解的精度和计算耗时，网格数选定为M=100,N=240，时间步长控制为Δt=0.1。图2 反映了一组参数下的径向速度三维分布，从图中可以看出，数值结果在开始时刻、原点附近出现了轻微震荡，但是随着时间和空间的迭代而趋于稳定和收敛。这个结果证明了原点的奇异性被消除，并且选择的计算网格是合适的。特别指出，本文提出的有限差分格式是条件稳定的，当松弛时间λ的值超过一定范围时，数值解不稳定。

图2 三维径向速度分布图（λ=0.1，α=0.6）

3 结果讨论

图3 至图5 分别给出了不同分数阶导数对r=1 处旋转充分发展的速度分量分布影响。图3表示不同α下的径向速度分布，每一条速度曲线都有一个最大值，在原点附近快速增加且当远离原点时缓慢减少到零。不同速度曲线的最大值几乎相等，但是极值点随着α的增加沿着z轴向上移动。值得注意的是，速度曲线两两相交，意味着随着流动的发展，分数阶导数参数对速度的影响会截然相反。开始流动时，速度随着α的增加而下降，但是经过一段时期的发展，速度随着α的增加而增加。这是因为分数阶粘弹性流体表现出延迟反应，所以对于较大的α，速度发展得更慢。α=0 代表了牛顿流体，通过与之对比，当α增大时，速度边界层变厚，表明了分数阶粘弹性流体具有剪切增稠的特性。

图3 不同α 下的径向速度分布

图4 不同α 下的角向速度分布

图5 不同α 下的轴向速度分布

图4 给出了不同α下的角向速度分布。每一条角向速度曲线是单调下降的，随着α的增加，角向速度增加而边界层的厚度也变得更大。角向速度对小α的值更加敏感，这表现为分布曲线之间的间距随着α的增加而减小。结果表明，流体的粘弹性特性增加时会导致流动阻力增加，所以牛顿流体（α=0）的动量传递更迅速。图5 给出了不同α下的轴向速度分布，负号的含义表示轴向速度的方向与z轴的正向相反。轴向速度在圆盘附近几乎无变化，但是远离圆盘时随着α的增加而增加。

4 小结

本文讨论了旋转圆盘上粘弹性流体的三维非稳态流动问题，在本构关系中引入了分数阶Maxwell 模型，建立了具有多项时间分数阶导数的非线性控制方程组，采用半移位的网格系统解决了原点的奇异性，提出了条件稳定的隐式差分格式。本文建立的分数阶控制方程组具有强非线性和耦合性，不易得到差分方程的收敛误差精度和稳定性的证明，这些问题将在下一步的研究中解决。