The energetic implications of using deforming reference descriptions to simulate the motion of incompressible, Newtonian fluids

Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

www.elsevier.com/locate/cma

The energetic implications of using deforming reference descriptions to simulate the motion of incompressible, Newtonian ¯uids S.J. Childs Department of Pure and Applied Mathematics, Rhodes University, Grahamstown, 6140, South Africa Received 05 March 1998

Abstract In this work the issue of whether key energetic properties (nonlinear, exponential-type dissipation in the absence of forcing and long-term stability under conditions of time dependent loading) are automatically inherited by deforming reference descriptions is resolved. These properties are intrinsic to real ¯ows and the conventional Navier±Stokes equations. A completely general reference description of an incompressible, Newtonian ¯uid, which reconciles the di€erences between opposing schools of thought in the literature is derived for the purposes of this investigation. The work subsequently focusses on establishing a class of time discretisations which inherit these self-same energetic properties, irrespective of the time increment employed. The ®ndings of this analysis have profound consequences for the use of certain classes of ®nite di€erence schemes in the context of deforming references. It is signi®cant that many algorithms presently in use do not automatically inherit the fundamental qualitative features of the dynamics. An `updated' approach as a means of avoiding ever burgeoning deformation gradients and a still further simpli®ed implementation are further topics explored. Ó 1999 Elsevier Science S.A. All rights reserved. Keywords: Energy conservation; Incompressible, Newtonian ¯uid; Completely general reference description; Arbitrary Lagrangian Eulerian; ALE; Rigid body in a ¯uid; Free surface; Finite elements; New Poincare inequality

1. Introduction Descriptions of ¯uid motion are conventionally based on the principles of conservation of mass and linear momentum. One might hope that all such descriptions would accordingly exhibit key energetic properties (nonlinear, exponential-type dissipation in the absence of forcing and long-term stability under conditions of time dependent loading) consistant with the principle of energy conservation. These properties are intrinsic to real ¯ows and the conventional, Eulerian Navier±Stokes equations. A completely general reference description of an incompressible, Newtonian ¯uid, which reconciles the di€erences between the so-called arbitrary Lagrangian Eulerian (ALE) formulation of Hughes et al. [4] (deformation gradients absent) and that of Soulaimani et al. [9] (deformation gradients present, but use is problematic), is derived for the purposes of this investigation. The implications of the resulting description are investigated in the context of energy conservation in a similar, but broader, approach to that taken by others (eg. Simo and Armero [8]) for the conventional, Eulerian Navier±Stokes equations. The work subsequently focusses on establishing a class of time discretisations which inherit these selfsame energetic properties irrespective of the time increment employed. The ®ndings of this analysis have

E-mail address: [email protected] (S.J. Childs). 0045-7825/99/$ - see front matter Ó 1999 Elsevier Science S.A. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 9 9 ) 0 0 0 7 6 - 6

220

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

profound consequences for the use of certain classes of di€erence schemes in the context of deforming references. It is signi®cant that many algorithms presently in use do not automatically inherit the fundamental qualitative features of the dynamics. An `updated' approach as a means of avoiding ever burgeoning deformation gradients which arise from the accumulated step-wise deformation of meshes and a still further simpli®ed implementation are further topics explored. The main conclusions of this work are based on a new inequality and a number of lemmas. These lemmas are mainly concerned with the new convective term. The new inequality is used in place of where the Poincare±Friedrichs inequality might otherwise have limited the analysis. This analysis is extended in that non-zero boundaries, so-called free boundaries and time-dependent loads are considered.

2. A completely general reference The implementation of most numerical time integration schemes would be problematic were a conventional Eulerian 2 description of ¯uid motion to be used in instances involving deforming domains. The reason is that most numerical time integration schemes require successive function evaluation at ®xed spatial locations. On the other hand meshes rapidly snarl when purely Lagrangian 3 descriptions are used. Eulerian and Lagrangian references are just two, speci®c examples of an unlimited number of con®gurations over which to de®ne ®elds used to describe the dynamics of deforming continua. They are both special cases of a more general reference description, a description in which the referential con®guration is deformed at will. A deforming ®nite element mesh would be a good example of just such a deforming reference in practice. The transformation to the completely general reference involves coordinates where used as spatial variables only and the resultant description is therefore inertial in the same way as Lagrangian descriptions are. 2.1. Notation Consider a material body which occupies a domain X at time t. The material domain, X0 , is that corresponding to time t ˆ t0 (the reference time, t0 , is conventionally, but not always, zero). A third con®g~ which is chosen arbitrarily is also de®ned for the purposes of this work. The three domains are uration, X, related in the sense that points in one domain may be obtained as one-to-one invertible maps from points in another. x; t†; t†, can be de®ned in terms of the doFor any general function f …x; t†, a function, f~…~ x; t†  f …k …~ mains and one-to-one, invertible mappings illustrated in Fig. 1. Similarly, f0 …x0 ; t†  f …k…x0 ; t†; t† can be de®ned. This notation can be generalised for the component-wise de®nition of higher order tensors. The key to understanding much of this work lies possibly in adopting a component-wise de®ned notation. f is not e and div In contrast to the function notation just established, the de®nition of the operators r based on r and div. They are instead the referential counterparts, that is e ˆ o r o~ x

and

f ˆ o ‡ o ‡ o : div o~x1 o~x2 o~x3

The notation A : B is used to denote the matrix inner product Aij Bij throughout this work, h
;
iL2 …
† denotes the L2 inner product and k
kL2 …
† the L2 norm. 2.2. Some general results for functions de®ned on the three domains Three important results are necessary for the derivation of the completely general reference description and these are presented below. 2 3

Eulerian or spatial descriptions are in terms of ®elds de®ned over the current con®guration. Lagrangian or material descriptions are made in terms of ®elds de®ned over a reference (material) con®guration.

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

221

Fig. 1. Schematic diagram of domains and mappings used in a completely general reference description.

2.2.1. The material derivative in terms of a completely general reference ~ in terms of a completely general, reference is The material derivative of any vector ®eld v i @~ v e h ~ ÿ1 ~ ref † : vÿv ‡ r~ v F …~ ot

…1†

~ is the deformation gradient given by ~ref is the velocity of the reference deformation, and F where v 

~ x† ˆ ok : F…~ o~ x This result (taken from Hughes et al. [4]) is obtained by recalling that the material derivative (total derivative) is the derivative with respect to time in the material con®guration. Thus D~ vi o ˆ f~ vi …~k…x0 ; t†; t†g Dt ot o~ vi o~ vi @ k~j ‡ : ˆ ot o~xj ot

…2†

A more practical expression is needed for okej =ot (the velocity as perceived in the distorting reference). This can be obtained by considering kk …x0 ; t† ˆ kk …~k…x0 ; t†; t†; (see Fig. 1) so that okk okk okk ok~j ˆ ‡ ot x0 fixed ot x~ fixed o~xj ot

222

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

or ok~j o~xj ˆ ot oxk

okk ot x0

@k ÿ k ot x~ fixed

! : fixed

Substituting this expression into Eq. (2), the desired, suitably practicable result is obtained. 2.2.2. An element of area in terms of a distorting reference The second important result can be recalled from general continuum mechanics. Consider an element of area, size dA, with an outward unitm normal n. Then ~ J~dA~ ~ ÿt N n dA ˆ F

…3†

