A finite volume approach for unsteady viscoelastic fluid flows

A finite volume, time‐marching for solving time‐dependent viscoelastic flow in two space dimensions for Oldroyd‐B and Phan Thien–Tanner fluids, is presented. A non‐uniform staggered grid system is used. The conservation and constitutive equations are solved using the finite volume method with an upwind scheme for the viscoelastic stresses and an hybrid scheme for the velocities. To calculate the pressure field, the semi‐implicit method for the pressure linked equation revised method is used. The discretized equations are solved sequentially, using the tridiagonal matrix algorithm solver with under‐relaxation. In both, the full approximation storage multigrid algorithm is used to speed up the convergence rate. Simulations of viscoelastic flows in four‐to‐one abrupt plane contraction are carried out. We will study the behaviour at the entrance corner of the four‐to‐one planar abrupt contraction. Using this solver, we show convergence up to a Weissenberg number We of 20 for the Oldroyd‐B model. No limiting Weissenberg number is observed even though a Phan Thien–Tanner model is used. Several numerical results are presented. Smooth and stable solutions are obtained for high Weissenberg number. Copyright © 2002 John Wiley & Sons, Ltd.     

A conservative semi‐implicit method for coupled surface–subsurface flows in regional scale

A semi‐implicit method for coupled surface–subsurface flows in regional scale is proposed and analyzed. The flow domain is assumed to have a small vertical scale as compared with the horizontal extents. Thus, after hydrostatic approximation, the simplified governing equations are derived from the Reynolds averaged Navier–Stokes equations for the surface flow and from the Darcy's law for the subsurface flow. A conservative free‐surface equation is derived from a vertical integral of the incompressibility condition and extends to the whole water column including both, the surface and the subsurface, wet domains. Numerically, the horizontal domain is covered by an unstructured orthogonal grid that may include subgrid specifications. Along the vertical direction a simple z‐layer discretization is adopted. Semi‐implicit finite difference equations for velocities and a finite volume approximation for the free‐surface equation are derived in such a fashion that, after simple manipulation, the resulting discrete free‐surface equation yields a single, well‐posed, mildly nonlinear system. This system is efficiently solved by a nested Newton‐type iterative method that yields simultaneously the pressure and a non‐negative fluid volume throughout the computational grid. The time‐step size is not restricted by stability conditions dictated by friction or surface wave speed. The resulting algorithm is simple, extremely efficient, and very accurate. Exact mass conservation is assured also in presence of wetting and drying dynamics, in pressurized flow conditions, and during free‐surface transition through the interface. A few examples illustrate the model applicability and demonstrate the effectiveness of the proposed algorithm. Copyright © 2015 John Wiley & Sons, Ltd.     

A robust multi‐grid pressure‐based algorithm for multi‐fluid flow at all speeds

This paper reports on the implementation and testing, within a full non‐linear multi‐grid environment, of a new pressure‐based algorithm for the prediction of multi‐fluid flow at all speeds. The algorithm is part of the mass conservation‐based algorithms (MCBA) group in which the pressure correction equation is derived from overall mass conservation. The performance of the new method is assessed by solving a series of two‐dimensional two‐fluid flow test problems varying from turbulent low Mach number to supersonic flows, and from very low to high fluid density ratios. Solutions are generated for several grid sizes using the single grid (SG), the prolongation grid (PG), and the full non‐linear multi‐grid (FMG) methods. The main outcomes of this study are: (i) a clear demonstration of the ability of the FMG method to tackle the added non‐linearity of multi‐fluid flows, which is manifested through the performance jump observed when using the non‐linear multi‐grid approach as compared to the SG and PG methods; (ii) the extension of the FMG method to predict turbulent multi‐fluid flows at all speeds. The convergence history plots and CPU‐times presented indicate that the FMG method is far more efficient than the PG method and accelerates the convergence rate over the SG method, for the problems solved and the grids used, by a factor reaching a value as high as 15. Copyright © 2003 John Wiley & Sons, Ltd.     

Coupling of finite element method and discontinuous Galerkin method to simulate viscoelastic flows

