STABILITY OF SOME EXPLICIT DIFFERENCE SCHEMES FOR FLUID-STRUCTURE INTERACTION PROBLEMS N . NEISHLOS~. M. ISRAELIand Y, KIVITY Department of Computer Science, Technion-Israel Institute of Technology, Haiia, Israel (Received6 Map 1980) Ah&act---The spectral stability theory of initial boundary value explicite finite-difkcncc schemes is used to develop a stability analysis method for problems of fluid-structure interaction. By this analysis it is t due to the interaction between the structure and fluid stability restrictions on the time step may shown be mo $ severe than commonly assumed. Four schemes of practical interest are analyzed in detail. The validity of the stability analysis is tested by somatic the effects of und~ater explosion on a submarine. The ~mp~ation~ results corroborate the prediction of the anafysis concerning the stabitity boundary.
1.lNTRODUCTlON
linearized shallow shell equations adequate for stability analysis. We consider an orthogonal curvilinear coordinate system (x’i , x;, x’&, where coordinates x’, and x; coincide with the principal curvature directions of the shell midsurface and x; is normal to the shell surface. The fluid extends from x’r=O. We suppose that for a given shell region the principal CttrVBtUrCS ?cl and K~ are constant and the diagonal components Aiand A, of the metric tensor arc equal to one. Then in various cases of practical interest, the use of simplifying assumptions followed by an appropriate orthogonal expansion leads to the following one dimensional model [2,4,6] :
In the explicit numerical integration of fluid-structure interaction problems there may be stability restrictions
on the time step which are more severe than commonly assumed. Denote by At, the maximum time step for stable integration in a fluid, by At, the corre5ponding time step in a structure, and by At the time step for stable integration in the coupled problem. Then it was found [I, 2] computationally that At is smaller than either At, or At, in certain difference schemes. The stability analysis currently available [3, 43 do not exhibit this restriction because of incomplete treatment of the interaction at the interface. In this work we consider a fluid ennui interacting with a thin shell at the boundary. The intention effects at the interface, namely (a) the coupling via the pressure, (b) the normal velocity continuity are included in the analysis while sliding between the fluid and shell is allowed. Three classes of schemes are possible. They are defined by the magnitude of the time step At required for stable integration of the coupled problem: At
d% ~+cu%=
*
-p
(2.1)
x3=0
Here w is the normal displacement into the fluid, p is the hydrodynamic pressure, t is the time, w, p, t and xj(j= 1.2.3)arenondimensionalized by L, p,u$L/a,,L respectively, where L=p&ppr, and pJ, a,, h are the density, sound velocity and thickness of the shell, pf and of are the density and sound velocity of the fluid. The frequency w of the shell depends on the mode of deformation (for stability analysis, normally the extensional mode), type of the shell theory and on the boundary conditions. Thus w is a faction of wave numbers in the reduced directions, and the parameters
= min (A”,, At&, for schemes of dass I
i We present and discuss four explicit difference schemes which demonstrate these three classes. Only schemes of class I are free from reduction of stability due to the interaction between the fluid and the shell. It should be noted here that the stability restrictions may be overcome by inte~~g the structure imp~~tly, but ,only with considerable additional ~rnpu~tio~l expense [S]. Thus the explicit approach remains attractive, especially for three dimensional problems.
pdpf,
a$lQJ.
K,L. K2L.
The flmd motion in the acoustic approximation is governed by the dimensionless wave equation in the (x,, x2, x,) coordinate system
2. EQUATIONS OF MOTION
The structural-hydrodynamic problem consists of solving the equations of structure and fluid together with contact, initial and boundary conditions. We make these assumptions: the fluid is acoustic and the structure is a thin shell. Furthermore, we use
where H,=A,(1+LK,X3).
H2=A2(1+LJc2x3~
Hr=l.
Considering the fluid motion in a thin layer near the shell surface we can assume Lqx3 Q 1 (i = 1,2). Using the previous assumptions on the shell metric and the above mentioned orthogonal expansion we
tFor correspondence: Computer Science Dept., Technion-Israel Institute of Technology. I-iaifa, Israel. 97
H.
9s
NEJSHLOSPI~
obtain from (2.2) the following: (2.3) where k is a funcuon of the wave numbers in the reduced direction and of the parameters tirL and K&. The contact condition at the interface is: (2.4) Equations (2.1), (2.3) and (2.4) together with initial and boundary conditions will be the basic model for our fluid-structure interaction problem. More details on the model and Justification of the assumptions may be found in 121. 3. NUMERICAL FORMULATION
We approximate the eqns (2.1), (2.3) and (2.4) by the following explicit finite difference scheme: 1 D;fl=(G/-k’)d; r3a) (Df + C0%“= -p$
(3b)
Dfw”= -G&j
(3c)
where D:, D:, 0: are the difference operators for second time derivatives with time step At. G: and G, are difference operators for second and first spatial derivatives with step Ax; j and n are space and time reference indices; j equal to zero denotes the shell surface. In the present work we consider a nondissipative case in which D!= 0: = 0,’ = D2. where D2 is the central difference operator. This choice gives rise to a time scheme, where the amplitude remains constant. Furthermore, the operator G: is centered and G, is any operator involving two or three spatial points. The solution procedure for the scheme (3) is: (1) Solve the shell and contact equations (3b) and (3c) for @*l and pz using the known displacements w”-‘, wn and pressures p; or p; and pi for the two or three point operators respectively. (2) Solve the fluid equation (3a) for d;” using the previouspressures~and~*‘(j=1,2,...). The schemes (3) with different types of G, will be the basic model for the present stability analysis. 4. ~ATHEMA~~AL
APPROACH
The system of eqns (3) may be considered as an approximation to an initial boundary value problem, where the eqn (3a) describes the Ruid (interior) and the eqns (3b. c) describe the shell (boundary). By this approach the spectral stability theory of initial boundary value finite difference schemes may be applied. The state of the art of this theory may be foupd in [7]. For stability analysis the su~titution
The functions dfs, At), g(r, AX), gC{r,AX) correspond to the operators D, G. 6, following the substitution above, and are given by 4 d2(s*At)=(s_
1
1)2&2 *
dr-
~_ AN=,r_-1,2AX2
13.2)
The following theorem is proved in [Z] : Tlrrorem The system of eqns (3) is stable if and only tf the following hold: (i) the equation h(r, Ax)=f‘(r. AX)
(4.3)
where h=g2 -k2
and
f=02’-
l-g,
(4.4)
has no roots r with Re(ri cOand im(r) #O, tiil the eqn (4.3) has no real roots m the Interval I -4 1 +4/k2AX2)“2, - l), (iii) the time step At must satisfy the following condition At Smin (Atf, At, At,),
(4.5)
where At, At, and At, are equal to Ati=2Ax(4+
kZAx2)-112
At,=2o-‘,
(4.61 (4.7)
2]fl- 1!2, over all roots of equation (4.4) in the interval ( - I, 0] if there are no roots of equation (4.4) in the interval (4.8) ( - LO]. In the following section we shah use this theorem to analyze the stability of the system (3) for d&rent contact operators G,. 5. STABILITY ANALYSlS
Let us consider four difference schemes where the operators G, and the corresponding functions gCand f are given in Table 1. In the following we prove that conditions (i) and (ii) of the theorem are satisfied for all operators d&ned in Table 1. The proofs are based upon a separation principle in the complex plane. Lemma 1. The functions fwhtch are defined in Table 1 and the function h satisfy the following condition: sign Im{h} = sign Im( -f}
(5.1)
for all Ax> 0 and r such that Re{r} ~0. Proof Let r=a + ib. and taken for example, case (ck wheref(r, Ax)=[2r/tr2-IlAx-2r]. We find:
and sign Im{fJ=sign Since, a<0 by assumption, (5.1) is proved. Similarly for the other three cases. Lemma 2. For real r, the functions f and A satisfy:
is made for the scalers p and w. The system (3) takes the form &+__k2 pcw2
_!?L-* 1-s
f
and
h>O.
(5.2)
for all Ax>0 and r in the interval (-(1+4/k’Ax’)‘;‘. - 1).
Stabihty of some explicit difference schemes for fluid-structure interaction problems
99
Table 1.
1
2 (r-I)Ax-2
f ,
2( r-2) (r-t )2Ax-2(r-2)
2r (r2-( )Ax-2r
2 (rtl)Ax-2
Proof. For r restricted as above we find: t
h(r, Ax) =
4 - X’>[I
-k2=0.
(rZ--
On the other hand it follows from (4.4) thatf is negative when g, is negative, but all gCin Table 1 are negative if r c - 1, and (5.2) is proved. Now we consider three classes I, II and III of schemes: Class I: Schemes where the interaction between the @id and the shell do,5 not reduce the &i&j .for all Ax Lemma 3. For functionsfdefined by schemes (a) and (b) from the Table 1 and for function h the following holds for all Ax and O&r> -1: AhtigAt,
(5.3)
6
ZiX
EX
AX
Fig. 1. Stability boundaries for the three classes of schemes. I, II and III (UJ> k).
Proef For r and Ax restricted as above f~ctions~are continuous and decrease monotoni~ly for all x from -w2 at r- - 1 to f(0, Ax)> -02 at r=O. Therefore, if the eqn (4.3) has a root ro, then 2
At, a
,;=A[,
2/w
Shell
----__ I-
2/K
but if the eqn (4.3) has no roots then Ati= At,
zj
by definition. and (5.3) is proved. To summarize. by Lemmas 1,2 and 3 we have two forward, lirst and second order schemes for the contact condition (3c), where the interaction does not reduce the stability, the stability is defined completely by the stability of the fluid and of the shell separately. The value of & separating the regions of fluid dominated stability as opposed to shell dominated stability is given by ~xx2(t$-k2)-1’2.
if
w>k.
-----
I Class Tt
---
tklss m
-
Lx
CIOSS
Ax
Fig. 2. Stability boundaries for the three classes of schemes, I1 II. and III (w
(5.4)
The typical stability boundaries for schemes of class 1 are given in Figs. 1 and 2 by full curve. Class II: Schemes where the stability is reduced fbr ull Ax as a result of the interaction between fluid and shell To this class belong the scheme with G, defined by (c) from Table 1. This case may be important for practical purposes being a second order scheme with a small truncation error.
Lemma 4. For function f. defined by the scheme (c) from the Table 1. and for function h the following holds forallsandO%rr -1
At,
(5.5)
Proof? Thefunction/(r,Ax)= 2rco2[(r2- 1)Ax -2r]- i decreases monoto~~ly from -co2 at r = - t to - CC at r=r_, =( l-41 + Ag)/Ax. On the other hand the
function h(r, Ax) increases monotomcally from - x at r= - 1 Co 4/(At# at PO. It follows that there is a single root, r,,of eqn (4.3) in the interval ( - 1.01 where jflrQ* Ax)1= f&roVAx)! is larger than both jf’( - 1, Ax)/ and @(O,Axf(. Consequently the lemma is proved. To estimate the effect of the interaction on stability we investigate some limitrng cases. We find [2] assymptotic expressions for Ar, as follows:
f’
I
Fig. 3 Comparison between mtiuenc_e oE the mteracuon, Q. and its esttmate Q +
J
Class Jli: Schemes where the stability is reduced for some Ax as a result of the interaction betweenjluul and shell
Lemma 5, For functionf’defined by scheme Id) from Table 1, and for function h, the following hoid for oar> -1
P 1 ,_+ ... , V,~Ax”I ifw=k.
Atiamin (At,. At,), for Ax
We see that Ar, approach At, and At, when Ax approaches 0 and c respectively. The typical stability boundary is given on Figs. 1 and 2 by the dashed curve. The influence of the interaction may be measured by the quantity : Q(Ax.k, ml=
min (At, Ar&-A&, min (At,. AtS) *
15.8)
An important practical case is w> k. which corresponds to a shell having higher sound velocity than the fluid. Then the region of the greatest reduction of stability is thg vicinity of Ax defined in f5.41. fn this case the inikence of the interaction on a specil?c scheme may be estimated by: _ At,-At, (5.9) Qfk’ @= At, ar=~ Using (5y and the assymptotic expressions and (5.7b) Q may be approximated by: (i(k, cu)-&,(k.
m)~(d-kz)-“Z.
66) (5.10)
In [2] we considered the response of a steel cyiindrical she%, submerged in water, to a transverse shock wave. For this case the extensional frequency is given by m=ak
where I=@&~
t 23.3 for steef shell and water), ~~~~rn~e~~a~~a~~nu~~ and R is the radius of the cylin$er. The estimate (5.10) for Q reduces to
k=~~/~~~;~~s~e
&=k-l(&l)-l?
where x is the single real root in the interval f& 21 of the equation ,f‘(O,Ax)= MO,Ax).
(5.13)
Proof: For Ax<2 the function t(r. Ax)=Zwi f(r+ l)Ax-21 decreases monotonicaiiy from -w2 at r = -1 ro f(O, Ax)< -& at r-0, while the function hlr. Ax) increases monotonically from - .x) at r = - i to h(0, Ax) at r - 0. Therefore, there exists a root r0 of eqn (4.3) in the interval ( - 1.01, if ti0, Axl<.h(O, Axt.
(5.14)
It may be easily shown ~~a~~~lly, for instancei,that eqn (5.13) always has a single real soMion Ax ~2, thus the inequality (5.14) holdrJ for Axa Ax. For Ax> 2 the functionf has a vertical asymptotic at O~=r,~=(2--A~/dxt~--t, and there exists a root r0 Such that if(rof\>wr. C sequentty, the lemma is proved. The typical stability bounckies for schemes of class III are presented at Figs. 1 and 2 by dotted curve. For an important ease kAhx= 2 (i.e. a square mesh in the f&id) the solutlon of the eqn 15.I31 is --;I;;=q1+ 14-w~t-“2. (5.151 For o--* co, G-0, and the stability behaviour of the scheme of class IIf approaches the stability b&aviour of the schemes of class II. ~syrnpto~&anaIy~ may be produced for the schemes of class III, as we have done for the schemes of class II
PI.
(5.11)
In Fig. 3 we present the curve (z, according to f5.111 arid the curve Q which is obtained n~e~~lly from @e compiete analysis. It is obscrve$ that for k>O.35 QA @vcs a very good estimate for Q, whereas f& k co.35 QA gives an overconservative estimate. In [l] it was found computationally for a problem using the scheme (c) that the time step Ar for stable int~ti~~ of the coupied problem mutt satisfy At cm min (Atp At,). This is in agreement with the above analvsis.
(5.12)
6.hPPtICATIONS
validity of the stability analysis was testedwith a nonlinear Lagrangian two-~en~o~~ program DISCO ES]. The problem considered is that of a submerged cylindrical shell subjected to a step shock wave, simulating the efFects of underwater explosions on a submarine. The code employs the scheme (a) of TabIe 1. The computational results corroborate the predictions of the analysis concerning the stability boundary. As an example, for a 5Omm-thick: steel shell cylinder with a diameter of 1Om for a 2MPa pressure wave: (acoustic The
Stability of some explicit difference schemes for fluid-structure interaction problems
l ’
-As*@AX
/
l
.. /
/AssAX
101
between the computed stability boundary of a twodimensional problem and the estimated stability boundary of a one-dimensionaI analysis, as illustrated in Section 6. These studies will be presented in a forthcoming publication. REFERENCES
AX,
mm
Fig. 4. Comparison between the computed and estimated stability boundaries, for a submerged cylindrical shell subjected to a step shock wave. region) the resuits are illustrated in Fig. 4, where ~~_2dsisatypicalmesh~iztintheshelfandAtisthe maximum allowed time step. 7. FURTHRR STUDIES
Using the same approach, we are able to produce a stability analysis, involving two dynamical equations for the shell without reduction to the one shell equation model (2.1). We iind that the results of such twodimensional anaiysis are simiiar to results which are presultai in present work. That means we can estimate the stability hehaviour of the numerical schemes by analyzing a onedimensional modd. The efkctiveness of this approach depends only on the estimate of the
highest shell frqucnce, that which is determining the numerical stability. This is coniirmed by the agreement
1. J. D. Bitzberger, Two Dimcnsio~I Analysis OJ’F&dStructure lnteracdon &y Method of Finite-LJi@rence~y~ru#~fc Ram The Fuel Tank Problem, NTIS, AD-783839 (1979). 2. H. Neishlos Dynamic Fluid-Structure Interaction, Ph.D. Thesis, Computer Science Dept., Technion-IIT, Haifa, Israel (1980). 3. T. Belytschko, H. J. Yen and R. Mullen, Mixed methods for time integration. Comput. Meth. Appl. Mech. Engng 17/W, 259-275 (1979). 4. K. C. Park, C. A. Feiippa and J. A. DeRuntz, Stabilization of staggered solution procedures for fl~d~~u~ure interaction analysis. In Corn~~ac~o~l wefts for Fluidhi~twe hFe?action ~oblems (Edited by T. ~elytschko et ol.l Vol. 26, DD.95-124 A.S.M.E. AD&& Mechanics Symposia Ser&‘AMD (1977). ** 5. T. Belytschko and R. Mullen, Stability of explicit-implicit mesh partitions in time integration. Int. J. Numer. Meth. Engng 12, 1575-1586 (1978). 6. K. For&erg, Influence of boundary conditions on the model chara&ristics of thin cylind&al shells. Al AA J. 2, 215&2157 (1964). 7. K. W. Morton, I~ti~-v~ue problems by ~~te~ff~~ and other methods. In The &ate of the-Art in Numericnf Analysis(Edited by D. A. H. Jacob). DD. 699-745. Acsdcmic.Press. New k’ork (1977). ” _* 8. Y. Kivity and D. Tzur, Study of projectile formation in lined shallow cavity charges. Proc. 4th Int. Symp. on Ballistics, Monterey, California (1972).