1.   引　言

黏弹性流体作为非牛顿流体力学领域里的一类重要研究对象, 广泛存在于自然界和工业生产过程中, 如石油输送和聚合物加工成型等. 为了提高产品质量和生产效率, 对黏弹性流体复杂流变特性的相关研究已成为工业生产过程中一个重要的环节. 目前, 基于网格的数值方法, 如有限差分法 (finite difference method, FDM)、有限体积法和有限元等, 已被用于该类问题的数值模拟. Viezel等^[1]使用FDM模拟了Oldroyd-B液滴碰撞问题. Li等^[2]通过流体体积法模拟了黏弹性液滴在低Weber数下的流动. 尽管基于网格的数值方法在模拟黏弹性流体方面取得了成功, 但其自身依赖于网格生成的质量. 此外, 对于含自由面的黏弹性流动模拟, 需借助额外的界面追踪技术, 实施过程相对繁琐.

近年来, 无网格粒子法^[3]引起了广大学者的关注以及重视. 光滑粒子动力学 (smoothed particle hydrodynamics, SPH)方法是一种典型的无网格粒子法, 1977年由Gingold, Monaghan^[4]和Lucy^[5]提出. 与基于网格的数值方法相比, SPH方法具有无网格特性、自适应特性、Lagrangian特性和易于编程等优势, 因此该方法已被成功用于黏性流^[6]、动画^[7]、传热^[8]、多相流^[9]、爆炸^[10]和流固耦合^[11]等领域中. 然而, 传统SPH方法也存在一些缺点, 如计算精度低、稳定性差和边界不易处理等. 针对这些问题, 近年来很多学者已提出了多种不同的改进方案. 例如, 为了提高计算精度, Liu等^[12]用积分再生核思想发展了再生核粒子方法(reproducing kernel particle method, RKPM); Liu等^[13]推导了具有较高计算精度的有限粒子法(finite particle method, FPM); Fang等^[14]提出了模拟黏性不可压缩流动的高精度SPH方法. 为了改善SPH方法的数值稳定性, Yang等^[15]提出了一种新的核函数, 克服了黏性液滴形成过程中的张力不稳定; Antuono等^[16]在密度方程中引入密度耗散项, 以解决流动过程中出现的压力不稳定; Lyu等^[17]发展了粒子迁移技术, 保证流动过程中流体粒子的规则分布. 为了施加边界条件, Monaghan等^[18]提出了只施加一层粒子在边界的边界排斥力法; Morris等^[19]发展了布置多层粒子在边界的虚拟粒子法; Liu等^[20]结合排斥力法和虚拟粒子法的优势, 建立了两者相耦合的动力学边界处理技术. 虽然这些改进方法已被成功提出并得到一定的应用, 但它们在黏弹性流体力学领域中的应用并不多见.

对于黏弹性流体的SPH模拟, Fang等^[21]将SPH方法推广到含自由面的Oldroyd-B黏弹性液滴撞击固壁问题的模拟中. Hashemi等^[22]研究了两个相互作用的固体颗粒在Oldroyd-B剪切流场下的迁移运动, 分析了Deborah数对固体颗粒运动轨迹的影响. Xu等^[23]提出了改进SPH方法模拟基于Oldroyd-B模型的黏弹性液滴撞击固壁和挤出胀大问题. Ozgen等^[24]利用SPH方法捕捉了非连续剪切增稠流体特有的流动现象. Vahabi和Kamkari^[25]基于SPH方法数值模拟了气泡在Giesekus黏弹性溶液中的上升和变形. King和Lind^[26]将对数构象公式耦合到SPH方法模拟了高Weissenberg数下的黏弹性圆柱绕流问题.

eXtended Pom-Pom (XPP)模型由Verbeeten等^[27]在Pom-Pom模型的基础上基于支链聚合物的分子理论推导得到. 与传统的简单模型(例如Oldroyd-B, 上随体Maxwell模型)相比, XPP模型能够更合理地描述聚合物溶液的剪切和拉伸行为, 且在一定程度上克服了应力奇点问题. 此外大量实验^[27]表明, XPP模型能够合理地表征黏弹性流体的剪切稀化. 根据文献调研结果, 目前运用SPH方法模拟XPP黏弹性流体的文献报道并不多见, 因此本文重点围绕基于XPP模型的黏弹性流动问题展开研究.

本文首先对传统SPH方法进行改进. 特别地, 采用一种不含核导数计算的核梯度修正离散格式以提高计算精度. 为了防止粒子穿透固壁, 提出一种增强型的边界处理技术. 为了消除张力不稳定性, 将人工应力耦合到动量守恒方程中. 随后, 运用改进SPH方法数值模拟了基于XPP模型的黏弹性Poiseuille流和黏弹性液滴撞击固壁问题, 通过与解析解或其他方法解的比较和对数值收敛性的评价验证了改进SPH方法的有效性和优势, 并在此基础上, 深入分析了Reyonlds数、Weissenberg数、溶剂黏度比、各向异性参数、松弛时间比和分子链臂数等流变参数对流动过程的影响.

2.   控制方程

在Lagranigian坐标系下, 二维等温、黏弹性流体的控制方程可写为

$$\frac{{{\text{d}}\rho }}{{{\text{d}}t}} = - \rho \frac{{\partial {v^\beta }}}{{\partial {x^\beta }}} ,$$

(1)

$$\frac{{{\text{d}}{v^\alpha }}}{{{\text{d}}t}} = \frac{1}{\rho }\frac{{\partial {\sigma ^{\alpha \beta }}}}{{\partial {x^\beta }}} + {F^\alpha } ,$$

(2)

其中, $\rho$是流体密度, t 是时间, ${x^\beta }$是空间坐标, ${v^\beta }$为速度, ${{\sigma} ^{\alpha \beta }}$为总应力张量, ${F^\alpha }$是作用在流体上的外力, $\dfrac{{\text{d}}}{{{\text{d}}t}}$表示物质导数, 即$\dfrac{{\text{d}}}{{{\text{d}}t}} = \dfrac{\partial }{{\partial t}} + {v^\beta }\dfrac{\partial }{{\partial {x^\beta }}}$. (2)式中的总应力张量${{\sigma }^{\alpha \beta }}$通常可分解为压力$p$、黏性应力${\tau} _{\text{s}}^{\alpha \beta }$和弹性应力${\tau} _{\text{p}}^{\alpha \beta }$三者之和, 即

$${{\sigma }^{\alpha \beta }} = - p{{\delta} ^{\alpha \beta }} + {\tau} _{\text{s}}^{\alpha \beta } + {\tau} _{\text{p}}^{\alpha \beta } ,$$

