Viscoelastic Flow through an Axisymmetric Contraction Using the Grid-by-Grid Inversion Method

Abstract

The newly developed algorithm called the grid-by-grid inversion method is a very convenient method for converting an existing computer code for Newtonian flow simulations to that for viscoelastic flow simulations. In this method, the hyperbolic constitutive equation is split such that the term for the convective transport of stress tensor is treated as a source which is updated iteratively. This allows the stress tensors at each grid point to be expressed in terms of velocity gradient tensor at the same location, and the set of stress tensor components is found after inverting a small matrix at each grid point. To corroborate the robustness and accuracy of the grid-by-grid inversion method, we apply it to the 4 : 1 axisymmetric contraction problem. This algorithm is found to be robust and yields accurate results as compared with other finite volume methods. Any commercial CFD packages for Newtonian flow simulations can be easily converted to those for viscoelastic fluids exploiting the grid-by-grid inversion method.

1. Introduction

Contrary to the techniques of computational fluid dynamics for Newtonian fluids, the numerical algorithms for viscoelastic flows are not so matured. The hyperbolic nature of the constitutive equation incurs peculiar flow phenomena such as rod-climbing and extrudate swell as well as causes difficulties in numerical simulation [1–3]. Since the momentum balance equation is elliptic in steady state and parabolic in unsteady state, the complete set for viscoelastic flows is a mixed type, hyperbolic-elliptic, or hyperbolic-parabolic. This situation is difficult to treat numerically since it is difficult to devise a numerical algorithm that works for mixed systems. Another difficulty associated with the hyperbolic constitutive equation is the choice of appropriate boundary conditions for the stress field at the boundary of computational domain. One can impose nonslip boundary condition on the walls for the velocity field but there exist no apparent or natural boundary conditions for the stress components at the wall. Various numerical techniques for solving viscoelastic flows, such as finite volume methods, finite element methods, and spectral methods, are well documented in the references cited [4–9].

In the present investigation, a newly developed algorithm called the grid-by-grid inversion method [10] is employed to solve the viscoelastic flows through an axisymmetric contraction. Figure 1 shows the flow geometry for the axisymmetric contraction. Viscoelastic fluid passes from one circular tube into the other tube of smaller radius and generates a complex flow having a strong shear near the walls and uniaxial extension along the centerline. The existence of strong shearing zones and uniaxial extension zone makes this flow geometry a good test bed for numerical algorithms of viscoelastic fluids. Many investigators of rheology adopt this flow geometry as an important benchmark problem, especially the 4 : 1 contraction geometry [11]. We shall solve the 4 : 1 axisymmetric contraction problem employing the newly developed “grid-by-grid inversion method” [10], which is implemented based on a finite volume method [10].

Figure 1 

The 4 : 1 axisymmetric contraction.

The hyperbolic constitutive equations of viscoelastic fluids have a nonlocal character because of the term representing convective transport of stress tensor. If this term is assumed to be known, the constitutive equation becomes local and the stress tensor is easily evaluated for a given velocity gradient tensor at the same location. The six stress tensor components for the cases of a three-dimensional flow are found after inverting a six by six matrix at each grid point and are substituted into the Navier-Stokes equation as a source term. In this way, the numerical solution of viscoelastic flows becomes as straightforward as that of Newtonian fluids. We call this algorithm the grid-by-grid inversion method since the viscoelastic stress tensor is obtained by the grid-by-grid inversion of a matrix equation at each grid point. This algorithm can easily be implemented using finite volume methods, finite element methods and spectral methods. When applied to the 4 : 1 axisymmetric contraction problem, it is found that the grid-by-grid inversion method yields accurate results efficiently in comparison with numerical results of traditional algorithms.

2. Governing Equations and the Grid-by-Grid Inversion Method

