Viscoelastic Flow Calculations Research Articles

A parallel adaptive viscoelastic flow solver with template based dynamic mesh refinement

A parallel adaptive mesh refinement algorithm has been incorporated into the side-centered finite volume method [Sahin, A stable unstructured finite volume method for parallel large-scale viscoelastic fluid flow calculations. J. Non-Newtonian Fluid Mech., 166 (2011) 779–791] in order to obtain highly accurate numerical results for viscoelastic fluid flow problems. The present recursive mesh refinement algorithm is based on a conformal refinement of unstructured quadrilateral/hexahedral elements with templates based on 1:3 refinement of edges. In order to transfer cell-centered data between source and target meshes, a second-order conservative interpolation (remapping) technique similar to the work of Menon and Schmidt [Supermesh construction for conservative interpolation on unstructured meshes: An extension to cell–centered finite-volume variables. Comput. Methods Appl. Mech. Eng., 200 (2011), 2797–2804] are employed and the approach has been extended for side-centered data. The proposed framework has been applied to the classical benchmark problem of an Oldroyd-B fluid past a confined circular cylinder in a rectangular channel and a sphere falling in a circular tube. The calculations confirm that high accuracy can be achieved with the present adaptive mesh refinement.

Parallel large-scale numerical simulations of purely-elastic instabilities behind a confined circular cylinder in a rectangular channel

The parallel large-scale unstructured finite volume method proposed in [Sahin, A stable unstructured finite volume method for parallel large-scale viscoelastic fluid flow calculations, J. Non-Newtonian Fluid Mech. 166 (2011) 779–791] has been applied to investigate the three-dimensional creeping flow of an Oldroyd-B fluid past a confined circular cylinder in a rectangular channel at relatively high Weissenberg numbers. The numerical method is based on side-centered finite volume method where the velocity vector components are defined at the mid-point of each cell face, while the pressure term and the extra stress tensor are defined at element centroids. The present arrangement of the primitive variables leads to a stable numerical scheme and it does not require any ad-hoc modifications in order to enhance the pressure–velocity–stress coupling. The combination of the present numerical method with the log-conformation representation proposed in [R. Fattal, R. Kupferman, Constitutive laws for the matrix-logarithm of the conformation tensor, J. Non-Newtonian Fluid Mech. 123 (2004) 281–285] and the geometric non-nested multilevel preconditioner for the Stokes system have enabled us to simulate large-scale viscoelastic fluid flow problems on highly parallel machines. The calculations are presented for an Oldroyd-B fluid past a confined circular cylinder in a rectangular channel at relatively high Weissenberg numbers and compared to those obtained for Newtonian fluids. The present numerical calculations reveal three-dimensional purely-elastic instabilities in the wake of a confined single cylinder which is in accord with the earlier experimental results in the literature. In addition, the flow field is found out to be no longer symmetric in the wake of the cylinder at high Weissenberg numbers.

A stable unstructured finite volume method for parallel large-scale viscoelastic fluid flow calculations

A new stable unstructured finite volume method is presented for parallel large-scale simulation of viscoelastic fluid flows. The numerical method is based on the side-centered finite volume method where the velocity vector components are defined at the mid-point of each cell face, while the pressure term and the extra stress tensor are defined at element centroids. The present arrangement of the primitive variables leads to a stable numerical scheme and it does not require any ad-hoc modifications in order to enhance the pressure–velocity–stress coupling. The log-conformation representation proposed in [R. Fattal, R. Kupferman, Constitutive laws for the matrix–logarithm of the conformation tensor, J. Non-Newtonian Fluid Mech. 123 (2004) 281–285] has been implemented in order improve the limiting Weissenberg numbers in the proposed finite volume method. The time stepping algorithm used decouples the calculation of the polymeric stress by solution of a hyperbolic constitutive equation from the evolution of the velocity and pressure fields by solution of a generalized Stokes problem. The resulting algebraic linear systems are solved using the FGMRES( m) Krylov iterative method with the restricted additive Schwarz preconditioner for the extra stress tensor and the geometric non-nested multilevel preconditioner for the Stokes system. The implementation of the preconditioned iterative solvers is based on the PETSc library for improving the efficiency of the parallel code. The present numerical algorithm is validated for the Kovasznay flow, the flow of an Oldroyd-B fluid past a confined circular cylinder in a channel and the three-dimensional flow of an Oldroyd-B fluid around a rigid sphere falling in a cylindrical tube. Parallel large-scale calculations are presented up to 523,094 quadrilateral elements in two-dimension and 1,190,376 hexahedral elements in three-dimension.

