An anisotropic h-adaptive finite element method for compressible Navier-Stokes equations

Computer methods in applied mechanics and englneering Comput. Methods Appl. Mech. Engrg. 146 (1997) 231-252

ELSEVIER

An anisotropic h-adaptive finite element method for compressible Navier-Stokes equations W. Rachowicz Section of Applied Mathematics, Technical University of Cracow, ~1. Warszawska 24, 31-155 Krakbw. Poland

Received 19 March 1996; revised 27 September 1996

Abstract

Application of the finite element method with significantly stretched elements to solve compressible Navier-Stokes equations is presented. Stretched elements are introduced in an adaptive process as a result of selective directional subsecting of elements. An adaptive strategy of such anisotropic h-type refinements was developed. It is based on minimization of the interpolation errors, reduction of the residual error indicators and an adaptive recovery of the viscous fluxes with the possibility to control their accuracy. An iterative solver with the convergence characteristics independent of the aspect ratio of elements is presented. The solver is the GMRES algorithm with the domain decomposition Schwarz-type preconditioner. The proposed techniques are used to solve a benchmark problem, a hypersonic high Reynolds number flow with strong shock-boundary layer interaction.

1. Introduction

Evaluation of viscous fluxes, the skin friction and the heat flux coefficients, is one of the most challenging tasks in the numerical analysis of compressible Navier-Stokes equations. ‘While obtaining a relatively accurate solution of the problem of a viscous flow inside the computational domain is of the same level of difficulty as solving an inviscid problem, an accurate simulation of viscous fluxes along the solid wall boundaries requires significantly finer approximation. This is caused by the conflict of scales: the viscous phenomena develop in very thin boundary layers where the solution is characterized by large gradients. The aerothermal loads, the viscous fluxes, are proportional to these gradients. As evaluation of the aerothermal loads is a primary goal of any simulation of a viscous flow (if not the only one), the accurate approximation of solutions inside boundary layers is an inseparable aspect of any analysis of a viscous flow. In this work we present an attempt to use an anisotropic h-adaptive finite element method [l] to solve the compressible Navier-Stokes equations. By ‘anisotropic h-adaptive’ we understand an approach in which the quadrilateral elements can be subsected in one of two possible directions. By performing a sequence of such refinements in elements approximating boundary layers we expect to approximate in the most optimal way the solutions which vary rapidly in the direction perpendicular to the boundary and are relatively smooth along the wall. We present an adaptive strategy which allows one to decide which elements should be subsected to reach a prescribed level of accuracy. The strategy is based on reducing an error indicator of derivatives of the solution and the requirement to reduce the contribution of artificial dissipation stabilizing the scheme. Since anisotropic refinements provide a possibility to generate a wider class of meshes than isotropic refinements (when elements are divided into four elements-sons) they allow one to approximate more optimally solutions not only in boundary layers but in the whole computational domain, that was demonstrated in [l]. We take advantage of this feature of the method developing an anisotropic adaptive strategy 00457825/97/$17.00 @ 1997 Elsevier Science S.A. All rights reserved SOO45-7825(96)01230-3

PII

232

W Rachowicz/Comput.

Methods Appl. Mech. Engrg. 146 (1997) 231-252

for improving the solution in the whole computational domain. The approach is directed to reduce the residual of the solution and, at the same time, to take advantage of the anisotropic refinements to better adjust the mesh to variation of the solution in different directions. Finally, we address the critical problem of solving the system of linear equations resulting from the finite element discretization with the significantly stretched elements. Such elements raise a common criticism as they result in ill-conditioned systems which are hard if not impossible to solve with iterative solvers. As a remedy of this problem we present a version of the overlapping domain decomposition preconditioner introduced in [2] in the context of symmetric elliptic equations discretized with stretched elements. We show that the performance of the preconditioner combined with the GMRES algorithm is uneffected by stretched elements in the case of flow problems formulations. The techniques proposed in the work are illustrated with a solution of a benchmark problem, a hypersonic flow with large Reynolds number.

2. Compressible

Navier-Stokes

We consider compressible U,, + Fj(U).i

equations

and the finite element

Navier-Stokes

method

equations in conservative

discretization

variables:

= FY(U, VU),;

(1)

In (1) and throughout the paper commas denote differentiation with respect to time and the spatial coordinates, (.),, = 3/dt, (.),; = d/dXi, and the Einstein summation convention over repeated indices is used. U = {Q,ml, m2, e}T denotes a vector of conservative variables, density Q, momentum per unit volume mi and the total energy per unit volume e. J’,(U) are the Eulerian fluxes: P; = {m;, nrrm;/Q + &;p,

%mi/e+

&;p,

(e

+p)W/@}T,

i

=

1,2

(2)

where p denotes the thermodynamical pressure. p = (y - 1)~. with L = e - 1/2(mt + rn:)/e being the internal energy per unit volume. y the ratio of the specific heats at constant pressure C, and at constant volume C,, y = 1.4. Fy(U, VU),I’ = 1,2 are the viscous fluxes: Fy = (0, TIN,TZ~,Ti?;uj + LJ;}~~

(3)

qi = -KO,i iS the heat flux, Lli = mile where rij are the viscous stresses, rij = /.L(u,, + Uj,; - Z!/LJSijUk,k), denotes the velocities and 0 = L/(&,) is the temperature. The viscosity coefficient p and the heat conductivity K are functions of the temperature 0 (the Sutherland law). As viscous fluxes are linear functions of VU they can be expressed as Fy = Kij(U)U,j, with Kij = aFy/aU,j. We discretize compressible Navier-Stokes equations using the version of the SUPG method for the formulation in conservation variables [3,4]. We assume that the computational domain R is divided into a mesh of quadrilateral elements K. The mesh defines the standard finite element space of bilinear shape functions Vh and we assume that any u E Vh satisfies homogeneous Dirichlet boundary conditions. We look for the solution U E Vh + ug such that:

J 0

J J +AiU,i -(Kiju.j),i] CI.X =J J~aw:Ak[U,r VdWTU.i dX

WT(U., + Ail-J.0 dx + II WT(KijU,j) dx + C

K

K

WTHdS,

+c K

VW E Vh

no,

K

where u 0 and H are the data specifying the inhomogeneous Dirichlet and Neumann boundary conditions, the latter being prescribed on 80~ c 6’0. A; are the Eulerian jacobians, A, = dFi/dU and coefficients r, and vd are given as (AU.ilTA,‘AjU.j

mW=l,2(hIPil) To

=

2(c+ju

.pI)

’

(S,.iU,i)TAi’ b.jU,j

(5)

W Rachowicz/Comput.

where A0 is the symmetrizer

of the Navier-Stokes

with H = -Q ln(L/eY) being the entropy

233

Methods Appl. Mech. Engrg. 145 (1997) 231-32

equations.

, c In (5) p is . a unit vector parallel to V ]]U ]]11,,,

function.

is the speed of sound, and h, are sizes of an element, h; = (~f&,(FJx;j~3&)‘)‘~‘.

The map x,(<,) is the