We consider incompressible isothermal flows of viscoelastic fluids. The governing equations may be written in dimensionless variables as follows:Re𝜕𝐯∇⋅𝐯=0,(2.1)𝜆𝜕𝝉𝝈𝜕=𝑡−+𝑃𝐯𝐈⋅+∇2𝐯𝛽𝐃=+∇𝝉,⋅𝝈,((22..23))𝜕𝑡+𝐯⋅∇𝝉=2(1−𝛽)𝐃+𝜆(∇𝐯)^𝑇⋅𝝉+𝝉⋅∇𝐯−𝝉.(2.4) In the above equations, 𝑃 is pressure, 𝐃 is the rate of deformation tensor, 𝝉 is the viscoelastic part of the stress tensor. The parameter 𝛽 is the ratio of the retardation and relaxation times of the fluid. Equation (2.4) is the constitutive equation of the Oldroyd-B model [9], 𝜆 is the dimensionless relaxation time or the Deborah number, and Re is the Reynolds number. The superscript 𝑇 in (2.4) is the transpose. The Reynolds number and the dimensionless Deborah number are defined byRe=𝜌𝑈𝐿𝜂,𝜆=𝜆¹𝑈𝐿,(2.5) where 𝜌 is the density, 𝜆₁ the dimensional relaxation time, 𝑈 the characteristic speed, and 𝐿 is the characteristic length. The characteristic velocity 𝑈 and the characteristic length 𝐿 are taken as the average velocity in the downstream tube and the radius of the downstream tube, respectively. To compare with the results of other investigators [9, 11], we take the value of 𝛽 to be 1/9.

Next, we consider the grid-by-grid inversion method as applied to the set of equations (2.1)–(2.4). Contrary to the Newtonian fluids where the stress field 𝝉 depends on the velocity gradient tensor ∇𝐯 locally, the stress depends on the velocity gradient tensor nonlocally due to the fact that constitutive equation is hyperbolic partial differential equation. The convective transport of the stress tensor, represented by 𝐯⋅∇𝝉 in (2.4), takes care of the memory effect of 𝝉 in its dependence on ∇𝐯 and causes the functional dependence of 𝝉 on ∇𝐯 to be nonlocal. In the grid-by-grid inversion method, 𝐯⋅∇𝝉 is assumed to be a known source term which shall be updated iteratively in the process. Then the constitutive equation (2.4) can be converted such that 𝝉 is represented as a local function of ∇𝐯 as in the case of Newtonian fluids. This arrangement renders the numerical solution of viscoelastic fluids as straightforward as that of Newtonian fluids. The steps for obtaining the velocity and the stress at time step 𝑛+1, 𝐯^𝑛+1 and 𝝉^𝑛+1, using known velocity and stress at time step 𝑛,𝐯^𝑛 and 𝝉^𝑛, proceed in an iterative manner as follows. First, define𝐂≡𝐯^{𝑛+1(𝑖𝑡)}⋅∇𝝉^{𝑛+1(𝑖𝑡)},(2.6) where the superscript 𝑛+1(𝑖𝑡) indicates variable at time step 𝑛+1 in the 𝑖𝑡th iteration. Discretizing (2.4) in time implicitly and representing 𝐃 in terms of the velocity gradient ∇𝐯, we find the following local matrix equations defined at each grid point.Δ𝜆𝑡+1.0𝝉^{𝑛+1(𝑖𝑡+1)}−𝜆(∇𝐯)^{𝑇𝑛+1(𝑖𝑡)}⋅𝝉^{𝑛+1(𝑖𝑡+1)}−𝜆𝝉^{𝑛+1(𝑖𝑡+1)}⋅(∇𝐯)^{𝑛+1(𝑖𝑡)}=𝜆Δ𝑡𝝉^𝑛−𝜆𝐂+(1−𝛽)(∇𝐯)^{𝑛+1(𝑖𝑡)}+(∇𝐯)^{𝑇𝑛+1(𝑖𝑡)}.(2.7) Once 𝐂 is evaluated using variables at 𝑛+1(𝑖𝑡), (2.7) can be solved for 𝝉^{𝑛+1(𝑖𝑡+1)} by inverting a six-by-six matrix at each grid point for the six independent stress components for the case of three-dimensional flows. Equation (2.7) is also solved at the boundary grid points to find the boundary stress field. This suggests a natural method of imposing boundary conditions for the hyperbolic constitutive equations at all boundaries of the computational domain. This method had been employed to impose outflow boundary conditions by Papanastasiou et al. [12] and is called the open boundary condition. Using 𝝉^{𝑛+1(𝑖𝑡+1)} obtained from (2.7), the velocity at the 𝑛+1 th time step in the 𝑖𝑡+1th iteration, 𝐯^{𝑛+1(𝑖𝑡+1)}, is found by solving (2.1)–(2.3) as follows:Re𝐯_{𝑛+1(𝑖𝑡+1)}∇⋅𝐯^{𝑛+1(𝑖𝑡+1)}=0,(2.8)−𝐯^𝑛Δ𝑡+𝐯^{𝑛+1(𝑖𝑡+1)}⋅∇𝐯^{𝑛+1(𝑖𝑡+1)}=−∇𝑃^{𝑛+1(𝑖𝑡+1)}+2𝛽∇²𝐯^{𝑛+1(𝑖𝑡+1)}+∇⋅𝝉^{𝑛+1(𝑖𝑡+1)}.(2.9) Although (2.8)-(2.9) can be solved using any numerical methods for incompressible Navier-Stokes equations, we employ a finite volume method based on the SIMPLE algorithm [13, 14]. To stabilize the numerical scheme for large values of 𝜆, (2.6) is evaluated using one of the various upwind schemes. In the present work, we adopt a higher-order upwind scheme, MinMod [14], to evaluate 𝐂.