~ denote the respective analogous size and outward unit normal of this element of area in the where dA~ and N referential con®guration (N so as to remain consistent with the notation, since N~i 6ˆ ni in this case) and ~ This result is demonstrated in most popular textbooks on continuum mechanics (eg. Lai et al. J~ ˆ det F. [6]). _ 0 ˆ J0 div v 2.2.3. The kinematic result J The material derivative of the Jacobian J0 is given by the relation _ 0 ˆ J0 div v; J where J0 is de®ned as follows:   ok J0  det : ox0 This result is demonstrated in most popular textbooks on continuum mechanics (eg. Lai et al. [6]). 2.3. Derivation of the completely general equation One way in which to derive a completely general reference description of an incompressible, Newtonian ¯uid is to start with the balance laws in global (integral) form, and to make the necessary substitutions in these integrals. The desired numerical implementation (similar to the one which has been thoroughly investigated and found to be stable) is then obtained. 2.3.1. Conservation of mass Let X…t† be an arbitrary sub-volume of material. The principle of conservation of mass states that Z d q dX ˆ 0 …rate of change of mass with time ˆ 0† dt X…t† Z d q J dX0 ˆ 0 …reformulating in terms of the material configuration; X0 † dt X0 0 0 Z o fq0 J0 g dX0 ˆ 0 …since limits are not time dependent in the material configuration† X ot Z  0  _ 0 ‡ q_ 0 J0 dX0 ˆ 0 …by the chain rule† q0 J …4† X0 Z   _ 0 ˆ J0 div v† q_ ‡ q div v dX ˆ 0 …using the kinematic result J X…t† ! Z o~ vi o~xj ~ ~ q_ ‡ q J dX ˆ 0 …reformulating in terms of the distorting referential o~xj oxi ~ X…t† ~ configuration; X…t††

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

  ev : F ~ ÿt J~ ˆ 0 ) q_ ‡ qr~

223

…integrand must be zero since the volume was arbitrary†:

Thus, for a material of constant, non-zero density, ev : F ~ ÿt ˆ 0 r~

J~ 6ˆ 0

since

…mappings are one-to-one and invertible†:

Notice also that Eq. (4) implies o fq J g ˆ 0; ot 0 0

…5†

since the volume was arbitrary and the integrand must therefore be zero. 2.3.2. Conservation of linear momentum (and mass) The principle of conservation of linear momentum for an arbitrary volume of material X…t† with boundary C…t† states that d dt

Z

Z X…t†

qv dX ˆ

Z

X…t†

qb dX ‡

C…t†

rn dA

…6†

where q is density, b is the body force per unit mass, r is the stress, n the outward unit normal to the boundary and v is the velocity. The term on the lefthand side can be rewritten as follows: d dt

Z X…t†

qv dX ˆ

d dt

Z

Z ˆ

X0

X0

o fq v0 J0 g dX0 ot 0

X0

Z ˆ

X…t†

Z ˆ

…reformulating in terms of the material configuration; X0 †

Z  ˆ

q0 v0 J0 dX0

~ X…t†

…since limits are not time dependent in the

material configuration†  ov0 o q0 J0 ‡ v0 fq0 J0 g dX0 ot ot

qv_ dX

…the second term above is zero as aconsequence of Eq:…5††

~ ~_ J~dX qv

…reformulating in terms of the distorting referential

~ configuration; X† ! Z i o~ v e h ~ ÿ1 ref ~ ~ † J~dX ˆ q vÿv ‡ r~ v F …~ ot ~ X…t†

…using result …1††

where v_ denotes the material derivative of v. The surface integral becomes Z

Z C…t†

rn dA ˆ

~ C…t†

Z ˆ

~ X…t†

~ J~dA~ …reformulating in terms of a distorting reference using result …3†† ~ ÿt N ~F r f f~ ~ ~ ÿt J~g dX div rF

…by the divergence theorem†:

224

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

Finally, the term involving body force becomes Z Z ~ …reformulating in terms of a distorting reference†: qb dX ˆ q~ bJ~dX ~ X…t†

X…t†

Substituting these expressions into Eq. (6), remembering that the volume used in the argument was arbitrary and that the entire integrand must therefore be zero, the conservation principles of linear momentum and mass may be written in primitive form as ! o~ v e ~ ÿ1 fP ~ ~ ref † J~ ˆ q~ vÿv bJ~ ‡ div …7† q ‡ r~ vF …~ ot and ev : F ~ ÿt ˆ 0; r~

…8†

~ is the Piola±Kircho€ stress tensor of the ®rst kind, P ~ ˆr ~ ÿt J~. In terms of the constitutive relation, ~F where P r ˆ ÿpI ‡ 2lD, for a Newtonian ¯uid,     t   t  e vF ~ ÿ1 ~ ÿt J~ since D ~ ˆ ÿ pI ‡ l r~ e vF ~ ÿ1 ‡ r~ ~ ˆ 1 rv f ‡ rv f : F P 2 The derivation of a variational formulation is along similar lines as that for the Navier±Stokes equations (the purely Eulerian description). For a ¯uid of constant density, the variational formulation Z Z h i o~ v~ ~ ev F ~ ÿ1 …~ ~ ~
J dX ‡ q w ~
r~ ~ ref † J~dX vÿv q w ot ~ ~ X X Z Z Z Z ~ ‡ p~r~ ~N ~ dC ~ w† : D…~ ~ v†J~dX ~ ‡q w ew : F ~ ÿt J~dX ~ ÿ 2l D…~ ~ ~
~ ~P bJ~dX …9† ˆq w ~ X

Z ~ X

~ X

~ ˆ0 ev : F ~ ÿt dX q~r~

~ X

~ C

…10†

~ are respectively the arbitrary pressure and velocity of the variational formuis obtained, where q~ and w lation. 2.4. Reconciling the di€erent schools of thought Eqs. (7) and (8) are the completely general reference description of an incompressible, Newtonian ¯uid. They reduce to the so-called ALE equations of Hughes et al. [4] for an instant in which spatial and referential con®gurations coincide. These simpli®ed equations should, however, not be implemented where the implementation requires evaluation about more than one point within each time step (see Section 5 for a further, indepth explanation). Under such circumstances the equations of Hughes et al. are an arbitrary Lagrangian Eulerian description in the very true sense (this is not surprising considering the equations have their origins in the arbitrarily, either Lagrangian or Eulerian programmes of Hirt et al. [3]). This fact is further borne out in observing that key energetic properties, consistant with the principle of energy conservation, are not automatically inherited by the equations of Hughes et al. in the context of more general references. ~ ÿ1 J~ ~F The momentum equations of Soulaimani et al. [9] are ¯awed as a result of the mistaken belief that r is the Piola±Kircho€ stress tensor of the ®rst kind (Soulaimani et al. [9] p. 268). Yet another problem is illustrated by rewriting the conventional incompressibility condition using the chain rule. The new incompressibility condition which arises is most certainly o~ vi o~xj ˆ0 o~xj oxi

and not

@~ vi o~xi ˆ 0: o~xj oxj

Further errors arising (e.g. J^ omitted in the ®rst term on the right hand side of the momentum equation, Eq. (10) on p. 268 of Soulaimani et al. [9]) make the use of these equations problematic.

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

225

