Analysis of complex viscoelastic flows using a finite element method

Abstract

The field of computational fluid mechanics of viscoelastic flows has been well explored in the three decades since its inception. Still, even with the vast amount of work detailed in the literature, much remains to be done towards the improvement of models of viscoelastic fluids and the improvement of the numerical methods used to solve the set of governing equations. The work contained in this document is concentrated in the latter of these areas. The main goal of this body of work is to develop a robust, efficient simulation package to model three-dimensional viscoelastic flows. In order to accomplish this goal, improvements to the numerical methods and equation formulation were necessary to help reduce the overall size of the equation set used to describe viscoelastic flows in three-dimensional geometries. In order to test their viability for use in reducing the overall size of the problem, concepts involving changing the formulation of the equations and the numerical methods used to find the solution to the equations were first implemented and analyzed in a previously developed two-dimensional finite element simulation package. Implementation and analysis is discussed of a formulation change involving decoupling the calculation of the velocity gradient interpolant equation and the momentum and mass continuity equations in the DEVSS-G formulation.

(cont.) Two different decoupled methods for computing the velocity gradient, one using a global least squares approximation and the other a local patch algorithm, are explored. While both methods reduce to the true velocity gradient with mesh refinement, the patch algorithm is shown to require significantly more mesh refinement than the global least squares approximation to order to attain equivalent refinement of the solution. Comparison of the two methods taking into account the additional refinement requirements of the local patch algorithm makes clear the superiority of the decoupled global least squares approximation for calculation of the velocity gradient interpolant. The versatility and robustness of the decoupled form of the DEVSS-G equations are demonstrated through the addition and modification of the evolution equations describing the stress of the polymer as well as new physical quantities of the flow. A time-dependent, free-surface finite element method is developed in which an evolution equation derived from the kinematic boundary condition is used to describe the height of the free surface as a function of time.

(cont.) This new evolution equation is incorporated into the decoupled formulation by simply adding an additional step to the time integration to evaluate the change in the height of the surface during the current timestep and then updating the element locations in the deformable region of the mesh. Application of the new equation in this manner requires no knowledge of the direct dependence of the system on changes in the new quantity, allowing for quick and easy implementation. Incorporation of more advanced constitutive equations is used as further example of the utility of the decoupled form of the DEVSS-G equations. For most continuum based constitutive equations, the dependence of the equations on the flow variables can be expressed explicitly, allowing for the coupled set of equations to be solved with Newton's method. However, the dependence of the stress on the flow cannot be explicitly written for more advanced constitutive equations such as those derived from kinetic theory or those employing Brownian dynamics, greatly hindering the performance of Newton's method in locating the solution to the system. As an illustrative example, incorporation into the decoupled equation formulation of the closed form of the Adaptive-Length-Scale model (ALS-C) is presented.

(cont.) Simulations are presented capturing for the first time the pressure drop enhancement with increasing viscoelasticity of the model of the flow of a Boger fluid in the 4:1:4 axisymmetric contraction-expansion geometry observed experimentally (Rothstein et al., 2001). Simulations of the flow of a 4-mode FENE-P model fluid within the geometry are also presented. Though its dependence on the flow field can be expressed analytically, the cost of computation using multimode models is typically prohibitive when using fully coupled equation sets as the overall problem size grows considerably with the addition of each new mode. Incorporation of the 4-mode model within the decoupled equation formulation adds relatively little computational cost to the overall calculation. Employing the formulation and numerical methods developed herein, a new three-dimensional finite element package is described for simulating confined viscoelastic flows. To make the package more robust, a number of different boundary conditions are included for modeling different geometries used in polymer processing. To help reduce the burden associated with mesh refinement in three-dimensional meshes, a commercial meshing package utilizing o-grid refinement for localization of refinement is employed.

(cont.) Furthermore, to allow for computation of the large equation sets typically associated with three-dimensional geometries, a parallel implementation of the three-dimensional simulation package is developed based on the two-dimensional parallel method developed by Caola et al. ((Caola et al., 2001), (Caola et al., 2002)). Simulation results demonstrating the accuracy and performance of the method are presented. As a test of the robustness of the three-dimensional method, simulations of the flow of Newtonian and Oldroyd-B fluids through a periodic, linear array of cylinders are presented. Comparisons with previous calculations for the Oldroyd-B flow in an infinitely wide domain with no variations in the direction of the width show the same trend in the drag force on the cylinder with increasing viscoelasticity as well as in the size and shape of the vortices formed in the gap between the cylinders. The study of this flow includes effects of modeling the cross section of the flow as an infinite domain with no variation in the direction of the width, an infinite domain of periodic computational width, an infinite domain of periodic computational width and a symmetric flow above and below the cylinders, and a bounded domain with solid walls located 4 cylinder radii apart.