SIMULATIONS OF NON-ISOTHERMAL VISCOELASTIC COMPLEX FLOWS BY IMPROVED SPH METHOD
-
摘要: 非等温黏弹性流体广泛存在于自然界和工业生产中, 准确预测黏弹性流体的非等温流动机理和复杂流变特性有着重要的应用价值. 文章提出一种改进的光滑粒子流体动力学(smoothed particle hydrodynamics, SPH)方法对非等温黏弹性复杂流动进行了数值模拟, 其中流体的黏弹特性通过eXtended Pom-Pom本构模型来表征. 为了提高模拟结果的精度, 采用了一种核函数梯度的修正算法; 为了灵活地施加边界条件, 发展了边界粒子和虚拟粒子相联合的边界处理方法; 为了消除流动过程中的拉伸不稳定性, 施加了粒子迁移技术. 运用改进SPH方法数值模拟了液滴撞击固壁和F型腔注塑成型问题, 通过与Basilisk软件得到的结果进行比较验证了改进SPH方法求解非等温黏弹性流体的有效性. 通过利用不同粒子初始间距进行计算, 评价了改进SPH方法的数值收敛性. 研究了非等温流动相较于等温流动的不同流动特征, 深入分析了不同热流变参数对流动过程的影响. 数值结果表明, 文章提出的改进SPH方法可稳定、准确地描述非等温黏弹性复杂流动的传热机理、复杂流变特性和自由面变化特性.Abstract: Non-isothermal viscoelastic fluid flow phenomena widely exist in nature and industrial productions, such as oil reservoir engineering, injection molding, etc. These flows generally exhibit a non-isothermal state. Accurate prediction of non-isothermal flow mechanism and complex rheological properties of viscoelastic fluid has important engineering application value. In this paper, an improved smoothed particle hydrodynamics (SPH) method is proposed for the numerical simulation of non-isothermal viscoelastic complex flow, in which the viscoelastic properties of the fluid are characterized by the eXtended Pom-Pom constitutive model. To improve the accuracy of simulation results, a kernel function gradient correction algorithm is adopted. To enforce the boundary conditions flexibly, a boundary treatment method combining boundary particles and virtual particles is developed. To eliminate the tensile instability in the flow process, the particle migration technology is applied. The improved SPH method is used to numerically simulate the impact of a droplet on the solid wall and injection molding of an F-shaped cavity. The effectiveness of the improved SPH method in solving the non-isothermal viscoelastic fluid is verified by comparing the SPH results with those obtained by the Basilisk software. Good agreement between these two numerical solutions is achieved. The numerical convergence of the improved SPH method is evaluated by using several different initial particle spacings. The different flow characteristics of non-isothermal flow compared with isothermal flow are investigated. It is found that the introduction of temperature leads to stronger contraction behavior of droplet. The influences of some different thermal rheological parameters such as the Péclet number, the Reyonlds number, the Weissenberg number, the solvent viscosity ratio, the anisotropy parameter, the relaxation time ratio and the molecular chain arm number on the flow process are deeply analyzed. The numerical results show that the improved SPH method proposed in this paper can accurately and stably describe the heat transfer mechanism, complex rheological properties, and free surface variation characteristics of non-isothermal viscoelastic fluid.
-
引言
黏弹性流体流动现象广泛存在于自然界和工业生产中, 如油藏工程、注塑成型等. 这些流动一般呈现非等温状态, 因此有必要考虑温度模型[1]. 准确预测黏弹性流体的非等温流动机理有着极为重要的工程应用价值, 但由于其多物理场耦合的复杂性和黏弹流变特性的复杂性, 实验和解析求解都较困难, 因此数值模拟是一个很好的选择.
目前, 基于网格的数值方法已被应用于此类问题的数值模拟中. Li等[2]发展了浸入边界法模拟非等温黏弹性熔体充填过程. Fernandes等[3]基于对数构象张量法提出了求解非等温黏弹性层流的分块耦合算法. Moreno等[4]利用变分多尺度有限元公式对非等温黏弹性圆柱绕流进行了模拟与分析. 虽然网格类方法在模拟黏弹性自由表面流方面取得了成功, 但它们需要借助额外的界面追踪技术, 实施过程较为复杂.
光滑粒子流体动力学(smoothed particle hydrodynamics, SPH)[5-6]方法是一种拉格朗日型的无网格数值方法. 与基于网格的数值方法相比, SPH方法具有粒子特性、自适应特性、拉格朗日特性和易于编程等优势, 因此近年来已被用于冲击[7]、爆炸[8]、流固耦合[9-10]、多相流[11]等的数值模拟中. 然而, 传统SPH方法也存在模拟结果精度低、边界条件不易施加和易出现拉伸不稳定性的缺点. 为了提高模拟结果的精度, Liu等[12]推导具有较高计算精度的有限粒子方法; Xue等[13]发展基于黎曼SPH和再生核粒子法的耦合算法; Ng等[14]提出模拟流固耦合的改进SPH-VCPM方法. 为了灵活地施加边界条件, Monaghan等[15]提出只施加一层粒子在固壁上的排斥力法; Morris等[16]发展了布置多层粒子在固壁上的虚拟粒子法; Liu等[17]建立排斥力法和虚拟粒子法相耦合的动力学边界方法; Yildiz等[18]提出多重边界切线法; Chen等[19]发展适用于入水问题的改进边界处理. 为了解决流动过程中出现的拉伸不稳定性, Yang等[20]引入一种新的核函数; Sun等[21]提出多分辨率δ + SPH方法, 解决了高Reynolds数下圆柱绕流等问题的数值不稳定; Chalk等[22]设计了基于应力节点的SPH修正算法, 克服了滑坡的张力不稳定性; Lyu等[23]发展了增强的粒子迁移技术, 消除了溃坝和液面晃荡等问题的粒子分布不均. 虽然这些改进方法已被成功提出并得到应用, 但其在黏弹性流体流动问题中的应用并不多见.
对于黏弹性流体流动问题的SPH模拟, Rafiee等[24]利用不可压SPH方法研究了基于Oldroyd-B模型的黏弹性自由表面流问题. 杨波等[25]运用SPH方法数值模拟了基于Phan-Thien-Tanner本构模型的黏弹性泊肃叶流动和液滴撞击固壁面, 分析了Reynolds数、Weissenberg数等对流动过程的影响. Xu等[26]发展一种改进SPH方法用于模拟黏弹性液滴撞击固壁面和挤出胀大问题. King等[27]将对数构象公式与应力分裂技术相结合, 提出模拟高Weissenberg数的不可压缩SPH算法. Duque-Daza等[28]将SPH与流体流动模拟的等效弹性势能相结合, 建立黏弹性流体建模与计算的简化方法. Vahabi等[29]用弱可压缩SPH方法研究Oldroyd-B黏弹性溶液中两个初始圆形的气泡在上升过程中的相互作用. Moinfar等[30]分析Giesekus液滴在简单剪切流中的形变. 值得注意的是, 以上研究工作大多基于Oldroyd-B, Giesekus等模型, 且是在等温条件下进行求解的.
本文侧重于对基于eXtended Pom-Pom (XPP)模型的黏弹性流体进行非等温数值模拟. 一方面, XPP模型是从支化聚合物分子理论推导得出的, 因此相较于Oldroyd-B等模型, 能够更合理地描述黏弹性流体的剪切稀化行为. 另一方面, 随着高聚物成型技术的发展, 黏弹性流体流变特性的研究已成为工业生产中的一个重要环节, 而工业生产中的流动通常都是非等温流动, 因此有必要研究非等温黏弹性流动. 对于非等温黏弹性流动的数值模拟, 涉及流体控制方程与本构方程和温度方程的耦合求解以及自由界面的精确捕捉, 因此相关研究较少, 而本文极大地丰富了SPH方法在非等温非牛顿流体力学领域中的应用.
本文首先对传统SPH方法进行改进. 为了提高模拟结果的精度, 采用一种核函数梯度的修正算法; 为了灵活地施加边界条件, 发展了边界粒子和虚拟粒子相联合的边界处理方法; 为了消除黏弹性流体的拉伸不稳定性, 施加粒子迁移技术. 运用改进SPH方法数值模拟了液滴撞击固壁和F型腔注塑成型问题, 通过与Basilisk软件结果进行比较验证了改进SPH方法求解非等温黏弹性流体的有效性. 利用不同粒子初始间距进行计算, 评价改进SPH方法的数值收敛性. 研究非等温流动相较于等温流动的不同流动特征, 分析不同热流变参数对流动过程的影响.
1. 控制方程
非等温黏弹性流体的控制方程包括连续性方程、动量守恒方程和温度方程. 其在拉格朗日坐标系下可分别表示为[31]
$$\qquad\qquad \frac{{{\rm{d}}\rho }}{{{\rm{d}}t}} = - \rho \nabla \cdot {\boldsymbol{u}} $$ (1) $$\qquad\qquad \frac{{{\rm{d}}{\boldsymbol{u}}}}{{{\rm{d}}t}} = \frac{1}{\rho }\nabla \cdot {\boldsymbol{\sigma }} + {\boldsymbol{F}} $$ (2) $$\qquad\qquad \frac{{{\rm{d}}T}}{{{\rm{d}}t}} = \frac{1}{{\rho {c_p}}}\nabla \cdot \left( {\kappa \nabla T} \right) $$ (3) 式中,
${\rm{d}}()/{\rm{d}}t$ 为物质导数,$\rho $ 为流体密度,$t$ 为时间,${\boldsymbol{u}}$ 为速度,${\boldsymbol{\sigma }}$ 为总应力张量,${\boldsymbol{F}}$ 为外力,$T$ 为温度,${c_p}$ 为比热容,$\kappa $ 为导热系数.总应力张量
${\boldsymbol{\sigma }}$ 可分解为各向同性压力、溶剂贡献和聚合物贡献之和$$ {\boldsymbol{\sigma }} = - p{\boldsymbol{I}} + 2{\eta _s}{\boldsymbol{d}} + {\boldsymbol{\tau }} $$ (4) 其中,
${\boldsymbol{I}}$ 为单位矩阵,$ {\eta _s} $ 为溶剂黏度,${\boldsymbol{\tau }}$ 表示弹性偏应力张量,${\boldsymbol{d}}$ 为形变率张量$$ {\boldsymbol{d}} = \frac{1}{2}\left[ {\nabla {\boldsymbol{u}} + {{\left( {\nabla {\boldsymbol{u}}} \right)}^{\rm{T}}}} \right] $$ (5) 1.1 非等温XPP本构方程
为了封闭控制方程, 需结合关于弹性偏应力张量
${\boldsymbol{\tau }}$ 的本构方程. 本文考虑非等温XPP模型, 其本构方程为[32]$$ \begin{split} & f\left( {\lambda {\text{ }},{\text{ }}{\boldsymbol{\tau }}} \right){\boldsymbol{\tau }} + {\lambda _1}\mathop {\boldsymbol{\tau }}\limits^\nabla + {G_0}\left[ {f\left( {\lambda {\text{ }},{\text{ }}{\boldsymbol{\tau }}} \right) - 1} \right]{\boldsymbol{I}}+ \\ &\qquad \frac{\alpha }{{{G_0}}}{\boldsymbol{\tau }} \cdot {\boldsymbol{\tau }} = 2{\lambda _1}{G_0}{\boldsymbol{d}} \end{split} $$ (6) 其中, 上随体导数
$$ \mathop {\boldsymbol{\tau }}\limits^\nabla = \frac{{{\rm{d}}{\boldsymbol{\tau }}}}{{{\rm{d}}t}} - {\boldsymbol{\tau }} \cdot \left( {\nabla {\boldsymbol{u}}} \right) - {\left( {\nabla {\boldsymbol{u}}} \right)^{\rm{T}}} \cdot {\boldsymbol{\tau }} $$ (7) 由式(6)和式(7), 可得
$$ \begin{split} & \frac{{{\rm{d}}{\boldsymbol{\tau }}}}{{{\rm{d}}t}} = {\boldsymbol{\tau }} \cdot \left( {\nabla {\boldsymbol{u}}} \right) + {\left( {\nabla {\boldsymbol{u}}} \right)^{\rm{T}}} \cdot {\boldsymbol{\tau }}-\frac{{f\left( {\lambda {\text{ }},{\text{ }}{\boldsymbol{\tau }}} \right)}}{{{\lambda _1}}}{\boldsymbol{\tau }} - \\ &\qquad \frac{{{G_0}}}{{{\lambda _1}}}\left[ {f\left( {\lambda {\text{ }},{\text{ }}{\boldsymbol{\tau }}} \right) - 1} \right]{\boldsymbol{I}} - \frac{\alpha }{{{\lambda _1}{G_0}}}{\boldsymbol{\tau }} \cdot {\boldsymbol{\tau }} + 2{G_0}{\boldsymbol{d}} \end{split}$$ (8) 其中
$$ f\left( {\lambda {\text{ }},{\text{ }}{\boldsymbol{\tau }}} \right) = 2\frac{{{\lambda _1}}}{{{\lambda _2}}}{{\rm{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] $$ (9) 式中,
$\lambda = \sqrt {1 + \dfrac{1}{{3{G_0}}}{\text{tr}}\left( {\boldsymbol{\tau }} \right)}$ 表示分子链拉伸量,${G_0}$ 表示线性松弛模量,${\lambda _1}$ 和${\lambda _2}$ 分别表示分子链的取向松弛和拉伸松弛时间,$\alpha $ 为各向异性流变参数,${Q_0}{\text{ = }}2/Q$ , 其中Q表示分子链臂数. 引入$\gamma = {\lambda _2}/{\lambda _1}$ , 表示拉伸松弛与取向松弛时间之比. 定义${\eta _p} = {G_0}{\lambda _1}$ , 表示聚合物黏度, 于是流体总黏度可表示为溶剂黏度和聚合物黏度之和, 即$$ \eta = {\eta _s} + {\eta _p} $$ (10) 本文考虑非等温黏弹性流动, 因此流体黏度和聚合物分子链的松弛时间会受到温度变化的影响. 对于流体黏度和松弛时间的温度依赖, 本文依据时温等效原理进行建模和计算, 并选用如下的Reynolds指数模型进行求解[33]
$$\qquad\qquad {\eta _s}\left( T \right){\text{ = }}{\eta _s}{{\rm{e}}^{ - \varphi \left( {T - {T_0}} \right)}} $$ (11) $$\qquad\qquad {\eta _p}\left( T \right){\text{ = }}{\eta _p}{{\rm{e}}^{ - \varphi \left( {T - {T_0}} \right)}} $$ (12) $$\qquad\qquad {\lambda _1}\left( T \right) = {\lambda _1}{{\rm{e}}^{ - \varphi \left( {T - {T_0}} \right)}} $$ (13) $$\qquad\qquad {\lambda _2}\left( T \right) = {\lambda _2}{{\rm{e}}^{ - \varphi \left( {T - {T_0}} \right)}} $$ (14) 其中,
$\varphi $ 为温度依赖系数,${T_0}$ 为流体的参考温度. 另外, 为了更好地描述非等温流动物理和复杂流变特性, 引入以下无量纲数$$ \begin{split} & Pe = \frac{{\rho {c_p}VL}}{\kappa },\;Re = \frac{{\rho VL}}{\eta } \\ & Wi = \frac{{{\lambda _1}V}}{L},\;\beta = \frac{{{\eta _s}}}{\eta } \end{split} $$ (15) 其中, Péclet数Pe表征对流作用与扩散作用的相对大小, Reynolds数Re是运动物体上的惯性力与黏性力之间的比率, Weissenberg数Wi比较了弹性力与黏性力, 通常由流体的应力松弛时间与应变速率的乘积关系给出, 溶剂黏度比
$\beta $ 表示黏弹性流动中牛顿溶剂黏度与流体总黏度的比值, V和L分别表示黏弹性流动的特征速度和特征长度.1.2 状态方程
在SPH方法中, 有两种常用的求解控制方程的方法: 一种是不可压缩SPH 方法[34], 其中流体压力通过Poisson方程隐式计算; 另一种是弱可压缩SPH方法[26], 其中压力根据状态方程显式计算. 本文采用后者, 使用的状态方程为[26]
$$ p(\rho ) = {c^2}(\rho - {\rho _0}) $$ (16) 其中,
${\rho _0}$ 为流体参考密度,$c$ 为声速. 为了保证弱可压缩流体的流动行为足够接近不可压缩流动行为, 声速$c$ 通常选取为不低于流动特征速度的10倍左右.2. 光滑粒子动力学方法
2.1 SPH离散
在SPH方法中, 对流体控制方程有多种不同的离散方案[35]. 本文对连续性方程和动量守恒方程采用的SPH离散格式分别为
$$ \frac{{{\rm{d}}{\rho _i}}}{{{\rm{d}}t}} = \sum\limits_j {{m_j}\left( {v_i^\beta - v_j^\beta } \right)\frac{{\partial {W_{ij}}}}{{\partial x_i^\beta }}} $$ (17) $$ \frac{{{\rm{d}}v_i^\alpha }}{{{\rm{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 } $$ (18) 其中
$$ \sigma _i^{\alpha \beta } = - p{\delta ^{\alpha \beta }} + {\eta _s}\left( {k_i^{\alpha \beta } + k_i^{\beta \alpha }} \right) + \tau _{p,\;i}^{\alpha \beta } $$ (19) $$ 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 }}} $$ (20) 式中, i和j表示粒子编号, m表示粒子质量,
${W_{ij}} = W(\left| {{x_i} - {x_j}} \right|,h)$ 是核函数,$ h $ 是核函数影响区域的光滑长度. 对于核函数, 本文均采用5次样条函数进行计算.对于温度方程, 其等号右端涉及拉普拉斯算子的离散化. 借鉴Shao等[36]对黏性项拉普拉斯算子离散的思想, 本文对温度方程采用如下的SPH离散格式
$$ {\left( {\frac{{{\rm{d}}T}}{{{\rm{d}}t}}} \right)_i} = \frac{\kappa }{{{c_p}}}\sum\limits_j {\frac{{8{m_j}}}{{{{({\rho _i} + {\rho _j})}^2}}}} \cdot \frac{{{{\boldsymbol{r}}_{ij}} \cdot {\nabla _i}{W_{ij}}}}{{{{\left| {{{\boldsymbol{r}}_{ij}}} \right|}^2}}}({T_i} - {T_j}) $$ (21) 对于XPP本构方程, 本文采用的SPH离散格式为
$$ \begin{split} & {\left( {\frac{{{\rm{d}}\tau _p^{\alpha \beta }}}{{{\rm{d}}t}}} \right)_i} = k_i^{\alpha \gamma }\tau _{p,\;i}^{\gamma \beta } + k_i^{\beta \gamma }\tau _{p,\;i}^{\gamma \alpha } - \frac{1}{{{\lambda _1}}}f\left( {\lambda ,\;{\boldsymbol{\tau }}} \right)\tau _{p,\;i}^{\alpha \beta } - \\ &\qquad \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 _{p,\;i}^{\alpha \gamma } \cdot \tau _{p,\;i}^{\gamma \beta }} \right) + \\ &\qquad {G_0}\left( {k_i^{\alpha \beta } + k_i^{\beta \alpha }} \right) \end{split} $$ (22) 另外, 传统SPH方法的计算精度低, 限制了其向更复杂流动问题的应用. 为了提高SPH方法的计算精度, 本文采用了如下的核函数梯度修正算法[26,37]
$$ \left( {\begin{array}{*{20}{c}} {\dfrac{{\mathop {\partial {W_{ij}}}\limits^ \sim }}{{\partial {x_i}}}} \\ {\dfrac{{\mathop {\partial {W_{ij}}}\limits^ \sim }}{{\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) $$ (23) 其中, 核函数梯度修正矩阵M为
$$ {\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) $$ (24) 关于核函数梯度修正算法的更多内容及提高模拟结果精度的相关验证, 可参阅文献[26, 37].
2.2 边界处理
固壁边界的处理是SPH计算实施过程的一个重要部分, 其处理是否得当, 关乎 SPH模拟的效率和稳定性. 为了灵活地施加边界条件, 发展了边界粒子和虚拟粒子相联合的边界处理方法. 对于直边固壁, 在固壁边界上以粒子初始间距
$\Delta x$ 布置一层边界粒子, 在固壁边界外以粒子初始间距$\Delta x$ 布置3层虚拟粒子. 虚拟粒子看作是边界粒子在边界外法向上的垂直延伸. 边界粒子和虚拟粒子的温度在计算过程中均保持不变. 边界粒子不再对邻近边界处的流体粒子施加排斥力; 相反, 其参与控制方程中温度、速度、压力和弹性应力的计算. 计算过程中, 边界粒子的密度和位置保持不变. 速度设置为0, 以施加无滑移边界条件. 压力和弹性应力通过对相邻流体粒子进行Shepard插值得到$$ {A_i} = \frac{{\displaystyle\sum\limits_j {\frac{{{m_j}}}{{{\rho _j}}}{A_j}{W_{ij}}} }}{{\displaystyle\sum\limits_j {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}} }} $$ (25) 式中, A表示压力或弹性应力,
$i$ 为边界粒子,$j$ 为与粒子$i$ 相邻的流体粒子. 由于虚拟粒子是边界粒子在边界外法向上的垂直延伸, 因此每个虚拟粒子均与其边界法向上的边界粒子相对应. 为满足压力和弹性应力的纽曼边界条件, 虚拟粒子的压力和弹性应力设置与边界法向上对应的边界粒子值相同. 速度设置为0, 以施加无滑移边界条件.2.3 粒子迁移技术
进行含自由面的黏弹性复杂流动的SPH模拟, 极易出现流动的拉伸不稳定. 为了解决该问题, 施加了粒子迁移技术[38]. 该技术的思想是在控制方程求解之后, 对粒子的位移进行微小地迁移, 并通过Taylor级数展开对相关流体变量进行修正[38]
$$ {\phi _i}^\prime = {\phi _i} + \delta {{\boldsymbol{r}}_{ii'}} \cdot {(\nabla \phi )_i} + O(\delta {\boldsymbol{r}}_{ii'}^2) $$ (26) 其中
$i$ 和$i'$ 分别是粒子迁移前和迁移后的位移,$\delta {{\boldsymbol{r}}_{ii'}}$ 为粒子迁移前后位移的向量.粒子迁移向量由菲克扩散定律[38]控制, 即
$$ J = - D'\nabla C $$ (27) 其中,
$J$ 是扩散通量,$D'$ 是扩散系数,$C$ 是粒子浓度, 定义为$$ {C_i} = \sum\limits_j {\frac{{{m_j}}}{{{\rho _j}}}{W_{ij}}} $$ (28) 假定扩散通量
$J$ 与粒子迁移速度${{\boldsymbol{v}}_s}$ 成正比, 则粒子迁移向量$\delta {{\boldsymbol{r}}_i} = {{\boldsymbol{v}}_s}\Delta t$ . 于是$$ \delta {{\boldsymbol{r}}_i} \propto - {D_i}^\prime \nabla {C_i}\Delta t $$ (29) 对于扩散系数, 其上限可由对流扩散方程的Von Neumann稳定性分析得到
$$ {D_i}^\prime \leqslant \frac{1}{2}\frac{{{h^2}}}{{\Delta t}} $$ (30) 其中,
$\Delta t$ 表示时间步长. 时间步长$\Delta t$ 需满足CFL条件$$ \Delta t \leqslant 0.1\frac{h}{{\left| {{{\boldsymbol{v}}_i}} \right|}} $$ (31) 其中,
$ \left| {{{\boldsymbol{v}}_i}} \right| $ 表示第i个粒子的速度. 本文选用$\Delta t = 0.1\dfrac{h}{{\left| {{{\boldsymbol{v}}_i}} \right|}}$ , 并将其代入式(30)有$$ {D_i}^\prime \leqslant 5h\left| {{{\boldsymbol{v}}_i}} \right| $$ (32) 由式(29)和式(32), 可得粒子迁移向量[38]
$$ \delta {{\boldsymbol{r}}_i} = - 5h\left| {{{\boldsymbol{v}}_i}} \right|\nabla {C_i}\Delta t $$ (33) 对于自由表面流, 还需处理自由表面粒子的扩散. 为此, 需要首先检测出自由表面粒子, 算法如下: (1) 通过密度法
$ \rho _i^* < 0.97{\rho _0} $ 检测出潜在的自由表面粒子$i$ , 其中$ \rho _i^* $ 可由密度求和公式求得; (2) 对于在第1步标记的潜在的自由表面粒子$i$ , 进一步判断其法向向量的扇形区域是否存在邻近粒子. 如果存在, 则粒子$i$ 被视为内部粒子; 否则, 粒子$i$ 被视为最终的自由表面粒子. 依据以上两个步骤, 只有最外一层的流体粒子被检测出来, 被视为自由表面粒子. 这些自由表面粒子代表了流体的自由表面. 一旦自由表面粒子被检测出来, 可通过菲克修正扩散定律[38]来控制自由表面粒子的迁移向量$$ \delta {{\boldsymbol{r}}_i} \propto - {D_i}^\prime ({\boldsymbol{I}} - {{\boldsymbol{n}}_i} \otimes {{\boldsymbol{n}}_i})\nabla {C_i}\Delta t $$ (34) 其中,
$ {{\boldsymbol{n}}_i} $ 是第i个粒子的单位法向.由式(32)和式(34), 可得自由表面粒子的迁移向量为[38]
$$ \delta {{\boldsymbol{r}}_i} = - 5h\left| {{{\boldsymbol{v}}_i}} \right|({\boldsymbol{I}} - {{\boldsymbol{n}}_i} \otimes {{\boldsymbol{n}}_i})\nabla {C_i}\Delta t $$ (35) 通过对流体内部和自由表面粒子的位置进行微小地位移, 可有效地消除流动过程中出现的拉伸不稳定性. 值得注意的是, 虽然粒子迁移技术已被Lind等[38]提出, 并得到一定的应用[23,39-40], 但绝大多数工作研究无黏性流或牛顿流. 目前, 将粒子迁移技术应用到黏弹性复杂流动问题的文献报道并不多见. 因此, 本文研究工作可看作是粒子迁移技术应用到复杂非等温黏弹性自由表面流问题的成功尝试.
2.4 时间积分
由于蛙跳积分法具有存储需求量低、计算效率高的优点, 因此本文选用该方法进行时间积分. 关于该积分格式的详细内容, 可参阅文献[35]. 另外, 为了保持数值稳定性, 时间步长
$ \Delta t $ 需要满足以下条件$$ \Delta t \leqslant \min \left( {\frac{h}{c},{\text{ }}\mathop {{\text{min}}}\limits_i \sqrt {\frac{h}{{{F_i}}}} ,{\text{ }}0.125\frac{{{h^2}}}{{{\upsilon _0}}}} \right) $$ (36) 其中,
$ {F_i} $ 表示每单位质量的力,$ {\upsilon _0} = \eta /\rho $ 表示流体的运动黏度.3. 数值算例
3.1 液滴撞击固壁
本文选取的第1个数值算例为基于XPP本构模型的非等温黏弹性液滴撞击固壁问题, 该问题含有复杂的自由面. 图1给出了该问题的计算模型. 初始时刻, 液滴为圆形, 直径
${d_0} = 0.02\;{{\rm{m}}}$ . 液滴到固壁的高度为0.03 m. 液滴具有初始下落速度$V = 1\;{{\rm{m}}} \cdot {{{\rm{s}}} ^{ - 1}}$ . 液滴在重力作用下加速下落, 并撞击固壁发生铺展变形. 液滴密度$\rho = 1000\;{{\rm{kg}}} \cdot {{{\rm{m}}} ^{ - 3}}$ , 取向松弛时间${\lambda _1} = 0.02\;{{\rm{s}}}$ . 液滴总黏度$\eta = 4\;{{\rm{Pa}}} \cdot {{{\rm{s}}} ^{ - 1}}$ , 其中溶剂黏度${\eta _s} = 0.4\;{\text{Pa}} \cdot {{\text{s}}^{ - 1}}$ 和聚合物黏度${\eta _{{p}}} = 3.6\;{{\rm{Pa}}} \cdot {{\rm{s}}^{ - 1}}$ . 液滴温度${T_{\rm{f}}} = 443\;{\text{K}}$ , 固壁温度${T_{\rm{w}}} = 313\;{\text{K}}$ . 比热容${c_p} = 1\;{\text{J}} \cdot {\text{k}}{{\text{g}}^{ - 1}} \cdot {{\text{K}}^{ - 1}}$ , 导热系数$\kappa = 1\;{\text{W}} \cdot {{\text{m}}^{ - 1}} \cdot {{\text{K}}^{ - 1}}$ . 温度依赖系数$\varphi = 0.002\;{{\text{K}}^{ - 1}}$ . 分别选取液滴直径${d_0}$ 和下降速度$V$ 作为特征长度和特征速度, 则Péclet数Pe = 20, Reynolds数Re = 5, Weissenberg 数Wi = 1, 溶剂黏度比β = 0.1. 其他XPP流变参数选取如下: α = 0.01, γ = 0.9, Q = 4. 粒子初始间距设置为Δx = 0.0002, 对应于7845个流体粒子, 301个边界粒子和903个虚拟粒子. 核函数采用5次样条函数, 光滑长度$h = 1.5\Delta x$ . 声速$c = 12.5\;{{\rm{m}}} \cdot {{{\rm{s}}} ^{ - 1}}$ , 时间步长$\Delta t = 1.0 \times {10^{ - 6}}\;{{\rm{s}}}$ .![]()
图 1 液滴冲击固体表面的计算模型Figure 1. Computational model of droplet impinging solid wall surface图2给出液滴在撞击固壁之后6个时刻
${t^*}$ = 1.5, 1.85, 2.35, 2.55, 3.85, 4.85的SPH模拟结果, 其中${t^*} = tV/{d_0}$ 表示液滴铺展的无量纲时间. 可以看出, 液滴在撞击固壁后的流动过程分为3个阶段. (1) 快速铺展阶段(${t^*}$ = 1.5, 1.85, 2.35). 由于液滴具有初始向下的速度且受重力作用的影响, 因此液滴在该阶段的快速铺展主要受惯性力所驱动. 在该阶段, 液滴左半部分流速向左, 右半部分流速向右; 宽度增加, 而高度下降. 同时, 由于液滴底部区域最先接触冷固壁, 故液滴底部温度首先降低. (2) 收缩阶段(${t^*}$ = 2.55, 3.85). 液滴在铺展到最大宽度后, 受弹性力的作用而逐渐收缩. 在该阶段, 液滴的速度方向发生了转变, 即左半部分流速向右, 右半部分流速向左. 同时, 由于液滴向冷固壁的传热, 液滴底部区域的温度已冷却到和冷固壁相同. (3) 再铺展阶段(${t^*}$ = 4.85). 液滴在收缩到最小宽度后, 由于能量的耗散和受重力的作用, 再次缓慢铺展. 在该阶段, 液滴持续向冷固壁传热, 导致整个液滴温度的下降, 并逐渐冷却到冷固壁的温度.![]()
图 2 运用改进SPH方法模拟得到的液滴冲击固体表面在不同时刻的温度分布Figure 2. The temperature distribution of droplet impinging solid wall surface at different times obtained by improved SPH method图3给出液滴的铺展宽度随时间的变化情况. 为了便于分析, 这里也用液滴直径d0对液滴铺展宽度dT 进行无量纲处理. 从图中可以清晰地观察到液滴撞击固壁之后的3个流动阶段, 其中在
${t^*}$ ≈ 2.250时铺展到最大宽度(约2.225), 在${t^*}$ ≈ 3.955时收缩到最小宽度(约1.775). 另外, 为了验证本文改进SPH方法模拟非等温黏弹性复杂流动的有效性, 本文特别利用了Basilisk流体仿真软件[41]对该问题进行了数值模拟. 该软件基于有限体积法求解控制方程, 而自由界面的追踪通过流体体积法来实现. 图3给出了利用本文改进SPH方法和Basilisk软件得到的液滴铺展宽度随时间变化情况的比较. 不难看出, 利用两种方法模拟得到的结果是基本吻合的, 均经历了快速铺展、收缩和再铺展的阶段, 从而表明了本文改进SPH方法模拟非等温黏弹性复杂流动是准确有效的. 同时, 也需注意, 两种数值结果在液滴的收缩阶段有一定的差异, 可能归因于: Basilisk软件使用Poisson方程求解压力, 而本文考虑弱可压缩流体, 使用状态方程求解压力.![]()
图 3 利用改进SPH方法和Basilisk软件得到的液滴铺展宽度随时间变化的比较Figure 3. Comparison of the time changes of the droplet spread width obtained by the improved SPH method and Basilisk software液滴撞击固壁问题含有复杂的自由面, 而对于含自由面的非等温黏弹性流体的数值模拟, 若直接使用SPH方法进行模拟, 极易产生拉伸不稳定性, 并最终导致计算的崩溃(见图4(a)). 利用本文粒子迁移技术可有效地解决该问题, 并得到稳定准确的结果(见图4(b)), 这表明了SPH方法中施加粒子迁移技术的有效性. 此外, 液滴不仅具有初始下落的速度, 而且受重力的作用加速下落, 因此液滴会对固壁边界产生较大的冲击力. 然而, 在模拟过程中, 并未观察到任何粒子穿透固壁边界现象的发生, 表明本文发展的边界粒子和虚拟粒子相联合的边界处理方法能够有效地处理固壁边界. 总之, 相较于网格类方法, 本文为非等温黏弹性流体的数值模拟提供了一种新型的、稳定精确的数值方法.
![]()
图 4 利用不同方法得到的液滴在t* = 2.5时刻的模拟结果Figure 4. The simulation result of droplet at t* = 2.5 obtained by different methods为了观察温度的引入对液滴的流动物理是否有影响, 图5给出了等温液滴和非等温液滴铺展宽度的对比, 其中等温液滴在程序中只需设置温度依赖系数φ = 0. 可以看出, 在
${t^*}$ < 2.8时刻, 等温液滴和非等温液滴的铺展宽度变化几乎相同. 这表明: 在流动的早期阶段, 温度的引入并未对液滴的流动物理产生明显的影响. 然而, 随着时间的推移, 非等温液滴表现出了较等温液滴更强的收缩行为. 特别是, 非等温液滴的最小收缩宽度(约1.775)明显低于等温液滴的最小收缩宽度(约1.825). 在等温液滴和非等温液滴均收缩到最小宽度后, 它们均再次铺展. 然而, 非等温液滴的再铺展宽度略低于等温液滴的再铺展宽度.![]()
图 5 温度依赖系数φ对液滴铺展宽度随时间变化的影响Figure 5. Influence of the temperature dependency coefficient φ on the time change of the droplet spread width为了研究温度依赖系数φ对液滴流动物理的影响, 图5给出了3个不同φ (φ = 0.002, 0.005, 0.01)取值下液滴的铺展宽度随时间变化的比较. 可以看出, φ越大, 对液滴流动过程的影响越大. 虽然在早期阶段, 不同φ取值下的液滴铺展宽度几乎相同; 但在
${t^*}$ > 2.0, 三者的不同之处愈加明显. 特别是, 液滴的最大铺展宽度随着φ的增大而减小. φ = 0.010的液滴最大铺展宽度约为2.105, 明显低于φ = 0.005, 0.002时的液滴最大铺展宽度(约2.160和2.225). 另外, φ越大, 液滴的最小收缩宽度也越小. 对于φ = 0.002, 0.005, 0.010, 其最小收缩宽度分别约为1.775, 1.655, 1.575. 由于液滴较小的收缩宽度, 液滴再铺展阶段的铺展宽度也因此变小. 以上现象可解释为: 在流动的早期阶段, 虽然有液滴向冷固壁的传热, 但该阶段的传热并不极大地改变液滴的物理属性. 早期阶段的液滴主要受惯性力驱动而快速铺展因此, 即使有液滴向冷固壁的传热, 温度对液滴流动物理的影响是可以忽略的. 然而, 随着时间的推移, 液滴底部温度逐渐冷却到冷固壁温度. 在这种情况下, 根据时温等效原理式(11)和式(12), 液滴底部区域的“有效”黏度逐渐增大, 并最终影响到液滴的流动物理. 因此, φ越大, 液滴的最大铺展宽度越小. 另外, 在流动后期, 仍有液滴向冷固壁的持续传热, 液滴最终冷却到和冷固壁温度相同. 因此, φ越大, 液滴的最小收缩宽度也越小. 在${t^*}$ > 7.0后, φ = 0.01的液滴铺展较其他两种情况缓慢, 因此再铺展宽度相应地变小.下面讨论其他各流变参数对液滴铺展宽度随时间的变化的影响.
(1) Re的影响
Re是描述流体流动的一个重要参数, 它比较了惯性力和黏性力. Re对液滴铺展宽度的影响如图6(a)所示. 可以看出, Re越大, 液滴的最大铺展宽度越大, 达到最大铺展宽度所需的时间相对滞后. 对于Re = 3, 5, 7, 液滴的最大铺展宽度分别约为1.935, 2.225, 2.395, 达到最大铺展宽度所需的时间分别约为2.105, 2.250和2.415. 这是因为: Re越大, 液滴惯性力越大, 液滴撞击固壁后有更快的铺展行为. 另外, 随着Re的增大, 液滴的最小收缩宽度相应地增加. 不同Re得到的液滴最小收缩宽度分别约为 1.475, 1.775和 2.015. Re越大, 液滴的收缩行为越弱. 这是因为, Re越大, 液滴弹性力相对越小.
![]()
图 6 不同流变参数下液滴铺展宽度随时间的变化Figure 6. Time change of the droplet spread width under different rheological parameters(2) Wi的影响
Wi是衡量黏弹性流体流动的一个重要参数, 它比较了弹性力和黏性力. 图6(b)描绘了Wi对液滴铺展宽度的影响. 可以看出, 随着Wi 的增加, 液滴达到的最大铺展宽度显著增加, 达到最大铺展宽度所需的时间相对滞后. 对于 Wi = 0.5, 1, 3, 5, 液滴达到的最大铺展宽度分别约为1.985, 2.225, 2.545和2.720, 达到最大铺展宽度的时间分别约为 2.150, 2.250, 2.600和2.755. 另外, 随着Wi的增大, 液滴的最小收缩宽度也相应增加. 不同 Wi 得到的液滴最小收缩宽度分别约为 1.745, 1.775, 1.955和 2.120. 这是因为, XPP模型可合理地描述黏弹性流体的剪切稀化行为; 随着Wi 的增大, 黏弹性流体的剪切稀化增强, 导致液滴的铺展更快, 因此液滴的铺展宽度相应增加.
(3) β的影响
β表示黏弹性流动中牛顿溶剂黏度与流体总黏度的比值, 它对液滴铺展宽度的影响如图6(c)所示. 可以看出, β越大, 液滴的收缩行为越弱, 但最终铺展宽度相差不大. 这是因为, 随着β 的增大, 牛顿溶剂黏度增加, 聚合物黏度降低, 导致液滴的弹性力变小. 在这种情形下, 液滴表现出类似于牛顿流体的铺展行为. 另外, 改变β只改变牛顿溶剂黏度占总黏度的比值, 而不改变总黏度的大小, 因此液滴的最终铺展宽度相差不大.
(4) α的影响
α是黏弹性流动中控制各向异性拖曳力相关的流变参数, 它对液滴铺展宽度的影响如图6(d)所示. 可以看出, 改变α 并不影响液滴达到的最大铺展宽度, 但在一定程度上影响液滴的收缩行为, 即α越大, 液滴的收缩行为越弱, 液滴的最小收缩宽度越大. 对于 α = 0, 0.01, 0.1, 0.3, 液滴的最小收缩宽度分别约为1.770, 1.775, 1.800和 1.825. 由于液滴收缩行为变弱, 液滴的最后铺展宽度随着α的增大而增大.
(5) γ的影响
γ表示聚合物分子链的拉伸松弛时间和取向松弛时间的比值. 图6(e)给出了γ对液滴铺展宽度的影响. γ越小, 液滴的最大铺展宽度越大, 达到最大铺展宽度所需的时间也相对滞后. 对于γ = 0.9, 0.7, 0.5, 0.3, 液滴的最大铺展宽度分别约为 2.225, 2.235, 2.240 和 2.300, 达到最大铺展宽度的时间分别约为 2.250, 2.265, 2.300和 2.485. 另外, γ在一定程度上也影响着液滴的收缩行为, 即γ越小, 液滴的最小收缩宽度越大, 这是因为: 降低γ导致液滴的拉伸松弛时间变短, 聚合物分子链得不到有效的拉伸, 因此液滴的弹性力变小. 由于变弱的液滴收缩行为, 液滴的最后铺展宽度随着γ的减小而增大.
(6) Q的影响
分子链臂数Q对液滴铺展宽度的影响如图6(f)所示. 可以看出, Q越大, 对流动行为的影响越弱. 另外, Q越小, 液滴的收缩行为越弱, 达到的最终铺展宽度越大. 对于Q = 8, 4, 2, 1, 液滴的最小收缩宽度分别约为1.755, 1.775, 1.805和 1.875.
3.2 F型腔注塑成型
本文选取的第2个算例为基于XPP本构模型的非等温黏弹性F型腔注塑成型问题, 该问题也含有复杂的自由面. 图7给出了该问题的计算模型. F型腔的右侧有两个注射口, 宽度均为
$L = 0.1\;{{\rm{m}}}$ . 熔体注射速度$V = 1\;{{\rm{m}}} \cdot {{{\rm{s}}} ^{ - 1}}$ . 熔体密度$\rho = 1000\;{{\rm{kg}}} \cdot {{{\rm{m}}} ^{ - 3}}$ , 取向松弛时间${\lambda _1} = 0.1\;{{\rm{s}}}$ . 熔体总黏度$\eta = 20\;{{\rm{Pa}}} \cdot {{{\rm{s}}} ^{ - 1}}$ , 其中溶剂黏度${\eta _{{s}}} = 10\;{\text{Pa}} \cdot {{\text{s}}^{ - 1}}$ 和聚合物黏度${\eta _{{p}}} = 10\;{{\rm{Pa}}} \cdot {{\rm{s}}^{ - 1}}$ . 熔体初始温度${T_{\rm{f}}} = 443\;{\text{K}}$ , 模壁温度${T_{\rm{w}}} =$ 313 K. 比热容${c_p} = 10\;{\text{J}} \cdot {\text{k}}{{\text{g}}^{ - 1}} \cdot {{\text{K}}^{ - 1}}$ , 导热系数$\kappa = 10\;{\text{W}} \cdot {{\text{m}}^{ - 1}} \cdot {{\text{K}}^{ - 1}}$ , 温度依赖系数$\varphi = 0.001\;{{\text{K}}^{ - 1}}$ . 分别选取注射口宽度L和注射速度$V$ 作为特征长度和特征速度, 则无量纲Péclet数Pe = 100, Reynolds数Re = 5, Weissenberg 数Wi = 1, 溶剂黏度比β = 0.5. 其他XPP流变参数选取如下: α = 0.01, γ = 0.9, Q = 4. 粒子初始间距设置为Δx = 0.002, 对应于1102个边界粒子和3360个虚拟粒子, 同时使用了27386个流体粒子充满整个型腔. 核函数采用5次样条函数, 光滑长度$h = 1.5\Delta x$ . 声速$c = 10\;{{\rm{m}}} \cdot {{{\rm{s}}} ^{ - 1}}$ , 时间步长$\Delta t = 1.0 \times {10^{ - 5}}\;{{\rm{s}}}$ . 对于该算例, 本文在联想计算机(Intel(R) Core(TM) i7-10700 CPU @ 3.60 GHz)模拟7万个时间步, 所耗时间约为4.6 h.![]()
图 7 F型腔注塑成型问题的计算模型Figure 7. Computational model of injection molding process for an F-shaped cavity图8给出了运用改进SPH方法模拟得到的F型腔注塑成型在4个不同时刻t = 0.32 s, 0.44 s, 0.56 s, 0.68 s的温度分布. 在早期阶段, 两注射口的熔体沿各自管道由右向左流动. 熔体温度沿水平管道上下对称. 此时, 已发生熔体向冷固壁的传热, 因此靠近管道上下模壁处的熔体温度逐渐降低. 由于下管道较短, 因此下注射口的熔体率先到达左侧垂直固壁(t = 0.32 s), 并沿该固壁分成上下两股熔体. 其中, 上股熔体向上流动, 而下股熔体向下流动. 在此阶段, 同样伴有熔体向左侧模壁的传热, 因此靠近左侧模壁处的熔体温度随之降低. 随着时间的推移, 上注射口的熔体也很快到达左侧垂直固壁(t = 0.44 s), 并沿左上270°转角向下流动进入垂直管道. 靠近转角处熔体的温度明显低于中心区域处熔体的温度. 不久之后, 上注射口的向下流动熔体和下注射口流入垂直管道的上股熔体相遇(t = 0.56 s), 形成熔接痕. 随着熔体的持续注入, F型腔的上半部分被完全充满(t = 0.68 s). 两注射口的熔体汇成一股熔体并沿垂直管道向下流动; 该股熔体最终充满整个腔体. 对于该问题的SPH模拟结果, 并未见到流体的拉伸不稳定性和粒子穿透固壁边界现象的发生, 这进一步体现了本文提出的改进SPH方法模拟非等温黏弹性复杂流动的优势.
![]()
图 8 运用改进SPH方法模拟得到的F型腔注塑成型在不同时刻的温度分布Figure 8. The temperature distribution of F-shaped cavity injection molding at different times obtained by improved SPH method为了研究Péclet数Pe对成型过程中传热特性的影响, 图9给出了利用4组不同Pe (Pe = 50, 100, 200, 400)得到的F型腔注塑成型在t = 0.6 s时刻的温度分布的比较. 可以看出, Péclet数越大, 腔体的热边界层越薄. 该现象可解释为: Péclet数表示流动中对流速率与扩散速率之比. Péclet数越大, 比热容越大, 熔体热扩散越慢. 在这种情况下, 熔体会携带更多的热量进入腔体, 因此腔体的热边界层越薄.
![]()
图 9 Péclet数对F型腔注塑成型在t = 0.6 s时刻的温度分布的影响Figure 9. Effect of Péclet number on the temperature distribution at t = 0.6 s in F-shaped cavity injection molding对于注塑成型, 应力对高聚物制品的质量有决定性的影响. 因此, 图10进一步给出了运用改进SPH方法模拟得到的F型腔注塑成型在4个不同时刻t = 0.32 s, 0.44 s, 0.56 s, 0.68 s的第一法向应力差分布. 可以看出, 在早期阶段, 熔体第一法向应力差沿水平管道上下对称. 越靠近管道中心区域, 第一法向应力差越小; 越靠近管道上下模壁, 第一法向应力差越大. 在t = 0.44 s时刻, 下注射口熔体流入左侧垂直管道. 垂直管道内熔体的第一法向应力差由正值逐渐变为负值(见垂直管道内熔体的中心区域处). 在t = 0.56 s时刻, 上注射口熔体沿左上转角流入垂直管道. 该股熔体在垂直管道内的第一法向应力差也由正值逐渐变为负值. 在t = 0.68 s时刻, 腔体上半部分被完全充满, 但可清晰地观察到熔接痕处的第一法向应力差较大. 由于较大的第一法向应力差, 高聚物制品在该位置处的强度较低. 此外, 腔体右下90°转角处的第一法向应力差分布较为复杂, 其中转角上模壁附近熔体的第一法向应力差为正值, 左模壁附近熔体的第一法向应力差为负值.
![]()
图 10 运用改进SPH方法模拟得到的F型腔注塑成型在不同时刻的第一法向应力差分布Figure 10. The first normal stress difference distribution of F-shaped cavity injection molding at different times obtained by improved SPH method对于该问题, 由于熔体流变特性(剪切变稀)的复杂性和多物理场(温度、应力、压力等)耦合的复杂性, 目前并未见到基于网格的方法对其进行数值模拟的文献报道, 因此本文通过评价SPH数值解的收敛性来验证SPH方法的有效性. 特别增加了粒子初始间距分别为Δx = 0.004和0.00125的数值试验, 而其它参数均保持不变. 图11给出了利用这2组不同粒子初始间距和原始粒子初始间距(Δx = 0.002)得到的点A处第一法向应力差随时间变化的比较, 其中点A到腔体右下90°转角的高度为0.01 m(如图7所示). 可以看出, 在t < 0.22 s, 点A处的第一法向应力差一直为0, 表明熔体在0.22 s之前还未到达该点位置. 在0.22 s < t < 0.35 s, 随着熔体流经点A, 第一法向应力差逐渐增大, 并于t ~ 0.35 s达到最大值. 在0.35 s < t < 0.60 s, 达到最大值的第一法向应力差呈小幅震荡的水平状态. 在t > 0.60 s, 第一法向应力差略有轻微下降. 通过3组不同粒子初始间距得到的第一法向应力差随时间变化的结果基本一致, 表明了本文改进SPH方法模拟非等温黏弹性注塑成型问题是收敛的. 此外, 为了定量地评价数值收敛性, 进一步引入
${l_2}$ 范数误差![]()
图 11 利用不同粒子间距Δx得到的点A处第一法 向应力差N1随时间变化的比较Figure 11. Comparison of the time change of the first normal stress difference N1 at point A obtained by different particle spacings Δx$$ {{E}} = \sqrt {\dfrac{{\displaystyle\sum {{{\left( {{R_{\Delta x = 0.001\;25}} - {R_{{\rm{SPH}}}}} \right)}^2}} }}{{\displaystyle\sum {{{\left( {{R_{\Delta x = 0.001\;25}}} \right)}^2}} }}} $$ (37) 其中, 以Δx = 0.00125数值解作为基准解. 经计算, 利用Δx = 0.004, 0.002得到的SPH解相较于基准解的
${l_2}$ 范数误差分别为2.69%和0.95%. 误差随着粒子的加密而下降, 进一步表明了本文改进SPH方法模拟非等温黏弹性注塑成型问题是收敛的.下面讨论其它各流变参数对成型过程中应力随时间变化的影响.
(1) Re的影响
图12(a)展示了Re对点A处第一法向应力差N1随时间变化的影响. 可以看出, 随着Re的增大, 第一法向应力差逐渐减小. 这是因为: Re越大, 流体总黏度越低, 其中的聚合物黏度也越低, 因此计算得到的第一法向应力差越小.
![]()
图 12 不同流变参数下点A处的第一法向应力差N1随时间的变化Figure 12. Time change of the first normal stress difference N1 at point A under different rheological parameters(2) Wi的影响
图12(b)描述了Wi对点A处第一法向应力差N1随时间变化的影响. Wi越大, 第一法向应力差越小. 这是因为: XPP模型可合理地描述黏弹性流体的剪切变稀行为; 当Wi较大时, 黏弹性流体发生了明显的剪切变稀, 因此计算得到的第一法向应力差变小.
(3) β的影响
图12(c)给出了β对点A处第一法向应力差N1随时间变化的影响. β越大, 第一法向应力差越小. 这是因为: 随着β的增大, 牛顿溶剂黏度逐渐增加. 由于流体总黏度不变, 导致其中的聚合物黏度降低, 因此计算得到的第一法向应力差变小.
(4) α的影响
α对点A处第一法向应力差N1随时间变化的影响如图12(d)所示. 可以看出, 利用不同α得到的第一法向应力差均几乎相同, 这表明α并非影响成型过程中应力的主要因素.
(5) γ的影响
图12(e)描述了γ对点A处第一法向应力差N1随时间变化的影响. γ越大, 第一法向应力差越大. 这是因为: γ越大, 聚合物熔体中分子链的拉伸松弛时间越长, 熔体需要较长的时间来达到分子链拉伸的稳态, 因此计算得到的第一法向应力差较大.
(6) Q的影响
Q表示聚合物分子链的臂数, 它对点A处第一法向应力差N1随时间变化的影响如图12(f)所示. 可以看出, Q越大, 第一法向应力差越大.
4. 结论
本文提出了一种改进的SPH方法对基于XPP模型的非等温黏弹性复杂流动进行了数值模拟. 为了提高模拟结果的精度, 采用了一种核函数梯度的修正算法; 为了灵活地施加边界条件, 发展了边界粒子和虚拟粒子相联合的边界处理方法; 为了消除流动过程中的拉伸不稳定性, 施加了粒子迁移技术. 运用改进SPH方法数值模拟了液滴撞击固壁和F型腔注塑成型问题, 本文所得结论如下.
(1) 对于非等温黏弹性复杂流动, 利用本文提出的改进SPH方法可有效地解决流动过程中出现的拉伸不稳定性, 且能够有效地防止粒子穿透固壁边界, 因此本文为非等温黏弹性流体的模拟提供了一种新型的、稳定准确的数值方法.
(2) 对于非等温黏弹性液滴撞击固壁, 通过比较SPH结果和Basilisk流体仿真软件得到的结果, 验证了改进SPH方法模拟非等温黏弹性复杂流动的有效性. 此外, 温度的引入会导致液滴表现出更强的收缩行为. Re和Wi越大, 液滴铺展宽度越大. β和α越大, 液滴收缩行为越弱. γ和Q越大, 液滴收缩行为越强.
(3) 对于非等温黏弹性F型腔注塑成型, 通过利用不同粒子初始间距进行计算, 表明了改进SPH方法模拟非等温黏弹性复杂流动是数值收敛的. 此外, Pe越大, 腔体的热边界层越薄; Re, Wi和β 越大, 腔体的第一法向应力差越小; γ和Q越大, 第一法向应力差越大; α对第一法向应力差的影响可以忽略.
(4) 本文提出的改进SPH方法能够稳定、准确地描述非等温黏弹性流体的传热机理、复杂流变特性和自由面变化特性.
-
![]()
图 1 液滴冲击固体表面的计算模型
Figure 1. Computational model of droplet impinging solid wall surface
![]()
图 2 运用改进SPH方法模拟得到的液滴冲击固体表面在不同时刻的温度分布
Figure 2. The temperature distribution of droplet impinging solid wall surface at different times obtained by improved SPH method
![]()
图 3 利用改进SPH方法和Basilisk软件得到的液滴铺展宽度随时间变化的比较
Figure 3. Comparison of the time changes of the droplet spread width obtained by the improved SPH method and Basilisk software
![]()
图 4 利用不同方法得到的液滴在t* = 2.5时刻的模拟结果
Figure 4. The simulation result of droplet at t* = 2.5 obtained by different methods
![]()
图 5 温度依赖系数φ对液滴铺展宽度随时间变化的影响
Figure 5. Influence of the temperature dependency coefficient φ on the time change of the droplet spread width
![]()
图 6 不同流变参数下液滴铺展宽度随时间的变化
Figure 6. Time change of the droplet spread width under different rheological parameters
![]()
图 7 F型腔注塑成型问题的计算模型
Figure 7. Computational model of injection molding process for an F-shaped cavity
![]()
图 8 运用改进SPH方法模拟得到的F型腔注塑成型在不同时刻的温度分布
Figure 8. The temperature distribution of F-shaped cavity injection molding at different times obtained by improved SPH method
![]()
图 9 Péclet数对F型腔注塑成型在t = 0.6 s时刻的温度分布的影响
Figure 9. Effect of Péclet number on the temperature distribution at t = 0.6 s in F-shaped cavity injection molding
![]()
图 10 运用改进SPH方法模拟得到的F型腔注塑成型在不同时刻的第一法向应力差分布
Figure 10. The first normal stress difference distribution of F-shaped cavity injection molding at different times obtained by improved SPH method
![]()
图 11 利用不同粒子间距Δx得到的点A处第一法 向应力差N1随时间变化的比较
Figure 11. Comparison of the time change of the first normal stress difference N1 at point A obtained by different particle spacings Δx
-
[1] Peters GWM, Baaijens FPT. Modelling of non-isothermal viscoelastic flows. Journal of Non-Newtonian Fluid Mechanics, 1997, 68(2-3): 205-224 doi: 10.1016/S0377-0257(96)01511-X
[2] Li Q, Qu F. A level set based immersed boundary method for simulation of non-isothermal viscoelastic melt filling process. Chinese Journal of Chemical Engineering, 2021, 32: 119-133 doi: 10.1016/j.cjche.2020.09.057
[3] Fernandes C. A fully implicit log-conformation tensor coupled algorithm for the solution of incompressible non-isothermal viscoelastic flows. Polymers, 2022, 14(19): 4099 doi: 10.3390/polym14194099
[4] Moreno L, Codina R, Baiges J. Numerical simulation of non-isothermal viscoelastic fluid flows using a VMS stabilized finite element formulation. Journal of Non-Newtonian Fluid Mechanics, 2021, 296: 104640 doi: 10.1016/j.jnnfm.2021.104640
[5] Gingold RA, Monaghan JJ. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society, 1977, 181(3): 375-389 doi: 10.1093/mnras/181.3.375
[6] Lucy LB. A numerical approach to the testing of the fission hypothesis. The Astronomical Journal, 1977, 82: 1013-1024 doi: 10.1086/112164
[7] 王璐, 徐绯, 杨扬. 完全拉格朗日 SPH 在冲击问题中的改进和应用. 力学学报, 2022, 54(12): 3297-3309 (Wang Lu, Xu Fei, Yang Yang. Improvement of the total Lagrangian SPH and its application in impact problems. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(12): 3297-3309 (in Chinese) [8] Liu MB, Zhang ZL, Feng DL. A density-adaptive SPH method with kernel gradient correction for modeling explosive welding. Computational Mechanics, 2017, 60(3): 513-529 doi: 10.1007/s00466-017-1420-5
[9] 王平平, 张阿漫, 彭玉祥等. 近场水下爆炸瞬态强非线性流固耦合无网格数值模拟研究. 力学学报, 2022, 54(8): 2194-2209 (Wang Pingping, Zhang Aman, Peng Yuxiang, et al. Numerical simulation of transient strongly-nonlinear fluid-structure interaction in near-field underwater explosion based on meshless method. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(8): 2194-2209 (in Chinese) [10] 姚学昊, 陈丁, 武立伟等. 流固耦合破坏分析的多分辨率 PD-SPH 方法. 力学学报, 2022, 54(12): 3333-3343 (Yao Xuehao, Chen Ding, Wu Liwei, et al. A multi-resolution PD-SPH coupling approach for structural failure under fluid-structure interaction. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(12): 3333-3343 (in Chinese) [11] He F, Zhang H, Huang C, et al. A stable SPH model with large CFL numbers for multi-phase flows with large density ratios. Journal of Computational Physics, 2022, 453: 110944 doi: 10.1016/j.jcp.2022.110944
[12] Liu MB, Liu GR. Restoring particle consistency in smoothed particle hydrodynamics. Applied Numerical Mathematics, 2006, 56(1): 19-36 doi: 10.1016/j.apnum.2005.02.012
[13] Xue B, Wang SP, Peng YX, et al. A novel coupled Riemann SPH–RKPM model for the simulation of weakly compressible fluid–structure interaction problems. Ocean Engineering, 2022, 266: 112447 doi: 10.1016/j.oceaneng.2022.112447
[14] Ng KC, Alexiadis A, Ng YL. An improved particle method for simulating fluid-structure interactions: The multi-resolution SPH-VCPM approach. Ocean Engineering, 2022, 247: 110779 doi: 10.1016/j.oceaneng.2022.110779
[15] Monaghan JJ, Kajtar JB. SPH particle boundary forces for arbitrary boundaries. Computer Physics Communications, 2009, 180(10): 1811-1820 doi: 10.1016/j.cpc.2009.05.008
[16] Morris JP, Fox PJ, Zhu Y. Modeling low Reynolds number incompressible flows using SPH. Journal of Computational Physics, 1997, 136(1): 214-226 doi: 10.1006/jcph.1997.5776
[17] Liu MB, Shao JR, Chang JZ. On the treatment of solid boundary in smoothed particle hydrodynamics. Science China Technological Sciences, 2012, 55(1): 244-254 doi: 10.1007/s11431-011-4663-y
[18] Yildiz M, Rook RA, Suleman A. SPH with the multiple boundary tangent method. International Journal for Numerical Methods in Engineering, 2009, 77(10): 1416-1438 doi: 10.1002/nme.2458
[19] Chen C, Zhang AM, Chen JQ, et al. SPH simulations of water entry problems using an improved boundary treatment. Ocean Engineering, 2021, 238: 109679 doi: 10.1016/j.oceaneng.2021.109679
[20] Yang X, Liu M, Peng S. Smoothed particle hydrodynamics modeling of viscous liquid drop without tensile instability. Computers & Fluids, 2014, 92: 199-208
[21] Sun PN, Colagrossi A, Marrone S, et al. Multi-resolution Delta-plus-SPH with tensile instability control: Towards high Reynolds number flows. Computer Physics Communications, 2018, 224: 63-80 doi: 10.1016/j.cpc.2017.11.016
[22] Chalk CM, Pastor M, Peakall J, et al. Stress-Particle smoothed particle hydrodynamics: An application to the failure and post-failure behaviour of slopes. Computer Methods in Applied Mechanics and Engineering, 2020, 366: 113034 doi: 10.1016/j.cma.2020.113034
[23] Lyu HG, Sun PN. Further enhancement of the particle shifting technique: Towards better volume conservation and particle distribution in SPH simulations of violent free-surface flows. Applied Mathematical Modelling, 2022, 101: 214-238 doi: 10.1016/j.apm.2021.08.014
[24] Rafiee A, Manzari MT, Hosseini M. An incompressible SPH method for simulation of unsteady viscoelastic free-surface flows. International Journal of Non-Linear Mechanics, 2007, 42(10): 1210-1223 doi: 10.1016/j.ijnonlinmec.2007.09.006
[25] 杨波, 欧阳洁, 蒋涛等. PTT 黏弹性流体的光滑粒子动力学方法模拟. 力学学报, 2011, 43(4): 667-673 (Yang Bo, Ouyang Jie, Jiang Tao, et al. Numerical simulation of the viscoelastic flows for PTT model by the SPH method. Chinese Journal of Theoretical and Applied Mechanics, 2011, 43(4): 667-673 (in Chinese) [26] Xu X, Deng XL. An improved weakly compressible SPH method for simulating free surface flows of viscous and viscoelastic fluids. Computer Physics Communications, 2016, 201: 43-62 doi: 10.1016/j.cpc.2015.12.016
[27] King JRC, Lind SJ. High weissenberg number simulations with incompressible smoothed particle hydrodynamics and the log-conformation formulation. Journal of Non-Newtonian Fluid Mechanics, 2021, 293: 104556 doi: 10.1016/j.jnnfm.2021.104556
[28] Duque-Daza C, Alexiadis A. A simplified framework for modelling viscoelastic fluids in discrete multiphysics. Chem Engineering, 2021, 5(3): 61
[29] Vahabi M, Hadavandmirzaei H, Kamkari B, et al. Interaction of a pair of in-line bubbles ascending in an Oldroyd-B liquid: A numerical study. European Journal of Mechanics-B/Fluids, 2021, 85: 413-429 doi: 10.1016/j.euromechflu.2020.11.004
[30] Moinfar Z, Vahabi S, Vahabi M. Numerical simulation of drop deformation under simple shear flow of Giesekus fluids by SPH. International Journal of Numerical Methods for Heat & Fluid Flow, 2023, 33(1): 263-281
[31] Meburger S, Niethammer M, Bothe D, et al. Numerical simulation of non-isothermal viscoelastic flows at high weissenberg numbers using a finite volume method on general unstructured meshes. Journal of Non-Newtonian Fluid Mechanics, 2021, 287: 104451 doi: 10.1016/j.jnnfm.2020.104451
[32] Verbeeten WMH, Peters GWM, Baaijens FPT. Differential constitutive equations for polymer melts: The extended Pom–Pom model. Journal of Rheology, 2001, 45(4): 823-843 doi: 10.1122/1.1380426
[33] Tanner RI. Engineering Rheology. OUP Oxford, 2000
[34] O'connor J, Domínguez JM, Rogers BD, et al. Eulerian incompressible smoothed particle hydrodynamics on multiple GPUs. Computer Physics Communications, 2022, 273: 108263 doi: 10.1016/j.cpc.2021.108263
[35] Liu GR, Liu MB. Smoothed Particle Hydrodynamics: A Meshfree Particle Method. Singapore: World Scientific, 2003
[36] Shao S, Lo EYM. Incompressible SPH method for simulating newtonian and non-newtonian flows with a free surface. Advances in Water Resources, 2003, 26(7): 787-800 doi: 10.1016/S0309-1708(03)00030-7
[37] Bonet J, Lok TSL. Variational and momentum preservation aspects of smooth particle hydrodynamic formulations. Computer Methods in Applied Mechanics and Engineering, 1999, 180(1-2): 97-115 doi: 10.1016/S0045-7825(99)00051-1
[38] Lind SJ, Xu R, Stansby PK, et al. Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves. Journal of Computational Physics, 2012, 231(4): 1499-1523 doi: 10.1016/j.jcp.2011.10.027
[39] Wang PP, Meng ZF, Zhang AM, et al. Improved particle shifting technology and optimized free-surface detection method for free-surface flows in smoothed particle hydrodynamics. Computer Methods in Applied Mechanics and Engineering, 2019, 357: 112580 doi: 10.1016/j.cma.2019.112580
[40] Yang L, Rakhsha M, Hu W, et al. A consistent multiphase flow model with a generalized particle shifting scheme resolved via incompressible SPH. Journal of Computational Physics, 2022, 458: 111079 doi: 10.1016/j.jcp.2022.111079
[41] Popinet S. An accurate adaptive solver for surface-tension-driven interfacial flows. Journal of Computational Physics, 2009, 228(16): 5838-5866 doi: 10.1016/j.jcp.2009.04.042
下载:
计量
- 文章访问数: 734
- HTML全文浏览量: 254
- PDF下载量: 136