(3)

式中, ${{\delta} ^{\alpha \beta }}$为Kronecker函数. 黏性应力$\tau _{\text{s}}^{\alpha \beta }$可按牛顿流体模型进行计算:

$${\tau} _{\text{s}}^{\alpha \beta } = 2{\eta _{\text{s}}}{{d}^{\alpha \beta }} ,$$

(4)

其中, ${\eta _{\text{s}}}$是牛顿溶剂黏度; ${{d}^{\alpha \beta }}$是形变率张量,

$${{d}^{\alpha \beta }} = \frac{1}{2}\left(\frac{{\partial {v^\alpha }}}{{\partial {x^\beta }}} + \frac{{\partial {v^\beta }}}{{\partial {x^\alpha }}}\right) .$$

(5)

2.1   XPP本构方程

为了封闭控制方程, 需附加一个与弹性应力$\tau _{\text{p}}^{\alpha \beta }$相关的本构方程. 本文考虑XPP模型, 其本构方程为

$$\begin{split} & f\left( {\lambda ,\;{\boldsymbol{\tau}} } \right)\tau _{\text{p}}^{\alpha \beta } + {\lambda _1}\mathop {\tau _{\text{p}}^{\alpha \beta }}\limits^\nabla + {G_0}\left[ {f\left( {\lambda ,\;{\boldsymbol{\tau}} } \right) - 1} \right]{\delta ^{\alpha \beta }} \\ & + \frac{\alpha }{{{G_0}}}\left( {\tau _{\text{p}}^{\alpha \gamma } \cdot \tau _{\text{p}}^{\gamma \beta }} \right) = 2{\lambda _1}{G_0}{d^{\alpha \beta }} , \end{split}$$

(6)

式中,

$$\begin{split} f\left( {\lambda ,\;{\boldsymbol{\tau}} } \right) =& 2\frac{{{\lambda _1}}}{{{\lambda _2}}}{{\text{e}}^{{Q_0}\left( {\lambda - 1} \right)}}\left( {1 - \frac{1}{\lambda }} \right) \\ & +\frac{1}{{{\lambda ^2}}}\left[ {1 - \frac{\alpha }{{3G_0^2}}{\text{tr}}\left( {{\boldsymbol{\tau}} \cdot {\boldsymbol{\tau}} } \right)} \right] ; \end{split}$$

(7)

τ为弹性应力的矩阵形式, $\tau _{\text{p}}^{\mathop {\alpha \beta }\limits^\nabla }$是$\tau _{\text{p}}^{\alpha \beta }$的上随体导数, 即

$$\tau _{\text{p}}^{\mathop {\alpha \beta }\limits^\nabla } = \frac{{{\text{d}}\tau _{\text{p}}^{\alpha \beta }}}{{{\text{d}}t}} - \frac{{\partial {v^\alpha }}}{{\partial {x^\gamma }}}\tau _{\text{p}}^{\gamma \beta } - \frac{{\partial {v^\beta }}}{{\partial {x^\gamma }}}\tau _{\text{p}}^{\alpha \gamma };$$

(8)

${\lambda _1}$和${\lambda _2}$分别表示聚合物分子链的取向松弛和拉伸松弛时间; $\gamma = {\lambda _2}/{\lambda _1}$表示拉伸松弛与取向松弛时间之比; ${G_0}$表示线性松弛模量; ${\eta _{\text{p}}} = {G_0}{\lambda _1}$表示聚合物黏度, $\eta = {\eta _{\text{s}}} + {\eta _{\text{p}}}$表示流体总黏度; $\alpha$为控制各向异性拖曳力的流变参数; ${Q_0}$反比于分子链臂数Q, 即${Q_0}{\text{ = }}2/Q$; 分子链拉伸量$\lambda$的表达式为

$$\lambda = \sqrt {1 + \frac{1}{{3{G_0}}}{\text{tr}}\left( {\boldsymbol{\tau}} \right)} .$$

(9)

由(6)式—(9)式, 可得

$$\begin{split} \frac{\text{d}{\tau }_{\text{p}}^{\alpha \beta }}{\text{d}t}= &\frac{\partial {v}^{\alpha }}{\partial {x}^{\gamma }}{\tau }_{\text{p}}^{\gamma \beta }+\frac{\partial {v}^{\beta }}{\partial {x}^{\gamma }}{\tau }_{\text{p}}^{\alpha \gamma }-\frac{1}{{\lambda }_{1}}f\left(\lambda ,{\boldsymbol{\tau}} \right){{{\tau}} }_{\text{p}}^{\alpha \beta }\\ & -\frac{1}{{\lambda }_{1}}{G}_{0}\left[f\left(\lambda ,{\boldsymbol{\tau}} \right)-1\right]{\delta }^{\alpha \beta } \\ & -\frac{\alpha }{{\lambda }_{1}{G}_{0}}\left({{\tau} }_{\text{p}}^{\alpha \gamma }\cdot {{\tau} }_{\text{p}}^{\gamma \beta }\right) +2{G}_{0}{d}^{\alpha \beta }. \\[-10pt] \end{split}$$

(10)

显然, 当$\alpha = 0$和$f\left( {\lambda , \;{\boldsymbol{\tau }}} \right) = 1$时, XPP本构方程将退化为Oldroyd-B本构方程.

2.2   状态方程

在SPH框架内, 有两种常用的求解控制方程的方法: 一种是不可压SPH方法^[28], 其中流体压力通过Possion方程隐式计算; 另一种是弱可压SPH方法^[23], 其中压力根据状态方程显式计算. 本文采用后者, 使用的状态方程为

$$p(\rho ) = {c^2}(\rho - {\rho _0}) ,$$

(11)

式中, $c$是声速, ${\rho _0}$是初始流体密度. 在弱可压SPH方法中, 马赫数应小于0.1. 基于此, 声速$c$应取为不低于流动特征速度值的10倍, 以保证弱可压流体的流动行为足够接近不可压流体的行为. 声速$c$越大, 模拟所用时间步长$\Delta t$越小.

3.   SPH方法

3.1   SPH离散

对流体控制方程可进行多种不同的SPH离散. 对质量、动量守恒方程和XPP本构方程采用的SPH离散为

$$\frac{{{\text{d}}{\rho _i}}}{{{\text{d}}t}} = \sum\limits_j {{m_j}\left( {v_i^\beta - v_j^\beta } \right)\frac{{\partial {W_{ij}}}}{{\partial x_i^\beta }}} ,$$

(12)