3. The energetic implications of a deforming reference The e€ect of quantities parameterising reference deformation on key energetic properties ± nonlinear, exponential-type dissipation in the absence of forcing and long-term stability under conditions of time dependent loading ± is investigated in this section. These properties, K…v† 6 K…vjt0 † eÿ2mCt

and

lim sup K…v† 6

t!1

M2 2m2 C 2

respectively (where K ˆ 12 qkvk2L2 …X† is the total kinetic energy), are intrinsic to real ¯ows and the conventional, Eulerian Navier±Stokes equations (see Temam [10], [11], Constantin and Foias [1], Simo and ~ ref on the afore mentioned aspects of conservation of the quantity Armero [8], in this regard). The e€ect of v

2 1 1 ~ J~2 q v

2~ 2 L …X† is essentially what is being investigated, with a view to establishing a set of conditions under which the discrete approximation can reasonably be expected to inherit these self-same energetic properties. One might anticipate key energetic properties to be manifest only in instances involving a ®xed contributing mass of material, whether its boundaries be dynamic, or not. An analysis of this nature only makes sense in the context of a constant volume of ¯uid which, for simplicity, will have material limits. Inequalities of the Poincare±Friedrichs type are a key feature of any stability analysis of this nature. Gradient containing L2 terms need to be re-expressed in terms of energy. In the case of a `no slip' (v ˆ 0) condition on the entire boundary the situation is straightforward, in that it is possible to use the standard Poincare±Friedrichs inequality: there exists a constant C1 > 0 such that n

kvkL2 6 C1 krvkL2 for all v 2 ‰H01 …X†Š : The use of the classical Poincare±Friedrichs inequality is otherwise identi®ed as a major limitation, even in the conventional Navier±Stokes related analyses. The Poincare±Friedrichs inequality is only applicable in very limited instances where the value for the entire boundary is stipulated to be identically zero. For boundary conditions of a more general nature, such as those encountered in this study, in which parts of the boundary may be either a free surface or subject to traction conditions, a more suitable inequality is required (notice that subtracting a boundary velocity and analysing the resulting equation is not feasible as the equations are nonlinear). The Poincare±Friedrichs inequality does, furthermore, not hold on subdomains of the domain in question and the constant is not optimal. Further investigation (communication [7]) reveals a similar result, the so-called Poincare±Morrey inequality, holds providing the function attains a value of zero somewhere on the boundary. The Poincare± Morrey inequality states that a constant C2 > 0 exists such that n

kvkL2 6 C2 krvkL2 for all v 2 ‰H01 …X†Š : The proof of the Poincare±Morrey inequality is, however, similar to that of one of Korn's inequalities (see Kikuchi and Oden [5]). In particluar, it is non-constructive, by contradiction and the constant cannot therefore be determined as part of the proof. Viewed in this light the forthcoming inequality amounts to a speci®cation of the hypothetical constant in the Poincare±Morrey inequality for domains of a particular geometry. The particular types of geometry considered are those that arise in problems involving the motion of rigid bodies such as pebbles on the sea bed; thus a free surface is present, and the domain may be multiply connected. Inequality 1 (A New `Poincar e' Inequality). Suppose v is continuous and differentiable to first order and that v attains a maximum absolute value, c, on an included, finite neighbourhood of minimum radius Rmin about a point xorigin (as depicted in Fig. 2). If X is a bounded, star-shaped (about a point xorigin ) 4 domain in R3 , then 4

By which is meant that every point in the domain can be reached by a straight line from xorigin that does not pass outside of X.

226

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

Fig. 2. A ®nite neighbourhood of minimum radius Rmin about a point xorigin .



1 …Rmax ÿ Rmin †…R3max ÿ R3min † 2 kvkL2 …X† 6 krvkL2 …X† ‡ kckL2 …X† ; 3Rmax Rmin where Rmax is the distance to the farthest point in X from xorigin . Proof. Consider the change to spherical coordinates ; r sin h sin / ÿ xorigin ; r cos h ÿ xorigin † vi …r; h; /† ˆ vi …r sin h cos / ÿ xorigin 1 2 3 centred on xorigin . Suppose the radial limits of the domain and neighbourhood are denoted Rb …h; /† and Ra …h; /† respectively. By the fundamental theorem of integral calculus !2 Z r  2 o vi …n; h; /† dn vi …r; h; /† ÿ vi jRa …h;/† ˆ Ra …h;/† or Z ˆ

r Ra …h;/†

Z 6

r Ra …h;/†

Z 6

Rmax Rmin

1 o vi …n; h; /† dn n n or

1 dn n2 1 dn n2

Z

Z

…Rmax ÿ Rmin † ˆ Rmax Rmin ˆ

where

Vi …h; /† ˆ

Z

r Ra …h;/† Rb …h;/†

Ra …h;/†

Z

!2

o vi …n; h; /† or

!2

o vi …n; h; /† @r

Rb …h;/† Ra …h;/†

n2 dn

…by Schwarz inequality†

!2

o vi …n; h; /† or

 n2 dn …for r 2 X† !2 n2 dn

…Rmax ÿ Rmin †  Vi …h; /†; Rmax Rmin

Rb …h;/†

Ra …h;/†

o vi …n; h; /† or

!2 n2 dn:

 outside the neighbourhood (angular extent being Integrating this result over that part of X Ha …/† 6 h 6 Hb …/† and Ua 6 / 6 Ub )

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

Z

Ub Ua

Z

Hb …/†

Ha …/†

Z

Rb …h;/† 

Ra …h;/†

vi …r; h; /† ÿ vi jRa …h;/†

6

…Rmax ÿ Rmin † Rmax Rmin

…Rmax ÿ Rmin † 6 Rmax Rmin

Z Z

2

Ub Ua Ub Ua

r2 sin h dr dh d/

Z

Ha …/†

Z

Z

Hb …/†

Hb …/†

Ha …/†

6

‡

…Rmax ÿ

3Rmax Rmin

r2

1 sin2 h

@ vi o/

ÿ

Rb …h;/†

Ra …h;/†

Vi …h; /†r2 sin h dr dh d/

Z  Vi …h; /†

…Rmax ÿ Rmin †…R3max ÿ R3min † 6 3Rmax Rmin Rmin †…R3max

227

R3min †

Z

Ub

Z

Ub

Rmin

Hb …/†

Ha …/†

Ua

Z

Rmax

Z

Hb …/†

Ha …/†

Ua

Z

 2

r dr sin h dh d/ Rb …h;/†

Ra …h;/†

Z

Rb …h;/†

Ra …h;/†

o vi @r

!2

2 vi 4 o @r

r2 sin h dr dh d/ !2

1 ‡ 2 r

@ vi @h

!2

!2 3 5 r2 sin h dr dh d/

…Rmax ÿ Rmin †…R3max ÿ R3min † ˆ 3Rmax Rmin

Z

Ub

Ua

Z

Hb …/†

Ha …/†

Z

Rb …h;/†

Ra …h;/†



   r vi
r vi r2 sin h dr dh d/:

Changing back to the original rectangular coordinates and de®ning vjbndry to be a radially constant  jRa …h;/† for r ˆ Ra …h; /†, function throughout X which takes the values of v Z  Z 2 …Rmax ÿ Rmin †…R3max ÿ R3min † vi …x† ÿ vi jbndry dX 6 …rvi …x††
…rvi …x††dX 3Rmax Rmin X X where X is X excluding the neighbourhood. Summing over i, Z  Z    …Rmax ÿ Rmin †…R3max ÿ R3min † v ÿ vjbndry
v ÿ vjbndry dX 6 …rv† : …rv†dX: 3Rmax Rmin X X Making use of either the Cauchy±Schwarz or triangle inequality,  2

