Finite Element Method for Viscoelastic Fluid Flow

In Simcenter STAR-CCM+, the numerical algorithm for capturing the complex flow field of non-Newtonian viscoelastic fluids is based on the finite element method. The discretization of a closed-form non-linear differential constitutive equation is described exemplarily by means of the Oldroyd-B model and can be readily extended to other differential constitutive equations.

The governing set of equations consisting of continuity equation, momentum equations, and constitutive equation is non-dimensionalized. The coupled system of equations is discretized by using a particular mixed finite element method, which employs three stabilization techniques. These stabilization techniques address issues with the stress-velocity coupling, the non-linear convection terms, and the incompressibility constraint of the fluid.

In steady-state and under isothermal conditions, the governing equations for an incompressible viscoelastic fluid flow can be written in non-dimensional form:

Continuity Equation

$$\nabla \cdot v=0$$

(1014)

Momentum Equation

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

(1015)

Constitutive Equation

To describe the numerical procedure, the viscoelastic constitutive equation is represented by the Oldroyd B model:

$$\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)

where:

$v$	non-dimensional velocity
$p$	non-dimensional pressure
$T$	non-dimensional viscous stress tensor
$\mathrm{Re}=\frac{\rho v_c{l}_c}{{\mu }_c}$	Reynolds number
$v_c$	characteristic velocity
$l_c$	characteristic length
${\mu }_c$	characteristic dynamic viscosity
${\mathrm{Wi}}_i=\frac{{\lambda }_i{v}_c}{l_c}$	Weissenberg number that is associated with the $i$th mode
${\lambda }_i$	relaxation time of the $i$th mode
${\zeta }_i={\mu }_i/{\mu }_c$	dimensionless viscosity of the $i$th mode
${\mu }_i$	viscosity that is associated with the $i$th mode
${\mu }_s$	solvent viscosity

The above equations are non-dimensionalized by $v\sim v_c,l\sim l_c,p,T\sim {\mu }_c{v}_c/l_c$.

In this form, the above system of equations yields a double saddle-point problem, because the velocity and pressure do not appear in the momentum and continuity equations explicitly. Although the velocity appears explicitly only in the momentum equation through the solvent viscosity term, the contribution of this term is often negligible for viscoelastic fluids. To overcome this difficulty, the Discrete Elastic Viscous Split Stress (DEVSS) method proposed by Guénette and Fortin [175] is used. The method defines an additional tensorial unknown $d$ :

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

(1018)

where $\alpha$ is an arbitrary constant.

This substitution results in the following $(v,p,T_i^p,d)$ formulation:

$$\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)

Variational Formulation

To obtain the variational formulation, the above set of equations is multiplied with the weight functions $({\tilde{p}}_h,{\tilde{v}}_h,{\tilde{T}}_h,{\tilde{d}}_h)$ and then integrated over the volume domain $\Omega$ with the enclosing boundary $\mathrm{\Gamma }$ and normal $n_{\mathrm{\Gamma }}$ .

By using the standard Galerkin approach for the system of $(p_h,v_h,T_{i,h}^p,d_h)$ and applying the divergence theorem, the system of equations is given by:

$${\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)

where $H_0^k\left(\mathrm{\Omega }\right)$ denotes the Sobolev space consisting of square integrable functions for which all derivatives up to order $k$ are square integrable, with the magnitude zero at the boundaries with essential boundary conditions and ${\mathbf{n}}_{\mathrm{\Gamma }}$ is the unit normal vector pointing outward to the boundary $\mathrm{\Gamma }$ enclosing the domain $\mathrm{\Omega }$ .

The continuous space domain is discretized into a finite number of elements, which are interconnected at the vertices.

Discrete Elastic Viscous Split Stress with Traceless Gradient (DEVSS/TG)

In any incompressible flow, the gradient must remain traceless, i.e. $\text{tr}\left(\nabla \mathbf{v}\right)=\nabla \cdot \mathbf{v}=0$ . However, the discretized velocity field computed according to the finite element method is not precisely traceless. In particular, the trace of the interpolated velocity gradient becomes extremely large in regions where there is a sudden change in velocity; this may lead to instability of the numerical simulation. See Pasquali & Scriven [200].

To guarantee the divergence-free property of the velocity gradient, Pasquali & Scriven modified the Discrete Elastic Viscous Split Stress (DEVSS) of Guénette & Fortin [175] so that the velocity gradient becomes traceless and thus it is called DEVSS/TG. The rate of deformation tensor $\mathbf{D}\left(\mathbf{v}\right)$ in Eqn. (1022) is defined as $\mathbf{D}\left(\mathbf{v}\right)=\frac{1}{2}(\nabla \mathbf{v}+\nabla {\mathbf{v}}^{\text{T}})$ . To make the velocity gradient $\nabla \mathbf{v}$ traceless, replace $\nabla \mathbf{v}$ with:

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