$$\frac{{{\text{d}}v_i^\alpha }}{{{\text{d}}t}} = \sum\limits_j {{m_j}\left( {\frac{{\sigma _i^{\alpha \beta }}}{{\rho _i^2}} + \frac{{\sigma _j^{\alpha \beta }}}{{\rho _j^2}}} \right)\frac{{\partial {W_{ij}}}}{{\partial x_i^\beta }} + F_i^\alpha } ,$$

(13)

$$\begin{split} {\left(\frac{\text{d}{\tau }_{\text{p}}^{\alpha \beta }}{\text{d}t}\right)}_{i}=& {k}_{i}^{\alpha \gamma }{\tau }_{\text{p},i}^{\gamma \beta }+{k}_{i}^{\beta \gamma }{\tau }_{\text{p},i}^{\gamma \alpha }-\frac{1}{{\lambda }_{1}}f\left(\lambda ,{\boldsymbol{\tau}} \right){\tau }_{\text{p},i}^{\alpha \beta }\\ & -\frac{1}{{\lambda }_{1}}{G}_{0}\left[f\left(\lambda ,{\boldsymbol{\tau}} \right)-1\right]{\delta }^{\alpha \beta }\\ & -\frac{\alpha }{{\lambda }_{1}{G}_{0}}\left({\tau }_{\text{p},i}^{\alpha \gamma }\cdot {\tau }_{\text{p},i}^{\gamma \beta }\right) \\ & +{G}_{0}\left({k}_{i}^{\alpha \beta }+{k}_{i}^{\beta \alpha }\right), \\[-15pt] \end{split}$$

(14)

其中,

$$\sigma _i^{\alpha \beta } = - p{\delta ^{\alpha \beta }} + {\eta _{\text{s}}}\left( {k_i^{\alpha \beta } + k_i^{\beta \alpha }} \right) + \tau _{{\text{p}},i}^{\alpha \beta } ,$$

(15)

$$k_i^{\alpha \beta } = {\left( {\frac{{\partial {v^\alpha }}}{{\partial {x^\beta }}}} \right)_i} = \sum\limits_j {\frac{{{m_j}}}{{{\rho _j}}}\left( {v_j^\alpha - v_i^\alpha } \right)\frac{{\partial {W_{ij}}}}{{\partial x_i^\beta }}} ,$$

(16)

i和j表示粒子编号, m表示粒子质量, ${W_{ij}} = W(\left| {{x_i} - {x_j}} \right|, h)$是核函数, $h$是核函数影响区域的光滑长度.

传统SPH方法的计算精度低, 限制了其向更复杂流动问题的应用. 为了提高计算精度, 许多学者已提出多种修正的SPH算法, 如RKPM^[12]和FPM^[13]等. 对SPH离散(12)式—(16)式中的核函数梯度进行修正, 即利用如下不含核导数计算的修正矩阵^[29]:

$${\boldsymbol{M}} = \left( {\begin{array}{*{20}{c}} {\displaystyle\sum\limits_j {\frac{{{m_j}}}{{{\rho _j}}}{x_{ji}}{x_{ji}}{W_{ij}}} }&{\displaystyle\sum\limits_j {\frac{{{m_j}}}{{{\rho _j}}}{y_{ji}}{x_{ji}}{W_{ij}}} } \\ {\displaystyle\sum\limits_j {\frac{{{m_j}}}{{{\rho _j}}}{x_{ji}}{y_{ji}}{W_{ij}}} }&{\displaystyle\sum\limits_j {\frac{{{m_j}}}{{{\rho _j}}}{y_{ji}}{y_{ji}}{W_{ij}}} } \end{array}} \right) ,$$

(17)

将原核函数梯度修正为

$$\left( {\begin{array}{*{20}{c}} {\dfrac{{\partial {{\tilde W}_{ij}}}}{{\partial {x_i}}}} \\ {\dfrac{{\partial {{\tilde W}_{ij}}}}{{\partial {y_i}}}} \end{array}} \right) = {{\boldsymbol{M}}^{ - 1}}\left( {\begin{array}{*{20}{c}} {\dfrac{{\partial {W_{ij}}}}{{\partial {x_i}}}} \\ {\dfrac{{\partial {W_{ij}}}}{{\partial {y_i}}}} \end{array}} \right) .$$

(18)

利用新的核梯度(18)式进行计算, 可以提高SPH方法的计算精度. 另外, 值得注意的是, 本文核梯度修正矩阵不含核函数一阶导数的计算, 因此对核函数的适用性更高, 进行相应模拟的稳定性也更好. 关于改进SPH离散的更多内容, 可参阅文献[29].

3.2   边界处理

边界处理是数值模拟的一个重要部分, 其处理是否得当关乎SPH模拟的稳定性和精度. 目前, 较为广泛使用的边界方法有排斥力法^[18]和虚拟粒子法^[19]等. 排斥力法只施加一层粒子在边界, 程序设计简单, 但守恒性较差. 虚拟粒子法布置多层粒子在边界, 解决了邻近边界处粒子被边界截断的问题, 但虚拟粒子的速度、压力等各物理量需通过构造伪粒子进行Shepard插值得到. 对于复杂边界, 伪粒子的构造较为繁琐.

本文提出一种由边界粒子和虚拟粒子两类粒子组成的增强型边界处理技术. 首先, 在边界上以粒子初始间距$\Delta x$布置一层边界粒子. 与排斥力法^[18]不同的是, 边界粒子不再对邻近边界处的流体粒子施加排斥力. 相反, 边界粒子参与控制方程中速度、压力和弹性应力的计算, 以有效防止粒子穿透边界. 计算过程中, 边界粒子的密度和位置保持不变. 速度设置为0, 以施加无滑移边界条件. 压力和弹性应力通过对相邻流体粒子进行Shepard插值得到:

$${A_i} = \frac{{\displaystyle\sum\nolimits_j {\frac{{{m_j}}}{{{\rho _j}}}{A_j}{W_{ij}}} }}{{\displaystyle\sum\nolimits_j {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}} }} ,$$

(19)

式中, A表示压力或弹性应力, $i$为边界粒子, $j$为与粒子$i$相邻的流体粒子.

其次, 为解决邻近边界处粒子被边界截断的问题, 在边界外以粒子初始间距$\Delta x$规则布置三层虚拟粒子. 与边界粒子类似, 虚拟粒子的密度和位置在计算过程中也保持不变. 但与虚拟粒子法^[19]不同的是, 虚拟粒子的各物理量不再通过构造伪粒子进行Shepard插值得到. 本文每个虚拟粒子均与边界法向上的一个边界粒子相对应. 为满足压力和弹性应力的纽曼边界条件, 虚拟粒子的压力和弹性应力设置与边界法向上对应的边界粒子值相同. 另外, 所有虚拟粒子的速度设置为0, 以满足无滑移边界条件.