The case of 𝛽=0 requires a special deliberation before employing standard algorithms for incompressible Navier-Stokes equations in the grid-by-grid inversion technique. When 𝛽=0, we decompose the total stress 𝝈 into pure viscous and elastic parts −𝑃𝐈+2𝐃 and 𝚺, respectively, such that𝝈=−𝑃𝐈+2𝐃+𝚺,(2.10) where𝚺≡𝝉−2𝐃.(2.11) Then the momentum and constitutive equations for the case of 𝛽=0 may be written asRe𝜕𝐯∇⋅𝐯=0,(2.12)𝜕𝑡+𝐯⋅∇𝐯=−∇𝑃+∇²𝜆𝜕𝚺𝐯+∇⋅𝚺,(2.13)𝜕𝑡+𝐯⋅∇𝚺−(∇𝐯)^𝑇⋅𝚺−𝚺⋅(∇𝐯)+𝚺=−2𝜆𝜕𝐃𝜕𝑡+𝐯⋅∇𝐃−(∇𝐯)^𝑇⋅𝐃−𝐃⋅∇𝐯.(2.14) This decomposition is called the EVSS formulation and was proposed by Perera and Walters [15]. The implementation of grid-by-grid inversion method to the set of (2.12)–(2.14) is straightforward. First, we evaluate𝐄≡𝐯^{𝑛+1(𝑖𝑡)}⋅∇𝚺^{𝑛+1(𝑖𝑡)}.(2.15) Next, convert (2.14) to the following local matrix equations for 𝚺 decoupled at each grid point:Δ𝜆𝑡+1.0𝚺^{𝑛+1(𝑖𝑡+1)}−𝜆(∇𝐯)^{𝑇𝑛+1(𝑖𝑡)}⋅𝚺^{𝑛+1(𝑖𝑡+1)}−𝜆𝚺^{𝑛+1(𝑖𝑡+1)}⋅(∇𝐯)^{𝑛+1(𝑖𝑡)}=𝜆Δ𝑡𝚺^𝑛−𝜆𝐄−2𝜆1Δ𝑡𝐃^{𝑛+1(𝑖𝑡)}−1Δ𝑡𝐃^𝑛+𝐯⋅∇𝐃−(∇𝐯)^𝑇⋅𝐃−𝐃⋅∇𝐯^{𝑛+1(𝑖𝑡)}.(2.16) Then, (2.12) and (2.13) are discretized as follows: Re𝐯_{𝑛+1(𝑖𝑡+1)}∇⋅𝐯^{𝑛+1(𝑖𝑡+1)}(2.17)−𝐯^𝑛Δ𝑡+𝐯^{𝑛+1(𝑖𝑡+1)}⋅∇𝐯^{𝑛+1(𝑖𝑡+1)}=−𝑃^{𝑛+1(𝑖𝑡+1)}+∇²𝐯^{𝑛+1(𝑖𝑡+1)}+∇⋅𝚺^{𝑛+1(𝑖𝑡+1)}.(2.18) The structure of (2.16) is the same as that of (2.7) and can be solved for 𝚺^{𝑛+1(𝑖𝑡+1)} by inverting a six-by-six matrix at each grid point for a three-dimensional flow once (∇𝐯)^{𝑛+1(𝑖𝑡)} is given. After obtaining 𝚺^{𝑛+1(𝑖𝑡+1)}, 𝐯^{𝑛+1(𝑖𝑡+1)} is found by solving (2.17)-(2.18) using the SIMPLE algorithm [13]. The second-order derivative terms appearing in 𝐯⋅∇𝐃 of (2.14) are evaluated using a one-sided finite difference formula at the boundaries.

