p-Version least squares finite element formulation for axisymmetric incompressible non-Newtonian fluid flow

+__li!B

2s

Computer methods in applied mechanics and engineering Comput. Methods Appl. Mech. Engrg. 113 (1994) 271-300

EtsEVIER

p-Version least squares finite element formulation for axisymmetric incompressible non-Newtonian fluid flow Nathan B. Edgar, Karan S. Surana* The University of Kansas, Department of Mechanical Engineering,

3013 Learned

Hall, Lawrence,

KS 66045, USA

Received 16 July 1992 Revised manuscript received 16 June 1993

Abstract The paper presents a p-version least squares finite element formulation (LSFEF) for axjsymmetric incompressible nonNewtonian Auid flow. The dimensionless form of the equations describing the fluid motion are cast into a set of first order coupled partial differential equations involving non-Newtonian stresses as auxiliary variables. The pressure, velocities (primary variables) and the stresses (auxiliary variables) are interpolated using equal order, Co continuity, p-version hierarchical approximation functions. The least squares functional (or error functional) is constructed using the system of coupled first order non-linear partial differential equations (integrated sum of squares of the errors resulting from the individual equations for the entire discretization) without linearization, approximations or assumptions. The minimization of this least squares error functional results in finding a solution vector {S } for which the partial derivatives of the error ~nctional with respect to the nodal degrees of freedom (6) becomes zero. This is accomplished by using Newton’s method with a line search. The paper presents an implementation of the power law model for non-Newtonian viscosity. Numerical examples are presented to demonstrate the convergence characteristics and the accuracy of the proposed formulation.

The finite element method has emerged as one of the most general and versatile numerical simulation techniques for systems described by differential or partial differential equations. Even .though its application in computational fluid dynamics (CFD) is rather recent (20 years), the method has already established its superiority over the finite difference and control volume based techniques. The variational, Galerkin, and more recently collocation and least squares methods have been applied in CFD. The variational methods provide the best approximation to the exact solution of the associated variational problem. However, the variational form in general is difficult (or may even be impossible) to construct for non-linear partial differential equations describing fluid motion. Galerkin, collocation and least squares processes are special cases of the weighted residual approach and provide possible alternatives to variational approach. Galerkin based formulations when applied to CFD using pressure and velocities as variables (p~mitive variables) lead to many difficulties, a summa~ of which can be found in [l]. The published literature on the finite element formulations and solution methods is enormous. In the following we provide a brief discussion of some pertinent literature. Hughes et al. [2] in 1979 presented a finite element formulation of incompressible viscous flows using penalty function formulation. The paper presents a discussion of Galerkin and ‘upwind’ treatments of the convective terms. Brooks and Hughes [3] in 1982 presented streamline u~wind/Petrov-Galerkin finite element form~ations with particular emphasis on the incompressible Navier-Stakes equations. In

* Corresponding

author.

0045-7825 /94/$~.~

SSDI

@ 1994 Elsevier Science B.V. All rights reserved

0045-7825(93)EOl21-N

272

N.3.

Edgar. KS.

Surana ! Compuf.

Methods Appl.

Mech. Engrg. 113 (1994) 2li-300