…Rmax ÿ Rmin †…R3max ÿ R3min †

2 6 krvkL2 …X † ; kvkL2 …X † ÿ vjbndry 2 L …X † 3Rmax Rmin and remembering that sup j vjRa …h;/† j 6 c,  1 …Rmax ÿ Rmin †…R3max ÿ R3min † 2 kvkL2 …X † 6 krvkL2 …X † ‡ kckL2 …X † : 3Rmax Rmin Consider the terms kvkL2 and kckL2 . Comparing these terms under circumstances of sup jvj 6 c leads to the conclusion that the inequality holds over the neighbourhood and that the inequality is therefore una€ected when the domain of integration is extended to include the neighbourhood. Of course, the radial extension of vjbndry can be used in place of c in instances where inclusion of the neighbourhood is not required. This inequality is similar to the Poincare±Friedrichs inequality when c ˆ 0, but is extended to a geometrical subclass of domains which have free and partly non-zero boundaries. It has a further advantage in that the constant is an order of magnitude more optimal when used under the `no slip' Poincare±Friedrichs condition (under such conditions the domain can always be deconstructed into a number of subdomains in which Rmin ˆ 13 Rmax ). The Poincare±Friedrichs inequality is a special case of the above inequality. The necessary lemma (below) follows naturally from the above inequality.

228

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

Lemma 1 (Deviatoric stress term energy). The kinetic energy satisfies the bound

2

C ~ 1

~ ~ 2 K…~ v† 6 D…~ v† J

2 ~ ; where C is related to the constant in Inequality 1; C > 0: q L …X† Proof. If, in particular, vjbndry ˆ 0 in Inequality 1, krvkL2 …X† p C 2 2 1 C 2kvkL2 …X† 6 kD…v†kL2 …X† : kvkL2 …X† 6

(The relationship between D and rv arises in the context of the original equations involving div r. It is because ÿ
Dij;j ˆ 12 vi;jj ‡ vj;ij ÿ
ˆ 12 vi;jj ‡ vj;ji …changing the order of differentiation† ˆ 12 vi;jj

…divv ˆ 0 by incompressibility†;

~ assuming, of course, that v is continuous and di€erentiable to ®rst order.) Rewriting in terms of X

2

C ~ ~ v†J~12 K…~ v† 6 D…~

2~: q L …X† The following lemma is vital to the deforming reference analysis in particular. It will form the basis to the next lemma and another concerned with the time discrete analysis. Lemma 2 (Basic to Lemmas 3 and 5). The relation D E D E D   E e u† F ~ ÿ1 w ew : F ~ ÿt ; v ~ ÿ1 w e v† F ~ J~ ~ J~ ~; …r~ ~u r~ ~u; …r~ ~J~ ˆ ÿ v ÿ 2 2 ~ L …X†

~ L …X†

~ L2 …X†

is valid for

n o ~ ÿt N ~
w ~ : ~ 2W ˆ w ~ :w ~ ˆ 0 or F ~ ˆ 0 on C w

e v† F ~ ÿ1 w ~ J~: Proof. Consider ~u
…r~ ~ k J~ ˆ ÿ~ ~ k J~ ÿ u~i v~i …F~jkÿ1 w ~ k J~†; j ‡ …~ ~ k J~†;j u~i v~i;j F~jkÿ1 w ui v~i F~jkÿ1 w ui;j v~i F~jkÿ1 w ~ k J~†;j by the product rule. In the terms arising from …F~jkÿ1 w o o~xj ÿ1 F~jk; j ˆ o~xj oxk o @~xj ˆ oxk o~xj ˆ0

and …order of differentiation ~…x; t† interchangeable for x continuous†

  o ox o~xj J~; j F~jkÿ1 ˆ det o~xj x oxk o~ o ox ˆ det oxk o~ x ˆ 0:

The latter result becomes apparent when considering that all terms which arise from the di€erentiation of the determinant contain factors of the form o2 xi =oxj o~xk . These vanish as continuity once again a€ords the order of di€erentiation interchangeable. Thus ~ k J~ ˆ ÿ~ ~ k J~ ÿ u~i v~i F~jkÿ1 w ~k;j J~ ‡ …~ ~ k J~†;j : ui;j v~i F~jkÿ1 w ui v~i F~jkÿ1 w u~i v~i;j F~jkÿ1 w

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

229

~ and applying the divergence theorem, Integrating over the domain X D E D E D   E e u† F ~ ÿ1 w ew : F ~ ÿt ; v e v† F ~ ÿ1 w ~ J~ ~ J~ ~; …r~ ~u r~ ~u; …r~ ~ J~ ˆ ÿ v ÿ ~ ~ ~ L2 …X† L2 …X† L2 …X† D   E ~
w ~ ÿt N ~ J~ ~ F ‡ ~ u; v : 2

…11†

~ L …C†

The condition of this lemma dictates the manner in which the reference must deform to ensure that the equation will inherit the desired energetic properties. This lemma is crucial to the deforming reference analysis. The lemma immediately below will facilitate the elimination of the convective energy rate in the forthcoming analysis. Lemma 3 (Convective energy rate). The relation * + D   E ~ 1 @ J ÿ1 ref e v† F ~ ~; v ~ ~ÿv ~ J~ ~; …r~ v ˆÿ q v ÿq v ~ L2 …X† 2 @t

~ L2 …X†

is valid in instances where a purely Lagrangian description is used to track free boundaries and/or boundaries are of a fixed impermeable type. ~ÿv ~ref vanProof. In instances where a purely Lagrangian description is used to track free boundaries, v ÿt ~
v ~ N ~ at ®xed impermeable boundaries. The condition at the boundary for Lemma 2 is ishes, as does F therefore satis®ed. Thus the term D  D  D    E  E   E e v† F ~ ÿ1 v e v e v† F ~ ÿ1 v ~ ÿt ; v ~ÿv ~ ref J~ ~ÿv ~ ref J~ ~; …r~ ~ r ~ÿv ~ref : F ~; …r~ ~ J~ ˆ q v ‡ q v ÿq v ~ ~ ~ L2 …X† L2 …X† L2 …X† D   E D   E e v† F ~ ÿ1 v e vref : F ~ ÿt ; v ~ÿv ~ ref J~ ~; …r~ ~ r~ ~ J~ ˆq v ÿq v ~ L2 …X†

~ L2 …X†

…by incompressibility† 1 D  e ref ~ ÿt  ~E ~ r~ ~J ˆ ÿ q v v : F ;v ~ L2 …X† 2 * + 1 @ J~ ~; v ~ ˆ ÿ q v ; 2 @t 2 ~ L …X†

~ ÿt

e vref : F ) in the same vein as J _ 0 ˆ J0 div v (the kinematic result used since oJ~=ot ˆ J~ div vref (which is J~r~ earlier). This lemma concludes the preliminaries required for the deforming reference energy analysis. 3.1. Exponential dissipation in the absence of forcing The issue of whether nonlinear, exponential-type dissipation in the absence of forcing is a property intrinsic to the deforming reference description is resolved as follows. Theorem 1. (Exponential dissipation in the absence of forcing) A sufficient condition for the completely general referential description to inherit nonlinear, exponential type energy dissipation ~ v† 6 K…~ ~ vj † eÿ2mCt K…~ t0 1 ~ v†  1 qk~ (where K…~ vJ~2 k2L2 …X† ~ ) in the absence of forcing (an intrinsic feature of real flows and the conventional, 2 Eulerian Navier±Stokes equations) is that the reference moves in a purely Lagrangian fashion at free boundaries.

~ for w ~ Proof. The ®rst step towards formulating an expression involving the kinetic energy is to substitute v in the variational momentum equation, Eq. (9). Then

230

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

*

+ o~ v~ ~; J q v @t

D E e v; F ~ ÿt J~ ˆ p~r~ ~ L2 …X†

~ L2 …X†

D E ~; ~ ‡q v bJ~

~ L2 …X†

D E ~ v†; D…~ ~ v†J~ ÿ 2l D…~

D E ~N ~ ~; P ‡ v

~ L2 …C†

~ L2 …X†

D   E e v† F ~ ÿ1 v ~ÿv ~ref J~ ~ ; …r~ ÿq v

~ L2 …X†

:

…12†

The term containing the pressure, that is D E ev : F ~ ÿt J~ ; p~r~ 2 ~ L …X†

vanishes as a result of incompressibility (Eq. (8)). The order of integration and di€erentiation are interchangeable (limits are time-independent in the reference which tracks the free boundary perfectly ± a description which becomes fully Lagrangian at boundaries was stipulated). Eq. (12) can be rewritten 0 1 * +

2

D   E ~

1 2 1 @d @ J 1 ref ~

e v†F ~ ÿ1 v ~ A ~ ~ 2 2 ~ ~ ~ ~ ~ ~ r~ ÿ v ÿ q v ; … ÿ v ; v ˆ ÿ 2l D J q v J J 2

2~ ~ L …X† 2 dt 2 ~ @t L …X†

~ L2 …X†

L …X†

D E ~ ; ~bJ~ ‡q v

~ L2 …X†

D E ~N ~ ~; P ‡ v

~ L2 …C†

as a result. The conditions of Lemma 3 are also satis®ed for a description which becomes fully Lagrangian at free boundaries and an expression

D E D E ~ v†

1 2 dK…~ ~ ~N ~ ~ J~2 ~ ~ ~ ‡ q v ; b J ‡ v ; P ˆ ÿ 2l D

2~ ~ ~ L2 …X† L2 …C† dt L …X†

2

~12 ~J is therefore obtained, where K~ ˆ 12 q v is the total kinetic energy. Using Lemma 1 ~ L2 …X† D E D E ~ v† dK…~ ~N ~ ~ v† ‡ q v ~; ~ ~; P 6 ÿ 2mC K…~ bJ~ 2 ‡ v : …13† ~ ~ L …X† L2 …C† dt Eq. (13) has a solution of the form ~ vjt † eÿ2mCt K~ 6 K…~ 0 ~N ~ ˆ 0), providing a purely Lagrangian description is used at in the absence of forcing (`no forcing' ) ~ bˆP free boundaries. A nonlinear, exponential-type energy dissipation in the absence of forcing is therefore an intrinsic property of the completely general reference description. This contractive ¯ow property is also an intrinsic property of the conventional Navier±Stokes equations. 3.2. Long-term stability under conditions of time-dependent loading The formulation of suitable load and free surface bounds is necessary before the issue of long-term stability (L2 -stability) under conditions of time-dependent loading can be resolved. The following lemma facilitates the formulation of load and free surface bounds. Lemma 4 (Force, free surface bounds). The inequality D E ~; ~bJ~ q v

