# A New Method to Simulate Free Surface Flows for Viscoelastic Fluid

Authors: Yu Cao; Xiao-Guang Ren; Xiao-Wei Guo; Miao Wang; Qian Wang; Xin-Hai Xu; Xue-Jun Yang
Journal: Advances in Materials Science and Engineering
 (2015)

---

## Abstract

Free surface flows arise in a variety of engineering applications. To predict the dynamic characteristics of such problems, specific numerical methods are required to accurately capture the shape of free surface. This paper proposed a new method which combined the Arbitrary Lagrangian-Eulerian (ALE) technique with the Finite Volume Method (FVM) to simulate the time-dependent viscoelastic free surface flows. Based on an open source CFD toolbox called OpenFOAM, we designed an ALE-FVM free surface simulation platform. In the meantime, the die-swell flow had been investigated with our proposed platform to make a further analysis of free surface phenomenon. The results validated the correctness and effectiveness of the proposed method for free surface simulation in both Newtonian fluid and viscoelastic fluid.


---

## Body

## 1. Introduction

In recent years, there is a growing interest in free surface simulations (e.g., die-swell and bubble growing). Die-swell, as a typical phenomenon seen in viscoelastic fluids, is widely investigated. Numerical simulations of extrusion flow need to address free surface problems. Numerical schemes for free surface treatment can be classified into the fixed grid (interface capturing) and the moving grid (interface tracking) methods. Fixed grid methods capture the interface in a Eulerian description, including Volume of Fluid (VOF) method [1], Level Set method [2], and Marker and Cell (MAC) method [3], in which the interface is reconstructed through numerical schemes in these methods. Though these methods could address sharp interface and large deformation issues well, the position of the free surface is usually not precise enough. On the contrary, moving grid methods explicitly track the position of interface by computational grids which guarantees the precision of the interface position, but the quality of the mesh might be reduced during mesh moving in large deformation cases. Arbitrary Lagrangian-Eulerian (ALE) Method [4] is one of the moving grid methods. It combines the advantage of the Lagrangian method in tracking interface and the efficiency of Eulerian method in solving governing equations. Therefore, it has been widely used in wave-like free surface simulations.ALE method is usually combined with Finite Element Method (FEM) [5], in which the traction boundary condition on the interface could be treated in a neat and exact way by integrating the force in the momentum equation by parts. However, with the development of Finite Volume Method (FVM), FVM system is more attractive in fluid simulations for its good property in local conservation. Therefore, combining ALE method with FVM is capable of representing the position of free surface by the computational mesh directly, while benefitting from a good conservation property with FVM. Tuković and Jasak [6] developed a platform to simulate Newtonian fluid free surface problems by using the ALE method in FVM system. In order to get the accurate velocity on the surface, a new surface mesh is coupled with the traditional mesh, but it increases the computational cost and complexity. Moreover, it could not tackle with viscoelastic fluids. The greatest challenge of doing this is setting proper boundary conditions on the free surface.In this paper, a novel numerical method is proposed for 2D viscoelastic extrusion flow simulations. To the best of the authors’ knowledge, this paper is the first implementation for combining ALE technique and FVM in studying time-dependent viscoelastic extrusion simulations, which extends the application fields of ALE method. The main contribution of this paper may be summarized as follows:(i)
We addressed the challenges in coupling the ALE with FVM method to simulate time-dependent viscoelastic extrusion flows. The obstacle of setting proper free surface boundary conditions is well addressed by deriving the Dirichlet pressure boundary condition and introducing a balance force imposed on the free surface. The proposed method is capable of dealing with both Newtonian fluids and viscoelastic fluids.
(ii)
We implemented a coupled ALE-FVM simulation platform based on OpenFOAM (Open Source Field Operation and Manipulation) [7]. The novel numerical method is validated by comprehensively comparing with previous literature. Compared to Tuković and Jasak’s work [6], our approach significantly reduced the storage and computation cost by abandoning the additional surface mesh. Moreover, based on OpenFOAM, our free surface numerical solver could be easily scaled to massive processor cores for large-scale parallel simulations.The rest of the paper is organized as follows. Details of the proposed method are presented in Section2. The numerical results and discussion are reported in Section 3. Conclusions are drawn in Section 4.


## 2. Methodology

In this section, the mathematical formulation of Oldroyd-B model and numerical algorithms used in the paper will be presented.


### 2.1. Mathematical Model

The governing equations for isothermal and incompressible viscoelastic fluid flows include the mass (continuity) conservation equation, momentum equation, and constitutive equation. The forms of those equations for the Oldroyd-B model are given by(1)
ρ
∂
u
∂
t
+
u
-
w
·
∇
u
=
-
∇
p
+
∇
·
τ
+
η
s
Δ
u
,
∇
·
u
=
0
,
τ
+
λ
τ
∇
=
2
η
p
D
,where ρ is the fluid density, u is the fluid velocity, w is the velocity of the cell, p is the pressure, and τ is the polymeric contribution for stress tensor. η
s and η
p are the solvent viscosity and polymer viscosity, respectively. λ is the characteristic relaxation time, and D is the rate of deformation tensor, which can be expressed as(2)
D
=
1
2
∇
u
+
∇
u
T
.The upper-convected time derivative of the polymer stress tensor is defined as(3)
τ
∇
=
∂
τ
∂
t
+
u
-
w
·
∇
τ
-
∇
u
T
·
τ
-
τ
·
∇
u
.For the convenience of comparison with others’ works, variables above are usually converted into their dimensionless form. New dimensionless variable Reynolds number Re and Weissenberg numberW
i are introduced by scaling the following variables as(4)
x
′
=
x
L
,
y
′
=
y
L
,
u
′
=
u
U
,
p
′
=
p
L
η
U
,
τ
′
=
τ
L
η
U
,
t
′
=
L
U
t
,
R
e
=
ρ
U
L
η
,
Wi
=
λ
U
L
,
β
=
η
s
ηin which L is the characteristic length (the half width of the die in this work), U is the characteristic velocity (the average velocity in this work), and η is the sum of the η
s and η
p; namely, total viscosity η
=
η
s
+
η
p. β is the viscosity ratio parameter.


### 2.2. Numerical Algorithm

#### 2.2.1. Implementation of Free Surface Boundary Condition