Comparison of Different Formulations for the Numerical Calculation of Unsteady Incompressible Viscoelastic Fluid Flow

In this paper we compare different methods currently used in the sta- bilization of numerical simulations of time-dependent viscoelastic fluid flows de- scribed with the Oldroyd-B and related models. The methods under consideration, based on the separation of newtonian-like components from the stress tensor, are applied to a finite volume analysis of two simple benchmark problems (the plane Poiseuille startup and pulsated flows), for which analytical solutions are known. The relative performances of each method are evaluated regarding stability, accu- racy and efficiency.

From atomistic trajectories to primitive paths to tube models: linking atomistic simulations with the reptation theory of polymer dynamics

Three different approaches have been reported in the last few years for identifying topological constraints and generating ensembles of primitive paths (PPs) in entangled, multi-chain polymeric systems. In addition to providing predictions for the static (statistical) properties of the underlying entanglement network, such a development has opened up the way to interfacing atomistic simulation data with reptation, admittedly the most successful phenomenological theory of polymer dynamics for entangled systems. The link between atomistic molecular dynamics simulation results and reptation theory is achieved by geometrically constructing the effective tube around each primitive chain and then documenting chain motion in terms of a curvilinear diffusion inside the effective tube around the coarse-grained chain contour. The outcome of such a topological and dynamical mapping is the computation of the function ψ(s,t), the probability that a segment s of the primitive chain remains inside the initial tube after time t. The purpose of this article is to discuss the new computational developments and to describe how they can be used to advance and improve our current understanding of polymer melt viscoelasticity. We emphasize in particular the opportunity opened today to bring together three different approaches to polymer dynamics (in addition to acquiring reliable experimental data): atomistic simulations, mesoscopic entanglement networks, and tube models. By developing methodologies that can consistently map the results of accurate computer models and simulations of polymer structure and dynamics onto theoretical treatments based on phenomenological concepts (that sometimes defy precise definition), our hope is that we can get a deeper understanding of the predominant relaxation mechanisms in entangled polymers, both in the linear and nonlinear regime, and succeed in our effort to encode them in the form of a constitutive equation for large-scale viscoelastic flow calculations.

A new, convenient way of imposing open-flow boundary conditions in two- and three-dimensional viscoelastic flows

Open flow boundaries are often present in viscoelastic flow calculations; their presence is not dictated by the physics of the problem, but rather by the need of truncating the computational domain. Viscoelastic liquids flowing in complex two- and three-dimensional domains are normally modeled by hyperbolic transport equations for the viscoelastic stress or conformation tensor, v ⋅ ∇ S = F ( ∇ v , S ) − G ( S ) , where S is the stress or conformation tensor, v is the velocity, and ∇ denotes gradient in space. In steady flows, the streamlines are the characteristics of these hyperbolic equations and boundary conditions on S are necessary where the liquid enters the flow domain. Open flow boundaries are almost always located in regions of fully-developed, rectilinear flow. Traditionally, several methods have been used to prescribe inflow conditions; each of them has one or more drawbacks in terms of applicability to general models, computational expense and complexity, and inability to deal with unknown flowrates, inflow-outflow boundaries, or unstructured meshes. Here, we propose a new, general way of imposing inflow boundary conditions based on solving the coupled algebraic equations of fully developed flow at the inflow, i.e., solving the equation F ( ∇ v , S ) − G ( S ) = 0 at the inflow coupled with the flow inside the domain. The equation holds because v ⋅ ∇ S ≡ 0 in fully developed rectilinear flow. Imposing the inflow boundary condition in this fashion is fully general and does not require additional programming in solvers based on finite elements, spectral methods, and finite differences, whether or not Newton’s method is used for solving the nonlinear algebraic equations arising from the discretized partial differential equation. We test this method and find excellent agreement with analytical results in combined Poiseuille and Couette flow of Oldroyd-B and Giesekus liquids in 2-D and 3-D channels and annuli, for which analytical expressions of the velocity and conformation (elastic stress) fields are available. We demonstrate that the new method yields shorter or equal upstream lengths than traditional ones.