The overall solution procedure for the grid-by-grid inversion method may be summarized as follows.(1)𝐯^𝑛 and 𝝉^𝑛 have obtained in the previous time step 𝑛.

The procedure for the time step 𝑛+1 begins as follows.(2) Assume 𝐯^{𝑛+1(𝑖𝑡)} and 𝝉^{𝑛+1(𝑖𝑡)}. For the first iteration (𝑖𝑡=1), 𝐯^{𝑛+1(𝑖𝑡)}=𝐯^𝑛 and 𝝉^{𝑛+1(𝑖𝑡)}=𝝉^𝑛.(3) Evaluate 𝐂 of (2.6) (𝛽≠0) or 𝐄 of (2.15) (𝛽=0) using an upwind scheme.(4) Using 𝐯^{𝑛+1(𝑖𝑡)}, solve (2.7) (𝛽≠0) or (2.16) (𝛽=0) for 𝝉^{𝑛+1(𝑖𝑡+1)} by inverting a six-by-six matrix at each grid point including the boundary grids.(5) Using 𝝉^{𝑛+1(𝑖𝑡+1)},

(𝛽≠0) solve the momentum equation, (2.9), and continuity equation to find 𝐯^{𝑛+1(𝑖𝑡+1)}.

(𝛽=0) solve the momentum equation, (2.18), and continuity equation to find 𝐯^{𝑛+1(𝑖𝑡+1)}.(6) Convergence check for 𝝉^𝑛+1 and 𝐯^𝑛+1. If not converged, go to step (2). Otherwise, update the time step and go to step (1).

Usually convergence is attained in two or three iterations. The novelty of the grid-by-grid inversion method is its easiness of numerical implementation. As noted in the above procedure, one can easily convert any existing computer code for Newtonian flow simulations to that for viscoelastic flow simulation by adding a subroutine that solves the viscoelastic constitutive equation using the grid-by-grid inversion method, summarized as steps (3)~(4), to evaluate the viscoelastic source terms in the Navier-Stokes equation, ∇⋅𝝉. Adding a source term to the Navier-Stokes equation is an easy procedure whether we employ a finite volume method or a finite element method. In the subroutine for the viscoelastic constitutive equation, one inverts a small matrix at each grid point, which can be performed cheaply. Therefore, one can easily convert any commercial CFD package for Newtonian fluid flows to that for viscoelastic fluids employing the grid-by-grid inversion method. Its robustness and accuracy are corroborated in the next section, where the grid-by-grid inversion method is employed to solve the 4 : 1 viscoelastic axisymmetric contraction problem. Although Phillips and Williams [16] solve the viscoelastic constitutive equation by converting small matrix equations in their semi-Lagrange method, it is very difficult to convert an existing Newtonian code to a viscoelastic code using the semi-Lagrange method since it treats the momentum and the constitutive equation simultaneously.

3. Viscoelastic Flow through an Axisymmetric Contraction