The most difficult problems in die-swell simulation are the implementation of free surface boundary condition. Boundary conditions imposed on the free surface could be divided into two categories, the kinematic boundary condition and the dynamic boundary condition. The kinematic boundary condition requires that the normal component of the fluid velocity equals the normal component of the interfacial grid velocity, as denoted by(5)
u
·
n
=
w
·
n
.The equation above means the relative speed between the fluid and the grid perpendicular to the interface is zero. In OpenFOAM, the velocity of fluid is stored at the cell center, which means it would be hard to get the velocity of interfacial grid directly. However, the flux of the fluid through each cell face could be obtained easily. Therefore, the velocity of interfacial face could be calculated according to the net mass flux on the surface. A detailed description about the implementation of the kinematic boundary condition can be found in [6, 8]. After moving the interfacial nodes at the interface, the internal nodes need to be moved according to the solution of solving a Laplace equation:(6)
∇
·
Γ
∇
u
=
0
,where Γ is the diffusion coefficient. The solver chosen to move nodes in this work is Laplace Face Decomposition [9], which could maintain a good quality of the mesh in large deformation cases. And the diffusion coefficient chosen is inversely proportional to the square of distance from the free surface. Regenerating the internal grids in this way could form lines of equal potential and decrease the nonorthogonality which in turn increases the quality of the mesh.The dynamic boundary condition is also called the traction boundary condition. Ignoring the surface tension, it means the traction exerted from the fluid to the atmosphere should be equal to the atmospheric forces acting on the fluid. Without loss of generality, the atmosphere pressure is set asp
0
=
0. So the boundary condition is given by(7)
σ
·
n
=
-
p
I
+
2
β
D
+
τ
·
n
=
-
p
I
+
T
·
n
=
0
,where σ is the total stress tensor, namely, the Cauchy stress tensor across the free surface, and T
=
2
β
D
+
τ is the extra stress tensor. When solving equations in OpenFOAM, PISO algorithm requires the term -
∇
p and the term η
s
Δ
u exists separately in the equation. However, the term σ is a combination of p, τ
, and D. Therefore, this boundary condition is hard to satisfy directly. It has been known that -
∇
p
+
∇
·
τ
+
η
s
Δ
u
=
∇
·
σ, and in FVM systems, every term in the equation should be integrated over the cell volume V. Therefore, the right side of (1) is actually solved as ∫
V
(
-
∇
p
+
∇
·
τ
+
η
s
Δ
u
)
d
V
=
∫
V
(
∇
·
σ
)
d
V. According to Gauss’s theorem, the spatial derivative term could be converted to integrals over the cell surface S bounding the volume. Then, we have(8)
∫
V
∇
·
σ
d
V
=
∫
S
d
S
·
σ
=
∑
a
σ
a
·
n
a
S
a
=
∑
a
≠
b
σ
a
·
n
a
S
a
+
σ
b
·
n
b
S
b
,where S is the face area vector and n is the unit normal vector of face area. a and b represent faces around a cell and free surface faces, respectively.The dynamic boundary condition requiresσ
b
·
n
b
=
0; however, it would not be the case if any of the boundary conditions is improper. Noting that the boundary condition would only influence the calculation of the cell that owns the boundary patch, we introduce a balance force σ
∗
=
-
σ on the free surface (it is zero elsewhere including the internal field and other boundary fields) to compensate the influence of nonzero boundary force. By adding the term ∇
·
σ
∗ into the original momentum equation, the right side of (1) can be solved as(9)
∫
V
-
∇
p
+
∇
·
τ
+
η
s
Δ
u
+
∇
·
σ
∗
d
V
=
∫
V
∇
·
σ
+
∇
·
σ
∗
d
V
=
∫
S
d
S
·
σ
+
σ
∗
=
∑
a
σ
a
·
n
a
S
a
-
σ
b
·
n
b
S
b
=
∑
a
≠
b
σ
a
·
n
a
S
a
.It can be seen that the influence of the total stress imposed on the free surface is offset in this way, which means the traction boundary condition is satisfied equivalently. The momentum equation is actually solved in the form:(10)
ρ
∂
u
∂
t
+
u
-
w
·
∇
u
=
-
∇
p
+
∇
·
τ
+
η
s
Δ
u
+
∇
·
σ
∗
Γ
f
,where Γ
f indicates the term is only applied on the free surface boundary. Since the SIMPLE algorithm is used to iteratively solve the linear momentum equation and the continuity equation, a proper pressure boundary condition is required during the step of solving the Poisson equation. According to the force balance in the normal direction on the interface, Dirichlet boundary condition for the pressure is given by(11)
p
=
n
·
T
·
n
.It is important to note that the pressure condition here is just a derived quantity, which has arisen only because of the need to solve the Poisson equation for p in isolation. While the traction boundary condition is the real physical quantity, it must be satisfied on the fluid. Provided that p and u are solved jointly, no pressure boundary condition is needed anymore [10].


#### 2.2.2. DEVSS Method for Numerical Stability

In order to enhance the stability of the linear momentum equation, especially for highW
i number flow, the discrete elastic split stress (DEVSS) numerical strategy proposed by Matallah et al. [11] is used. The momentum equation can be rewritten as(12)
ρ
∂
u
∂
t
+
u
-
w
·
∇
u
=
-
∇
p
+
∇
·
Σ
p
+
η
t
Δ
u
+
∇
·
σ
∗
Γ
f
,where η
t
=
η
s
+
α and ∇
·
Σ
p is the added artificial viscosity ratio. The term ∇
·
Σ
p
=
∇
·
τ
-
α
Δ
u is treated as a source term. A typical choice of α is η
p. This strategy greatly increases numerical stability and can be applied in modelling various constitutive models.


#### 2.2.3. The Simulation Procedure

Apart from the treatment of the free surface boundary conditions, the simulation procedure is similar to [6] as shown below:(1)
For a new time stept
n, the velocity field u
∗, the pressure field p
∗, the polymer stress field τ
∗, and the fluid mass fluxes m
˙
f are obtained from the previous time step. The initial net mass fluxes through the interface faces φ
n are calculated according to the previous interface position by subtracting the fluid mass fluxes by the volumetric face fluxes V
˙
f.
(2)
Move the interface mesh points so that the net mass flux through the interface could be compensated and the mass conservation law is satisfied.
(3)
According to the displacement of the interface mesh points, the position of internal mesh points is obtained by solving the Laplacian equation.
(4)
Beginning of the outer iteration loop:
(a)
getting the new estimated polymer stress tensorτ
∗
∗ by solving the specified constitutive equation with previous velocity field u
∗;
(b)
setting up the pressure boundary condition and estimating the balance forceσ
∗;
(c)
calculating the latest velocity fieldu
∗
∗ and the pressure filed p
∗
∗ with SIMPLE algorithm;
(d)
moving the interface mesh points similar to step (2) until the prescribed accuracy is reached.
(5)
If the programme has not reached the final time, return to step (1).Step (4)(d) is performed again in order to avoid the operation of step (3), which is a time-consuming step. Such procedure is reasonable provided that the displacement of the interface in one time step is smaller than the thickness of the nearest interface cell. And this condition could be easily satisfied by decreasing the time interval.


#### 2.2.4. Other General Numerical Schemes Used in the Simulation

In addition to the approaches described above, underrelaxation algorithm is also used to enhance the stability of the numerical simulation, where the relaxation parameters forp and U are 0.3 and 0.5, respectively.Backward Euler method with second order is used for time derivatives. In terms of the linear system solver, PCG (Preconditioned Conjugate Gradient method) with AMG (Algebraic Multi-Grid method) preconditioning is used for pressure, while BiCGstab (Bi-Conjugate Gradient Stabilized method) with Diagonal Incomplete Lower-Upper (DILU) preconditioning is used for velocity and stress [12, 13]. Absolute tolerances are 1.0
×
1
0
-
6 for the velocity and 1.0
×
1
0
-
7 for pressure and stress, respectively.


## 2.1. Mathematical Model

