ASIA-PACIFIC JOURNAL OF CHEMICAL ENGINEERING Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 Published online 11 September 2009 in Wiley InterScience (www.interscience.wiley.com) DOI:10.1002/apj.384 Research Article Influence of the Weißenberg number on the stability of Oldroyd kind fluids N. Scurtu1 * and E. Bänsch2 1 Department of Aerodynamics and Fluid Mechanics, Brandenburg University of Technology Cottbus, Siemens-Halske-Ring 14, 03046 Cottbus, Germany 2 Department of Applied Mathematics III, Friedrich-Alexander University of Erlangen-Nürnberg, Haberstr. 2, 91058 Erlangen, Germany ABSTRACT: This paper is concerned with nonlinear rheological fluids of Oldroyd type. We present a (formal) stability analysis of the corresponding system of equations, showing stability limits on the Weißenberg number in certain cases. To this end, we proceed in several steps, thus separating the possible sources for instabilities. First, a spectral analysis of the linearized Oldroyd system is presented. Then, the influence of the βa -term on the stability of the constitutive stress equation and of the full Oldroyd system is examined. Moreover, because this stability analysis is of formal and linear nature, we complement it by numerical simulations for the system showing that the upper limit of the Weißenberg number found by the stability analysis is fairly sharp. We thereby try to shed some light on the high Weißenberg number problem, that is, the problem why in certain cases there seem to exist no solutions to the Oldroyd problem for large Weißenberg numbers. 2009 Curtin University of Technology and John Wiley & Sons, Ltd. KEYWORDS: CFD; non-Newtonian fluids; viscoelastic flow simulation; Oldroyd fluid; stability analysis INTRODUCTION The equations governing the flow of viscoelastic nonNewtonian fluids have been intensively studied for instance because of their industrial applications in, for example, polymer processing. Besides so called generalized Newtonian fluids, where the viscosity is a function of the strain, more sophisticated models introduce a material constitutive equation for taking into account memory effects of such polymeric materials. In 1950, Oldroyd proposed a differential equation for modeling viscoelastic fluids.[1] During the last few decades, significant progress has been made in the development of numerical algorithms for stable and accurate solution of these type of viscoelastic flow problems, see Refs [2–7]. For incompressible fluids, the constitutive assumptions in combination with the equations of motion taking into account elastic and memory effects lead to a highly nonlinear system of partial differential equations. For Oldroyd fluids, the flow is described by the following dimensionless system of equations of hyperbolic–parabolic–elliptic type: *Correspondence to: N. Scurtu, Department of Aerodynamics and Fluid Mechanics, Brandenburg University of Technology Cottbus, Siemens-Halske-Ring 14, 03046 Cottbus, Germany. E-mail: scurtu@tu-cottbus.de 2009 Curtin University of Technology and John Wiley & Sons, Ltd. ∂τ + (u · ∇)τ + β (τ, ∇u) We a ∂t +τ − 2αD(u) = 0, ∂u Re + (u · ∇)u ∂t −2(1 − α)divD(u) − div τ + ∇p = f , div u = 0 (1) The unknown fields are the velocity u, the pressure p, and the symmetric stress tensor τ . In order to fulfill the frame indifference principle, an objective time derivative is introduced in the constitutive equation by adding the so called β-term to the material time derivative: βa (T , ∇u) := TW − WT − a(DT + TD) (2) where D = D(u) = (∇u + ∇ T u)/2 and W = W (U ) = (∇u − ∇ T u)/2 are the rate of deformation tensor and the spin tensor, respectively. The Oldroyd system is to be solved in a domain ⊂ Rd , d ≥ 2, supplemented by appropriate initial and boundary conditions. For the velocity field, the boundary condition may be of, for example, Dirichlet type on the whole boundary = ∂, whereas for the stress field Dirichlet conditions have to impose on the inflow boundary part − only. 658 N. SCURTU AND E. BÄNSCH Asia-Pacific Journal of Chemical Engineering The Oldroyd system is characterized by four dimensionless parameters, the Reynolds number Re > 0, the Weißenberg number We > 0, defined as the ratio of a relaxation time of the polymer over a typical time scale in the flow, the rate of viscoelastic viscosities α ∈ [0, 1] and a ∈ [−1, 1]. The slip parameter a was interpreted by Johnson and Segalman[8] as a measure of the nonaffinity of polymer deformation, that is, the fractional stretch of the polymeric material with respect to the stretch of the flow field. In the limit of We = 0, the system reduces to the Navier–Stokes equations, corresponding to the purely viscous fluid, whereas the limit case α = 1 corresponds to the purely elastic Maxwell fluid. The particular cases of the Oldroyd fluid with a lower and upper convective derivatives, respectively, are referred to as Oldroyd-A (a = −1) and Oldroyd-B (a = 1) fluid, see Ref.[1]. The most critical parameter of the system is the the Weißenberg number. Although small values of the Weißenberg (describing slightly elastic liquids) can be handled well, there are unsolved problems for high Weißenberg numbers corresponding to a large relaxation and highly elastic fluids. The problems associated with such cases of highly elastic viscoelastic fluids, where We is large, are consequently called the high Weißenberg number problem. An upper limit for Wecr , beyond which no solution exists, is found in all published works applying numerical techniques to differential and integral models, see Refs [2,3,7,9,10]. Minor changes in the constitutive equation and/or the algorithm employed can lead to higher limit values of We. As We approaches the critical value Wecr , it is often observed that spurious oscillations appear in the field variables; the stress components are then more severely affected than the velocity components, yielding large and erroneous stress gradients. The spurious oscillations have no physical background, and their wavelength depends on the mesh used for the discretization. Mesh refinement and corner strategies affect the critical conditions for breakdown, but it is difficult to discern an overall consistent trend in published works. For a number of benchmark problems in case of steady flows, the reviews[2,10,11] showed that the limits in the maximum attainable Weißenberg number still exist, independently of the numerical method used for solving the Oldroyd problem. Concerning the existence of solutions for the Oldroyd system with suitable boundary conditions, Renardy[12] has summarized the known results that are basically four types of existence results: For the time-dependent case, existence of solutions locally in time and for small data, globally in time has been proved in Banach spaces by FernándezCara et al .[13] For the Oldroyd-B model without the convective term in the constitutive equation, existence of solution in a convex polygonal domain in the stationary case and for small relaxation times has been proved by Picasso and Rappaz.[14] Mathematical existence of a solution in a Banach space was proved by Bonito et al.[15] for the time-dependent system but also without the convective term. For the case of a = 0, Lions and Masmoudi[16] have proved the global existence of a weak solution that satisfies energy identities in the two-dimensional case. The main aim of this paper is to present a (formal) stability analysis of the Oldroyd system and thereby shed some light on the question about the mechanism that leads to an upper limit of We. Upon performing a linearized stability analysis, we will show that the system can become unstable at supercritical Weißenberg numbers. Because this analysis is of formal nature, we complement it by numerical simulations. We computationally investigate different situations, where instabilities show up, and try to classify the origin of the arising instabilities. Our findings imply that no accurate numerical scheme should overshoot the limit of Wecr . Simulations that go much beyond this limit should be estimated with care. Although we mainly focus on the Oldroyd fluid model, our results may be extended to other nonNewtonian fluids like the Giesekus model[17] or the PTT model.[18] The rest of this paper is organized as follows. Section Stability Analysis contains a formal spectral analysis of the linearized Oldroyd system. For the full nonlinear system, stability limits on the Weißenberg number will be presented. In Section Numerical Study, we perform numerical simulations on three test cases for the validation of the stability limits found in Section Stability Analysis. We conclude the paper by a summary in Section Conclusions. (1) results on existence locally in time for initial value problems, when the solution is a small perturbation of the initial data; (2) results on global time existence and asymptotic decay if the initial conditions are small perturbations of the rest state; In this section, we study stability properties of the Oldroyd system by linear perturbation technique. We proceed in several steps. First, we examine the case, when neglecting all nonlinear terms, like the convective term and the β-term. Then we analyze the stress transport equation alone, but with the β-term included. 2009 Curtin University of Technology and John Wiley & Sons, Ltd. (3) results on existence of steady flows that are small perturbation of the rest state; (4) results on existence of steady flows that are small perturbation of the Newtonian flow. STABILITY ANALYSIS Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering INFLUENCE OF WEIßENBERG NUMBER Next, we linearize the full system at a given state (τ 0 , u 0 , p 0 ), but still without convection. Finally, we study the full system. It turns out that this procedure is adequate, as we identify the possible mechanism leading to instability separately. Our finding is that the possible source for instabilities is the β-term, as it is a term of 0-th order with an indefinite sign. Convection may transport this instabilities. The limitation of our analysis of course is that we cannot investigate strong nonlinear interactions. However, numerical examples in the next section support our theoretical findings. Let us mention that without loss of generality we restrict ourselves to the case of two space dimensions. Stability analysis for the linearized Oldroyd system We start by studying stability properties of the Oldroyd system (1) linearized at u ≡ 0, that is, by neglecting the convective terms and the β-term: Problem 1 Find (τ , u, p) fulfilling summed up. Multiplying the second equation of Oldyfour by ξj and using the third equality of Oldyfour, the pressure can be written in the form: p̂ = τ (x , t) = R2 exp (−i ξ · x )τ̂ (ξ, t) dξ, R2 ∂ τ̂mj ˆ We ∂t + τ̂mj + i α(ξm ûj + ξj ûm ) = fs , ∂ û ξξ ξ Re ∂tj +i ξm τ̂mj +(1 − α)|ξ |2 ûj −i j m2 s τ̂sm = ĝj , (7) |ξ | ξm ûm = 0 where ĝ = fˆ + i ξ · (ξ · fˆ )/|ξ |2 . For simplicity, without loss of generality, in the sequel we will treat the case fs = 0 and f = 0 only. Considering now the vector of unknowns U := (τ̂xx , τ̂xy , τ̂yy , ûx , ûy ) , system (7) can be put in the form: ∂U + AU = 0, with ∂t 1 1 B I3 We We A = 2 1 C (1 − α)|ξ | I 2 Re Re exp (−i ξ · x )û(ξ, t) dξ, R2 exp (−i ξ · x )p̂(ξ, t) dξ (8) where the submatrices B and C are defined as follows 2ξ1 0 B = i α ξ2 ξ1 0 2ξ2 i ξ1 ξ22 −ξ2 z −ξ1 ξ22 C = (9) 2 ξ12 ξ2 |ξ |2 −ξ1 ξ2 ξ1 z p(x , t) = (6) Here, |ξ | is the length of the wave vector and z = ξ12 − ξ22 . The eigenvalues of the matrix A are given by: u(x , t) = |ξ |2 Thus, by eliminating the pressure field in Eqn 5, one obtains the following system for the unknowns τ̂ and û ∂τ We ∂t + τ − 2αD(u) = fs , ∂u (3) Re ∂t − 2(1 − α)divD(u) − div τ + ∇p = f , div u = 0 with appropriate initial and boundary conditions. The right-hand side fs in the stress transport equation is introduced for the matter of generality. We (formally) perform a Fourier transform of this system and its unknowns in R2 : i ξ · fˆ + ξ · (ξ · τ̂ ) λe4,5 This transformation leads to the following system: ∂ τ̂mj We ∂t + τ̂mj + i α(ξm ûj + ξj ûm ) = fˆs , (5) ∂ û Re ∂tj + i ξm τ̂mj + (1 − α)|ξ |2 ûj − i ξj p̂ = fˆj , ξm ûm = 0 where fˆs and fˆ are the corresponding Fourier transform of the right-hand sides fs and f , m, j ∈ {1, 2}, respectively, and repeated indices are understood to be 2009 Curtin University of Technology and John Wiley & Sons, Ltd. 1 (1 − α)|ξ |2 , λe3 = , We Re √ 1 1 (1 − α)|ξ |2 + ± e = 2 Re We λe1,2 = (4) (10) with e = 1 (1 − α)|ξ |2 − Re We 2 − 4α|ξ |2 ReWe (11) The first three eigenvalues are clearly real and positive, so for the stability of the system (8) it remains to show that the other two eigenvalues λ4 , λ5 Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj 659 660 N. SCURTU AND E. BÄNSCH Asia-Pacific Journal of Chemical Engineering have positive real parts. To do so, two cases can be distinguished. First, we observe that for e < 0, the real part R(λe4,5 ) = (1 − α)|ξ |2 /Re + 1/We is positive and thus we have stability in this case. Secondly, if e ≥ 0 then λ4,5 ∈ R and one can easily check that these two eigenvalues are positive. Thus, the linearized Oldroyd system is unconditionally stable. Influence of the β-term on the stability of the constitutive equation The next step in examining the stability properties of the Oldroyd system is to study the influence of the βterm. To this end, we first look at the stress problem alone. Problem 2 Given a solenoidal velocity field u, solve in × R+ the following equation for the stress tensor τ 2α 1 ∂τ + (u · ∇)τ + βa (τ, ∇u) + τ = D ∂t We We (12) supplemented by the boundary condition τ|− = τ on the inflow boundary and the initial condition τ|t =0 = τ0 . Neglecting the convective term (u · ∇)τ , the stress equation with given divergence-free velocity field reduces to a simple initial value problem. By representstress tensor in the compact form τ = T ing the symmetric τxx τxy τyy , the stress equation can be written as: ∂τ + Aτ = R, A = ∂t We −1 − 2aD 2(−aDxy + Wyx ) 0 xx −1 −(aDxy + Wyx ) We −aDxy + Wyx 0 −2(aDxy + Wyx ) We −1 + 2aDxx (13) with R some right-hand side depending on u and D = Dxx Dxy 0 Dxy ,W = −Dxx Wyx −Wyx 0 (14) Because of the incompressibility condition (1)3 , the component Dyy of the rate of deformation tensor is equal to −Dxx . The three eigenvalues of matrix A depend on the velocity field and the Weißenberg number in the following way: λs = We −1 , 2 + D2 ) − W 2 λ± = We −1 ± 2 a 2 (Dxx xy yx In the case of a = 0, this condition leads to an upper stability limit for the Weißenberg number: −1 2 + D2 ) − W 2 We cr = R 2 a 2 (Dxx xy yx (16) which is caused by the β-term. More specifically, in the case of a = ±1, this stability bound implies the following: if their exists any point x in the computational domain, where the velocity field satisfies det(∇u(x )) < 0, then the stress equation can be stable only up to Wecr = (2 | det(∇u)|)−1 (17) This limit is confirmed by our numerical tests in Section on Numerical Study, see also Ref.[19]. At this point, we remark that the restriction imposed in Ref.[14] on the relaxation time for the test example used there is nothing else but the stability restriction det(∇u(x )) > 0 found in this work. In the case of a = 0, the radicand in Eqn (15) is always negative, so no stability limit for the Weißenberg number exists. This result is in agreement with the conclusion of Jourdain et al.[20] , who have proven formally exponential convergence to equilibrium for the case a = 0 of Eqn (12). Constitutive equation with β-term and convection Stability limits can also exist for the full stress constitutive equation with convection. However, this situation is more involved. In this case, we have to look at the possible development of instabilities along streamlines, on which the stress equation can be again treated as an initialboundary problem. It turns out that along streamlines, which leave the domain in a finite time, instabilities cannot blow up, but may rather disturb the solution with at most finite amplitude. This is because these perturbations have only a finite time to grow during traversing of streamlines through a region, where the stability condition is violated. However, if there exists a stagnation point of the flow, where det(∇u) < 0 (in the case of a = ±1), or if closed streamlines intersect such a region, then instabilities along these streamlines arise for Weißenberg numbers beyond the critical value. (15) Influence of the nonlinearity β on the stability of the Oldroyd system The stable regime is again given by the condition that the real parts of the eigenvalues of matrix A are positive. In this subsection, the full Oldroyd system but without the convective terms is considered: 2009 Curtin University of Technology and John Wiley & Sons, Ltd. Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering INFLUENCE OF WEIßENBERG NUMBER Problem 3 Solve in × R+ the following system of equations ∂U + A5 U = 0, with ∂t 1 B + B̃ A We A5 = 1 1 − α |ξ |2 I 2 Re C Re 1 2α ∂τ ∂t + βa (τ, ∇u) + We τ − Re divD = 0 2(1 − α) divD − 1 div τ + 1 ∇p = f , (18) ∂u ∂t − Re Re Re div u = 0 with boundary condition u| = u , and initial conditions u|t =0 = u0 , τ|t =0 = τ0 . The matrices B and C are those defined in Eqn (9), the matrix A is analogous to Eqn (13), but with D and W replaced by D(u 0 ) and W (u 0 ), and B̃ is defined by: 0 0 ξ1 − (1 + a)τxy ξ2 −2aτxx 0 (1 − a)τxy ξ1 λξ = u(x , t) = u (x ) + u(x , t), 0 (19) Neglecting higher order terms, this splitting leads to the following linear system of equations for the perturbations: Problem 4 Solve in × R+ for τ , u and p, fulfilling the following system of equations: ∂τ 1 τ + βa (τ , ∇u 0 ) + βa (τ 0 , ∇u) + We ∂t 2α D(u) = 0, − We ∂u − 2(1 − α) divD(u) − 1 div τ + 1 ∇p = 0, ∂t Re Re Re div u = 0 (20) with boundary condition u | = 0 and initial conditions u |t =0 = u 0 and τ |t =0 = τ 0 . We would now like to perform a Fourier transform again like in the previous section. But this is only feasible, if τ 0 , ∇u 0 , which act like coefficients for the unknowns in Eqn (20), would be constant in space, which is an unrealistic assumption in general. However, in performing a Fourier transform one might think of not too oscillatory functions τ 0 , ∇u 0 or local properties. Thus, consider now formally the Fourier transforms (τ̂ , û, p̂) of the unknown fields (τ , u, p). Removing again the pressure from the set of equation like in Eqn (7), Eqn (20) can be written as: 2009 Curtin University of Technology and John Wiley & Sons, Ltd. (22) The case a = 1 In the case of a = 1, for the Oldroyd-B model, one eigenvalue of matrix A5 is: τ (x , t) = τ 0 (x ) + τ (x , t), p(x , t) = p 0 (x ) + p(x , t) 1 [(1 − a)τ 0 − (1 + a)τ 0 ]ξ yy xx 1 2 0 0 −2aτyy ξ2 − (1 + a)τxy ξ1 0 0 B̃ = i 1 [(1 − a)τxx − (1 + a)τyy ]ξ2 2 0 (1 − a)τxy ξ2 Let (τ 0 , u 0 , p 0 ) be a stationary solution of Problem 3, and (τ , u, p) be a small perturbation of the stationary solution such that (21) 1−α 2 |ξ | Re (23) and the other four eigenvalues can be obtained from the solutions λ of the polynomial equation: 1−α 2 |ξ | − λ Re 2 2 1 1 0 2 × 4 det (∇u ) + |ξ | + −λ −λ We We 1 −λ We α |ξ |2 0 0 0 |ξ |2 − ξ1 2 τxx + 2ξ1 ξ2 τxy + ξ22 τyy We Re 1 + − λ 2Dxx ξ1 ξ2 − Dxy z + Wyx |ξ |2 We 2 0 0 0 × −ξ1 ξ2 τxx + z τxy + ξ1 ξ2 τyy Re 4α − [Dxx z + 2Dxy ξ1 ξ2 ]2 ReWe 4 0 [Dxx z + 2Dxy ξ1 ξ2 ] [Dxy |ξ |2 − Wyx z ]τxy + Re 0 + [Dxx ξ1 + (Dxy + Wyx )ξ2 ]ξ1 τxx 0 = 0 (24) + [−Dxx ξ2 + (Dxy − Wyx )ξ1 ]ξ2 τyy × with z = ξ12 − ξ22 as in Eqn (9). It is interesting to observe that the left-hand side of this equation comprises only terms of order two and four in |ξ |. For small wave vectors, the fourth-order terms can be Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj 661 662 N. SCURTU AND E. BÄNSCH Asia-Pacific Journal of Chemical Engineering neglected compared with the second-order terms. Thus, one obtains a reduced equation 1 −λ λ We 4 det (∇u ) + 0 1 −λ We 2 = 0 (25) which possess the solutions λ0 = 0 and λs , λ± as in the pure stress constitutive equation (15). This means, the stability properties of system (20) for case a = ±1 are, as expected, not better than those of the pure stress constitutive equation (12). The case a = 0 In the case of a = 0, the matrix A5 has the two eigenvalues λξ and λs as above. The other three eigenvalues are solutions of the following equation, which is of third order with respect to λ: 2 1 1−α 2 −λ |ξ | − λ We Re 1 1−α 2 2 + 2Wyx |ξ | − λ + −λ Re We α |ξ |2 1 z − τxy ξ1 ξ2 + (τyy − τxx ) We 2 4 Re 1 2 1 + Wyx τxy z + (τyy − τxx )ξ1 ξ2 =0 Re (26) It seems difficult to give explicit solutions for Eqns (24) and (26) for general wave vectors. However, by fixing the direction of the wave vector to, for example, ξ = (ξ1 , 0), ξ = (0, ξ2 ) or ξ = (ξ1 , ±ξ1 ), one can simplify the analysis of the eigenvalues for a concrete problem for which the velocity and stress fields are known. This will be carried out for the examples in Section on Numerical Study. Influence of the convective terms on the stability of the Oldroyd system Applying the splitting Eqn (19) to the full Oldroyd system (1), after neglecting higher order terms, the following linearized problem for the perturbed fields is obtained: Problem 5 Solve in × R+ for τ , u and p the following system of equations ∂τ + (u 0 · ∇)τ + (u · ∇)τ 0 + β (τ , ∇u 0 ) a ∂t 1 2α 0 + βa (τ , ∇u) + We τ − We D(u) = 0, ∂u + (u 0 · ∇)u + (u · ∇)u 0 − 2(1 − α) divD(u) ∂t Re 1 div τ + 1 ∇p = 0, − Re Re div u = 0 2009 Curtin University of Technology and John Wiley & Sons, Ltd. with boundary condition u | = 0 and initial conditions u |t =0 = u 0 and τ |t =0 = τ 0 . Introducing again the Fourier transforms of the perturbed fields, Problem 5 can be written as an initial value problem of type (21) with the matrix A5 now replaced by Â5 = A5 − i (ux0 ξ1 + uy0 ξ2 ) I5 + D5 (27) where A5 is defined by Eqn (21) and D5 = 0 0 0 ∂τxx ∂x 0 B̂ , B̂ = ∂τxy D̂ ∂x 0 ∂τyy ∂x 0 ∂ux ∂x D̂ = 0 ∂uy ∂x 0 ∂τxx ∂y 0 ∂τxy , ∂y 0 ∂τyy ∂y ∂ux0 ∂y (28) ∂uy0 ∂y The convective terms (u 0 · ∇)τ and (u 0 · ∇)u contribute to the matrix Â5 only through the purely imaginary term i (ux0 ξ1 + uy0 ξ2 ) times identity I5 and thus does not influence the real parts of the eigenvalues. The other terms arising from linearizing the convective terms, namely (u · ∇)τ 0 and (u · ∇)u 0 collected in D5 , however, may influence the real parts of the eigenvalues. Thus, the real parts of the eigenvalues of matrix Â5 may differ from the real parts of the eigenvalues of matrix A5 . The results of the previous stability analysis are summarized in Table 1. For the case a = ±1, the critical Weißenberg number Wecr is (neglecting fourth order terms) Wecr = (2 | det(∇u 0 )|)−1 , computed from Eqn (25). By P we denote the set of points x ∈ that are either stagnation points for the flow or are located on streamlines that form closed curves in the domain . For the case a = 0, the critical Weißenberg number Wecr0 is defined as the minimum value of We for which Eqn (26) has a solution with negative real part in some point of the domain for some wave vector. By We we include the effect of matrix D5 , which we consider as a perturbation to the eigenvalues without this term. For a concrete given example one may compute the exact eigenvalues of matrix Â5 , though. NUMERICAL STUDY In this section, we use numerical simulations to validate our formal stability analysis of the previous section. To this end we will define three test examples. The numerical method we are using was developed in Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering INFLUENCE OF WEIßENBERG NUMBER Table 1. Results of the stability analysis for the different cases.P is the set of points that are either stagnation points for the flow or are located on streamlines that form closed curves in the domain . We0cr is defined as the minimum value of We for which Eqn (26) has a solution with negative real part. By We we include the effect of matrix D5 , which we consider as a perturbation to the eigenvalues without this term . Linear system Stress eqn./wo conv. Stress eqn./w conv. Oldroyd/wo conv. Full Oldroyd a=0 Uncond. stable Uncond. stable Uncond. stable Unstable if We > Wecr0 Unstable if We > Wecr0 + We a = ±1 Uncond. stable Unstable if det(∇u 0 (x )) < 0, x ∈ and We > Wecr Unstable if det(∇u 0 (x )) < 0, x ∈ P and We > Wecr Unstable if det(∇u 0 (x )) < 0, x ∈ and We > Wecr Unstable if det(∇u 0 (x )) < 0, x ∈ P and We > Wecr + We Ref.[19]. Here, we just sketch the basic ideas. For further reference, see the citations above. The method is based on a splitting strategy via the fractional step θ -scheme for time discretization. The latter one goes back to Bristeau et al.[21] . By this algorithm, the three major numerical difficulties of the Oldroyd system are decoupled: the treatment of the nonlinearity in the momentum equation (given by the velocity transport term), the solenoidal condition, and the stress transport term in the constitutive equation. By the operator splitting, one reduces the Oldroyd system to three considerable simpler subproblems: a linear self-adjoint Stokes like problem, a nonlinear Burgers problem, and a transport problem for the stress tensor. The hyperbolic stress constitutive equation is discretized in space by the discontinuous Galerkin method with piecewise linear functions. For the velocity and pressure fields, we have chosen the LBB-stable TaylorHood element. This choice of finite elements is justified by an inf-sup condition, see Ref. [22]. The algorithm was implemented via the finite element package Alberta, see Ref. [23]. Numerical tests were performed for different refinement levels of the underlying mesh starting from a macro triangulation of the domain consisting of two triangles. For refinement level n, the mesh width h is then h = 2−n . All simulations were run as time evolutions. They were started with an interpolant of an exact, steady state solution as initial condition. Then we followed the time evolution of the discrete solution. In case of stable solutions, they converged to a steady discrete state not too far off the exact solution. In case of instabilities, the solutions departed from this exact solution and blew up for time t → ∞. To further assess the quality of the numerical solutions, we introduce the experimental order of convergence (EOC), defined by EOC := 1 Errh , log 2 Errh/2 Errh = lim vex − vh (·, t) t→∞ (29) 2009 Curtin University of Technology and John Wiley & Sons, Ltd. Here, Errh is the error of an exact steady solution vex and the numerical solution vh corresponding to a mesh with width h. v may stand for u, p, or τ . The EOC gives the computational rate of convergence for the steady discrete solution to the exact one with respect to the mesh size h. For the pressure and velocity fields, we chose the L2 - and H 1 -norm, respectively, and for the stress components we considered the L2 -norm. Note that for the choice of elements used in our computations, that is, globally continuous, piecewise quadratics for the velocity, globally continuous, piecewise linears for the pressure, and discontinuous, piecewise linears for the stresses, the optimal order is 2 (in case of the H 1 -norm) and 3 (for the L2 -norm) for the velocity, 1 (in case of the H 1 -norm) and 2 (for the L2 -norm) for the pressure and 2 for the stresses. Example 1 We choose the right-hand sides fs and f in such a way that the exact velocity, stress components, and pressure are: uex = sin( π2 x ) sin( π2 y) , π x ) cos( π y) , cos( 2 2 xx , xy, yy = pex = cos( π2 x ) sin( π2 y) τex (30) The reason for this choice is that the eigenvalues from Section on Stability Analysis, which can be computed explicitly here, are subcritical in some parts of the domain, whereas supercritical in others. The computational domain is the unit square, = [0, 1]2 . We will now study the behavior of the computational solutions in a similar order like for the stability analysis. We start with the first example, where one might encounter instability. Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj 663 664 N. SCURTU AND E. BÄNSCH Stress equation without convective term for a=1 In the case of a = 1, the eigenvalue λ− from Eqn (15) for this example is given by 1 We π π π π − π cos2 ( x ) sin2 ( y) − sin2 ( x ) cos2 ( y) (31) 2 2 2 2 λ− = For increasing values of We, this eigenvalue is the first passing through zero and this takes first place in the corner point with coordinates (x , y) = (0, 1) at the critical Weißenberg number Wecr = 1/π ≈ 0.3183. Beyond this critical state, the part of the domain, where the eigenvalue λ− is negative, becomes increasingly larger. The respective eigenvalues and eigenvectors corresponding to the A-matrix for this example problem in the corner point (0, 1) are: 1 1 1 , λ+ = + π, λ− = − π, We We We e0 = [0 1 0], e+ = [0 0 1], e− = [1 0 0] (32) λ0 = Consequently, the eigenvector corresponding to the eigenvalue λ− is the unit vector in the x direction, that Asia-Pacific Journal of Chemical Engineering is, λ− affects only the stress component τxx . Thus, when λ− becomes negative, then the component τxx of the stress tensor will first be affected by the instability. The instability in the other stress components τxy and τyy appear much later. This behavior is clearly confirmed by the computational results. Figure 1 shows a perturbation emanating from point (0, 1), and the part of the domain with instability is broadening with increasing Weißenberg number. We note that this phenomenon is essentially independent of the mesh width. Stress equation with convection for a = 1 As just seen, the stress equation without convection and a = 1 admits stable solutions only for subcritical Weißenberg numbers. For the stress transport equation with convection, the stress is transported along streamlines that are depicted for the present example in Fig. 2. The inflow boundary for Example 1 is the line y = 0 and x = 1 is the outflow boundary. All interior points of the domain lie on streamlines that leave the domain, thus for these points the stress field is unperturbed, see Fig. 2, which shows a very different behavior than for the above case without convection, Section on Stress Equation without Convective Term for a = 1. Figure 1. Isolines of the solution τxx of the stress equation without convective term for a = 1 and mesh size h = 2−6 . Stable solution for We = 0.3 and unstable solutions for We ∈ {0.32, 0.4, 0.5} at time t = 10. For t → ∞, the numerical solutions for the three last cases blew up. This figure is available in colour online at www.apjChemEng.com. Figure 2. Velocity field (left) and isolines for the steady solution τxx (right) of the stress equation with convection for a = 1 and mesh size h = 2−6 at We ∈ {0.4, 0.5, 0.6}. The stress field is only slightly perturbed on the top of the domain. This figure is available in colour online at www.apjChemEng.com. 2009 Curtin University of Technology and John Wiley & Sons, Ltd. Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering INFLUENCE OF WEIßENBERG NUMBER Oldroyd system without β-term, without convection Just like in the stability analysis from Section on Stability Analysis for the Linearized Oldroyd System, for the Oldroyd system without the β-term no stability limits exists. Without showing corresponding results we note that the computations confirm convergence to the steady state solution for any value of the Weißenberg number. Oldroyd system with β-term, without convection for a = 1 In the case of the Oldroyd system without convection for a = 1, the influence of the β-term is similar to the stress equation without convection for a = 1 with given velocity field, see Section on Influence of the β-term on the Stability of the Constitutive Equation. Because for the pure stress equation the critical state appears at the corner (0, 1) first, we take a look at the eigenvalues and eigenvectors of the A5 -matrix from Eqn (21) of the Oldroyd system at this point. At the corner point (0, 1) the respective velocity field, the rate of deformation tensor, the spin tensor, and the stresses calculated from Eqn (30) are: ux = uy = 0, Dxx = π , 2 π Dxy = 0, Wyx = 0, det(∇u) = − , 2 τxx = τxy = τyy = 1 (33) One eigenvalue of matrix A5 (21) is λξ (defined in Eqn (23)) and is always positive. Computing the other four eigenvalues from Eqn (24) explicitly is not a simple task in general. Therefore, we make a particular, but hopefully representative choice of the wave vectors, for example, ξ = (ξ1 , 0) or ξ = (0, ξ2 ), where ξ1 and ξ2 are arbitrary. For this choice, it is easy to compute the 2009 Curtin University of Technology and John Wiley & Sons, Ltd. exact eigenvalues of the perturbed solution at the corner point (0, 1) corresponding to the above wave vectors. The eigenvalues determined by Eqn (24) are given by: 1 ± π, We |ξ |2 1 1 (1 − α) + = 2 Re We λ1,2 = λ3,4 ! 2 " 2 2 " |ξ | 1 |ξ | 1 ±# (1 − α) + −1 −4 . Re We We Re The real parts of the eigenvalues λ1 and λ3 are always positive, the real part of λ4 is positive as long as We < 1, but the real part of the eigenvalue λ2 becomes negative already at Wecr = 1/π . The eigenvector corresponding to the eigenvalue λ2 is given by e2 = [1 , 0 , 0 , 0 , 0]. Thus, similar as in the analysis of the pure stress equation without the convective term, the stress component τxx is first affected by instabilities, when λ2 becomes negative. h L2 In Fig. 3, the time evolution of the error τxx − τxx is shown for two subcritical values of the Weißenberg number We ∈ {0.3, 0.31} and for the supercritical value We = 0.32. For subcritical values, the time evolution converged to the stationary (discrete) state, but for supercritical values no convergence was achieved and the computations eventually blew up. The stress component τxx is plotted in Fig. 4 for We = 0.32 along the diagonal y = 1 − x of the square domain. One can observe the instability appearing in the region near the corner and its growth in time. 5 x 10−3 We = 0.32 We = 0.31 We = 0.3 4.5 ||τxx − τhxx||L2 − error norm In contrast to the interior points, we see some perturbations at the line y = 1. The solution of the stress transport equation with convection is steady but perturbed. The perturbation exists only along the streamline y = 1 in the upper element layer and grows with growing Weißenberg number. In the following subsections, we describe the results of simulations for the Oldroyd system both without and with β-term, as well as without and with convective terms. The numerical runs were carried out for different values of the Weißenberg number, while keeping the other parameters constant at the values Re = 1.0 and α = 0.89. 4 3.5 3 2.5 2 1.5 1 0.5 0 0 20 40 60 80 100 t − time h 2 Figure 3. Time evolution of the error τxx − τxx L for the Oldroyd-B system without the stress transport term for h = 2−6 at different We. This figure is available in colour online at www.apjChemEng.com. Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj 665 666 N. SCURTU AND E. BÄNSCH Asia-Pacific Journal of Chemical Engineering Full Oldroyd-B system, a = 1 In the case of the full Oldroyd system with a = 1, we analyze the eigenvalues of matrix Â5 (defined in Eqn (27)) for the particular choice of the wave vectors of the perturbed solution, ξ = (ξ1 , 0) and ξ = (0, ξ2 ), as before. For the wave vector ξ = (ξ1 , 0), the eigenvalues of matrix Â5 (defined in Eqn (27)) are: instable regime for each We, whereas for large wave vectors the stable regime is determined by the condition $ % 1 π Re π 2 , 1− |ξ | > Re and We < min 2 π 2 |ξ |2 (37) For the wave vector ξ = (0, ξ2 ) of the perturbed solution, the eigenvalues of matrix Â5 are given by 1−α 2 π 1 |ξ | + , λ1,2 = ± π, Re 2 We ! " 2 2 2 2 " |ξ | |ξ | 1 π 1 π |ξ | π 1 1 1 + − ± # (1 − α) + − −1 − = (1 − α) −4 (34) 2 Re We 2 Re We 2 We Re 2 We λξ = λ3,4 1−α 2 π |ξ | − , λ1,2 = Re 2 ! " 2 2 " |ξ |2 |ξ |2 1 1 1 π 1 π |ξ | # (1 − α) = (1 − α) −4 + + ± + + −1 + 2 Re We 2 Re We 2 We Re λξ = λ3,4 1 ± π, We π 1 (38) 2 We The real parts of the eigenvalues λξ , λ1 and λ3 are always positive, while the real part of λ2 becomes negative already at We = 1/π , whereas the real part of the eigenvalue λ4 is positive as long as 2 |ξ | π 1 1 −1 − >0 (35) We Re 2 We Again, the eigenvalue λξ is positive only for large wave vectors satisfying |ξ |2 > π Re/2(1 − α). The real part of the eigenvalues λ1 and λ3 are always positive, the real part of λ2 becomes negative already at We = 1/π , but the real part of the eigenvalue λ4 is positive as long as: For small wave vectors, if |ξ |2 < Re π/2, the real part of the eigenvalue λ4 is negative for each value of We, which interestingly means that we have a long wave instability. We note that as in our example Re = 1 and thus Re π/2 ≈ 1.57 and = [0, 1]2 , these long wave instabilities do not show up in our computations. For large wave vectors, if |ξ |2 > Re π/2, the real part of the eigenvalue λ4 is positive only for (39) We < 1 − π Re 2 |ξ |2 (36) Thus, the eigenvalue analysis corresponding to the wave vector ξ = (ξ1 , 0) shows that for small wave vectors, |ξ |2 < Re π/2, the Oldroyd-B system is in the 2009 Curtin University of Technology and John Wiley & Sons, Ltd. 2 |ξ | π 1 1 −1 + >0 We Re 2 We and accordingly for values of the Weißenberg number that satisfy π Re (40) We < 1 + 2 |ξ |2 Therefore, the eigenvalue analysis corresponding to the wave vector ξ = (0, ξ2 ) yields that the stable regime is given by the condition $ % 1 π Re π Re and We < min , 1+ |ξ |2 > 2 1−α π 2 |ξ |2 (41) From conditions (37) and (41), one can conclude that the solutions are stable only for large wave vectors and Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering INFLUENCE OF WEIßENBERG NUMBER 1.2 4 t = 50 t = 80 exact solution 0.8 0.6 error norm τxx − stress component 1 0.4 0.2 3.5 ||τxx − τhxx||L2 3 ||τxy − τhxy||L2 2.5 ||τyy − τhyy||L2 ||u − uh||L2 2 1.5 1 0 −0.2 x 10-3 0.5 0 0.2 0.4 0.6 0.8 1 0 position along y = 1x 0 5 10 15 20 25 t−time Numerical values of τxx for the Oldroyd-B system without the stress transport term along the line y = 1 − x at We = 0.32. This figure is available in colour online at www.apjChemEng.com. Figure 4. More precisely, for small wave vectors Wecr = 0 and for large wave vectors Wecr = 1/π if α ≥ 1/π and Wecr = α if α < 1/π . Similarly as in the case of the pure stress equation with given velocity field, the corner point (0, 1) may cause instability, which is propagated along the streamline traversing this point. The numerical simulations show stable steady solutions up to We ≈ 0.4, which is slightly above the critical value Wecr . Beyond this value, however, for instance for We ∈ {0.5, 0.6}, the perturbation in the τxx -component, emanating from the corner point, is growing along the line y = 1 and eventually blowing up. Because the Oldroyd system is a nonlinear coupled system for the unknown fields, the instability in τxx eventually also affects the other components of the solution, which can be seen in Fig. 5. Here, the growth of the error norms in time is displayed for We = 0.5 and We = 0.6 for the mesh with size h = 2−6 and time step size t = 0.01. The stress component τxx is plotted in Fig. 6 at different time instants for We = 0.6 along the line y = 1 − h/2. One can see the deviation of the numerical values of τxx from the exact solution near the streamline passing the stagnation point. Neither further mesh refining nor reduction of the time step led to better convergence results, confirming that this phenomenon is inherent in the equation. To show that the perturbation in the upper element layer is only caused by the instability arising at the 2009 Curtin University of Technology and John Wiley & Sons, Ltd. Oldroyd-B system for h = 2−6 at We = 0.6. This figure is available in colour online at www.apjChemEng.com. 1 0.5 τxx − stress component small Weißenberg numbers, which satisfy We < 1 ifα ≥ 1 π Re π π and (42) |ξ |2 > We < 1 − π2 Re2 ifα < π1 2 1−α |ξ | Figure 5. Time evolution of the error L2 -norms for the full 0 −0.5 −1 −1.5 −2 t = 15 t = 25 exact solution −2.5 −3 0 0.2 0.4 0.6 0.8 1 x − position along y = 0.9961 Numerical values of τxx for the OldroydB system along the line y = 1 − h/2 for h = 2−7 at We = 0.6. This figure is available in colour online at www.apjChemEng.com. Figure 6. stagnation point (x , y) = (0, 1), and is not influenced by some boundary effect, calculations were performed also in the domain = [−1, 1] × [0, 2]. In Fig. 7, the velocity field and the time evolution of the stress component τxx are presented for We = 0.6 for a mesh with size h = 2−5 and time step width t = 0.01. Stress transport problem for a = 0 Now, let us clarify what happens in the case of a = 0. For the pure stress transport, the eigenvalues (15) are in the stable regime for each value of the Weißenberg number. Consequently, no stability limit exists in this case. Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj 667 668 N. SCURTU AND E. BÄNSCH Asia-Pacific Journal of Chemical Engineering τxx at t = 30 u at t = 100 τxx at t = 100 Figure 7. Velocity field and propagation of the τxx perturbation at We = 0.6 for the full Oldroyd-B system; = [−1, 1] × [0, 2]. This figure is available in colour online at www.apjChemEng.com. Table 2. Errors and EOCs for Example 1 for a = 0 and We = 1000; ‘Ref.’ is the refinement level of the mesh. xy yy Ref. τ xx − τhxx L2 EOC τ xy − τh L2 EOC τ yy − τh L2 EOC 4 5 6 7 6.272069e-04 1.464079e-04 3.497572e-05 8.504148e-06 2.10 2.07 2.04 – 5.865086e-04 1.410928e-04 3.429843e-05 8.415914e-06 2.06 2.04 2.03 – 5.677649e-04 1.382838e-04 3.388190e-05 8.354178e-06 2.04 2.03 2.02 – This fact is reflected in Table 2, which contains EOC tests for We = 1000. The computations converged to a steady state for any value of the Weißenberg number. Oldroyd system without convection for a = 0 Although the stress transport equation for a = 0 is unconditionally stable, in the Oldroyd system without convection instabilities may arise, as it can be seen from the eigenvalue analysis of matrix A5 defined in Eqn (21). Like above, we simplify the explicit analysis by fixing the wave vector to, for example, ξ = (ξ1 , 0) and ξ = (0, ξ2 ), where the components ξ1 and ξ2 are arbitrary. At the corner (x , y) = (0, 1), where relations (33) hold, Eqn (26) corresponding to the wave vector ξ = (ξ1 , 0) or ξ = (0, ξ2 ), admits a solution equal to 1/We and the other two solutions satisfy 1−α 2 1 1 |ξ |2 2 + |ξ | + = 0 (43) λ −λ We Re We Re The two solutions of Eqn (43) & λ1,2 = b ± 1 b= 2 1 |ξ |2 , b2 − We Re 1 1−α 2 + |ξ | We Re For the wave vectors ξ = (ξ1 , ±ξ1 ), Eqn (26) admits a solution equal to 1/We and the other two solutions satisfy the quadratic equation 2 1−α 2 |ξ | 1 1 + |ξ | + ∓1 = 0 λ −λ We Re We Re (45) with roots 2 & 2 |ξ | 1 ∓1 , =b+ − We Re & 2 |ξ | 1 ± 2 λ2 = b − b − ∓1 We Re λ± 1 b2 (46) − The solutions λ− 1 and λ2 corresponding to the wave vector ξ = (ξ1 , −ξ1 ) are always positive, but those corresponding to the wave vector ξ = (ξ1 , ξ1 ), namely, + λ+ 1 and λ2 , are both positive only as long as We ≤ 1.0. The eigenvector of matrix A5 corresponding to the eigenvalue λ+ 2 from Eqn (46) for the wave vector ξ = (ξ1 , ξ1 ), is given by −1 1 α + −1 − λ2 = 2i ξ1 0 We We −1 1 α + −1 − λ2 −1 1 −2i ξ1 We We e2+ (44) (47) are always positive. 2009 Curtin University of Technology and John Wiley & Sons, Ltd. Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering Thus, the perturbation will show up first in the stress componentsτxx and τyy . By means of the previous analytical considerations one can explain the instabilities arising in the Oldroyd system, without and with the convective terms, for a = 0, when We is increasing. In the case of the Oldroyd system without convective terms, convergence to a steady state was obtained computationally up to Wecr0 = 1.0. At We = 1.1 convergence was achieved only for large mesh sizes, for example, h ∈ {2−3 , 2−4 }, but not for finer meshes. For supercritical Weißenberg numbers, the perturbation appeared first in the stress components τxx and τyy , and slightly later in the τxy component, which is in agreement with the form of the eigenvector e2+ in Eqn (47). In Fig. 8, the growth of the perturbation in the stress component τxx at We = 1.1 is shown at different times. The perturbation is propagated in all directions of the computational domain, emanating from the corner point (x , y) = (0, 1). INFLUENCE OF WEIßENBERG NUMBER Full Oldroyd system for a = 0 Let us turn to the case of the full Oldroyd system with a = 0 and Â5 defined by Eqn (27). The eigenvalue analysis corresponding to the wave vectors ξ = (ξ1 , 0) and ξ = (0, ξ2 ) shows that only for large wave vectors satisfying |ξ |2 > π Re/2 the Oldroyd system is in the stable regime for each We. For the wave vectors ξ = (ξ1 , ±ξ1 ), two eigenvalues of matrix Â5 are λs (15), and the other three satisfy the cubic equation 1−α 2 π 1 −λ |ξ | + − λ We Re 2 1−α 2 π 1−α 2 |ξ | − − λ + |ξ | − λ Re 2 Re |ξ |2 α ∓1 =0 (48) We Re The numerical simulations show convergence up to We ≈ 1.4. Likewise, as for the full Oldroyd-B problem, Figure 8. Example 1, propagation of the perturbation in τxx for the Oldroyd system without the convective terms for a = 0 and We = 1.1. This figure is available in colour online at www.apjChemEng.com. Figure 9. Example 1, propagation of the perturbation in the pressure for the full Oldroyd system for a = 0 and We = 1.5, h = 2−6 . This figure is available in colour online at www.apjChemEng.com. 2009 Curtin University of Technology and John Wiley & Sons, Ltd. Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj 669 N. SCURTU AND E. BÄNSCH Asia-Pacific Journal of Chemical Engineering 200 −5 −6 180 −6 −6 −5 −3 −7 −7 −5 0 0 −3 −5 −6 −6 −5 −7 −7 −3 160 140 −6 −3 −5 −6 0 6 80 10 8 6 4 2 0 −2 We = 0.1 We = 0.109 We = 0.125 We = 0.2 We = 0.5 0 − − 3 −65 40 As a second example for the full Oldroyd system, we define = [−1, 1]2 and the following exact solution: u = (4 (2 − r 2 ) y , −4 (2 − r 2 ) x ), ex 2 r = x + y 2, (49) xx , xy, yy = 1 + exp(−2) − exp(−r 2 ), τex pex = exp(−r 2 ) −8 0 7 Example 2 −6 −6 −3 100 0 −5 0 120 60 −4 −7 −6 −7 −5 in the present case, because of the convective terms, the velocity and stress fields are transported along the streamlines. Thus again, along the streamlines that leave the computational domain no instabilities arise but on the streamline containing the stagnation point, the negative eigenvalue λ+ 2 defined in Eqn (46) and corresponding to the wave vector ξ = (ξ1 , ξ1 ), gives rise to instabilities. In Fig. 9, the growth of the perturbation in the pressure field at We = 1.5 is shown at different times. It is obvious that the perturbation appears first in the upper element layer on the boundary streamline downstream the stagnation point. real(λ−) 670 0.2 0.4 0.6 0.8 1 1.2 1.4 radial position Figure 10. Example 2, eigenvalue λ− on the line y = x, x ≥ 0 at different We. This figure is available in colour online at www.apjChemEng.com. 20 − 6 −5 −7 20 40 −3 0 −7 60 80 0 −3 −5 6 6 100 120 −5 −6 −7 140 160 −7 −6−5 180 200 Example 2, isolines of Rλ− at We = 0.5. This figure is available in colour online at www.apjChemEng.com. Figure 11. Again the right-hand sides were chosen in such a way that the above functions are solutions of the Oldroyd–B system. One key feature of this solution is the fact that the streamlines consist of circular arcs or closed circles. The critical Weißenberg number defined by relation (16) and corresponding to the velocity field (49) is approximatively Wecr ≈ 0.109. Figure 10 shows the real part of the eigenvalue λ− on the diagonal line y = x and x ≥ 0 for different values of We. For We < 0.125 the points, where Rλ− < 0, are located only on streamlines that leave the domain. Therefore, only for We > 0.125, one can expect instabilities to appear. This critical value of the Weißenberg number is available for the stress transport equation with a = 1. However, for the full Oldroyd-B system, like as for Example 1, the stability limits can be different. Simulations of the Oldroyd-B system for this example still converged at We = 0.3, and gave good EOC values, as it is shown in Table 3. Table 3. EOC tests for Example 2 applied to the Oldroyd-B system for We = 0.3, t = 0.001. Ref. 2 3 4 5 Ref. 2 3 4 5 p − ph H 1 2.530e-01 1.149e-01 5.720e-02 2.841e-02 τ xx − τhxx L2 1.355e-02 2.625e-03 6.107e-04 1.531e-04 EOC 1.13 1.00 1.00 – EOC 2.36 2.10 1.99 – u − uh L2 9.840e-03 1.211e-03 1.562e-04 1.886e-05 xy τ xy − τh L2 1.240e-02 2.359e-03 5.562e-04 1.450e-04 2009 Curtin University of Technology and John Wiley & Sons, Ltd. EOC 3.02 2.95 3.04 – EOC 2.39 2.08 1.93 – u − uh H 1 2.296e-01 5.476e-02 1.402e-02 3.356e-03 yy τ yy − τh L2 1.427e-02 2.484e-03 5.652e-04 1.439e-04 EOC 2.06 1.96 2.06 – EOC 2.52 2.13 1.97 – Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering INFLUENCE OF WEIßENBERG NUMBER t = 9.0 u t = 11.5 Example 2: velocity field and time propagation of the τxx perturbation in the full Oldroyd-B system for We = 0.5. This figure is available in colour online at www.apjChemEng.com. Figure 12. τxx u τxy Figure 13. Example 3: velocity field and isolines of the stress components τxx and τxy at We = 100. This figure is available in colour online at www.apjChemEng.com. Table 4. EOC tests for Example 3 for the full Oldroyd-B system for We = 100 and t = 10−3 . Ref 3 4 5 6 7 Ref 3 4 5 6 7 p − ph H 1 6.836e+00 2.877e+00 1.020e+00 3.932e-01 1.709e-01 τ xx − τhxx L2 7.378e-01 1.078e-01 2.222e-02 4.158e-03 9.685e-04 EOC 1.25 1.49 1.37 1.20 – EOC 2.77 2.28 2.41 2.10 – u − uh L2 6.476e-02 1.009e-02 1.882e-03 2.738e-04 3.444e-05 xy τ xy − τh L2 4.268e-01 9.977e-02 1.773e-02 4.012e-03 1.009e-03 By increasing the Weißenberg number to 0.5, the system does not converge any more. In Fig. 12, one can see the time evolution of the perturbation in the stress component τxx on mesh with size h = 2−5 . The perturbation grows and is propagated along the streamlines, in the region of the negative eigenvalue (compare with Fig. 11), where the streamlines are closed curves. On the streamlines that leave the computational domain, no perturbation appear, even though some of them traverse the region of unstable regime. But along the closed streamlines that lie in the unstable region, the perturbation is propagated and grows in time. 2009 Curtin University of Technology and John Wiley & Sons, Ltd. EOC 2.68 2.42 2.78 2.99 – EOC 2.09 2.49 2.14 2.04 – u − uh H 1 1.240e+00 4.077e-01 1.071e-01 2.620e-02 6.711e-03 yy τ yy − τh L2 6.474e-01 1.507e-01 2.148e-02 4.316e-03 9.730e-04 EOC 1.61 1.93 2.03 1.96 – EOC 2.10 2.81 2.32 2.09 – Example 3 This third example is taken from Ervin and Ntasin.[24] The computational domain is an L-shaped domain given by = [−1, 1]2 \ [0, 1]2 . The velocity, stress and pressure fields are defined by y − 0.1 x − 0.1 u= − , r r r = (x − 0.1)2 + (y − 0.1)2 , τ = 2 α D(u), p = (2 − x − y)2 (50) Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj 671 672 N. SCURTU AND E. BÄNSCH The velocity and stress fields have singularities just exterior to the domain near the point (0, 0). The crucial feature of this example is the fact that the velocity field satisfies det(∇u) = 0 pointwise everywhere. According to Eqn (15), the real part of the eigenvalue λ− , corresponding to the stress equation, is therefore positive in the whole computational domain, and thus no stability limit for the Weißenberg number exists. Moreover, the streamlines neither comprise a stagnation point nor are closed curves in the whole domain . Hence, theoretically, one expects no instabilities for the Oldroyd problem. Indeed, the simulations converged to the steady state for every value of the Weißenberg number we tested, see Fig. 13. In Table 4, corresponding EOCs are presented for the Oldroyd-B model for We = 100. Asia-Pacific Journal of Chemical Engineering for arbitrary values of the Weißenberg number (for a = ±1). The eigenvalue analysis showed that even in the case of a = 0 instabilities can occur. The numerical tests confirmed this stability limits. For the full Oldroyd-B system, at least those stability limits arising in the pure stress constitutive equation, exist, if the flow domain comprises a stagnation point in the region where det(∇u) < 0, or if the streamlines are closed curves that intersect such a region. The numerical simulations confirmed this stability limit derived by the eigenvalue analysis in principle, however with a slightly larger value of the critical Weißenberg number. The numerical examples given here are only relevant test cases. Application of the obtained results to real flows on classical benchmarks, such as the flow in a 4:1 contraction, will be presented in a next paper. CONCLUSIONS REFERENCES In this paper, we studied the stability of the Oldroyd system by both a formal linear spectral analysis and numerical simulations. Our findings are summarized as follows. The pure stress equation without the β-term as well as the linearized Oldroyd system are both unconditionally stable. The possible source for instability is the β-term, because it is a zero-order term with indefinite sign. By linear stability analysis, we found conditions to characterize the region of stability, both for the pure stress equation and for the full Oldroyd system. These stability limits were fairly well confirmed by the simulations. Considering the pure stress equation without the stress transport term and for a given stationary velocity field, in the case of a = 0, a stability limit can exist. More precisely, if their exists a point or a region in the computational domain where the velocity field fulfills det(∇u) < 0, then for√a = ±1 the stress equation is stable until Wecr = (2 − det(∇u))−1 . In the case of a = 0, no stability limit was found. Considering the stress constitutive equation with given (stationary) velocity field, the stress is transported along the streamlines. Along streamlines, which leave the computational domain, no instability arises. However, if there exists a stagnation point of the flow in the region where det(∇u) < 0 (in the case a = ±1), or if the streamlines are closed curves that intersect such a region, then instability of the stress components along the streamlines arises for Weißenberg numbers above the critical value. For a = 0 the pure stress equation is not affected, but in the full Oldroyd system stability limits exist. For the Oldroyd system without the stress transport term, in the case a = 0, one detects the same stability limit Wecr . However, if the velocity field fulfills det(∇u) ≥ 0 everywhere, the system remains stable 2009 Curtin University of Technology and John Wiley & Sons, Ltd. [1] J.G. Oldroyd. Proc. R. Soc. London, 1950; A 200, 523–541. [2] F.P.T. Baaijens. J. Non-Newtonian Fluid Mech., 1998; 79, 361–385. [3] C. Bodart, M.J. Crochet. Theoret. Comput. Fluid Dynamics, 1993; 5, 57–75. [4] R. Guénette, M. Fortin. J. Non-Newtonian Fluid Mech., 1995; 60, 27–52. [5] D. Sandri. Comput. Meth. Appl. Mech. Eng., 2002; 191, 5045–5065. [6] P. Saramito. Math. Modell. Num. Anal., 1994; 28, 1–34. [7] P. Wapperon, M.F. Webster. Comput. Meth. Appl. Mech. Eng., 1999; 180(3-4), 281–304. [8] M. Johnson, D. Segalman. J. Non-Newtonian Fluid Mech., 1977; 2, 255–270. [9] M.A. Alves, P.J. Oliveira, F.T. Pinho. J. Non-Newtonian Fluid Mech., 2003; 110, 45–75. [10] R. Keunings. J. Non-Newtonian Fluid Mech., 1986; 20, 209–226. [11] B. Caswell. J. Non-Newtonian Fluid Mech., 1996; 62, 99–110. [12] M. Renardy. Mathematical Analysis of Viscoelastic Flows, CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM 2000. [13] E. Fernández-Cara, F. Guillén, R.R. Ortega. Mathematical modeling and analysis of viscoelastic fluids of the Oldroyd kind, Handbook of Numerical Analysis, Vol VIII, NorthHolland, Amsterdam 2002; pp.543–661. [14] M. Picasso, J. Rappaz. ESIAM: Math. Modell. Num. Anal., 2001; 35(5), 879–897. [15] A. Bonito, P. Clément, M. Picasso. Numer. Math., 2007; 107, 213–255. [16] P.L. Lions, N. Masmoudi. China Ann. Math., 2000; 21B(2), 131–146. [17] H. Giesekus. Phänomenologische Rheologie, Springer, Berlin, Heidelberg 1994. [18] N. Phan-Thien, R.I. Tanner. J. Non-Newtonian Fluid Mech., 1977; 2, 353–365. [19] N. Scurtu. Stability Analysis and Numerical Simulation of NonNewtonian Fluids of Oldroyd Kind. PhD Thesis. FriedrichAlexander-Universität Erlangen-Nürnberg 2005. [20] B. Jourdain, C. Le Bris, T. Leliévre, F. Otto. J. Arch. Rational Mech. Anal., 2006; 181, 97–148. [21] M.O. Bristeau, R. Glowinski, J. Periaux. Comput. Phys. Rep., 1987; 6, 73–187. [22] M. Fortin, A. Fortin. J. Non-Newtonian Fluid Mech., 1989; 32, 295–310. [23] A. Schmidt, K.G. Siebert.Design of adaptive finite element software: The finite element toolbox ALBERTA. Springer LNCSE Series 42 2005. Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering [24] V.J. Ervin, L.N. Ntasin. SIAM Numer. Meth. Diff. Equations, 2004; 21, 297–322. [25] E. Bänsch. Acta Math. Univ. Comenianae, 1998; LXVII, 101–114. 2009 Curtin University of Technology and John Wiley & Sons, Ltd. INFLUENCE OF WEIßENBERG NUMBER [26] R.B. Bird, R.C. Armstrong, O. Hassager. Dynamics of Polymeric Liquids, John Wiley and Sons, New York 1977. Asia-Pac. J. Chem. Eng. 2010; 5: 657–673 DOI: 10.1002/apj 673

1/--страниц