粘弹性流体流的有限元方法

在 Simcenter STAR-CCM+ 中，用于捕捉非牛顿粘弹性流体的复杂流场的数值算法基于有限元方法。对闭合形式的非线性微分本构方程离散化的描述可借鉴 Oldroyd-B 模型，并且该方程可随时扩展为其他微分本构方程。

由连续性方程、动量方程和本构方程组成的控制方程组无量纲。耦合方程组使用特定的混合有限元方法（采用三种稳定化方法）进行离散化。这些稳定化方法使用应力-速度耦合、非线性对流项和流体的不可压缩性约束解决问题。

在稳态和等温条件下，可以无量纲形式编写不可压缩粘弹性流体流的控制方程：

连续性方程

$$\nabla \cdot v=0$$

(1014)

动量方程

$$\mathrm{Re}(v\cdot \nabla v)=-\nabla p+\nabla \cdot T$$

(1015)

本构方程

为描述数值计算过程，需要通过 Oldroyd B 模型表示粘弹性本构方程：

$$\mathrm{Wi}{\stackrel{▿}{T}}_i^p+T_i^p=2{\zeta }_{i}D$$

(1016)

$$T=2{\zeta }_{s}D+{\displaystyle \sum _{i=1}^N{T}_i^p}$$

(1017)

其中：

$v$	无量纲速度
$p$	无量纲压力
$T$	无量纲粘性应力张量
$\mathrm{Re}=\frac{\rho v_c{l}_c}{{\mu }_c}$	雷诺数
$v_c$	特征速度
$l_c$	特征长度
${\mu }_c$	特征动力粘度
${\mathrm{Wi}}_i=\frac{{\lambda }_i{v}_c}{l_c}$	与第 $i$ 模式关联的 Weissenberg 数
${\lambda }_i$	第 $i$ 模式的松弛时间
${\zeta }_i={\mu }_i/{\mu }_c$	第 $i$ 模式的无量纲粘度
${\mu }_i$	与第 $i$ 模式关联的粘度
${\mu }_s$	溶剂粘度

上述方程通过 $v\sim v_c,l\sim l_c,p,T\sim {\mu }_c{v}_c/l_c$ 无量纲化。

采用这种形式时，由于速度和压力没有显式出现在动量和连续性方程中，因此上述方程组将产生双鞍点问题。虽然速度通过溶剂粘度项仅显式出现在动量方程中，但是对于粘弹性流体，此项的贡献通常可忽略不计。要解决此问题，需要使用 Guénette 和 Fortin [175] 提出的离散弹粘性分割应力 (DEVSS) 方法。该方法定义了额外的张量未知量 $d$ ：

$$T\to 2{\zeta }_{s}D(v)+{\displaystyle \sum _{i=1}^N{T}_i^p}-2\alpha d+2\alpha D(v)$$

(1018)

其中，$\alpha$ 为任意常数。

此替代将生成以下 $(v,p,T_i^p,d)$ 公式：

$$\nabla \cdot v=0$$

(1019)

$$\text{Re}(v\cdot \nabla v)+\nabla p-(2{\zeta }_s+2\alpha )\nabla \cdot D(v)-{\displaystyle \sum _{i=1}^N\nabla \cdot T_i^p}+2\alpha \nabla \cdot d=0$$

(1020)

$$\begin{array}{c}{\mathrm{W}\mathrm{i}}_i{T}_i^p+T_i^p-2{\xi }_{i}D(v)=0\end{array}$$

(1021)

$$\mathbf{d}-\mathbf{D}\begin{array}{c}\left(\mathbf{v}\right)=0\end{array}$$

(1022)

变分公式

要获取变分公式，需要将上述方程组乘以加权函数 $({\tilde{p}}_h,{\tilde{v}}_h,{\tilde{T}}_h,{\tilde{d}}_h)$ ，然后对具有闭合边界 $\mathrm{\Gamma }$ 和法线 $n_{\mathrm{\Gamma }}$ 的体积域 $\Omega$ 进行积分。

通过将标准 Galerkin 方法用于 $(p_h,v_h,T_{i,h}^p,d_h)$ 的方程组并应用散度定理，方程组由以下公式给出：