In this paper, a numerical method, which is about the coupling of continuous and discontinuous Galerkin method based on the splitting scheme, is presented for the calculation of viscoelastic flows of the Oldroyd‐B fluid. The momentum equation is discretized in time by using the Adams‐Bashforth second‐order algorithm, and then decoupled via the splitting approach. Considering the Oldroyd‐B constitutive equation, the second‐order Runge‐Kutta approach is selected to complete the temporal discretization. As for the spatial discretizations, the fundamental purpose is to make the best of finite element method (FEM) and discontinuous Galerkin (DG) method to handle different types of equations. Specifically speaking, for the subequations, FEM is chosen to treat the Poisson and Helmholtz equations, and DG is employed to deal with the nonlinear convective term. In addition, because of the hyperbolic nature, DG is also utilized to discretize the Oldroyd‐B constitutive equation spatially. This coupled method avoids resorting to extra stabilization technique occurred in standard FEM framework even for moderately high values of Weissenberg number and also reduces the complexity compared with unified DG scheme. The Oldroyd‐B model is applied to investigate several typical and challenging benchmarks, such as the 4:1 planar contraction flow and the lid‐driven cavity flow, with a wide range of Weissenberg number to illustrate the feasibility, robustness, and validity of our coupled method.     

Solving viscoelastic free surface flows of a second‐order fluid using a marker‐and‐cell approach

This work is concerned with the numerical simulation of two‐dimensional viscoelastic free surface flows of a second‐order fluid. The governing equations are solved by a finite difference technique based on the marker‐and‐cell philosophy. A staggered grid is employed and marker particles are used to represent the fluid free surface. Full details for the approximation of the free surface stress conditions are given. The resultant code is validated and convergence is demonstrated. Numerical simulations of the extrudate swell and flow through a planar 4:1 contraction for various values of the Deborah number are presented. Copyright © 2006 John Wiley & Sons, Ltd.     

A finite element thermally coupled flow formulation for phase‐change problems

A finite element, thermally coupled incompressible flow formulation considering phase‐change effects is presented. This formulation accounts for natural convection, temperature‐dependent material properties and isothermal and non‐isothermal phase‐change models. In this context, the full Navier–Stokes equations are solved using a generalized streamline operator (GSO) technique. The highly non‐linear phase‐change effects are treated with a temperature‐based algorithm, which provides stability and convergence of the numerical solution. The Boussinesq approximation is used in order to consider the temperature‐dependent density variation. Furthermore, the numerical solution of the coupled problem is approached with a staggered incremental‐iterative solution scheme, such that the convergence criteria are written in terms of the residual vectors. Finally, this formulation is used for the solutions of solidification and melting problems validating some numerical results with other existing solutions obtained with different methodologies. Copyright © 2000 John Wiley & Sons, Ltd.     

Die-swell,splashing drop and a numerical technique for solving the Oldroyd B model for axisymmetric free surface flows

This work deals with the development of a numerical method for simulating viscoelastic axisymmetric free surface flow of an Oldroyd B fluid. A novel formulation is developed for the computation of the non-Newtonian extra-stress components on rigid boundaries and on the symmetry axis. The full free surface stress conditions are employed. The resulting governing equations are solved by finite differences on a Marker-and-cell (MAC) type grid. Validation is provided by simulating a pipe flow problem. The classical die-swell problem is solved and swelling ratios are provided. The height of the splash caused by a falling liquid drop for various Reynolds and Weissenberg numbers is then studied, and the height of the splash is shown to diminish with increasing viscoelasticity.     

A coupled surface‐subsurface model for hydrostatic flows under saturated and variably saturated conditions

In this paper, the governing differential equations for hydrostatic surface‐subsurface flows are derived from the Richards and from the Navier‐Stokes equations. A vertically integrated continuity equation is formulated to account for both surface and subsurface flows under saturated and variable saturated conditions. Numerically, the horizontal domain is covered by an unstructured orthogonal grid that may include subgrid specifications. Along the vertical direction, a simple z‐layer discretization is adopted. Semi‐implicit finite difference equations for velocities, and a finite volume approximation for the vertically integrated continuity equation, are derived in such a fashion that, after simple manipulation, the resulting discrete pressure equation can be assembled into a single, two‐dimensional, mildly nonlinear system. This system is solved by a nested Newton‐type method, which yields simultaneously the (hydrostatic) pressure and a nonnegative fluid volume throughout the computational grid. The resulting algorithm is relatively simple, extremely efficient, and very accurate. Stability, convergence, and exact mass conservation are assured throughout also in presence of wetting and drying, in variable saturated conditions, and during flow transition through the soil interface. A few examples illustrate the model applicability and demonstrate the effectiveness of the proposed algorithm.     

Application of differential quadrature method to solve entry flow of viscoelastic second‐order fluid

