Effect of fluid elasticity on the numerical stability of high-resolution schemes for high shearing contraction flows using OpenFOAM

Theoretical & Applied Mechanics Letters (

)

–

Contents lists available at ScienceDirect

Theoretical & Applied Mechanics Letters journal homepage: www.elsevier.com/locate/taml

Letter

Effect of fluid elasticity on the numerical stability of high-resolution schemes for high shearing contraction flows using OpenFOAM T. Chourushi Department of Mechanical Engineering, Indian Institute of Technology Indore, Simrol, Indore - 453 552, India

article

info

Article history: Received 12 November 2016 Accepted 28 December 2016 Available online xxxx *This article belongs to the Fluid Mechanics Keywords: High resolution schemes (HRS) Viscoelastic fluid Contraction flows Elasticity number OpenFOAM

a b s t r a c t Viscoelastic fluids due to their non-linear nature play an important role in process and polymer industries. These non-linear characteristics of fluid, influence final outcome of the product. Such processes though look simple are numerically challenging to study, due to the loss of numerical stability. Over the years, various methodologies have been developed to overcome this numerical limitation. In spite of this, numerical solutions are considered distant from accuracy, as first-order upwind-differencing scheme (UDS) is often employed for improving the stability of algorithm. To elude this effect, some works been reported in the past, where high-resolution-schemes (HRS) were employed and Deborah number was varied. However, these works are limited to creeping flows and do not detail any information on the numerical stability of HRS. Hence, this article presents the numerical study of high shearing contraction flows, where stability of HRS are addressed in reference to fluid elasticity. Results suggest that all HRS show some order of undue oscillations in flow variable profiles, measured along vertical lines placed near contraction region in the upstream section of domain, at varied elasticity number E ≈ 5. Furthermore, by E, a clear relationship between numerical stability of HRS and E was obtained, which states that the order of undue oscillations in flow variable profiles is directly proportional to E. © 2017 The Author(s). Published by Elsevier Ltd on behalf of The Chinese Society of Theoretical and Applied Mechanics. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

1. Introduction Viscoelastic fluids are special category of non-Newtonian fluids which show viscous, elastic and non-linear responses to the flows. Such characteristics make these fluids very useful in daily routine works such as laundry liquids, glues, gels, oils, sauces, paints, etc. Also, these fluids are abundantly used in process and polymer industries. Typical industrial application includes extrusion and mould flow processes. In these processes, final outcome of the product depends significantly on the property and pattern of flow [1]. As these processes offer abrupt contraction to the flow of fluid, it gives rise to stress singularity near the re-entrant corner. At this point, the stress term rises exponentially to very high value. Hence, this problem is widely studied in past, both experimentally and numerically [2,3]. With the advancement in conventional rheometers, experimental works offer comparatively an easy prediction of flow properties [4]. However, this problem is numerically difficult to study due to the non-linear nature of constitutive equations [5]. Also, these works are subjected to loss of numerical stability at high Deborah numbers, commonly referred as high Deborah number E-mail addresses: [email protected], [email protected].

problem (HDNP) [6]. To overcome this elevated issue of HDNP, numerous numerical methodologies have been put forward [7–14]. Amid this, some of the famous methodologies which can be applied to finite volume technique are both-side-diffusion (BSD) [11], log-conformation representation (LCR) [12], positive definiteness preserving scheme (PDPS) [13] and square-root conformation representation (SRCR) [14]. In the past, Jovani et. al [15] also presented split-stress approach in reference to OpenFOAM which resembles the BSD approach. This approach promotes numerical stability of the algorithm by adding additional diffusive terms on both sides of the momentum equation. While, other approaches promote numerical stability of the algorithm, by preserving positive definiteness of conformation tensor on discrete level [16]. As numerical stability of the algorithm is often considered important, accuracy of the results are usually overlooked for viscoelastic fluids. This problem was first addressed by Alves et al. [10], in which first-order upwind-differencing scheme (UDS) and second-order numerical schemes were compared. Broadly, accuracy of numerical solution depends heavily on the order of convective schemes [17]. These schemes are classified according to its polynomial relationship of control volume grid points. Numerical schemes relying only on the immediate neighbouring control volume grid point are referred as lower-order schemes. However,

http://dx.doi.org/10.1016/j.taml.2017.01.005 2095-0349/© 2017 The Author(s). Published by Elsevier Ltd on behalf of The Chinese Society of Theoretical and Applied Mechanics. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

Please cite this article in press as: T. Chourushi, Effect of fluid elasticity on the numerical stability of high-resolution schemes for high shearing contraction flows using OpenFOAM, Theoretical & Applied Mechanics Letters (2017), http://dx.doi.org/10.1016/j.taml.2017.01.005.

2