A fully explicit characteristic based split (CBS) scheme for viscoelastic flow calculations

AbstractIn this paper, an explicit characteristic based split (CBS) scheme is proposed for the numerical solution of incompressible viscoelastic flow equations. The scheme proposed is free from simultaneous solution to the matrices arising from the finite element discretization of the governing equations. The experience gained from the solution of Newtonian fluid dynamics problems has been applied to the solution of viscoelastic flows. The Oldroyd‐B model has been employed to solve two benchmark problems of viscoelastic flow. They are viscoelastic flow past a circular cylinder and viscoelastic flow through planar contraction geometry. The results show that the solutions obtained are stable for the Weissenberg or Deborah number range studied in this paper. The solutions obtained at lower Weissenberg or Deborah numbers are accurate and agree excellently with the majority of available numerical data. However at higher Weissenberg or Deborah numbers, results show some sign of negative influence of the artificial dissipation added to the discrete constitutive equations. Copyright © 2004 John Wiley Sons, Ltd.

A first error indicator for micro–macro simulations of viscoelastic fluids having closed-form constitutive equations

An a posteriori error indicator for viscoelastic flow calculations is proposed by the author for a recently introduced method. This new method is derived from the Brownian configuration field approach and can be applied to models for which closed-form expressions already exist. Numerical results are presented in the context of the spectral element method to show the reliability of the error indicator. The benchmark problem of the flow of an Oldroyd-B fluid past a cylinder in a channel is used and we show that this new method is always more accurate than using a constitutive equation.

Upwinding with deferred correction (UPDC): an effective implementation of higher-order convection schemes for implicit finite volume methods

A conceptually simple and computationally inexpensive implementation of higher-order convection interpolation schemes for stable and faster convergence of the implicit finite volume method (FVM) is evaluated for multi-dimensional viscoelastic flow calculations as well as convection-dominated Newtonian flow calculations. Based on the deferred-correction method, this effective formulation consistently satisfies the four rules that guarantee bounded numerical solutions. In addition, artificial numerical diffusion, commonly found in numerical solutions with the first-order upwinding convection scheme, can be effectively controlled by consistently using a third-order quadratic upstream interpolation over the entire computational domain. Extensive tests performed for the steady-state wall-driven square enclosure flow of a Newtonian fluid at Reynolds number ( Re) up to 5000 demonstrate that the formulation is stable, accurate, and converges well. The accuracy of the formulation for multi-dimensional viscoelastic flow calculations is evaluated by solving the Oldroyd-B equations for the problem of a fully two-dimensional steady frictionless plane flow, and comparing the numerical results with the exact solution available. With the first-order upwinding scheme, the results show that, for hyperbolic constitutive models, regardless of the grid size, there always exists artificial diffusion whenever the flow streamlines are not closely aligned with the grid lines. With the new implementation, artificial numerical diffusion is effectively controlled and the results demonstrate third-order accuracy.

How accurate is your solution?: Error indicators for viscoelastic flow calculations

The present work is an extension of earlier results of Owens [R.G. Owens, Comput. Methods Appl. Mech. Eng. 164 (1998) 375–395] and is an attempt to provide some theoretical undergirding to the question of appropriate error indicators for numerical solutions of flows of viscoelastic fluids having a differential constitutive equation. In particular, it is shown that the local elemental residual for the elastic stresses only accounts for the so-called stress cell error and as a consequence is, in general, an inadequate measure of the local error. An improved error indicator is then proposed which takes account of the transmitted error. Acknowledging that the stress errors form the major contribution to our error indicator for sufficiently large Deborah numbers we have devised a method of computing the error indicator on an element-by-element basis. Numerical results are presented to show how the approximate error and the ‘exact’ error obtained by calculating the difference between the numerical solution and a reference calculation decrease with mesh refinement. As an illustration of the use of our error indicator we proceed to describe an adaptive spectral element method for flow of an Oldroyd-B fluid past a sphere in a tube. The use of consistent upwinding is shown to result in greater accuracy and stability than is possible with a Galerkin approach. The results are compared with those in the literature.