~ L2 …X†

! !

2

2

2

2

1

1 mC 1



~ ~2 ~~ ~ J~2 ~ q q 6 ‡

v

bJ 2 ~ ‡ PN L2 …C†

2 ~ ‡ v ~ ~ ~ L2 …C† L2 …C† 2 2mC L …X† L …X†

D E ~N ~ ~; P ‡ v

holds where mC is a constant, mC > 0.

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

231

Proof. In terms of the Cauchy±Schwarz inequality,



D E

1

1

v

~ ~2 ~2 ~ ~; ~bJ~ 6 J v

bJ 2 ~

~ L2 …X† ~ L2 …X† L …X†

2

2

mC ~1 1 1

~ 2 ~ v J 6 ‡ bJ~2



2 ~ for mC > 0 2 2mC ~ L2 …X† L …X† by Young's inequality. Similarly,

D E mC 1

2

~ ~ 2 ~N ~ ~ ~; P 6 ‡ v P N v

2 ~ for mC > 0: ~ ~ L2 …C† L2 …C† L …C† 2 2mC This done, the mathematical machinery necessary to the long-term stability analysis is in place. Theorem 2 (Long-term stability). A sufficient condition for the completely general referential description to inherit the property of long-term stability 2 ~ v† 6 M lim sup K…~ t!1 2m2 C 2

under conditions of time-dependent loading (an intrinsic feature of real flows and the Navier±Stokes equations), where this time-dependent loading and the speed of the free surface is bounded in such a way that

2

2

1

~ ~ 2 2 2 ~bJ~2 ~ ‡ P N ‡ m C 6 M 2; q

v

2 2 ~ L2 …X†

~ L …C†

~ L …C†

is that the description becomes purely Lagrangian at free boundaries. Proof. Using Lemma 4 in Eq. (13), then applying the above bound, 2 ~ v† dK…~ ~ v† 6 M : ‡ mC K…~ 2mC dt

Using the Gronwall lemma (see Hirsch and Smale [2]) leads to the di€erential inequality ~ v† dK…~ M 2 ÿmCt e ; 6 2mC dt which, when solved, yields 2 ÿ
~ v† 6 eÿmCt K…~ ~ vjtˆt † ‡ 1 ÿ eÿmCt M : K…~ 0 2m2 C 2

This in turn implies 2 ~ v† 6 M : lim sup K…~ t!1 2m2 C 2

The preceding analyses lead to natural notions of nonlinear dissipation in the absence of forcing and long-term stability under conditions of time-dependent loading for the analytic problem. 4. The energetic implications of the time discretisation This section is concerned with establishing a class of time discretisations which inherit the self-same energetic properties (nonlinear dissipation in the absence of forcing and long-term stability under conditions of time dependent loading) as the analytic problem, irrespective of the time increment employed. In this section a generalised, Euler di€erence time-stepping scheme for the completely general reference equation is formulated and the energetic implications investigated in a similar vein as the analytic equations in the previous section.

232

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