transformation of the reference element l? = [-1, l]* into the actual element K. and &(x,) in (5) is its inverse. In (4) the first two integrals correspond to the Galerkin formulation of the problem while the remaining integrals are the SUPG stabilizing terms. For the reasons which we discuss later, we simplify the second of the stabilizing integrals by dropping out U,, term. Then we discretize (4) in time using a version of the Euler method obtaining the following time integration scheme: WTU"+'

= / . f‘J

dx +

At + JW~(Kij J J

TaAiAj + vdS;j)Uy’ dx

f2

WTA;U;;. dx + At

WTU”dx-At

R

TaWlAk(K;jUn),; dx + At

R

J

WTH dS

212,

where At is the time step, U” E U(k At) and coefficients Ai, Kij, vd and r, are evaluated for U”. A similar formulation can be obtained by a modification of the Taylor-Galerkin method as described for instance in [15]. In this approach the contribution WTAiAjUy’ in (7) is proportional to At2/2 instead of At T. Such a choice makes the scheme second-order accurate in time. On the other hand, however, if the steady state solutions are sought (U”+’ = U”), the factor At2/2 makes the stabilization term At-dependent, and. practically, absent for small At. Formulation (7) does not suffer from this drawback. The motivation for neglecting U., in the stabilizing integral in (4) is its influence on the convergence properties of the iterative solver which we discuss in detail in Section 4. Note that taking into account U,, in deriving (7) would result in the presence of the first-order term J,, r0W$ Un+’ d\: which would dominate the second-order terms (proportional to At) especially in the case of h-adaptive meshes, when the global value At may be much smaller then the local value of 7,. Such a strong dominance of unsymmetric terms in the formulation requires a special treatment when using iterative algorithms to solve the discrete equations. Namely, it requires the presence of a global, involving all degrees-of-freedom. step in the iterations [14], which may be not an optimal approach (see discussion in Section 4). On the other hand. neglecting U,l in the stabilizing part of (4) does not affect steady state solutions. An alternative version of the SUPG method based on the entropy variables leads to a symmetric formulation. Choosing the version based on conservation variables we follow the suggestions expressed in [3] and [4] that the results obtained with the two methods are practically identical while the formulation using conservation variables is essentially simpler. 3. An anisotropic h-adaptive

strategy

Establishing the technique of modifying meshes which would lead to achieving an accurate solution at the possibly lowest computational cost is a difficult task. The most popular strategy consists in refining the elements with the biggest values of a certain kind of an error indicator. Such a technique is based on the intuitively obvious condition that obtaining an accurate solution requires elimination of elements with large errors. A quantitative justification of this refinement strategy is based on the assumption that the error associated with an element depends on the element size as f(x)h” with f(x) being the solution dependent function, a = const, and that we want to minimize the sum of such errors with the constraint that the number of elements is fixed [5]. Despite the fact that the assumed relation between the local error and h can be accepted practically only for interpolation errors, the technique breaking elements with large errors and thus leading to equidistribution of the errors is commonly used. In constructing an anisotropic h-adaptive strategy for viscous flows we encounter the following problems:

234

W Rachowicz/Comput.

Methods

Appl.

Mech. Engrg. 146 (1997) 231-252

(1) The strategy should lead to the reduction of the residual corresponding to the stabilized formulation of the Navier-Stokes equations and, at the same time, to the reduction of the influence of the regularizing terms on the solution. (2) We must be able to decide not only which elements are to be refined but also whether directional or isotropic refinement should be applied. (3) The adaptive technique should provide the means to control the accuracy of the solution in boundary layers, which is the primary goal of the analysis of the viscous problem. These goals should be realized in such a way that obtaining the final solution would require the least possible computational effort. Construction of the strategy satisfying conditions l-3 faces, however, a few pessimistic facts which indicate the limitations we must take into account. The reduction of the residual below a certain low threshold value is computationally unrealistic in the regions like shocks. There is no relation between the anisotropic reduction of an element size and the reduction of the error. Finally, the problem of the error estimation of pointwise values, like the fluxes along the solid wall boundary, has been solved thus far only for a very limited class of problems (e.g. via post-processing for some elliptic boundary value problems, see [6]), and it is open for flow problems. Still, despite these theoretical limitations, computational experience indicates that construction of adaptive meshes may result in very substantial improvements of solutions at relatively low cost even if it is not based on mathematically rigorous arguments. The adaptive strategy we propose is based on the following principles: (1) The adaptive process consists of two steps: (i) A number of adaptations of the mesh throughout the whole computational domain. (ii) A sequence of mesh adaptions limited to boundary layers leading to the recovery of accurate viscous fluxes. Returning back to step (i) after performing step (ii) is frequently impossible as (ii) introduces extremely small elements which may essentially reduce the timestep At and thus make any further evolution of the solution in the domain prohibitively expensive. (2) One possible way to perform the refinements inside the computational domain is to base the adaption on the criterion of minimization of the interpolation error for a fixed number of elements [5,1], i.e. to adjust the mesh to the solution from the point of view of the possibility to interpolate this solution. An advantage of this approach is that it selects the optimal direction of breaking of the element (as it was shown in [l] and discussed in Section 3.2). The essential drawback is that decisions about the refinements are based on the estimates of the second-order derivatives of the solution, which one should consider much less reliable than the undoubtedly true indication of the error--the residual. These arguments are the basis of the following philosophy of mesh adaption: the decision which elements are subsected is based on residual error indicators-the elements with the largest indicators are broken; the direction of subsecting follows the principles of mesh adjustment for the interpolation. (3) To estimate the accuracy of the viscous fluxes along the solid wall boundary we use the error indicator based on the interpolation principles. The adaption consists in subsecting elements adjacent to the boundary and updating the solution until the value of the indicator drops below a given level of the error. The procedure can be augmented by additional refinements reducing the error indicators and artificial dissipation in the remaining elements of the boundary layer. The details of the components of the strategy outlined above are presented next.

3.1. Anisotropic interpolation estimates in the L2- and HI-norms We consider optimization of meshes via h-anisotropic refinements in the sense of the ability to interpolate a given solution with the minimum error. Unlike in the case of the second-order elliptic equations where the most natural measure of the error is the H’-norm, it is not obvious which is the most adequate norm to measure the accuracy of solutions to the Navier-Stokes equations. We postpone the discussion of this problem and in this section we obtain interpolation results for stretched elements for the L2-norm of the error and we recall the analogous results for the HI-seminorm obtained in [l].

W Rachowicz/Comput.

Methods

Appl.

Mech. Engrg. 146 (1997) 231-252

235

Anisotropic h-adaptive meshes are obtained by a sequence of refinements of quadrilateral elements constituting the initial mesh. We assume that the initial mesh is quasiuniform, i.e. that the following conditions hold: RK ff~>ClH,

HK

2 C,

for all the initial elements K,

(8)