The entry flow of viscoelastic second‐order fluid between two parallel plates is discussed. The governing equations of vorticity and the streamfunction are expanded with respect to a small parameter that characterizes the elasticity of the fluid by means of the standard perturbation method. By using the differential quadrature method with only a few grid points, high‐accurate numerical solutions are obtained. The numerical results show a lot of the features of a viscoelastic second‐order fluid. Copyright © 1999 John Wiley & Sons, Ltd.     

Spectral element method for viscoelastic flows in a planar contraction channel

A new algorithm, which combines the spectral element method with elastic viscous splitting stress (EVSS) method, has been developed for viscoelastic fluid flows in a planar contraction channel. The system of spectral element approximations to the velocity, pressure, extra stress and the rate of deformation variables is solved by a preconditioned conjugate gradient method based on the Uzawa iteration procedure. The numerical approach is implemented on a planar four‐to‐one contraction channel for a fluid governed by an Oldroyd‐B constitutive equation. The behaviour of the Oldroyd‐B fluids in the contraction channel is investigated with various Weissenberg numbers. It is shown that numerical solutions obtained here agree well with experimental measurements and other numerical predictions. Copyright © 2003 John Wiley & Sons, Ltd.     

A robust algorithm for computing fluid flows on highly non‐smooth staggered grids

A new method for computing the fluid flow in complex geometries using highly non‐smooth and non‐orthogonal staggered grid is presented. In a context of the SIMPLE algorithm, pressure and physical tangential velocity components are used as dependent variables in momentum equations. To reduce the sensitivity of the curvature terms in response to coordinate line orientation change, these terms are exclusively computed using Cartesian velocity components in momentum equations. The method is then used to solve some fairly complicated 2‐D and 3‐D flow field using highly non‐smooth grids. The accuracy of results on rough grids (with sharp grid line orientation change and non‐uniformity) was found to be high and the agreement with previous experimental and numerical results was quite good. Copyright © 2008 John Wiley & Sons, Ltd.     

Parallel domain decomposition method for finite element approximation of 3D steady state non‐Newtonian fluids

We introduce a stabilized finite element method for the 3D non‐Newtonian Navier–Stokes equations and a parallel domain decomposition method for solving the sparse system of nonlinear equations arising from the discretization. Non‐Newtonian flow problems are, generally speaking, more challenging than Newtonian flows because the nonlinearities are not only in the convection term but also in the viscosity term, which depends on the shear rate. Many good iterative methods and preconditioning techniques that work well for the Newtonian flows do not work well for the non‐Newtonian flows. We employ a Galerkin/least squares finite element method, with stabilization parameters adjusted to count the non‐Newtonian effect, to discretize the equations, and the resulting highly nonlinear system of equations is solved by a Newton–Krylov–Schwarz algorithm. In this study, we apply the proposed method to some inelastic power‐law fluid flows through the eccentric annuli with inner cylinder rotation and investigate the robustness of the method with respect to some physical parameters, including the power‐law index and the Reynolds number ratios. We then report the superlinear speedup achieved by the domain decomposition algorithm on a computer with up to 512 processors. Copyright © 2015 John Wiley & Sons, Ltd.     

Modeling and simulation of the viscoelastic fluid mold filling process by level set method

A model for unifying a viscoelastic fluid and a Newtonian fluid is established, in which the governing equations for the viscoelastic fluid and the Newtonian fluid are successfully united into a system of generalized Navier–Stokes equations. A level set method is set up to solve the model for capturing the moving interface in the mold filling process. The physical governing equations are solved by the finite volume method on a non-staggered grid and the interpolation technique on the collocated grid is used for the pressure-velocity and the stress-velocity decoupling problems. The level set and its reinitialization equation are solved by the finite difference method, in which the spatial derivatives are discretized by the 5th-order Weighted Essentially Non-Oscillatory (WENO) scheme, and the temporal derivatives are discretized by the 3rd-order Total Variation Diminishing Runge–Kutta (TVD-R–K) scheme. The validity of the method is verified by some benchmark problems. Then a simulation of viscoelastic fluid mold filling process is pursued with the method. The moving interface and all the information of the physical quantities during the injection process are captured. The die swelling phenomenon is found in the simulation. The influences of elasticity and viscosity on the physical quantities such as stresses etc. in the mold filling process are analyzed. Numerical results show that elastic characteristics such as the stretch and die swelling etc. reinforce accordingly as Weissenberg number increases. Pressures increase continuously in the mold filling process and the pressure maintains the maximum value at the inlet. Injection velocity is proportional to injection pressure. A higher viscosity leads to a higher pressure distribution, that is, the pressure decreases as Reynolds number increases.     

