基于分子动力学方法模拟聚合物稀溶液喷泉流动

Simulation of Fountain Flow Based on Molecular Dynamics Method

摘要: 本文基于分子动力学原理利用FENE哑铃模型对聚合物稀溶液流动前沿喷泉效应开展数值模拟。利用欧拉法求解力学模型在简单剪切流场下的构型方程和哑铃分子的位形方程，获取示踪流线和流动前沿哑铃分布，计算应力场，分析流变演化规律，研究温度和剪切速率等参数对该模型的影响，分析FENE哑铃分子模型模拟喷泉效应的有效性。结果表明，在喷泉效应的作用下，随着剪切速率的升高哑铃分子的伸展变大，流动前沿聚合物应力会出现应力过冲，随后趋于稳定的复杂变化，哑铃分子沿喷泉流线取向。

Abstract: In this paper, numerical simulations of the fountain effect on the flow front of polymer dilute solu-tions are carried out based on molecular dynamics principles and the FENE dumbbell model. The Euler method is used to solve the constitutive equations of the mechanical model for the simple shear flow field and the dislocation equations of the dumbbell molecules, and subsequently obtain the tracer flow lines and the dumbbell distribution of the flow front. The results are used to calcu-late the stress field, analyze the rheological evolution law, study the effect of temperature and shear rate and other parameters on the model and analyze the effectiveness of the FENE dumbbell molec-ular model for the fountain effect. The study shows that due to the fountain effect, the dumbbell molecules stretch more with increasing shear rate, the polymer stress at the flow front will show a complex change of stress overshoot and then stabilize, and the dumbbell molecules are oriented along the fountain flow line.

1. 引言

大多数聚合物溶液都属于黏弹性流体，对黏弹性流体的研究关键在于确定聚合物的微观信息，采用传统的连续介质力学方法无法得到聚合物分子微观运动的信息，为了得到聚合物分子微观运动信息，可以采用聚合物分子力学模型来探究高聚物体系的流变性能，从而揭示分子的微观结构与形态，以及和宏观流变性能之间的内在联系。

在流体充填过程中，高聚物溶液流动存在两个不同的区域：主流区和前沿流区，在这两个不同的流动区域内，聚合物的流动行为有很大区别，主流区近乎是一维流动，在前沿区域的流动较为复杂，高聚物溶液由于受到剪切流动的影响，沿着流动方向排列，不断强迫流体从中心流向模壁，呈现喷泉流动形态。喷泉流动对制品的微观结构、机械和热物理特性具有非常重要的影响。

传统的控制体、拉格朗日欧拉法、水平集法无法模拟前沿的喷泉流动，需要对前沿流体进行专门描述。Rose [1] 首先用流体置换方法定义并模拟了毛细管中界面演化。E. Mitsoulis [2] 对平面和轴对称几何体中的牛顿流体进行参数化研究，全面分析了各种流体力学参数对喷泉流的影响，极大提高了喷泉流研究水平。文艳 [3] 建立了短纤维增强聚合物熔体充模过程的三维模型来描述熔体前沿的喷泉效应，实现了动态追踪纤维运动情况和取向情况，较好地反映了三维状态下熔体前沿的喷泉效应。Bechara等人 [4] 则通过采用有限体积模型和机械模型相结合的方法模拟了喷泉流动对纤维取向和分布的影响，从传统方法扩展到与机械方法相结合。Papathanasiou等人 [5] 研究了前进自由表面附近的喷泉流对填充流纤维取向的影响，系统地介绍和讨论了关于该主题的实验和计算结果。Liu等人 [6] 基于Rolie-Poly模型探究了入口速度、聚合物分子量和温度对聚苯乙烯熔体前沿喷泉演化过程的影响，模拟了聚苯乙烯的黏弹性行为。Borzenko [7] [8] [9] 通过Ostwald-de Waele流变方程研究了液体介质的假塑性和膨胀特性对圆管填充过程中非牛顿流体的喷泉流动运动影响，解释了传统方法中接触线奇异性的消除 [10] 。通过PTT-XPP模型预测了喷泉流的不稳定性，比较了流动中的弹性效应，验证了临界Weissenberg数减少对喷泉流动性的影响。Jong等人 [11] 采用可视化模具设计研究了气体反压和模具温度对喷泉流动效果的影响，推进了喷泉流动可视化研究的进展。Velmaat等人 [12] 通过使用无网格方法模拟了喷泉效应。Ren等人 [13] 采用改进的SPH方法追踪了型腔中粒子运动以验证熔体前沿的喷泉流动效应。