The governing equations for isothermal and incompressible viscoelastic fluid flows include the mass (continuity) conservation equation, momentum equation, and constitutive equation. The forms of those equations for the Oldroyd-B model are given by(1)
ρ
∂
u
∂
t
+
u
-
w
·
∇
u
=
-
∇
p
+
∇
·
τ
+
η
s
Δ
u
,
∇
·
u
=
0
,
τ
+
λ
τ
∇
=
2
η
p
D
,where ρ is the fluid density, u is the fluid velocity, w is the velocity of the cell, p is the pressure, and τ is the polymeric contribution for stress tensor. η
s and η
p are the solvent viscosity and polymer viscosity, respectively. λ is the characteristic relaxation time, and D is the rate of deformation tensor, which can be expressed as(2)
D
=
1
2
∇
u
+
∇
u
T
.The upper-convected time derivative of the polymer stress tensor is defined as(3)
τ
∇
=
∂
τ
∂
t
+
u
-
w
·
∇
τ
-
∇
u
T
·
τ
-
τ
·
∇
u
.For the convenience of comparison with others’ works, variables above are usually converted into their dimensionless form. New dimensionless variable Reynolds number Re and Weissenberg numberW
i are introduced by scaling the following variables as(4)
x
′
=
x
L
,
y
′
=
y
L
,
u
′
=
u
U
,
p
′
=
p
L
η
U
,
τ
′
=
τ
L
η
U
,
t
′
=
L
U
t
,
R
e
=
ρ
U
L
η
,
Wi
=
λ
U
L
,
β
=
η
s
ηin which L is the characteristic length (the half width of the die in this work), U is the characteristic velocity (the average velocity in this work), and η is the sum of the η
s and η
p; namely, total viscosity η
=
η
s
+
η
p. β is the viscosity ratio parameter.


## 2.2. Numerical Algorithm

### 2.2.1. Implementation of Free Surface Boundary Condition

The most difficult problems in die-swell simulation are the implementation of free surface boundary condition. Boundary conditions imposed on the free surface could be divided into two categories, the kinematic boundary condition and the dynamic boundary condition. The kinematic boundary condition requires that the normal component of the fluid velocity equals the normal component of the interfacial grid velocity, as denoted by(5)
u
·
n
=
w
·
n
.The equation above means the relative speed between the fluid and the grid perpendicular to the interface is zero. In OpenFOAM, the velocity of fluid is stored at the cell center, which means it would be hard to get the velocity of interfacial grid directly. However, the flux of the fluid through each cell face could be obtained easily. Therefore, the velocity of interfacial face could be calculated according to the net mass flux on the surface. A detailed description about the implementation of the kinematic boundary condition can be found in [6, 8]. After moving the interfacial nodes at the interface, the internal nodes need to be moved according to the solution of solving a Laplace equation:(6)
∇
·
Γ
∇
u
=
0
,where Γ is the diffusion coefficient. The solver chosen to move nodes in this work is Laplace Face Decomposition [9], which could maintain a good quality of the mesh in large deformation cases. And the diffusion coefficient chosen is inversely proportional to the square of distance from the free surface. Regenerating the internal grids in this way could form lines of equal potential and decrease the nonorthogonality which in turn increases the quality of the mesh.The dynamic boundary condition is also called the traction boundary condition. Ignoring the surface tension, it means the traction exerted from the fluid to the atmosphere should be equal to the atmospheric forces acting on the fluid. Without loss of generality, the atmosphere pressure is set asp
0
=
0. So the boundary condition is given by(7)
σ
·
n
=
-
p
I
+
2
β
D
+
τ
·
n
=
-
p
I
+
T
·
n
=
0
,where σ is the total stress tensor, namely, the Cauchy stress tensor across the free surface, and T
=
2
β
D
+
τ is the extra stress tensor. When solving equations in OpenFOAM, PISO algorithm requires the term -
∇
p and the term η
s
Δ
u exists separately in the equation. However, the term σ is a combination of p, τ
, and D. Therefore, this boundary condition is hard to satisfy directly. It has been known that -
∇
p
+
∇
·
τ
+
η
s
Δ
u
=
∇
·
σ, and in FVM systems, every term in the equation should be integrated over the cell volume V. Therefore, the right side of (1) is actually solved as ∫
V
(
-
∇
p
+
∇
·
τ
+
η
s
Δ
u
)
d
V
=
∫
V
(
∇
·
σ
)
d
V. According to Gauss’s theorem, the spatial derivative term could be converted to integrals over the cell surface S bounding the volume. Then, we have(8)
∫
V
∇
·
σ
d
V
=
∫
S
d
S
·
σ
=
∑
a
σ
a
·
n
a
S
a
=
∑
a
≠
b
σ
a
·
n
a
S
a
+
σ
b
·
n
b
S
b
,where S is the face area vector and n is the unit normal vector of face area. a and b represent faces around a cell and free surface faces, respectively.The dynamic boundary condition requiresσ
b
·
n
b
=
0; however, it would not be the case if any of the boundary conditions is improper. Noting that the boundary condition would only influence the calculation of the cell that owns the boundary patch, we introduce a balance force σ
∗
=
-
σ on the free surface (it is zero elsewhere including the internal field and other boundary fields) to compensate the influence of nonzero boundary force. By adding the term ∇
·
σ
∗ into the original momentum equation, the right side of (1) can be solved as(9)
∫
V
-
∇
p
+
∇
·
τ
+
η
s
Δ
u
+
∇
·
σ
∗
d
V
=
∫
V
∇
·
σ
+
∇
·
σ
∗
d
V
=
∫
S
d
S
·
σ
+
σ
∗
=
∑
a
σ
a
·
n
a
S
a
-
σ
b
·
n
b
S
b
=
∑
a
≠
b
σ
a
·
n
a
S
a
.It can be seen that the influence of the total stress imposed on the free surface is offset in this way, which means the traction boundary condition is satisfied equivalently. The momentum equation is actually solved in the form:(10)
ρ
∂
u
∂
t
+
u
-
w
·
∇
u
=
-
∇
p
+
∇
·
τ
+
η
s
Δ
u
+
∇
·
σ
∗
Γ
f
,where Γ
f indicates the term is only applied on the free surface boundary. Since the SIMPLE algorithm is used to iteratively solve the linear momentum equation and the continuity equation, a proper pressure boundary condition is required during the step of solving the Poisson equation. According to the force balance in the normal direction on the interface, Dirichlet boundary condition for the pressure is given by(11)
p
=
n
·
T
·
n
.It is important to note that the pressure condition here is just a derived quantity, which has arisen only because of the need to solve the Poisson equation for p in isolation. While the traction boundary condition is the real physical quantity, it must be satisfied on the fluid. Provided that p and u are solved jointly, no pressure boundary condition is needed anymore [10].


### 2.2.2. DEVSS Method for Numerical Stability

In order to enhance the stability of the linear momentum equation, especially for highW
i number flow, the discrete elastic split stress (DEVSS) numerical strategy proposed by Matallah et al. [11] is used. The momentum equation can be rewritten as(12)
ρ
∂
u
∂
t
+
u
-
w
·
∇
u
=
-
∇
p
+
∇
·
Σ
p
+
η
t
Δ
u
+
∇
·
σ
∗
Γ
f
,where η
t
=
η
s
+
α and ∇
·
Σ
p is the added artificial viscosity ratio. The term ∇
·
Σ
p
=
∇
·
τ
-
α
Δ
u is treated as a source term. A typical choice of α is η
p. This strategy greatly increases numerical stability and can be applied in modelling various constitutive models.


### 2.2.3. The Simulation Procedure

Apart from the treatment of the free surface boundary conditions, the simulation procedure is similar to [6] as shown below:(1)
For a new time stept
n, the velocity field u
∗, the pressure field p
∗, the polymer stress field τ
∗, and the fluid mass fluxes m
˙
f are obtained from the previous time step. The initial net mass fluxes through the interface faces φ
n are calculated according to the previous interface position by subtracting the fluid mass fluxes by the volumetric face fluxes V
˙
f.
(2)
Move the interface mesh points so that the net mass flux through the interface could be compensated and the mass conservation law is satisfied.
(3)
According to the displacement of the interface mesh points, the position of internal mesh points is obtained by solving the Laplacian equation.
(4)
Beginning of the outer iteration loop:
(a)
getting the new estimated polymer stress tensorτ
∗
∗ by solving the specified constitutive equation with previous velocity field u
∗;
(b)
setting up the pressure boundary condition and estimating the balance forceσ
∗;
(c)
calculating the latest velocity fieldu
∗
∗ and the pressure filed p
∗
∗ with SIMPLE algorithm;
(d)
moving the interface mesh points similar to step (2) until the prescribed accuracy is reached.
(5)
If the programme has not reached the final time, return to step (1).Step (4)(d) is performed again in order to avoid the operation of step (3), which is a time-consuming step. Such procedure is reasonable provided that the displacement of the interface in one time step is smaller than the thickness of the nearest interface cell. And this condition could be easily satisfied by decreasing the time interval.