This stability analysis is inspired by the approach of others to schemes for the conventional Navier± Stokes equations. The desirability of the attributes identi®ed as key energetic properties is recognised and they have been used as a benchmark in the analysis of various of the conventional, Eulerian Navier±Stokes schemes by a host of authors. Related work on the conventional, Eulerian Navier±Stokes equations can be found in a variety of references, for example Temam [11] and Simo and Armero [8]. The analyses presented here are extended, not only in the sense that they deal with the completely general reference equation, but also in that non-zero boundaries, so-called free boundaries and time-dependent loads are able to be taken into account (the former two as a consequence of Inequality 1). The ®ndings of this work have profound consequences for the implementation of the deforming reference equations. It is signi®cant that many algorithms used for long-term simulation do not automatically inherit the fundamental qualitative features of the dynamics. 4.1. A generalised time-stepping scheme An expression for a generalised Euler di€erence time-stepping scheme can be formulated by introducing an `intermediate' velocity ~jt‡Dt ‡ …1 ÿ a† v ~ jt ~ n‡a  a v v

for

a 2 ‰0; 1Š

…14†

~jt and v ~jt‡Dt are the solutions at times t and t ‡ Dt to the variational momentum equation (Eq. (9)) where v respectively, Dt being the time step. It is in this way that a generalised time-discrete approximation of the momentum equation E D E D E qD e w; F ~ ÿt J~n‡a ~ w†; D…~ ~ vn‡a †J~n‡a ~ ; …~ ~n †J~n‡a r~ ~ ˆ p ÿ 2l D…~ w vn‡1 ÿ v n‡a ~ n‡a † ~ n‡a † ~ n‡a † L2 …X L2 …X L2 …X Dt D   E e vn‡a †F ~ ÿ1 v ~ ~ ; …r~ ~ ref ÿq w n‡a Jn‡a n‡a ~ n‡a ÿ v 2 D E ~; ~ ‡q w bn‡a J~n‡a

~ n‡a † L …X

~ n‡a † L2 …X

D E ~ n‡a N ~ n‡a ~; P ‡ w

~ n‡a † L2 …C

…15†

~ n‡a , is derived, where h : iL2 …X~ n‡a † denotes the L2 inner product over the deforming domain at time t ‡ aDt. C ~ n‡a , P ~ n‡a , and ~ ~n‡a , J~n‡a , D bn‡a are likewise de®ned to be the relevant quantities evaluated at time t ‡ aDt. F It will presently become clear that it makes sense to perform the analyses for the time-discrete equation in the context of divergence free rates of reference deformation only. This is since relevant energy terms are not readilly recovered from the time-discrete equations for deforming references in general. This investigation is accordingly restricted to a subclass of reference deformations in which `reference volume' is conserved. This is for reasons of expedience alone and the subclass of deformations is thought to be representative. Assumption 1. The assumptions J~n ˆ J~n‡a and J~n‡1 ˆ J~n‡a are made so that



1 2

1 2 1 1 2

~ vn † ˆ q v ~ n J~2 ~ n J~ ˆ q v K…~ 2 n L2 …X~ n † 2 n‡a L2 …X~ n‡a † and

2

2

1 1 1 1 2 2

~ ~ ~ ~ n‡1 Jn‡1 ~ n‡1 Jn‡a ˆ q v K…~ vn‡1 † ˆ q v

2~ 2 2 ~ n‡1 † L2 …X L …Xn‡a † (by Eq. (14) and since the volume of material over which integration is being performed is constant). ~ ~ Remark. Notice that Jn‡1DtÿJn ˆ J~ div vref n‡a ; the discrete form of

div vref n‡a ˆ 0

oJ~ ˆ ot

J~ div vref , can consequently be rewritten as

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

233

under the conditions of the above assumption. It is for the practical expedience afforded by Assumption 1 alone that this analysis is limited to instances in which div vref n‡a ˆ 0. The following lemma will establish that the rate of energy change associated with the convective term vanishes as a result of the assumption. Lemma 5 (Discrete convective energy rate). The discrete convective term D   E ref ~ ÿ1 v e vn‡a †F ~n‡a ~ ~ ~ ; …r~ ÿ v J ÿq w n‡a n‡a n‡a 2 ~ n‡a † L …X

vanishes under circumstances of div vref n‡a ˆ 0 and a purely Lagrangian description is used at free boundaries (alternatively boundaries are of the fixed, impermeable type). D E ~ ÿ1 w e F ~ J~ Proof. The operator
; …r
† is skew-symmetric for ~ L2 …X† o n ~
w e w† : F ~ ÿt ˆ 0 on X; ~ w ~ ÿt N ~ ~ ˆ 0 on C ~ 2W ˆ w ~ : …r~ ~ ˆ 0 or F w ~ n‡a ÿ v ~ ref by Eq. (11). In instances where a purely Lagrangian description is used to track free boundaries v n‡a ÿt ~ ~ ~ n‡a vanishes. The condition at the boundary is vanishes. At ®xed, impermeable boundaries Fn‡a Nn‡a
v therefore satis®ed, under all of the afore mentioned circumstances. Apply the stipulated condition ~ ÿt ˆ 0 and set w e vref : F ~ ˆv ~n‡a ÿ v ~ref r~ n‡a etc. Remark. Recall that in the investigation of the analytic problem, a term arising from the manipulation of the acceleration containing term (the term containing the rate of change of the Jacobian) cancelled with the convective energy. It is therefore not surprising that assumptions pertaining to the acceleration containing term (in particular to the rate of change of the Jacobian) in the discrete problem will, once made, also be necessary for the corresponding discrete convective energy term to vanish (reffering to the div vref ˆ 0 condition of Lemma 5). This is a good prognosis for the energetic behaviour of the discrete problem in circumstances of reference deformations excluded by Assumption 1. This concludes the preliminaries required for the analysis of the time-discrete equation. 4.2. Nonlinear dissipation in the absence of forcing The following analysis establishes a class of time-stepping schemes which exhibit nonlinear dissipation in the absence of forcing regardless of the time increment employed. Theorem 3 (Nonlinear dissipation in the absence of forcing). Suppose that the description is pure Lagrangian at any free boundaries and that the deformation rate of the reference is divergence free. A sufficient condition for the kinetic energy associated with the generalised class of time-stepping schemes to decay nonlinearly

2

1 2

~ ~ ~ ~ vn † 6 ÿ Dt 2l D…~ vn‡a †Jn‡a K…~ vn‡1 † ÿ K…~

~ n‡a † L2 …X

in the absence of forcing and irrespective of the time increment employed, is that the scheme is as, or more, implicit than central difference. That is a P 12: ~n‡a in terms of Eq. (14) and subtracting, the result ~ n‡1 and v Proof. Expressing the intermediate velocities v 2   ÿ
~ n‡a ˆ a ÿ 12 v ~ n‡1 ÿ v ~n ‡ v ~ n‡1 …16† v 2 is obtained. The ®rst step towards formulating an expression involving the kinetic energy of the gener~ n‡a . By further substituting alised time stepping-scheme (15) is to replace the arbitrary vector, w, with v

234

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

Eq. (16) into Eq. (15) and eliminating the pressure containing term in a similar manner to that in Theorem 1, an expression involving the di€erence in kinetic energy over the duration of a single time step is obtained. ~ÿv ~ ref vanishes in instances where a purely Lagrangian description is used to track free The vector v ~ ~ ÿt N ~n‡a vanishes where boundary conditions are of a ®xed impermeable boundaries. The quantity F n‡a n‡a
v type. The condition at the boundary for Lemma 5 is therefore satis®ed. Incompressibility and a restriction on reference deformations to those for which div vref n‡a is zero ensure that the remaining Lemma 5 condition is satis®ed. The equation,

2

 1

2

1 ÿ
 1 2 2

~ ~ ~ ~ ~ ~ ~ vn † ˆ ÿ q a ÿ 2 vn‡1 ÿ vn Jn‡a ÿ Dt 2l D…~ vn‡a †Jn‡a K…~ vn‡1 † ÿ K…~

D E ~ n‡a ; ~ ‡ Dtq v bn‡a J~n‡a

~ n‡a † L2 …X

~ n‡a † L2 …X

