EuropeanCongressonComputationalMethodsinAppliedSciencesandEngineering ECCOMAS 2004 P.Neittaanm¨aki,T.Rossi,S.Korotov,E.On˜ate,J.P´eriaux,andD.Kn¨orzer(eds.) Jyva¨skyl¨a,24–28July2004 STABILIZED FINITE ELEMENT METHODS OF GLS TYPE FOR MAXWELL-B AND OLDROYD-B VISCOELASTIC FLUIDS Marek Behr(cid:1), Dhruv Arora†, and Matteo Pasquali‡ (cid:1)Lehrstuhl fu¨r Numerische Mechanik Technische Universit¨at Mu¨nchen, 85747 Garching, Germany e-mail: [email protected], web page: http://www.lnm.mw.tum.de † Department of Mechanical Engineering and Materials Science Rice University, MS 321, 6100 Main St, Houston, TX 77005, USA e-mail: [email protected] ‡ Department of Chemical Engineering Rice University, MS 362, 6100 Main St, Houston, TX 77005, USA e-mail: [email protected] Key words: Stabilized FEM, Viscoelastic Fluids, Computational Methods. Abstract. We evaluate the stabilized three-field stress-velocity-pressure Galerkin/Least- Squares finite element formulation for viscoelastic fluids, using a benchmark problem of Oldroyd-B flow past a cylinder at various Weissenberg numbers. To address the issue of weak consistency exhibited by low-order velocity interpolations in the context of stabilized formulations, we also employ velocity gradient recovery, and study how such an approach affects the quality ofcomputedresults. We showthat characteristic flowquantities obtained with the new formulation are in good agreement standard DAVSS and DEVSS results, while the cost of fully-implicit velocity gradient computation may be in some cases avoided. 1 Marek Behr, Dhruv Arora and Matteo Pasquali 1 INTRODUCTION Numerical modeling of viscoelastic fluids of rate type is fraught with difficult issues. The viscoelastic stress, or the conformation tensor, must be represented as an additional flow unknown. The advective nature of the constitutive equations, and the interaction of multiple discrete unknown fields (viscoelastic stress, velocity and pressure) both present obstacles to numerical method development. Yet, as has been the case with Newtonian fluid modeling before, these obstacles are being gradually overcome. In the finite element arena, on which we concentrate our attention, adding the Streamline-Upwind/Petrov- Galerkin (SUPG) terms [1,2] to the traditional Galerkin formulation was instrumental in overcoming the difficulties associated with the advectiveterms present in the constitutive equation. The Discontinuous-Galerkin (DG) [3]approach provided similarbenefits, albeit with a marked increase in the complexity of the implementation. Dealing with instabilities related to the interplay of discrete interpolations of the vis- coelastic stress and the velocity fields has been even more challenging. Compatibility conditions on these interpolations have been since formulated, and appear to be well un- derstood [4–6]. These requirements are analogous to the LBB condition [7] for mixed formulations of incompressible flow. In order to satisfy such conditions, complex combi- nations of interpolation functions were needed, such as the 4 ×4 stress sub-element [2] complementing quadratic velocity interpolation. In time, alternative methods were developed which admitted simpler, more efficient, and easier to implement,equal-order interpolations of the viscoelastic stress and velocity. Recentreviews by Baaijens [8] and Keunings [9,10] outline the developmentof the Elastic ViscousStress Splitting(EVSS[11]) approaches, inwhichthe originalGalerkin,SUPG,or DG formulation is modified and new variables or terms are introduced. The state-of-the- art methods are the so-called Discrete EVSS [12,13] methods (DEVSS-G/SUPG [14,15] and DEVSS-G/DG [16]), which include the following features: • the velocity gradient is an independent variable approximated with continuous in- terpolation functions; • the viscous stress is split into two contributions, associated with the continuous velocity gradient and with the discontinuous gradient of the velocity field; • the constitutive equation of the viscoelastic stress (or the conformation) is trans- formedintoweakformwiththestreamline-upwind(SUPG)ordiscontinuousGalerkin (DG) method. Baaijens[8]pointsout thatevenintherelativelysimpletwo-dimensionalflowinanabrupt contraction, mesh-converged results cannot be obtained downstream of corner (location of the singularity) with the DEVSS-G/SUPG and DEVSS-G/DG methods. The common characteristicsof all the EVSS-based methods is the addition to the weak formofthemomentumequationofanelliptictermwhichrepresentsthedifferencebetween 2 Marek Behr, Dhruv Arora and Matteo Pasquali the continuous and discontinuous velocity gradient, and vanishes with increasing mesh resolution. A similar stabilizing term arises naturally when the Petrov-Galerkin approach is replaced with a Least-Squares (LS) method, as was done by Carey’s group [17,18]. Yet, the LS approach carries the usual penalty of high condition numbers of the resulting equation system matrices, and stringent continuity requirements for the interpolation functions (although the orders of interpolation are not higher than those required by the EVSS, DEVSS, and DEVSS-G methods). Moreover, the divergence of the stress is not integrated by parts in the LS form of the momentum equation; thus, the traction at the boundary does not appear naturally in the weak form and imposing boundary conditions at free surfaces and deformable boundaries is difficult. One approach that has been quite successful in circumventing the LBB condition in the case of the Navier-Stokes equations has been the Galerkin/Least-Squares (GLS) method [19]. Here, the stabilizing least-squares form of the governing equations is added to the Galerkin weak form, in the element interiors. The resulting formulation recovers SUPG terms,and also includespressure-stabilizingtermsthat alleviatethe needto satisfy the LBB condition. Moreover,the traction term is present in the GLS formulation, which allows imposing free-surface boundary conditions naturally. This particular approach has been used by Behr and Franca [20] to design and analyze a Galerkin/Least-Squares variational formulation of Navier-Stokes equations involving viscous stress, velocity and pressure as the primary variables, without the usual restrictions on the interpolation function spaces. The GLS approach was extended by Behr [21] to the Maxwell-B and Oldroyd-B constitutive models—hereafter referred to as three-field GLS formulation, or GLS3. The method showed a number of desirable characteristics: • SUPG stabilization terms in the constitutive equation were obtained immediately from the GLS terms; • equal-order interpolations for all flow field variables (viscoelasticstress, velocityand pressure) were admissible; • hence, the implementation was straightforward and the computational cost was modest. The formulation has been then used to compute a benchmark problem of flow in a con- traction. The present work represents an effort to resume and complete the evaluation of the GLS3 approach, and to improve upon it by developing several variants of the method. We will begin by introducing the equations of motion for the Oldroyd-B fluid and its Maxwell-B limit case in Section 2. In Section 3, we present the three-field stress-velocity- pressure Galerkin/Least-Squares finite element formulation. We comment on the weak consistency of this formulation when used with low-order interpolation functions, and introduce two variants that improve the consistency of the method in such situations. In Section 4, we present an example of flow past a cylinder placed in a slit, and compare obtained drag coefficient at various Weissenberg numbers with published results. 3 Marek Behr, Dhruv Arora and Matteo Pasquali 2 PROBLEM STATEMENT We consider a viscous, incompressiblefluid occupying at an instant t ∈ (0,T) a bound- ed region Ωt ⊂ Rnsd, with boundary Γt, where nsd is the number of space dimensions. The velocity and pressure, u(x,t) and p(x,t), are governed by the momentum and mass balance equations: (cid:1) (cid:2) ∂u ρ +u·∇u−f −∇·σ = 0 on Ω ∀t ∈ (0,T), (1) t ∂t ∇·u = 0 on Ω ∀t ∈ (0,T), (2) t where ρ is the fluid density, assumed to be constant, and f(x,t) is an external, e.g., gravitational, force field. The closure is obtained with a constitutive equation relating the stress tensor σ to velocity and pressure fields. Both the Dirichlet and Neumann-type boundary conditions are taken into account, represented as: u = g on (Γ ) , (3) t g n·σ = h on (Γ ) , (4) t h where (Γ ) and (Γ ) are complementary subsets of the boundary Γ . The vector sub- t g t h t scripts signify that this decomposition of Γ may be different for each of the velocity t components. The initial condition consists of a divergence-free velocity field specified over the entire domain: u(x,0) = u , ∇·u = 0 on Ω . (5) 0 0 0 Viscoelastic fluids exhibit dependence of the stress not only on the instantaneous rate of strain, but also on thestrain history. For theupper-convectedMaxwellfluid, also called the Maxwell-B fluid, the constitutive equation acquires an evolutionary character: σ = −pI+T, (cid:1) T+λT = 2µε(u), (6) (cid:1) where the T denotes an upper-convected derivative: (cid:3) (cid:4) (cid:1) ∂T T = +u·∇T− ∇uT+T(∇u)T , (7) ∂t and the rate-of-strain tensor is defined as: (cid:5) (cid:6) 1 ε(u) = ∇u+(∇u)T . (8) 2 Note that, in the case of steady flows considered in the remainder of this article, the time derivative in (7) is dropped, and domain Ω is replaced simply by a constant region Ω. t 4 Marek Behr, Dhruv Arora and Matteo Pasquali The upper-convected Oldroydmodel, also known as the Oldroyd-B model, includesthe Newtonian and Maxwellmodels, and coversthe cases in which an elasticfluid obeying the Maxwell relation is mixedwith a fluid governed by a Newtonian law. This corresponds to a situation in which an elastic polymer with viscosity µ is dissolved in a viscous solvent 1 with viscosity µ : 2 σ = −pI+T, T = T +T , 1 2 (cid:1) T +λT = 2µ ε(u), 1 1 1 T = 2µ ε(u), 2 2 µ +µ = µ. (9) 1 2 The Maxwell fluid is extremely difficult to handle numerically, partially because of the convective character of the stress evolution equation (6). The difficulty is lessened considerably by even a small addition of the Newtonian solvent in an Oldroyd-B fluid. Even so, one may expect the problems normally associated with an advective systems to arise when the relevant constitutive equation is discretized. 3 THREE-FIELD GLS FORMULATION Due to the non-trivial nature of the constitutiveequation, the components of the extra stress haveto be treatedas additional degreesof freedom,complementingthe velocityand pressure fields. A mixed formulation is naturally extended to accommodate the added equation and unknown. However, an arbitrary choice of the stress interpolation often leads to failure, as documented in [2]; possible remedies have been discussed in Section 1. Also, the values of the Weissenberg number for which a Galerkin formulation would remain convergent are extremely small, as seen, e.g., in [22]. The remedy follows the path of the various upwinding methods developed for the advection-diffusion equation, or its advective limit. In [2] both the consistent SUPG and an inconsistent Streamline Upwind(SU) methods are considered. The SU approach is found to stabilizethe Galerkin method sufficiently, but the various coupling effects render the potentially more accurate SUPG method ineffective. Similar conclusion is reached in [3], where a consistent up- winding inherent in the Lesaint-Raviart method has to be augmented by an inconsistent SU addition. We recall here the stress-velocity-pressure formulation introduced in [20] for Stokes equations, and generalizedto Maxwell-Band Oldroyd-B flows in [21]. The n = n (n + tc sd sd 1)/2 independent tensor components of the extra stress T are necessarily treated as 1 additional unknowns, and equation (9) enters the variational formulation directly. We drop the subscript from T , as this is the only stress component entering the formulation 1 explicitly. The interpolation function spaces for the velocity, pressure and extra stress 5 Marek Behr, Dhruv Arora and Matteo Pasquali tensor are given as: (cid:7) (cid:8) (cid:9) (cid:10) Sh = uh | uh ∈ H1h(Ω) nsd ,uh =. gh on Γ , (10) u (cid:7) (cid:8) (cid:9) g(cid:10) Vh = uh | uh ∈ H1h(Ω) nsd ,uh =. 0 on Γ , (11) u (cid:7) (cid:10) g Sh = Vh = ph | ph ∈ H1h(Ω) , (12) p p (cid:7) (cid:8) (cid:9) (cid:10) Sh = Vh = Th | Th ∈ H1h(Ω) ntc . (13) T T In the λ > 0 case, the spaces Sh and Vh must also account for the essential boundary T T conditions for the extra stress at the inflow boundary of the domain. The three-field Galerkin/Least-Squares, or GLS3, velocity-pressure-stress formulation for Oldroyd-B fluid is written as follows: find uh ∈ Sh, ph ∈ Sh and Th ∈ Sh such that: u p T (cid:11) (cid:11) (cid:11) (cid:5) (cid:6) wh ·ρ uh ·∇uh −f dΩ− ∇·whphdΩ+ ε(wh) : ThdΩ Ω (cid:11) (cid:11) Ω (cid:11)Ω + 2µ ε(wh) : ε(uh)− wh ·hhdΓdΩ+ qh∇·uhdΩ 2 (cid:11)Ω Γ(cid:11)h (cid:11) Ω 1 λ (cid:1) + Sh : ThdΩ+ Sh : ThdΩ− Sh : ε(uh)dΩ 2µ 2µ 1 Ω 1 Ω Ω (cid:11) (cid:12)nel 1 (cid:8) (cid:5) (cid:6) (cid:9) + τ ρ uh·∇wh +∇qh −∇·Sh −2µ ∇·ε(wh) MOMρ 2 e=1 Ωe (cid:8) (cid:5) (cid:6) (cid:9) · ρ uh·∇uh −f +∇ph −∇·Th −2µ ∇·ε(uh) dΩ 2 (cid:11) (cid:13) (cid:14) (cid:12)n el 1 λ (cid:1) + τ 2µ Sh + Sh −ε(wh) CONS 1 2µ 2µ Ωe 1 1 e=1 (cid:13) (cid:14) 1 λ (cid:1) : Th + Th −ε(uh) dΩ 2µ 2µ 1 1 (cid:11) (cid:12)n el + τ ∇·whρ∇·uhdΩ = 0, ∀wh∈ Vh, ∀qh∈ Vh, ∀Sh∈ Vh. (14) CONT u p T e=1 Ωe The stabilization parameters τ and τ follow standard definitions given, e.g., MOM CONT in [23]. The parameter τ is taken here as: CONS (cid:15) (cid:3) (cid:4) max 1, he , λ|uh| > 1 τ = 2λ|uh|2 2 (15) CONS max(1,he), λ|uh| ≤ 1 2 Remark 1 The addition of the least-squares form of the momentum equation, i.e., the τ -term MOM in(14), stabilizesthemethod against well-knownadverseeffects of under-diffusivityof the 6 Marek Behr, Dhruv Arora and Matteo Pasquali Galerkin discretization of the momentum equation at high element-level Peclet numbers, and a possible lack of compatibility between velocity and pressure function spaces. Remark 2 The least-squares form of the continuity equation, i.e., the τ -term in (14), improves CONT the convergence of non-linear solvers at high Reynolds numbers. Remark 3 The least-squares form of the constitutive equation, i.e. the τ -term in (14), stabilizes CONS the method against two further causes of spurious numerical oscillations: the under- diffusivity of the Galerkin discretization of the constitutive equation at high Weissenberg numbers, and a possible lack of compatibilitybetween velocityand stress function spaces. As shown in [20], this term allows for arbitrary combinations of interpolation functions for velocity and the extra stress, which do not have to satisfy a specific LBB condition. The combination of the stabilization terms circumventing the inf-sup conditions pro- videsabsolutefreedominselectingtheinterpolationfunctionspaces, allowinginparticular convenient equal-order interpolation for velocity, pressure and the stress. In the example that follows, a piecewise bi-linear interpolation is in fact used for all three fields (so called Q1Q1Q1 element). Although simpletoimplementand computationallyefficient,linearorbi-linearvelocity interpolation suffers from one important drawback. In the presence of the Newtonian solvent (Oldroyd-B fluid), the second derivatives of velocity in the stabilization terms are not adequately represented—they are either zero in the linear case, or small and unrelated to actual second derivative of velocity field in the bi-linear case. The GLS3 formulation, normally considered to be consistent, is only weakly consistent in such cases as pointed out in [24]. Solutions to this problem in general involveobtaining a continuous velocity gradient, with well-defined derivatives inside element domain. To this end, one canemployeitherof twostrategies: asimplifiedrecoveryof nodal velocitygradientusinga least-squares approach, or introduction of a continuous velocitygradient field, resulting in a four-field velocity-pressure-stress-gradient stabilized formulation. The former approach is taken here; the latter approach leads to a GLS4 formulation which will be explored in a future article. The GLS4 bears similarities to DEVSS approach, although the extra velocity gradient field is used only for the proper representation of second derivatives of velocity in the case of Oldroyd-B fluid, i.e., to improve accuracy, whereas in DEVSS, such field is necessary for bypassing the compatibility condition on velocity and stress interpolations, i.e., to prevent a breakdown of discretization. The least-squares recoveryproceduresolves, in a decoupledmanner, the weak problem: given uh—piecewise linear or bi-linear computed velocity field—find Lh ∈ Sh such that: L (cid:11) (cid:11) Kh : LhdΩ = Kh : ε(uh) ∀Kh∈ Vh, (16) L Ω Ω 7 Marek Behr, Dhruv Arora and Matteo Pasquali where the function spaces: (cid:16) (cid:17) (cid:8) (cid:9) Sh = Vh = Lh | Lh ∈ H1h(Ω) nsd×nsd , (17) L L are by definition composed of continuous functions, with non-trivial derivatives in the element interior. The recovered values Lh are then used to form the residual of (14) at the next iteration of the non-linear iteration procedure, namely, in the computation of the τ -terms involving second derivatives of velocity. The mass matrix of the equation MOM system originating from (16) can be eitherlumped, for increased computational efficiency, or solvedwitha director iterativesolverfor increasedaccuracy. Depending on this choice, we arrive at two GLS3 variants: GLS3-L Base GLS3 method with recovery of continuous velocity gradient using lumped mass matrix, GLS3-M Base GLS3 method withrecoveryof continuousvelocitygradientusing consistent mass matrix. 4 NUMERICAL EXAMPLE Flow of Oldroyd-B fluid past a circular cylinder placed between parallel fixed plates is a standard benchmark problem for two-dimensional viscoelastic flow simulation. The results for the ratio of cylinder diameter to slit width of 1/8 have been reported by Sun et al. [16] using DAVSS-G/DG formulation, and by Pasquali [25] using DEVSS-G/SUPG formulation, and are reported here for the GLS3 formulation variants. The flow domain is shown in Figure 1. The center of the cylinder is assumed to be located at (0,0), the radius of the cylinder is taken as R = 1, and the slit half-width as h = 8. The distance from the cylinder center to the inflow boundary Γ and the outflow i boundary Γ is 10 and 20, respectively. A parabolic flow profile is prescribed on the inflow o and outflow boundaries: u = 1.5 (Q/h) (1−x2/h2), 1 2 u = 0, 2 (cid:5) (cid:6) T = 2 λ µ −1.5x /h2 2, 11 (cid:5) 1 2(cid:6) T = µ −1.5x /h2 , 12 1 2 T = 0. (18) 22 with flow rate per unit thickness Q = 8 reflecting the value used in [16,25]. A no-slip condition was prescribed at the cylinder wall Γ and the slit wall Γ. Finally, a symmetry c s condition was applied at the upstream and downstream symmetry lines Γ and Γ . u d The fluid and solvent viscosities are µ = 0.41 and µ = 0.59, respectively, giving the 1 2 ratio of solvent viscosity to total viscosity of 0.59. The relaxation time λ is varied to arrive at the desired Weissenberg number: Ws = Qλ/hR. (19) 8 Marek Behr, Dhruv Arora and Matteo Pasquali Γ s Γ h Γ i o R Γ Γ Γ u c d Figure 1: Flow past a circular cylinder: domain description. A semi-structured mesh with 25,025 nodes and 24,619 bi-linear elements covers the domain, as shown in Figure 2. The size of the smallest elements, i.e., the ones adjacent to the cylinder,is approximately 0.010×0.018 in the radial and circumferentialdirection, respectively. Figure 2: Flow past a circular cylinder: finite element mesh. A characteristic quantity for this flow field is the drag force exerted by the fluid on the cylinder, and its variation with the Weissenberg number. The drag force is computed using the integral: (cid:11) F = −2 e ·σh ·ndΓ, (20) d 1 Γc where e is the horizontal unit vector and n is the unit normal at the cylinder surface 1 pointing out of the flow domain. Note that this line integral includes contributions from the pressure ph, viscoelastic stress Th, and the viscous stress 2µ ε(uh). The drag for 2 Weissenberg numbers 0.0 to 2.0 is tabulated in Table 1 and shown in Figure 3. The agreement between GLS3-M and results of Sun et al. [16] is excellent up to Weissenberg number of 1, while GLS3-L under-predicts the drag by less than 1%. For Weissenberg numbers higher than 1, the GLS3-M and GLS3-L formulations depart further from Sun et al. [16] results, especially for Ws > 1.8, where the difference reaches 2.5%. Note that Ws (cid:5) 2.0 is the limit of convergence for many methods when applied to this problem, including that of Sun et al. [16]. 9 Marek Behr, Dhruv Arora and Matteo Pasquali Ws Sun et al. GLS3-L GLS3-M 0.0000 15.722 15.6720 15.7004 0.1000 15.705 15.6748 15.7062 0.2000 15.691 15.6648 15.6960 0.3000 15.682 15.6516 15.6822 0.4000 15.668 15.6384 15.6684 0.5000 15.662 15.6292 15.6584 0.6000 15.652 15.6270 15.6556 0.7000 15.665 15.6346 15.6624 0.8000 15.693 15.6536 15.6808 0.9000 15.728 15.6852 15.7122 1.0000 15.774 15.7306 15.7570 1.1000 15.843 15.7898 15.8160 1.2000 15.924 15.8626 15.8886 1.3000 16.007 15.9502 15.9764 1.4000 16.120 16.0492 16.0754 1.5000 16.235 16.1628 16.1892 1.6000 16.360 16.2854 16.3120 1.7000 16.530 16.4158 16.4428 1.8000 16.689 16.5472 16.5754 1.9000 16.884 16.6630 16.6910 2.0000 17.148 16.7608 16.7912 Table 1: Flow past a circular cylinder: drag as a function of Weissenberg number. 10