### 2.2.4. Other General Numerical Schemes Used in the Simulation

In addition to the approaches described above, underrelaxation algorithm is also used to enhance the stability of the numerical simulation, where the relaxation parameters forp and U are 0.3 and 0.5, respectively.Backward Euler method with second order is used for time derivatives. In terms of the linear system solver, PCG (Preconditioned Conjugate Gradient method) with AMG (Algebraic Multi-Grid method) preconditioning is used for pressure, while BiCGstab (Bi-Conjugate Gradient Stabilized method) with Diagonal Incomplete Lower-Upper (DILU) preconditioning is used for velocity and stress [12, 13]. Absolute tolerances are 1.0
×
1
0
-
6 for the velocity and 1.0
×
1
0
-
7 for pressure and stress, respectively.


## 2.2.1. Implementation of Free Surface Boundary Condition

The most difficult problems in die-swell simulation are the implementation of free surface boundary condition. Boundary conditions imposed on the free surface could be divided into two categories, the kinematic boundary condition and the dynamic boundary condition. The kinematic boundary condition requires that the normal component of the fluid velocity equals the normal component of the interfacial grid velocity, as denoted by(5)
u
·
n
=
w
·
n
.The equation above means the relative speed between the fluid and the grid perpendicular to the interface is zero. In OpenFOAM, the velocity of fluid is stored at the cell center, which means it would be hard to get the velocity of interfacial grid directly. However, the flux of the fluid through each cell face could be obtained easily. Therefore, the velocity of interfacial face could be calculated according to the net mass flux on the surface. A detailed description about the implementation of the kinematic boundary condition can be found in [6, 8]. After moving the interfacial nodes at the interface, the internal nodes need to be moved according to the solution of solving a Laplace equation:(6)
∇
·
Γ
∇
u
=
0
,where Γ is the diffusion coefficient. The solver chosen to move nodes in this work is Laplace Face Decomposition [9], which could maintain a good quality of the mesh in large deformation cases. And the diffusion coefficient chosen is inversely proportional to the square of distance from the free surface. Regenerating the internal grids in this way could form lines of equal potential and decrease the nonorthogonality which in turn increases the quality of the mesh.The dynamic boundary condition is also called the traction boundary condition. Ignoring the surface tension, it means the traction exerted from the fluid to the atmosphere should be equal to the atmospheric forces acting on the fluid. Without loss of generality, the atmosphere pressure is set asp
0
=
0. So the boundary condition is given by(7)
σ
·
n
=
-
p
I
+
2
β
D
+
τ
·
n
=
-
p
I
+
T
·
n
=
0
,where σ is the total stress tensor, namely, the Cauchy stress tensor across the free surface, and T
=
2
β
D
+
τ is the extra stress tensor. When solving equations in OpenFOAM, PISO algorithm requires the term -
∇
p and the term η
s
Δ
u exists separately in the equation. However, the term σ is a combination of p, τ
, and D. Therefore, this boundary condition is hard to satisfy directly. It has been known that -
∇
p
+
∇
·
τ
+
η
s
Δ
u
=
∇
·
σ, and in FVM systems, every term in the equation should be integrated over the cell volume V. Therefore, the right side of (1) is actually solved as ∫
V
(
-
∇
p
+
∇
·
τ
+
η
s
Δ
u
)
d
V
=
∫
V
(
∇
·
σ
)
d
V. According to Gauss’s theorem, the spatial derivative term could be converted to integrals over the cell surface S bounding the volume. Then, we have(8)
∫
V
∇
·
σ
d
V
=
∫
S
d
S
·
σ
=
∑
a
σ
a
·
n
a
S
a
=
∑
a
≠
b
σ
a
·
n
a
S
a
+
σ
b
·
n
b
S
b
,where S is the face area vector and n is the unit normal vector of face area. a and b represent faces around a cell and free surface faces, respectively.The dynamic boundary condition requiresσ
b
·
n
b
=
0; however, it would not be the case if any of the boundary conditions is improper. Noting that the boundary condition would only influence the calculation of the cell that owns the boundary patch, we introduce a balance force σ
∗
=
-
σ on the free surface (it is zero elsewhere including the internal field and other boundary fields) to compensate the influence of nonzero boundary force. By adding the term ∇
·
σ
∗ into the original momentum equation, the right side of (1) can be solved as(9)
∫
V
-
∇
p
+
∇
·
τ
+
η
s
Δ
u
+
∇
·
σ
∗
d
V
=
∫
V
∇
·
σ
+
∇
·
σ
∗
d
V
=
∫
S
d
S
·
σ
+
σ
∗
=
∑
a
σ
a
·
n
a
S
a
-
σ
b
·
n
b
S
b
=
∑
a
≠
b
σ
a
·
n
a
S
a
.It can be seen that the influence of the total stress imposed on the free surface is offset in this way, which means the traction boundary condition is satisfied equivalently. The momentum equation is actually solved in the form:(10)
ρ
∂
u
∂
t
+
u
-
w
·
∇
u
=
-
∇
p
+
∇
·
τ
+
η
s
Δ
u
+
∇
·
σ
∗
Γ
f
,where Γ
f indicates the term is only applied on the free surface boundary. Since the SIMPLE algorithm is used to iteratively solve the linear momentum equation and the continuity equation, a proper pressure boundary condition is required during the step of solving the Poisson equation. According to the force balance in the normal direction on the interface, Dirichlet boundary condition for the pressure is given by(11)
p
=
n
·
T
·
n
.It is important to note that the pressure condition here is just a derived quantity, which has arisen only because of the need to solve the Poisson equation for p in isolation. While the traction boundary condition is the real physical quantity, it must be satisfied on the fluid. Provided that p and u are solved jointly, no pressure boundary condition is needed anymore [10].


## 2.2.2. DEVSS Method for Numerical Stability

In order to enhance the stability of the linear momentum equation, especially for highW
i number flow, the discrete elastic split stress (DEVSS) numerical strategy proposed by Matallah et al. [11] is used. The momentum equation can be rewritten as(12)
ρ
∂
u
∂
t
+
u
-
w
·
∇
u
=
-
∇
p
+
∇
·
Σ
p
+
η
t
Δ
u
+
∇
·
σ
∗
Γ
f
,where η
t
=
η
s
+
α and ∇
·
Σ
p is the added artificial viscosity ratio. The term ∇
·
Σ
p
=
∇
·
τ
-
α
Δ
u is treated as a source term. A typical choice of α is η
p. This strategy greatly increases numerical stability and can be applied in modelling various constitutive models.


## 2.2.3. The Simulation Procedure