相较于传统方法模拟喷泉流动可能存在忽略和近似的缺点，分子动力学方法可以准确的模拟分子运动和相互作用，并且时间尺度较小，可以提供详细的数据。FENE哑铃分子模型是一种用于模拟聚合物分子运动的分子动力学方法，它建立在哑铃分子模型的基础上，通过引入额外的弹性能量来模拟聚合物分子之间的相互作用，认为分子内部的相互作用力可以由各种弹簧力或约束条件来描述 [14] ，具有效率高、精度高和灵活性好的特点。本文根据聚合物溶液流变过程的特点，提出基于分子Brown动力学原理的FENE哑铃分子模型，对聚合物稀溶液在成型流动过程存在的喷泉效应进行模拟分析求解。利用求解的构型信息和位形信息计算前沿应力的变化和演化趋势，以及各参数对模型的影响，通过喷泉流线分析验证模型的适用性。

2. 理论基础

有限可伸展非线性弹性哑铃模型(Finitely Extensible Nonlinear Elastic, FENE)，由两个小珠和一根具有一定长度的弹簧组成，如图1所示，是模拟聚合物分子链之间的相互作用的分子动力学模型，适用于聚合物稀溶液的数值模拟研究，珠子代表着溶液中相互作用的质点，不计质量的弹簧代表分子链的约束， $r_1,r_2$ 表示珠子的位形向量， $Q$ 表示珠子之间的哑铃连接向量， $r_c$ 表示哑铃分子链的质量中心。刚性珠之间的弹簧约束力代表着熵弹性的影响，是和哑铃构型向量直接相关的力。在FENE哑铃分子模型中，根据Warner弹性定律给出的有限可伸展非线性弹簧的有效约束力 $F_i^s$ 表示为 [15]

Figure 1. Schematic diagram of the dumbbell molecular model

图1. 哑铃分子模型示意图

$F_i^s=\frac{H_s{Q}_i}{1-{\left(Q_i/Q_0\right)}^2}$ (1)

其中 $Q_i=r_{i+1}-r_i$ 为弹簧的连接向量， $H_s$ 为弹簧常数， $Q_0$ 表示每根弹簧的最大可伸长长度。

高聚物分子的流变性能主要受到黏性拖拽、弹簧约束作用、布朗力作用、分子内部相互作用、排除体积、分子链缠结作用等的影响，忽略流动过程中的分子内部相互作用力、排除体积和分子缠结的作用，同时考虑到前沿流动的“喷泉”效应，对哑铃分子进行受力分析时需要计算剪切应力对喷泉流线的影响，前沿流动如图2所示。

忽略分子链的惯性效应和质量力，动态链的控制方程可以由作用在刚性珠子上的外部力总和为零的方程获得 [16]

$F_i^{drag}+F_i^B+F_i^s+F_i^{\tau }=0$ (2)

上式中的四个力分别表示为液体阻力(包括表面张力)，弹簧约束力，布朗力和剪切应力。

珠子在溶剂中运动时，液体阻力是由于分子链在流场中运动时与溶剂分子之间存在相对速度差而产生的摩擦力，以及流动前沿自由表面所产生的表面张力，刚性珠受到的Stokes阻力作用

Figure 2. Schematic diagram of the frontier area

图2. 前沿区域示意图

$F_i^{drag}=-\xi \left({\stackrel{\dot{}}{r}}_i-v\left(r_i\right)\right)$ (3)