(1027)

As a result :

$$\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)

Then Eqn. (1022) can be re-written as:

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

(1029)

Non-Isothermal Flow

For non-isothermal flow, the energy equation for uniform conductivity in the non-dimensional form is:

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

(1030)

where:

• $T$ is the non-dimensional fluid temperature, scaled with the characteristic temperature of the domain $T_c$ .
• $\mathbf{T}$ is the non-dimensional total stress tensor.
• $\text{Pe}$ and $\text{Br}$ are the Péclet number and Brinkman numbers respectively:
  $$\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}$$
  and:
  • $v_c$ is the characteristic velocity.
  • $k$ is the fluid conductivity.
  • $C_p$ is the specific heat capacity.
  • $T_c$ is the characteristic temperature of the domain.

To obtain a variational form, multiply Eqn. (1030) by the weight function ${\tilde{T}}_h$ and then integrate over the volume domain $\mathrm{\Omega }$ with the enclosing boundary $\mathrm{\Gamma }$ and normal ${\mathbf{n}}_{\mathrm{\Gamma }}$ . By using the standard Galerkin approach and applying the divergence theorem, Eqn. (1030) can be written as:

$$\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)

where ${\mathbf{e}}_n$ is the dimensionless unit normal.

In non-isothermal flow, Eqn. (1031) is solved with the continuity, momentum and constitutive equation in a coupled way.

Maximum and Minimum Temperature

The transient non-isothermal simulation of polymers is often complex, since the conductivity of the polymers is usually small and the Péclet number is very high . This often leads to numerical instability next to walls, where there is a very thin thermal boundary layer. The numerical instability appears as sharp spikes in temperature by the walls. The temperature spike is often significant (sometimes result in a local negative temperature), causing the simulation to collapse immediately. The magnitude of this spike is limited by introducing minimum and maximum values for temperature. If the temperature is greater (or lesser) than the specified $T_{\text{max}}$ (or $T_{\text{min}}$ ), a sink (or source) term is added locally to the energy equation, to cool down (or heat up) the local temperature. The term is in the following form:

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

(1032)

where $\alpha$ is a small non-dimensional value determining the degree of tolerance and $\mathrm{\Delta }t$ is the time-step. Plus (or minus) accounts for cooling (or heating). To identify the cell in the domain where this term should be locally incorporated, the local temperature at any Gaussian point is tested according to the following smooth function:

$$\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)

where $ϵ=$ 0.01. The above functions generate zero for all points for which the temperature is below (above) $T_{\text{max}}$ ( $T_{\text{min}}$ ).

Source Terms for Momentum and Energy

Introducing a body force ${\mathbf{f}}_b$ to the momentum equation adds a quadrature in the following form to the left-hand side of 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)

Similarly, introducing an energy source $q$ to the energy equation adds the following term to the left-hand side of 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)

Passive Scalars in Viscoelastic Fluid Flow

Passive scalars are treated as components of the viscous fluid that can be advected by the flow field and diffused according to Fick's law, but occur in concentrations so small that their presence does not alter the flow field. With multiple passive scalars, there is no interaction between the components. The distribution of each passive scalar in the system is described by an advection-diffusion equation:

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

(1036)

where:

• $C_i$ is the non-dimensional distribution of the $i$ th species in the system.
• ${\text{Pe}}_i$ is the Péclet number, ${Pe}_i=v_c{l}_c/D_i$ .
• $D_i$ is the molecular diffusivity of the $i$ th species that (for small molecule species) is obtained by the Stokes-Einstein relation.
• $t$ is non-dimensional physical time scaled with the advection time scale $t_c\sim a/v_c$ .
• $\mathbf{v}$ is the non-dimensional velocity vector.

This equation is one-way coupled to the hydrodynamic equations.

To obtain a variational formulation for passive scalars, multiply Eqn. (1036) by the weight function ${\tilde{C}}_i$ and then integrate over the volume domain $\mathrm{\Omega }$ . This gives:

$${\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)

Integrate the last term by parts and use the divergence theorem to derive:

$${\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)

The passive scalar equation does not modify the flow field, but it cannot be solved independently, since a velocity field must be known in order to solve the above discretized equation simultaneously with the discretized hydrodynamic equations.