Apart from the treatment of the free surface boundary conditions, the simulation procedure is similar to [6] as shown below:(1)
For a new time stept
n, the velocity field u
∗, the pressure field p
∗, the polymer stress field τ
∗, and the fluid mass fluxes m
˙
f are obtained from the previous time step. The initial net mass fluxes through the interface faces φ
n are calculated according to the previous interface position by subtracting the fluid mass fluxes by the volumetric face fluxes V
˙
f.
(2)
Move the interface mesh points so that the net mass flux through the interface could be compensated and the mass conservation law is satisfied.
(3)
According to the displacement of the interface mesh points, the position of internal mesh points is obtained by solving the Laplacian equation.
(4)
Beginning of the outer iteration loop:
(a)
getting the new estimated polymer stress tensorτ
∗
∗ by solving the specified constitutive equation with previous velocity field u
∗;
(b)
setting up the pressure boundary condition and estimating the balance forceσ
∗;
(c)
calculating the latest velocity fieldu
∗
∗ and the pressure filed p
∗
∗ with SIMPLE algorithm;
(d)
moving the interface mesh points similar to step (2) until the prescribed accuracy is reached.
(5)
If the programme has not reached the final time, return to step (1).Step (4)(d) is performed again in order to avoid the operation of step (3), which is a time-consuming step. Such procedure is reasonable provided that the displacement of the interface in one time step is smaller than the thickness of the nearest interface cell. And this condition could be easily satisfied by decreasing the time interval.


## 2.2.4. Other General Numerical Schemes Used in the Simulation

In addition to the approaches described above, underrelaxation algorithm is also used to enhance the stability of the numerical simulation, where the relaxation parameters forp and U are 0.3 and 0.5, respectively.Backward Euler method with second order is used for time derivatives. In terms of the linear system solver, PCG (Preconditioned Conjugate Gradient method) with AMG (Algebraic Multi-Grid method) preconditioning is used for pressure, while BiCGstab (Bi-Conjugate Gradient Stabilized method) with Diagonal Incomplete Lower-Upper (DILU) preconditioning is used for velocity and stress [12, 13]. Absolute tolerances are 1.0
×
1
0
-
6 for the velocity and 1.0
×
1
0
-
7 for pressure and stress, respectively.


## 3. Results and Discussion

### 3.1. Basic Geometric Shape and Boundary Conditions

A sketch of the geometry and the boundary conditions for modelling the Newtonian and Oldroyd-B fluid is shown in Figure1. The half width of the die L is chosen as the characteristic length. The lengths of the upstream and the downstream are set as L
1 and L
2, respectively. They should be long enough to ensure the flow before entry of the die exit and further downstream at the outlet is fully developed. In this paper, length settings are L
=
1 and L
1
=
L
2
=
8.Figure 1
Die-swell flow configurations.As the pressure, the velocity, and the polymer stress are solved in a coupled way in our system, we usually give either all the correct boundary conditions or any one of them on the boundary. In order to avoid unnecessary numerical failure near the inlet, the conditions of fully developed flow are imposed on the inlet boundary. The velocity alongx direction at the inlet is u
i
n
=
1.5
(
1
-
y
2
), while the velocity along y direction is zero. Therefore, the average velocity of the inflow is u
¯
=
1, and the shear rate near the wall is γ
˙
w
=
3. The polymer stress tensor at the inlet for the Oldroyd-B model is set according to the analytical solution. No-slip boundary condition is set at the wall. The symmetry boundary condition is set at the centreline for all fields. Symmetry boundary condition is provided by OpenFoam itself. It means there is neither convection flux nor diffusion flux across the plane; that is, the vertical velocity and the shear stress component are zero (u
x
=
0 and σ
x
y
=
0). We assumed that the flow at the outlet has been already fully developed; therefore, the pressure is set as zero at the outlet. Note that the pressure here is just a reference pressure which is used to solve the governing equations. In terms of the free surface, the pressure is initially set as fixed value zero, and it would change dynamically during the simulation (see (11)). Without special declaration, zero gradient boundary condition is set for the remaining boundaries.


### 3.2. Validation of Computational Model with Newtonian Fluid

By settingR
e
=
0.01, W
i
=
0, β
=
1, and ɛ
=
0, Oldroyd-B model degrades to Newtonian fluid. Under isothermal condition and neglecting the effects of gravity and surface tension, the swelling ratio, defined as the maximum diameter of the extrudate divided by the diameter of the die, is 1.193 from our simulation. It is in an excellent agreement with the results of Mitsoulis et al. [14] and Georgiou and Boudouvis [15], which is around 1.19.


### 3.3. Validation of Computational Model with Viscoelastic Fluid Using Oldroyd-B Model