相较于传统的排斥力法^[18]和虚拟粒子法^[19], 本文提出的增强型边界处理技术不仅程序设计简单, 而且避免了构造伪粒子进行Shepard插值的繁琐操作. 数值试验表明, 本文提出的增强型边界处理技术能够获得精确稳定的计算结果; 即便对于冲击力较大的黏弹性液滴撞击固壁问题, 也能够有效地防止流体粒子穿透边界.

3.3   人工应力

当用SPH方法模拟固体弹性变形问题时, 会产生张力不稳定性. 在进行含自由面的黏弹性流动问题模拟时, 由于弹性力的作用, 也会产生张力不稳定性, 即粒子在运动过程中形成团块并最终导致模拟中断. 对于该问题, 目前已有一些解决方法, 如采用新核函数法^[15]、粒子位移修正法^[30]等. 本文将Monaghan^[31]和Gray等^[32]提出的人工应力耦合到动量守恒方程中, 以消除张力不稳定性. 人工应力的基本思想是在一对相邻的粒子之间引入一个小的短程排斥力, 以防止两个粒子靠的太近而形成团块, 其表达式为

$$S_{ij}^{\alpha \beta } = f_{ij}^n(R_i^{\alpha \beta } + R_j^{\alpha \beta }) ,$$

(20)

其中, $n = W(0, h)/W(\Delta x, h)$为指数型因子, ${f_{ij}} = {W_{ij}}/W(\Delta x, h)$为人工应力系数, ${R^{\alpha \beta }}$为人工应力分量, $\Delta x$为粒子初始间距.

二维问题中, 人工应力张量${R^{\alpha \beta }}$的分量可通过以下方式转换得到:

$$\left\{\begin{aligned} & {R}^{xx}={{R}^{\prime }}^{xx}{\rm{cos}}^{2}\theta +{{R}^{\prime }}^{yy}{\rm{sin}}^{2}\theta \text{, }\\ & {R}^{yy}={{R}^{\prime }}^{xx}{\rm{sin}}^{2}\theta +{{R}^{\prime }}^{yy}{\rm{cos}}^{2}\theta \text{, }\\ & {R}^{xy}={\rm{sin}}\theta {\rm{cos}}\theta ({{R}^{\prime }}^{xx}-{{R}^{\prime }}^{yy}), \end{aligned} \right.$$

(21)

其中, ${R'^{xx}}$和${R'^{yy}}$分别是沿$x'$轴和$y'$轴人工应力张量的对角分量. 坐标系$(x', y')$相对于$(x, y)$旋转$\theta$角度, 即

$$\theta = \frac{1}{2}{\tan ^{ - 1}}\left(\frac{{2{\sigma ^{xy}}}}{{{\sigma ^{xx}} - {\sigma ^{yy}}}}\right) .$$

(22)

${R'^{xx}}$的表达式为

$${{R}^{\prime }}^{xx}=\left\{\begin{array}{l}-\varepsilon \dfrac{{{\sigma }^{\prime }}^{xx}}{{\rho }^{2}},\text{ }\text{ }{{\sigma }^{\prime }}^{xx} \gt 0\text{, }\\ 0,\qquad\quad\;\; \text{otherwise}\text{, }\end{array}\right.$$

(23)

其中$\varepsilon$是一个介于0和1之间的人工应力参数, ${R'^{yy}}$也是类似的形式. 对角分量${\sigma '^{xx}}$和${\sigma '^{yy}}$可通过以下关系式计算:

$${\sigma '^{xx}} = {\cos ^2}\theta {\sigma ^{xx}} + 2\cos \theta \sin \theta {\sigma ^{xy}} + {\sin ^2}\theta {\sigma ^{yy}} ,$$

(24)

$${\sigma '^{yy}} = {\sin ^2}\theta {\sigma ^{xx}} - 2\cos \theta \sin \theta {\sigma ^{xy}} + {\cos ^2}\theta {\sigma ^{yy}} .$$

(25)

于是, 通过耦合人工应力, 动量方程(13)的SPH离散形式修正为

$$\frac{{{\text{d}}v_i^\alpha }}{{{\text{d}}t}} = \sum\limits_j {{m_j}\left(\frac{{\sigma _i^{\alpha \beta }}}{{\rho _i^2}} + \frac{{\sigma _j^{\alpha \beta }}}{{\rho _j^2}} + S_{ij}^{\alpha \beta }\right)\frac{{\partial {W_{ij}}}}{{\partial x_i^\beta }} + F_i^\alpha } .$$

(26)

3.4   时间积分

为了求解控制方程, 采用蛙跳公式进行时间积分. 令${{\boldsymbol{X}}_i}$表示未知变量$\left( {{\rho _i}, {{\boldsymbol{u}}_i}, {{\boldsymbol{\tau}} _i}} \right)$的向量, ${{\boldsymbol{B}}_i}$表示所对应方程的右端向量, 那么蛙跳公式可概括如下:

在第一个时间步长${t_0}$结束后, ${{\boldsymbol{X}}_i}$向前推进半个时间步, ${{\boldsymbol{r}}_i}$向前推进一个时间步:

$$\left\{\begin{aligned} & t={t}_{0}+\Delta t\text{ }\text{, }\text{ }\text{ }\text{ }\\ & {{\boldsymbol{X}}}_{i}\left({t}_{0}+\Delta t/2\right)={{\boldsymbol{X}}}_{i}\left({t}_{0}\right)+{{\boldsymbol{B}}}_{i}\left({t}_{0}\right)\frac{\Delta t}{2}\text{, }\text{ }\\ & {{\boldsymbol{r}}}_{i}\left({t}_{0}+\Delta t\right)={{\boldsymbol{r}}}_{i}\left({t}_{0}\right)+{{\boldsymbol{u}}}_{i}\left({t}_{0}+\Delta t/2\right)\Delta t. \end{aligned} \right.$$

(27)

为了实现每个时间步的一致性, 在每个时间步的开端, 各粒子的速度向前推进半个时间步, 从而保证和位移的一致性:

$${{\boldsymbol{X}}_i}\left( t \right) = {{\boldsymbol{X}}_i}\left( {t - \Delta t/2} \right) + {{\boldsymbol{B}}_i}\left( {t - \Delta t} \right)\frac{{\Delta t}}{2} .$$

(28)

在随后时间步的末端, 粒子的密度、速度和位移按照标准的蛙跳公式进行推进:

$$\left\{\begin{aligned} & t=t+\Delta t\text{ }\text{, }\text{ }\text{ }\text{ }\\ & {{\boldsymbol{X}}}_{i}\left(t+\Delta t/2\right)={{\boldsymbol{X}}}_{i}\left(t-\Delta t/2\right)+{{\boldsymbol{B}}}_{i}\left(t\right)\Delta t\text{, }\\ & {{\boldsymbol{r}}}_{i}\left(t+\Delta t\right)={{\boldsymbol{r}}}_{i}\left(t\right)+{{\boldsymbol{u}}}_{i}\left(t+\Delta t/2\right)\Delta t. \text{ }\end{aligned}\right.$$

(29)

此外, 为了保持数值稳定性, 时间步长$\Delta t$必须满足以下条件:

$$\Delta t \leqslant \min \left( {\frac{h}{c},{\text{ }}\mathop {{\min}}\limits_i \sqrt {\frac{h}{{{F_i}}}} ,{\text{ }}0.125\frac{{{h^2}}}{{{\upsilon _0}}}} \right) ,$$

(30)

其中, ${F_i}$表示每单位质量的力, ${\upsilon _0} = \eta /\rho$表示流体的运动黏度.

4.   数值算例

4.1   黏弹性Poiseuille流

目前, 尚未见到基于XPP模型对黏弹性Poiseuille流进行SPH模拟研究的相关文献报道, 因此首先选取该问题作为数值模拟的第一个算例, 并分析XPP模型中各流变参数对黏弹性流动过程的影响. 在模拟过程中, 为了更好地描述黏弹性流动, 引入无量纲Reynolds数

$$Re = {{\rho VL}}/{\eta }\text{, }$$

(31)

无量纲Weissenberg数

$$Wi = {{{\lambda _1}V}}/{L}\text{, }$$

(32)

以及无量纲溶剂黏度比

$$\beta = {{{\eta _{\rm{s}}}}}/{\eta }.$$

(33)

这里V和L分别表示流体流动的特征速度和特征长度. Wi比较了弹性力与黏性力, 通常由流体的应力松弛时间与具体的加工时间的关系给出.

图1给出了黏弹性Poiseuille流的几何模型和初始状态. 当$t = 0$时, 流体位于两块水平固定且间距为L的无限长平板之间, 流体初始速度为0. 当$t \gt 0$时, 流体受到水平外力F的作用开始流动. 模拟中, 沿x轴方向采用周期边界条件, 故计算区域可取为${L_x} \times {L_y} = 0.5 \times 1$. 水平外力$F = 1.0$, 粒子密度$\rho = 1$, 各流变参数取值如下: Re = 2, Wi = 1,β = 0.1, α = 0.01, γ = 0.9, Q = 4, 其中特征长度L = 1, 特征速度$V = 0.25\;{{\rm{m}}} \cdot {{{\rm{s}}} ^{ - 1}}$. 粒子初始间距为$\Delta x = 0.02\;{{\rm{m}}}$, 对应于1274个流体粒子、64个边界粒子和192个虚拟粒子. 核函数采用5次样条函数, 光滑长度$h = 0.9\Delta x$. 声速$c = 10 V$, 时间步长$\Delta t = 1.0 \times {10^{ - 5}}$. 对于声速$c$, 还选取较大的$20 V$和$40 V$进行数值实验, 所得模拟结果和$10 V$的模拟结果几乎相同, 表明声速$c$取值为特征速度的10倍, 已能够保证弱可压流体的流动行为充分接近不可压流体的行为.

图2(a)给出了利用本文改进SPH方法模拟基于XPP模型的黏弹性Poiseuille流关于y轴的速度分布. 可以看出, 黏弹性Poiseuille流的速度关于y = 0对称; 且由于受到弹性力作用的影响, 速度有明显的过冲现象, 即速度在t = 0.8时过冲到最大, 随后在t = 2.1时降低至最小, 最终在t = 8时达到稳态. 空间点1—3 (如图1所示)的速度随时间的变化情况如图2(b)所示. 很明显, 3个空间点均展现出了速度过冲现象; 且空间点离边界越近, 速度值越小. 此外, 由于空间点1和3关于中心点2对称, 因此它们的速度值相同. 从图2(c)可知, 空间点2的弹性剪切应力一直为0, 而空间点1和3的弹性剪切应力大小相等, 方向相反. 显然, 空间点离边界越近, 弹性剪切应力值越大.

4.1.1   改进SPH方法的有效性和优势

目前, 基于XPP模型的黏弹性Poiseulle流尚没有瞬态解析解. 但当$f\left( {\lambda , \;{{\rm{{\boldsymbol{\tau}}}} _i}} \right) = 1$和$\alpha = 0$时, XPP本构方程将退化为Oldroyd-B本构方程, 而基于Oldroyd-B模型的黏弹性Poiseuille流具有瞬态解析解^[33]. 因此, 为了验证本文改进SPH方法模拟XPP黏弹性流体的准确性和可靠性, 增加了$f\left( {\lambda , \;{{\rm{{\boldsymbol{\tau}}}} _i}} \right) = 1$和$\alpha = 0$的数值试验, 而其他参数保持不变. 将得到的u数值解与基于Oldroyd-B模型的u解析解进行比较, 如图3所示, SPH结果与解析解吻合很好, 从而验证了本文提出的增强型边界处理技术是合理的, 以及本文改进的SPH方法模拟XPP黏弹性流体是有效的.

为了展现本文改进SPH方法的优势, 这里利用传统SPH方法对$f\left( {\lambda , \;{{\rm{{\boldsymbol{\tau}}}} _i}} \right) = 1$和$\alpha = 0$的数值实验进行了测试, 计算得到的u数值解如图3所示. 同时, 为了更好地展示不同曲线标志的差异性, 还增加了特定区域的局部放大图. 可以看出, 相比于传统SPH方法, 利用本文改进SPH方法得到的u数值解更接近解析解. 另外, 为了定量地体现本文改进SPH方法的优势, 进一步引入L-2范数误差:

$$\varDelta =\sqrt{\frac{{\displaystyle \sum {\left({R}_{解析解}-{R}_{数值解}\right)}^{2}}}{{\displaystyle \sum {\left({R}_{解析解}\right)}^{2}}}} ,$$

(34)

其中, R代表SPH数值解或解析解. 图4给出了利用本文改进SPH方法和传统SPH方法得到的u数值解与解析解的L-2范数误差随时间变化的比较. 很明显, 相比于传统SPH方法, 利用本文改进SPH方法得到的u数值解与解析解的L-2范数误差更小, 这表明了本文改进SPH方法相较于传统SPH方法具有更高的计算精度.

另外, 为了评价本文SPH方法模拟XPP黏弹性流体的数值收敛性, 以图2为例(Re = 2, Wi = 1, β = 0.1, α = 0.01, γ = 0.9, Q = 4), 本文特别增加了粒子初始间距分别为$\Delta x = 0.025$, 0.0125和0.01的三个数值实验, 而其他参数均保持不变. 所有数值试验的时间步长均设置为$\Delta t = 1.0 \times {10^{ - 5}}$, 以保证数值稳定性. 图5(a)和图5(b)分别给出了利用不同粒子初始间距$\Delta x$得到的点2处的速度和点1处的弹性剪切应力随时间的变化情况. 可以看出, 利用不同粒子初始间距$\Delta x$得到的点2处的速度和点1处的弹性剪切应力随时间的变化均相同, 从而表明了本文SPH方法具有良好的数值收敛性.

4.1.2   XPP流变参数对流动过程的影响

下面讨论XPP本构模型中各流变参数对流动过程的影响. 由于空间点离边界越近, 速度值越小, 弹性剪切应力值越大, 故讨论各流变参数对速度的影响时, 只考虑中心点2; 当讨论各流变参数对弹性剪切应力的影响时, 只考虑中心点1.

1) Re的影响

Re对速度u和弹性剪切应力τ_xy的影响如图6(a)和图7(a)所示. 由图6(a)可见, Re越大, 速度过冲达到的最值越大, 达到最值所需的时间越长. 对于Re = 1, 2, 3, 4和5, 速度过冲达到的最值分别为0.385, 0.575, 0.710, 0.890和1.000, 达到最值所需的时间分别为0.510, 0.575, 1.105, 1.255和1.450. 另外, 能观察到速度过冲现象随着Re的增大而减弱. 当Re = 2时, 有2次明显的速度过冲现象; 但当Re增大到4以上时, 却只观察到了1次, 这是因为Re越大, 流体惯性力越大, 弹性力相对变小. Re对速度的稳态值也有较大的影响. Re越大, 速度达到的稳态值越大. Re = 1时的速度稳态值为0.115, 明显低于Re = 4和5时的速度稳态值(分别为0.625和0.920). 对于弹性剪切应力: Re越大, 弹性剪切应力过冲现象越弱, 但过冲达到的最值越大, 达到最值所需的时间越长. 此外, Re对弹性剪切应力稳态值的影响不大, 表明Re并非影响弹性力的主要因素.

2) Wi的影响

Wi是黏弹性流动的一个重要参数. 图6(b)和图7(b)分别描述了Wi对速度u和弹性剪切应力τ_xy的影响. 由图6(b)可见, Wi越大, 速度过冲达到的最值越大, 达到最值所需的时间越长. 对于Wi = 0.1, 0.2, 0.5, 1.0, 2.0, 3.0, 4.0和5.0, 速度过冲达到的最值分别为0.220, 0.310, 0.425, 0.530, 0.715, 0.815, 0.900和0.985, 达到最值所需的时间分别为0.355, 0.410, 0.520, 0.750, 0.955, 1.115, 1.425和1.655. 另外, 速度过冲现象随着Wi的增大先增强后减弱. 从图6(b)可以看出, Wi = 1时速度过冲现象最为显著, 发生有2次. 低于或高于该值, 速度过冲现象均减弱. 当Wi 降到0.2以下时, 仅观察到1次, 这是因为Wi越小, 弹性力越小. 当Wi 增到4以上时, 也仅观察到1次, 这是因为XPP模型可合理地描述流体的剪切变稀行为; 当Wi较大时, 流体发生了剪切变稀. Wi对速度稳态值也有显著的影响, 即Wi越大, 速度达到的稳态值越大. 由图7(b)可见, 由于剪切变稀, 随着Wi的增大, 弹性剪切应力的稳态值反而变小, 但过冲达到的最值先增大后减小, 达到最值所需的时间越长.

3) β的影响

β表示黏弹性流动中牛顿溶剂黏度与总黏度的比值. 图6(c)和图7(c)分别给出了β对速度u和弹性剪切应力τ_xy的影响. 由图6(c)可见, 随着β的增大, 速度过冲现象逐渐变弱. 这是因为β越大, 牛顿溶剂黏度越高, 流体越接近于牛顿流体. 对于β = 0.1, 0.3, 0.5, 0.7和0.9, 速度过冲达到的最值分别为0.565, 0.430, 0.365, 0.310和0.235. 不同β值达到最值所需的时间差异不大, 稳态值也近乎相同. 由图7(c)可见, β极大地影响着弹性剪切应力的稳态值. β越小, 弹性剪切应力的稳态值越大, 这是因为β越小, 聚合物黏度越高, 弹性力越大. 可见, β是影响弹性力大小的一个主要因素.

4) α的影响

α是黏弹性溶液中与各向异性拖曳力相关的一个流变参数, 它对速度u和剪切应力τ_xy的影响如图6(d)和图7(d)所示. 由图6(d)可见, α并不影响速度过冲达到的最值和达到最值所需的时间, 但对速度达到的稳态值产生一定的影响, 即α越大, 速度稳态值越大. 由图7(d)可见, α对弹性剪切应力过冲行为和稳态值的影响均不大.

5) γ的影响

图6(e)和图7(e)分别描述了γ对速度u和弹性剪切应力τ_xy的影响. 可以看出, γ并不影响速度过冲达到的最值, 但在一定程度上影响弹性剪切应力过冲达到的最值. γ越大, 弹性剪切应力过冲达到的最值越大. 这是因为γ越大, 流体拉伸松弛时间越长. γ对速度达到的稳态值有一定的影响, 即γ越大, 速度稳态值越小. γ也轻微影响弹性剪切应力达到的稳态值. 随着γ的增大, 弹性剪切应力稳态值轻微地增大.

6) Q的影响

Q表示聚合物分子链的臂数. Q对速度u和剪切应力τ_xy的影响如图6(f)和图7(f)所示. 由图6(f)可见, Q并不影响速度过冲达到的最值和达到最值所需的时间, 但它轻微影响速度达到的稳态值, 即Q越大, 速度稳态值轻微地变小. 另外, 由图7(f)可知, Q也轻微地影响弹性剪切应力的过冲行为, 但其对应力稳态值的影响并不大.

4.2   黏弹性液滴撞击固壁

接下来, 研究基于XPP模型的黏弹性液滴撞击固壁问题, 该问题含有复杂的自由面. 值得注意的是, 对于含自由面的黏弹性流体的数值模拟, 若直接使用SPH方法进行模拟, 则会产生张力不稳定性(见图8), 并最终导致模拟过程的中断. 为了解决该问题, 我们特别采用了第3.3节介绍的人工应力.

图9给出了液滴撞击固壁问题的计算模型. 初始时刻的液滴为圆形, 直径${d_0} = 0.02\;{{\rm{m}}}$, 圆心位于(0 m, 0.04 m). 固壁几何尺寸取为$\left\{ (x, y): - 0.3 \leqslant x \leqslant 0.3, \;y = 0 \right\}$. 当t = 0时, 液滴具有初始下落速度$V = 1\;{{\rm{m}}} \cdot {{{\rm{s}}} ^{ - 1}}$; 当t > 0时, 液滴在重力加速度g的作用下开始加速下落, 然后撞击固壁, 发生铺展变形. 对于该问题, 其铺展变形主要受初始下降速度和重力作用的驱动, 表面张力对其影响很小^[21], 因此本文并不考虑表面张力.

模拟中, 液滴密度$\rho = 1000\;{{\rm{kg}}} {\cdot} {{{\rm{m}}} ^{ - 3}}$, 取向松弛时间${\lambda _1} = 0.02{\text{ }}{{\rm{s}}}$. 液滴总黏度$\eta = 4{\text{ }}{{\rm{Pa}}} {\cdot} {{{\rm{s}}} ^{ - 1}}$, 其中牛顿溶剂黏度${\eta _{\text{s}}} = 0.4{\text{ Pa}} {\cdot} {{\text{s}}^{ - 1}}$, 聚合物黏度${\eta _{\text{p}}} = 3.6{\text{ }}{{\rm{Pa}}} {\cdot} {{\rm{s}}^{ - 1}}$. 分别选取液滴直径${d_0}$和下降速度$V$作为特征长度和特征速度, 则无量纲Reynolds数Re = 5, Weissenberg 数Wi = 1, 溶剂黏度比β = 0.1. 其他XPP流变参数为: α = 0.01, γ = 0.9, Q = 4. 为了解决张力不稳定性, 人工应力参数$\varepsilon = 0.5$. 粒子初始间距设置为Δx = 0.0002, 对应于7845个流体粒子、301个边界粒子和903个虚拟粒子. 核函数采用5次样条函数, 光滑长度$h = 1.5\Delta x$. 声速$c = 12.5{\text{ }}{{\rm{m}}} {\cdot} {{{\rm{s}}} ^{ - 1}}$, 时间步长$\Delta t = 1.0 \times {10^{ - 6}}{\text{ }}{{\rm{s}}}$. 对于声速$c$, 分别取$c = 20{\text{ }}{{\rm{m}}} {\cdot} {{{\rm{s}}} ^{ - 1}}$和$40{\text{ }}{{\rm{m}}} {\cdot} {{{\rm{s}}} ^{ - 1}}$进行数值实验, 所得模拟结果和$12.5{\text{ }}{{\rm{m}}} {\cdot} {{{\rm{s}}} ^{ - 1}}$的模拟结果没有显著差异, 表明声速$c$取值为$12.5{\text{ }}{{\rm{m}}} {\cdot} {{{\rm{s}}} ^{ - 1}}$, 已能够保证弱可压流体的流动行为充分接近不可压流体的行为. 对于该算例, 本文在联想深腾6800型高性能服务器上, 使用1个CPU (Intel Xeon Platinum 8273 @ 3.0 GHz)模拟20 万个时间步, 所耗时间约为5.3 h.

图10给出了利用本文改进SPH方法模拟基于XPP模型的黏弹性液滴在撞击固壁之后6个不同时刻的结果. 值得注意的是, 本文选用无量纲时间$T = tV/{d_0}$来记录液滴的铺展历程. 可以看出, 液滴在撞击固壁之后的流动过程分为三个阶段: 第一阶段, 液滴在撞击固壁后迅速铺展. 该阶段的液滴虽具有初始向下的速度, 但在撞击固壁之后, 液滴左半部分流速向左, 右半部分流速向右(见T = 1.55和1.85). 第二阶段, 液滴铺展到最大宽度后, 受弹性力的作用逐渐开始收缩. 在该阶段, 液滴的速度方向发生了转变, 即液滴左半部分流速向右, 右半部分流速向左(见T = 2.35, 2.85和3.45). 第三阶段, 液滴收缩到最小宽度后, 由于能量的耗散和重力的影响, 再次缓慢铺展(见T = 5.00).

图11给出了XPP液滴的铺展宽度随时间的变化. 为了便于分析, 这里也用液滴直径d₀对液滴铺展宽度d(T)进行无量纲处理. 从图11可清晰地观察到液滴撞击固壁之后的三个流动阶段, 其中液滴在T约为2.405时达到最大铺展宽度(～2.255); 在T约为4.015时收缩到最小宽度(～1.875). 对于该算例, 由于引入了人工应力, 因此为了验证本文改进SPH方法模拟含自由面的黏弹性流动问题的有效性, 这里也展示了利用FDM方法^[34]模拟该问题的数值结果(见文献[34]中的图9(a)), 其中自由界面的追踪通过标记单元法来实现. 从图11可以看出, 利用本文改进SPH方法得到的结果和FDM解是基本符合的, 从而表明改进SPH方法在解决张力不稳定性的同时, 能够得到可靠的模拟结果. 值得注意的是, 两种数值结果在液滴收缩阶段有微小的差异, 可能归因于: FDM方法使用Possion方程求解压力, 而本文考虑弱可压缩流体, 使用状态方程求解压力. 为了评价本文改进SPH方法的数值收敛性, 还特别增加了粒子初始间距分别为Δx = 0.0001, 0.00025和0.0004的数值试验, 而其他参数均保持不变. 得到的XPP液滴铺展宽度随时间的变化曲线如图11所示. 可以看出, 通过3组不同粒子初始间距和原始粒子初始间距(Δx = 0.0002)得到的液滴铺展宽度结果基本一致, 从而验证了本文改进SPH方法模拟黏弹性自由表面流是数值收敛的.

对于该问题, 由于液滴不仅具有初始下落速度, 而且受重力的作用加速下落, 因此液滴会对固壁边界产生较大的冲击力. 然而, 由于液滴较大的黏度($\eta = 4{\text{ }}{\rm Pa} {\cdot} {{\rm s} ^{ - 1}}$), 因此在撞击壁面过程中并没有出现明显的小液滴飞溅. 模拟过程中, 也并未发现任何流体粒子流出固壁边界(见图10), 表明本文提出的增强型边界处理技术能够有效地防止流体粒子穿透边界. 此外, SPH结果与FDM解^[34]是基本吻合的(见图11), 这进一步表明了本文提出的增强型边界处理技术是合理的. 总之, 相较于传统的排斥力法和虚拟粒子法, 本文为SPH方法的边界处理提供了一种合理的、可供替代的边界处理方案, 且该方案具有程序设计简单、避免构造伪粒子进行Shepard插值繁琐操作的优势.

下面讨论XPP本构模型中各流变参数对液滴铺展宽度随时间的变化的影响.

7) Re的影响

Re对液滴铺展宽度的影响如图12(a)所示. 可以看到, 随着Re的增大, 液滴的最大铺展宽度显著增加, 达到最大铺展宽度所需的时间相对滞后. 对于Re = 2, 5和10, 液滴的最大铺展宽度分别约为1.735, 2.215和2.705, 达到最大铺展宽度所需的时间分别约为1.850, 2.005和2.495. 这是因为Re越大, 液滴惯性力越大, 导致液滴撞击固壁后有更快的铺展行为. 另外, Re越大, 液滴的收缩行为越弱. 从图12(a)可知, Re = 2的收缩行为明显强于Re = 10的收缩行为. 这是因为Re越大, 液滴弹性力相对越小. 由于Re越大, 液滴的铺展越快, 导致液滴的最小收缩宽度和最后铺展宽度也相应地更大.

8) Wi的影响