On the discrete EVSS method

This paper is a completion of a previous paper written by two of the authors [M. Fortin, R Guénette, R. Pierre, Numerical analysis of the modified EVSS method, Comput. Methods. Appl. Mech. Engrg. 143 (1997) 79–95] concerning the DEVSS method, which is a discrete variant of the popular elastic-viscous-split-stress (EVSS) method introduced by [Rajagopalan, R.A. Brown, R.C. Armstrong, Finite element methods for calculation of steady viscoelastic flow using constitutive equations with a Newtonian viscosity, J. Non-Newtonian Fluid Mech. 36 (1990) 159–199] within the context of viscoelastic fluids. The method uses the rate of deformation tensor as an additional variable. In this paper, it is shown that a simple modification of the theoretical framework introduced in the first paper [M. Fortin, R Guénette, R. Pierre, Numerical analysis of the modified EVSS method, Comput. Methods. Appl. Mech. Engrg. 143 (1997) 79–95] allows the treatment of multi-mode constitutive equations of viscoelastic flows as well as that of the velocity gradient introduced as the additional variable instead of the rate of deformation tensor. Moreover, this general framework based on a generalized Brezzi–Babuska theory, allows different discretizations for the stress and the additional variable. The paper presents a general criteria insuring the stability of the method. Three particular mixed finite elements satisfying this criteria are studied. An error analysis is performed on Stokes and viscoelastic flows.

A posteriori error estimates for spectral element solutions to viscoelastic flow problems

In this paper we develop an a posteriori error estimate for viscoelastic flow calculations using spectral elements. The error estimate is compared with the ‘exact’ error obtained by computing the difference between the numerical solution and a reference calculation for the benchmark problem of viscoelastic flow past a rigid sphere in a tube.

High weissenberg number simulation of an annular extrudate swell using the differential type constitutive equation

AbstractAnnular extrudate swell simulations at high Weissenberg numbers were made using a differential type constitutive equation. The streamline‐upwinding method with a sub‐element for extra stress components, which is called SU4 × 4, is one of the best mixed finite element methods for computation of viscoelastic flows. Planar and capillary extrudate swell calculations at high Weissenberg numbers (We > 1000) were accomplished by SU4 × 4. However, annular extrudate swell simulations at high We by SU4 × 4 were not successful. The calculated We was less than about 4. A new calculation technique using a Newton‐Raphson discretization of the equation of motion was developed. This technique is called a “new under‐relaxation method.” The calculated We of annular extrudate swell simulation by the new under‐relaxation method with SU4 × 4 was about 6∼250 times larger than those by SU4 × 4. Reasonable calculation results were obtained in an annular flow and a capillary extrudate swell by this method, and the reliability and the utility of the new under‐relaxation method are shown. It is now possible to consider the swell shapes of annular extrudate under industrially useful conditions. The calculated swelling ratios were also compared with experimental ones.

Brownian configuration fields and variance reduced CONNFFESSIT

Stochastic simulation techniques, such as Brownian dynamics, provide us an extremely powerful tool for solving the usually nonlinear equations describing polymer dynamics in solutions and melts [1]. However, the most challenging problems (e.g. the investigation of the universal behaviour of long polymer chains, or the flow calculation based on stochastic simulation techniques) involve a very large number of degrees of freedom and hence require an enomous amount of computer time. In order to solve such problems on currently available computers it is therefore necessary to develop strategies to drastically suppress the level of the fluctuations in the simulations. The purpose of this note is to show that the recently proposed concept of Brownian configuration fields [2] in viscoelastic flow calculations can be regarded as an extremely powerful extension of variance reduction techniques based on parallel process simulation.

2‐D time‐dependent viscoelastic flow calculations using connffessit