Based on K-BKZ models, Tanner [16] derived an analytical formula to account for the viscous and elastic contributions to the swell ratio as(13)
S
R
=
0.19
+
1
+
1
12
N
1
σ
x
y
w
2
1
/
4
,where N
1
=
σ
x
x
-
σ
y
y is the first normal stress difference and σ
x
y is the shear stress. Both of them are evaluated at the wall. Introducing the parameter of recoverable shear R
s
=
(
N
1
/
2
σ
x
y
)
w, the swelling ratio can be estimated as(14)
S
R
=
0.19
+
1
+
1
3
R
s
2
1
/
4
.For the Oldroyd-B model, the recoverable shear of a fully developed Poiseuille flow is equal to(15)
R
s
=
2
Wi
1
-
β
γ
˙
2
2
γ
˙
w
=
1
-
β
Wi
γ
˙
w
.Since in our caseγ
˙
w
=
3 and β
=
0.1, the recoverable shear is R
s
=
2.7
W
i.Mesh convergence was tested atW
i
=
0.1 with seven meshes, which are listed in Table 1. Δ
x
m
i
n
/
L and Δ
y
m
i
n
/
L are the relative minimum length along x-axis and y-axis, respectively. These cells appear near the exit of the die. Through the comparison of all these seven meshes, it is found that the value of Δ
x
m
i
n
/
L plays a dominant role when determining the welling ratio. Furthermore, if Δ
y
m
i
n
/
L
≤
0.01 and Δ
x
m
i
n
/
L
≤
0.025, the relative error of the swelling ratio is below 0.2%. The solution convergence is reached by using M6. The time interval is Δ
t
=
1
0
-
6 for M7 and Δ
t
=
1
0
-
5 for other cases (except for M1).Table 1
Mesh parameters for convergence validation with Wi = 0.1.
Name
Size
Δ
x
m
i
n
/
L
Δ
y
m
i
n
/
L
S
R
M1
160
∗
30
0.1
0.02
1.1677
M2
320
∗
30
0.05
0.02
1.1727
M3
160
∗
60
0.1
0.01
1.1691
M4
640
∗
30
0.025
0.02
1.1789
M5
160
∗
120
0.1
0.005
1.1697
M6
360
∗
60
0.025
0.01
1.1770
M7
400
∗
60
0.0125
0.01
1.1787The value of pressure and polymer stress near the singularity is illustrated in Figure2, where x is in the range of [
-
0.4,0.4
], along with different mesh densities (M4, M6) at W
i
=
0.1. It is shown that the two curves are very close to each other, which indicates that the given mesh density is sufficient. In terms of the tendency, the pressure gradually declines at first and then experiences a steep fall followed by a mild elevation before the exit. Subsequently, it dramatically drops from positive to negative and approaches zero afterwards. As shown in Figures 2(b) and 2(c), there is a highest point for τ
x
x and lowest point for τ
y
y, where both points are around the exit. Furthermore, the curve of τ
y
y is sharper than the one of τ
x
x. As x increases, the value of τ
x
y decreases at first, reaches the bottom, and sharply grows to an apogee, before it declines to zero gradually. In terms of the trend of the profiles as well as the maximum and the minimum values near the singularity, the results of the proposed method are close to those of the spectral element method with 3rd order polynomial adopted in Claus work [17], which again demonstrates that the current mesh is capable of providing reliable results. The small deviation in terms of the positions of the maximum value and the minimum value are attributed to the fact that the velocity is stored on the cell center in this paper and on the cell vertex in [17], respectively.Figure 2
Pressure and components of polymer stress near the singularity atW
i
=
0.1.
(a)
p
(b)
τ
x
x
(c)
τ
y
y
(d)
τ
x
yIn general, with the increase ofW
i, the profile of components of the polymer stress should be similar but with a greater value. However, the figure for τ
y
y is much different from the intuition. Figure 3 shows that the nadir for W
i
=
0.1 and W
i
=
0.2 which is negative becomes the peak for W
i
=
0.5 which is positive. The force along y direction caused by the polymer stress components τ
y
y is n
y
τ
y
y, where n
y is vertical component of the standard interfacial normal. Thus, the more likely n
y is approaching to one (corresponding to a small swelling ratio with a gentle curve), the greater the influence of τ
y
y would be. As the force increases, the extrusion is growing larger.Figure 3
τ
y
y near the singularity for various W
i.The swelling ratio is given in Figure4 with a comparison with others’ works. Note that, in their works, the definition of relative parameters could be different. For example, Russo and Phillips [18] and Tomé et al. [19] used the full width as the characteristic length. We just altered them to be the same one in our work. Moreover, Crochet and Keunings [20] and Tomé et al. [19] took the Reynolds number as 0 and 0.5, respectively. Russo and Phillips [18] tested two situations with R
e
=
0.001 and R
e
=
0.5 and the results shed light on the fact that no significant change was achieved by varying the Reynolds number when Re is small enough. All of these three papers took β
=
1
/
9 in their similar simulations. Our simulation is kind of similar to Crochet and Keunings [20] and Tomé et al. [19] who varied the relaxation time λ in different cases; therefore, the trend of the curves is the same. However, Russo and Phillips [18] changed the input velocity, that is, both shear rate near the wall γ
˙
w and the Weissenberg number W
i are changed. In his case, the recoverable shear rate is far different from the one in our simulation even with the same W
i. That is the reason the trend of his curve seems different from others. The swelling ratio of all these curves shown increases with the increase of the Weissenberg number. From Figure 4, our results are above Crochet and Keunings [20] and Tomé et al. [19] but below Russo and Phillips [18]. With the increase of W
i number, the elastic effect plays a dominant role in determining the die-swell. In addition, slight drop of swelling ratio at low W
i number illustrated in Figure 4 is due to the negative value of τ
y
y and was also observed in [20].Figure 4
Swelling ratio with the increase ofW
i.The steady state profiles of the free surface at variousW
i numbers are shown in Figure 5. It can be seen that a larger W
i value would result in an increased slope of the curve. This is because the first normal stress at the exit of the channel for bigger W
i number is larger than the smaller ones. Such phenomenon also indicates the main factor of the swelling is the value of the first normal stress near the exit of the die. Under higher W
i number, the hexahedron cell is deformed greatly which increases the mesh nonorthogonality. The maximum mesh nonorthogonality has reached as high as 40 at W
i
=
1.0 and the threshold in OpenFOAM is usually set as 70. Though the nonorthogonality value could be reduced by mesh refinement, the computational time would increase. The maximum swelling ratio obtained from our simulation is about 2.26 at W
i
=
2.5 using Mesh M1.Figure 5
Free surface profile for variousW
i for β
=
0.1.Contour plots of velocity componentsU
x for W
i
=
0.1, 0.5, 0.8, 1.0 and U
y for W
i
=
0.5, 1.0 are displayed in Figure 6. By increasing W
i number, the horizontal components of the velocity U
x increase along y directions. The contour line of U
x
=
1.4 expands upward and forward. The vertical components of the velocity U
y also increase due to larger elastic effects which results in a larger swelling ratio.Figure 6
Contour plots of velocity componentsU
x (a)–(d) and U
y (e)-(f) for a range of W
i.
(a)
W
i
=
0.1
(b)
W
i
=
0.5
(c)
W
i
=
0.8
(d)
W
i
=
1
(e)
W
i
=
0.5
(f)
W
i
=
1The profile of the first normal stress differenceN
1 just before the exit (at x
=
-
0.04, and the exit is at x
=
0) of the die for a large range of W
i is shown in Figure 7. These lines are sampled along y-axis at fixed x. Apparently, the maximum value of each line increases as W
i becomes larger, which is consistent with (15). All of these four lines begin with zero at the middle of the die (y
=
0). As for W
i
=
1.0, it goes up gradually until it approaches the wall, where it grows instantly to the peak. Subsequently, it drops dramatically. The profiles of other W
i values smaller than 1.0 are something like shifting the curve of W
i
=
1.0 downwards to the right. In detail, there is no peak in the case of W
i
=
0.1. Generally speaking, a larger W
i results in a higher peak, which indicates the first normal stress dominates the extrusion process.Figure 7
First normal stress difference atx
=
-
0.04 for various W
i with β
=
0.1.Decreasing the viscosity ratioβ from 1 to 0 means gradually increasing the contribution of elastic stress and decreasing the effects of purely viscous stress. When β reaches zero, the model degrades to the upper-convected Maxwell (UCM) model. In order to investigate the role of β in die-swell simulations, simulations under β
=
0.5 and β
=
0.9 are carried out. Figure 8 illustrates the trend that the swelling ratio would increase as W
i number increases, regardless of the value of β. Moreover, the result departs from Tanners theory obviously. Considering an extreme case, as β approaches one, the swelling ratio almost remains a constant even for a large value of W
i according to Tanners theory. However, this does not match the experimental result. For the case β
=
0.9, W
i
=
1 and the case β
=
0.1, W
i
=
1, the recoverable shear R
s are 3 and 2.7, respectively. According to (15), the swelling ratio of the former one should be larger than the later one, while the fact is contrary. Russo and Phillips [18] attribute the discrepancy to the fact that Tanners formula derived from a KBK-Z integral constitutive equation only considered a single relaxation time; however, the Oldroyd-B model owns a retardation time in addition to the relaxation time. This result might be related to the fact that the elongational response for β
=
0.1 is much stronger than the one for β
=
0.9 as illustrated in Figure 9. Thus, the extensional response also plays an important role in swelling process and should be taken into account when predicting the final swelling ratio. The simulation results in [21] demonstrate that flows with a sufficiently high elongational viscosity are able to swell to an apparently large state even at a relatively low flow rate, ignoring their performance in simple shear flows. Nevertheless, decreasing the viscosity ratio would lead to an increase in the recoverable shear and, in turn, increases the value of swelling ration according to (15). By comparing the value in Figures 4 and 8, it is reasonable to say that Tanners theory still makes some sense. Apart from that, with a careful comparison, it is found that when the viscosity ratio β is very large (for the case of β
=
0.9), which represents a very dilute fluid, the swelling ratio suffers a significant drop. In comparison, when β is relatively small (for the case of β
=
0.5), the swelling ratio is much less sensitive. This phenomena was also observed in [22]. Note that because of different definitions of the viscosity ratio, β in this paper is equivalent to 1
-
β in their works.Figure 8
Swelling ratio with the increase ofW
i for β
=
0.5 and β
=
0.9.Figure 9
Planar extensional viscosity for Oldroyd-B model with variousβ.As can be seen from (13), the value of the shear stress would not change once the shear rate is determined. Consequently, the first normal stress difference is the only variable which determines the swelling ratio. Figures 7 and 10 demonstrate that a high proportion of polymeric stress contribution would result in a larger value of N
1, which would provide greater power to swell. According to the analytical solution of the Oldroyd-B model, smaller β would result in a higher value of the first normal stress. Therefore, the first normal stress plays a vital role in determining the behaviour of extrusion flows and could be used to predict the final results of swelling ratio.Figure 10
First normal stress difference atx
=
-
0.04 for various W
i with β
=
0.5 and β
=
0.9.
(a)
β
=
0.5
(b)
β
=
0.9