上式中， $\xi$ 表示为液体的阻力系数， ${\stackrel{\dot{}}{r}}_i$ 表示为刚性小珠的速度， $v\left(r_i\right)$ 表示第i个珠子处高分子溶剂的速度。

在剪切流场中会产生指向模壁方向的剪切应力，引起向模壁方向的径向运动，剪切应力表示为：

$F_i^{\tau }=\mu \frac{\text{d}{v}_x}{\text{d}y}=\mu \frac{\Delta x}{\Delta t\Delta y}$ (4)

布朗力是引起高分子溶剂分子对高分子链在短时间内的迅速碰撞的原因，满足

$\langle F_i^B\left(t\right)\rangle =0$ (5)

$\langle F_i^B\left(t\right)F_j^B\left(t+\Delta t\right)\rangle =6k_{B}T{\delta }_{ij}\zeta \delta \left(\Delta t\right)\approx \frac{2k_{B}T{\delta }_{ij}\zeta }{\Delta t}$ (6)

上式中 $k_B$ 表示玻尔兹曼常数，T绝对温度， $\delta$ 表示delta函数， $\langle \cdots \rangle$ 表示为系综平均，在时间范围 $\text{d}t$ 内可以用(6)中的最后一项表示。

将公式(1)、(3)、(4)、(6)代入方程(2)并整理后，可以得到以下控制珠子位置演化的方程：

$\frac{\text{d}{r}_i}{\text{d}t}=v\left(r_i\right)+\frac{1}{\xi }\left(F_i^B+F_i^s+F_i^{\tau }\right)$ (7)

$\text{d}{r}_i=\left[v\left(r_i\right)+\frac{F_i^s+F_i^{\tau }}{\xi }\right]\text{d}t+{\left[\frac{2k_{B}T\text{d}t}{\zeta }\right]}^{1/2}\cdot n$ (8)

上式中， $n$ 表示维纳过程的随机二维向量， $n$ 的每一分量取高斯分布随机数，其均值为0，方差为1，在计算中可以用处于[−1, 1]之间的一个随机数代替，这样既可以节省计算时间，同时又对解的收敛性不会产生很大的影响。

以及单个哑铃分子的质心方程：

$r_c=\frac{r_1+r_2}{2}$ (9)

刚性小珠的位形演化方程相减获得哑铃的位形控制方程：

$\text{d}Q=\left[v\left(Q\right)+\eta \frac{\frac{\Delta r_{2x}^{}}{\Delta r_{2y}^{}}-\frac{\Delta r_{1x}^{}}{\Delta r_{1y}^{}}}{\xi \Delta t}-\frac{2}{\xi }\frac{H_{s}Q}{1-{\left(Q/Q_0\right)}^2}\right]\text{d}t+{\left[\frac{2k_{B}T\text{d}t}{\xi }\right]}^{1/2}\cdot \left(n_2-n_1\right)$ (10)

其中单个哑铃分子的与流动方向的取向角可以表示为：

$\theta =\arctan \left(\frac{Q_y}{Q_x}\right)$ (11)

上式中 $x,y$ 分别表示流动方向和速度梯度方向。

聚合物溶液的偏应力张量由溶剂分子和聚合物分子两部分贡献所得，聚合物应力张量表示为：

${\tau }_P=-n_0\langle Q{F}^E\rangle +n_0{k}_{B}T\delta$ (12)

上式中 $n_0$ 表示单位体积内聚合物分子数目， $\delta$ 表示为单位张量。

3. 数值方法

3.1. 位形计算

对于FENE哑铃分子模型，可分别引入时间常数 $\frac{\xi }{4H_{s}t}$ ，距离常数 $\sqrt{kT/H}$ ，和一个无量纲参数 $b=H_s{Q}_0^2/k_{B}T$ 。

