A New Method to Simulate Free Surface Flows for Viscoelastic Fluid

Free surface flows arise in a variety of engineering applications. To predict the dynamic characteristics of such problems, specific numericalmethods are required to accurately capture the shape of free surface.This paper proposed a newmethodwhich 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 proposedmethod for free surface simulation in both Newtonian fluid and viscoelastic fluid.


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 2 Advances in Materials Science and Engineering 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 timedependent 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 Section 2. The numerical results and discussion are reported in Section 3. Conclusions are drawn in Section 4.

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

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 where  is the fluid density, u is the fluid velocity, w is the velocity of the cell,  is the pressure, and  is the polymeric contribution for stress tensor.  and   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 The upper-convected time derivative of the polymer stress tensor is defined as For the convenience of comparison with others' works, variables above are usually converted into their dimensionless form.
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: 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 as  0 = 0.So the boundary condition is given by  ⋅ n = (−I + 2D + ) ⋅ n = (−I + T) ⋅ n = 0, (7) where  is the total stress tensor, namely, the Cauchy stress tensor across the free surface, and T = 2D +  is the extra stress tensor.When solving equations in OpenFOAM, PISO algorithm requires the term −∇ and the term   Δu exists separately in the equation.However, the term  is a combination of , , and D. Therefore, this boundary condition is hard to satisfy directly.It has been known that −∇ + ∇ ⋅  +   Δ = ∇ ⋅ , and in FVM systems, every term in the equation should be integrated over the cell volume .Therefore, the right side of ( 1) is actually solved as ∫  (−∇+∇⋅ +  Δ) = ∫  (∇⋅).According to Gauss's theorem, the spatial derivative term could be converted to integrals over the cell surface  bounding the volume.Then, we have where S is the face area vector and n is the unit normal vector of face area. and  represent faces around a cell and free surface faces, respectively.The dynamic boundary condition requires   ⋅ n  = 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 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: where Γ  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 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  in isolation.While the traction boundary condition is the real physical quantity, it must be satisfied on the fluid.Provided that  and u are solved jointly, no pressure boundary condition is needed anymore [10].

