c
° TU Delft, The Netherlands, 2006
EXPLICIT AND SEMI-IMPLICIT CHARACTERISTIC
BASED SPLIT (CBS) SCHEMES FOR VISCOELASTIC
FLOW CALCULATIONS
C-B. Liu∗, P. Nithiarasu†
University of Wales Swansea, School of Engineering Swansea SA2 8PP, UK
∗e-mail: C.B.Liu@swansea.ac.uk † e-mail: P.Nithiarasu@Swansea.ac.uk
web page: http://www.engineering.swan.ac.uk/
Key words: CBS scheme, artificial compressibility(AC), semi-implicit, Oldroyd-B, Maxwell, PTT, viscoelastic flow
Abstract. A fully explicit characteristic based split (CBS) based artificial
compressibil-ity (AC) method and a semi-implicit CBS scheme, for viscoelastic flow over a circular cylinder placed in a rectangular channel are presented in this paper. In the explicit form the pressure equation is solved explicitly by introducing an artificial compressibility para-meter, while a pressure Poisson equation is implicitly solved in the semi-implicit form. Three different constitutive models have been employed to investigate flow past a circu-lar cylinder. They are the Oldroyd-B, upper convected Maxwell (UCM) and simplified Phan-Thien/Tanner (PTT) models. The main difference from the previously published results is that in the present study the artificial dissipation is not switched on until it is absolutely essential. No loss of convergence to steady state was observed in any of the results presented in this paper. Comparison of the present numerical solution with other available numerical data shows that the CBS algorithm is in excellent agreement at lower Deborah numbers. However, at higher Deborah numbers, the present results differ sub-stantially from other numerical results. The two presented schemes give almost identical drag force values and the solutions presented demonstrate that both schemes are suitable for viscoelastic flow calculations.
1 INTRODUCTION
Numerical solution of such equations are often subjected to strong extra stress oscillations, especially at higher Deborah numbers. It is also obvious from the available literature that obtaining convergence to steady state is extremely difficult at higher Deborah numbers.
In order to avoid numerical instability due to high elasticity in viscoelastic flows, most of papers suggest splitting the extra stress into elastic and viscous parts of the poly-meric contributions of the constitutive equation. This approach together with interpo-lation of the velocity-gradient/vorticity or re-structured momentum equation are used for better accuracy. Several methods have been proposed and being widely used in the finite element, finite volume or finite difference literature. These include elastic viscous split stress (EVSS)1,2,3,4 discrete elastic viscous split stress (DEVSS)4,5, explicitly elliptic
momentum equation (EEME)6,7, adaptive viscoelastic stress splitting (AVSS)8,9, EVSS
gradient (EVSS-G)10,11, DEVSS gradient (DEVSS-G/DG)11,12,13, discrete adaptive
vis-coelastic stress splitting gradient (DAVSS-G/DG)14,15, DEVSS vorticity (DEVSS-ω)16,
DAVSS-ω16.
Many of the reported work using the Oldroyd-B model for flow over a circular cylinder, at the higher Deborah numbers, suffer from loss of convergence to steady state on very fine meshes. Phan-Thien and Dou3,16,17 identified that the solution breaks down due to
the stress oscillations in the wake along with large values of primary (first) normal stress difference at the front stagnation point of a cylinder. Fan4 suggested that the solutions
resulted from the higher Deborah numbers may be numerical artifacts. Fan also indicated lack of convergence of the solution beyond 0.7. It was shown by Caola et al.12using the fine
finite element mesh and DEVSS-G based discontinuous Galerkin method, that reaching a steady state beyond a Deborah number of unity is difficult. This is due to the linear increase in the extra stress at the rear stagnation point. The proposed CBS scheme is expected to reduce the stress oscillations at higher Deborah numbers due to the inherent higher order stabilization of convective terms and fractional step based stabilization for pressure variables.
Over the last ten years, a unified computational fluid dynamics solver, referred to as the CBS algorithm, for viscous compressible and incompressible flow calculations has been investigated by Zienkiewicz, Nithiarasu and co-workers 18−23. It is based on a reference
particle following the characteristic direction in the time-space domain to establish the self-adjoint equation leading to temporal semi-discrete formulation. By using a Taylor series expansion, consistent, extra second-order time terms are obtained in a non-moving frame. These terms in addition to increasing the stabilization, also provides second order accuracy in time. On the other hand, in order to circumvent BB or LBB restrictions, the CBS scheme uses fractional step method. Both fully explicit24,25 and semi-implicit
forms26−28 were presented and developed in the past for Newtonian fluids. The extension
of the CBS scheme in its explicit form was previously carried out using the Oldroyd-B model18. In the present paper, we provide details of the extension of the CBS scheme to
using standard simplified, explicit characteristic Galerkin method29.
Three different constitutive models, the Oldroyd-B, upper convected Maxwell (UCM) and Phan-Thien/Tanner (PTT)30 models, have been employed in this work to model
viscoelastic flows. Both the fully explicit and semi-implicit forms are studied and pre-sented for a wide range of Deborah numbers. Section 2 gives the governing equations for the three different viscoelastic fluid models and CBS scheme. In Section 3, qualitative and quantitative results are shown over wide range of Deborah numbers. Finally, the conclusions derived from the present research are described in Section 4.
2 MATHEMATICAL FORMULATIONS
2.1 Governing equations
The non-dimensional isothermal fluid dynamics equations for viscoelastic flows in con-servative form can be written as
∂W ∂t + ∂Fj ∂xj +∂Gj ∂xj = 0 (1)
where the increase time-rate vector
W = (ρ, ρu1, ρu2)T, (2)
the convective acceleration vector Fj = Ã ρuj, δ1jp Re + ρu1uj, δ2jp Re + ρu2uj !T , (3)
and the viscous diffusion vector Gj = µ 0, 1 Re ³ −τ1 1j − τ1j2 ´ , 1 Re ³ −τ1 2j− τ2j2 ´¶T (4) and the Newtonian deviatoric stress tensor
τ1 ij = (1 − α) Ã ∂ui ∂xj + ∂uj ∂xi − 2 3 ∂uk ∂xk δij ! (5) with the non-Newtonian deviatoric stress tensor (τ2
ij = τije + τijv), which is given by the
constitutive equation of extra stress, i.e.
τv ij = α Ã ∂ui ∂xj + ∂uj ∂xi − 2 3 ∂uk ∂xk δij ! (7) In the above governing equations, p is the pressure, ρ is the density, uj are the velocity
components, δij is the Kronecker delta, τije is the polymer-contributed to elastic stress,
τv
ij is the polymer-contributed to viscous stress, α = ηm0/η0 in which ηm0 represents the
polymer-contributed viscosity and η0 represents the zero shear rate viscosity, Re is the
Reynolds number and De is the Deborah number defined as
Re = ρ∞u∞L
η0
; De = λu∞
L (8)
where subscript ∞ indicates a free stream value, L is a characteristic length, λ is the relaxation time and η0 = ηn+ηm0in which ηnrepresents the Newtonian dynamic viscosity.
Apparently, when ε = 0 and 0 < α < 1, the constitutive equation describes the Oldroyd-B model and the upper convected Maxwell (UCM) model is obtained when ε = 0 and
α = 1. Both models explicitly represent the elastic stress (Equation (6)) and the viscous
stress (Equation (7)) and contribute to polymeric solutions as the last term disappears in Equation (6). It should also be noted that the PTT model is obtained by substituting
α = 1 and ε 6= 0.
2.2 Characteristic based split (CBS) formulation
The CBS scheme is based on the characteristic-Galerkin procedure and a fractional step method. It is widely used to carry out fluid dynamic calculations. In this paper, explicit CBS-AC and semi-implicit CBS schemes have been employed to solve viscoelastic flow. Following the intermediate momentum at the first step, the pressure expression is given from the mass conservation at step2. The velocity field is corrected at step3. Finally, the constitutive equations are solved at the fourth step. The semi-discrete form of these steps are written as
Step1: Intermediate momentum ∆U? j = Uj?− Ujn = ∆t " − ∂ ∂xk (ukUj) + 1 Re ∂τ1 ij ∂xi + 1 Re ∂τ2 ij ∂xi #n + (∆t) 2 2 ( um ∂ ∂xm " ∂ ∂xk (ukUj) #)n (9) where Un
Step2: Pressure 1 Re µ1 c2 ¶n ∆p ≈ 1 Re à 1 β2 !n (pn+1− pn) = − ∆t " ∂Un j ∂xj + θ1 ∂∆U? j ∂xj − ∆tθ1 Re à ∂2pn ∂xj∂xj + θ2 ∂2∆p ∂xj∂xj !# (10) where c is the speed of sound which assumes density changes are related to pressure changes for small compressibility or elastic deformability and approaches infinity for in-compressible flows. However, a wave speed c is reasonably replaced by an artificial com-pressibility parameter β to calculate incompressible flow solution. For the explicit CBS-AC scheme, β depends on convective velocity, diffusive velocity and viscoelastic velocity fields.
Step3: Momentum correction
∆Uj = Ujn+1− Ujn= ∆Uj?− ∆t Re ∂p ∂xj n+θ2 (11) where 0.5 ≤ θ1 ≤ 1 and θ2 = 0 is in explicit forms. 0.5 ≤ θ1 ≤ 1 and 0.5 ≤ θ2 ≤ 1 is in
semi-implicit forms. In this work, θ1 = 1 and θ2 = 0 are employed with the fully explicit
scheme whilst θ1 = 1 and θ2 = 1 chosen for the semi-implicit scheme.
Uj = NuU˜j; ∆Uj = Nu∆ ˜Uj; ∆Uj? = Nu∆ ˜U?j; uj = Nuu˜j; p = Npp;˜
∆p = Np∆˜p; τij2 = Nτ2τ˜ij2 (13)
where N are the shape functions and˜indicates a nodal quantity.
Applying the standard Galerkin approximation with the divergence theorem, we get the following weak forms, i.e.
Step1: Weak form of intermediate momentum
Z ΩN T u∆Uj?dΩ = ∆t " − Z ΩN T u ∂ ∂xk (ukUj)dΩ − 1 Re Z Ω ∂NT u ∂xi (τ1 ij + τij2)dΩ #n + ∆t2 2 "Z Ω ∂ ∂xm (umNTu) Ã − ∂ ∂xk (ukUj) ! dΩ #n + ∆t ·Z ΓN T utddΓ ¸n (14) In the above equation td = [τij1/Re]n indicates the part of the traction corresponding
to the Newtonian stresses only and n are the components of the outward normal to the boundaries. As the pressure term is completely removed from the first step, we have only deviatoric stresses part of the traction left in the equation.
Step2: Weak form of pressure equation
explicit scheme Z ΩN T p 1 Re à 1 β2 !n ∆pdΩ = −∆t Z ΩN T p ∂ ∂xj Un jdΩ − ∆t Z ΓN T p à ∆U? j − ∆t Re ∂p ∂xj n! njdΓ + ∆t Z Ω ∂NT p ∂xj à ∆U? j − ∆t Re ∂p ∂xj n! dΩ (15)
In the above equation, pressure and ∆U∗
i terms are integrated by parts and nj are the
components of the outward normal to the boundaries.
Semi-implicit scheme Z ΓN T p 1 Re ∂p ∂xj n+1 njdΓ − Z Ω ∂NT p ∂xj 1 Re ∂p ∂xj n+1 dΩ = 1 ∆t Z ΩN T p ∂U? j ∂xj dΩ (16)
Here, pressure term is integrated by parts.
explicit scheme Z ΩN T u∆UjdΩ = Z ΩN T u∆Uj?dΩ + ∆t Z Ω ∂NT u ∂xj 1 Rep ndΩ − ∆tZ ΓN T utpdΓ (17) Semi-implicit scheme Z ΩN T u∆UjdΩ = Z ΩN T u∆Uj?dΩ + ∆t Z Ω ∂NT u ∂xj 1 Rep n+1dΩ − ∆tZ ΓN T utpdΓ (18)
In the above two equations tp= [(pn+ θ2∆p)/Re]n only indicates the part of the traction
corresponding to the pressure which was removed from step1. It is simply ignored and assumed to be zero as the full traction is prescribed and employed in step1.
Step4: Weak form of constitutive equation
Z ΩN T τ2∆τij2dΩ = ∆t − Z ΩN T τ2 ∂ ∂xk ³ ukτij2 ´ dΩ − Z ΩN T τ2 α + εDe trhτ2 ij i αDe τ 2 ijdΩ n + ∆t "Z ΩN T τ2 à τ2 ik ∂uj ∂xk + τ2 jk ∂ui ∂xk ! dΩ #n + ∆t " α De Z ΩN T τ2 à ∂ui ∂xj +∂uj ∂xi − 2 3 ∂uk ∂uk δij ! dΩ #n + (∆t) 2 2 "Z Ω ∂ ∂xm (umNTτ2) à − ∂ ∂xk ³ ukτij2 ´! dΩ #n − (∆t)2 2 Z Ω ∂ ∂xm (umNTτ2) α + εDe trhτ2 ij i αDe τ 2 ijdΩ n − (∆t)2 2 "Z Ω ∂ ∂xm (umNTτ2) à τ2 ik ∂uj ∂xk + τ2 jk ∂ui ∂xk ! dΩ #n − (∆t)2 2 "Z Ω ∂ ∂xm (umNTτ2) α De à ∂ui ∂xj +∂uj ∂xi −2 3 ∂uk ∂uk δij ! dΩ #n (19) It should be noted that the convection stabilizing treatment indicates last four terms by the characteristic-Galerkin procedure.
Once the finite element approach is worked out by the above weak forms, the matrix expression is convenient to establish, i.e.
Step1: Intermediate momentum ∆ ˜U? = −Mu−1∆t
h
(CuU + K˜ τu + G˜ τ2τ˜2− fu) − ∆t(KuU)˜
in
Step2: Pressure
(Mp+ ∆t2θ1θ2H)∆˜p = ∆t[GuU˜n+ θ1Gu∆ ˜U?− ∆tθ1H˜p − fp]n (21)
Step3: Momentum correction
∆ ˜U = ∆ ˜U?− M−1u ∆thGpT(˜pn+ θ2∆˜p)
i
(22) Step 4: Constitutive equation
∆ ˜τ2 = −M τ2−1∆t h Cτ2τ˜2+ Kτ2τ˜2− Kτvu − D˜ τu˜ in + HOTn (23)
where the velocity, pressure and extra stress variables are approximated using equal in-terpolation functions at all computational points in the domain and are given as
Mu= Z ΩNu TN udΩ; Kτ = Z ΩB T(1 − α) Re µ Io− 2 3mm T ¶ BdΩ; Cu = Z ΩNu T(∇T(uN u))dΩ; Gτ2 = Z Ω 1 Re(∇Nu) TN τ2dΩ; Ku = − 1 2 Z Ω(∇ T(uN u))T(∇T(uNu))dΩ; fu = Z ΓNu Tt ddΓ; Mp = Z ΩNp T à 1 Reβ2 !n NpdΩ; H = Z Ω 1 Re(∇Np) T∇N pdΩ; Gu = Z Ω(∇Np) TN udΩ; fp = ∆t Z ΓNp T · NuU˜n+ θ1 µ ∆ ˜U?− ∆t Re∇p n+θ2 ¶¸ nTdΓ; Gp = Z Ω 1 Re(∇Np) TN udΩ; Cτ2 = Z ΩNτ 2T(∇T(uNτ2))dΩ; Mτ2 = Z ΩNτ 2TNτ2dΩ; Kτ2 = α + εDe trhτ2 ij i αDe Mτ2; Kτv = Z Ω α DeB TI 0mIdΩ; Dτ = Z ΩNτ 2Tτ˜2∇NudΩ; Gu = Z Ω(∇Np) TN udΩ (24)
In the above the strain shape function matrix B is given as
B = SNu (25)
where S is an strain matrix operator derived from the deviatoric stress and strain relations. For a two dimensional case, it can be written as
m = [1, 1, 0]T (27) mI = [1, 1, 1]T (28) and Io = 2 2 1 (29)
3 RESULTS AND DISCUSSION
3.1 Steady state convergence
The steady state convergence criteria is based on the L2 norm of the velocity field. Here
the norm of difference between time step n + 1 and n velocities is normalized by the Euclidean norm of the velocity at time step n + 1. The L2 norm is given as
kek2 = · PN N i=1 ³ kukn+1 i − kukni ´2¸1/2 · PN N i=1 ³ kukn+1 i ´2¸1/2 (30)
where NN are the number of nodes. The above tolerance was reduced to a value of at least 10−7 to assume a steady state.
28R 12R 2R 4R 2R Θ r Θ e u1 e x1 ex2 e
Figure 1: Viscoelastic flow over a circular cylinder placed in a channel.
3.2 Flow over a circular cylinder placed in a channel
(a) (b)
(c) (d)
(e)
(f)
Figure 2: Viscoelastic flow over a circular cylinder. (a) Unstructured mesh1 (∆d = 0.097); (b) Unstruc-tured mesh2 (∆d = 0.041); (c) UnstrucUnstruc-tured mesh3 (∆d = 0.008); (d) Hybrid mesh4 (∆d = 0.001); (e) Unstructured mesh3. Nodes: 27451; Elements: 53516; (f) Hybrid mesh4. Nodes: 35506; Elements: 69918.
inlet to exit section is 28R. The geometry and boundary conditions are shown in Figure 1.
Figure 2(a)-(e) shows the four meshes in the vicinity of the cylinder. All meshes are refined close to the wall. The element size near the cylinder wall are respectively 0.097, 0.041, 0.008 and 0.001 for the four meshes. The hybrid mesh4 consists of 22 structured layers to capture high elastic stress contribution around the cylinder surface. The mesh of the full domain for the unstructured mesh3 and hybrid mesh4 are shown in Figure 2(f) and (g) respectively.
Mesh Typical element size Drag force ∆d D Unstructured mesh1 0.097 119.23 Unstructured mesh2 0.041 116.83 Unstructured mesh3 0.008 114.96 Hybrid mesh4 0.001 117.54
Table 1: Meshes convergence at Re = 0 and De = 3 by using explicit form with the Oldroyd-B model.
to be zero which assumes convective terms of the momentum equation are removed. The boundary conditions in non-dimensional form are given below3,16,17. The inlet and
exit velocity profiles and extra stress are
u1 = 1.5 Ã 1 −x22 4 ! (31) u2 = 0 (32) τ2 11 = 9 8(De)x 2 2 (33) τ2 12= − 3 4x2 (34) and τ2
22 = 0. On the channel solid walls and cylinder surface, no slip conditions are
assumed. The pressure is zero at exit for semi-implicit scheme. The extra stresses on the channel walls for Oldroyd-B fluid are given as
τ2 11= 0.82De à ∂u1 ∂x2 !2 (35) τ2 12 = 0.41 ∂u1 ∂x2 (36) and τ22 = 0. In cylindrical coordinates, the extra stresses on the cylinder surface are
and τrr = 0. It should be noted, however, that no need to explicitly prescribe the extra
stresses on the channel walls and on cylinder surface for UCM and PTT models.
The detail of mesh convergence at a Deborah number of 3 using the Oldroyd-B model are given in Table 1. It shows that the drag forces converge consistently on the unstruc-tured meshes. However, the hybrid mesh with finer elements show a higher drag force than the converged solution on the unstructured mesh.
The drag force around the circular cylinder per unit length is calculated in non-dimensional form as Fx = Z 2π 0 h (−p + τ111 + τ112)cosθ + (τ121 + τ122)sinθiRdθ (39) (a) (b) (c) (d) (e) (f)
Figure 3: Oldroyd-B fluid flow over a circular cylinder at De = 1 using both fully explicit scheme (left)
and the semi-implicit scheme (right). (a) τ2
11contours. τ112min = -4.07, τ 2 11max = 191.59; (b) τ 2 11contours. τ2 11min = -0.81, τ 2 11max = 159.80; (c) τ 2 12 contours. τ122min = -38.77, τ 2 12max = 42.86; (d) τ 2 12 contours. τ2 12min = -45.73, τ 2 12max = 48.75; (e) τ 2 22 contours. τ222min = -0.76, τ 2 22max = 19.63; (f) τ 2 22contours. τ222min = -0.38, τ2 22max = 21.28.
Figure 3 shows the extra stress contours at De = 1 from explicit and semi-implicit forms using the Oldroyd-B model. As seen, the large normal stress in the x1 direction occurs on
the top and bottom cylinder surfaces. Both schemes give some minor spatial oscillations of τ2
11 in the wake, but the semi-implicit scheme is more stable than explicit scheme.
110 120 130 140 150 160 170 0 0.5 1 1.5 2 2.5 3 Drag force, D Deborah number, De CBS scheme; Oldroyd-B fluid
Explicit; Hybrid mesh4 Semi-implicit; Hybrid mesh4 SMART; Alves et al.[30] DAVSS-G/DG; Sun et al.[14] GLS-MIX1; Fan et al.[4] DAVSS-ω; Dou et al.[16] Extrapolated; Dou et al.[16] DEVSS-G/DG; Caola et al.[12]
(a) 76 78 80 82 84 86 88 90 92 94 0 0.5 1 1.5 2 2.5 3
Pressure component of drag force, D
p
Deborah number, De CBS scheme; Oldroyd-B fluid
Explicit; Hybrid mesh4 Semi-implicit; Hybrid mesh4 DEVSS-G/DG; Caola et al.[12]
(b) 32 34 36 38 40 42 44 0 0.5 1 1.5 2 2.5 3
Viscoelastic stress component of drag force, D
s
Deborah number, De CBS scheme; Oldroyd-B fluid
Explicit; Hybrid mesh4 Semi-implicit; Hybrid mesh4 DEVSS-G/DG; Caola et al.[12]
(c)
Figure 4: Oldroyd-B fluid flow over a circular cylinder. (a) Drag force, D, distribution; (b) Pressure com-ponent of drag force, Dp, distribution; (c) Viscoelastic stress comcom-ponent of drag force, Ds, distribution.
(a) (b)
(c) (d)
Figure 5: Oldroyd-B fluid flow over a circular cylinder at De = 3 using both fully explicit scheme (left)
and the semi-implicit scheme (right). (a) u1contours. u1min = 0.0, u1max = 3.17; (b) u1 contours. u1min
= 0.0, u1max = 3.15; (c) u2contours. u2min= -0.97, u2max = 1.07; (d) u2contours. u2min = -0.94, u2max
= 1.04.
Figure 4(a) shows the comparison of the variation of the drag force at different Deborah numbers. The results show an excellent comparison with the literature up to a Deborah number of unity. However, beyond this limit the drag force did not increase dramatically like the ones published previously. The pressure and viscoelastic stress (including New-tonian and non-NewNew-tonian stresses) components of drag force are also depicted in Figure 4(b) and (c).
Figure 5 shows contours of horizontal velocity component u1, vertical velocity
compo-nent u2 and pressure at De = 3 using both the variations of the CBS schemes. As seen
Explicit Semi-implicit
De (Mesh4) (Mesh4) Caola et al.[12] Sun et al.[14] Fan et al.[4] Dou et al.[16] Alves et al.[31] 1.0 117.58 117.15 117.80 119.32 118.49 118.78 118.69 1.1 117.38 117.20 118.00 120.39 119.69 1.2 117.27 117.25 121.65 120.50 1.3 117.16 117.30 123.07 121.54 1.4 117.21 117.35 124.66 122.86 1.5 117.26 117.39 126.32 1.6 117.27 117.45 128.01 1.7 117.30 117.48 129.67 1.8 117.33 117.49 131.37 1.85 132.30 1.9 117.36 117.54 2.0 117.38 117.58 2.1 117.41 117.62 2.2 117.43 117.66 2.3 117.46 117.70 2.4 117.47 117.75 2.5 117.49 117.80 2.6 117.51 117.84 2.7 117.51 117.88 2.8 117.54 117.92 2.9 117.54 117.95 3.0 117.54 118.00
Table 2: Drag force calculations by using CBS schemes with the Oldroyd-B model at higher Deborah numbers.
The extra stress contours at the Deborah number of three are shown in Figure 6. As seen, the minor spatial oscillations of τ2
11 still exists in the wake. Some oscillations are
also noted in τ2
12 contours. When the Deborah number is increased, the maximum extra
stress components τ2
11, τ222 and τ122 , becomes larger due to a very thin shear boundary layer
around the cylinder surface. Table 2 presents the drag force values for different Deborah numbers up to a value of three.
The drag forces from UCM and PTT models are compared with Alves et al.31, Fan et
al.4 and Phan-Thien et al.3’ in Tables 3 and 4.
4 CONCLUSIONS
Phan-Explicit Semi-implicit
De (Mesh3) (Mesh4) (Mesh3) (Mesh4) Alves et al.[31] Fan et al.[4] Phan-Thien et al.(M4)[3] 0.0 131.79 132.22 131.84 132.29 132.34 132.36 131.91 0.00164 131.55 0.0162 131.24 0.01 132.13 0.025 131.66 131.81 130.80 131.65 131.82 0.0323 130.02 0.05 130.68 131.27 130.78 0.075 129.14 129.27 129.27 129.14 0.1 128.11 128.17 127.30 127.15 127.22 127.42 0.129 125.31 0.2 118.96 116.50 118.24 116.58 117.65 117.83 0.248 115.14 0.3 111.70 109.50 109.14 107.70 108.52 108.68 0.323 109.72 0.4 102.66 101.13 101.29 101.43 0.437 104.31
Table 3: Drag force calculations by using CBS schemes with the UCM model.
Thien/Tanner (PTT) models were studied in this paper. One of the major problems of Oldrod-B model is the convergence to steady state beyond a Deborah number of unity. In this paper we have proved that the CBS scheme gives steady state results up to a Deborah number of three. The major conclusions derived are:
• The CBS scheme is validated for Oldroyd-B, UCM and PTT fluids over a wide range
Deborah numbers. It is clear that most of the models agree excellently at Deborah numbers below unity.
• At higher Deborah numbers, CBS scheme gives converged solution but the
agree-ment with other numerical solution is poor. This could be due to the lack of con-vergence reported by other researchers.
• The CBS scheme can work for UCM and PTT fluids over a circular cylinder even
though the majority part of the scheme are explicit in nature.
• At lower Deborah numbers, the explicit scheme is faster than the semi-implicit
scheme. However, at higher Deborah numbers, the semi-implicit scheme is faster. This conclusion could be different if the extra artificial dissipation is switched on.
• The lack of extra dissipation results in some minor extra stress oscillations. Further
Explicit Semi-implicit Phan-Thien et al.[3] De (Mesh3) (Mesh4) (Mesh3) M4 M5 0.0 131.79 132.22 131.84 131.91 131.60 0.00164 131.49 131.20 0.0162 130.28 130.20 0.025 128.69 128.79 0.0323 126.59 126.55 0.05 121.40 119.10 121.50 0.0647 116.31 116.13 0.075 113.44 111.27 0.1 109.09 107.00 106.14 0.129 98.38 98.04 0.162 91.47 91.19 0.2 84.65 83.02 84.73 0.248 77.80 77.20 0.3 71.85 70.44 73.40 0.323 69.44 69.35 0.4 64.55 63.26 0.5 59.88 58.67 0.6 55.54 54.38 0.611 51.18
Table 4: Drag force calculations by using CBS schemes with the PTT model.
5 Acknowledgements
This research is funded by the Engineering and Physical Sciences Research Council (EPSRC) grant no. EP/C515498/1.
REFERENCES
[1] M.G.N. Perera and K. Walters, Long-range memory effects in flows involving abrupt changes in geometry: Part I: flows associated with I-shaped and T-shaped geometries.
Journal of Non-Newtonian Fluid Mechanics, 2, 49–81 (1977).
[2] D. Rajagopalan, R.C. Armstrong and R.A. Brown, Finite element methdos for cal-culation of steady, viscoelastic flow using constitutive equations with a Newtonian viscosity. Journal of Non-Newtonian Fluid Mechanics, 36, 159–192 (1990).
[3] N. Phan-Thien and H.S. Dou, Viscoelastic flow past a cylinder: drag coefficient.
(a) (b)
(c) (d)
(e) (f)
(g) (h)
Figure 6: Oldroyd-B fluid flow over a circular cylinder at De = 3 using both fully explicit scheme (left)
and the semi-implicit scheme (right). (a) Primary normal stress difference contours. N1min = -55.84,
N1max = 647.35; (b) Primary normal stress difference contours. N1min = -24.11, N1max = 593.37; (c)
τ2 11contours. τ112min = -5.42, τ 2 11max = 646.92; (d) τ 2 11 contours. τ112min = -0.82, τ 2 11max = 594.15; (e) τ 2 12 contours. τ2 12min = -112.56, τ 2 12max = 124.39; (f) τ 2 12 contours. τ122min = -137.52, τ 2 12max = 147.96; (g) τ 2 22 contours. τ2 22min = -0.75, τ 2 22max = 62.26; (h) τ 2 22 contours. τ222min = -0.38, τ 2 22max = 56.49.
[4] Y. Fan, R.I. Tanner and N. Phan-Thien, Galerkin/least-square finite-element meth-ods for steady viscoelastic flows. Journal of Non-Newtonian Fluid Mechanics, 84, 233–256 (1999).
[5] R. Gu´enette and M. Fortin, A new mixed finite element method for computing viscoelastic flows. Journal of Non-Newtonian Fluid Mechanics, 60, 27–52 (1995). [6] R.C. King, M.R. Apelian, R.C. Armstrong and R.A. Brown, Numerically stable finite
element techniques for viscoelastic calculations in smooth and singular geometries.
[7] D. Rajagopalan, R.C. Armstrong and R.A. Brown, Calculation of steady viscoelastic flow using a multimode Maxwell model: application of the explicitly elliptic momen-tum equation (EEME) formulation. Journal of Non-Newtonian Fluid Mechanics, 36, 135–157 (1990).
[8] J. Sun, N. Phan-Thien and R.I. Tanner, An adaptive viscoelastic stress splitting scheme and its applications: AVSS/SI and AVSS/SUPG. Journal of Non-Newtonian
Fluid Mechanics, 65, 75–91 (1996).
[9] V. Warichet and V. Legat, Adaptive high-order prediction of the drag correction factor for the upper-convected Maxwell fluid. Journal of Non-Newtonian Fluid
Me-chanics, 73, 95–114 (1997).
[10] K. Foteinopoulou, V.G. Mavrantzas and J. Tsamopoulos Numerical simulation of bubble growth in Newtonian and viscoelastic filaments undergoing stretching. Journal
of Non-Newtonian Fluid Mechanics, 122, 177–200 (1997).
[11] A.W. Liu, D.E. Bornside, R.C. Armstrong and R.A. Brown, Viscoelastic flow of poly-mer solutions around a periodic, linear array of cylinders: comparisons of predictions for microstructure and flow fields. Journal of Non-Newtonian Fluid Mechanics, 77, 153–190 (1998).
[12] A.E. Caola, Y.L. Joo, R.C. Armstrong and R.A. Brown, Highly parallel time integra-tion of viscoelastic flows. Journal of Non-Newtonian Fluid Mechanics, 100, 191–216 (2001).
[13] F.P.T. Baaijens, An iterative solver for the DEVSS/DG method with application to smooth and non-smooth flows of the upper convected Maxwell fluid. Journal of
Non-Newtonian Fluid Mechanics, 75, 119–138 (1998).
[14] J. Sun, M.D. Smith, R.C. Armstrong and R.A. Brown, Finite element method for viscoelastic flows based on the discrete adaptive viscoelastic stress splitting and the discontinuous Galerkin method: DAVSS-G/DG. Journal of Non-Newtonian Fluid
Mechanics, 86, 281–307 (1999).
[15] X.-J. Fan, N. Phan-Thien and R. Zheng, Simulation of fibre suspension flows by the Brownian configuration field method. Journal of Non-Newtonian Fluid Mechanics, 84, 257–274 (1999).
[16] H.S. Dou and N. Phan-Thien, The flow of an Oldroyd-B fluid past a cylinder in a channel: adaptive viscosity vorticity (DAVSS-ω) formulation. Journal of
[17] H.-S. Dou and N. Phan-Thien, Parallelisation of an unstructured finite volume code with PVM: viscoelastic flow around a cylinder. Journal of Non-Newtonian Fluid
Mechanics, 77, 21–51 (1998).
[18] P. Nithiarasu, A fully explicit characteristic based split (CBS) scheme for viscoelastic flow calculations. International Journal for Numerical Methods in Engineering, 60, 949–978 (2004).
[19] O.C. Zienkiewicz and R. Codina, A general algorithm for compressible and incom-pressible flow –part I. the split, characteristic-based scheme. International Journal
for Numerical Methods in Fluids, 20, 869–885 (1995).
[20] O.C. Zienkiewicz and P. Nithiarasu, A universal algorithm for fluid dynamics. The Characteristic Based Split (CBS) procedure. Some tests on stability and boundary conditions. Archives of Mechanics, 52, 857–887 (2000).
[21] N. Massarotti, P. Nithiarasu and O.C. Zienkiewicz, Characteristic-Based-Split (CBS) algorithm for incompressible flow problems with heat transfer. International Journal
of Numerical Methods for Heat & Fluid Flow 8, 969–990 (2001).
[22] P. Nithiarasu, J.S. Mathur, N.P. Weatherill and K. Morgan, Three dimensional incompressible flow calculations using the Characteristic Based Split (CBS) scheme.
International Journal for Numerical Methods in Fluids, 44, 1207–1229 (2004).
[23] P. Nithiarasu, R. Codina and O.C. Zienkiewicz, The Characteristic Based Split (CBS) Scheme − a unified approach to fluid dynamics. International Journal for
Numerical Methods in Engineering, (In press 2006).
[24] O.C. Zienkiewicz, B.V.K. Satya Sai, K. Morgan, R. Codina and M. V´azquez, A general algorithm for compressible and incompressible flow –part II. tests on the explicit form. International Journal for Numerical Methods in Fluids, 20, 887–913 (1995).
[25] P. Nithiarasu, An efficient Artificial Compressibility (AC) scheme based on Char-acteristic Based Split (CBS) method for incompressible flows. International Journal
for Numerical Methods in Engineering, 56, 1815–1845 (2003).
[26] O.C. Zienkiewicz, P. Nithiarasu, R. Codina, M. V´azquez and P. Ortiz, The characteristic-based-split procedure: an efficient and accurate algorithm for fluid problems. International Journal for Numerical Methods in Fluids, 31, 359–396 (1999).
[28] R. Codina, M. V´azquez and O.C. Zienkiewicz, A general algorithm for compressible and incompressible flow –part III. The semi-implicit form. International Journal for
Numerical Methods in Fluids, 27, 13–32 (1998).
[29] O.C. Zienkiewicz, R.L. Taylor and P. Nithiarasu, The Finite Element Method for
Fluid Dynamics, 5th Edition, Elsevier, 2006.
[30] N. Phan-Thien and R.I. Tanner, A new constitutive equation derived from network theory. Journal of Non-Newtonian Fluid Mechanics, 2, 353–365 (1977).
[31] M.A. Alves, F.T. Pinho and P.J. Oliveira, The flow of viscoelastic fluids past a cylinder: finite-volume high-resolution methods. Journal of Non-Newtonian Fluid
