Advances in Materials Science and Engineering

Volume 2015, Article ID 159831, 10 pages

http://dx.doi.org/10.1155/2015/159831

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

^{1}College of Computer, National University of Defense Technology, Changsha, Hunan 410073, China^{2}State Key Laboratory of High Performance Computing, National University of Defense Technology, Changsha, Hunan 410073, China^{3}School of Chemical Engineering and Analytical Science, Manchester Institute of Biotechnology, University of Manchester, Manchester M13 9PL, UK

Received 11 September 2015; Revised 1 December 2015; Accepted 13 December 2015

Academic Editor: Olanrewaju Ojo

Copyright © 2015 Yu Cao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.

#### 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 Section 2. 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 bywhere is the fluid density, is the fluid velocity, 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 is the rate of deformation tensor, which can be expressed asThe 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. New dimensionless variable Reynolds number Re and Weissenberg number are introduced by scaling the following variables asin which is the characteristic length (the half width of the die in this work), 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.

##### 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 byThe 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 . So the boundary condition is given bywhere is the total stress tensor, namely, the Cauchy stress tensor across the free surface, and is the extra stress tensor. When solving equations in OpenFOAM, PISO algorithm requires the term and the term exists separately in the equation. However, the term is a combination of , **,** and . 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 havewhere is the face area vector and is the unit normal vector of face area. and represent faces around a cell and free surface faces, respectively.

The dynamic boundary condition requires ; 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 asIt 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 byIt 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 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 high number flow, the discrete elastic split stress (DEVSS) numerical strategy proposed by Matallah et al. [11] is used. The momentum equation can be rewritten aswhere and is the added artificial viscosity ratio. The term 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.

###### 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 , 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 .(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 ;(b)setting up the pressure boundary condition and estimating the balance force ;(c)calculating the latest velocity field and the pressure filed 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 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 for the velocity and 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 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 and , 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 and .