A fixed‐grid b‐spline finite element technique for fluid–structure interaction

We present a fixed‐grid finite element technique for fluid–structure interaction problems involving incompressible viscous flows and thin structures. The flow equations are discretised with isoparametric b‐spline basis functions defined on a logically Cartesian grid. In addition, the previously proposed subdivision‐stabilisation technique is used to ensure inf–sup stability. The beam equations are discretised with b‐splines and the shell equations with subdivision basis functions, both leading to a rotation‐free formulation. The interface conditions between the fluid and the structure are enforced with the Nitsche technique. The resulting coupled system of equations is solved with a Dirichlet–Robin partitioning scheme, and the fluid equations are solved with a pressure–correction method. Auxiliary techniques employed for improving numerical robustness include the level‐set based implicit representation of the structure interface on the fluid grid, a cut‐cell integration algorithm based on marching tetrahedra and the conservative data transfer between the fluid and structure discretisations. A number of verification and validation examples, primarily motivated by animal locomotion in air or water, demonstrate the robustness and efficiency of our approach. Copyright © 2013 John Wiley & Sons, Ltd.     

Parallel finite element simulation of mooring forces on floating objects

The coupling between the equations governing the free‐surface flows, the six degrees of freedom non‐linear rigid body dynamics, the linear elasticity equations for mesh‐moving and the cables has resulted in a fluid‐structure interaction technology capable of simulating mooring forces on floating objects. The finite element solution strategy is based on a combination approach derived from fixed‐mesh and moving‐mesh techniques. Here, the free‐surface flow simulations are based on the Navier–Stokes equations written for two incompressible fluids where the impact of one fluid on the other one is extremely small. An interface function with two distinct values is used to locate the position of the free‐surface. The stabilized finite element formulations are written and integrated in an arbitrary Lagrangian–Eulerian domain. This allows us to handle the motion of the time dependent geometries. Forces and momentums exerted on the floating object by both water and hawsers are calculated and used to update the position of the floating object in time. In the mesh moving scheme, we assume that the computational domain is made of elastic materials. The linear elasticity equations are solved to obtain the displacements for each computational node. The non‐linear rigid body dynamics equations are coupled with the governing equations of fluid flow and are solved simultaneously to update the position of the floating object. The numerical examples includes a 3D simulation of water waves impacting on a moored floating box and a model boat and simulation of floating object under water constrained with a cable. Copyright © 2003 John Wiley & Sons, Ltd.     

Pure axial flow of viscoelastic fluids in rectangular microchannels under combined effects of electro-osmosis and hydrodynamics

This paper presents an analysis of the combined electro-osmotic and pressure-driven axial flows of viscoelastic fluids in a rectangular microchannel with arbitrary aspect ratios. The rheological behavior of the fluid is described by the complete form of Phan-Thien–Tanner (PTT) model with the Gordon–Schowalter convected derivative which covers the upper convected Maxwell, Johnson–Segalman and FENE-P models. Our numerical simulation is based on the computation of 2D Poisson–Boltzmann, Cauchy momentum and PTT constitutive equations. The solution of these governing nonlinear coupled set of equations is obtained by using the second-order central finite difference method in a non-uniform grid system and is verified against 1D analytical solution of the velocity profile with less than 0.06% relative error. Also, a parametric study is carried out to investigate the effect of channel aspect ratio (width to height), wall zeta potential and the Debye–Hückel parameter on 2D velocity profile, volumetric flow rate and the Poiseuille number in the mixed EO/PD flows of viscoelastic fluids with different Weissenberg numbers. Our results show that, for low channel aspect ratios, the previous 1D analytical models underestimate the velocity profile at the channel half-width centerline in the case of favorable pressure gradients and overestimate it in the case of adverse pressure gradients. The results reveal that the inapplicability of the Debye–Hückel approximation at high zeta potentials is more significant for higher Weissenberg number fluids. Also, it is found that, under the specified values of electrokinetic parameters, there is a threshold for velocity scale ratio in which the Poiseuille number is approximately independent of channel aspect ratio.     

Numerical simulation of free‐surface flow using the level‐set method with global mass correction