图12(b)展示了Wi对液滴铺展宽度的影响. 可以看出, 增加Wi导致液滴达到的最大铺展宽度显著增加, 达到最大铺展宽度所需的时间也相对滞后. 对于Wi = 0.5, 1, 2, 4和8, 液滴达到的最大铺展宽度分别约为2.005, 2.215, 2.455, 2.700和2.875, 达到最大铺展宽度的时间分别约为2.015, 2.355, 2.500, 2.765和2.955. 另外, 随着Wi的增大, 液滴的最小收缩宽度相应增加. 不同Wi得到的液滴最小收缩宽度分别约为1.805, 1.845, 1.975, 2.110和2.400. 这是因为, 随着Wi的增大, 黏弹性流体的剪切稀化增强, 导致液滴的铺展更快, 因此液滴的铺展宽度相应增加.

9) β的影响

β对液滴铺展宽度的影响如图12(c)所示, β越大, 液滴的收缩行为越弱, 但最终铺展宽度相差不大. 特别是, 对于β = 0.9时的液滴已观察不到收缩现象. 这是因为随着β的增大, 牛顿溶剂黏度增加, 聚合物黏度降低, 导致液滴的弹性力变小. 在这种情形下, 液滴表现出类似于牛顿流体的铺展行为. 另外, 改变β只改变牛顿溶剂黏度占总黏度的比, 而不改变总黏度的大小, 因此液滴的最终铺展宽度相差不大.