In this section, we solve the flow of an Oldroyd-B fluid through a 4 : 1 axisymmetric contraction geometry using the grid-by-grid inversion method. The schematic representation of the flow geometry is depicted in Figure 1. The viscoelastic flow through the circular contraction geometry has served as a standard benchmark problem for numerical algorithms for viscoelastic fluids. The presence of a singularity in the entry-flow geometry as well as regions of high shear and extension in the flow have been a significant challenge for a long time, and many investigators have struggled to develop accurate and robust algorithms for viscoelastic flows through contraction geometries [17, 18]. The contraction ratio 4 : 1 is the standard one in this benchmark problem. A good review of old works may be Trebotich et al. [18]. Recently, Phillips and Williams [11, 16] used a finite volume method to solve the flow of an Oldroyd-B fluid through planar and axisymmetric contractions. Among these two geometries, the axisymmetric contraction yields more dramatic viscoelastic effects. In the present investigation, we will solve the flow of Oldroyd-B fluid through the axisymmetric geometry using the grid-by-grid inversion method based on the SIMPLE algorithm to corroborate the accuracy and robustness of the grid-by-grid inversion method. When solving (2.8)-(2.9) or (2.17)-(2.18) using the SIMPLE algorithm, we adopt a collocated mesh and the pressure oscillation induced by the collocated mesh is eliminated using the Rhie-Chow method [19], and a higher-order upwind scheme, MinMod [14], is used to treat the inertia force term in the momentum equation and the convective transport of stress tensor in the constitutive equation. For the 4 : 1 axisymmetric contraction geometry depicted in Figure 1, the number of relevant components of stress tensor is four because the ∇⋅𝝉 term in (2.9) contains a circular stress component 𝜏^𝜃𝜃 as follows:(∇⋅𝝉)_𝑟=1𝑟𝜕𝜕𝑟(𝑟𝜏^𝑟𝑟)+𝜕𝜕𝑧𝜏^𝑟𝑧−𝜏^𝜃𝜃(∇⋅𝝉)_𝑧𝑟,=1𝑟𝜕𝜕𝑟(𝑟𝜏^𝑟𝑧)+𝜕𝜕𝑧𝜏^𝑧𝑧.(3.1) Although the relevant velocity components for axisymmetric flows are (𝑣^𝑟,𝑣^𝑧), it is necessary to find the four stress components, 𝜏^𝑟𝑟, 𝜏^𝑟𝑧, 𝜏^𝑧𝑧, and 𝜏^𝜃𝜃 which appear in (3.1). The constitutive equation (2.4) may be written for an axisymmetric geometry as follows:𝜆𝜕𝜏^𝑟𝑟𝜕𝑡+𝑣^𝑟𝜕𝜏^𝑟𝑟𝜕𝑟+𝑣^𝑧𝜕𝜏^𝑟𝑟𝜕𝑧=2(1−𝛽)𝜕𝜕𝑣𝑟^𝑟+𝜆2𝜕𝜕𝑣𝑟^𝑟𝜏^𝑟𝑟+2𝜕𝑣^𝑟𝜕𝑧𝜏^𝑟𝑧−𝜏^𝑟𝑟𝜆𝜕𝜏_𝑟𝑧,(3.2)𝜕𝑡+𝑣^𝑟𝜕𝜏^𝑟𝑧𝜕𝑟+𝑣^𝑧𝜕𝜏^𝑟𝑧=(1−𝛽)𝜕𝑣_𝑧𝜕𝑧𝜕𝑟+𝜕𝜕𝑣𝑧^𝑟+𝜆𝜕𝜕𝑣𝑟^𝑟𝜏^𝑟𝑧+𝜕𝑣^𝑟𝜕𝑧𝜏^𝑧𝑧+𝜏^𝑟𝑟𝜕𝑣^𝑧𝜕𝑟+𝜏^𝑟𝑧𝜕𝑣^𝑧𝜕𝑧−𝜏^𝑟𝑧𝜆𝜕𝜏_𝑧𝑧,(3.3)𝜕𝑡+𝑣^𝑟𝜕𝜏^𝑧𝑧𝜕𝑟+𝑣^𝑧𝜕𝜏^𝑧𝑧𝜕𝑧=2(1−𝛽)𝜕𝜕𝑣𝑧^𝑧+𝜆2𝜕𝜕𝑣𝑟^𝑧𝜏^𝑟𝑧+2𝜕𝑣^𝑧𝜕𝑧𝜏^𝑧𝑧−𝜏^𝑧𝑧𝜆𝜕𝜏_𝜃𝜃,(3.4)𝜕𝑡+𝑣^𝑟𝜕𝜏^𝜃𝜃𝜕𝑟+𝑣^𝑧𝜕𝜏^𝜃𝜃𝜕𝑧=2(1−𝛽)𝑣𝑟^𝑟+𝜆2𝑣𝑟^𝑟𝜏^𝜃𝜃−𝜏^𝜃𝜃.(3.5) It is to be noted that equations for (𝜏^𝑟𝑟,𝜏^𝑟𝑧,𝜏^𝑧𝑧) are coupled with each other, while that for 𝜏^𝜃𝜃 is decoupled. The relevant components of 𝐂 of (2.6) are𝐶^𝐼𝐽=𝑣^𝑟𝜕𝜏^𝐼𝐽𝜕𝑟+𝑣^𝑧𝜕𝜏^𝐼𝐽𝜕𝑧,(3.6) where 𝐼𝐽=𝑟𝑟,𝑟𝑧,𝑧𝑧,𝜃𝜃. Assuming known values of 𝐶^𝐼𝐽, (3.2)–(3.4) may be represented as a set of matrix equations at each grid point including the boundary grids as follows:⎡⎢⎢⎢⎢⎢⎣𝜆Δ𝑡+−1𝜆−𝜕𝑣2_𝑧𝜆𝜕𝜕𝑣𝑟^𝑟−2𝜆𝜕𝜕𝑣𝑧^𝑟00𝜕𝑟Δ𝜆𝑡+1−−𝜆2𝜆𝜕𝜕𝜕𝑣𝑣𝑟^𝑟_𝑧+𝜕𝜕𝑣𝑧^𝑧−𝜆𝜕𝜕𝑣𝑧^𝑟𝜕𝑟Δ𝜆𝑡+1−2𝜆𝜕𝜕𝑣𝑧^𝑧⎤⎥⎥⎥⎥⎥⎦^{𝑛+1(𝑖𝑡)𝑖𝑘}⎡⎢⎢⎢⎢⎣𝜏^𝑟𝑟𝜏^𝑟𝑧𝜏^𝑧𝑧⎤⎥⎥⎥⎥⎦^{𝑛+1(𝑖𝑡+1)𝑖𝑘}=⎡⎢⎢⎢⎢⎢⎣𝜆𝜏^{𝑟𝑟(𝑛)}Δ𝑡−𝜆𝐶^𝑟𝑟+2(1−𝛽)𝜕𝑣^𝑟𝜆𝜏_{𝑟𝑧(𝑛)}𝜕𝑟Δ𝑡−𝜆𝐶^𝑟𝑧+(1−𝛽)𝜕𝑣^𝑧𝜆𝜕𝑟+𝜕𝜕𝑣𝑧^𝑟Δ𝑡𝜏^{𝑧𝑧(𝑛)}−𝜆𝐶^𝑧𝑧+2(1−𝛽)𝜕𝑣^𝑧𝜕𝑧⎤⎥⎥⎥⎥⎥⎦^{𝑛+1(𝑖𝑡)𝑖𝑘}.(3.7) On the other hand, from (3.5), 𝜏^{𝜃𝜃𝑖𝑘} can be found at each grid point within the computational domain and on the boundary as follows:𝜏^𝜃𝜃^{𝑛+1(𝑖𝑡+1)𝑖𝑘}=𝜆Δ𝑡𝜏^𝜃𝜃^𝑛−𝜆𝐶^𝜃𝜃+2(1−𝛽)𝑣^𝑟𝑟^{𝑛+1(𝑖𝑡)𝑖𝑘}Δ𝜆𝑡+1−2𝜆𝑣𝑟^𝑟^{𝑛+1(𝑖𝑡)𝑖𝑘}.(3.8) Once (𝜏^{𝑟𝑟𝑖𝑘},𝜏^{𝑟𝑧𝑖𝑘},𝜏^{𝑧𝑧𝑖𝑘},𝜏^{𝜃𝜃𝑖𝑘})^{𝑛+1(𝑖𝑡+1)} are obtained from (3.7)-(3.8), they are used to solve (2.8)-(2.9) to find (𝑣^𝑟,𝑣^𝑧)^{𝑛+1(𝑖𝑡+1)𝑖𝑘} in the 4 : 1 axisymmetric contraction geometry.

