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 .
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 , while the velocity along direction is zero. Therefore, the average velocity of the inflow is , and the shear rate near the wall is . 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 ( and ). 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 setting , , , and , 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 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 .
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 aswhere 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 , 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 and , the recoverable shear is .
Mesh convergence was tested at with seven meshes, which are listed in Table 1. and 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 plays a dominant role when determining the welling ratio. Furthermore, if and , the relative error of the swelling ratio is below 0.2%. The solution convergence is reached by using M6. The time interval is for M7 and 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 , along with different mesh densities (M4, M6) at . 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.
(a)
(b)
(c)
(d)
In general, with the increase of , 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 and which is negative becomes the peak for 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 and 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 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 are changed. In his case, the recoverable shear rate is far different from the one in our simulation even with the same . 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 number, the elastic effect plays a dominant role in determining the die-swell. In addition, slight drop of swelling ratio at low 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 numbers are shown in Figure 5. It can be seen that a larger 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 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 number, the hexahedron cell is deformed greatly which increases the mesh nonorthogonality. The maximum mesh nonorthogonality has reached as high as 40 at 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 using Mesh M1.
Contour plots of velocity components for , 0.5, 0.8, 1.0 and for , 1.0 are displayed in Figure 6. By increasing number, the horizontal components of the velocity increase along directions. The contour line of expands upward and forward. The vertical components of the velocity also increase due to larger elastic effects which results in a larger swelling ratio.
(a)
(b)
(c)
(d)
(e)
(f)
The profile of the first normal stress difference just before the exit (at , and the exit is at ) of the die for a large range of is shown in Figure 7. These lines are sampled along -axis at fixed . Apparently, the maximum value of each line increases as becomes larger, which is consistent with (15). All of these four lines begin with zero at the middle of the die (). As for , it goes up gradually until it approaches the wall, where it grows instantly to the peak. Subsequently, it drops dramatically. The profiles of other values smaller than 1.0 are something like shifting the curve of downwards to the right. In detail, there is no peak in the case of . Generally speaking, a larger 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 die-swell simulations, simulations under and are carried out. Figure 8 illustrates the trend that the swelling ratio would increase as 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 according to Tanners theory. However, this does not match the experimental result. For the case , and the case , , 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 is much stronger than the one for 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 ), which represents a very dilute fluid, the swelling ratio suffers a significant drop. In comparison, when is relatively small (for the case of ), 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 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 , 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.
(a)
(b)
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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
Thanks are due to the computational resources provided by the University of Manchester, partial funding from the National Natural Science Foundation of China under Grants no. 61221491, no. 61303071, and no. 61120106005, and Open fund (no. 201503-01 and no. 201503-02) from State Key Laboratory of High Performance Computing.