1986, Hughes et al. [4] presented a Petrov-Galerkin finite element formulation of Stokes problem accommodating equal order Co interpolations which circumvented the BabuSka-Brezzi condition. Hughes and Franca (51 in 1987 presented a finite element formulation for Stokes flow with various well posed boundary conditions. The formulation produced symmetric matrices and the authors claimed convergence for all velocity/pressure spaces. The formulations described in [2-51 address many techniques of overcoming the difficulties described in [l] when using Galerkin based formulations using velocities and pressure as variables. Hirsch and Warzee [6] presented a finite element formulation based on the Galerkin method for axisymmetric flow computations in a turbo-machine. The equations for the meridional through-bow were formulated as a quasi-ha~onic non-linear equation. The non-linear equations resulting from the Galerkin process were solved iteratively using under-relaxation. Numerical simulation results were compared with experimental data. The finite element solution of viscous, incompressible, Newtonian jet flows with surface tension using the Galerkin method was presented by Reddy and Tanner [7]. In their work the resulting system of non-linear algebraic equations was solved using the method of successive substitution or the Picard method. A Galerkin based finite element formulation for an axisymmetric, Newtonian jet of incompressible viscous fluid with surface tension to study die-swell was presented by Omodei [8]. The numerical results were also compared with experimental results. Yoseph et al. [9] presented a finite element formulation of the Navier-Stokes equations in rotating flow using (1) Galerkin’s approach with application of Green’s theorem to velocity terms, (2) application of Green’s theorem to velocity and pressure terms and (3) by introducing a penalty parameter in the second approach. A number of numerical studies are presented to compare these approaches. Donea et al. (101 presented an explanation of the significance of discontinuous pressure at the interelement boundaries in a Galerkin finite element formulation of the Navier-Stokes equations when using 4-noded isoparametric elements. They suggested an alternate scheme using 4-noded non-conforming elements, with strict satisfaction of the tangential velocity boundary condition which produced chequerboard free and correct pressure values everywhere. Carey and Krishnan [11] utilized a continuation technique in Reynolds number to obtain the solution of non-linear algebraic equations resulting from penalty finite element approximation of Navier-Stokes equations. Numerical results were presented for the driven cavity problem. The approach was also extended to arc length continuation to treat problems with singular points on the continuation path. The least squares finite element method (LSFEM) overcomes many difficulties present in other formulation procedures and offers a more general and assumption free platform for treating problems described by non-linear coupled partial differential equations, The first successful application of the LSFEM to linear problems described by self-adjoint (elastic torsion of prismatic bar), non-self-adjoint, biharmonic (plane stress) differential equation and plane steady state flow of a viscous incompressible fluid with negligible inertia terms (without convective terms) was presented by Lynn and Arya [12]. In each case, the linear differential equation or the system of coupled linear differential equations were reduced to a coupled system of first order partial differential equations by introducing auxiliary variables. The least squares finite element formulation was then constructed for the system of first order coupled partial differential equations using equal order interpolation for both primary and auxiliary variables. In 1974, Lynn and Arya [13] presented a weighted discrete or element-wise least squares minimization of differential equation residuals for linear problems. In this paper: (1) the equations were cast as a system of first order coupled partial differential equations using auxiliary variables and (2) the smoothing of discontinuous trial functions was employed. The latter essentialfy relaxes the continuity requirements and yields non-confo~irig finite elements. Zienkiewicz and Owen 1141 presented least squares finite element formulation for elasto-statics problems using the approach of Lynn and Arya [12] and advocated the use of ‘reduced integration’ as a means of improving the accuracy of the computed solution. The first application of the LSFEM to non-linear differential equations was presented by Lynn [15] for laminar boundary layer flows using the linearized form of the Blasius non-finear differential equation for incompressible viscous flows afong a plate with similarity property. Eason and Monte [16] presented a solution of the non-linear boundary value problems by the discrete least squares finite element method. In this formulation, the feast squares integrals are replaced by the sum of squares of the residuals at individual points and the numerical solution is obtained using Powell’s method [17].

N.B.

Edgar,

K.S. Surana

I Comput.

Methods Appl. Mech. Engrg.

113 (1994) 271-300

213

Polk and Lynn [18] presented a LSFEF for one-dimensional unsteady gas dynamics using a space-time finite element approach in which an approximate least squares functional was constructed by integrating the sum of squares of the residuals in space and time using a linearized form of the equations. A least squares finite element formulation for linear problems of mixed type was given by Fix and Gunzburger [ 191. Fletcher [20], in 1979, presented a primitive variable (using velocities and pressure) LSFEF for inviscid compressible flow. A unique feature of this formulation was the representation of groups of variables rather than single variables. Bristeau et al. [21] presented considerable mathematical details on the formulation and numerical solution of fluid dynamics problems using a LSFEM and a conjugate gradient solution method. The least squares based space-time finite element formulations were also presented by Nguyen and Reynen [22,23] for the convection-diffusion and Burgers’ equations. The authors begin the formulation procedures with a true least squares functional in space and time assuming space and time effects to be coupled. However, in the process of obtaining element properties, assumptions and approximations were made which altered the exact nature of the least squares functional in space and time. First, we note that both for the convection-diffusion and Burgers’ equations, the expressions for E (error) involve second derivatives of the dependent variable. Thus a linear finite element approximation for the dependent variable in space is unsuitable and will destroy the true form of the least squares functional. It could be argued that 6E (variation of E) may be treated like a weight function wj (which is a function of space and time) and that the least squares method is a special case of the Petrov-Galerkin technique. Then of course any suitable form for wj is permissible. However, the formulation resulting for any other wj except 6E will not represent a true least squares space-time finite element formulation. Second, the authors perform integration by parts for some terms involving second derivatives and drop the others due to linear approximation. Last, it is important to point out that the authors are able to make this procedure a marching procedure in time only by introducing A4i (change in the dependent variable 4 over the time step). This is only possible for a linear approximation in time. For higher order approximations in time, the procedure presented in [22,23] cannot be utilized as a time-stepping procedure. In conclusion, the formulations presented in [22,23] are not true least squares space-time formulations and the time stepping is possible only for linear approximation in time. Tabarrok and Saghir [24] presented a mixed formulation using a SLFEM for two-dimensional incompressible flows. The authors utilized the ‘true’ form of the least squares functional in the minimization process. The conditions resulting from the least squares minimizations were arranged in the form of a system of simultaneous algebraic equations such that all non-linear terms were moved to the right side making the coefficient matrix linear (independent of field variables) and symmetric. Since the authors did not use auxiliary variables to reduce continuity requirements on the approximation functions, the formulation required the use of C’ elements. Aziz et al. [25] presented a weighted LSFEF for elliptic systems of partial differential equations along with error analysis and some numerical examples. Jiang and Carey [26] presented a LSFEF for the Laplace equation. The resulting equations were solved using an element-by-element conjugate gradient solution procedure. The authors also presented an adaptive mesh refinement strategy based on mean square element error functionals (element error functional/ element area). In another paper, Carey and Jiang [27] presented LSFEF for a linear first order equation and elliptic system described by a pair of first order equations. The authors advocated the use of a preconditioned conjugate gradient procedure as the solution method. In 1987, Carey and Jiang [28] presented a LSFEF for non-linear partial differential equations. Unfortunately, the equations were linearized by lagging coefficients (see [28, Eq. (15)]) which represented an approximation to the true stationary condition (as admitted by the authors). Thus, this formulation does not deal with the true least squares functional that corresponds to the actual non-linear partial differential equations at hand. In this paper the authors also advocated a preconditioned conjugate gradient method for computing the numerical solution. Jiang and Carey [29] also presented a least squares finite element formulation for non-linear hyperbolic problems. In this work, the equations are first linearized over the time step At=t”+‘t”. A backward difference approximation is then used to obtain an implicit time-difference problem for which the least squares finite element approximation is constructed in space using the functional Z = ], (R* + P(aR/ax)*) da, 0 s /3 c 1 instead of Z = ], R* dR where R represents the

274

N.B.

Edgar,

KS.

Surana

I Comput. Methods Appl. Mech. Engrg.

113 (1994) 2716?cNI

residuals resulting from the linearized time discretized equations. Numerical examples include inviscid Burgers’ equation, isothermal flow in a nozzle, and a shock tube problem. In a series of papers, Kececioglu and Rubinsky [30-32] presented a mixed variable continuously deforming finite element method based on a least squares approach for evolution problems. These papers discuss theoretical aspects, their application to problems of phase-change in a porous media, and numerical implementation. The authors advocate a true space-time coupled least squares minimization approach in theory, but, the effects of space and time are actually decoupled in the formulations presented in the papers. The authors first discretized temporal derivatives using difference equations in time and then constructed LSFEF in space using the decoupled equations. Jiang and Chang [33] presented a LSFEF for Stokes flow. The equations of fluid motion were recast using pressure, velocity and vorticity (an auxiliary variable: {w} = curl{u}) as variables which yields a system of ‘linear’ coupled partial differential equations for which the LSFEF can be easily constructed without any approximations. In another paper, Jiang and Povinelli [34] presented applications of the LSFEM to various one- and two-dimensional flows. Unfortunately, only the linearized form of the equations of fluid motion were used in constru~ing the least squares error functional to be minimized. Thus these formulations only represent an approximation to the true least squares error functional which must be derived using actual non-linear equations of fluid motion. From the published literature, we note that even though the LSFEF has been advocated as a general tool for CFD, a general presentation of the formulation procedure for non-linear problems applicable to all classes of non-linear problem regardless of the nature and the severity of the non-linearities is rather hard to find. Almost all published works (except 1241) start the formulation procedure with a true least squares error functional which corresponds to the actual non-linear differential equations, but when it comes down to the details of the formulation for specific problems, the partial differential equations are linearized before constructing the error functional to be minimized. All such formulations represent only an approximation to the true least squares error functional which must be minimized. The h, p and h-p versions of the finite element method are well known procedures for improving the accuracy of the solution. The p-version of finite element method possesses superior convergence characteristics compared to h-version [35]. For problems involving singularities, the h-p version is necessary ta achieve the fastest rate of convergence [35]. Jiang and Sonnad [36] applied the p-version least squares finite element method to incompressible fluid flow using a vorticity approach and the approximation functions based on Legendre polynomials. Winterscheidt and Surana [ 1,37,38] have reported p-version least squares finite element formulations for convection-diffusion, Burger’s equation and two-dimensional incompressible fluid flow. In their work, the equations are recast into a set of first order coupled partial differential equations by introducing auxiliary variables. The least s uares finite element formulation is constructed for this set of first order equations using equal order C?i, p-version hierarchical approximation functions for both primary and auxiliary variables without linearizing the equations or introducing any other assumptions or approximations. The vast majority of the published literature in finite element techniques for CFD is concentrated in the area of Newtonian lluid how using low order element approximations. In a recent paper, Bell and Surana [39] presented a p-version least squares formulation for two-dimensional incompressible non-Newtonian fluid flow by using equal order Co, p-version hierarchical approximation functions for both primary and auxiliary variables. In Bell and Surana’s work [39], also, the least squares error functional to be minimized was constructed using actual non-linear partial differential equations without assumptions or approximations. In this paper, we present a p-version squares finite element formulation for axisymmet~c, incompressible, non-Newtonian fluid flow. The dimensionless form of the equations of axisymmetric fluid motion in the cylindrical coordinate system are recast into a set of first order coupled partial differential equations using non-Newtonian stresses as auxiliary variables for which a p-version least squares finite element formulation is constructed. The least squares error functional used in the minimization process is the true least squares functional corresponding to the actual non-linear equations describing the non-Newtonian axisymmetric flow of a power law fluid. The pressure, velocities and the stresses are interpolated over an element ,using equal order, Co, p-version hierarchical approximation functions. A solution vector {a}, which satisfies the conditions resulting from the least squares minimization process,

N.B.

Edgar,

K.S. Surana

I Comput. Methods Appl. Mech. Engrg.

113 (1994) 271-300

275

is found by using Newton’s method with a line search. Numerical examples are presented for the power law model. The results obtained from the present formulation are compared with analytical solutions and the results reported in the literature.

2. Least squares finite element formulation The following is a summary of the p-version described by non-linear differential equations. functional I’ for an element e is defined as I' = i jj=, 0’

(E;)’ dR

least squares finite element formulation for problems For a system of N differential equations, the error

,

where Ey, i = 1,2, . . . , N are the errors which result when the finite element approximation to the true solution is substituted into the differential equations. For a finite element mesh consisting of NE elements, the total error functional I for the entire mesh is obtained by summing up the element error functionals I’.

t?=I

From Eq. (2), we note that Z is a function of nodal degrees of freedom (6) and thus minimization would require that we differentiate Z with respect to (6) and set it to zero.

of Z

(3) where

(4) We note that {g} in Eq. (3) is indeed the variation of Z and thus minimization 6Z= F 61’ = {g} = (0)

)

of Z would require (5)

e=l

where

with (7) Eqs. (3) or (5) represent the true least squares minimization statement. For a (6) to minimize I, {g} must become a null vector. The least squares minimization statement ends here. In other words, in the least squares process we minimize Z given by Eq. (2) which results in finding a (6 } that makes {g} a null vector. In the least squares statement presented above there are no assumptions or approximations and the non-linear terms are treated properly (see derivation presented later), for example terms like u au/ax are differentiated with respect to {S} using the product rule while formulating {g} in Eq. (3) as opposed to linearizing them by lagging u. It is important to note that finding a vector (8) which satisfies {g} = (0) is a separate issue and has nothing to do with the least squares error functional and its minimization. It may also be noted that the least squares minimization principle stated above is exact and is applicable to any non-linear system of partial

N.B.

276

Edgar,

KS.

Surana

i Comput. Methods Appl. Me&.

differential equations regardless of the nature of non-linearities the non-linear terms.

Engrg.

113 (1994) 2lI-300

and it does not require linearization

of

3. Solution procedure We consider problems described by non-linear differential partial equations. For such problems, {g} is a non-linear function of {S} and our objective is to find a {S} which would make {g} a null vector. At this point several different strategies can be adopted. In the first approach, {g} = (0) may be recast in the following forms:

In Eq. (Sa), the coefficient matrix [K] is non-symmetric and in general its coefficients are non-linear functions of (6). In Eq. (8b), the coefficient matrix [K] is symmetric and is independent of (6). In this form, all (6) dependent terms are contained in the vector {f}. In the third form given by Eq. (8c), the coefficient matrix [KJ may or may not be symmetric (depending upon which terms have been transferred to the right hand side vector {f}), but both [#I and (f} are functions of {a}. The solution vector (6) may be calculated by operating on Eq. (8) iteratively. In the second scheme, we operate on {g} directly. The element by element conjugate gradient (including preconditioned conjugate gradient) method, and methods based on the Taylor series expansion of {g} (Newton’s method) fall into this second category. Here we utilize and present the details of a solution method based on the second approach. Since the actual solution (6) which makes {g) a null vector is not known, we assume a starting solution {S,} for which in general we have {g({%J)) + (01 * Let {A??} be the correction M&J

(9) to {S,} such that

+ (As))} = (01.

(101

We expand (10) in Taylor series about {S,} and retain only the linear terms in (A.6): M&J

+ {Aa))) = {g({&J)} +#

From (ll),

we can write

Substituting

from (3) into (12), we obtain

(sol {Aa) = (01 -

(11)

(13) Let (14) The coefficient matrix [H”] is called the element Hessian matrix. Using the notation in Eq. (14), we can rewrite Eqs. (12) or (13),

EQo, (ASI = -M{4J)I

,

(15)

N.B.

Edgar,

KS.

Surana

I Comput. Methods Appl. Mech. Engrg.

113 (1994) 271-300

271

where

[HI=

;l . WI

(16)

Using Eq. (4) we can write the following expression for the element Hessian matrix [H”]: (17) From Eq. (17) we note that the element Hessian matrix [H”] is symmetric. The improved value of the solution is given by (6 > = {V + {Aa > .

(18)

At this point we note that if the incremental change {AS} is too large, updating the solution vector (8,) according to Eq. (18) may result in an increase rather than a decrease in the error functional 1. Thus we update the solution vector {a,,} using {6}={6,}+a{A~}.

(19)

The scalar (Yis selected to minimize the error functional I. In our numerical computations we consider the range of CYas 0 < (Y< 2. In this range, (Yis incremented in steps of 0.1 and for each value of (Ythe error functional Z is calculated. The lowest value of Z and a value on either side of the lowest value are used to construct a parabolic fit from which a value of (Yis calculated which minimizes I. This value of CY is used in Eq. (19) to obtain the updated vector {S} for the next iteration. This procedure described above is termed ‘Newton’s method with a line search’. When (Y= 1, it reduces to the classical Newton’s method. In the present work, we use a starting vector {S,} = (0) an d (Y= 1 for the first iteration, the line search is used for subsequent iterations. In all numerical studies conducted (including the ones presented in this paper), the actual range of (I! was 0 < a G 1.5. When the updated solution is in the very close neighborhood of the true solution, (Y becomes very close to 1 and eventually reaches 1 at convergence. For a given mesh and p-levels, iterations are performed until the error functional Z and each component of {g} are below a prespecified tolerance or cease to decrease below a certain value. More details on the number of iterations and the actual tolerances used are presented with the numerical examples. 3.1. Comments and discussion (1) Some discussion regarding the element Hessian matrix is in order first. Numerical studies have shown that the second term in the integrand of Eq. (17) which involves the second derivative of El with respect to {ae} can be neglected without affecting the convergence of the iterative solution procedure. In fact in many numerical studies the convergence of the iterative procedure is improved by doing so. In problems where convective terms are the only source of non-linearities (Burger’s equation, incompressible Newtonian fluid flow, etc.), obtaining explicit expressions for the second derivatives of EL with respect to {ae} is not too terribly complicated, however, in applications like non-Newtonian fluid flow where the generalized viscosity may be represented by power law behavior, obtaining expressions for this term is algebraically very tedious and in some cases may even become a formidable task. Additional computational time needed in computing [H’] due to inclusion of this term depends upon the degree of non-linearity. For two-dimensional incompressible Newtonian flow, the increase in computational time is only a small fraction; but, for non-Newtonian fluid flow applications, inclusion of this term may significantly increase the computational time for [He]. Thus, in summary, including the second term in the computation of [H’] increases the complexity of programming, increases computational effort and does not improve the convergence of the iterative process. Even though at this point we are unable to offer a mathematical justification for not including the second term in the computation of [H’], the benefits described here amply justify computation of [H”] with the second term. (2) It should be noted that this approximation (described in item (1) above) in the computation of

278

N.B. Edgar, KS. Surana ! Comput. Methods Appl. Mech. Engrg, 11.3 (IW)

271-300

the element Hessian matrix only alters the search direction in the iterative procedure and in no way affects the least squares error functional and the minimization procedure described in the earlier section. (3) We note that {g} = (0) 1s . only a necessary condition for the minimization of I. The sufficient condition to ensure that we have the global minima of I is that [H] (given by Eq. (16)) must be positive definite for a (6) which satisfies {g} = (0). That is, we must show that {a*}’ [H( {6})]{6*} > 0 for an arbitrary vector (8 *}. Numerically this is rather easy to compute and demonstrate (even though we have not done so in this paper). Rather than proving or numerically showing that [H({6})] is positive definite, we present a different reasoning and state another condition which ensures that a (8) which satisfies {g} = {0} indeed minimizes I. First we note that Z represents the summation over all of the elements of the mesh of the sum of squares of the errors resulting from the individual differential equations integrated over each element. Let A represent the threshold value for numerically computed zero (usually A = 10-6-10-12 depending upon the word size of the machine) i.e., all computed quantities that are less than or equal to A in magnitude are considered zero. In our iterative procedure, the solution is considered converged to the correct solution when the magnitude of each component of (g} is less than or equal to A and when Z 6 A. We note that Z < A or I = 0 is the minimum value of Z we are seeking since it corresponds to zero error in each differential equation at each integration point of each element of the mesh. Thus, in conclusion, when we find a {S} for which {g} = (0) and Z = 0 (both in the sense of threshold A) we have indeed found a (6) for which Z is minimum because Z = 0 is the minima we are seeking. Thus {g} = {0} and Z = 0 are necessary and sufficient conditions for the minima of 1. {g) z (0) resulting in Z > 0 indicates either lack of mesh refinement or inadequate p-levels or in some cases both. (4) The p-version LSFEF, when applied to the system of non-linear partial differential equation results in finding a (6) for which {g} must become a null vector. But, since {g} may be a complex non-linear function of {S}, we apply Newton’s method with a line search to find this {S}. In many pubhshed papers, the non-linear differential equations are linearized at the onset before constructing the least squares functional [6,9,13]. Such approaches result in an approximation to the true least squares error functional. Clearly, as just demonstrated, such linearization is totally unwarranted as far as the least squares minimization functional is concerned. Another important comment regarding the linearization of differential equations is that this can only be done in simple cases (like convective terms in momentum equations). In more complex situations like non-Newtonian fluid flow, the linearization approaches are similar to those presented and advocated in [6,9,13] and others may result in undesirable approximation or error. We argue that, first, the linearization of the differential equations is an unnecessary approximation which can be easily avoided and secondly the least squares functional corresponding to the linearized differential equations is not the true functional which corresponds to the actual non-linear equations. Any numerical scheme based on this approximate functional will result in an approximation in the calcmated solution as well. In other words, a {S} which makes the error functional Z resulting from the linearized system zero may resuh in the true error functional (based on actual non-linear differential equations) Z > 0, indicating that this {S} is not the solution of the actual non-linear problem. (5) We also wish to point out that application of Newton’s method to find a (6) for which {g} becomes a null vector corresponding to the true least squares functional presented here is not the same as linearizing the differential equations and then constructing the Ieast squares functional followed by its minimization. As we stated earlier, construction of the Ieast squares functional and the associated minimization and the method to find a solution vector (6) that satisfies the conditions dictated by the least squares minimization principle are two separate issues. The point to note is that linearized differential equations result in a functional that does not correspond to the statement of the original problem. In our presentation we have demonstrated that true least squares functional and minimization principle can be constructed without introducing any approximations in the differential equations and that Newton’s method with a linear search can indeed be used to find the solution vector (6) for which {g} = (0) and Z = 0 which in fact is the minima of Z we are seeking. The approximation in the computation of the Hessian matrix only affects the search direction during iterations and has absolutely no effect on the exactness of the least squares functional and the minimization principle.

N.B.

Edgar,

K.S. Surana

I Comput. Methods Appl. Mech. Engrg.

113 (1994) 271-300

219

4. Equations of motion for axisymmetric fluid flow The continuity, momentum and constitutive equations (total stresses for inelastic non-Newtonian fluids) for steady-state axisymmetric incompressible fluid flow can be written in the following form using the cylindrical coordinate system: .

n

1

S+;+$=o, ~

aii

n

ali

u-g+lJ-g+a;-

ai)

at ai

(a,,+++,__ A ai

*

26

ti

i ( 7 >> = 0,

(20)

The non-Newtonian fluid viscosity ;I depends upon the flow field parameters. The general form of the constitutive equations in (20) is the same as that for the Newtonian fluids and is often referred to as generalized Newtonian fluid models. The specific form of +j is chosen to describe the behavior of the fluid of interest. When 4 is constant for the entire flow field, we have a Newtonian fluid. There are many non-Newtonian fluid models described in the literature [40,41]. In the present work we consider the power law model which is quite commonly used for many fluids. 4.1. Power law model for ij The power law model is one of the most commonly used models for non-Newtonian viscosity because of its ability to describe both shear thinning as well as shear thickening fluid behavior. For power law fluids, we define ,+j= @#n-U/Z

)

where I2 is the second invariant of the strain rate tensor. For axisymmetric following explicit form of j2:

(21)

flow, we can write the

(22) where n is a dimensionless quantity called ‘power index’ and h is called ‘consistency’. For n < 1, the fluid is shear thinning, for n > 1, the fluid is shear thickening, and for 12= 1, Eq. (21) reduces to Newtonian fluid behavior.

5. Dimensionless form of the equations of fluid motion It is possible to consider two different dimensionless forms of Eqs. (20)-(22) which arise as a result of two different choices of the characteristic scale for the stresses. If the characteristic viscous stress (&uO /D) is chosen to scale stresses, the Reynolds number appears in the momentum equations. If the characteristic kinetic energy (p,-,u%)is used to non-dimensionalize the stresses, the Reynolds number appears in the stress equations. The choice of the dimensionless form of these equations has a dramatic effect on the convergence of the iterative process used for solving the non-linear discretized equations of equilibrium resulting from the least squares process. A detailed discussion of these issues along with

280

N.B. Edgar,

K.S. Surana I Comput. Methods Appi. Mech. Engrg.

113 (iYY4) 271-300

a numerical example illustrating the effect of the two choices of non-dimensionalizing the stress equations can be found in a recent paper by Bell and Surana [39]. Details are not repeated here for the sake of brevity. Bell and Surana [39] showed that when the characteristic viscous stress (~u~/~) is used to non-dimensionalize the stresses, the resulting equations lead to slow or non-convergence of the iterative solution procedure or convergence to spurious solutions for a variety of problems, even for low to moderate Reynolds numbers. In the present work, we use characteristic kinetic energy to non-dimensionalize the stresses. The following dimensionless variables are used:

(23)

Using these dimensionless form:

variables, Eqs . (20)-( 22) can be expressed in the following dimensionless

$+:+$o, 14$+vg+dr-

u$+ug+?zT..-2vv() $

C3P (

aP

27l

~+++~_~

2.%+~+eg (

>

(

;

1)

=o,

=a, (241

=o,

TJ"(~+~)=o, 7zz-27&($) The dimensionless 7” = j$

=o. non-Newtonian

viscosity rl, described by the power law can be written as

(zZ)(n-1)‘2 )

(25)

1Z=2(!$>Z+2(~)2+2($)2+($+$)2,

(26)

where

where Re, is the Reynolds number for the power law fluid and is given by Re, = MJF”f)”

*

m0

(27)

6. LSFEF using the dimensionless form of the equations of motion Eqs. (24) represent a set of six coupled first order partial differential equations. The purpose of considering the set of coupled first order equations, Eq. (24), is to permit the use of Co element approximation functions [l, 37,381. When the finite element approximations

N.B.

Edgar,

K.S. Surana

I Comput.

Methods Appl. Mech. Engrg.

113 (1994) 271-300

281

auh uh+avh

Ee

I

ar

I-

r

a.2’

and 77”

=$

(#-I)/2

)

1:=2(g)2+2($-)2+2(+$)2+ 6.1.

(5+$).

(29)

Element approximation

If we consider the same order of approximation Gh =

for P, u, v, T,,, r,, and T,,, then we can write

[W@~ 7

(30)

where SDis the field variable, [N] is the hierarchical approximation function matrix and {CD}is a vector of the hierarchical nodal degrees of freedom for the field variable @. Now, by letting @ equal P, u, v, ‘rL) T,, and T,* we can write the element approximation for each of the seven variables. For example, etc. Here we consider a nine node quadrilateral element with curved u =[N] {u},vh=[N] {v},..., sides. The element is mapped into a two unit square in the 5, n natural coordinate space (Fig. 1). The

a

Nodes with hierarchical degrees of freedom

0

Non-hierarchical nodes

Fig. 1. A nine node p-version axisymmetric element.

hierarchical approximation functions and hierarchical nadal variables are constructed in the natural coordinate system. Details of this procedure can be found in [l? 371 and are omitted here for the sake of brevity. In addition to the errors given by Eqs. (28), the least squares fo~ulat~on presented in the previous section requires expressions for the derivatives of the errors with respect to nodal degrees of freedom (6 ‘} . If we define

then

(32) Using Eq. (30) with 6; = F, u, u, r,,, r,, and r,, and Eqs. (28), the required expression in Eq. (323 can be easily obtained:

(33)

(34)

I$$)‘=[PI. - [%[-g] +(~~~){~}t]? - [?“[$I +(~+~){~}‘1,,Ol,PI, INI] 7

(37)

(38)

Partial derivatives (33>-(38), can be derivatives of [N] Partial derivatives Eq. (29). Details

of the approximation functions [N] with respect to r and z, as needed in Eqs, easily obtained by using the Jacobian matrix of transformation [rz/&] and the partial with respect to 5 and 3. This procedure is well known, hence the details are omitted. of r), with respect to u and u for the power law model can be easily obtained by using are given in the following.

N.B.

6.2.

Edgar,

KS.

Surana

I Comput. Methods Appl. Mech. Engrg. 113 (1994) 271-300

283

Power law model for vu

2&

(2+)(z;)(-{-.&})

(39)

{-&} =e

(!+)(zy-w~{-&})

(40)

(-&}

where

(41) (42) The errors given by Eqs. (28) and the derivatives of the errors with respect to {ae} given by Eqs. (33)-(38) are used to compute the vector {g’} (Eq. (4)) and the Hessian matrix [H’] (Eq. (17)). The numerical values of the integrals are calculated using Gaussian quadrature with the (2p + 1) rule (exact integration) in both 6 and 77 directions where p is the order of the polynomial for the element approximation. A numerical integration procedure using a lesser number of integration points than those required for the exact integration is termed reduced integration. The order of Gaussian quadrature required for exact integration is determined based on the highest degree of the polynomials appearing in the integrands of the coefficients of [H’], {g’} and I’ without considering the terms resulting from the non-Newtonian viscosity 7, because 7, is a complicated function (as opposed to a simple algebraic polynomial) of the approximation functions. However, we have found that the (2p + 1) integration rule, which is exact for all other terms except those containing vu, is satisfactory for 0.15 s IZs 1.5. For 0.15 s n s 1.5, a higher order quadrature produced absolutely no change in the computed results. The system of algebraic equations in Eq. (15) is solved using wavefront method. The CPU time required in this method is a function of maximum wavefront and rms wavefront. We have used this procedure successfully to solve problems with maximum wavefront in excess of 3000 degrees of freedom. However, this solution method may not be the most efficient (in terms of storage and CPU time) method of solving Eq. (15) for meshes containing a large number of degrees of freedom. Here, we elected to use this method because it was easily available from other finite element subsystems in our finite element software system. The efficiency aspects of the formulation and solution methods are neither claimed nor discussed in this paper. These aspects are currently being investigated and will be presented in a separate paper.

7. Numerical

examples

Numerical examples are presented in this section to demonstrate the accuracy and convergence characteristics of the present p-version least squares formulation for constitutive behavior described by a power law model. The first example, annular couette flow, illustrates the accuracy of the formulation for crude meshes and power law index as low as 0.25. Discussion on mesh construction and the manner in which the mesh refinement is carried out is given. Also, the aspects of reduced versus exact integration are discussed. In Example 2, annular driven cavity, we also present a discussion on mesh construction and mesh refinement as well as a study of reduced versus exact integration. Example 3 presents results for a 2: 1 axisymmetric sudden expansion. Examples 2 and 3 demonstrate the ability of the formulation to accurately predict the behavior of the flows involving recirculations. For each of the three examples, comparisons with analytical, numerical and experimental results reported in the technical literature are given where appropriate. In the numerical examples, the exact integration rule is used to compute results unless stated

284

N.B.

Edgar,

K.S. Surana

I Camput. Methods Appl. Mech. Engrg.

otherwise. All numerical computations (64 bit word).

were performed

11.1 (1994) 271-30~~

on an HP 730 using double precision arithmetic

7.1. Example 1. Couette shear flow In this example, fully developed couette flow in an infinitely long annulus with a prescribed pressure gradient is investigated for power law flow indices ranging from 0.25 to 1.5. Fig. 2 shows a schematic of the two concentric tube configuration. The outer tube is stationary and the inner tube is moving in the negative z direction with a constant velocity of 3.0. The flow in the annulus between the tubes is fully developed and is pressure driven. Figs. 3(a) and (b) show a five element uniform mesh and a three element graded mesh with the specified boundary conditions. For n = 1.0 (Newtonian flow), the analytical solution can be obtained for the axial velocity profile and is given by V

v=ln(k)+

(PO-P,)R2 4mL

k2-1 o(-i;;ii;i-)1n(~)w(W)(1+(&)2)’

(43)

where k = Ri /R, is the annular geometry factor, y is the axial velocity of the inside tube, P,, P, are the pressure values at fore and aft ends of the annulus and L is the annulus length. Since the flow is fully developed and the quantities of interest are not varying in the axial direction, the p-level in the ,$’direction (axial direction) is chosen to be two and is kept fixed. The p-level in the 7 direction is progressively increased until the total error functional is below a prespecified tolerance. In designing a finite element mesh, the ultimate goal is to arrive at a mesh which yields most accurate results for the fewest total number of degrees of freedom. In context with the present formulation, we are looking for the lowest error functional value for the least number of total degrees of freedom. The

.

Ir

.lJJ P= 100

P=85

u=o

u=o

+

0.2

--

Ro

u=o

-

vi= -3

-

0.2 0.2 0.2 0.2

(a) Fig. 2. Schematic and boundary conditions for pressure driven annular couette flow (Example 1). Fig. 3. Finite element meshes for pressure driven annular couette flow (Example 1): (a) a five element uniform mesh; (b) a three element graded mesh.

N.B. Edgar, KS. Surana I Comput. Methods Appl. Mech. Engrg. 113 (1994) 271-300

285

uniform meshes, often lead to unnecessary degrees of freedom in many areas of the domain. The graded meshes on the other hand concentrate additional degrees of freedom in the areas of singularities where they are needed and thus result in meshes with substantially lower total degrees of freedom when compared with uniform meshes for a given accuracy. The LSFEF provides a ‘built in’ element error functional Z, (for an element e) which is an indication of the error in the solution (or accuracy) and can be used to guide the mesh refinement process. Using the element error functional I’, we define the mean square element error functional I’,,,,

where A’ is the area of element e. An optimally graded mesh would yield uniform values of II,, for all elements of the discretization. The finite elements of the mesh with values of Zks, greater than a threshold value (maybe preset or dynamically computed each time) are likely candidates for mesh refinement. The mesh design strategy used in this example is presented in the following. We start with an initial crude mesh consisting of three elements (uniform mesh) and obtain a solution for this mesh using p, = 2 and ps = 2 (fixed). The ZLs, values for the elements next to the walls were substantially higher (>O(lO-*)) than for the element in the middle indicating a need for mesh refinement in these two elements if p-levels are to be kept fixed. Uniform mesh refinement (elements were divided in half) was carried out in these two elements next to the walls and the solution was recomputed for p, = 2. For this mesh the I’,,, values were below 0( lo-*) for all elements of mesh indicating no need for further mesh refinement. The p-levels in the 7 direction were then increased for this graded mesh of five elements until the error functional Z and each component of {g} were below a convergence threshold value of lop4 or less. This mesh design procedure must be repeated for each value of the power law index. An optimally designed mesh for one value of ‘n’ may not be adequate or may even be overly refined for some other value of IZ. Results were also computed for uniformly refined meshes (for example, the five element uniform mesh of Fig. 3(b)). In general, as the power law index decreases, the performance of uniformly refined meshes begins to diminish. For values of at 3 1.0, uniform meshes performed well whereas, for IZ< 1.0, mesh grading was needed near the walls. Thus it is obvious that for 0.25 < n 6 1.5, there is no single optimal mesh. A five element uniform mesh produced good results for all values of it in the range 0.25 G IZ< 1.5 (even though it may be optimal for only one value of n or a much narrower range of n). The construction of a five element graded mesh prompted us to investigate the three element graded mesh of Fig. 3(a) where the elements near the walls are much smaller in size. The three element mesh of Fig. 3(a) yielded the same degree of accuracy as the five element uniform mesh, but for far fewer degrees of freedom. For all of the results reported for this example no more than lo-13 iterations were necessary for convergence tolerance O(10m4) for both {g} and I. Fig. 4(a) shows p-convergence of the error functional for various values of II for the three element graded mesh for Fig. 3(a). The plots of converged axial velocity for various values of it are shown in Figs. 5(a) and (b). Fig. 5(a) also shows a comparison between our computed conveyed axial velocity with the analytical solution for IZ= 1.0. For a power law index other than 1.0, a closed form analytical solution was not available. However, Hanks and Larson [42] provide values of the location of maximum velocity for various values of II and the annular geometry factor k. Table 1 shows a comparison between the results obtained from the present formulation and those reported by Hanks and Larson [42]. Another point of interest is the issue of reduced integration which refers to a lower order quadrature than that needed to exactly integrate the coefficients of the element matrices. The use of reduced integration in context with the least squares finite element formulation has been advocated by [14] as a means of improving accuracy. Even in more recent work dealing with LSFEF for CFD, reduced integration has been used for integrating the coefficients of element matrices. It was pointed out by Winterscheidt and Surana [37,38] (and possibly by others) in connection with p-version LSFEF for the convection diffusion problem and Burgers’ equation that the use of reduced integration in LSFEF leads to collocation which requires extremely refined meshes for accurate results. It was also pointed out that for coarse meshes, p-convergence of the error functional Z (recomputed using exact integration) is non-monotonic when reduced integration is employed and thus adaptive processes are not possible when reduced integration is used. In this example, we present reduced integration results for the three

N.B.

286

Edgar,

K.S. Sumna

i Comput.

lIx4mpla 1 lknllui4r cou4tt4 Thr44 Elamsnt m4ot

Methods Appl. Mech. Engrg.

Exrmpla 1 AnntIlar couatta Flow Thr44 Elamant Graded Yeah fft”dy$ IntegrationI

Flow

Graded

113 (1994) 271-300

Mad

Int4@4tionl

l

r(y02)25 QP,=2

0:ao A n = 0.75 0 P = 1.00 Q xl = 1.50

P, = 2

0 Xl =

-

---

_

XII

0 n l P 8 n 6, P

= = = = =

0.25 0.50 0.75 1.00 1.50

4

Pg = 5

Pq = 5

10’

10 =

10’

Dagreea of Freedom Fig. 4.

(a) p-Convergence of error fu~ctionai I, exact inte~ation integration (Example 1).

10’ Degress

(Example 1). (b) p-Convergence

10’ of Freedom

of error functional I, reduced

graded mesh of Fig, 3(a) which gives excellent results for all of values of n when integration rule is used. We choose a wide range of IZ values (0.25, 0.75 and 1.5). Fig. 4(b) ~-convergence of the error functional for reduced integration. We note that for II = 0.25 For values of YEclose to 1.0, reduced integration p-convergence of Z is non-monotonic. reasonable results but at the expense of twice as many iterations. As n deviated from ‘l’, element mesh with reduced integration produced spurious results indicating the need elements. element

7.2. Exa~~i~

2. Annular

10 d

the exact shows the and 0.75, produced the three for more

drive cavity

In this example, we consider an annular driven cavity for power law indices ranging from 0.15 to 1.50. Fig. 6 shows a schematic of the cavity and the boundary conditions. Since there are no published results for this problem, a two-dimensional cavity was simulated by making R, very large. For simulating two-dimensionai driven cavity, we consider Ri = 25, R, = 26 for Reynolds number of 100 and n = 1.0. For this do-dimensional cavity at Reynolds number 100, published numerical results of Ghia et al. [43] are available for comparison. We also consider this two-dimensional configuration at IZ= 1 and Re, = 100 for mesh design purposes. As demonstrated in Example 1, a single mesh will not perform optimally for all values of n for a fixed Reynolds number. Our strategy here would be to design a good mesh for this two-dimensional cavity at IZ= 1 and Re, = 100 and then use it for other values of n. If this mesh becomes inadequate for some values of n (which is expected) then additional mesh refinement can be performed based on Zks, values. An initial threshold of 0(10e2) is used for Zks, with a p-level of 2 ( pr = p, = 2) in designing the mesh. First we consider meshes. An initial uniform mesh of 16 elements (Fig. 7(a)) with p-level of 2 was used to compute the results.

N.B. Edgar, K.S. Surana f Comput. Methods Appl. Mech. Engrg. 113 (1994) 271-300

287

Example 1 Annular couette Flow Three Fknent Graded Meah LEsaet Intqrationl

eolu~on (xl = 1.0: --__

.!j 0.80

t t

z

0.0

-2.0

_----

-4.0 -8.0 Axial velocity v

n = 11= n = p =

1.50 1.00 0.76 0.50

-8.0

-10.0

0.80

* a B

Example 1 Annular CouetfA Flow Three Element Graded Yeah mad rntqrationl

0.80

ii! .--_

(npl==o$sP* = 5)

.

Fig. 5. (a) Converged axial velocity profiles for n = OS, 0.75, 1.0 and 1.5 (Example 1). (b) Converged axial velocity profile for n = 0.25 (Example 1).

Table 1 Radial location of maximum velocity for n = 0.3-1.0 and a comparison with Ref. [42] (Example 1) Power index of n

Radial distance to maximum velocity Ref. [42]

Present

0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.4991 0.5100 0.5189 OS262 OS324 0.5377 0.5422 0.5461

0.4990 0.5100 0.5190 OS260 0.53258 0.5375 0.5420 0.5460

1.0

N.B.

288

Edgar,

KS.

Surana

I Comput. Methods Appt. Me&. Engrg. 113 (1994) 271-300

T r

u=o.v=-1

u=c v=c

P=O-

. z

Fig. 6. Schematic and boundary conditions for annular driven cavity (Example 2).

I fl 1

0.25

1

0.25

[

0.25

f

0.25

i

I

1

I fl

i

I

17 -

-f-i

.-

0.25

j

k 0.25

,-

1

0.25

0.25

-4

/4-1.0

-

f=0.1,i=0.4,j=0.2,k=0.3,1=0.4 Cc)

T 1 1.0

e e

e e e e

IfIg

h

1

h

I 8 If1 -T -

T

-

h -

1.0

e e e e -

8

h g T -

1

l------*.0

------+

f=0.1,g=0.15,h=0.25 Cd)

Fig. 7. Finite element meshes for annular driven cavity (Example 2): (a) a 16 element uniform mesh; (b) a 100 element uniform mesh; (c) a 16 element graded mesh; (d) a 36 element graded mesh.

N.B.

Edgar,

K.S. Surana

I Comput.

Methods Appl. Mech. Engrg.

113 (1994) 271-300

289

For this mesh, the I’,,, values for the majority of the elements was 0(10m2) or lower except for the elements near the singularities (top two corners of the cavity) indicating a need for further mesh refinement in those areas. Successively more refined meshes were then constructed and for each mesh, Zks, values were examined until all elements of the mesh exhibited I’,,, values of O(lO-*) or lower. Fig. 7(b) shows the final 100 element uniform mesh arrived at by using this process. For this mesh, the p-levels were uniformly increased ( pr = p, = p for all elements of the mesh) until Z and {g} were below 26.0

25.6

.A

-

Example 2 Driven Cavity;

/’

mesh c R, = 25, R0 = 26

&I 25.6 -

Re, = 100, n = 1.0

a a

----

present. pc = p = 5 0 Ghia, Ghia,

ma shin t431

d

a

25.4 -

0.4

0.0

-0.2 -0.4 -0.6 Axial Velocity v et 2 = 0.5

-0.6

-1.0

Fig. 8. Converged axial velocity at z = 0.5 for large value of Ri (=25) and a comparison with Ref. (431 (Example 2).

Example 2 Driven Cavity; mesh (R, = 1. R, = 2)

p=2 X 0 A !I 0 x

n n n n n

= = = = =

c

0.25 0.50 0.75 1.00 1.50

p=5

2.6

Fig. 9. p-Convergence

3.0 3.3 LOG (degreea of freedom)

3.5

3.7

of the error functional I for n = 0.25-1.5 (Example 2)

a prespecified tolerance of 10T4. We note that for this 100 element uniform, mesh, the element size is 0.1 indicating that this is possibly the maximum element size capable of arresting singularities at the two corners of the cavity. With this addition& information, a much coarser graded mesh can be easily designed without much effort. The element sizes at the top two corners are kept at 0.2, The remaining portion of the cavity is uaiformfy divided into two columns of elements in the axial direction while generating three graded rows of elements in the radial direction as shown in Fig. 7(c). For this graded mesh of Fig. 7(c), p-levels were uniformly increased until the error functional I and the components of (g) were O(lO-“) or lower. Figure 8 shows a plot of the axial velocity u at z = 0.5 for the converged solution which is in excetXent agreement with the resu&s reported by Ebia et al. f43fS The r~m~i~~ng resu&s fur this driven catity problem are computed and presented for R, = 1, R, = 2, Re, = 100 and 0.15 d n G 1.5. First we use the 16 element graded mesh of Fig. 7(c) for various values of IZ and increase p-levels uniformly until convergence is achieved (I and {sj are 0(10S4) or lower). The 16 element graded mesh proved quite satisfactory for all values of n except 0.15. Fig, 9 shows p-convergence of the error functional for 0.25 c n .-E 1.5, For ra = 0.15, futiher mesh refinement was

= 1.00, meah c. p n = O.T& mesh c. p = 5 n = 0.50, mesh c. p = 8

n

u”,“______

--

----mm”m..*-~-

n

=

----

n

5

0.25, meah c. p = Q b.f5, mesh d, p = 4

-t

f.00 0.4

0.2

0.0

tiaz

-&2 vehG*

Fig. TO. Converged axial velocities u

-0.4 v

-0.8

-0.6

-LO

at. z - 0.5

at z = 0.5 for n = 0,15-1.5 (Example 2).

0.80

tn

0.1%

; II k %

0.06

d

,x ‘ij 2

ijry

-0.06

4

-0.18

-0.30 1.0

033

0.7

0.5

Axial

(34

0.2

Disbnce z

Fig. II. Conveqerl radial velocitias u at r = 1.5 for n =:0.15-1.5 (Exampla 2).

0.0

N.B.

Edgar,

K.S. Surana

291

I Comput. Methods Appl. Mech. Engrg. 113 (1994) 271-300

necessary to lower the error functional I. Fig. 7(d) shows a 36 element graded mesh. In this mesh also, the elements at the top two corners of the cavity are of size 0.1 and the remaining part of the cavity is divided uniformly into four columns of elements in the axial direction while generating five graded rows in the radial direction. For this graded mesh, p-levels were uniformly increased until I and {g} were O(10e4) or lower. A p-level of six was needed to accomplish this. Fig. 10 shows graphs of converged axial velocity u at z = 0.5 versus radius r for various values of n. Plots of converged radial velocity u at r = 1.5 versus axial distance z are shown in Fig. 11 for various values of n. Figs. 12 and 13 show some sample graphs of converged pressure at z = 0.5 and r = 1.5 for various values of 12. Sample plots of the streamlines for power law indices of 1.5, 1.0 and 0.25 are shown in Fig. 14. From these streamline plots we note that as IZ decreases, the location of the maximum intensity moves closer to the lid and to the right of the cavity while the recirculation zones in the bottom corners progressively become weaker (and will eventually cease to exist for very low values of n). For all values of n in the range 0.25 s IZc 1.5, Newton’s method converged in less than lo-13 iterations for a convergence tolerance of 0(10e4) for 2.00

1.80

1.60

1.40

n = 0.50, mesh c, p = 5 n = 0.25, mesh c. p = 5 n = 0.15, mesh d, p = 6

-------__._-..

1.20

--------

-0.04

0.01

Fig. 12.

-0.19

Converged pressures P at z = 0.5 for n = 0.15-1.5 (Example 2).

r

I Example 2: Driven Cavity z P n P P P

1.0

-0.14

-0.09 Pressure P at z = 0.5

0.8

0.7

0.5 Axial Dbtnnce

==‘;,~ ‘= = = = =

1.00, 0.75, 0.50, 0.25, 0.15,

0.4

= 2) meah meah meah meah meah

c, c, c, c, d,

p p p p p

= = = = =

5 5 5 5 tl

0.2

2

Fig. 13. Converged pressures P at r = 1.5 for n = 0.15-1.0 (Example 2).

292

N.B. Edgar,

l

**

x rr,

=

=

K.S. Surana

I Comput. Methods Appl. Mech. Engrg.

-0.03427 0.55@2E-06

Re. = 100,n = 0.25

113 (1994) 271-300

m *A

= -0.1364

x *_

=

Re, =

100,n = 1.00

0.21343-04

\

x *-

Re,

=

=

0.90373-04

100,n = 1.50

Fig. 14. Streamline plots for n = 0.25, II = 1.0, and n = 1.5 (Example 2).

both Z and {g}. However, for n = 0.15, mesh d was necessary (at a p-level of 5) and the procedure required 34 iterations for a convergence tolerance of O(10-4). As in Example 1, here we also investigate the use of reduced integration. For this study we only consider the 16 element graded mesh (mesh c) of Fig. 7(c) and the 100 element uniform mesh (mesh b) of Fig. 7(b). Figs. 15 and 16 show p-convergence of the error functional for reduced integration for these two meshes and a comparison with the exact integration (for some sample values of TI). Even

N.B.

Edgar,

KS.

Surana

I Comput. Methods Appl. Mech. Engrg.

1.0 [

0.0

Exemple 2: Driven Cavity 16 element 2reded meeh 0 n = 0.75 exe& inte2ration (2p+1) an= 0.75 reduced intqration (p+l) an= 0.25 exact inte2retton (2p+l) n n= 0.25 reduced intqration (p+l)

-

-1.0

-

-2.0

-

-3.0

-

-4.0

-

113 (1994) 271-300

= i

=5

’ ’

-5.0

2.0

Fig. 15. p-Convergence

’

’ ’

’

2.0

’

’

’

’

3.0 LOG (de2reee

’

’

’

’

of the error functional I, reduced integration, 1.0

0.0

0 0 B l

-

-1.0

-

-2.0

-

-3.0

-

-4.0

-

’

’

3.3 of freedom)

’

’

’

’

3.5

3.7

16 element graded mesh (Example 2).

Exemple 2: Driven Cavity 100 element uniform meeh n = 0.75, exact inteqetion ( +l) n = 0.7s reduced inta2retion “ipp+l) n = 0.25 exect inte2retion’ (2p+l) n = 0.25 reduced intepatoin (p+l)

t: s”

’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’ ’

-5.0

3.3

Fig. 16. p-Convergence

3.5

3.7 3.0 4.1 LOG (de2reee of freedom)

4.3

4.5

of the error functional I, 100 element uniform mesh (Example 2).

293

294

N.B. Edgar, KS. Surana i Comput.MethodsAppl. Mech. Engrg. 113 (1994) 271-300

though graphs of Z versus degrees of freedom for reduced integration are monotonic in this case, the error functional values obtained with reduced integration are much higher than those obtained with exact integration. Fig. 17 shows plots of the axial velocity at z = 0.5 versus radius r. As discussed in Example 1, reduced integration will perform well for very refined meshes, this is true in this case as well. From Fig. 17 we note that the solution for the 100 element uniform mesh with reduced integration matches the solution of mesh c (16 element graded mesh) with exact integration perfectly. However, the solution obtained from the 16 element graded mesh with reduced integration is grossly in error. 7.3. Example 3. A 2: 1 ax~yrnrne~~~ sudden expansion A 2:l (Z?,/R,) axisymmetric sudden expansion is investigated in this example for power law index ranging from 0.5 to 1.5 for a Reynolds number of 10. Fig. 18(a) shows the schematic of the problem, and the boundary conditions. A 24 element graded finite element mesh is shown in Fig. 18(b). The fully developed velocity profiles were computed for each power law index considered here and were imposed at the inlet. The p-levels in the 5 and n directions for each element of the mesh were uniformly increased until the error functional, Z, and each component of {g} were 0(10m4) or less. The graded 2.00

Exa~~pls 2: Driven Cavity

1.00

0.4

0.6

0.2

-0.4 0.0 -0.2 Axial Velocity P at a = 0.5

-0.6

-0.8

-1.0

Fig. 17. Converged axial velocity at z = 0.5 for n = 1.0, exact and reduced integration (Example 2). 1

r u=v=o a

P=O U=O

2 u=o “=* .

7h 0

+=%I!! a

l

.

-

40

A-+

Z u=o Tn=

0

* specified inlet velocity (a)

I

a

lblbl

e

I

a=0.9

f

I

b=O.l

g

e=0.4

I

h

f= 1.5 g= IO.0 h=28.0 @>

Fig. 18. A 2: 1 axisymmetric sudden expansion (Example 3). (a) Schematic and boundary conditions for 2: 1 axisymmetric sudden expansion; (b) a 24 element graded mesh.

I

N.B.

Edgar,

K.S.

Surana

I Comput.

Methods

Appl.

Mech.

Engrg.

113 (1994)

271-300

295

finite element mesh presented in Fig. 18(b) was arrived at using the procedures similar to those described for Examples 1 and 2. Initially a coarse mesh was constructed and Z&, values were used to guide the mesh refinement process. The number of iterations needed for convergence were 13 in all cases reported here except n = 1.0 which required 11 iterations. As a check on the accuracy of our imputations, the reattachment length was computed for Newtonian fluid flow (n = 1.0) for Reynolds numbers of 25, 50,7.5 and 100 using the mesh of Fig. 18(b).

J

8.00

t

2

3.00

/ -

Example 3 Ashymm6tric Suddan &par&on Newtonian Fluid (n = 1.0) 0 MacaSno at.aI.. erp6r. A Soott and Nina. FEN * Napolitano &al. FDM 0 pre6ent formulation Scott and Nina, flttad IXUW

0.00 25.0

0.0

50.0 hyIiOld6

Fig. 19. Reattachment

75.0 Numb6r i&e,

100.0

125.0

length versus Reynolds number (Example 3).

-2.00

-2.40

-2.60 3 ‘;?’ s -3.20

-3.60

fl

-4.00 2.8

Fig. 20. p-Convergence

3.7 3.1 3.4 LOc(degtse6 of freedom)

4.0

of the error functional 2 for n =0.5, 0.75, 1.0 and 1.5 (Example 3).

296

N.B.

Edgar,

KS.

Surana

I Comput. Methods Appl. Mech. Engrg.

113 (1994) 271-300

The Reynolds number was based on a unit average inlet velocity and an inlet radius of 1.0. The values of reattachment length obtained from the present formulation is compared with the experimental results of Macagno and Hung [44], the finite element numerical prediction of Scott and Mirza [45], and the numerical finite difference solution of Napolitano et al. [46]. From Fig. 19, we note that our results show excellent agreement with both experimental and numerical results reported in the literature. Fig. 20 shows p-convergence of the error functional for various values of it. Figs. 21 and 22 show plots of the converged axial velocity profiles at various axial stations along the length of the expansion for n = 1.5 and IZ= 0.5. Plots of the wall shear stress at r = R, versus axial length z for various values of IZ are shown in Fig. 23. The transition of 7rz from negative to positive indicates flow reattachment and the end of the recirculation zone. It may be noted that as n decreases, the reattachment length also decreases. Figs. 24 and 25 show graphs of the converged pressure at r = R, versus axial length and the converged pressure at the centerline versus axial length for various values of n. 2.00

1.60

*

-_.-_. ---_._.--. ________ _---

1.20

-_-

ij p:

0

Sudden Expandon) Fk, = 10, n = 1.5 PI = P, = 5 fzYcd 2 = -1.2 z = -2.2 z = -3.0 Z = -4.2 outlet

0.80

0.40

0.00 0.25

-0.25

-0.75 Md

-1.25 Velooity v

-1.75

-2.25

Fig. 21. Converged axial velocity profiles at various axial locations for n = 1.5 (Example 3). 2.00 Example 3 (Msymmetric Sudden Expansion) Rem = 10, P = 0.5 PI = p_ = 5

1.60

*

m-q. -.-_--.. ____--__ m-m ---.

1.20

9 a e a

0 UldJtiCd Inlet z =

-1.2 z = -2.2 z = -3.9 z = -6.0 outlet

0.80

0.40

0.00 0.25

-0.25

-0.75 -1.25 Axial Velodty v

-1.75

-2.25

Fig. 22. Converged axial velocity profiles at various axial locations for n = 0.5 (Example 3).

N. 3. Edgar,

KS.

Swam

i Comput. methods Appl. Mech. Engrg.

113 (1994) 271-300

297

0.10

0.08

d II

0.00

*

3

r 1 ii!

2

ii

0.04

0.02

0.00

-0.02

-0.04 -1.0

-2.5

-4.0 Axial Distance x

-5.5

-7.0

Fig. 23. Converged shear stress T,, at r = R, versus axial distance z for n = 0.5, 0.75, 1.0 and 1.5 (Exampfe 3). 0.00

-1.00

rr’ II L %

RI

-2.00

[ :: PI -3.00

-4.00 -1.0

-14.0 Axial Distance z

-27.0

-40.0

Fig. 24. Converged pressures P at r = R, versus axial distance t for n = 0.5, 0.75, 1.0 and 1.5 (Example 3).

8. Summary and concIusions An exact p-version least squares minimization procedure has been presented for coupled non-linear partial differential equations. In the development of the least squares error functional, actual non-linear partial differential equations are utilized without linearizing the non-linear terms or introducing any other approximations. This procedure has been applied to the specific case of axisymmetric, incompressible, non-Newtonian, steady state flow described by the Navier-Stokes equations. The power law model has been used to characterize the fluid viscosity. For non-linear problems, the least squares minimization principle requires finding a (8) which makes {g} (a vector of partial derivatives of the error functional with respect to the global nodal degrees of freedom) a null vector. This is accomplished by utilizing Newton’s method with a line search. Three numerical examples are presented to demonstrate the accuracy and convergence characteristics of the method. The following specific conclusions can be drawn based on this research.

298

N.B.

Edgar,

K.S. Surana

I Comput. Methods Appl. Mech. Engrg.

113 (1994) 271-300

Example 3 (Aqvmmetric Sudden Expamion)

-1.90 -

-3.50[ 0.0

’

’

’

’ -10.0

’

’

’

’

’

-20.0 AxiaJDirtancez

’

’

’ -30.0

’

’

’

1

-40.0

Fig. 25. Converged pressures P at r = 0 versus axial distance z for n = 0.5, 0.75, 1.0 and 1.5 (Example 3).

The least squares minimization procedure presented here is applicable to any set of linear or non-linear differential equations and will result in a true least squares functional free of assumptions and approximations. Linearization of the non-linear differential equations before constructing the least squares functional unnecessarily destroys the true least squares functional that corresponds to the actual problem at hand. The construction of the least squares functional and finding a solution vector (6) that satisfies the necessary and sufficient conditions for minimization are two unrelated issues. Assumptions or approximations employed in devising a solution procedure to find such a {S} have absolutely no effect on the exactness of the least squares functional to be minimized. The condition {g({S})} = (0) for a solution vector {S} only represents the necessary conditions for the global minima of I. The sufficient condition is given by {6*}‘[H( {6})]{S*} > 0 for an arbitrary vector {S *}. Since I = 0 is the minima we are seeking, a solution vector (8) which satisfies Z = 0 will automatically satisfy the sufficiency condition of the positive definiteness of the Hessian matrix [HI. Thus for the procedure presented here, the sufficient condition given by (6 *}‘[H( {a})] (6 *} > 0 can be replaced with Z = 0. Making an approximation in the computation of the Hessian matrix [H’] (neglecting the term containing the second derivative of I’ with respect to (6’)) only alters the search directions during the Newton’s iteration procedure and has absolutely no effect on the exactness of the least squares minimization functional. This approximation has beneficial effects of simplifying programming, speeding up computations and accelerating the convergence of the Newton’s method with line search. The least squares formulation automatically provides a measure of the solution error through the error functional. The element error functional values are very useful in adaptive mesh refinement or adaptive p-level changes. The total error functional is a monotonic function of the total degrees of freedom as the p-levels are increased for a fixed mesh. The numerical results are also presented for reduced integration (Examples 1 and 2) using the (p + 1) rule and are compared with those obtained using exact integration ((2~ + 1) rule). p-Convergence of the error functional Z is non-monotonic in the case of Example 1 for n = 0.25 and 0.75 when the three element graded mesh is used. In Example 2, even though the p-convergence of the error functional Z for reduced integration is monotonic (for the meshes considered), the actual numerical values of Z are much higher than those for the exact integration. Furthermore, we note that for the meshes for which exact integration produces converged and good results, reduced integration fails to do so. The numerical examples clearly demonstrate the effectiveness of the formulation and the solution procedure for a wide range of power law index in obtaining accurate and converged solutions. From

N.B.

Edgar,

K.S. Surana

I Comput.

Methods Appl. Mech. Engrg.

113 (1994) 271-300

299

p-convergence graphs, we note that the convergence rates are pretty much independent of the power law index. The computational results agree very well with both experimental measurements and other numerical solutions considered to be quite accurate. In summary, the p-version least squares formulation presented here produces excellent results, has good convergence characteristics and provides an accurate tool for numerical simulation of axisymmetric, incompressible, non-Newtonian steady state fluid flow.

Acknowledgments

The computing facilities provided by the Computational Mechanics Laboratory (CML) of the Department of Mechanical Engineering are gratefully acknowledged. The financial support provided to the first author by the Department of Mechanical Engineering through a graduate teaching assistantship and by the School of Engineering through the Strobe1 Fellowship is also acknowledged.

References [l] D. Winterscheidt and K.S. Surana, p-Version least squares finite element formulation for two dimensional incompressible fluid flow, Internat. J. Numer. Methods Fluids (1993) in press. [2] T.J.R. Hughes, W.K. Liu and A. Brooks, Finite element analysis of incompressible viscous flows by the penalty functional formulation, J. Comput. Phys. 30 (1979) l-60. [3] A.N. Brooks and T.J.R. Hughes, Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199-259. [4] T.J.R. Hughes, L.P. Franca and M. Balestra, A new finite element formulation for computational fluid dynamics: V. Circumventing the BabuSka-Brezzi condition: a stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations, Comput. Methods Appl. Mech. Engrg. 59 (1986) 85-99. [5] T.J.R. Hughes and L.P. Franca, A new finite element formulation for computational fluid dynamics: VII. The Stokes problem with various well-posed boundary conditions: symmetric formulations that converge for all velocity/pressure spaces, Comput. Methods Appl. Mech. Engrg. 65 (1987) 85-96. [6] C.H. Hirsch and G. Warzee, A finite element method for the axisymmetric flow computation in a turbomachine, Intemat. J. Numer. Methods Engrg. 10 (1976) 93-113. [7] K.R. Reddy and R.L. Tanner, Finite element solution of viscous jet flows with surface tension, Comput. Fluids 6 (1978) 83-91. [8] Bernard J. Omodei, On the die-swell of an axisymmetric Newtonian jet, Comput. & Fluids 8 (1980) 275-289. (91 P. Bar-Yoseph, J.J. Bleck and A. Solan, Finite element solution of the Navier-Stokes equations in rotating flows, Internat. J. Numer. Methods Engrg. 17 (1981) 1123-1146. [lo] J. Donea, S. Giuliani and K. Morgan, The significance of Chequerboarding in the Galerkin finite element solution of the Navier-Stokes equations, Intemat. J. Numer. Methods Engrg. 17 (1981) 790-795. [ll] G.F. Carey and R. Krishnan, Continuation techniques for a penalty approximation of the Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 48 (1985) 265-282. [12] P.P. Lynn and S.K. Arya, Use of the least squares criterion in the finite element formulation, Internat. J. Numer. Methods Engrg. 6 (1973) 75-88. [13] P.P. Lynn and S.K. Arya, Finite elements formulated by the weighted discrete least squares method, Internat. J. Numer. Methods Engrg. 8 (1974) 71-90. [14] O.C. Zienkiewicz and D.R.J. Owen, Least squares finite element for elasto-static problems, use of ‘reduced’ integration, Internat. J. Numer. Methods Engrg. 8 (1974) 341-358. [15] P.P. Lynn, Least squares finite element analysis of laminar boundary layer flows, Internat. J. Numer. Methods Engrg. 8 (1974) 865-876. 1161 E.D. Eason and C.D. Monte, Jr., Solution of non-linear boundary value problems by discrete least squares, Internat. J. Numer. Methods Engrg. 11 (1977) 641-652. [17] M.J.D. Powell, A method for minimizing a sum of squares of non-linear functions without calculating derivatives, Comput. J. 7 (1965) 303-307. [18] J.F. Polk and P.P. Lynn, A least squares finite element approach to unsteady gas dynamics, Internat. J. Numer. Methods Engrg. 12 (1978) 3-10. [19] G.J. Fix and M.D. Gunzburger, On least squares approximations to indefinite problems of the mixed type, Internat. J. Numer. Methods Engrg. 12 (1978) 453-469.

300

N.B.

[ZO] C.A.J. Fletcher,

Edgar,

K.S. Surana

I Comput. Methods Appl. Mech. Engrg. 113 (1994) 271-300

A primitive variable finite element formulation for inviscid, compressible flow, J. Comput. Phys. 33 (1979) 301-312. [21] M.O. Bristeau, 0. Pironneau, R. Glowinski, J. Periaux and P. Perrier, On the numerical solution of non-linear problems in fluid dynamics by least squares and finite element methods (1) least squares formulation and conjugate gradient solution of the continuous problems, Comput. Methods Appl. Mech. Engrg. 17118 (1979) 619-657. [22] H. Nguyen and J. Reynen, A space-time least-squares finite element scheme for advection-diffusion equations, Comput. Methods Appl. Mech. Engrg. 42 (1984) 331-342. [23] H. Nguyen and J. Reynen, A space-time finite element approach to Burgers’ equation, in: C. Taylor, E. Hinton and D.R.J. Owen, eds., Numerical Methods for Non-Linear Problems, Vol. 2 (Pineridge Press, Swansea, UK, 1984) 718-728. 1241 B. Tabarrok and M. Ziad Saghir, A new mixed formulation for 2D incompressible flow, Comput. Methods Appl. Mech. Engrg. 43 (1984) 81-102. [25] A.K. Aziz, R.B. Kellogg and A.B. Stephens, Least squares method for elliptic systems, Math. Comp. 44 (1985) 53-70. [26] B.N. Jiang and G.F. Carey, Adaptive refinement for least squares finite elements with element-by-element conjugate gradient solution, Internat. J. Numer. Methods Engrg. 24 (1987) 569-580. ]27] G.F. Carey and B.N. Jiang, Least squares finite element method and preconditioned conjugate gradient solution, Internat. J. Numer. Methods Engrg. 24 (1987) 1283-1296. [28] G.F. Carey and B.-N. Jiang, Nonlinar preconditioned conjugate gradient and least-squares finite elements, Comput. Methods Appi. Mech. Engrg. 62 (1987) 145-154. [29] B.N. Jiang and G.F. Carey, A stable least-squares finite element method for non-linear hyperbolic problems, Intemat. J. Numer. Methods Fluids 8 (1988) 933-942. [30] I. Kececioglu and B. Rubinsky, A mixed-variable continuously deforming finite element method for parabolic evolution problems, Part I: The variational formulation for single evolution equation, Internat. J. Numer. Methods Engrg. 28 (1989) 2583-2607. 1311 I. Kececioglu and B. Rubinsky, A mixed-variable continuously deforming finite element method for parabolic evolution problems. Part II: The coupled problem of phase change in porous media, Internat. J. Numer. Methods Engrg. 28 (1989) 2609-2634. [32] I. Kececioglu and B. Rubinsky, A mixed-variable continuously deforming finite element method for parabolic evolution problems. Part III: Numerical implementation and computational results, Internat. J. Numer. Methods Engrg. 28 (1989) 27152760. [33] B.-N. Jiang and C.-L. Chang, Least-squares finite elements for the Stokes problem, Comput. Methods Appl. Mech. Engrg. 78 (1990) 297-311. f34] B.-N. Jiang and L.A. Povinelli, Least-squares finite element method for fluid dynamics, Comput. Methods Appl. Mech. Engrg. 81 (1990) 13-37. [35] I. BabuSka, O.C. Zienkiewicz, J. Gago and E.R. de Al. Oliveira, Accuracy Estimates and Adaptive Refinements in Finite Element Computations (Wiley, New York, 1986). [36] B.N. Jiang and V. Sonnad, Least-squares solution of incompressible Navier-Stokes equations with the p-version of finite elements, NASA Technical Memorandum 105203, ICMP-91-14, 1991. 1371 D. Winterscheidt and K.S. Surana, p-Version least squares finite element formulation for convection-diffusion problems, Internat. J. Numer. Methods Engrg. 36 (1993) 111-133. [38] D. Winterscheidt and K.S. Surana, p-Version least squares finite element formulation of Burgers’ equation, Internat. J. Numer. Methods Engrg., in press, 1391 B.C. Bell and K.S. Surana, p-Version least squares finite element formulation for two-dimensional incompressible non-Newtonian isothermal and non-isothermal fluid flow, Internat. J. Numer. Methods Fluids, in press. [40] R.B. Bird, W.E. Stewart and E.N. Lightfoot, Transport Phenomena (Wiley, New York, 1960). [41] R.B. Bird, R.C. Armstrong and 0. Hassager, Dynamics of Polymetric Liquids, Vol. 1 (Wiley, New York, 1987). [42] R.W. Hanks and K.M. Larsen, Ind. Engrg. Chem., Fundam. 18 (1979) 33-35. [43] U. Ghia, K.N. Ghia and C.T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comput. Phys. 48 (1982) 387-411. [44] E.O. Macagno and T. Hung, Computational and experimental study of a captive annular eddy, J. Fluid Mech. 28 (1967) 43-64. 1451 P.S. Scott and F.A. Mirza, A finite element analysis of laminar flows through planar arid ~isymmetric abrupt expansions, Comput. & Fluids 14 (1986) 423-432. [46] M. Napolitano and P. Cinella, A numerical study of planar and axially-symmetric sudden expansion flows, Comput. & Fluids 17 (1989) 185-193.