10) α的影响

各向异性流变参数α对液滴铺展宽度的影响如图12(d)所示. 可以观察到, 改变α并不影响液滴达到的最大铺展宽度, 但在一定程度上影响液滴的收缩行为, 即α越大, 液滴的收缩行为越弱, 液滴的最小收缩宽度因此相应地增加. 对于α = 0, 0.01, 0.1, 0.3和0.5, 液滴的最小收缩宽度分别约为1.835, 1.837, 1.875, 1.925和1.995. 由于变弱的液滴收缩行为, 液滴的最后铺展宽度随着α的增大而增大.

11) γ的影响

如图12(e)所示, 降低γ导致液滴的最大铺展宽度显著增加, 达到最大铺展宽度所需的时间也相对滞后. 对于γ = 0.9, 0.7, 0.5, 0.3, 0.2和0.1, 液滴的最大铺展宽度分别约为2.215, 2.235, 2.265, 2.335, 2.400和2.495, 达到最大铺展宽度的时间分别约为2.355, 2.357, 2.475, 2.500, 2.515和2.535. 另外, γ在一定程度上也影响着液滴的收缩行为, 即γ越小, 液滴的最小收缩宽度越大, 这是因为降低γ导致液滴的拉伸松弛时间变短, 聚合物分子链得不到有效的拉伸, 因此液滴的弹性力变小, 导致液滴的收缩行为变弱. 由于变弱的液滴收缩行为, 液滴的最后铺展宽度随着γ的减小而增大.