4. Results

In this section, we compare the results from the grid-by-grid inversion method with those of Phillips and Williams [11] to corroborate the accuracy and robustness of the grid-by-grid inversion method. To make the comparison meaningful, we adopt the same flow geometry and inlet conditions as those employed by Phillips and Williams [11]. Namely, the contraction ratio is 4 : 1, the length of the larger channel is 27𝑅, and that of smaller channel is 49𝑅 when 𝑅 is the radius of the smaller channel (cf. Figure 1). At the inlet, the parabolic Poiseuille flow is adopted as Phillips and Williams [11] have done.(-) inlet velocity𝑣^𝑧(𝑟)=16416−𝑟²,𝑣^𝑟=0,0≤𝑟≤4,(4.1) where 𝑟≡𝑟^/𝑅 and 𝑟′ is the dimensional radial distance. Contrary to Phillips and Williams [11], we do not impose inlet boundary conditions for the stress field, since the grid-by-grid inversion method does not require them and simply solve (3.7)-(3.8) at the inlet grid points to find the stress field at the inlet location. This is also an important merit of the grid-by-grid inversion method over other algorithms for the viscoelastic flows. Although Papanastasiou et al. [12] employed this method of imposing open outflow boundary conditions previously and named it open boundary condition method or free boundary condition method, it appears very naturally in the grid-by-grid inversion method.