T. Chourushi / Theoretical & Applied Mechanics Letters (

schemes relying on large number of control volume grid points are mentioned as higher-order (HO) schemes. These unbounded schemes though accurate, give rise to undue oscillations of numerical solutions for convection dominated flows. These unwanted oscillations can be suppressed by imposing any of these conditions: (1) total variation diminishing (TVD) [18], (2) convection boundedness criterion (CBC) [19]. To further improve the convergence and stability of numerical schemes, a special class of bounded schemes are used which are commonly known as high-resolution schemes (HRS) [19–23]. Over the last decade, some works were reported using HRS for viscoelastic fluids [6,15,16,24–26]. In broad terms, these works assessed the numerical stability of algorithm for different Deborah numbers. Also, these studies are limited to low shearing flows or creeping flows. Consequently, insufficient information is available on the numerical stability of high shearing flows. This article therefore presents a study of high shearing flows using the existing viscoelastic fluid flow solver [15] and HRS, wherein elasticity number is varied, in order to assess the numerical stability of HRS. In the present context of this paper, HRS are implemented using the blending factor (γ ) method [21] in OpenFOAM. The article is presented as follows: The governing equations for isothermal, incompressible and viscoelastic fluid are detailed in Sect. 2.1, followed by a brief description of BSD approach, numerical discretization and numerical algorithm in Sects. 2.2, 2.3 and 2.4, respectively. In Sect. 2.3, an overview of normalized variable approach and HRS are presented. Later in Sect. 3, flow geometry and mesh characteristics of planar contraction is outlined. Thereafter, under the results and discussion section, a comparative study of undue oscillations of flow variable profiles for HRS are detailed and assessment of numerical stability of HRS for varied elasticity number (E) are also presented, respectively. Finally, conclusions are reported based on these results which details the numerical stability of HRS in reference to E. 2. Numerical methodology 2.1. Governing equations The governing mass and momentum conservation equations for isothermal, incompressible and viscoelastic fluid in vector form are written as:

∇ · u = 0, (1) ∂ (ρ u) + ∇ · (ρ uu) = −∇ p + ∇ · τ, (2) ∂t where u refers to the velocity vector, ρ refers to the density, p refers to the pressure and τ refers to the total stress. This total stress can be decomposed into viscous and elastic parts:

τ = τS + τP,

(3)

where τ S denotes the solvent (viscous) contribution and τ P denotes the polymeric (elastic) contribution. The viscous term τ S is further detailed as under:

[ ] τ S = ηS ∇ u + (∇ u)T .

(4)

In this equation, ηS refers to the solvent viscosity. And, the polymeric term τ P referred in Eq. (3), depends upon the viscoelastic fluid model. As the simplest viscoelastic models (UCM and Oldroyd-B) does not predict the shear-thinning characteristics of the fluid, so the single-mode Giesekus model [27] is being considered for formulation of constitutive equation. Moreover, this model enables decent prediction of first and second stress differences for the mobility factor α in the range of 0 < α < 0.5. ∇

τ P + λτ P + α

[ ] λ (τ P · τ P ) = ηP ∇ u + (∇ u)T , ηP

(5)

)

–

where ηP represents the polymeric contribution of shear viscosity, λ represents the relaxation time, α represents the mobility factor and (τ P · τ P ) represents the non-linear term, which corresponds to qualitative description of viscoelastic properties [27–29]. And, ∇

τ P refers to the upper convected derivative which is defined as under [30]: ∇

τP =

Dτ P Dt

− (∇ u)T · τ P − τ P · ∇ u.

(6)

By substituting Eq. (6) in Eq. (5), following constitutive equation is obtained:

α δτ P τP + + ∇ · (uτ P ) + (τ P · τ P ) λ δt ηP ] ηP [ = ∇ u + (∇ u)T + (∇ u)T · τ P + τ P · ∇ u. λ

(7)

Thus, by substituting Eqs. (3), (4) and (7) into Eq. (2), give rise to the complete governing momentum equation for viscoelastic fluids.

[ ] ∂ (ρ u) + ∇ · (ρ uu) − ∇ · ηS ∇ u + (∇ u)T = −∇ p + ∇ · τ P . (8) ∂t These equations though complete in all aspects, are subjected to the loss of numerical stability when solvent viscosity (ηS ≈ 0) reduces to negligibly small value. To overcome this shortcoming, both-side-diffusion (BSD) approach is considered in the present study [11]. 2.2. Both-side-diffusion approach In momentum equation, non-linear hyperbolic stress term (τ P ) is calculated based on the constitutive equation. When ηS → 0, elliptic diffusive term on the left-hand-side (LHS) of Eq. (8), reduces to zero. This absence of an explicit diffusive term makes these equations difficult to converge for high shearing flows. In order to overcome this limitation, an additional diffusion term is added on both sides of momentum equation. Thus, Eq. (8) changes to:

[ ] ∂ (ρ u) + ∇ · (ρ uu) − ∇ · (ηS + ηP ) ∇ u + (∇ u)T ∂t [ ] = −∇ p − ∇ · ηP ∇ u + (∇ u)T + ∇ · τ P .

(9)

This approach of adding diffusive terms is commonly referred as BSD approach. In above equation, all terms on the LHS are implicitly treated and remaining right-hand-side (RHS) terms are explicitly added as a source. In case of numerical convergence, these added terms cancel each other. 2.3. Numerical discretization In this section, the governing flow equation (referred in Eq. (9)) are integrated over a control volume and gauss divergence theorem is applied over it, to give a set of discretized equations [17]. In this equation, unsteady, time-dependent terms are discretized using the Crank Nicolson scheme and convective terms of momentum and constitutive equations are discretized using the high resolution scheme (HRS). While, remaining terms are discretized using the central difference scheme. All flow variables are stored using the non-staggered, collocated grid arrangement of OpenFOAM. And localize normalized variable approach [21] is being considered for formulation of HRS. 2.3.1. Normalized variable approach Figure 1 shows the actual plot of convected variables and distances of the control volume, along positive flow direction. These variables and distances are normalized using the localize normalized variable approach for any arbitrary control volume mesh.

Please cite this article in press as: T. Chourushi, Effect of fluid elasticity on the numerical stability of high-resolution schemes for high shearing contraction flows using OpenFOAM, Theoretical & Applied Mechanics Letters (2017), http://dx.doi.org/10.1016/j.taml.2017.01.005.

T. Chourushi / Theoretical & Applied Mechanics Letters (

a

)

–

3

b

Fig. 1. Convected variables along the positive flow direction: (a) actual plot of convected variables and distances, and (b) arbitrary control volume mesh.

Wherein, variable φ of central grid point C of the control volume is bounded between its adjoining interfaces f − and f + as follows:

φˆ C =

φC − φf− φf+ − φf−

=1−

φf+ − φC φf+ − φf−

.

+

fx = (fD) /(CD) .

(11)

Based on Eqs. (10) and (11), the numerator term (φf+ − φC ) is given as under

φf+ − φC = fx (φD − φC ) = fx

φD − φC (xD − xC ) , xD − xC

(12)

where xD and xC represents distance for downwind (D) and central (C ) grid points of the control volume, respectively. Hence, Eq. (12) reduces to:

φD − φC = (∇φ)f · dˆ . xD − xC

(13)

By substituting Eq. (13) in Eq. (12) we have,

φf+ − φC = fx (∇φ)f · dˆ (xD − xC ) .

(14)

Here, dˆ refers to the unit vector, as depicted in Fig. 1(b). Similarly, denominator term of Eq. (10) can be written as:

φf+ − φf− =

φf+ − φf− ( +

−

xf − xf

φˆ C = 1 −

(10)

Here, φf and φf are interpolated control volume face variables. For any arbitrary control volume mesh (referred in Fig. 1(b)), interpolation factor is written as −

Eq. (18) can be transformed as,

− − x+ = (∇φ)C · dˆ x+ . f − xf f − xf

)

(

)

(15)

φD − φC 2(∇φ)C · d

.

(20)

Thus, Eq. (20) represents localized normalized convected variable for the central grid point (φˆ C ) of control volume. For any numerical scheme, its functional relationship is formulated in reference to φˆ C , as follows

( ) φˆ f = f φˆ C , xˆ C , xˆ f .

(21)

Here, xˆ C and xˆ f represent normalized distances for central grid point and face of the control volume, respectively. 2.3.2. High resolution schemes Table 1 shows the normalized functional relationships of HRS considered in this work. Figure 2 details the plot of these schemes in NVD. These normalized convected face variable φˆ f is then unnormalized using the blending factor (γ ), which is calculated in reference to HRS. Since, far-upwind convected variable (φU ) is not considered for normalization, unnormalized face variable φf is reconstructed as below

φf = [1 − γ (1 − fx )] φC + γ (1 − fx ) φD .

(22)

The other control volume faces are calculated using the similar procedure. These convected face variables φf are then substituted for convective terms, to give the discretized governing equation. As, these variables are constituted from neighbouring grid points of the control volume face, it generates compact implicit tri-diagonal system of linear equations.

By substituting Eqs. (14) and (15) in Eq. (10) we have,

φˆ C = 1 −

φf+ − φC φf+ − φf−

=1−

fx (∇φ)f · dˆ (xD − xC )

). −