## 3.1. Basic Geometric Shape and Boundary Conditions

A sketch of the geometry and the boundary conditions for modelling the Newtonian and Oldroyd-B fluid is shown in Figure1. The half width of the die L is chosen as the characteristic length. The lengths of the upstream and the downstream are set as L
1 and L
2, respectively. They should be long enough to ensure the flow before entry of the die exit and further downstream at the outlet is fully developed. In this paper, length settings are L
=
1 and L
1
=
L
2
=
8.Figure 1
Die-swell flow configurations.As the pressure, the velocity, and the polymer stress are solved in a coupled way in our system, we usually give either all the correct boundary conditions or any one of them on the boundary. In order to avoid unnecessary numerical failure near the inlet, the conditions of fully developed flow are imposed on the inlet boundary. The velocity alongx direction at the inlet is u
i
n
=
1.5
(
1
-
y
2
), while the velocity along y direction is zero. Therefore, the average velocity of the inflow is u
¯
=
1, and the shear rate near the wall is γ
˙
w
=
3. The polymer stress tensor at the inlet for the Oldroyd-B model is set according to the analytical solution. No-slip boundary condition is set at the wall. The symmetry boundary condition is set at the centreline for all fields. Symmetry boundary condition is provided by OpenFoam itself. It means there is neither convection flux nor diffusion flux across the plane; that is, the vertical velocity and the shear stress component are zero (u
x
=
0 and σ
x
y
=
0). We assumed that the flow at the outlet has been already fully developed; therefore, the pressure is set as zero at the outlet. Note that the pressure here is just a reference pressure which is used to solve the governing equations. In terms of the free surface, the pressure is initially set as fixed value zero, and it would change dynamically during the simulation (see (11)). Without special declaration, zero gradient boundary condition is set for the remaining boundaries.


## 3.2. Validation of Computational Model with Newtonian Fluid

By settingR
e
=
0.01, W
i
=
0, β
=
1, and ɛ
=
0, Oldroyd-B model degrades to Newtonian fluid. Under isothermal condition and neglecting the effects of gravity and surface tension, the swelling ratio, defined as the maximum diameter of the extrudate divided by the diameter of the die, is 1.193 from our simulation. It is in an excellent agreement with the results of Mitsoulis et al. [14] and Georgiou and Boudouvis [15], which is around 1.19.


## 3.3. Validation of Computational Model with Viscoelastic Fluid Using Oldroyd-B Model