对于控制哑铃的位置演化方程，对公式(8)和公式(10)两端同除距离常数 $\sqrt{kT/H}$ ，然后对等式右侧再乘时间常数 ${\lambda }_H=\frac{\xi }{4H_{s}t}$ ，完成无量纲化的处理，经过无量纲化得到 $\text{d}\bar{Q}$ ， $\text{d}{\bar{r}}_1$ ， $\text{d}{\bar{r}}_2$ 。为了方便起见，仍然将 $\bar{Q}$ 记为 $Q$ ， $\bar{r}$ 记为 $r$ ，如公式(13)、(14)、(15)所示。

$\text{d}Q=\left[{\lambda }_H\cdot v\left(Q\right)+\eta \frac{\frac{\Delta r_{2x}}{\Delta r_{2y}}-\frac{\Delta r_{1x}}{\Delta r_{1y}}}{\Delta t\sqrt{kT/H}}-\frac{1}{2}\frac{Q}{1-Q^2/b}\right]\text{d}t+{\left[\frac{\text{d}t}{2}\right]}^{1/2}\cdot \left(n_2-n_1\right)$ (13)

$\text{d}{r}_1=\left[{\lambda }_{H}v\left(r_1\right)+\frac{1}{4}\frac{Q}{1-Q^2/b}+\eta \frac{\Delta r_{1x}}{\sqrt{kT/H}\Delta t\cdot \Delta r_{1y}}\right]\text{d}t+{\left[\frac{\text{d}t}{2}\right]}^{1/2}\cdot n$ (14)

$\text{d}{r}_2=\left[{\lambda }_{H}v\left(r_2\right)-\frac{1}{4}\frac{Q}{1-Q^2/b}+\eta \frac{\Delta r_{2x}}{\sqrt{kT/H}\Delta t\cdot \Delta r_{2y}}\right]\text{d}t+{\left[\frac{\text{d}t}{2}\right]}^{1/2}\cdot n$ (15)

分析得到哑铃分子的水平构形方程为：

$Q_x\left(t+\Delta t\right)=Q_x\left(t\right)+\left({\lambda }_H\cdot v\left(Q\right)-\frac{1}{2}\cdot \frac{Q_x}{1-Q^2/b}\right)\cdot \text{d}t+{\left(\frac{\text{d}t}{2}\right)}^{1/2}\cdot \left(n_2-n_1\right)$ (16)

竖直方向构形方程为：

$Q_y\left(t+\Delta t\right)=Q_y\left(t\right)+\eta \frac{\frac{\Delta r_{2x}}{\Delta r_{2y}}-\frac{\Delta r_{1x}}{\Delta r_{1y}}}{\sqrt{kT/H}}-\frac{1}{2}\cdot \frac{Q_y\cdot \text{d}t}{1-Q^2/b}+{\left(\frac{\text{d}t}{2}\right)}^{1/2}\cdot \left(n_2-n_1\right)$ (17)

材料黏度无量纲化表示为 [14] ：

$\frac{\eta \left(\gamma \right)-{\eta }_s}{n_0{k}_{B}T{\lambda }_H}=-\frac{1}{{\lambda }_H\gamma }\frac{{\tau }_{p,xy}}{n_0{k}_{B}T}$ (18)

哑铃初始位形在 $\left[0,\sqrt{b/2}\right]$ 范围内随机获取得到， $r_1$ 初始位置沿中心面往模壁方向在x = 0处随机获取，然后通过计算 $r_2=r_1+Q$ 得到珠子的位置矢量 $r_2$ ，首次计算时剪切应力取 $F=\eta \stackrel{\dot{}}{\gamma }$ ，之后按照公式(16)和公式(17)迭代运算，代入 $Q$ 计算下一步的构型方程 $Q$ ，从而得到所有时刻的构型向量和位形向量。

3.2. 材料参数与计算方案

本文根据流体填充过程简单剪切流场中流体前沿喷泉效应，假设在前沿喷泉效应区域内稀溶液不可压缩，忽略溶质分子的质量力和惯性力，忽略流动过程中的热对流和热传导，采用欧拉显式差分法模拟流场及流动前沿。哑铃分子模型中材料计算模拟参数如下表1所示。