AbstractTwo‐dimensional time‐dependent calculations for a molecular model of finite extensibility in the journal‐bearing geometry are presented. The flow is considered to be incompressible and isothermal. The momentum conservation equation is integrated using a time‐marching procedure in which local ensembles of dumbbells act as stress calculators. The calculations are based on the Calculation of non‐Newtonian flows: finite elements and stochastic simulation technique (CONNFFESSIT) and combine deterministic (finite elements) and stochastic techniques to advance the velocity and stress fields in time. The ability of CONNFFESSIT to treat models for which no closed‐form constitutive equation can be derived is illustrated by performing calculations using FENE dumbbells. Significant differences in the stress field between the true FENE and the linearized FENE‐P are found, especially during the inception period. Steady‐state kinematics are, however, identical within error bars for both FENE and FENE‐P and for the Newtonian fluid. The essential algorithm of 2‐D CONNFFESSIT is detailed, as well as experience gathered from its parallel and vector versions.

A comparative study of the discontinuous Galerkin and continuous SUPG finite element methods for computation of viscoelastic flows

The performance of the discontinuous Galerkin (DG) method and the continuous streamline upwind Petrov-Galerkin (SUPG) method in the computation of steady viscoelastic flows have been compared in this paper. For the continuous SUPG method it is necessary to include the elastic-viscous-split-stress (EVSS) formulation to satisfy the compatibility requirement on the velocity-stress interpolation (SUPG/ EVSS), while for the discontinuous Galerkin method both the pure mixed formulation (DG/MIX) and the EVSS formulation (DG/EVSS) have been tested. We first examine a benchmark smooth problem, i.e. the steady flow of an upper convected Maxwell (UCM) fluid around a sphere falling along the axis of a cylindrical tube, by using high-order finite element approximation. DG/MIX is found not to be a stable discretization in this problem, its limiting Weissenberg number decreases with increasing the interpolation order; on the other hand, DG/EVSS is an accurate and stable algorithm, its solution converges as the approximation order increases (the so-called p-convergence) and it produces almost identical stress boundary layers to that of the SUPG/EVSS algorithm at a Weissenberg number as high as 2.0. To explore benefits of the discontinuous interpolation for the stress in the DG method, we further examine two benchmark problems with stress singularities by using low-order interpolations on sufficiently fine meshes; they are the axisymmetric steady flows of the modified Chilcott and Rallison (MCR) fluid through a four-to-one abrupt contraction and in a tube with the stick-to-slip transition in the wall velocity boundary condition. In the contraction problem the DG/EVSS method is found to be a significant improvement over the SUPG/EVSS method, quality of the solutions is better for DG/EVSS and a much higher Weissenberg number limit is observed. In the stick-slip problem both DG/EVSS and SUPG/EVSS are very robust, they confront no limit of the Weissenberg number and predict asymptotic stress profiles as the number becomes sufficiently large. In both problems the pure mixed method DG/MIX produces severe stress oscillations near the singularities and fails to converge at medium Weissenberg numbers.

Effect of artificial stress diffusivity on the stability of numerical calculations and the flow dynamics of time-dependent viscoelastic flows

In this work, we investigate the effect of the introduction of a stress diffusive term into the classical Oldroyd-B constitutive equation on the numerical stability of time-dependent viscoelastic flow calculations. The channel Poiseuille flow at Re ⪢ 1 and O(1) We is chosen as a test problem. Through a linear stability analysis, we demonstrate that the introduction of a small amount of (dimensionless) diffusivity, typically of the order 10 −3, does not affect the critical eigenmodes of the viscoelastic Orr-Sommerfeld problem appreciably. However, a diffusive term of that magnitude is shown to have a significant influence on the singular eigenmodes of the classical Oldroyd-B model, associated with the continuum spectra. A finite amplitude perturbation is constructed as a linear superposition of the eigenvectors corresponding to the most unstable eigenvalues of the problem. This is superimposed on the steady Poiseuille flow solution to provide the initial conditions for time-dependent simulations. The numerical algorithm involves a fully spectral spatial discretization and a semi-implicit second order integration in time. For the Oldroyd-B fluid, depending on the magnitude of the initial perturbation, numerical instabilities set in at relatively short times while the components of the conformation tensor increase monotonically in magnitude with time. Introduction of a diffusive term into this model is shown to stabilize the calculations remarkably, and for a three-dimensional simulation with Re = 5000 and We = 1, no instabilities were observed even at very large times. The effect of the magnitude of the diffusivity on the stability and the flow dynamics is addressed through a direct comparison of the results with those obtained for the Oldroyd-B model.