Based on K-BKZ models, Tanner [16] derived an analytical formula to account for the viscous and elastic contributions to the swell ratio as(13)
S
R
=
0.19
+
1
+
1
12
N
1
σ
x
y
w
2
1
/
4
,where N
1
=
σ
x
x
-
σ
y
y is the first normal stress difference and σ
x
y is the shear stress. Both of them are evaluated at the wall. Introducing the parameter of recoverable shear R
s
=
(
N
1
/
2
σ
x
y
)
w, the swelling ratio can be estimated as(14)
S
R
=
0.19
+
1
+
1
3
R
s
2
1
/
4
.For the Oldroyd-B model, the recoverable shear of a fully developed Poiseuille flow is equal to(15)
R
s
=
2
Wi
1
-
β
γ
˙
2
2
γ
˙
w
=
1
-
β
Wi
γ
˙
w
.Since in our caseγ
˙
w
=
3 and β
=
0.1, the recoverable shear is R
s
=
2.7
W
i.Mesh convergence was tested atW
i
=
0.1 with seven meshes, which are listed in Table 1. Δ
x
m
i
n
/
L and Δ
y
m
i
n
/
L are the relative minimum length along x-axis and y-axis, respectively. These cells appear near the exit of the die. Through the comparison of all these seven meshes, it is found that the value of Δ
x
m
i
n
/
L plays a dominant role when determining the welling ratio. Furthermore, if Δ
y
m
i
n
/
L
≤
0.01 and Δ
x
m
i
n
/
L
≤
0.025, the relative error of the swelling ratio is below 0.2%. The solution convergence is reached by using M6. The time interval is Δ
t
=
1
0
-
6 for M7 and Δ
t
=
1
0
-
5 for other cases (except for M1).Table 1
Mesh parameters for convergence validation with Wi = 0.1.
Name
Size
Δ
x
m
i
n
/
L
Δ
y
m
i
n
/
L
S
R
M1
160
∗
30
0.1
0.02
1.1677
M2
320
∗
30
0.05
0.02
1.1727
M3
160
∗
60
0.1
0.01
1.1691
M4
640
∗
30
0.025
0.02
1.1789
M5
160
∗
120
0.1
0.005
1.1697
M6
360
∗
60
0.025
0.01
1.1770
M7
400
∗
60
0.0125
0.01
1.1787The value of pressure and polymer stress near the singularity is illustrated in Figure2, where x is in the range of [
-
0.4,0.4
], along with different mesh densities (M4, M6) at W
i
=
0.1. It is shown that the two curves are very close to each other, which indicates that the given mesh density is sufficient. In terms of the tendency, the pressure gradually declines at first and then experiences a steep fall followed by a mild elevation before the exit. Subsequently, it dramatically drops from positive to negative and approaches zero afterwards. As shown in Figures 2(b) and 2(c), there is a highest point for τ
x
x and lowest point for τ
y
y, where both points are around the exit. Furthermore, the curve of τ
y
y is sharper than the one of τ
x
x. As x increases, the value of τ
x
y decreases at first, reaches the bottom, and sharply grows to an apogee, before it declines to zero gradually. In terms of the trend of the profiles as well as the maximum and the minimum values near the singularity, the results of the proposed method are close to those of the spectral element method with 3rd order polynomial adopted in Claus work [17], which again demonstrates that the current mesh is capable of providing reliable results. The small deviation in terms of the positions of the maximum value and the minimum value are attributed to the fact that the velocity is stored on the cell center in this paper and on the cell vertex in [17], respectively.Figure 2
Pressure and components of polymer stress near the singularity atW
i
=
0.1.
(a)
p
(b)
τ
x
x
(c)
τ
y
y
(d)
τ
x
yIn general, with the increase ofW
i, the profile of components of the polymer stress should be similar but with a greater value. However, the figure for τ
y
y is much different from the intuition. Figure 3 shows that the nadir for W
i
=
0.1 and W
i
=
0.2 which is negative becomes the peak for W
i
=
0.5 which is positive. The force along y direction caused by the polymer stress components τ
y
y is n
y
τ
y
y, where n
y is vertical component of the standard interfacial normal. Thus, the more likely n
y is approaching to one (corresponding to a small swelling ratio with a gentle curve), the greater the influence of τ
y
y would be. As the force increases, the extrusion is growing larger.Figure 3
τ
y
y near the singularity for various W
i.The swelling ratio is given in Figure4 with a comparison with others’ works. Note that, in their works, the definition of relative parameters could be different. For example, Russo and Phillips [18] and Tomé et al. [19] used the full width as the characteristic length. We just altered them to be the same one in our work. Moreover, Crochet and Keunings [20] and Tomé et al. [19] took the Reynolds number as 0 and 0.5, respectively. Russo and Phillips [18] tested two situations with R
e
=
0.001 and R
e
=
0.5 and the results shed light on the fact that no significant change was achieved by varying the Reynolds number when Re is small enough. All of these three papers took β
=
1
/
9 in their similar simulations. Our simulation is kind of similar to Crochet and Keunings [20] and Tomé et al. [19] who varied the relaxation time λ in different cases; therefore, the trend of the curves is the same. However, Russo and Phillips [18] changed the input velocity, that is, both shear rate near the wall γ
˙
w and the Weissenberg number W
i are changed. In his case, the recoverable shear rate is far different from the one in our simulation even with the same W
i. That is the reason the trend of his curve seems different from others. The swelling ratio of all these curves shown increases with the increase of the Weissenberg number. From Figure 4, our results are above Crochet and Keunings [20] and Tomé et al. [19] but below Russo and Phillips [18]. With the increase of W
i number, the elastic effect plays a dominant role in determining the die-swell. In addition, slight drop of swelling ratio at low W
i number illustrated in Figure 4 is due to the negative value of τ
y
y and was also observed in [20].Figure 4
Swelling ratio with the increase ofW
i.The steady state profiles of the free surface at variousW
i numbers are shown in Figure 5. It can be seen that a larger W
i value would result in an increased slope of the curve. This is because the first normal stress at the exit of the channel for bigger W
i number is larger than the smaller ones. Such phenomenon also indicates the main factor of the swelling is the value of the first normal stress near the exit of the die. Under higher W
i number, the hexahedron cell is deformed greatly which increases the mesh nonorthogonality. The maximum mesh nonorthogonality has reached as high as 40 at W
i
=
1.0 and the threshold in OpenFOAM is usually set as 70. Though the nonorthogonality value could be reduced by mesh refinement, the computational time would increase. The maximum swelling ratio obtained from our simulation is about 2.26 at W
i
=
2.5 using Mesh M1.Figure 5
Free surface profile for variousW
i for β
=
0.1.Contour plots of velocity componentsU
x for W
i
=
0.1, 0.5, 0.8, 1.0 and U
y for W
i
=
0.5, 1.0 are displayed in Figure 6. By increasing W
i number, the horizontal components of the velocity U
x increase along y directions. The contour line of U
x
=
1.4 expands upward and forward. The vertical components of the velocity U
y also increase due to larger elastic effects which results in a larger swelling ratio.Figure 6
Contour plots of velocity componentsU
x (a)–(d) and U
y (e)-(f) for a range of W
i.
(a)
W
i
=
0.1
(b)
W
i
=
0.5
(c)
W
i
=
0.8
(d)
W
i
=
1
(e)
W
i
=
0.5
(f)
W
i
=
1The profile of the first normal stress differenceN
1 just before the exit (at x
=
-
0.04, and the exit is at x
=
0) of the die for a large range of W
i is shown in Figure 7. These lines are sampled along y-axis at fixed x. Apparently, the maximum value of each line increases as W
i becomes larger, which is consistent with (15). All of these four lines begin with zero at the middle of the die (y
=
0). As for W
i
=
1.0, it goes up gradually until it approaches the wall, where it grows instantly to the peak. Subsequently, it drops dramatically. The profiles of other W
i values smaller than 1.0 are something like shifting the curve of W
i
=
1.0 downwards to the right. In detail, there is no peak in the case of W
i
=
0.1. Generally speaking, a larger W
i results in a higher peak, which indicates the first normal stress dominates the extrusion process.Figure 7
First normal stress difference atx
=
-
0.04 for various W
i with β
=
0.1.Decreasing the viscosity ratioβ from 1 to 0 means gradually increasing the contribution of elastic stress and decreasing the effects of purely viscous stress. When β reaches zero, the model degrades to the upper-convected Maxwell (UCM) model. In order to investigate the role of β in die-swell simulations, simulations under β
=
0.5 and β
=
0.9 are carried out. Figure 8 illustrates the trend that the swelling ratio would increase as W
i number increases, regardless of the value of β. Moreover, the result departs from Tanners theory obviously. Considering an extreme case, as β approaches one, the swelling ratio almost remains a constant even for a large value of W
i according to Tanners theory. However, this does not match the experimental result. For the case β
=
0.9, W
i
=
1 and the case β
=
0.1, W
i
=
1, the recoverable shear R
s are 3 and 2.7, respectively. According to (15), the swelling ratio of the former one should be larger than the later one, while the fact is contrary. Russo and Phillips [18] attribute the discrepancy to the fact that Tanners formula derived from a KBK-Z integral constitutive equation only considered a single relaxation time; however, the Oldroyd-B model owns a retardation time in addition to the relaxation time. This result might be related to the fact that the elongational response for β
=
0.1 is much stronger than the one for β
=
0.9 as illustrated in Figure 9. Thus, the extensional response also plays an important role in swelling process and should be taken into account when predicting the final swelling ratio. The simulation results in [21] demonstrate that flows with a sufficiently high elongational viscosity are able to swell to an apparently large state even at a relatively low flow rate, ignoring their performance in simple shear flows. Nevertheless, decreasing the viscosity ratio would lead to an increase in the recoverable shear and, in turn, increases the value of swelling ration according to (15). By comparing the value in Figures 4 and 8, it is reasonable to say that Tanners theory still makes some sense. Apart from that, with a careful comparison, it is found that when the viscosity ratio β is very large (for the case of β
=
0.9), which represents a very dilute fluid, the swelling ratio suffers a significant drop. In comparison, when β is relatively small (for the case of β
=
0.5), the swelling ratio is much less sensitive. This phenomena was also observed in [22]. Note that because of different definitions of the viscosity ratio, β in this paper is equivalent to 1
-
β in their works.Figure 8
Swelling ratio with the increase ofW
i for β
=
0.5 and β
=
0.9.Figure 9
Planar extensional viscosity for Oldroyd-B model with variousβ.As can be seen from (13), the value of the shear stress would not change once the shear rate is determined. Consequently, the first normal stress difference is the only variable which determines the swelling ratio. Figures 7 and 10 demonstrate that a high proportion of polymeric stress contribution would result in a larger value of N
1, which would provide greater power to swell. According to the analytical solution of the Oldroyd-B model, smaller β would result in a higher value of the first normal stress. Therefore, the first normal stress plays a vital role in determining the behaviour of extrusion flows and could be used to predict the final results of swelling ratio.Figure 10
First normal stress difference atx
=
-
0.04 for various W
i with β
=
0.5 and β
=
0.9.
(a)
β
=
0.5
(b)
β
=
0.9


## 4. Conclusion

In this work, a new method to track the free surface for viscoelastic fluids was proposed. The obstacle of setting a proper free surface boundary condition was well resolved by deriving a Dirichlet pressure boundary condition and introducing a balance force. The proposed numerical approach was implemented opon OpenFOAM, which is widely used in CFD simulations [23, 24]. Comparing with previous literature, the results showed that our solutions were close to those with other methods (e.g., FEM, Spectral Element Method, and MAC), which validated the correctness and effectiveness of the proposed numerical method. Moreover, the results showed that the swelling ratio would decline gradually with the increase of the viscosity ratio. In the meantime, brief analysis of the influence of the viscosity ratio was performed. Moreover, providing appropriate constitutive equations, investigation of die-swell behaviour for new polymeric materials can be performed based on our new platform. In addition, since the FVM systems are suitable for parallelization, the free surface simulations in our platform would be easily scaled to massive processors for parallel computing.



---

*Source: 159831-2015-12-31.xml*