D E ~ n‡a ~ n‡a N ~ n‡a ; P ‡ Dt v

~ n‡a † L2 …X

~ n‡a † L2 …C

;

is then obtained. Since it is assumed that there is no forcing,

2  1

2

1
 1 2 2 ~

~ ~ ~ ~ ~n‡1 ÿ v ~n Jn‡a vn † 6 ÿ q…a ÿ 2 v ÿ Dt2l D…~ vn‡a †Jn‡a K…~ vn‡1 † ÿ K…~

~ n‡a † L2 …X

~ n‡a † L2 …X

…17†

:

Thus the kinetic energy inherent to the algorithmic ¯ow decreases nonlinearly in the absence of forcing, irrespective of the time increment employed and for arbitrary initial conditions provided that a P 12

and

div vref n‡a ˆ 0:

The former requirement translates directly into one specifying the use of schemes as, or more, implicit than central di€erence. Only for descriptions which become fully Lagrangian at free boundaries can it be guaranteed that energy will not be arti®cially introduced by way of the reference. Remark. Notice (by Lemma 1) that for a ˆ 12 an identical rate of energy decay ~ vn † ~ vn‡1 † ÿ K…~ K…~ ~ vn‡a † 6 ÿ 2mC K…~ Dt is obtained for the discrete approximation as was obtained for the equations. 4.3. Long-term stability under conditions of time-dependent loading This second part of the time-discrete analysis establishes a class of time stepping schemes which exhibit long-term stability under conditions of time dependent loading irrespective of the time increment employed. The following lemma is necessary to the analysis and is concerned with devising a bound for the energy at an intermediate point in terms of energy values at either end of the time step. Lemma 6 (Intermediate point energy). The following bound applies   ~ vn † ~ vn‡1 † ‡ …1 ÿ a† 1 ÿ a ÿ a K…~ ~ vn‡a †  a…a ÿ c ‡ ac†K…~ K…~ c where c is some constant, c > 0. Proof. By Young's inequality



1

1 2 2

v

~ ~ ~ ~ v J J

n‡1 n‡a 2 ~ n n‡a L …Xn‡a †

~ n‡a † L2 …X

6

2  c 

1 2

v

~ ~ J n‡1 n‡a

2 2

~ n‡a † L …X

 ‡

2 

1 1 2

v

~ ~ J n n‡a

2c ~ n‡a † L2 …X

…18†

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

235

~ vn‡a † explicitly in terms of the `intermediate' velocity de®nition, Eq. (14), leads to for c > 0. Writing K…~ D E ~ vn‡a † ˆ a2 K…~ ~ vn‡1 † ‡ …1 ÿ a†2 K…~ ~ vn † ‡ 2a…1 ÿ a† v ~ n‡1 ; v ~n J~n‡a K…~ 2 ~ n‡a † L …X



1 2 ~ 2 ~ 2

~ ~ n‡1 Jn‡a P a K…~ vn‡1 † ‡ …1 ÿ a† K…~ vn † ÿ 2a…1 ÿ a† v

~ n‡a † L2 …X



1 2

v

~

~ n Jn‡a

~ n‡a † L2 …X

h i ~ vn‡1 † ‡ …1 ÿ a† …1 ÿ a† ÿ a K…~ ~ vn † P a‰a ÿ …1 ÿ a†cŠK…~ c using Eq. (18). The optimal choice of the constant c is established farther on. The following theorem establishes a class of time-stepping schemes which exhibit long-term stability under conditions of time-dependent loading regardless of the time increment employed. Theorem 4. (Long-term stability) Suppose that the description is pure Lagrangian at any free boundaries and the rate at which the reference is deformed is divergence free. A sufficient condition for the algorithmic flow to exhibit long-term stability under conditions of time-dependent loading (intrinsic to real flows and the Navier± Stokes equations), assuming this time-dependent loading and the speed of the free surface is bounded in such a way that

2



2

1

~ ~ 2

2 2 2 ~

~ ~ v N ‡ P ‡ m C q bn‡a Jn‡a



2 ~ 6 M2 n‡a n‡a n‡a

2~ ~ n‡a † L2 …C L …Cn‡a † L …Xn‡a †

is a > 12: Proof. Substituting Lemma 4 and Lemma 1 into Eq. (17), applying the above bound and choosing a P one obtains

1 2

2 ~ vn † ~ vn‡1 † ÿ K…~ K…~ ~ vn‡a † 6 M : ‡ mC K…~ 2mC Dt

From this point on the argument used is identical to that of Simo and Armero [8] for the conventional, Eulerian Navier±Stokes equations. Substitution of Lemma 6 leads to a recurrence relation, ~ vn‡1 † 6 K…~

1 ÿ mC…1 ÿ a†…1 ÿ a ÿ ac†Dt M 2 Dt ~ vn † ‡ : K…~ 2mC ‰1 ‡ mCa…a ÿ 1 ‡ ac†DtŠ 1 ‡ mCa…a ÿ c ‡ ac†Dt

Using this recurrence relation to take cognisance of the energy over all time steps, ~ vn‡1 † 6 K…~



1 ÿ mC…1 ÿ a†…1 ÿ a ÿ ac†Dt 1 ‡ mCa…a ÿ c ‡ ac†Dt

n

~ v0 † K…~

k nÿ1  X …1 ÿ mC…1 ÿ a†…1 ÿ a ÿ ac†Dt† M 2 Dt ‡ 2mC ‰1 ‡ mCa…a ÿ c ‡ ac†DtŠ kˆ0 1 ‡ mCa…a ÿ c ‡ ac†Dt is obtained. An in®nite geometric series which converges so that  ÿ1 …1 ÿ mC…1 ÿ a†…1 ÿ a ÿ ac†Dt† M 2 Dt 1ÿ 2mC ‰1 ‡ mCa…a ÿ c ‡ ac†DtŠ 1 ‡ mCa…a ÿ c ‡ ac†Dt M2   ˆ 2mC mCa…a ÿ c ‡ ac† ‡ mC…1 ÿ a†…1 ÿ a ÿ ac†

~ vn‡1 † 6 lim sup K…~

n!1

…19†

236

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

results, providing the absolute ratio of the series is less than unity. That is 1 ÿ mC…1 ÿ a†…1 ÿ a ÿ ac†Dt 1 ‡ mCa…a ÿ c ‡ ac†Dt < 1: Therefore either

 a ÿ1 ÿ mCa…a ÿ c ‡ ac†Dt < 1 ÿ mC…1 ÿ a† 1 ÿ a ÿ Dt c

or

 a 1 ÿ mC…1 ÿ a† 1 ÿ a ÿ Dt < 1 ‡ mCa…a ÿ c ‡ ac†Dt c

…20†

in order for the bound to exist. Notice, furthermore, that for this desired convergence to be unconditional (regardless of the time increment employed) requires a ÿ c ‡ ac P 0: The denominator in the series ratio might otherwise vanish for some value of Dt. For a 2 (21) together imply

1 2



…21†

; 1 Eq. (20) and

…1 ÿ a† a < c6 a …1 ÿ a† which in its turn implies …1 ÿ a† a < : a …1 ÿ a† The choice of the parameter a > 12 therefore leads to an in®nite geometric series which forms the desired upper bound. The minimum value of this bound occurs for c chosen according to inf

…1ÿa† a a
1 1 ˆ : mCa…a ÿ c ‡ ac† ‡ mC…1 ÿ a†…1 ÿ a ÿ ac† mC…2a ÿ 1†2

The value of this upper bound, which occurs for the choice of the parameter a > 12, is then ~ vn‡1 † 6 lim sup K…~

n!1

M2 2m2 C 2 …2a ÿ 1†2

:

In this way one arrives at a class of algorithms which are unconditionally (irrespective of the time increment employed) stable. Remark. Notice that for a ˆ 1 one obtains an identical energy bound for the discrete approximation as was obtained for the equations. 5. An `updated' approach and a simpli®ed implementation Ever burgeoning deformation gradients accumulate for a straight forward implementation of the equations. Using an `updated' approach is one way of coping with this otherwise rather daunting prospect. An `updated' approach is the result of a little, well-worthwhile lateral thinking. An `updated' approach amounts to choosing a new referential con®guration during each time step. In the case of time stepping schemes based about a single instant (e.g. the generalised class of Euler di€erence schemes investigated in Section 4) a considerably simpli®ed implementation can further be achieved by a particularly appropriate choice of con®gurations. Making the choice of a referential con®guration which coincides with the spatial con®guration at the instant about which the time stepping scheme is based allows the deformation gradient to be omitted altogether (the deformation gradient is

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