Letter to the Editor: Mesh refinement limits in viscoelastic flow calculations?

First Page

A new mixed finite element method for viscoelastic flows governed by differential constitutive equations

A new mixed finite element method is presented that has improved numerical stability and has similar numerical accuracy compared to the EVSS/FEM developed by Rajagopalan et al. (J. Non-Newtonian Fluid Mech. 36 (1990) 159–192). The new method, denoted EVSS-G/FEM, is based on separate bilinear interpolation of all relevant components of the velocity gradient tensor. Lagrangian bilinear approximations of this variable are formed by least-squares interpolation of the derivatives of a Lagrangian biquadratic approximation to the velocity field that is generated by solution of the momentum and continuity equations using a standard mixed formulation for velocity and pressure. Elastic components of the stress are computed by solving the constitutive equation using either the streamline upwinding (SU) or streamline upwind Petrov-Galerkin (SUPG) methods with bilinear interpolation. The motivation for the interpolation of the velocity gradients and for the stress interpolation is to force the compatibility between these variables that is needed to satisfy differential viscoelastic constitutive equations in the limit where the velocity field vanishes, as it does near solid boundaries. The enhanced numerical stability of the EVSS-G/FEM over the EVSS/FEM formulation is demonstrated for the calculation of the linear stability of planar Couette flow of a UCM fluid, where closed form results exist for the linear stability of the steady-state flow; hence, any instability detected in the numerical solution is a result of either the finite element approximation or the time integration method. The second possibility is eliminated by using a fully implicit time integration method. Calculations with the EVSS-G/FEM are stable for De>100, whereas calculations with the EVSS/FEM become numerically unstable for De>5. The EVSS-G formulation is combined with finite element discretizations that incorporate local adaptive mesh refinement for increasing the accuracy of the calculations in regions where the stress and velocity gradients change rapidly. The accuracy of the EVSS-G/FEM is demonstrated by steady-state calculations for the flow between eccentric rotating cylinders, flow through a wavy-walled tube, and for flow through a square array of cylinders; each has been used before as a test problem for the accuracy of viscoelastic flow calculations. As in previous calculations, discretization of the constitutive equation by SUPG gives superior accuracy at high values of De compared to the application of SU.

A finite volume approach for calculation of viscoelastic flow through an abrupt axisymmetric contraction

This paper reports the first convergent numerical algorithm for the steady inertialess flow of an Upper-Convected Maxwell (UCM) fluid through a four-to-one abrupt axisymmetric contraction. The Finite Volume Method (FVM) is adopted along with a stream function-vorticity approach in the Elastic Viscous Split Stress (EVSS) form with a first-order upwind approximation applied to the convective terms of the stress constitutive equation. A staggered volume discretization of the flow variables eliminates the stress singularity from the computational domain without loosing any flow physics. The volume integrals of the governing equations over the flow domain result in a system of nonlinear algebraic equations that are solved iteratively by a semi-implicit line-to-line method with a pseudo-trasient term added to the stress constitutive equation. Computations of an UCM fluid using this method are carried out to a much higher value of the Deborah number ( De) than previous numerical simulations using the Finite Element Method (FEM). The solutions are found to be smooth, stable, and convergent with the finger stress tensor remaining posiitive-definite throughout the flow domain. Calculations are not performed above De = 6.25 because of the decreasing pseudo-time step constraint at higher elasticity. The finite volume algorithm approximates better solutions upon mesh refinement, demonstrates the smoothness and the mathematical well-posedness of this problem, and predicts a Newtonian-like flow structure near the singularity for an UCM fluid. Using this new method, the inertialess flow of an UCM fluid through 4:1 abrupt axisymmetric contraction, for the first time, produces larger corner vortices at higher values of De.