We have employed three sets of grid system to ensure grid convergence of numerical results. For the three sets of grids, that is, 14,580 (Grid A), 32,535 (Grid B), and 57,960 (Grid C), the extrastress component 𝜏^𝑧𝑧 is plotted when 𝜆=1.0 for Re=0.0 along 𝑟=1 in Figure 2(a), and along 𝑟=0.064 in Figure 2(b). For there three mesh systems, the numbers of cells along 𝑟 in the small tube are 30, 45, and 60, respectively. For Grid A, the minimum Δ𝑟 is 0.005 at the wall and it grows gradually until it becomes 5.1×10⁻² at the center. For Grid B, (Δ𝑟)_min is 0.003 and (Δ𝑟)_max is 3.4×10⁻² at the center, and for Grid C, (Δ𝑟)_min=0.002 at the wall and (Δ𝑟)_max=2.3×10⁻² at the center. The results of the present work are also compared with those of Phillips and Williams [11]. Though the grid number increases from 32,535 to 57,960, there is not appreciable change in 𝜏^𝑧𝑧. Therefore, we adopt Grid B in the subsequent computations. Figure 3 shows the extrastress component 𝜏^𝑧𝑧 at Re=0.0 for 𝜆=0.5,1.0, and 1.5 along 𝑟=1.0 (Figure 2(a)) and along 𝑟=0.064 (Figure 2(b)). The grid-converged data for 𝜏^𝑧𝑧 obtained from the grid-by-grid inversion method are somewhat smaller than the benchmark data of Phillips and Williams [11]. At the corner of the contractor, 𝜏^𝑧𝑧 has a sharp overshoot, which increases with respect to 𝜆, and it settles down to a downstream value rapidly. Near the center of the contracted channel (𝑟=0.064), a smaller overshoot of 𝜏^𝑧𝑧 appears, which increases as 𝜆 increases, and it settles down to a downstream value much slowly as compared to that at the wall (𝑟=1.0).

(a)

(b)

Figure 2 

Confirmation of grid convergence of the grid-by-grid inversion method when Re=0.0 and 𝜆=1.0 (a) along 𝑟=1.0 (b) along 𝑟=0.064.

(a)

(b)

Figure 3 