$${\displaystyle \underset{\Omega }{\int }\nabla \cdot v_h {\tilde{p}}_h d\Omega }=0, \forall {\tilde{p}}_h\in H_0^1(\Omega )$$

(1023)

$$\text{Re}{\displaystyle \underset{\Omega }{\int }(v_h\cdot \nabla v_h)\cdot {\tilde{v}}_h d\Omega }-{\displaystyle \underset{\Omega }{\int }{p}_h\nabla \cdot {\tilde{v}}_h d\Omega }+(2{\zeta }_s+2\alpha ){\displaystyle \underset{\Omega }{\int }D(v_h):D({\tilde{v}}_h) d\Omega }+{\displaystyle \sum _{m=1}^N{\displaystyle \underset{\Omega }{\int }{T}_i^p:D({\tilde{v}}_h) d\Omega }}-2\alpha {\displaystyle \underset{\Omega }{\int }d:D({\tilde{v}}_h) d\Omega }-{\displaystyle \underset{\Gamma }{\int }{\mathbf{n}}_{\mathrm{\Gamma }}\cdot (-p_{h}I+2{\zeta }_{s}D(v_h)+{\displaystyle \sum _{m=1}^N{T}_i^p}+2\alpha (D(v_h)-d))\cdot {\tilde{v}}_h d\Gamma }=0, \forall {\tilde{v}}_h\in {(H_0^1(\Omega ))}^3,$$

(1024)

$${\mathrm{We}}_i{\displaystyle \underset{\Omega }{\int }{\stackrel{▿}{T}}_{i,h}^p:{\mathbf{T}}_h d\Omega }+{\displaystyle \underset{\Omega }{\int }{T}_{i,h}^p:{\mathbf{T}}_h d\Omega }-2{\zeta }^m{\displaystyle \underset{\Omega }{\int }D(v_h):{\mathbf{T}}_h d\Omega }=0, \forall {\mathbf{T}}_h\in {(H_0^1(\Omega ))}^6$$

(1025)

$${\displaystyle \underset{\Omega }{\int }{d}_h:d_h d\Omega }-{\displaystyle \underset{\Omega }{\int }{d}_h:D(v_h) d\Omega }=0, \forall d_h\in {(H_0^1(\Omega ))}^6$$

(1026)

其中， $H_0^k\left(\mathrm{\Omega }\right)$ 表示由平方可积函数（高达 $k$ 阶的所有导数均为平方可积）组成的 Sobolev 空间，在具有重要边界条件的边界处的幅值为零，且 ${\mathbf{n}}_{\mathrm{\Gamma }}$ 为向外指向包围域 $\mathrm{\Omega }$ 的边界 $\mathrm{\Gamma }$ 的单位法向矢量。

连续空间域将离散化成在节点处互连的有限数量的单元。

具有无轨迹梯度的离散弹粘性分割应力 (DEVSS/TG)

在任何不可压缩流中，梯度必须保持无轨迹，即 $\text{tr}\left(\nabla \mathbf{v}\right)=\nabla \cdot \mathbf{v}=0$ 。但是，根据有限元方法计算的离散速度场并不是精确无迹的。特别是，在速度发生突然变化的区域中，插值速度梯度的轨迹变得非常大；这可能会导致数值模拟不稳定。请参见 Pasquali 和 Scriven [200]。

为了保证速度梯度的无发散性，Pasquali 和 Scriven 修改了 Guénette 和 Fortin [175] 的离散弹粘性分割应力 (DEVSS)，从而使速度梯度变为无轨迹，因此它被称为 DEVSS/TG。Eqn. (1022) 中的变形率张量 $\mathbf{D}\left(\mathbf{v}\right)$ 定义为 $\mathbf{D}\left(\mathbf{v}\right)=\frac{1}{2}(\nabla \mathbf{v}+\nabla {\mathbf{v}}^{\text{T}})$ 。要使速度梯度 $\nabla \mathbf{v}$ 变得无轨迹，将 $\nabla \mathbf{v}$ 替换为：

$$\nabla \mathbf{v}-\frac{\nabla \cdot \mathbf{v}}{\text{tr}\left(\mathbf{I}\right)}\mathbf{I}$$

(1027)

因此：