Table 1. Material parameter values

表1. 材料参数取值

聚合物稀溶液中的剪切流场的速度分量设置为 $v_x=\left(1-y/h\right)\cdot \stackrel{\dot{}}{\gamma }$ ， $v_y=0$ ，其中h为中心面到模壁的长度。考虑表面张力的影响，将溶液前沿部分流场速度分量设置为 $v_x=\left(1/2-y/h\right)\cdot \stackrel{\dot{}}{\gamma }$ ， $v_y=0$ ，并以此作为流动前沿的边界条件，考虑到流场的对称性，只选取上半部分进行计算。

由于方程无法通过直接求解得到解析解，用前文所述的欧拉数值方法求解。首先获取初始哑铃位形信息，利用哑铃初始位形求解公式(14)和(15)，然后求解得到初始状态信息，进而求解方程(19)、(20)，在每一次迭代运算中对哑铃分子的位形信息排序，最前沿部分的哑铃分子采用前沿流场进行迭代计算，其他哑铃分子按照原流场速度分量进行计算，从而获取所有连接向量以及小珠位形。

$r_1\left(t+\Delta t\right)=r_1\left(t\right)+\left[{\lambda }_{H}v\left(r_1\right)+\frac{1}{4}\frac{Q}{1-Q^2/b}+\eta \frac{\Delta r_{1x}}{\sqrt{kT/H}\Delta t\cdot \Delta r_{1y}}\right]\text{d}t+{\left(\text{d}t\right)}^{1/2}\cdot n$ (19)

$r_2\left(t+\Delta t\right)=r_2\left(t\right)+\left[{\lambda }_{H}v\left(r_2\right)+\frac{1}{4}\frac{Q}{1-Q^2/b}+\eta \frac{\Delta r_{2x}}{\sqrt{kT/H}\Delta t\cdot \Delta r_{2y}}\right]\text{d}t+{\left(\text{d}t\right)}^{1/2}\cdot n$ (20)

$Q=r_2-r_1$ (21)

本文哑铃分子总数为n = 10,000，无量纲时间步长设置为1e‒3。当计算时间步长不是特别小时，由于计算误差和随机力的作用，在模拟计算过程中可能会出现弹簧长度超出最大拉伸长度的不收敛情况，本文采用两种方面处理，一是提高哑铃弹簧的最大拉伸长度；二是当出现拉伸过长的情况时采用逐次缩小时间步长的方法 [17] ，将时间步长缩小N倍，对上一状态重新计算N步，如果仍然出现超长情况，则将N扩大，继续从前一状态计算，之后仍用原本的时间步长计算，这样既可以避免持续采用较短的时间步长，节省计算时间，同时也可以通过较短的时间步长控制向量的解在收敛区间之内，计算流程如图3所示。

Figure 3. Flow chart of simulation calculation

图3. 模拟计算流程图

在入口截面处，先获取每个时刻所有质心在截面处的哑铃分子信息，之后遍历获取所有附近的哑铃分子进行系综平均，计算得到聚合物应力。同时追踪单个哑铃在运动过程中的位形变化和分子链波动，对任意时刻在前沿处的哑铃进行系综平均求取流动前沿的应力变化。

4. 结果与讨论

通过对哑铃分子位形方程进行追踪，在 ${\lambda }_H\stackrel{\dot{}}{\gamma }=1$ 情况下随机选取9个哑铃分子追踪获得运动轨迹如图4所示，与E. Mitsoulis [2] 描述的喷泉流线宏观上一致，与理论中的流动前沿喷泉效应流线相同，随着哑铃分子的运动，主流方向速度逐渐减弱且获得径向速度，慢慢涌向壁面，趋近于中心面的哑铃分子运动轨迹更符合理想中的喷泉效应，并且涌动的喷泉更加明显。

Figure 4. Dumbbell molecular tracer streamlines

图4. 哑铃分子示踪流线

