Journal of Applied MathematicsVolume 2012 (2012), Article ID 420189, 15 pagesdoi:10.1155/2012/420189

Research Article

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

H. M. Park

Department of Chemical and Biomolecular Engineering, Sogang University, Seoul 121-742, Republic of Korea

Received 1 December 2011; Revised 5 January 2012; Accepted 15 January 2012

Academic Editor: M. F. El-Amin

Copyright © 2012 H. M. Park. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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: 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 is the Reynolds number. The superscript in (2.4) is the transpose. The Reynolds number and the dimensionless Deborah number are defined by 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 .

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 , and , using known velocity and stress at time step and , proceed in an iterative manner as follows. First, define where the superscript indicates variable at time step in the 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. Once is evaluated using variables at , (2.7) can be solved for 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 obtained from (2.7), the velocity at the th time step in the iteration, , is found by solving (2.1)–(2.3) as follows: 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 requires a special deliberation before employing standard algorithms for incompressible Navier-Stokes equations in the grid-by-grid inversion technique. When , we decompose the total stress into pure viscous and elastic parts and , respectively, such that where Then the momentum and constitutive equations for the case of may be written as 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 Next, convert (2.14) to the following local matrix equations for decoupled at each grid point: Then, (2.12) and (2.13) are discretized as follows: The structure of (2.16) is the same as that of (2.7) and can be solved for by inverting a six-by-six matrix at each grid point for a three-dimensional flow once is given. After obtaining , 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 begins as follows.(2) Assume and . For the first iteration (), = and =.(3) Evaluate of (2.6) () or of (2.15) () using an upwind scheme.(4) Using , solve (2.7) () or (2.16) () for by inverting a six-by-six matrix at each grid point including the boundary grids.(5) Using ,

() solve the momentum equation, (2.9), and continuity equation to find .

() solve the momentum equation, (2.18), and continuity equation to find .(6) Convergence check for and . 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: 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: 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 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: On the other hand, from (3.5), can be found at each grid point within the computational domain and on the boundary as follows: Once are obtained from (3.7)-(3.8), they are used to solve (2.8)-(2.9) to find 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 , and that of smaller channel is 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 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 for along in Figure 2(a), and along 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 at the center. For Grid B, is 0.003 and is at the center, and for Grid C, at the wall and 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 for , and 1.5 along (Figure 2(a)) and along (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 , 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 ().

(a)

(b)

Figure 2: Confirmation of grid convergence of the grid-by-grid inversion method when and (a) along (b) along .

(a)

(b)

Figure 3: Stress overshoot at various when (a) along (b) along .

Figure 4 shows comparison of streamlines for , 0.5, 1.0, and 1.5 obtained from the grid-by-grid inversion method with those from the Phillips and Williams [11] when . 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 . 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 , 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 and . For , the grid-by-grid inversion method is found to yield results similar to those of the benchmark data of Phillips and Williams [11]. For , 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 , we cannot compare the grid-by-grid inversion method with Oliveira et al. when . 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 . 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 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 , , , and when .

Figure 5: Comparison of streamlines for , , , and when .

Figure 6: Variation of the size of corner vortex with respect to when and .

Figure 7: Variation of the dimensionless vortex intensity with respect to when .

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 and .

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).

References

R. B. Bird, R. C. Amstrong, and O. Hassager, Dynamics of Polymeric Liquids, vol. 1, Wiley, New York, NY, USA, 1987.

R. R. Huilgol and N. Phan-Thien, Fluid Mechanics of Viscoelasticity, Elsevier, Amsterdam, The Netherlands, 1997.

D. D. Joseph, Fluid Dynamics of Viscoelastic Liquids, vol. 84 of Applied Mathematical Sciences, Springer, 1990.

M. A. Alves, F. T. Pinho, and P. J. Oliveira, “The flow of viscoelastic fluids past a cylinder: finite-volume high-resolutions method,” Journal of Non-Newtonian Fluid Mechanics, vol. 97, no. 2-3, pp. 207–232, 2001.

M. J. Crochet, A. R. Davis, and K. Walters, Numerical Simulation of Non-Newtonian Flow, Elsevier, Amsterdam, The Netherlands, 1984.

F. Yurun and M. J. Crochet, “High-order finite element methods for steady viscoelastic flows,” Journal of Non-Newtonian Fluid Mechanics, vol. 57, no. 2-3, pp. 283–311, 1995.

B. Khomami, K. K. Talwar, and H. K. Ganpule, “A comparative study of higher- and lower-order finite element techniques for computation of viscoelastic flows,” Journal of Rheology, vol. 38, no. 2, pp. 255–289, 1994.

M. S. N. Oliveira, P. J. Oliveira, F. T. Pinho, and M. A. Alves, “Effect of contraction ratio upon viscoelastic flow in contractions: the axisymmetric case,” Journal of Non-Newtonian Fluid Mechanics, vol. 147, no. 1-2, pp. 92–108, 2007.

R. G. Owens and T. N. Phillips, Computational Rheology, Imperial College Press, London, UK, 2002.

H. M. Park and J. Y. Lim, “A new numerical algorithm for viscoelastic fluid flows: the grid-by-grid inversion method,” Journal of Non-Newtonian Fluid Mechanics, vol. 165, no. 5-6, pp. 238–246, 2010.

T. N. Phillips and A. J. Williams, “Comparison of creeping and inertial flow of an Oldroyd B fluid through planar and axisymmetric contractions,” Journal of Non-Newtonian Fluid Mechanics, vol. 108, no. 1–3, pp. 25–47, 2002.

T. C. Papanastasiou, N. Malamataris, and K. Ellwood, “A new outflow boundary condition,” International Journal for Numerical Methods in Fluids, vol. 14, no. 5, pp. 587–608, 1992.

S. V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, NY, USA, 1980.

H. K. Versteeg and W. Malalasekera, An Introduction to Computational Fluid Dynamics, Pearson Education Limited, Essex, UK, 2007.

M. G. N. Perera and K. Walters, “Long-range memory effects in flows involving abrupt changes in geometry: part1. Flows associated with L-shaped and T-shaped geometries,” Journal of Non-Newtonian Fluid Mechanics, vol. 2, no. 2, pp. 191–204, 1977.

T. N. Phillips and A. J. Williams, “Viscoelastic flow through a planar contraction using a semi-Lagrangian finite volume method,” Journal of Non-Newtonian Fluid Mechanics, vol. 87, no. 2-3, pp. 215–246, 1999.

R. Keunings and M. J. Crochet, “Numerical simulation of the flow of a viscoelastic fluid through an abrupt contraction,” Journal of Non-Newtonian Fluid Mechanics, vol. 14, pp. 279–299, 1984.

D. Trebotich, P. Colella, and G. H. Miller, “A stable and convergent scheme for viscoelastic flow in contraction channels,” Journal of Computational Physics, vol. 205, no. 1, pp. 315–342, 2005.

C. M. Rhie and W. L. Chow, “Numerical study of the turbulent flow past an airfoil with trailing edge separation,” AIAA Journal, vol. 21, no. 11, pp. 1525–1532, 1983.