$$\mathbf{D}\left(\mathbf{v}\right)=\frac{\nabla \mathbf{v}+\nabla {\mathbf{v}}^{\text{T}}}{2}-\frac{\nabla \cdot \mathbf{v}}{\text{tr}\left(\mathbf{I}\right)}\mathbf{I}$$

(1028)

然后，Eqn. (1022) 可以重写为：

$$\mathbf{d}-\mathbf{D}+\frac{\nabla \cdot \mathbf{v}}{\text{tr}\left(\mathbf{I}\right)}\mathbf{I}=0$$

(1029)

非等温流

对于非等温流，无量纲形式的均匀导电率能量方程为：

$$\text{Pe}(\frac{\partial T}{\partial t}+\mathbf{v}\cdot \nabla T)-{\nabla }^{2}T-\text{Br}\mathbf{T}:\nabla \mathbf{v}=0$$

(1030)

其中：

• $T$ 为无量纲流体温度，使用域 $T_c$ 的特征温度进行缩放。
• $\mathbf{T}$ 为无量纲总应力张量。
• $\text{Pe}$ 和 $\text{Br}$ 分别是 Péclet 数和 Brinkman 数：
  $$\begin{array}{c}\text{Pe}=\frac{v_c{l}_c}{k/\rho C_p}\\ \text{Br}=\frac{{\mu }_c{v}_c^2}{k{T}_c}\end{array}$$
  且：
  • $v_c$ 为特征速度。
  • $k$ 为流体导电率。
  • $C_p$ 为比热容。
  • $T_c$ 为域的特征温度。

要获取变分形式，需要将 Eqn. (1030) 乘以加权函数 ${\tilde{T}}_h$ ，然后对具有闭合边界 $\mathrm{\Gamma }$ 和法线 ${\mathbf{n}}_{\mathrm{\Gamma }}$ 的体积域 $\mathrm{\Omega }$ 进行积分。通过使用标准 Galerkin 方法并应用散度定理，Eqn. (1030) 可以写为：

$$\begin{gathered} \underset{\mathrm{\Omega }}{\int }\text{Pe}{\tilde{T}}_h\frac{\partial T}{\partial t}d\mathrm{\Omega }+\underset{\mathrm{\Omega }}{\int }\text{Pe}({\mathbf{v}}_h\cdot \nabla T_h){\tilde{T}}_{h}d\mathrm{\Omega }+\underset{\mathrm{\Omega }}{\int }\nabla T_h\cdot {\tilde{T}}_{h}d\mathrm{\Omega }-\underset{\mathrm{\Omega }}{\int }\text{Br}({\mathbf{T}}_h:\nabla {\mathbf{v}}_h){\tilde{T}}_{h}d\mathrm{\Omega }-\underset{\mathrm{\Gamma }}{\int }(\nabla T_h\cdot {\mathbf{e}}_n){\tilde{T}}_{h}d\mathrm{\Gamma }=0, \\ \forall {\tilde{T}}_h\in H_0^1\left(\mathrm{\Omega }\right) \end{gathered}$$

(1031)

其中， ${\mathbf{e}}_n$ 为无量纲单位法向。

在非等温流中，Eqn. (1031) 通过连续性、动量和本构方程以耦合方式求解。

最高温度和最低温度

聚合物的瞬态非等温模拟通常很复杂，因为聚合物的导电率一般较小，Péclet 数非常高。这通常会导致具有非常薄的热边界层的壁面旁边的数值不稳定。数值不稳定显示为壁面温度的尖锐峰值。通常，温度峰值非常明显（有时会导致产生局部负温度），导致模拟立即崩溃。通过引入温度的最小值和最大值，限制此峰值的幅值。如果温度高于（或低于）指定的 $T_{\text{max}}$ （或 $T_{\text{min}}$ ），将在局部向能量方程添加能量汇（或源）项，以冷却（或加热）局部温度。该项采用以下形式：

$$\pm \frac{\alpha \rho C_p\mathrm{\Delta }T}{\mathrm{\Delta }t}$$

(1032)

其中， $\alpha$ 是一个较小的无量纲值，用于确定容差度，且 $\mathrm{\Delta }t$ 为时间步。加（或减）表示冷却（或加热）。要在应局部包含该项的域中标识网格单元，根据以下平滑函数测试任意高斯点的局部温度：