A new numerical method that couples the incompressible Navier–Stokes equations with the global mass correction level‐set method for simulating fluid problems with free surfaces and interfaces is presented in this paper. The finite volume method is used to discretize Navier–Stokes equations with the two‐step projection method on a staggered Cartesian grid. The free‐surface flow problem is solved on a fixed grid in which the free surface is captured by the zero level set. Mass conservation is improved significantly by applying a global mass correction scheme, in a novel combination with third‐order essentially non‐oscillatory schemes and a five stage Runge–Kutta method, to accomplish advection and re‐distancing of the level‐set function. The coupled solver is applied to simulate interface change and flow field in four benchmark test cases: (1) shear flow; (2) dam break; (3) travelling and reflection of solitary wave and (4) solitary wave over a submerged object. The computational results are in excellent agreement with theoretical predictions, experimental data and previous numerical simulations using a RANS‐VOF method. The simulations reveal some interesting free‐surface phenomena such as the free‐surface vortices, air entrapment and wave deformation over a submerged object. Copyright © 2009 John Wiley & Sons, Ltd.     

Numerical simulation of 2D lid‐driven cavity flow with CLEARER algorithm on extremely highly skewed grids at high Reynolds numbers

It is crucial to deal with the grid non‐orthogonality effectively in solving the flow in complex geometries, especially at high Reynolds numbers. In this study, the newly proposed Coupled and Linked Equations Algorithm Revised‐ER (CLEARER) algorithm is adopted to solve this problem successfully. In CLEARER algorithm the second relaxation factor is introduced in constructing the contravariant interface velocities, by setting it to a low value. CLEARER algorithm can overcome the severe grid non‐orthogonality and non‐linearity of equations effectively. After the numerical results with CLEARER are validated with the benchmark solutions, this algorithm is used to solve the lid‐driven flow in inclined cavity with inclination angles varying from 10 to 170^°, and Reynolds numbers varying from 5000 to 15 000. The streamlines and the centerline velocity distributions are provided in detail for all cases, which may offer some guidance for the study in this area. Copyright © 2010 John Wiley & Sons, Ltd.     

An extended finite element method for the simulation of particulate viscoelastic flows

We present an extended finite element method (XFEM) for the direct numerical simulation of the flow of viscoelastic fluids with suspended particles. For moving particle problems, we devise a temporary arbitrary Lagrangian–Eulerian (ALE) scheme which defines the mapping of field variables at previous time levels onto the computational mesh at the current time level. In this method, a regular mesh is used for the whole computational domain including both fluid and particles. A temporary ALE mesh is constructed separately and the computational mesh is kept unchanged throughout the whole computations. Particles are moving on a fixed Eulerian mesh without any need of re-meshing. For mesh refinements around the interface, we combine XFEM with the grid deformation method, in which nodal points are redistributed close to the interface while preserving the mesh topology. Our method is verified by comparing with the results of boundary fitted mesh problems combined with the conventional ALE scheme. The proposed method shows similar accuracy compared with boundary fitted mesh problems and superior accuracy compared with the fictitious domain method. If the grid deformation method is combined with XFEM, the required computational time is reduced significantly compared to uniform mesh refinements, while providing mesh convergent solutions. We apply the proposed method to the particle migration in rotating Couette flow of a Giesekus fluid. We investigate the effect of initial particle positions, the Weissenberg number, the mobility parameter of the Giesekus model and the particle size on the particle migration. We also show two-particle interactions in confined shear flow of a viscoelastic fluid. We find three different regimes of particle motions according to initial separations of particles.     

3D free‐surface flow computation using a RANSE/Fourier–Kochin coupling

A coupling method for numerical calculations of steady free‐surface flows around a body is presented. The fluid domain in the neighbourhood of the hull is divided into two overlapping zones. Viscous effects are taken in account near the hull using Reynolds‐averaged Navier–Stokes equations (RANSE), whereas potential flow provides the flow away from the hull. In the internal domain, RANSE are solved by a fully coupled velocity, pressure and free‐surface elevation method. In the external domain, potential‐flow theory with linearized free‐surface condition is used to provide boundary conditions to the RANSE solver. The Fourier–Kochin method based on the Fourier–Kochin formulation, which defines the velocity field in a potential‐flow region in terms of the velocity distribution at a boundary surface, is used for that purpose. Moreover, the free‐surface Green function satisfying this linearized free‐surface condition is used. Calculations have been successfully performed for steady ship‐waves past a serie 60 and then have demonstrated abilities of the present coupling algorithm. Copyright © 2003 John Wiley & Sons, Ltd.