Stress overshoot at various 𝜆 when Re=0.0 (a) along 𝑟=1.0 (b) along 𝑟=0.064.

Figure 4 shows comparison of streamlines for 𝜆=0.0, 0.5, 1.0, and 1.5 obtained from the grid-by-grid inversion method with those from the Phillips and Williams [11] when Re=0.0. The grid-by-grid inversion method yields accurate results as compared with the benchmark data for the range of 𝜆 values considered. The size of the corner vortex increases as 𝜆 increases. Figure 5 depicts comparison of streamlines at various values of 𝜆 when Re=1.0. The predictions of the grid-by-grid inversion method are also in good agreement with those of the benchmark data [11]. Comparing with the streamlines of Re=0.0, it is found that the vortex sizes are diminished due to the inertia effect [11]. The size of the corner vortex can be measured using a parameter 𝐿₁, as suggested by Phillips and Williams [11] and Oliveira et al.[8], which is the distance of the upstream separation point from the salient corner. Figure 6 shows the variation of 𝐿₁ with respect to 𝜆 for Re=0.0 and Re=1.0. For Re=1.0, the grid-by-grid inversion method is found to yield results similar to those of the benchmark data of Phillips and Williams [11]. For Re=0.0, the results of the grid-by-grid inversion method are compared with those of Oliveira et al. [8] as well as those of Phillips and Williams [11]. It is shown that both the grid-by-grid inversion and Phillips and Williams [11] yield 𝐿₁ somewhat larger than those of Oliveira et al. [8]. Since Oliveira et al. [8] do not consider the case Re=1.0, we cannot compare the grid-by-grid inversion method with Oliveira et al. when Re=1.0. To corroborate the accuracy of the grid-by-grid inversion method further, we compare the results for the vortex intensity Ψ_𝑅, which is defined as the flow rate in recirculation divided by inlet flow rate, in Figure 7 when Re=0.0. As in previous comparison, the grid-by-grid inversion method yields results similar to those of Phillips and Williams [11], while Oliveira et al. [8] predicts results somewhat smaller than the grid-by-grid inversion method. The grid-by-grid inversion method yields results for 𝜆 larger than 1.5 and Re larger than 1.0, which demonstrates the robustness of the grid-by-grid inversion method as compared to other finite volume methods or finite element methods [9, 11, 16, 18].

Figure 4 

Comparison of streamlines for 𝜆=0.0, 𝜆=0.5, 𝜆=1.0, and 𝜆=1.5 when Re=0.0.

Figure 5 

Comparison of streamlines for 𝜆=0.0, 𝜆=0.5, 𝜆=1.0, and 𝜆=1.5 when Re=1.0.

Figure 6 

Variation of the size of corner vortex 𝐿1 with respect to 𝜆 when Re=0.0 and Re=1.0.

Figure 7 

Variation of the dimensionless vortex intensity Ψ𝑅 with respect to 𝜆 when Re=0.0.

5. Conclusion

The grid-by-grid inversion method [10] has been applied to the 4 : 1 circular contraction problem in the present investigation. Viscoelastic flows through the contraction generate complex flows exhibiting strong shear and uniaxial extension, which is a good test bed for the robustness and accuracy of a new numerical algorithms. In the grid-by-grid inversion method [10], the hyperbolic constitutive equation is split such that the term for the convective transport of stress tensor is treated as a source, which is updated iteratively. This allows the stress tensor at each grid point to be expressed in terms of velocity gradient tensor at the same location, and the set of stress tensor components is found after inverting a small matrix at each grid point. The grid-by-grid inversion method can be implemented easily in any commercial CFD packages for Newtonian fluid flows to convert them to be used for viscoelastic fluid flows. The grid-by-grid inversion method is found to yield accurate results as compared with the benchmark data of Phillips and Williams [11]. It predicts accurately the variation of corner vortex size and the 𝜏^𝑧𝑧 overshoot with respect to 𝜆 at Re=0.0 and Re=1.0.

Acknowledgment

This paper was supported by the Human Resources Development of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) Grant funded by the Korea Government Ministry of Knowledge and Economy (no. 20114010203090).