- Email: [email protected]
Stability analysis using compressible viscoelastic formulation
J. Non-Newtonian Fluid Mech. 116 (2004) 471–477
Short communication
Stability analysis using compressible viscoelastic formulation R. Sureshkumar Department of Chemical Engineering, Materials Research Laboratory, Washington University, Saint Louis, MO 63130, USA Received 5 August 2003; received in revised form 9 November 2003
Abstract Compressible viscoelastic formulation (CVF) is used to perform stability analysis of planar and curvilinear flows. Incompressible eigensolutions are recovered from CVF in the limit of small compressibility. CVF allows for the derivation of a simple eigenvalue problem for stability analysis as compared to a generalized eigenvalue problem (GEVP) that results from the incompressible flow equations. Hence computational cost can be significantly reduced. A minimal CVF useful for multi-dimensional flows is presented. © 2003 Elsevier B.V. All rights reserved. Keywords: Viscoelastic; Compressible; Stability analysis
Stability analysis of viscoelastic flows is computationally demanding. The recent increase in computational power has helped reduce the CPU-time requirements and allows for the direct solution of larger eigenvalue problems. However, this is not adequate to overcome the computational requirements associated with the stability analysis of multi-dimensional viscoelastic flows. The difficulties range from the numerical resolution of the continuous part of the eigenspectrum to the isolation of spurious eigenvalues from the physically-relevant ones [1]. Knowledge of the composition of the eigenspectrum and the structure of the eigenfunctions in benchmark geometries for multi-dimensional flows could also prove to be useful in developing reliable algorithms for time dependent simulations. Linear stability analysis of incompressible viscoelastic flows based on normal mode analysis requires the solution of a generalized eigenvalue problem (GEVP). The direct solution of the GEVP is expensive even for relatively coarse meshes. Recently, Arora and Sureshkumar [2] proposed an algorithm based on submatrix-based transformation of the linearized equations (SubTLE) to help reduce the CPU-time associated with the solution of the GEVP. SubTLE has been implemented to enable global linear stability analysis of multidimensional viscoelastic flows [3]. However, this algorithm is applicable only in the creeping flow limit and requires the LU decomposition of a submatrix to transform the GEVP to a simple EVP of half the original size. Hence, it is desirable to develop algorithms that will allow for the derivation E-mail address: [email protected] (R. Sureshkumar). 0377-0257/$ – see front matter © 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2003.11.003
472
R. Sureshkumar / J. Non-Newtonian Fluid Mech. 116 (2004) 471–477
of a simple EVP directly from the governing equations, thereby obviating the need to perform expensive LU decompositions. Such an algorithm is developed and validated in this communication. We start with compressible viscoelastic formulation (CVF) described by the continuity equation, momentum conservation law and a differential constitutive model [4,5]. Without loss of generality, the upper convected Maxwell (UCM) constitutive model is considered in this work. We denote the length and velocity scales by L and U, respectively. The governing equations are made dimensionless by using L/U as the time scale, ps ≡ ηU/L as scale for pressure and stress where η denotes the fluid viscosity and 1/ps as the scale for compressibility. The continuity equation is given by: Dp + K∇ · u = 0, (1) Dt where p and u denote the pressure and velocity, respectively, K is the inverse of the isothermal compressibility and Dp ∂p ≡ + u · ∇p. Dt ∂t Eq. (1) is based on the p–density (ρ) relationship, p − p0 ρ = ρo 1 + , K
(1a)
(2)
valid in the limit |p − p0 |/K 1 [4]. The specific form of the p–ρ relationship is not of consequence in this analysis since compressibility is introduced as an artifact to retain the time derivative of pressure in Eq. (1). As shown below, a sufficiently large value of K may be used to recover the incompressible eigensolutions. The momentum equation is given by: Du 1 = [−∇p + ∇ · T + f ], (3) Dt Re where Re ≡ ρ0 UL/η is the Reynolds number, T is the viscoelastic stress tensor and f represents the additional terms arising from compressibility effects (see discussion following Eq. (4)). The evolution equation for T is given by: DT 1 = T · ∇u + (∇u)t · T − T ∇ · u − (4) [T − (∇u + (∇u)t )], Dt We where We ≡ λU/L is the Weissenberg number with λ denoting the fluid relaxation time. The term T ∇·u is due to Edwards and Beris [5] who accounted for the variability in the number density of polymer molecules caused by compressibility effects. For differential constitutive models of the UCM type, setting f = −∇(trace (T ))/3 in Eq. (3) allows for the Newtonian limit of f = ∇(∇·u)/3 to be recovered in the limit as We → 0 [4,5]. However, this choice of f will influence the steady state solutions appreciably except for homogeneous flows. Hence, we choose f = −(2/3)∇(∇·u). This choice, together with Eq. (4), allows for the recovery of the Newtonian limit as We → 0. In the limit as 1/K → 0, Eqs. (1), (3) and (4) reduce to the governing equations for the incompressible UCM fluid. Eqs. (1), (3) and (4) are supplemented by the no slip condition for velocity at solid boundaries. We first consider three shear flows: (i) plane Couette flow in the creeping flow regime for which analytical results exist for the eigenspectrum [6,7], (ii) plane Poiseuille flow and (iii) Taylor–Couette
R. Sureshkumar / J. Non-Newtonian Fluid Mech. 116 (2004) 471–477
473
flow, i.e. the flow generated in the gap of rotating coaxial cylinders. For cases (ii) and (iii) results are validated against those reported in [8–10], respectively. For these unidirectional shear flows, the steady states of the incompressible flow equations are also solution to Eqs. (1), (3) and (4). The velocity and length scales used for non-dimensionalization of the governing equations are as follows: the velocity of the top plate and the separation between plates for the plane Couette flow where only the top plate is moving; the centerline velocity and channel half height for the plane Poiseuille flow; and the velocity of the inner cylinder and the gap width for the Taylor–Couette flow. For the plane Couette and Poiseuille flows, the steady state is perturbed by two-dimensional normal mode perturbations of the type: q(x, y, t) = Real(φ(y)exp(iαx + σt)),
(5)
where q denotes any one of the flow variables, x and y denote the (steady) flow and wall-normal directions, respectively, and α and σ denote the wavenumber and the eigenvalue, respectively. For the Taylor–Couette flow, the most dangerous disturbance could be non-axisymmetric. Hence, the analysis is performed against three-dimensional perturbations of the type: q(r, θ, z, t) = Real(ϕ(r)exp(iαz + iξθ + σt)),
(5a)
where r, θ, and z denote the cylindrical coordinates with z-axis coinciding with the axis of the cylinders and α and ξ are the wavenumbers in the z and θ directions, respectively. Substitution of the normal mode perturbations into Eq. (1), (3) and (4) and linearization about the steady state result in a differential eigenvalue problem (DEVP) of the type Lv = σv where v denotes the vector of eigenfunctions of the flow variables p, u and T . The analogous DEVP for incompressible flow is of the form Lv = σMv where M is a singular operator due to the incompressibility constraint, i.e. the continuity equation does not contain explicit time derivate of pressure. We use a Chebyshev collocation algorithm to numerically discretize the eigenfunctions [2,8,10]. This results in a simple algebraic eigenvalue problem of dimension 6N in the case of plane Couette and Poiseuille flows and 10N in the case of the Taylor–Couette flow where N is the order of the Chebyshev expansion used to represent the eigenfunctions. The simple eigenvalue problem resulting from CVF is solved using a lapack subroutine zgeev.f. For comparison purposes and to estimate the CPU-time savings, the generalized eigenvalue problem resulting from the incompressible flow equations is also solved using the lapack subroutine zgegv.f. Several calculations were performed to validate the use of CVF. Typical results are reported in Figs. 1–3 for the plane Couette flow (Re = 0.1, We = 10), the plane Poiseuille flow (Re = 2310, We = 2.31) and the Taylor–Couette flow (Re = 0.001, We = 38.7) respectively. In all three cases, K = 106 . For comparison purposes, the results obtained from the incompressible formulation are also presented. For the plane Couette flow, the two Gorodtsov–Leonov eigenvalues [6,7] and the numerical approximation of the continuous eigenspectrum located at the −1/We line are correctly captured by CVF. Similar agreement with the incompressible formulation can be seen from Figs. 2 and 3 for the plane Poiseuille (see discussion associated with Fig. 1 in ref. [8] for details of the spectrum) and the Taylor–Couette flows (see discussion associated with Table 2 in ref. [10]) respectively. The leading eigenvalue(s) and the remainder of the spectrum are captured accurately by CVF. While using CVF in the purely elastic limit, a finite value of Re is required. For the plane Couette flow, the perturbation pressure p is given by: (σ + iαU)p = K∇ · u,
(6)
474
R. Sureshkumar / J. Non-Newtonian Fluid Mech. 116 (2004) 471–477
Fig. 1. Eigenspectrum of the plane Couette flow for Re = 0.01, We = 10, α = 1. CVF ( ), GEVP (䊐).
where u denotes the perturbation velocity and U = y is the steady state velocity profile. Hence the influence of compressibility in the right hand side of the momentum equation (i.e. from f = −(2/3)∇(∇·u) in Eq. (3)) is inversely proportional to ReK. This implies that K 1/Re to recover the incompressible solutions. This is consistent with numerical observations. For instance, when Re was decreased to 10−3 from 10−1 for the case reported in Fig. 1, K had to be increased from 106 to 108 to get comparable agreement with the incompressible flow solution. The compressibility effect appears in the constitutive equation through the term T ∇·u that contributes as T ss ∇·u in the eigenvalue problem where T ss denotes
Fig. 2. Eigenspectrum of the plane Poiseuille flow for Re = 2310, We = 2.31, α = 1.3. CVF ( ), GEVP (䊐).
R. Sureshkumar / J. Non-Newtonian Fluid Mech. 116 (2004) 471–477
475
Fig. 3. Eigenspectrum of the Taylor–Couette flow for Re = 0.001, We = 38.7, α = 5, ξ = 2, ratio of cylinder radii = 0.95, ratio of the rotational speeds of the inner to the outer cylinder = 0.5. CVF ( ), GEVP (䊐).
the steady state stress tensor. Since the (dimensionless) normal stress increases proportional to We in the Couette flow, Eq. (6) implies that We/K 1 for compressibility effects to be negligible. This is realizable if a sufficiently large value of K is used. Since in the limit as K → ∞ the influence of compressibility is negligible, it is tempting to set the compressibility corrections in the momentum and constitutive equations to zero. In order to assess the influence of such an approximation on the eigenvalues, the term f in Eq. (3) and the term T ∇·u in Eq. (4) as well as the convective term u·∇p in Eq. (1), were set to zero and the eigenspectrum was recomputed for each of the three cases reported in Figs. 1–3. As expected, the effect on the eigenvalues was imperceptible. This suggests the use of a minimal CVF in which the momentum and constitutive equations for the incompressible fluid are supplemented by the modified continuity equation ∂p/∂t + K∇ · u = 0. Minimal CVF could be useful for the stability analysis of multi-dimensional flows since the introduction of the time derivative of pressure does not alter the (incompressible) steady state solution. Its implementation in the solution of the EVP is easy. For instance, minimal CVF was applied to the computation of the eigenspectrum of the flow between eccentric cylinders of a UCM fluid [11,12]. The steady state was computed by using a Chebyshev–Fourier collocation method. Fig. 4 shows the leading eigenvalues of the spectrum obtained for the following condition where the inner cylinder rotates and the outer one is stationary: ratio of average gap width d to inner cylinder radius R1 = 0.0526, eccentricity ε = 0.2d and We ≡ λV1 /d = 11.41, where V1 is the velocity of the inner cylinder. This value of We is very near the critical condition for αd = 7.8 [11]. This represents approximately 50% decrease in critical We with respect to the Taylor–Couette flow realized for ε = 0. In CVF finite values of Re (=10−3 ) and K = 109 were used. The results for the incompressible flow were computed for Re ≡ 0 and by introducing a small value of solvent viscosity = 1% of the total viscosity using the SubTLE algorithm [2]. A mesh of 17 Chebyshev and 15 Fourier nodes is used. As seen from Fig. 4, the leading eigenvalues are captured accurately by CVF. The difference between the most dangerous eigenvalue obtained from the two methods is
476
R. Sureshkumar / J. Non-Newtonian Fluid Mech. 116 (2004) 471–477
Fig. 4. Eigenspectrum of the flow between eccentric cylinders for Re = 0.001, We = 11.4. CVF ( ), GEVP (䊐).
approximately one part in ten thousand. In addition to the compressibility effect, this difference includes those arising from the use of finite Re (in CVF) and finite solvent viscosity required for SubTLE [2]. CPU-time savings resulting from the implementation of CVF are significant. For the plane Couette flow the solution of a GEVP of dimension 774 took approximately 110 s on a DEC alpha station whereas CVF required only 30 s on the same machine. Such CPU-time savings are useful in reducing the turnover time for stability analysis that requires numerous calculations of the spectrum to determine the neutral stability conditions as well as parametric studies. For the eccentric cylinder flow, SubTLE and CVF algorithms required nearly the same amount of CPU-time. However, we note that SubTLE works only in the Re ≡ 0 limit. Moreover, CVF is more amenable for multi-dimensional flows where direct solution of the EVP may be infeasible and iterative procedures (e.g. Arnoldi method [8]) will have to be used to compute the leading eigenvalue. The bottleneck in employing iterative solvers is the step that involves the conversion of the GEVP into a simple eigenvalue problem on which the iterative procedure is used. This step requires LU decomposition of the dense coefficient matrix which is prohibitively time and memory intensive. Once fully tested and benchmarked, CVF could be employed for the stability analysis of multi-dimensional flows since it results in a simple EVP. The use of small values of compressibility may be justified based on physical grounds as well.
Acknowledgements Support from NSF through grant CTS 9874813 and the Donors of The Petroleum Research Fund administered by the ACS through grant 36098-AC9 are gratefully acknowledged. The author is thankful to Mr. Kartik Arora for useful discussions and providing the data for Fig. 4.
R. Sureshkumar / J. Non-Newtonian Fluid Mech. 116 (2004) 471–477
477
References [1] B. Sadanandan, R. Sureshkumar, Numerical eigenspectrum of non-viscometric viscoelastic flows: results for the periodic channel flow, J. Non-Newtonian Fluid Mech. 108 (2002) 143. [2] K. Arora, R. Sureshkumar, Efficient computation of the eigenspectrum of viscoelastic flows using submatrix-based transformation of the linearized equations (SubTLE), J. Non-Newtonian Fluid Mech. 104 (2002) 75. [3] B. Sadanandan, R. Sureshkumar, Global linear stability analysis of non-separated viscoelastic flow through a periodically constricted channel, submitted for publication. [4] F.R. Phelan Jr., M.F. Malone, H.H. Winter, A purely hyperbolic model for unsteady viscoelastic flow, J. Non-Newtonian Fluid Mech. 32 (1989) 197. [5] B.J. Edwards, A.N. Beris, Remarks concerning compressible viscoelastic fluid models, J. Non-Newtonian Fluid Mech. 36 (1990) 411. [6] V.A. Gorodtsov, A.I. Leonov, On a linear instability of a plane parallel Couette flow of viscoelastic fluid, J. Appl. Math. Mech. 31 (1967) 310. [7] M. Renardy, Y. Renardy, Linear stability of plane Couette flow of an upper convected Maxwell fluid, J. Non-Newtonian Fluid Mech. 22 (1986) 23. [8] R. Sureshkumar, A.N. Beris, Linear stability analysis of the viscoelastic Poiseuille flow using an Arnoldi orthogonalization algorithm, J. Non-Newtonian Fluid Mech. 56 (1995) 151. [9] H.J. Wilson, M. Renardy, Y. Renardy, Structure of the spectrum in zero Reynolds number shear flow of the UCM and Oldroyd-B liquids, J. Non-Newtonian Fluid Mech. 80 (1999) 251. [10] R. Sureshkumar, A.N. Beris, M. Avgousti, Non-axisymmetric subcritical bifurcations in viscoelastic Taylor–Couette flow, Proc. R. Soc. London Ser. A 447 (1994) 135. [11] A. Chawda, M. Avgousti, Stability of viscoelastic flow between eccentric rotating cylinders, J. Non-Newtonian Fluid Mech. 63 (1996) 97. [12] R. Sureshkumar, M. Avgousti, Influence of eccentricity on the stability of purely elastic Dean flow, J. Non-Newtonian Fluid Mech. 93 (2000) 61.
Short communication
Stability analysis using compressible viscoelastic formulation R. Sureshkumar Department of Chemical Engineering, Materials Research Laboratory, Washington University, Saint Louis, MO 63130, USA Received 5 August 2003; received in revised form 9 November 2003
Abstract Compressible viscoelastic formulation (CVF) is used to perform stability analysis of planar and curvilinear flows. Incompressible eigensolutions are recovered from CVF in the limit of small compressibility. CVF allows for the derivation of a simple eigenvalue problem for stability analysis as compared to a generalized eigenvalue problem (GEVP) that results from the incompressible flow equations. Hence computational cost can be significantly reduced. A minimal CVF useful for multi-dimensional flows is presented. © 2003 Elsevier B.V. All rights reserved. Keywords: Viscoelastic; Compressible; Stability analysis
Stability analysis of viscoelastic flows is computationally demanding. The recent increase in computational power has helped reduce the CPU-time requirements and allows for the direct solution of larger eigenvalue problems. However, this is not adequate to overcome the computational requirements associated with the stability analysis of multi-dimensional viscoelastic flows. The difficulties range from the numerical resolution of the continuous part of the eigenspectrum to the isolation of spurious eigenvalues from the physically-relevant ones [1]. Knowledge of the composition of the eigenspectrum and the structure of the eigenfunctions in benchmark geometries for multi-dimensional flows could also prove to be useful in developing reliable algorithms for time dependent simulations. Linear stability analysis of incompressible viscoelastic flows based on normal mode analysis requires the solution of a generalized eigenvalue problem (GEVP). The direct solution of the GEVP is expensive even for relatively coarse meshes. Recently, Arora and Sureshkumar [2] proposed an algorithm based on submatrix-based transformation of the linearized equations (SubTLE) to help reduce the CPU-time associated with the solution of the GEVP. SubTLE has been implemented to enable global linear stability analysis of multidimensional viscoelastic flows [3]. However, this algorithm is applicable only in the creeping flow limit and requires the LU decomposition of a submatrix to transform the GEVP to a simple EVP of half the original size. Hence, it is desirable to develop algorithms that will allow for the derivation E-mail address: [email protected] (R. Sureshkumar). 0377-0257/$ – see front matter © 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2003.11.003
472
R. Sureshkumar / J. Non-Newtonian Fluid Mech. 116 (2004) 471–477
of a simple EVP directly from the governing equations, thereby obviating the need to perform expensive LU decompositions. Such an algorithm is developed and validated in this communication. We start with compressible viscoelastic formulation (CVF) described by the continuity equation, momentum conservation law and a differential constitutive model [4,5]. Without loss of generality, the upper convected Maxwell (UCM) constitutive model is considered in this work. We denote the length and velocity scales by L and U, respectively. The governing equations are made dimensionless by using L/U as the time scale, ps ≡ ηU/L as scale for pressure and stress where η denotes the fluid viscosity and 1/ps as the scale for compressibility. The continuity equation is given by: Dp + K∇ · u = 0, (1) Dt where p and u denote the pressure and velocity, respectively, K is the inverse of the isothermal compressibility and Dp ∂p ≡ + u · ∇p. Dt ∂t Eq. (1) is based on the p–density (ρ) relationship, p − p0 ρ = ρo 1 + , K
(1a)
(2)
valid in the limit |p − p0 |/K 1 [4]. The specific form of the p–ρ relationship is not of consequence in this analysis since compressibility is introduced as an artifact to retain the time derivative of pressure in Eq. (1). As shown below, a sufficiently large value of K may be used to recover the incompressible eigensolutions. The momentum equation is given by: Du 1 = [−∇p + ∇ · T + f ], (3) Dt Re where Re ≡ ρ0 UL/η is the Reynolds number, T is the viscoelastic stress tensor and f represents the additional terms arising from compressibility effects (see discussion following Eq. (4)). The evolution equation for T is given by: DT 1 = T · ∇u + (∇u)t · T − T ∇ · u − (4) [T − (∇u + (∇u)t )], Dt We where We ≡ λU/L is the Weissenberg number with λ denoting the fluid relaxation time. The term T ∇·u is due to Edwards and Beris [5] who accounted for the variability in the number density of polymer molecules caused by compressibility effects. For differential constitutive models of the UCM type, setting f = −∇(trace (T ))/3 in Eq. (3) allows for the Newtonian limit of f = ∇(∇·u)/3 to be recovered in the limit as We → 0 [4,5]. However, this choice of f will influence the steady state solutions appreciably except for homogeneous flows. Hence, we choose f = −(2/3)∇(∇·u). This choice, together with Eq. (4), allows for the recovery of the Newtonian limit as We → 0. In the limit as 1/K → 0, Eqs. (1), (3) and (4) reduce to the governing equations for the incompressible UCM fluid. Eqs. (1), (3) and (4) are supplemented by the no slip condition for velocity at solid boundaries. We first consider three shear flows: (i) plane Couette flow in the creeping flow regime for which analytical results exist for the eigenspectrum [6,7], (ii) plane Poiseuille flow and (iii) Taylor–Couette
R. Sureshkumar / J. Non-Newtonian Fluid Mech. 116 (2004) 471–477
473
flow, i.e. the flow generated in the gap of rotating coaxial cylinders. For cases (ii) and (iii) results are validated against those reported in [8–10], respectively. For these unidirectional shear flows, the steady states of the incompressible flow equations are also solution to Eqs. (1), (3) and (4). The velocity and length scales used for non-dimensionalization of the governing equations are as follows: the velocity of the top plate and the separation between plates for the plane Couette flow where only the top plate is moving; the centerline velocity and channel half height for the plane Poiseuille flow; and the velocity of the inner cylinder and the gap width for the Taylor–Couette flow. For the plane Couette and Poiseuille flows, the steady state is perturbed by two-dimensional normal mode perturbations of the type: q(x, y, t) = Real(φ(y)exp(iαx + σt)),
(5)
where q denotes any one of the flow variables, x and y denote the (steady) flow and wall-normal directions, respectively, and α and σ denote the wavenumber and the eigenvalue, respectively. For the Taylor–Couette flow, the most dangerous disturbance could be non-axisymmetric. Hence, the analysis is performed against three-dimensional perturbations of the type: q(r, θ, z, t) = Real(ϕ(r)exp(iαz + iξθ + σt)),
(5a)
where r, θ, and z denote the cylindrical coordinates with z-axis coinciding with the axis of the cylinders and α and ξ are the wavenumbers in the z and θ directions, respectively. Substitution of the normal mode perturbations into Eq. (1), (3) and (4) and linearization about the steady state result in a differential eigenvalue problem (DEVP) of the type Lv = σv where v denotes the vector of eigenfunctions of the flow variables p, u and T . The analogous DEVP for incompressible flow is of the form Lv = σMv where M is a singular operator due to the incompressibility constraint, i.e. the continuity equation does not contain explicit time derivate of pressure. We use a Chebyshev collocation algorithm to numerically discretize the eigenfunctions [2,8,10]. This results in a simple algebraic eigenvalue problem of dimension 6N in the case of plane Couette and Poiseuille flows and 10N in the case of the Taylor–Couette flow where N is the order of the Chebyshev expansion used to represent the eigenfunctions. The simple eigenvalue problem resulting from CVF is solved using a lapack subroutine zgeev.f. For comparison purposes and to estimate the CPU-time savings, the generalized eigenvalue problem resulting from the incompressible flow equations is also solved using the lapack subroutine zgegv.f. Several calculations were performed to validate the use of CVF. Typical results are reported in Figs. 1–3 for the plane Couette flow (Re = 0.1, We = 10), the plane Poiseuille flow (Re = 2310, We = 2.31) and the Taylor–Couette flow (Re = 0.001, We = 38.7) respectively. In all three cases, K = 106 . For comparison purposes, the results obtained from the incompressible formulation are also presented. For the plane Couette flow, the two Gorodtsov–Leonov eigenvalues [6,7] and the numerical approximation of the continuous eigenspectrum located at the −1/We line are correctly captured by CVF. Similar agreement with the incompressible formulation can be seen from Figs. 2 and 3 for the plane Poiseuille (see discussion associated with Fig. 1 in ref. [8] for details of the spectrum) and the Taylor–Couette flows (see discussion associated with Table 2 in ref. [10]) respectively. The leading eigenvalue(s) and the remainder of the spectrum are captured accurately by CVF. While using CVF in the purely elastic limit, a finite value of Re is required. For the plane Couette flow, the perturbation pressure p is given by: (σ + iαU)p = K∇ · u,
(6)
474
R. Sureshkumar / J. Non-Newtonian Fluid Mech. 116 (2004) 471–477
Fig. 1. Eigenspectrum of the plane Couette flow for Re = 0.01, We = 10, α = 1. CVF ( ), GEVP (䊐).
where u denotes the perturbation velocity and U = y is the steady state velocity profile. Hence the influence of compressibility in the right hand side of the momentum equation (i.e. from f = −(2/3)∇(∇·u) in Eq. (3)) is inversely proportional to ReK. This implies that K 1/Re to recover the incompressible solutions. This is consistent with numerical observations. For instance, when Re was decreased to 10−3 from 10−1 for the case reported in Fig. 1, K had to be increased from 106 to 108 to get comparable agreement with the incompressible flow solution. The compressibility effect appears in the constitutive equation through the term T ∇·u that contributes as T ss ∇·u in the eigenvalue problem where T ss denotes
Fig. 2. Eigenspectrum of the plane Poiseuille flow for Re = 2310, We = 2.31, α = 1.3. CVF ( ), GEVP (䊐).
R. Sureshkumar / J. Non-Newtonian Fluid Mech. 116 (2004) 471–477
475
Fig. 3. Eigenspectrum of the Taylor–Couette flow for Re = 0.001, We = 38.7, α = 5, ξ = 2, ratio of cylinder radii = 0.95, ratio of the rotational speeds of the inner to the outer cylinder = 0.5. CVF ( ), GEVP (䊐).
the steady state stress tensor. Since the (dimensionless) normal stress increases proportional to We in the Couette flow, Eq. (6) implies that We/K 1 for compressibility effects to be negligible. This is realizable if a sufficiently large value of K is used. Since in the limit as K → ∞ the influence of compressibility is negligible, it is tempting to set the compressibility corrections in the momentum and constitutive equations to zero. In order to assess the influence of such an approximation on the eigenvalues, the term f in Eq. (3) and the term T ∇·u in Eq. (4) as well as the convective term u·∇p in Eq. (1), were set to zero and the eigenspectrum was recomputed for each of the three cases reported in Figs. 1–3. As expected, the effect on the eigenvalues was imperceptible. This suggests the use of a minimal CVF in which the momentum and constitutive equations for the incompressible fluid are supplemented by the modified continuity equation ∂p/∂t + K∇ · u = 0. Minimal CVF could be useful for the stability analysis of multi-dimensional flows since the introduction of the time derivative of pressure does not alter the (incompressible) steady state solution. Its implementation in the solution of the EVP is easy. For instance, minimal CVF was applied to the computation of the eigenspectrum of the flow between eccentric cylinders of a UCM fluid [11,12]. The steady state was computed by using a Chebyshev–Fourier collocation method. Fig. 4 shows the leading eigenvalues of the spectrum obtained for the following condition where the inner cylinder rotates and the outer one is stationary: ratio of average gap width d to inner cylinder radius R1 = 0.0526, eccentricity ε = 0.2d and We ≡ λV1 /d = 11.41, where V1 is the velocity of the inner cylinder. This value of We is very near the critical condition for αd = 7.8 [11]. This represents approximately 50% decrease in critical We with respect to the Taylor–Couette flow realized for ε = 0. In CVF finite values of Re (=10−3 ) and K = 109 were used. The results for the incompressible flow were computed for Re ≡ 0 and by introducing a small value of solvent viscosity = 1% of the total viscosity using the SubTLE algorithm [2]. A mesh of 17 Chebyshev and 15 Fourier nodes is used. As seen from Fig. 4, the leading eigenvalues are captured accurately by CVF. The difference between the most dangerous eigenvalue obtained from the two methods is
476
R. Sureshkumar / J. Non-Newtonian Fluid Mech. 116 (2004) 471–477
Fig. 4. Eigenspectrum of the flow between eccentric cylinders for Re = 0.001, We = 11.4. CVF ( ), GEVP (䊐).
approximately one part in ten thousand. In addition to the compressibility effect, this difference includes those arising from the use of finite Re (in CVF) and finite solvent viscosity required for SubTLE [2]. CPU-time savings resulting from the implementation of CVF are significant. For the plane Couette flow the solution of a GEVP of dimension 774 took approximately 110 s on a DEC alpha station whereas CVF required only 30 s on the same machine. Such CPU-time savings are useful in reducing the turnover time for stability analysis that requires numerous calculations of the spectrum to determine the neutral stability conditions as well as parametric studies. For the eccentric cylinder flow, SubTLE and CVF algorithms required nearly the same amount of CPU-time. However, we note that SubTLE works only in the Re ≡ 0 limit. Moreover, CVF is more amenable for multi-dimensional flows where direct solution of the EVP may be infeasible and iterative procedures (e.g. Arnoldi method [8]) will have to be used to compute the leading eigenvalue. The bottleneck in employing iterative solvers is the step that involves the conversion of the GEVP into a simple eigenvalue problem on which the iterative procedure is used. This step requires LU decomposition of the dense coefficient matrix which is prohibitively time and memory intensive. Once fully tested and benchmarked, CVF could be employed for the stability analysis of multi-dimensional flows since it results in a simple EVP. The use of small values of compressibility may be justified based on physical grounds as well.
Acknowledgements Support from NSF through grant CTS 9874813 and the Donors of The Petroleum Research Fund administered by the ACS through grant 36098-AC9 are gratefully acknowledged. The author is thankful to Mr. Kartik Arora for useful discussions and providing the data for Fig. 4.
R. Sureshkumar / J. Non-Newtonian Fluid Mech. 116 (2004) 471–477
477
References [1] B. Sadanandan, R. Sureshkumar, Numerical eigenspectrum of non-viscometric viscoelastic flows: results for the periodic channel flow, J. Non-Newtonian Fluid Mech. 108 (2002) 143. [2] K. Arora, R. Sureshkumar, Efficient computation of the eigenspectrum of viscoelastic flows using submatrix-based transformation of the linearized equations (SubTLE), J. Non-Newtonian Fluid Mech. 104 (2002) 75. [3] B. Sadanandan, R. Sureshkumar, Global linear stability analysis of non-separated viscoelastic flow through a periodically constricted channel, submitted for publication. [4] F.R. Phelan Jr., M.F. Malone, H.H. Winter, A purely hyperbolic model for unsteady viscoelastic flow, J. Non-Newtonian Fluid Mech. 32 (1989) 197. [5] B.J. Edwards, A.N. Beris, Remarks concerning compressible viscoelastic fluid models, J. Non-Newtonian Fluid Mech. 36 (1990) 411. [6] V.A. Gorodtsov, A.I. Leonov, On a linear instability of a plane parallel Couette flow of viscoelastic fluid, J. Appl. Math. Mech. 31 (1967) 310. [7] M. Renardy, Y. Renardy, Linear stability of plane Couette flow of an upper convected Maxwell fluid, J. Non-Newtonian Fluid Mech. 22 (1986) 23. [8] R. Sureshkumar, A.N. Beris, Linear stability analysis of the viscoelastic Poiseuille flow using an Arnoldi orthogonalization algorithm, J. Non-Newtonian Fluid Mech. 56 (1995) 151. [9] H.J. Wilson, M. Renardy, Y. Renardy, Structure of the spectrum in zero Reynolds number shear flow of the UCM and Oldroyd-B liquids, J. Non-Newtonian Fluid Mech. 80 (1999) 251. [10] R. Sureshkumar, A.N. Beris, M. Avgousti, Non-axisymmetric subcritical bifurcations in viscoelastic Taylor–Couette flow, Proc. R. Soc. London Ser. A 447 (1994) 135. [11] A. Chawda, M. Avgousti, Stability of viscoelastic flow between eccentric rotating cylinders, J. Non-Newtonian Fluid Mech. 63 (1996) 97. [12] R. Sureshkumar, M. Avgousti, Influence of eccentricity on the stability of purely elastic Dean flow, J. Non-Newtonian Fluid Mech. 93 (2000) 61.