12) Q的影响

图12(f)展示了聚合物分子链臂数Q对液滴铺展宽度的影响. 可以看出, 臂数越多, 对流动行为的影响越弱. 特别是, 当Q > 4时, 已观察不到其对液滴铺展宽度的影响. 另一方面, 臂数越少, 液滴的收缩行为越弱, 达到的最终铺展宽度越大. 对于单臂的情形(Q = 1), 得到的液滴铺展宽度明显高于Q > 4得到的液滴铺展宽度.

5.   结　论

本文提出一种改进SPH方法数值模拟了基于XPP模型的黏弹性流动. 为了提高计算精度, 采用了一种不含核导数计算的核梯度修正离散格式. 为了防止粒子穿透固壁, 提出了一种增强型的边界处理技术. 为了消除张力不稳定性, 将人工应力耦合到动量守恒方程中. 运用改进SPH方法数值模拟了基于XPP模型的黏弹性Poiseuille流和黏弹性液滴撞击固壁问题, 深入分析了不同流变参数对流动过程的影响, 结论如下.

1) 对于黏弹性Poiseuille流, 通过比较SPH数值解和解析解以及评价数值解的收敛性, 表明了本文改进SPH方法模拟XPP黏弹性流动问题是准确、有效的, 且相对于传统SPH方法具有更高的计算精度. 对于黏弹性液滴问题, 利用增强型的边界处理技术可有效地防止粒子穿透固壁; 借助人工应力可有效地消除张力不稳定性.

2) XPP流变参数对黏弹性Poiseuille流动过程有着重要的影响, 即Re, Wi和α越大, 速度稳态值越大; γ和Q越大, 速度稳态值越小; β越大, 速度过冲现象越弱, 但并不影响速度稳态值的大小. Wi和β越大, 弹性剪切应力稳态值越小; γ越大, 弹性剪切应力稳态值轻微增大; Re, α和Q对弹性剪切应力稳态值的影响不大.

3) XPP流变参数对黏弹性液滴撞击固壁流动过程有着重要的影响, 即Re和Wi越大, 液滴铺展越快, 铺展宽度越大. β越大, 液滴收缩行为越弱, 但并不影响液滴的最终铺展宽度. α越大, 液滴收缩行为越弱, 铺展宽度越大. γ越大, 液滴收缩行为越强, 铺展宽度越小. Q越大, 对液滴铺展宽度的影响越弱.

4) 本文改进SPH方法能够有效而准确地描述基于XPP模型的黏弹性流体的复杂流变特性和自由面变化特征.