(∇φ)C · dˆ x+ f − xf

(

2.4. Numerical algorithm (16)

Since, grid point C is located in middle of the control volume cell, following geometric relationship holds: − x+ f − xf

xD − xC

=2

x+ f − xC xD − xC

= 2fx .

(17)

AC U = H (U ) − ∇ p,

By substituting Eq. (17) in Eq. (16) we have,

φˆ C = 1 −

(∇φ)f · d (∇φ)f · dˆ =1− . 2(∇φ)C · d 2(∇φ)C · dˆ

(18)

As difference between convected variables of control volume grid points is,

φD − φC = (∇φ)f · d .

These governing equations can be written into semi-discretized form of momentum equation, which follows the spirit of Rhie– Chow interpolation [31,32] where pressure gradient term is not discretized at this stage. Thus, the semi-discretized form of momentum equation in reference to Jovani et al. [15] is as under,

(19)

(23)

where AC refers to the contribution of C , H (U ) comprises the contribution of neighbouring grid points and other added source terms. And, last term on the right hand side of equation denotes the pressure gradient. Furthermore, velocity field can be obtained by dividing Eq. (23) with AC , as follows: U =

H (U ) AC

−

1 AC

∇ p.

(24)

Please cite this article in press as: T. Chourushi, Effect of fluid elasticity on the numerical stability of high-resolution schemes for high shearing contraction flows using OpenFOAM, Theoretical & Applied Mechanics Letters (2017), http://dx.doi.org/10.1016/j.taml.2017.01.005.

4

T. Chourushi / Theoretical & Applied Mechanics Letters (

)

–

Table 1 Normalized functional relationships of MINMOD, GAMMA, SMART, WACEB and CUBISTA schemes in reference to uniform/non-uniform grids [19–23]. HR schemes

MINMOD

Uniform grids

φˆ f =

⎧3 ⎪ φˆ C , ⎪ ⎪ ⎨2 1

Non-uniform grids 0 < φˆ C < 1

φˆ C + , ⎪ ⎪ 2 ⎪ ⎩2 φˆ C ,

GAMMA

φˆ f =

1 2

1 2

φˆ f =

≤ φˆ C < 1

elsewhere

] ⎧[ 1 ⎪ 1+ (1 − φˆ C ) φˆ C , ⎪ ⎪ 2βm ⎨ 1

0 < φˆ C < βm

1

φˆ C + , ⎪ ⎪ 2 ⎪ ⎩2

βm ≤ φˆ C < 1

φˆ C ,

SMART

φˆ f =

⎧ ⎪ 3φˆ C , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪3 ⎨ 4

φˆ C +

3

,

⎪ ⎪ ⎪ ⎪ 1, ⎪ ⎪ ⎪ ⎩

φˆ C ,

WACEB

⎧ ⎪ 2φˆ C , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ 3 φˆ + 3 , C 8 φˆ f = 4 ⎪ ⎪ ⎪ ⎪ 1, ⎪ ⎪ ⎪ ⎩ φˆ C ,

CUBISTA

⎧7 ⎪ φˆ C , ⎪ ⎪ 4 ⎪ ⎪ ⎪3 ⎪ ⎨ φˆ + 3 , C 8 φˆ f = 4 ⎪ 3 ⎪1 ˆ ⎪ ⎪ φC + , ⎪ ⎪ 4 4 ⎪ ⎩ φˆ C ,

elsewhere 0 < φˆ C <

8

φˆ f =

1 6 5 6

≤ φˆ C ≤

1 6 5 6

φˆ f =

< φˆ C < 1

elsewhere 0 < φˆ C < 3 10 5 6

3 10 5

≤ φˆ C ≤

6

φˆ f

< φˆ C < 1

elsewhere 0 < φˆ C < 3 8 3 4

≤ φˆ C ≤

3 8 3 4

< φˆ C < 1

elsewhere

φˆ f

⎧ xˆ f ⎪ ⎪ φˆ C , ⎪ ⎪ ⎨ xˆ C

0 < φˆ C < xˆ C

1 − xˆ f

xˆ f − xˆ C

1 − xˆ f

xˆ f − xˆ C

φˆ C + , xˆ C ≤ φˆ C < 1 ⎪ ⎪ 1 − xˆ C 1 − xˆ C ⎪ ⎪ ⎩ˆ φC , elsewhere ⎧[ ] xˆ f − xˆ C ⎪ ⎪ 1+ (1 − φˆ C ) φˆ C , 0 < φˆ C < βm ⎪ ⎪ ˆ β (1 − x ) ⎨ m C φˆ C + , ⎪ ⎪ 1 − xˆ C 1 − xˆ C ⎪ ⎪ ⎩ˆ φC , ⎧ xˆ f (1 − 3xˆ C + 2xˆ f ) ⎪ ⎪ φˆ C , ⎪ ⎪ xˆ C (1 − xˆ C ) ⎪ ⎪ ⎪ ⎪ xˆ f (xˆ f − xˆ C ) ⎨ xˆ f (1 − xˆ f ) ˆ xˆ (1 − xˆ )

φC +

1 − xˆ

βm ≤ φˆ C < 1 elsewhere

,

C C C ⎪ ⎪ ⎪ ⎪ 1 , ⎪ ⎪ ⎪ ⎪ ⎩ˆ φC , ⎧ ⎪ ⎪ 2φˆ C , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ xˆ f (xˆ f − xˆ C ) ⎨ xˆ f (1 − xˆ f ) ˆ φ + , = xˆ C (1 − xˆ C ) C 1 − xˆ C ⎪ ⎪ ⎪ ⎪ 1, ⎪ ⎪ ⎪ ⎪ ⎩ˆ φC , ⎧[ ] xˆ f xˆ f − xˆ C ⎪ ⎪ φˆ C , 1+ ⎪ ⎪ 3(1 − xˆ C ) xˆ C ⎪ ⎪ ⎪ xˆ (1 − xˆ ) ⎪ ˆ ˆ f − xˆ C ) x ( f x f ⎨ f φˆ + , = xˆ C (1 − xˆ C ) C 1 − xˆ C ⎪ ⎪ 1 − xˆ f ⎪ ⎪1 − (1 − φˆ C ), ⎪ ⎪ ⎪ 2(1 − xˆ C ) ⎪ ⎩ˆ

φC ,

xˆ C 0 < φˆ C < 3 xˆ C xˆ C 3 xˆ C xˆ f

≤ φˆ C ≤

xˆ f

(1 + xˆ f − xˆ C )

(1 + xˆ f − xˆ C ) < φˆ C < 1

elsewhere xˆ C xˆ f (xˆ f − xˆ C )

0 < φˆ C <

2xˆ C (1 − xˆ C ) − xˆ f (1 − xˆ f ) xˆ C xˆ f (xˆ f − xˆ C ) xˆ C

≤ φˆ C ≤

2xˆ C (1 − xˆ C ) − xˆ f (1 − xˆ f ) xˆ C (1 + xˆ f − xˆ C ) < φˆ C < 1 xˆ f

xˆ f

(1 + xˆ f − xˆ C )

elsewhere 0 < φˆ C <

3 4

xˆ C

1 + 2(xˆ f − xˆ C ) xˆ C xˆ C ≤ φˆ C ≤ 4 2xˆ f − xˆ C 1 + 2(xˆ f − xˆ C ) xˆ C < φˆ C < 1 2xˆ f − xˆ C 3

elsewhere

Fig. 2. Plots of MINMOD, GAMMA, SMART, WACEB and CUBISTA schemes in the NVD.

Please cite this article in press as: T. Chourushi, Effect of fluid elasticity on the numerical stability of high-resolution schemes for high shearing contraction flows using OpenFOAM, Theoretical & Applied Mechanics Letters (2017), http://dx.doi.org/10.1016/j.taml.2017.01.005.

T. Chourushi / Theoretical & Applied Mechanics Letters (

)

–

5

By substituting Eq. (26) in Eq. (25), generates following pressure correction equation:

( ∇·

1 AC

) ∇p

( =∇· =

H (U )

)

[ AC ∑ ( H (U ) ) f

AC

] · Sf .

(27)

f

Once, pressure field is evaluated, the face flux for control volume is assembled to check any continuity errors. In this paper, SIMPLE algorithm [33] is adopted to obtain the solution of these equations. Steps of the algorithm are outlined below (for given initial conditions of velocity U, pressure p and stress τ fields): Fig. 3. Sketch of two-dimensional 4:1 planar contraction, wherein vertical and horizontal lines are marked with a and b letters, respectively.

Once the velocity field is calculated, its flux (F ) through the control volume face is assembled, as under: F =∇ ·U =

∑

Uf · Sf = 0.

(25)

f

Here, Sf refers to the area vector of cell-face and Uf refers to the velocity field of cell-face. These fields are linearly interpolated from variable values at control volume grid points.

[ Uf =

H (U ) AC

]

( − f

1 AC

)

(∇ p)f . f

(26)

1. Firstly, momentum equation is solved along each coordinate directions with the calculated or guessed values of pressure gradient at time, t = 0. 2. Secondly, H (U ) operator is evaluated using the predicted velocity U ∗ field to obtain the solution of pressure correction equation. (referred in Eq. (27)) 3. This new obtained pressure field p∗ is then substituted in the momentum equation to obtain corrected velocity U ∗∗ field. 4. Lastly, constitutive equations are solved using the corrected velocity field, to obtain new stress τ ∗ field. (refer to Eq. (7)) 5. Above steps are repeated, until all equations reach convergence. In these equations, asymmetric matrices are solved using the preconditioned conjugate gradient (PCG) solver and symmetric

Fig. 4. Residual plot of τxx equation required for the convergence of HRS in case of viscoelastic flow in a 4:1 planar contraction, at (a) E ≈ 1.8; (b) E ≈ 3.5; (c) E ≈ 5; (d) E ≈ 7.

Please cite this article in press as: T. Chourushi, Effect of fluid elasticity on the numerical stability of high-resolution schemes for high shearing contraction flows using OpenFOAM, Theoretical & Applied Mechanics Letters (2017), http://dx.doi.org/10.1016/j.taml.2017.01.005.

6

T. Chourushi / Theoretical & Applied Mechanics Letters (

Fig. 5. Plots of Ux velocity profiles measured along (a) vertical and (b) horizontal lines in the upstream section of 4:1 planar contraction for E ≈ 5.

matrices are solved using the preconditioned bi-conjugate gradient (PBiCG) solver [34]. 3. Flow geometry and mesh characteristics This section details the flow geometry and mesh characteristics for planar contraction geometry. Figure 3 shows the sketch of two-dimensional 4:1 planar contraction geometry. In this figure, both upstream and downstream lengths are taken sufficiently long to obtain fully developed flow. All lines in the upstream section of domain are numbered using positive subscript and lines in the downstream section of domain are numbered using negative subscript, while lines at the corner of domain are assigned with zero subscript. This geometry though looks simple, is numerically stringent to study, as it generates high stresses near the point of contraction. In previous work of Alves et al. [10], it was found that the numerical stability of HRS tends to decrease with the increase in mesh refinement for varying Deborah number. This work therefore considers the mesh characteristics similar to Jovani et al. [15], on account of its cost effectiveness and numerical accuracy. Accordingly, the domain is decomposed into hexahedral meshes of 20, 700 control volumes with smallest cell (1xmin /h = 0.0065×1ymin /h = 0.017) at the contraction region. All numerical simulations are performed using the volumetric flow-rate of Q = 252 cm3 ·s−1 , where shear-rate reaches up to the order of 48.4 s−1 .

)

–

Fig. 6. Plots of Ux velocity profiles measured along (a) vertical and (b) horizontal lines in the downstream section of 4:1 planar contraction for E ≈ 5.

4. Results and discussion Numerical work under this section is broadly classified as: (1) Comparison of flow variable profiles for HRS; (2) Evaluation of numerical stability of HRS for different elasticity number. 4.1. Comparison of flow variable profiles for high resolution schemes In this section, a comparative study for flow variable profiles of HRS is presented. All boundary conditions of the domain, are set in accordance with the flow characteristics. At the wall boundary, no-slip conditions are imposed for velocity field and zero-gradient conditions are applied for pressure and stress fields. And, at the inlet boundary, dirichlet condition is applied for velocity field which is set in reference to flow-rate and zero-gradient conditions are employed for pressure and stress fields. Also, at the outlet boundary, zero-gradient conditions are specified for velocity and stress fields, while pressure field is set to zero value. On the other hand, at the symmetric boundary, symmetric conditions are applied for velocity, pressure and stress fields. In all runs, underrelaxation factors for velocity, pressure and stress fields are considered as 0.5, 0.3 and 0.3, respectively. In order to attain numerical stability of time-dependent flows, run-time modifiable time-step

Please cite this article in press as: T. Chourushi, Effect of fluid elasticity on the numerical stability of high-resolution schemes for high shearing contraction flows using OpenFOAM, Theoretical & Applied Mechanics Letters (2017), http://dx.doi.org/10.1016/j.taml.2017.01.005.

T. Chourushi / Theoretical & Applied Mechanics Letters (

a

)

–

7

b

Fig. 7. Plots of shear-stress (τxy ) profiles measured along vertical and horizontal lines in the (a) upstream and (b) downstream sections of 4:1 planar contraction for E ≈ 5.

in reference to maximum Courant Friedrichs Lewy (CFL) number of 0.8 is considered. Flow parameters for single-mode Giesekus fluid are taken as, ρ = 803.87 kg·m−3 , α = 0.15, λ = 0.03 s, ηP = 1.422 Pa · s and ηS = 0.002 Pa · s. Thereby, retardation ratio and elasticity number of fluid are β = ηS /(ηS + ηP ) ≈ 0.001 and E = λ(ηS + ηP )/(ρ h2 ) = 5.18 ≈ 5, respectively. All simulations are run for sufficiently long time, until initial residual of stress equation falls below 10−6 . Fig. 4(c) shows residual plot of τxx equation required for the convergence of HRS, at E ≈ 5. As depicted in the figure, first-order UDS converges first, which is followed by MINMOD, GAMMA, and CUBISTA schemes. On the other hand, WACEB and SMART schemes failed to reach convergence, at this elasticity number. To further assess the numerical stability of HRS, flow variable profiles along vertical and horizontal lines (refer to Fig. 3) of the domain are detailed here. Figure 5(a) depicts the plot of Ux velocity profile measured along vertical lines in the upstream section of 4:1 planar contraction. In this figure, as the fluid moves toward contraction region, it gets accelerated in the centre while it decelerates toward the walls. Figure 5(b) also satisfy this behaviour, as outlined in plot of Ux velocity profiles measured along horizontal lines in the upstream section of domain. However, Ux velocity profiles measured along vertical and horizontal lines in the downstream section of domain show parabolic profiles, as detailed in Fig. 6. These results are in agreement with previous literatures [15,35]. Moreover, it has been found that GAMMA and CUBISTA show

minute undue oscillations in velocity profiles measured along vertical line placed nearer to contraction region in the upstream section of domain. As, shear stress (τxy ) relies on the velocity field, its profiles get affected. Figure 7 shows the plot of shear stress profiles measured along vertical and horizontal lines in the upstream and downstream sections of domain. In Fig. 7(a), shear stress value increases to maximum along vertical lines placed nearer to the re-entrant corner and thereafter it reduces to a constant value. While, these profiles constantly increase along vertical lines in the downstream section of domain, as outlined in Fig. 7(b). Broadly, HRS show some order of undue oscillations in flow variable profiles along vertical lines placed near the contraction region. And profiles along horizontal lines in the upstream and downstream sections of domain resemble this behaviour. Comparable results are obtained for first normal stress difference (τxx − τyy ) profiles. Figure 8 show its plot along vertical and horizontal lines in the upstream and downstream sections of domain. These results suggest that all HRS generate undue oscillations in flow variable profiles along vertical lines placed near contraction region in the upstream section of domain. To further evaluate the actual cause of this behaviour, numerical studies using different elasticity numbers [36] are presented next.

Please cite this article in press as: T. Chourushi, Effect of fluid elasticity on the numerical stability of high-resolution schemes for high shearing contraction flows using OpenFOAM, Theoretical & Applied Mechanics Letters (2017), http://dx.doi.org/10.1016/j.taml.2017.01.005.

8

T. Chourushi / Theoretical & Applied Mechanics Letters (

a

)

–

b

Fig. 8. Plots of first normal stress difference τxx − τyy profiles measured along vertical and horizontal lines in the (a) upstream and (b) downstream sections of 4:1 planar contraction for E ≈ 5.

(

)

4.2. Evaluation of numerical stability of HRS for different elasticity number Based on previous results, it is confirmed that all HRS show some order of undue oscillations in flow variable profiles for Re = 0.57 and De = 1.45, thereby β ≈ 0.001 and E ≈ 5. Also, during our numerical experiments it has been found that the order of undue oscillations in flow variable profiles of HRS remained comparatively same for Giesekus fluid, as β ranges from 0.001 to 0.1. Moreover, Rodd et al. [3] suggested that the inertia and elasticity of fluid plays an important role in contraction flows. Hence, we assessed the numerical stability of HRS against the varying elasticity number viz., E ≈ 1.8, 3.5, 5, 7. In this context, E is varied by adjusting only the λ value and all other parameters are considered similar to earlier section of the paper. All numerical simulations are run for sufficiently long period of time and are assumed to be converged when the initial residual of stress equation falls below 10−6 . Figure 4 shows the residual plot of τxx equation required for the convergence of HRS for different elasticity number. As depicted in Fig. 4(a), all HRS converges up to prescribed order of tolerance 10−6 . In broad terms, UDS converges first which is followed by MINMOD, GAMMA, and CUBISTA schemes. On the other hand, both WACEB and SMART schemes took long time for convergence. With further increase in elasticity number, numerical convergence of HRS became unattainable. When E ≈ 7, only UDS and MINMOD schemes reach convergence, as outlined in Fig. 4(d). In particular,

WACEB and SMART schemes which are dependent on the downwind line of NVD show loss of numerical stability at E ≈ 3.5, 5, 7. However, GAMMA and CUBISTA schemes remained numerically stable only up to the order of E ≈ 5. This suggests that only those HRS which are nearer to unconditionally stable UDS in the NVD (such as, MINMOD), possess better numerical stability than other HRS. To demonstrate this behaviour, we detailed the comparison of shear stress (τxy ) and normal stress difference (τxx − τyy ) profiles of HRS for varied elasticity numbers. Figure 9 depicts the plot of shear stress profiles for UDS and HRS measured along vertical line a1 (= − 0.37) in the upstream section of domain using different elasticity number (E). As depicted in the figure, UDS captures smoothly varying profiles for all E values, while MINMOD scheme captures minutely oscillating profiles. These oscillations further gain intensity for GAMMA and CUBISTA schemes. And in case of SMART and WACEB schemes large oscillating profiles are obtained, leading to the failure of numerical convergence at moderate elasticity numbers, that is, E ≥ 3.5. These oscillating profiles are more predominantly observed in first normal stress difference plots. Figure 10 shows the plot of these profiles for UDS and HRS measured along vertical line a1 (=− 0.37) in the upstream section of domain using varied E. In this figure, all HRS report some order of undue oscillations in first normal stress difference profiles and its intensity increase with the rise in E. At E ≥ 5, large undue oscillations are reported for MINMOD, GAMMA and CUBISTA schemes. While, both WACEB and SMART schemes capture large oscillating profiles for lower value of E(≈1.8), which

Please cite this article in press as: T. Chourushi, Effect of fluid elasticity on the numerical stability of high-resolution schemes for high shearing contraction flows using OpenFOAM, Theoretical & Applied Mechanics Letters (2017), http://dx.doi.org/10.1016/j.taml.2017.01.005.

T. Chourushi / Theoretical & Applied Mechanics Letters (

)

–

9

Fig. 9. Plots of shear stress τxy profiles measured along vertical line (a1 = −0.37) in the upstream section of 4:1 planar contraction for (a) UDS, (b) MINMOD, (c) GAMMA, (d) SMART, (e) WACEB and (f) CUBISTA schemes, using varied E ≈ 1.8, 3.5, 5, 7.

gives rise to the failure of numerical convergence. In fact, those HRS which are distant from the UDS in NVD, more likely results into loss of numerical stability for the increase in E. 5. Conclusions It is numerically challenging to study the behaviour of viscoelastic fluids, on account of its non-linearity which is dependent

upon the constitutive equation. These non-linear characteristics result into loss of numerical stability for HRS, at high shearing rates. To further enlighten these effects, the paper presented a comparative study of numerical stability of HRS for high shearing planar contraction flows. In the present context, localize normalized variable approach [21,32] is being adopted for implementation of HRS. And the existing viscoelastic fluid flow solver [15] is considered for numerical works. Results suggest that all HRS

Please cite this article in press as: T. Chourushi, Effect of fluid elasticity on the numerical stability of high-resolution schemes for high shearing contraction flows using OpenFOAM, Theoretical & Applied Mechanics Letters (2017), http://dx.doi.org/10.1016/j.taml.2017.01.005.

10

T. Chourushi / Theoretical & Applied Mechanics Letters (

a

b

c

d

e

f

)

–

Fig. 10. Plots of first normal stress difference τxx − τyy profiles measured along vertical line (a1 = −0.37) in the upstream section of 4:1 planar contraction for (a) UDS, (b) MINMOD, (c) GAMMA, (d) SMART, (e) WACEB and (f) CUBISTA schemes, using varied E ≈ 1.8, 3.5, 5, 7.

(

)

show some order of undue oscillations in flow variable profiles measured along vertical lines placed nearer to contraction region in the upstream section of domain. In general, order of undue oscillations in flow variable profiles of HRS is directly proportional to the elasticity number (E) of fluid and its intensity increases with

the rise in E. Furthermore, those HRS which are distant from the UDS in the NVD, loses numerical stability prior to other schemes. In this study, effect of fluid elasticity based on the numerical stability of HRS is addressed for different elasticity number E in reference to 4:1 planar contraction flows. In future, it is intended to

Please cite this article in press as: T. Chourushi, Effect of fluid elasticity on the numerical stability of high-resolution schemes for high shearing contraction flows using OpenFOAM, Theoretical & Applied Mechanics Letters (2017), http://dx.doi.org/10.1016/j.taml.2017.01.005.

T. Chourushi / Theoretical & Applied Mechanics Letters (

extend this work for different contraction ratios, viscoelastic fluid models and stabilization methods for viscoelastic fluids. Acknowledgements The author is grateful to the Indian Institute of Technology Indore (IITI) for providing with the research facility. Also, the author is thankful to Dr. S. Dhinakaran, Dr. D. Deshmukh, and Dr. S.S. Ahmad for their valuable support. References [1] X. Xingming, Z. Guoqun, Q. Shengxue, et al., Numerical simulation of viscoelastic extrudate swell through elliptical ring die, Chin. J. Chem. Eng. 19 (2011) 10–17. [2] S.A. White, A.D. Gotsis, D.G. Baird, Review of the entry flow problem: experimental and numerical, J. Non-Newton. Fluid Mech. 24 (1987) 121–160. [3] L.E. Rodd, T.P. Scott, D.V. Boger, et al., The inertio-elastic planar entry flow of low-viscosity elastic fluids in micro-fabricated geometries, J. Non-Newton. Fluid Mech. 129 (2005) 1–22. [4] K. Kang, L. James Lee, K.W. Koelling, High shear microfluidics and its application in rheological measurement, Exp. Fluids 38 (2005) 222–232. [5] P.J. Oliveira, F.T. Pinho, Plane contraction flows of upper convected Maxwell and Phan-Thien-Tanner fluids as predicted by a finite-volume method, J. NonNewton. Fluid Mech. 88 (1999) 63–88. [6] A.M. Afonso, P.J. Oliveira, F.T. Pinho, et al., Dynamics of high-Deborah-number entry flows: a numerical study, J. Fluid Mech. 677 (2011) 272–304. [7] P.J. Oliveira, F.T. Pinho, G.A. Pinto, Numerical simulation of non-linear elastic flows with a general collocated finite-volume method, J. Non-Newton. Fluid Mech. 79 (1998) 1–43. [8] J.H. Song, J.Y. Yoo, Numerical simulation of viscoelastic flow through a sudden contraction using a type dependent difference method, J. Non-Newton. Fluid Mech. 24 (1987) 221–243. [9] J.M. Marchal, M.J. Crochet, A new mixed finite element for calculating viscoelastic flow, J. Non-Newton. Fluid Mech. 26 (1987) 77–114. [10] M.A. Alves, F.T. Pinho, P.J. Oliveira, Effect of a high-resolution differencing scheme on finite-volume predictions of viscoelastic flows, J. Non-Newton. Fluid Mech. 93 (2000) 287–314. [11] L.J. Amoreira, P.J. Oliveira, Comparison of different formulations for the numerical calculation of unsteady incompressible viscoelastic fluid flow, Adv. Appl. Math. Mech. 2 (2010) 483–502. [12] R. Fattal, R. Kupferman, Time-dependent simulation of viscoelastic flows at high weissenberg number using the Log-Conformation Representation, J. NonNewton. Fluid Mech. 126 (2005) 23–37. [13] P.A. Stewart, N. Lay, M. Sussman, et al., An improved sharp interface method for viscoelastic and viscous two-phase flows, J. Sci. Comput. 35 (2008) 43–61. [14] N. Balci, B. Thomases, M. Renardy, et al., Symmetric Factorization of the Conformation Tensor in Viscoelastic Fluid Models, J. Non-Newton. Fluid Mech. 166 (2011) 546–553. [15] J.L. Favero, A.R. Secchi, N.S.M. Cardozo, et al., Viscoelastic flow analysis using the software OpenFOAM and differential constitutive equations, J. Non-Newton. Fluid Mech. 165 (2010) 1625–1636. [16] X. Chen, H. Marschall, M. Schäfer, et al., A comparison of stabilisation approaches for finite-volume simulation of viscoelastic fluid flow, Int. J. Comput. Fluid Dyn. 27 (2013) 229–250. [17] H.K. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, second ed., Pearson Education Limited, 2007. [18] P.K. Sweby, High resolution schemes using flux-limiters for hyperbolic conservation laws, Soc. Ind. Appl. Math. J. Numer. Anal. 21 (1984) 995–1011.

)

–

11

[19] P.H. Gaskell, A.K.C. Lau, Curvature-compensated convective transport: SMART, A new boundedness perserving transport algorithm, Internat. J. Numer. Methods Fluids 8 (1998) 617–641. [20] G.D. Van Albada, B. Van Leer, W.W. Roberts, A comparative study of computational methods in cosmic gas dynamics, Astron. Astrophys. 108 (1982) 76–84. [21] H. Jasak, H.G. Weller, A.D. Gosman, High Resolution NVD differencing scheme for arbitrarily unstructured meshes, Internat. J. Numer. Methods Fluids 31 (1999) 431–449. [22] B. Song, G.R. Liu, K.Y. Lam, et al., On a higher-order bounded discretization scheme, Internat. J. Numer. Methods Fluids 32 (2000) 881–897. [23] M.A. Alves, P.J. Oliveira, F.T. Pinho, A convergent and universally bounded interpolation scheme for the treatment of advection, Internat. J. Numer. Methods Fluids 41 (2003) 47–75. [24] M.A. Alves, P.J. Oliveira, F.T. Pinho, Benchmark solutions for the flow of OldroydB and PTT fluids in planar contractions, J. Non-Newton. Fluid Mech. 110 (2003) 45–75. [25] S.C. Xue, N. Phan-Thien, R.I. Tanner, Three dimensional numerical simulations of viscoelastic flows through planar contractions, J. Non-Newton. Fluid Mech. 74 (1998) 195–245. [26] F. Habla, M.W. Tan, W.J. Haßlberger, et al., Numerical simulation of the viscoelastic flow in a three-dimensional lid-driven cavity using the logconformation reformulation in OpenFOAM, J. Non-Newton. Fluid Mech. 212 (2014) 47–62. [27] H. Giesekus, A simple constitutive equation for polymer fluids based on the concept of deformation-dependent tensorial mobility, J. Non-Newton. Fluid Mech. 11 (1982) 69–109. [28] R.B. Bird, R.C. Armstorng, O. Hassager, Dyn. of pol. liq, in: Fluid Dyn, second ed., Wiley, New York, 1987, p. 1. [29] A.N. Beris, R.C. Armstrong, R.A. Brown, Finite element calculation of viscoelastic flow in a journal bearing: II. Moderate eccentricity, J. Non-Newton. Fluid Mech. 19 (1986) 323–347. [30] W.M. Christopher, Rheology: Principles Measurements and Applications, Wiley, New York, 1994. [31] C.M. Rhie, W.L. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA J. 21 (1983) 1525–1532. [32] H. Jasak, Error Analysis and Estimation in the Finite Volume Method with Applications to Fluid Flows (Ph.D. thesis), Imperial College, University of London, 1996. [33] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, in: Hemisphere Series on Comput. Meth. in Mech. and Therm. Sci., Washington DC, 1980. [34] H.G. Weller, G. Tabor, H. Jasak, et al., A tensorial approach to computational continuum mechanics using object-oriented techniques, J. Comput. Phys. 12 (1998) 620–631. [35] L.M. Quinzani, R.C. Armstrong, R.A. Brown, Birefringence and laser Doppler velocimetry (ldv) studies of viscoelastic flow through a planar contraction, J. Non-Newton. Fluid Mech. 52 (1994) 1–36. [36] J.M. Malheiro, P.J. Oliveira, F.T. Pinho, Planar time-dependent viscoelastic flow - pulsating flow and effect of elasticity number, in: III Conferência Nacional em Mecânica de Fluidos, Termodinâmica e Energia, MEFTE - BRAGANÇA 09, 2009, pp. 1–9.

Further Reading [1] F. Habla, A. Obermeier, O. Hinrichsen, Semi-implicit stress formulation for viscoelastic models: Application to three-dimensional contraction flows, J. Non-Newton. Fluid Mech. 199 (2013) 70–79.

Please cite this article in press as: T. Chourushi, Effect of fluid elasticity on the numerical stability of high-resolution schemes for high shearing contraction flows using OpenFOAM, Theoretical & Applied Mechanics Letters (2017), http://dx.doi.org/10.1016/j.taml.2017.01.005.