模拟的x = 0.05和x = 0.06处的应力如图5所示，应力随着时间变化先增长达到一个最大值，然后逐渐变为稳态值，并且在变化过程中由于哑铃分子运动造成分子数量变化，稳定之前存在一些波动，之后逐渐趋于稳定。x = 0.05处应力达到的波动最大值比x = 0.06处要大，但前者的稳定值更小，且前者的应力波动程度更大，同时达到稳定状态的时间也更短。

Figure 5. Cross-sectional stress fluctuation with time

图5. 截面应力随时间波动

图6记录了四个不同时刻部分哑铃分子在前沿的位形分布，随着哑铃分子的运动，前沿哑铃逐渐从与流场相似的斜面分布转变为接近于自由表面的前沿曲面分布，并且逐渐趋于稳定，在流动过程中哑铃分子链发生旋转、拉伸等变化，且随着时间的变化，哑铃分子链沿运动方向的取向更加明显，与Jong [11] 的仿真和实验结果相同，与Borzenko [8] 在圆管中的模拟也保持一致，验证了FENE哑铃模型对于聚合物流动前沿的模拟有效性，但与理想状态下的半圆形仍有差距，仍需要进一步深入研究。

对不同剪切速率下高聚物前沿分子的分子链随时间的波动伸展分析如图7所示，分子伸展随着哑铃分子在流场中运动逐渐伸长并趋向于平稳， ${\lambda }_H\stackrel{\dot{}}{\gamma }=4$ 时哑铃分子伸长度趋近于0.65， ${\lambda }_H\stackrel{\dot{}}{\gamma }=1$ 时哑铃分子伸长度趋向于0.4左右，在高剪切速率下哑铃伸长的稳定值更大，且变化过程中的波动明显，低剪切速率下哑铃达到稳定态所需的时间也更长。

(a) t = 0.05 s (b) t = 0.1 s (c) t = 0.15 s (d) t = 0.3 s

Figure 6. Frontside dumbbell molecular bit shape distribution

图6. 前沿哑铃分子位形分布

Figure 7. Dumbbell molecular chain elongation fluctuates with time

图7. 哑铃分子链伸长随时间波动

图8为聚合物黏度系数变化曲线图，从图中可以看出随着剪切速率的增大，聚合物黏度呈下降趋势，出现剪切变稀现象，并且随着温度的升高，分子间相互作用力减弱且布朗运动加剧，聚合物黏度会出现下降。对于不同剪切速率下对流动前沿聚合物应力的计算结果如图9所示，当剪切速率较大时，流动前沿聚合物应力较快增长达到最大值，然后逐渐变化为稳态值，同样存在高剪切速率下应力过冲的情况，与闵志宇 [14] 对聚合物分子流变性质的模拟结果基本一致。

Figure 8. Variation of shear viscosity with shear rate and temperature

图8. 剪切黏度随剪切速率和温度的变化

Figure 9. Variation of shear stress with time

图9. 剪切应力随时间的变化

5. 结论

本文基于FENE哑铃建立了聚合物稀溶液中前沿喷泉流动模拟方法，实现了对前沿喷泉效应的模拟，记录了哑铃的位形分布，得到了较为理想的喷泉流线，说明FENE哑铃分子模型适用于聚合物稀溶液前沿流动的模拟，且对流动过程中截面应力以及前沿应力进行了模拟分析，应力随着区域内分子量的变化也产生了较为复杂的应力变化曲线，截面处的应力受到前沿喷泉效应的影响会产生波动，之后趋于稳定。前沿剪切应力随着剪切速率的增大会出现应力过冲，且剪切速率越大，峰值越大，峰值区间越小，越靠近中心面位置的哑铃位形受到喷泉效应影响越大。前沿分子一开始沿主流动方向取向，伴随喷泉效应的发生，逐步沿喷泉流线方向取向，最后趋向于稳定值，分子链伸长增大逐渐趋于稳定，且前沿流动受到表面张力的影响较弱。

基金项目

国家自然科学基金(11672271、11432003)项目资助。