$$\begin{array}{c}\mathrm{\Delta }T=\frac{(T-T_{\text{max}})+\sqrt{{(T-T_{\text{max}})}^2+{ϵ}^2}}{2}\\ \mathrm{\Delta }T=\frac{-(T-T_{\text{min}})+\sqrt{{(T-T_{\text{min}})}^2+{ϵ}^2}}{2}\end{array}$$

(1033)

其中， $ϵ=$ 0.01。对于温度低于（高于） $T_{\text{max}}$ ( $T_{\text{min}}$ ) 的所有点，上述函数将生成 0。

动量和能量的源项

将体积力 ${\mathbf{f}}_b$ 引入动量方程会将以下形式的求积添加到 Eqn. (1024) 的左侧：

$$-\underset{\mathrm{\Omega }}{\int }{\mathbf{f}}_b\cdot {\tilde{\mathbf{v}}}_{h}d\mathrm{\Omega },\forall {\tilde{\mathbf{v}}}_h\in H_0^1\left(\mathrm{\Omega }\right)$$

(1034)

类似地，将能量源 $q$ 引入能量方程会将以下项添加到 Eqn. (1031) 的左侧：

$$-\underset{\mathrm{\Omega }}{\int }q{\tilde{T}}_{h}d\mathrm{\Omega },\forall \tilde{T}\in H_0^1\left(\mathrm{\Omega }\right)$$

(1035)

粘弹性流体流中的被动标量

被动标量被视为粘性流体的组分，这种流体可在流场作用下平流且可根据 Fick 定律扩散，但是它们出现时的浓度非常小，因此不会改变流场。存在多个被动标量时，组分之间没有任何相互作用。系统中每个被动标量的分布由平流-扩散方程描述：

$${\text{Pe}}_i(\frac{\partial C_i}{\partial t}+\mathbf{v}\cdot \nabla C_i)={\nabla }^2{C}_i$$

(1036)

其中：

• $C_i$ 为系统中第 个组分的无量纲分布。 $i$
• ${\text{Pe}}_i$ 为 Péclet 数 ${Pe}_i=v_c{l}_c/D_i$ 。
• $D_i$ 为第 个组分的分子扩散率，（对于小分子组分）通过 Stokes-Einstein 关系获得。 $i$
• $t$ 为无量纲物理时间，使用平流时间尺度 $t_c\sim a/v_c$ 。
• $\mathbf{v}$ 为无量纲速度矢量。

此方程单向耦合到流体动力方程。

要获取被动标量的变分公式，需要将 Eqn. (1036) 乘以加权函数 ${\tilde{C}}_i$ ，然后对体积域 $\mathrm{\Omega }$ 积分。由此得出：

$${\int }_{\mathrm{\Omega }}^{}{Pe}_i{\tilde{C}}_i\frac{\partial C_i}{\partial t}d\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}^{}{Pe}_i{\tilde{C}}_i(\mathbf{v}\cdot \nabla C_i)d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}^{}{\tilde{C}}_i{\nabla }^2{C}_{i}d\mathrm{\Omega }=0,\forall {\tilde{C}}_i\in H_0^1\left(\mathrm{\Omega }\right)$$

(1037)

分部积分最后一项，然后根据散度定理推导：

$${\int }_{\mathrm{\Omega }}^{}{Pe}_i{\tilde{C}}_i\frac{\partial C_i}{\partial t}d\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}^{}{Pe}_i{\tilde{C}}_i(\mathbf{v}\cdot \nabla C_i)d\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}^{}\nabla {\tilde{C}}_i\cdot \nabla C_{i}d\mathrm{\Omega }-{\int }_{\mathrm{\Gamma }}^{}{\tilde{C}}_i({\mathbf{n}}_{\mathrm{\Gamma }}\cdot \nabla C_i)d\mathrm{\Omega }=0,\forall {\tilde{C}}_i\in H_0^1\left(\mathrm{\Omega }\right)$$

(1038)

被动标量方程不会修改流场，但它无法独立求解，因为必须已知速度场才能同时求解上述离散方程与离散流体动力方程。