where C1 > 0, C, > 0 are constants, HK = dia(K), H = max K H K, RK is the radius of the largest circle contained in K. It is easy to see that the elements resulting from isotropic refinements also satisfy the condition of limited distortion (8)2 unlike the elements generated by anisotropic refinements: if the number of refinements in one of the directions exceeds the number of subsections in the other one the element becomes stretched and one has to reduce adequately constant C, in (8)2. The goal is to obtain interpolation results for elements which do not satisfy the quasiuniformity assumption (8)*. Consider first the situation when elements of the initial mesh are images of a reference element K = [0, 112through an affine mapping F (parallelogram elements). Then any element K is also an image of K through an affine mapping 4 (which is an adequately scaled F). Assume the elements are of the Lagrangian type of the order p. Denoting by uh the Lagrangian interpolant of u we have

(9) where z?= u o 4, i&, = uh o cj~,a(x,, x2)/8(&, 52) is the jacobian of 4. In the last step we used the classical interpolation result for a Lagrangian square element [7]. Denoting by k and 1 the versors of sides of element K and by hk, hl the sizes of its sides we can write (9) as IIU-UhIl;,K~c

[ h ,@+l)J,

(f&f$2

dx+h;(p+l)~

($)‘&I

(10)

A generalization of this result for curved elements is straightforward and merely notational. The problem is in a unique definition of sizes hk and hl for such elements. Having in mind, however, that we control rather levels of refinement than continuously varying sizes of elements, we can use purely formal definitions of hk, hl associating them with levels of refinements. Assume that the mapping F generating initial element K, K = F(k), is a homeomorphism satisfying the conditions: IIVF II,,k< CIH and IIVF-’

llrn,~< C2/H,

V initial elements K.

(11)

Let L be an element resulting from nl horizontal and 112vertical breakings. Then L is an image of a rectangle i c K’, L = F(i), the sizes of i are 2- ‘1 and 2-Q, respectively. The error estimate can be found as follows: 11 u -

uh

I&=

where we used (ll)l. thus:

h(ii- iih)2a,i;::;; d5 < CH2& -

Lh)2 dt

(12)

The last expression can be estimated using (10) (with K E J!,, hk = 2-“1, hl = 2-“*),

Using (11)~ and introducing

the notation hk = 2-“‘H, hl = 2-“2H we can write the estimate as

WI Rachowicz/Comput.

236

Methods Appl. Mech. Engrg. 146 (1997) 231-252

(14) where derivatives are calculated with respect to parameters & and t2 along curves F(., t2) and F(&, .), respectively. In the case of H’-seminorm the anisotropic estimate was obtained in [l]. The estimate is limited to parallelogram elements only and it is of the form:

where (Yis the angle between the sides of the element. If regularity of u admits the choice i = p + 1 then the second and the fourth expressions are of higher order and can be neglected leaving us with a formula similar to (14). 3.2. Optima&y conditions for h-adaptive meshes In this section we investigate the problem of constructing optimal meshes in the sense of minimizing the interpolation error with the number of elements being fixed. It was shown in [l] that minimization of the expression on the right-hand side of (15) (simplified by dropping out the higher-order terms) subject to the constraint of the fixed number of elements leads to the optimality criterion as the equidistribution of the element error indicators ‘&K, T)~,~: l+lcotaI

%+,K = 2

_

%,K =’

sin2 ff

2

“:J,(g)‘dx=A

l+Icot~I h;L(g)2dx=A 2

(16)

sin2 LY

throughout the computational domain. An identical argument results also in the equidistribution of the L2-error indicators:

applied to the L2-error estimate

(10)

(17)

(Minor modifications are required to cover curved elements). Conditions (16) or (17) imply the mesh adaption strategy which consists in calculating nk,$, ~L,Kfor all the elements and subsecting the elements with indicators exceeding a certain level, say p maxK (nk,J, q[,~), 0 < /I< 1, as such a process leads to the approximate equidistribution. Practical effectivity of this approach was demonstrated in [l]. The mesh adaption strategy considered in Section 3.4 takes advantage of the fact, that conditions (16) or (17) imply an optimal ratio hk/hl of sizes of elements. In the ideal situation, if hk and hl could assume any values, their ratio would be dictated by exact satisfaction of (16) or (17). This is, however, impossible: the element sizes can only be 2-“, n = 0,1, . . . fractions of the size of the initial element. Intuition suggests that the best ratio hk/hl should be selected such that nk,K/nl,K is as close as possible to 1. In fact the optimal value of ‘r&$/n[,$ must satisfy

(18)

?I! Rachowicz/Comput.

Methods Appl. Mech. Engrg. 146 (1997) 231-252

237

where CY= 2P+’ for L2-errors, cy = 2P for Hr-errors, as by a single directional breaking nk/q[ changes by a factor (Y*I . This quick conclusion is supported by a more rigorous reasoning. Assume that in some sufficiently small subdomain covered by s identical elements the integrands in (17) are nearly constant. Let 77/,> nl. Assume that t0 achieve the Optimal ratio r)k/v[ we perform m refinements reducing hk 2m times. Then isotropic refinements become optimal. As the total error reduction due to directional breaking is a = a2 we find that logarithms of the errors corresponding to anisotropic and isotropic refinements are the following functions of levels of refinements m, and mi: ya = ln(f$P~ + e,‘),

yj = In(P

. c)

(19)

where e: = ST;, ei = s$, c is a constant. We associate levels of refinements mar rni with the number of elements being produced N, = cr2”0, Ni = ~~2~~1,where cl, c2 are constants. Then the errors (19) are expressed as functions of logarithms of N,, Ni, n, = log, cr + m, and IZ~= log, c2 + 2mi: y,(n,) = ln(e~a’og~‘l+“~+ ei),

yi(ni) = lnc + (;-!!?!$!5)

m is the optimal number of directional refinements ycl(n) is less steep than yi(n) for II 3 m i.e. if:

y.(m)-y&-l)<:

if function yO(n) is steeper than yi(n) for n 6 m and

lna
Denoting by /3 = e~u’“szcl+m /ei the optimal ratio $/$

P+l

Pla+l‘

< uw

lna

6

(21) we can write (21) as

pa+l

(22)

P+l

which is equivalent to (18). 3.3. Residual error estimate In the case of Co-continuous approximation the residual corresponding to an inviscid flow is a piecewise polynomial possibly discontinuous function. If the viscous flow is considered, however, the situation is more complicated: the generalized second-order derivatives of the Co-continuous solution define the residual and the techniques devoloped for elliptic second-order equations must be applied. In this work we use the element residual method for symmetrizable elliptic problems presented in [8]. The method takes advantage of the approximate symmetrization of the bilinear form B(., .) corresponding to the weak statement (7). In the first step for every element K a small local boundary value problem is solved: find +K E V:,,(K) such that

WTakdJJnk

II

dS,

VM’ E V;,2(K)

(23)

where BK, L,y are restrictions of the bilinear and linear forms from the weak statement (7) to element is the space K, akl = At(TAkh + VdaklI +Kkl). [I u II s t an d s f or a jump of u across 8K and VE,2(K) of biquadratic element shape functions which vanish at vertices. Then, element error indicators are calculated as ei=

BK(#K,A~'(~K)

(24)

where A0 is the symmetrizer matrix for the Navier-Stokes equations defined in (6). We emphasise that the technique estimates the element contributions to the residual corresponding to the stabilized equations. Thus, the refinements must reduce not only the residual but the contributions of the regularizing terms as well.

238

H! Rachowicz/Comput.

Methods Appl. Mech. Engrg. 146 (1997) 231-252

3.4. Adaptive strategy Following the previous discussion in the first step of the adaptive solving of the flow problem we enhance the accuracy of the solution in the computational domain as a whole. The procedure consists of a number of mesh refinements followed by solving the problem on the new mesh. The decisions about subsecting of elements are based on values of residual error indicators eK while the choice of the direction of breaking results from the comparison of the directional error indicators qk,K and q[,$ and the requirement to reach the optimal value vk,K/~l,K of (18). At this point we must decide which norm, H’ or L2 is more appropriate to measure the interpolation error in flow problems. Taking into account the fact that the elliptic-like contribution in Navier-Stokes equations due to the true viscosity terms is practically negligible in most of the computational domain except the boundary layers we are more inclined to use the L2-norm as the basis for mesh optimization. An additional argument for this option we obtain by considering mesh optimization for shocks. If orientation of elements and a shock is given the optimal sizes of an element are such that the shock is parallel to the diagonal of the element. This is because it minimizes the size of the element in the direction n = (nl, n2) perpendicular to the shock with the number of elements along the shock being constant. If we consider interpolation of a model thin shock-like function f(nlxl + n2x2) then (a2f/ak2)2 = (k . t~)~f”~ and ($f/a12)2 = (I . r~)~f’~. Th us conditions (17) for optimization in the L2-norm and (16) for the H’seminorm imply the following relations for mesh parameters hk and hl: h;t(k . rQ4 = h;(/! . n)’ and hi(k . IZ)~ = hf(l . II)~

(2%

We observe that the first formula corresponds to the optimal situation with the shock being parallel to the diagonal when hklk . rz1= hllZ . nl. Therefore we assume that parameter cx = 2l+’ = 4 as for the L2-norm optimization and for p = 1 being the order of elements. Finally, we take into account that the solution is a vector-valued function. We decide to be rather conservative in introducing directional refinements and we break an element directionally only if this is indicated by the ratio qk/‘r& for all four components of the solution. The adaptive procedure includes also possible unrefinements as mesh modifications may cause significant shifts and changes in the distribution of errors. Combined together the algorithm of mesh adaption can be written as follows: (1) Converge to the steady state solution on a current mesh. (2) For all elements K calculate the residual error indicators eK and the directional interpolation indicators qi K, qfK, i = 1,. . .4, with i numbering the solution components. Find emax= maxK eK. Disrupt the brocess if the accuracy of the solution is sufficient. (3) Cluster back the elements resulting from previous refinements for which (CK ei)li2 6 alemax, where summation includes elements K resulting directly from a previous breaking of the common element-father. (4) Subsect elements K for which eK 3 (Y2emax.The way of breaking is chosen as follows: if

T&/$,K< a-l/2

if

~i,~Ir);,~ > a112

,

I

fori=l,...

fori=l,...

,4 refine directionally

,4 refine directionally

reducing h,,

(26)

reducing hk,

(27)

if none of the two conditions is satisfied-refine isotropically. (5) Go to 1. Above al, a2 are the user specified parameters, 0 < (Y~< CY~ < 1. A separate problem is estimation of the second-order derivatives of the piecewise bilinear solution. It can be performed by using appropriate finite difference quotients extended to elements-neighbors of the given element. Another option is to apply a certain post-processing returning continuous first derivatives of the solution, for instance discontinuous VU can be L2-projected onto the finite element space [9]. Unfortunately neither of these procedures is mathematically justified, nevertheless they are frequently used in practice.

W. Rachowicz/Comput.

Methods

Appl.

Mech. Engrg. 146 (1997) 231-252

239

3.5. Adaptive recovery of viscous fluxes

Extremely small values of the viscosity and conductivity constants p and K result in the fact that even if relatively accurate values of the solution of the flow problem are obtained the gradients of the solution in boundary layers can differ from the expected exact (experimental) values even by orders of magnitude. Accurate simulation of viscous fluxes (proportional to gradients) requires much finer approximation in boundary layers than is necessary to approximate U. The frequently used technique of obtaining reliable fluxes consists in generating a mesh which introduced a priori a thin layer of very fine elements along the solid wall boundaries. It requires estimation of the necessary characteristics of the layer which may be unavailable in general situations. In our approach we take advantage of the unique ability of the adaptive finite element method to adjust to the features of the solution and use it in boundary layers as well. We accept the assumption that by the global step of refinements sufficiently accurate values of the solution have been obtained. The goal is to adapt the mesh in boundary layers such that the required accuracy of gradients of the solution in this region is reached. Defining the problem in this way we want to specify the criteria which should be met to accept the solution in the boundary layer region. Or, to be more realistic, the criteria which if violated, would allow one to disqualify the solution. We assume that the errors in viscous fluxes, the skin friction and the heat flux along the solid wall boundary, are associated exclusively with the errors in gradients of the solution. The skin friction and the heat flux along the solid wall can be written as (28) and 86

K

q=Kdn=C,

_--_-

(29)

where IZis the unit vector normal to the boundary, s-the unit tangent vector, and the relations ui = mile, 8 = [e - (mi + m$/(2e)]/(CV~) and the conditions mi = 0, i = 1,2 along the wall were used. Therefore, the errors in r and q, AT and Aq, can be expressed by the errors of gradients of the conservative variables, A(&/%), A(ami/&z), A(&/&) as follows:

(30)

(31) The estimation of the accuracy of the derivatives of Q,mi, e with respect to n is based on the elementary one-dimensional interpolation estimate for linear approximation:

(32) resulting from the Taylor expansion u(h) - u(0) = h g

(0) + 5

2

(t),

for some 5 E [O,h]

where h denotes the size of the element measured along it, duh/dn = (uh(h) - uh(O))/h, and u denotes one of the components of U, uh its interpolant. Obviously, finite element solution does not coincide with the interpolant of the exact solution. Still, the large value of the indicator

rl”=2

h

a2u I 32 Im,()h] >>

(34)

is a convincing reason to claim a poor accuracy of du/dn. Therefore we accept 77”as an error indicator for au/an with d2u/dn2 being estimated by appropriate finite difference quotients of uh. The adaptive

240

W Rachowicz/Comput.

Methods Appl. Mech. Engrg. 146 (1997) 231-252

strategy based on the error indicator introduced above concerns the elements adjacent to the solid wall boundary. It consists of the following steps: (1) Converge to the steady state solution on a current mesh. (2) For all elements K adjacent to solid wall boundary estimate the errors of viscous fluxes ArK, AqK using formulas (30) and (31) with (35) Disrupt solving the problem if sufficient accuracy was achieved. (3) Refine elements K for which the relative errors of fluxes exceed the preassumed thresholds: AUK max(TK,

>a hf)

or

AqK max(qK7 qref)

>a ’

(36)

where a is the prescribed relative error (say a = O.Ol), qef, qref are the user specified reference values, rK, qK are the actual values of fluxes in K. (4) Go to 1. In general, the elements are subsected in the direction parallel to the boundary. We admit, however, isotropic refinements if condition (26) is not satisfied, with qk, 71 and (Ycorresponding to H’-seminorm as viscosity terms dominate in boundary layers. Since this procedure breaks only elements adjacent to the boundary, it generates the mesh with element sizes grading geometrically as we move away from the wall. Thus, even if errors ArK, AqK are reduced below the prescribed level there is no guarantee that the accuracy of the approximation in the remaining part of the boundary layer is sufficient. To help this problem we can extend the procedure by additional refinements leading to the reduction and equidistribution of the residual error indicators (24) and/or the interpolation error indicator for VU being the generalization of (34): (37) Obviously, such refinements are analogous to the previous modifications in the computational domain as a whole. This time, however, the errors in elements adjacent to dJ2 indicate the accuracy to be reached in the rest of the boundary layer. We mention also, that such a procedure requires some concrete definition of the boundary layer. Finally, the essential factor which may affect the reliability of fluxes is the presence of the artificial dissipation, as expressed by coefficients T, and vd in (4). Ideally, one would like to reduce them essentially below the natural viscosity coefficients. This could also be achieved by the adaptive procedure similar to the above. 4. An iterative solver A frequently raised critical argument against using significantly stretched elements is the ill-conditioning of the system of equations resulting from the discretization which destroys the convergence of the majority of iterative solvers. Thus, application of the suitable preconditioner capable of neutralizing this drawback is an important aspect of the successful application of the anisotropic h-adaptive finite element method. In [2] it was shown that the standard additive Schwarz method satisfies such requirements for symmetric elliptic boundary value problems. The analysis of the Schwarz methods as preconditioners was presented by Dryja and Widlund in [13] for the additive symmetric case, by Cai and Widlund in [14] and [ll] for the additive and multiplicative algorithms used to some nonsymmetric problems. The multiplicative methods for symmetric boundary value problems were analysed by Bramble et al. [12]. The successful application of Schwarz methods for solving symmetric elliptic problems on stretched meshes was the inspiration to use them to solve the nonsymmetric equations corresponding to discretization of flow problems. The algorithm can be summarized as follows:

W Rachowicz/Comput.

Methods Appl. Mech. Engrg. 146 (1997) 231-252

241

(1) We decompose the h-adaptive finite element mesh into N overlapping subdomains .ni. They are defined as supports of base shape functions corresponding to an auxiliary coarse mesh covering the computational domain 0. The coarse mesh is generated by performing a sequence of isotropic refinements of the elements of the initial mesh continued until the user specified relation between the sizes of the coarse mesh and the elements of the actual adaptive mesh is reached. For the details see [2]. (2) With each of the subdomains L!i we associate the local finite element space V: = Vh n Hi(Q). In addition, we define the space V: c Vh of bilinear shape functions associated with the coarse mesh. Subspaces V/ are used to define the additive or multiplicative Schwarz methods, which can be identified with the block-Jacobi method (if Q’s do not have common internal nodes) or with the Gauss-Seidel method. Formally the operators defining the two algorithms can be written as N

N

G = Z- C

RTA,‘RiA

and

G,Y= n (I - RTA,‘RiA)

I=0

(38)

i=O

where the notation is as follows: A is the global stiffness matrix, Ai, i = 1,. . . , N are the stiffness matrices associated with the subdomains Q. Matrices Ri extract from the global vector of unknowns the degrees-of-freedom associated with Qn,,while operators RT extend (by zeros) degrees of freedom of Q to the global vector. Analogously A0 is the stiffness matrix corresponding to V,h c Vh and &,R;f associate degrees-of-freedom of Vi with the global ones. Operators G, G, consist of a sequence of local projection operators Pi represented by matrices Z?TAly’Z?iA.They correspond to the following boundary value problems: find PiU E Vi” such that B(PiZl,vi) = B(u,Ui),

YUi E V/

(39)

(3) According to the idea of preconditioning by the standard iteration we use the conjugate gradientlike methods to solve the preconditioned systems M-‘(AU -f) = 0 with the preconditioner M-’ = (I - G)A-’ or M-’ = (I - G,)A-’ i.e. the systems: (I-G)U=k

or

(I-G,)U=k,

(40)

where f is the original load vector, k, k, are the modified right-hand side vectors and Z denotes the identity matrix. Since the formulation of flow problems is nonsymmetric we choose the Gauss-Seidel version of the preconditioner. It updates the solution more rapidly so it seems to be more efficient if no advantage of the symmetry is to be taken. The global step of the iterations improves the performance of the solver as it provides the exchange of the information in the whole 0. It requires, however, significant additional memory space. Thus application of this step may be considered a trade off between the speed of computations and the available size of the problem. Before we investigate the characteristics of G, corresponding to the bilinear form of (7): B(U,V)

=

s R

VTUdx+At

sR

VT (K;, + TaAiAj + vda;j) U, dx

(41)

let us point out the following features of B: (1) The quadratic form Uij = Kij + TaAiAj + v&iZ is nonsymmetric: $ # a;i. (2) Form aij is not necessarily uniformly positive definite. This is because matrix K of combined blocks Kij (i.e. VUTKVV = U,iKijV,j, YU,i, V,j) is of rank 5 in R8. Also, vectors Cf=, AiU,i may vanish for some U,i # 0. Moreover, coefficients T, and vd of stabilizing terms can be reduced significantly below the range of magnitude of coefficients Kij. These properties together with very rapid variation of r, and vd (depending upon U” and the finite element mesh) make the analysis of the preconditioner a difficult task. The lack of symmetry and uniform positive definiteness of aij does not allow us to use the results of Cai and Widlund [ll] which consider only the case with nonsymmetry caused by the presence of the first order terms. To our knowledge, the situation with nonsymmetric aij have not been investigated thus far. Within the existing theory we are

242

II! RachowiczlComput.

Methods Appl. Mech. Engrg. 146 (1997) 231-252

able to analyse the performance of the domain decomposition preconditioner for the alternative to (4) symmetric formulation of the SUPG method in the entropy variables. In this case the difficulty which remains is the lack of uniform positive definiteness of the second-order contributions. The discrete equations for the entropy variables formulation are obtained analogously to (7) and the corresponding bilinear form is similar to (41): B(U, V) =

I

D

VTAoUdx+At

IR

VT(Kij + ATTAj + vdSijAo)Uj dx

(42)

In (42) viscous matrices Kij and jacobians Ai correspond to viscous and Eulerian fluxes in the entropy variables formulation. The symmetric positive definite operator r and parameter vd > 0 are the functions of the finite element mesh and the current solution, for the details see [4]. The analysis of the multiplicative preconditioner for (42) can be based on the results of Bramble et al. [12]. In this work a general situation of solving the equation B(u, v) = (f, u),

vu E v

(43)

in a finite or infinite dimensional Hilbert space V with the inner product B(., .) is considered. vergence of the multiplicative algorithm based on the decomposition V = VI + . . . + V, (Vi of V), and the projections Pi (39) is shown providing that the projections satisfy the following there is Cc > 0 such that for every u E V there is a decomposition u = CL, ui, ui E Vi such

The consubspaces condition: that (44)

i=l

In the case (43) is a finite element formulation sR

of the elliptic boundary value problem:

ra,u,iu,jdr+Snbuudx=S,fUdr,

(45)

with (aij)zxz being a symmetric uniformly positive definite matrix and b > 0, condition (44) was proved in [13] with Co being independent of the discretization step h providing the mesh is quasiuniform (with V = Vh, Vi = V,h, i = 1,. . . ,N). The required decomposition in (44) can be defined using the coarse grid of subdomains Q. We consider a partition of unity on f2, 1 = CL, f3i(x), where @i(x) are the base shape functions corresponding to the coarse mesh. Clearly, 1 3 13,3 0, ](00, ]]_,a< CZZ-‘, i = 1,. . . , N, where Hi is the size of Q. The partition of unity is used to define the decomposition of any U E Vh: U=etYi,

Ui=Zh(@Uj,

i=l,...

,N

(46)

i=l

where Zh is the standard nodal interpolation operator. Property (44) is the basis for establishing the estimate of the largest eigenvalue of G, [12]: (47) where H = minZ!Zi, n = maxni with Eli being the number of subdomains Qj intersecting Q. In [2] the assumption of quasiuniformity of the mesh used in showing (44) was relaxed admitting the use of significantly stretched elements. Proving (44) for the flow problem formulation (42) requires relaxing the assumption of uniform positive definiteness of Uij. Our analysis of the preconditioner relies on the following properties of the matrices defining B(., .) in (42) which hold for any physically acceptable Bow field U”: (1) The symmetrizer A0 = A;f is uniformly bounded from above and it is uniformly positive definite. (2) Matrix r = 7T 3 0 and coefficient vd are uniformly bounded. (3) Matrix K defined by VUTKVV = UTKijV,j is symmetric semipositive definite, dim(ker(K)) = 3 in R8. For vectors W E (ker(K))l the quadratic form WTKW is uniformly positive definite, i.e. the smallest nonzero eigenvalue of K is uniformly bounded from zero. (In the formulation using

WI Rachowicz/Comput.

143

Methods Appl. Mech. Engrg. 146 (1997) 231-252

the natural variables, density, velocities and the temperature, K represents the linear elasticity and the Laplace operators with uniformly bounded coefficients.) We consider inequality (44) for each of bilinear forms constituting B(., .) separately. Inequality

CJ 4

i

Zh(0iU)TAo~/$3iU)

dx < C

Jn

UTAoUdx,

VU E Vh

results from elementary

properties of Zh. Consider the contribution of quadratic form Kij. Using the notation as above we need to show that

CJ

VUTKVU dx +

VZh(6iU)TKVZh(BiU) dX 6 C

Q

i

sn

with C being independent of the finite element mesh. Inequality of this type was shown in [2] in the case the first integral on the right-hand side is a norm in H,‘(n). The major difficulty now is that, as it was mentioned, the kernel of matrix K is nonempty. Thus it is not obvious whether the expression involving derivatives VU E ker(K) on the left-hand side of (49) can be estimated by an exclusively L2-norm of U. Moreover, the required estimate must be independent of, possibly unbounded, aspect ratio of elements. We begin with splitting each of the integrals over Q into contributions from elements K. In each K we express the integrand as Zh(&U) = &U + zh(pU), where 8i is the average value of 6, in K, /_?= 0, - 8i, which implies that

s

VZh(C3iU)TKVZh(0iU) dX < 2Gi2

K

sK

VUTKVU dx + 2

sK

VZh(pU)TKVZh(PU)

do

(50)

Our goal is to bound the second term on the right-hand side of (50) by the sum of JK VUTKVU dw and ]]u ]]$(K,* We simplify the discussion by considering only rectangular elements of sizes hi and h2 with the unrestricted value of hl/h2. Since Zh(pU) is a bilinear function its differentiation can be replaced by taking finite difference quatients. Denoting by U ,,, Ui2, U2i, U22 the nodal values U(O,O), U(0, h2), U(hl,O), U(hl, h2) in local coordinates such that K = [0, hl] x [0, h2] we find that VZh(pU) = Wi + W2 where U2l

-

Ull

h

+ Pll

W1 =

u12

-

u11

2

+ Ull

P21 -

2 W2 =

52) + U’“i

Ul2

P22 + Pl2

u12+u11

2

52

2

P12 + Pll

h2 u21

(1 _

2

hl

Pll

u22

(l-51)+

(1 _

-

u21

P22 + P21

2

h2 52) + u22;

u12

P22 -

hi

PI2

51

(51)

1

(2

ht

P12 h2 -Pll

(1 _

&)

+ u22

;

u21

(52)

p22 h, p21

g-1 I

& = xi/hi, the subsequent

rows of Wi and W2 correspond (1 _

h)

+ u22;

u12

52

(1 -

51) + (122;

u21

51

1

to aZh(pU)/axi,

i = 1,2. Similarily

(53)

(Let us point out that, in general, Wi is not co-linear with VU as the scaling coefficients pZl + fill and &2 + pi2 are the opposite numbers: their sum vanishes as Pii are the nodal values of the bilinear function of zero average value.) We observe that

J

VZ,(pU)TKVZ,(PU)

dx < 2

K

6 c

sK

W;KW2 dx + 2

11vP

ii;“(K)ii

u

sK Ilt?(K)

WTKWl +2

J

dx

K WTKWl

dX

(54)

244

W Rachowicz/Comput.

Consider the contribution

Methods Appl. Mech. Engrg. 146 (1997) 231-252

of IV, in (54). If K were uniformly positive definite then we could write

completing the proof. In the case K 3 0 the last step fails for U such that VU E ker(K). REMARK

meshes we can express WI as

1. In the case of quasiuniform

w

=

(-(U21

-

UllW

-

65) + (U22

-

Ul2)52)

$

$

(-(Ui2

-

~11)(1

-

51) + (U22

-

U21)51)

2

g

I

1 I

(56)

where derivatives of p are evaluated at the center of K. This, in view of the uniform bounds ct < ht /hZ < c2 for some 0 < cl < c2, implies that

J

w,Kwl

b <

c 11Ofi Iit-(K,II u Il:r(K,

(57)

K

which suffices to show (44). To proceed in a general case we must investigate in detail the properties of matrix K. Following [lo] Appendix A we find that ker(K) is spanned by the following vectors in RR: VI = (l,O,O, o,o,o,o,o,

U? = (O,O, o,o, l,O, O,O)T,

o)T,

-l,O, o)T

(58) where the first 4 components correspond to differentiation with respect to x1 the remaining to x2. Assume first that VU = clu;, i = 1,2. Then obviously WI = c2ui and (55) is trivially satisfied. Consider the case VU = cu3. Then, by comparing the third and the sixth rows of (53) (both vector-valued rows of (53) consist of 4 components) we find that _

u3

=

(O,O,

l,O,O,

4

(59)

for all &, t2 E [O,11, where superscripts denote the vector components 3 l421 -

UT,

=

u;2

-

2 92

14:2,

-u,,

2

-

U:l

and

=z.4;2-u;l

z4;1

of U. (59) is satisfied if _ _-

2 q2

-

Ufl

(60)

hl h2 The first two conditions mean that U is a linear (not a bilinear) function of x,, i.e. that VU is constant in K. Taking advantage of (60) we can write Wl as follows: 0 0

41hl

U:l

821 + Pli

(I

_

t21

+ P22 ;

PI2

52

1

0 n

0 3 _ %

P12fPll

(,_[,)+P22+P21

2 0

I

-(uf, - uf,) -;;

0

0

0

0

-1

El,

2 0

0

0

2

w, =

0

0

0 0

(Gl

-

u:,>

g

(h- i)

(52

-

;>

0 0

(61) where a = (~1, - u:t)(/&t + PI2 + p2t + p22)/(4hl) and where we used the fact that /3 is a bilinear function and derivatives of /3 are evaluated at the center of K. Since the first component of WI is in ker(K) the decomposition (61) implies that

W Rachowicz/Comput.

Methods Appl. Mech. Engrg. 146 (1997) 231-252

Finally, we collect (50) and (54), inequality (55) which holds if VU E ker(K)I, VU E ker(K) obtaining that

34.5

and the above results for

(63) where we used the bound ]I00, ]lm,o,6 C/Hi.

REMARK

1,. . , N the of the sizes of subdomains Hi (see [2]).

2.If the global, coarse grid subspace Vi is added to the local subspaces VI!, i =

decomposition

(44) can be found such that C0 is independent

Next, we analyse the contribution of the quadratic form VUTAVU := UTATTAiU,j to inequality (44). Form A is positive semidefinite. The analysis can be performed in a similar way as for K. This time, however, we cannot use the argument of decomposing VU into contributions from ker(A) and ker(A)l as the nonzero eigenvalues of A may be arbitrarily close to zero. We avoid this problem by using form K instead of A in the bounding expression. Consider vector Wi defined as in (51). Let A > 0, A > 0 be defined as A=

Then, for any U such that VU E ker(K)l

s

inf

1,

VUTKVU

dX

we can write

WTAW1 CI-Y< CA IIp llt=(K)ll VU llf,:~K~G C t

VUTKVU

II P II~LX(K,J K

K

(64)

VUc(ker(K))L

dx

(65)

which corresponds to (55). In boundary layers, where significantly stretched elements are introduced, the magnitude of the stabilizing terms ATrAj is comparable to that of Kija Therefore, A/A is not expected to be large. Consider VU E ker(K). For VU = ui , i = 1,2 given in (58) a direct calculation shows that estimate (55) with A replacing K is satisfied. Vector u3 of ker(K) is also in ker(A) therefore an estimate analogous to (62) holds. Finally, the isotropic contribution Yd&jZ in B(-, -) can be easily included in inequality (44).

5. A numerical

example

As an illustration of the methods developed in this work we present the analysis of a known benchmark problem: a hypersonic flow over a 20” compression corner. The problem was investigated experimentally by Holden [16] and numerically in a number of works (e.g. [17-191). The geometry of the problem and the h-adaptive finite element mesh are shown in Fig. 1. The statement of the problem is as follows: Mach number Reynolds number Inflow temperature Wall temperature

M = 14.1 Re = 72000 ft-* T, = 88.9”K Twali= 294.4 “K

The problem was solved using the anisotropic h-adaptive strategy. The initial mesh consisted of 23 x 11 elements. Four levels of adaptive refinements were introduced to obtain the solution with the acceptable accuracy inside the computational domain. The adaptive strategy of Section 3.4 based on the residual error estimate and directional refinements was applied. In the second step of the adaptive process the strategy of improving the solution in the boundary layers was applied (Section 3.5). Maximum of 11 additional levels of directional refinements were introduced making the total of 15 levels of subsecting of elements along the solid wall. Thus the minimum size of elements was 2-15h ‘v 1/32OOOh, where h = 0.01 is the vertical size of the initial mesh. In the process of solving the problem we encountered a few numerical difficulties:

W Rachowicz/Comput.

246

Methods Appl. Mech. Engrg. 146 (1997) 231-252

M = 14.1 Re = 72000ft-1

J

Fig. 1. Geometry

of the problem

and an h-adaptive

(1) The

mesh.

solution was loosing stability at the stagnation point. We helped this by adding artificial viscosity due to Lohner et al. [20] in elements directly surrounding the tip of the plate. (2) The size of the elements about the tip was significantly affecting the location of shocks. We stabilized this location by additional reduction of sizes of the elements in the vicinity of the tip. The reflected shock is especially strong in the shock-boundary layer interaction region with the (3) pressure ratio exceeding 100. We had to increase coefficient Vdof the shock capturing operator by 1.5 to prevent overshoots in the pressure resulting in nonphysical negative pressures. The increase of vd unfortunately deteriorated the solution in the rest of the computational domain. We attribute the insufficient value of Vd given by (5) to the fact that for significantly stretched elements the ratio of Vd to the element size in the direction perpendicular to the shock is ~6 times less than for undistorted elements. Fig. 2 presents some characteristics of the solution in the boundary layer region and the finite element ._ mesh displayed with the logarithmic resealing in the vertical direction, y’ = log,(y + E), with E = 2-lZh. Such resealing maps the original geometrically graded mesh into the mesh close to uniform in the vertical direction, making the presentation possible. The region being shown is l/8 of the initial element size wide and it covers a half of the inclined wall. Figs. 2 (a-h) show distributions of the following characteristics: (a) Local Mach number (to indicate the range of the boundary layer). (b) Residual error indicators (24). (c) Interpolation error indicator for VU, (37). (d) SUPG parameter T,, (5). (e,f) Parameter Vd (5) and the natural viscosity coefficient p both scaled in a fashion suggested in [4] to compare them with 7,:

PI@ v;l =

(c + ,:.

412

(notation as in (5)).

and

ii =

(c + 124 .

PI)2

W Rachowicz/Comput.

b

a

10**

247

Methods Appl. Mech. Engrg. 146 (1997) 231-252

-6

.OlG25

0.0

0.0

0.0 c

e Fig. 2. Characteristics

0.o

of the solution

(g,h) Local contributions natural dissipation:

in the shock-boundary

.!T

layer interaction

to the global entropy production

dUT al = G~A~A;‘A~~ I

I

and a2 = !!?A-‘K,i axi o

11

region

(details

in the text).

rates associated with the artificial and the

g

(67) I

In the adaption process we might perform additional refinements in boundary layer regions in order to reduce the SUPG stabilizing coefficients 7, and Vd. The second one turned out to be negligible if compared to natural viscosity p but it was not the case with r, which even exceeds ji. On the other hand the effects of the term r, and the natural viscosity on the solution as measured by the bilinear forms al and a2 (67) suggest a sufficient reduction of the influence of stabilizing mechanism on the solution. Finally, we observe that the residual error indicators and interpolation indicators are not sufficiently reduced inside the boundary layer if compared to the values obtained along the wall. In the adaption process, reduction of these indicators might be taken into account too. It is in a sense unexpected that despite neglecting to do this and focusing exclusively on the reduction of the error in fluxes along the wall, their values are satisfactorily accurate. Figs. 3-5 present the skin friction coefficient, the heat flux coefficient and the pressure coefficient compared with the experimental data of [16]. The viscous fluxes match the experiment even though the pressure is underestimated indicating the possibility of further adaptive improvements of the solution inside the computational domain. Fig. 6 shows the error indicators of the viscous fluxes, max(ArK/qef, AqK/qref) with qer = maxK TV, qref = maxK qK, obtained after performing 15 levels of refinements. Figs. 7 and 8 present the distribution of density Q and the vertical component of the velocity u. The total number of nodes in the final mesh was 6200, the number of elements reached 6700.

W Rachowiez/Comput.

248

Methods Appl. Mech. Engrg. 146 (1997) 231-252

Fig. 3. The skin friction

coefficient,

C,- = 27.

1.925 !

-0.525

1 0.

,

,

I 0.75

”

’

I

’

“I

1.75

Fig. 4. The heat flux coefficient,

/ 2.75

Ch = 29.

1

I 3.5

WI Rachowicz/Comput.

249

Methods Appl. Mech. Engrg. 146 (1997) 231-252

0.75

0.3

I

-0.15

I

0.

(

I

I

Fig. 5. The pressure

0.0525

(

I

I

I

1.75

0.75

coefficient,

(

I

(

2.15

3.5

2.15

3.5

C, = p/(l/2~,u,)

_

0.04125 _

0.02625 _

0.01125

0.

Ik 0.

0.75

Fig. 6. Error

1.75

estimate

for viscous

fluxes.

250

WI Rachowicz/Comput.

0.

Methods Appl. Mech. Engrg. 146 (1997) 231-252

I

I

a.15

1

17.5

Fig. 7. Distribution

0.01

28

of density

36.15

Q.

I

/

I

0.1225

0.235

0.3475

Fig. 8. Vertical

component

of velocity

u.

0.46

W Rachowicz/Comput.

6. Concluding

Methods Appl. Mech. Engrg. 146 (1997) 231-252

251

remarks

An anisotropic h-adaptive finite element method is ideally suited for solving compressible NavierStokes equations. First of all, it provides the perfect approximation in boundary layers, where solutions are close to one-dimensional. Moreover, the size of stretched elements in such regions is found automatically in an adaptive process in which the presumed accuracy of viscous fluxes is obtained. The capability of the anisotropic adaptive meshes to adjust to irregularities of the solution in the rest of the computational domain (especially to shocks) is also significantly better than for isotropic method. The anisotropic approach saves degrees-of-freedom especially in regions where directional features of the solution align with the mesh lines. The adaptive strategies we propose in the work are based on interpolation error indicators. They are not satisfactory in a sense that the low level of the estimated interpolation error does not guarantee a good accuracy of the solution (a drastic example: for a constant solution interpolation error indicators vanish). In the first, global step of mesh adaption we overcome this drawback by using the residual error indicators to decide which elements to subsect while interpolation indicators decide only about the direction of breaking. The adaption in boundary layers, however, must be considered in the following sense: we refine the elements to eliminate the detected inability to approximate large curvature of the solution as expressed by the indicators (32). Elimination of the drawbacks mentioned above can be anticipated in developing the residual error indicators which would take into account directional features of the solution and local error indicators for fluxes. Finally, we provide the iterative solver which makes application of the anisotropic h-adaptive meshes possible. The convergence of the solver is insensitive to the aspect ratio of the elements and it allows one to adjust the sizes of local problems being solved on subdomains to the rapidly changing density of elements. Application of the proposed techniques to solve a model benchmark problem, the high Reynolds number hypersonic flow with the shock-boundary layer interaction, confirmed their practical usefulness.

Acknowledgment

The support of this work by Komitet Badan Naukowych under Grant no 7T07A 010 08 is gratefully acknowledged.

(KBN) of the Polish Academy of Science

References [l] W. Rachowicz, An anisotropic h-type mesh-refinement strategy, Comput. Methods Appl. Mech. Engrg. 109 (1993) 169-181. [2] W. Rachowicz, An overlapping domain decomposition preconditioner for an anisotropic h-adaptive finite element method, Comput. Methods Appl. Mech. Engrg. 127 (1995) 269-292. [3] G.J. Le Beau. SE. Ray, S.K. Aliabadi and T.E. Tezduyar, SUPG finite element computation of compressible flows with the entropy and conservation variables formulations, Comput. Methods Appl. Mech. Engrg. 104 (1993) 397422. [4] SK. Aliabadi, S.E. Ray and T.E. Tezduyar, SUPG finite element computation of viscous compressible flows based on the conservation and entropy variables formulations, Comput. Mech. 11 (1993) 30&312. [S] L. Demkowicz, P Devloo and J.T. Oden, On an h-type mesh-refinement strategy based on minimization of interpolation errors, Comput. Methods Appl. Mech. Engrg. 53 (1985) 67-89. [6] I. BabuSka, Miller, The post-processing approach and other higher derivatives of the displacements,

in the finite element method-Part Int. J. Numer. Methods Engrg.

1: Calculation of displacements, 20 (1984) 1085-1109.

stresses

[7] PG. Ciarlet, The Finite Element Method for Elliptic Problems (North Holland, 1978). [S] J.T. Oden, L. Demkowicz, W. Rachowicz and T. Westermann, A posteriori error analysis in finite elements: the element residual method for symmetrizable problems with application to compressible Euler and Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. Special Issue: Reliability in Computational Mechanics, J.T. Oden, ed. 82 (1990) 183-204. [9] J.T. Oden and H.J. Brauchli, On the calculation of consistent stress distributions in finite element applications, Int. J. Numer. Methods Engrg. 3 (1971) 317-325. [lo] T.J.R. Hughes, L.P. Franca and M. Malet, A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the compressible Euler and Navier-Stokes equations and the second law of thermodynamics, Comput. Methods Appl. Mech. Engrg. 34 (1986) 223-234.

252

W Rachowicz/Comput.

Methods

Appl.

Mech. Engrg. 146 (1997) 231-252

[Ill X.-C. Cai and O.B. Widlund, Multiplicative Schwarz algorithms for some nonsymmetric and indefinite problems, SIAM J. Numer. Anal. 30 (1993) 936-952. [12] J.H. Bramble, J.E. Pasciak, J. Wang and J. Xu, Convergence estimates for product iterative methods with applications to domain decomposition, Math. Comput. 57 (195) (1991) 1-21. [13] M. Dryja and 0. Widlund, An additive variant of the Schwarz alternating method for the case of many subregions, Technical Report 339, Courant Institute of Mathematical Sciences, 1987. [14] X.-C. Cai and O.B. Widlund, Domain decomposition algorithms for indefinite elliptic problems, SIAM J. Sci. Stat. Comput. 13 (1992) 243-258. [15] W.W. Tworzydlo, J.T. Oden and E.A. Thornton, Adaptive implicit/explicit finite element method for compressible viscous flows, Comput. Methods Appl. Mech. Engrg. 95 (1992) 397-440. [16] M.S. Holden and J.R. Noselle, Theoretical and experimental studies of the shock wave-boundary layer interaction on compression surfaces in hypersonic flaw, CALSPAN Report ARL 70-0002, University of Buffalo Research Center, January 1970. [17] R.R. Thareja, R.K. Prabhu, K. Morgan, J. Peraire, J. Peiro and S. Soltani, Applications of an adaptive unstructured solution algorithm to the analysis of high speed flows, AIAA paper 90-0395, 1990. [18] E.A. Thornton, P Dechaumpai and G.R. Vemaganti, A finite element approach for prediction of aerothermal loads, AIAA paper 86-1050, 1986. [19] D.H. Rudy, J.L. Thomas, A. Kumar, PA. Gnoffo and S.R. Chakravarthy, A validation study of four Navier-Stokes codes for high speed flows, AIAA paper 89-1838, 1989. [20] R. Lohner, K. Morgan and T. Peraire, A simple extension to multidimensional problems of the artificial viscosity due to Lapidus, Comput. Appl. Numer. Methods 2 (1985) 141-147.