237

identity under such circumstances). For such implementations (those which require evaluation about a single point only) no error arises from the use of the equations cited in Hughes et al. [4],   @v ref ‡ rv…v ÿ v † ˆ qb ‡ div r q …22† @t div v ˆ 0:

…23†

These equations are not valid for any, arbitrary choice of reference or if the implementation requires the equation to be evaluated at more than one point within each time step (e.g. a Runge±Kutta or ®nite-element-in-time scheme). It is important to remember that in a discrete context the reference con®guration is ®xed for the duration of the entire time increment. Although the referential con®guration is hypothetical and can be chosen arbitrarily for each time step, once chosen it is static for the duration of the entire time ~ is de®ned by the deformation, step. Once the coincidence of con®gurations is ordained at a given instant, F both before and after, and must be consistant. There would seem to be no reason why one would wish to de®ne the deformation about a con®guration other than that at the instant about which the implementation is based (assuming the implementation used is indeed based about a single point e.g. a ®nite di€erence) thereby involving deformation gradients. Resolving the resulting diculties associated with the deformation gradients by means of a perturbation seems unnecessarily complicated in the light of the above reasoning.

6. Conclusions The correct equations, which describe the motion of an incompressible, Newtonian ¯uid and which are valid for a completely general range of reference deformations, are Eqs. (7) and (8). For implementations requiring the equations to be evaluated about a single instant within each time step only (e.g. ®nite differences), the deformation gradients may be assumed identity i.e. the equations of Hughes et al. [4] (Eq. 22 and 23) will suce. It is shown (as was hoped) that nonlinear, exponential-type dissipation in the absence of forcing and long-term stability under conditions of time dependent loading are properties automatically inherited by deforming reference descriptions. The single provisor is that such descriptions become fully Lagrangian at any moving boundaries (although further development of this work does allow the results to be extended to include a practically more realistic range of boundary descriptions). These properties are intrinsic to real ¯ows and the conventional, Eulerian Navier±Stokes equations. Relevant energy terms are not readily recovered from the time-discrete equations for deforming references in general. Only for divergence free rates of reference deformation which become fully Lagrangian at free boundaries could it consequently be guaranteed that energy would not be arti®cially introduced to the algorithmic ¯ow by way of the reference. The divergence free assumption was made for reasons of expedience alone and the limitations of the time-discrete analysis are consequently not expected to detract from the use of the method in any way. This is especially so when it is considered that, a term arising from the manipulation of the acceleration containing term (the term containing the rate of change of the Jacobian) cancelled with the convective energy in the investigation of the analytic problem and that assumptions pertaining to the acceleration containing term (in particular to the rate of change of the Jacobian) in the discrete problem were, once made, also necessary for the corresponding discrete convective energy term to vanish (re€ering to the div vref ˆ 0 condition of Lemma 5). If one were to be overly cautious on this basis one would be faced with the additional challenge of enforcing a fully Lagrangian description for nodes situated on any free boundaries, while deforming elements would be required to deform at a rate which is divergence free. Such a totally divergence free description may, however, not be possible. An alternative strategy would be to use a fully Lagrangian description. Both the purely Lagrangian and purely Eulerian ¯uid descriptions have divergence free rates of distortion. There are inherent problems with using certain classes of time-stepping schemes and the use of ®nite di€erence schemes more implicit than central di€erence is consequently advocated. Such di€erences exhibit the key energetic properties (nonlinear, exponential-type dissipation in the absence of forcing and long-term

238

S.J. Childs / Comput. Methods Appl. Mech. Engrg. 180 (1999) 219±238

stability under conditions of time dependent loading) irrespective of the time increment employed. A backward di€erence is the obvious choice. Calculations at time t ‡ aDt would require an intermediate mesh and associated quantities for instances in which a 6ˆ 1 (since a > 12). The author recommends a strategy in which a predominantly Eulerian description is used, where possible, for the bulk of the problem (from an eciency point of view) and the completely general reference description for the remainder is appropriate. Purely Eulerian descriptions have the advantage of a `one o€' ®nite element construction and involve none of the hazards of a badly distorted reference. Acknowledgements Daya Reddy is thanked for the loan of references and related advice. Grzegorz Luczonok and Ronald Becker are thanked for their respective opinions on the inequality. References [1] P. Constantin, C. Foias, Navier±Stokes Equations, University of Chicago Press, Chicago, IL, 1988. [2] M.W. Hirsch, Stephen Smale, Di€erential Equations Dynamical Systems and Linear Algebra, Academic Press, New York, 1974. [3] C.W. Hirt, A.A. Amsden, J.L. Cook, An arbitrary Lagrangian±Eulerian computing method of all speeds, J. Comput. Phys. 14 (1974) 227. [4] Thomas J.R. Hughes, Wing Kam Liu, Thomas K. Zimmerman, Lagrangian±Eulerian-®nite ®nite element formulation for incompressible viscous ¯ows, Comp. Meth. Appl. Mech. Engrg. 29 (1981) 329±349. [5] N. Kikuchi, J.T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, SIAM Studies in Applied Mathematics, SIAM, 1988. [6] W. Michael Lai, David Rubin, Erhard Krempl, Introduction to Continuum Mechanics, Pergamon Press, Oxford, 1978. [7] B.D. Reddy, University of Cape Town, South Africa, 1998, by communication. [8] J.C. Simo, F. Armero, Unconditional stability and long-term behaviour of tansient agorithms for the icompressible Navier±Sokes and Euler equations, Comp. Meth. Appl. Mech. Engrg. 111 (1993) 154. [9] A. Soulaimani, M. Fortin, G. Dhatt, Y. Ouellet, Finite element simulation of two and three- dimensional free surface ¯ows, Comp. Meth. Appl. Mech. Engrg. 86 (1990) 265±296. [10] R. Temam, Navier±Stokes equations and nonlinear functional analysis, CBMS-NSF Regional Conference Series in Applied Mathematics SIAM, Philadelphia, PA, 1983. [11] R. Temam, In®nite-dimensional Dynamical Systems in Mechanics and Physics, Number 68 in Applied Mathematical Sciences, Springer, Berlin, 1988.