DEVSS Method for Numerical Stability.
In order to enhance the stability of the linear momentum equation, especially for high Wi 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 where   =   +  and ∇ ⋅ Σ  is the added artificial viscosity ratio.The term ∇ ⋅ Σ  = ∇ ⋅  − Δu is treated as a source term.A typical choice of  is   .This strategy greatly increases numerical stability and can be applied in modelling various constitutive models.
Advances in Materials Science and Engineering 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 step   , the velocity field u * , the pressure field  * , the polymer stress field  * , and the fluid mass fluxes ṁ are obtained from the previous time step.The initial net mass fluxes through the interface faces   are calculated according to the previous interface position by subtracting the fluid mass fluxes by the volumetric face fluxes V  .(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.(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.

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 for  and  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 × 10 −6 for the velocity and 1.0 × 10 −7 for pressure and stress, respectively.modelling the Newtonian and Oldroyd-B fluid is shown in Figure 1.The half width of the die  is chosen as the characteristic length.The lengths of the upstream and the downstream are set as  1 and  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  = 1 and

Basic Geometric Shape and Boundary Conditions. A sketch of the geometry and the boundary conditions for
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 along  direction at the inlet is u in = 1.5(1 −  2 ), while the velocity along  direction is zero.Therefore, the average velocity of the inflow is u = 1, and the shear rate near the wall is γ  = 3.The polymer stress tensor at the inlet for the Oldroyd-B model is set according to the analytical solution.Noslip 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  = 0 and   = 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.

Validation of Computational Model with Newtonian Fluid.
By setting Re = 0.01, Wi = 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.

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 where  1 =   −   is the first normal stress difference and   is the shear stress.Both of them are evaluated at the wall.Introducing the parameter of recoverable shear   = ( 1 /2  )  , the swelling ratio can be estimated as For the Oldroyd-B model, the recoverable shear of a fully developed Poiseuille flow is equal to Since in our case γ  = 3 and  = 0.1, the recoverable shear is   = 2.7Wi.
Mesh convergence was tested at Wi = 0.1 with seven meshes, which are listed in Table 1.Δ min / and Δ min / are the relative minimum length along -axis and -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 Δ min / plays a dominant role when determining the welling ratio.Furthermore, if Δ min / ≤ 0.01 and Δ min / ≤ 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 Δ = 10 −6 for M7 and Δ = 10 −5 for other cases (except for M1).
The value of pressure and polymer stress near the singularity is illustrated in Figure 2, where  is in the range of [−0.4,0.4], along with different mesh densities (M4, M6) at Wi = 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   and lowest point for   , where both points are around the exit.Furthermore, the curve of   is sharper than the one of   .As  increases, the value of   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.
In general, with the increase of Wi, the profile of components of the polymer stress should be similar but with a greater value.However, the figure for   is much different from the intuition.Figure 3 shows that the nadir for Wi = 0.1 and Wi = 0.2 which is negative becomes the peak for Wi = 0.5 which is positive.The force along  direction caused by the polymer stress components   is     , where   is vertical component of the standard interfacial normal.Thus, the more likely   is approaching to one (corresponding to a small swelling ratio with a gentle curve), the greater the influence of   would be.As the force increases, the extrusion is growing larger.
The swelling ratio is given in Figure 4 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 Re = 0.001 and Re = 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 γ  and the Weissenberg number Wi are changed.In his case, the recoverable shear rate is far different from the one in our simulation even with the same Wi.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 Wi number, the elastic effect plays a dominant role in determining the die-swell.In addition, slight drop of swelling ratio at low Wi number illustrated in Figure 4 is due to the negative value of   and was also observed in [20].
The steady state profiles of the free surface at various Wi numbers are shown in Figure 5.It can be seen that a larger Wi 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 Wi 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 Wi number, the hexahedron cell is deformed     greatly which increases the mesh nonorthogonality.The maximum mesh nonorthogonality has reached as high as 40 at Wi = 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 Wi = 2.5 using Mesh M1.Contour plots of velocity components   for Wi = 0.1, 0.5, 0.8, 1.0 and   for Wi = 0.5, 1.0 are displayed in Figure 6.By increasing Wi number, the horizontal components of the velocity   increase along  directions.The contour line of   = 1.4 expands upward and forward.The vertical components of the velocity   also increase due to larger elastic effects which results in a larger swelling ratio.
The profile of the first normal stress difference  1 just before the exit (at  = −0.04,and the exit is at  = 0) of the die for a large range of Wi is shown in Figure 7.These lines are sampled along -axis at fixed .Apparently, the maximum value of each line increases as Wi becomes larger, which is consistent with (15).All of these four lines begin with zero at the middle of the die ( = 0).As for Wi = 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 Wi values smaller than 1.0 are something like shifting the curve of Wi = 1.0 downwards to the right.In detail, there is no peak in the case of Wi = 0.1.Generally speaking, a larger Wi results in a higher peak, which indicates the first normal stress dominates the extrusion process.
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 dieswell simulations, simulations under  = 0.5 and  = 0.9 are carried out.Figure 8 illustrates the trend that the swelling ratio would increase as Wi 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 Wi according to Tanners theory.However, this does not match the experimental result.For the case  = 0.9, Wi = 1 and the case  = 0.1, Wi = 1, the recoverable shear   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.
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  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.

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.

( 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 field u * * and the pressure filed  * * with SIMPLE algorithm; (d) moving the interface mesh points similar to step

Figure 4 :
Figure 4: Swelling ratio with the increase of Wi.

Figure 6 :
Figure 6: Contour plots of velocity components   (a)-(d) and   (e)-(f) for a range of Wi.

Figure 9 :
Figure 9: Planar extensional viscosity for Oldroyd-B model with various .
New dimensionless variable Reynolds number Re and Weissenberg number Wi are introduced by scaling the following variables as  is the characteristic velocity (the average velocity in this work), and  is the sum of the   and   ; namely, total viscosity  =   +   . is the viscosity ratio parameter.

Table 1 :
Mesh parameters for convergence validation with Wi = 0.1.