Fluid-structure interaction and ADINA

Compurers & Swucrures Vol Pnnted m Great Britam

17, No

5-5, pp

763-773,

1983

FLUID-STRUCTURE

0045-7949183 Pergamon

INTERACTION

$3 CKI+ 00 Press Ltd

AND ADINA

s. ZILLIACUS

David W. Taylor Naval Ship Research and Development Center, Bethesda, MD 20084. U.S.A Abstract-The theoretical basis of a fluid-structure interaction decoupling approximation for transient analysis of wave excited submerged structures is outlined. The approximation. called the Doubly Asymptotic Approximation (DAA) has been included in ADINA. Program modifications and additions are described. Results for two verification problems are given and discussed. It is concluded that the DAA augmented ADINA program gives satisfactory results

INTRODUCTION

The transient analysis of a wave loaded submerged structure is an interaction problem since loading depends upon response of the wet surface in contact with the fluid. At each time step in the analysis, that part of the pressure interaction loading depending on the movement of the wet surface must be computed[ 11. Most general-purpose structural codes have been written with intent of finding response to prescribed loading. Thus if one wishes to use a specific code such as ADINA, e.g. for fluid-structure interaction analysis, it must of necessity be augmented. It is assumed here that modeling of the exterior domain by fluid elements is prohibitively expensive. However, a decoupling scheme which concerns itself with fluid-structure boundary pressure only (the effect of the entire fluid domain on the wet surface) has been shown to be economical. It is the purpose of this paper to discuss the underlying concepts of such a scheme and its implementation into ADINA.

continuity

t +d’ 1 d= or = --‘i.vp= P ri;

P’ + PSR, P=

description

-1,” P

standard formulation, or modified formulation,

or p’+pS+pR,

standard formulation, or modified formulation,

wherep’(t) is (prescribed) pressure associated with an incident free field pressure wave; p SR= p s + p R;p s( t ) is scattered pressure due to interaction of p’ with S, for S rigid and fixed; and pR(u) is radiated and inertia pressure due to normal movement of S in the fluid. “Exact” analytical determination of fluid pressure can be effected by several techniques, among which are

The fluid in the exterior domain is assumed to be acoustic. Behavior is governed by the wave equation (WE) v2p = f j,

where V2 is the Laplacian operator; p is the transient fluid pressure; c is the sound speed in the fluid; and /i is the second Partial time derivative of D. Motion of the-structure is governed by the coupled equations of motion (EM) of the type

for an arbitrary point on S

where n^is the unit normal vector of S, positive into the fluid; p is the fluid density; C(t) is the acceleration component of S in the n^direction; and C’(t) is the component of particle acceleration normal to S associated with an incident pressure wave, positive on that part of S “illuminated” by the wave, negative in the shadow zone (standard formulation), zero in the shadow zone (modified formulation) (Fig. 1). The transient pressure p may be analyzed as follows. For a point in the fluid

FLUIMTRUCTURE INTERACTION-THEORY “Exact”

condition,

(1) Integral (time) Transform-Modal (spatial) Decomposition (ITMD) in cases where the surface Scan

D(u) = F(t, r) + G@), where D represents a spatial, temporal differential operator; F is the prescribed external forcing; G is the fluid-structure interaction forcing; and u(t, r) is the displacement component at time t, locationr. In addition to the natural or prescribed boundary conditions imposed upon the structure, continuity of normal fluid particle velocity and fluid-structure boundary normal velocity must be enforced since cavitation is prohibited by hypothesis. Thus the effective fluid particle acceleration component it normal to the fluid-structure CAS Vol 17. No 5/&J

boundary

Srructure

Spherical WaYe Fig.

surface S is given by the 763

1. Incident

(free-field) wave-structure geometry.

Plane Wave engulfment

164

s.

ZILLIACUS

be described by a separable coordinate system, for instance, by prolate spheroidal coordinates, for which simple separation constants may be found. (2) Application of Kirchhoffs Retarded Potential Integral Equation (KRPIE), which can be done numerically and be compatible with finite element analysis. KRPIE can be derived from the WE. Since p SR,ps, pR satisfy the WE and KRPIE individually, they may be solved for individually. In the standard formulation, scattering, radiation and inertia effects are combined, but are separated in the modified. Thus, in the latter, ps, may be computed for the entire time domain a priori since it depends only upon the incident wave and geometry of S. pR in the modified and pSR in the standard formulation, however, must be calculated interactively with the structure solution due to dependence upon the response of S. The continuity condition couples the EM with fluid response. It is incorporated into an integral relationship when KRPIE is used. examples of “classical” Two complete fluid-structure interaction solutions will be cited. First is the “exact” response of a spherical shell to a plane step wave obtained by Huang[2] using ITMD. After modal decomposition of displacement components in terms of Legendre functions of the colatitude angle coordinate, the coupled Lagrange equations of motion were Laplace transformed with respect to time. These two differential equations were thus converted to algebraic equations which could be readily combined into one. The pressure solution of the modal Laplace time transformed WE in spherical coordinates was obtained in terms of modified spherical Bessel functions. Multiplier coefficients in the solution were evaluated from the modal Laplace time transformed continuity condition. The Bessel functions were expanded as finite polynomials in terms of the Laplace transform parameter. Thus the transformed modal displacement components could be expressed as ratios of polynomials. Finally, a root finding program for the denominators of these ratios, Heaviside’s expansion theorem, and tabulated inverse Laplace transforms were used to obtain closed-form modal time response. Huang’s results gave a complete and detailed picture of the total transient vibratory response of the shell. The second example of a “classical” solution is the response of a submerged infinitely long cylinder to an incident transverse step pressure wave obtained by Geers[3] by means of ITMD supplemented by his Residual Potential technique. Displacement components and pressure were decomposed in terms of sines and cosines of the circumferential angle coordinate. The two coupled Lagrange equations were uncoupled by means of an approximation valid for the lower vibrational modes. The modal Laplace transformed WE yielded a solution for pressure in terms of integral order (cylindrical) modified Bessel functions of the second kind. Again, as in Huang’s spherical shell problem, the arbitrary coefficients in the solution were evaluated from the modal Laplace transformed continuity condition. The response histDetails of the derivation can be more easily explained for

the spherical problems followmg.

tory of the shell was now expressible in terms of sums of Laplace inverse transform integrals which could, however, be evaluated only by numerical techniques. At this juncture, Geers embarked upon a new course. He introduced his (first order modal) residual potential equation in lieu of the (second order) wave equation as the first step of a procedure designed to circumvent the need for numerical inversion. It should be noted that both fluid pressure or particle velocity can be expressed in terms of a velocity potential, so that the residual potential equation could equally well be written as a residual pressure equation. By combining the modal Lagrange equation, continuity condition and residual potential equation, Geers obtained a third order motion equation containing residual pressure on the right hand side. Next the modal residual potential pressure was expressed as the product of a function r,,(t) times a sum of terms containing modal incident wave pressure, shell acceleration and displacement. The functions r,(t), which in dimensionless form are independent of all physical parameters of the problem, closely resemble sums of exponentially decaying sinusoids, and can be tabulated in advance for the entire time domain. trapezoidal Finally, Geers employed and Runge-Kutta integration for simultaneous evaluation of the residual pressure and third order motion equations, respectively. The results were compared with various approximate fluid-structure decoupling schemes. Approximate description

In a subsequent paper[4] Geers, again using ITMD, extended the residual potential method from two to three dimensional analysis. He obtained the three-dimensional analogue r,,(t) of r,(t), in terms of r,(t). From this result he proceeded to find the three-dimensional cylindrical wave approximation d,,(t, a) = I$::(& PC

a) + t

where m, n are axial, circumferential

modal indices;

a, 1 are characteristic radial, axial dimensions; p = mnall; K, is the nth order modified Bessel function of the 2nd kind; and u,, pz are the mn th modal

normal particle velocity and interaction pressure. From this he obtained the corresponding matrix formulation for the decoupling approximation called the Doubly Asymptotic Approximation (DAA)t c=

$jSR+

M,-

IA~SR,

PC -

where MF is the fluid apparent added mass matrix, coupling normal motion of all discrete structure surface elements through the fluid; A is the diagonal matrix of effective surface areas; and r~ is normal effective fluid particle velocity vector. There are two special cases of DAA: (1) @ = 0 epR = 0, the case of an incident pres-

Fluid-structure sure wave scattered

765

Intelmactlon and ADINA

by a rigid, fixed body,

in which ps is caused by arrest of c’ at the structure surface. (2) ~‘=Oop’=ps=O, the case of no incident wave, but a moving-surface,

To arrive at the DAA formulation in this case, it is necessary to find the limiting values of the bracketed term for very large and small values of s, corresponding to early and late time behavior, respectively. For the early time limit, a descending series expansion of k,(s) in s is required, leading to the result a =-s,

a =-s-m

1 r; =@ =,c;“+MF-‘ApR,

C

C

L

in which pR is caused by surface motion only. The general case may be found by combining cases I and 2

for

“s %>I, C

from which

pmSR(l, a) = (standard

Insight into early and late time behavior may be gained by examination of the integrated form of DAA (for zero initial conditions): ?!SR= pc {p - M,-

a),

the PWA. For the late time limit, an ascending sion gives

series expan-

] = m + 1,

lim [

‘AcJ.‘~},

where

,+.O

whence gSR(‘) =

5

PmSR(f, a) =

‘pSR(r) dr. I-

For early time, gsR is’small approximation

pco,(t.

formulation).

so that the plane wave

the inertia approximation,

(PWA), Ref. [5], is applicable,

while later, the growth of qsR causes incompressible inertia effects to predominate. The DAA matrix formulation has the appearance of being quite general. It is asymptotically correct for early and late time. Thus, even though DAA was developed from an approximation to a solution for a specific geometry (finite cylinder with specified boundary conditions). it has been applied to diverse problems, one of which shall be described subsequently. In order to evaluate the adequacy of the DAA, Huang[6] returned to his classical sphere problem [2]. The relation between modal Laplace transformed effective particle acceleration z:,(s, a) and scattered plus radiated and inertia pressure pSR(s, a) on the spherical surface in contact with the fluid obtained by Huang is:

in which

pa

-o m+l

P SR= pee Huang

combined

c,(t, a),

A-‘M,.

these two limiting

lim [ r-large+ r-II

]=

is

+(m+l).

0

from which

pa 1.

m+l ti,(t, a) - -pmSR(t,

dmSR(r,a) = pc

[ the DAA. conditions)

or its integrated

pmSR(f, a) = pc where

cases

C tl,(t,

version

a)

(for zero initial

m+l a)

-

___

Pa

qmSR(f.

a)

1 ,

‘pmSR(r,a) dz,

qmSR(f,a) = s0

where m is the meridional modal index; a is the shell mean radius; s is the Laplace transform parameter; k, is the order m modified spherical Bessel function of the 3rd kind; and f’s, a) = Jo”fct, a) e-“‘dt, the Laplace transform of f(r, a).

Computation of fluid apparent added muss matrix Fluid added mass matrix MF is obtained by a surface integral method[7]. The structure enclosed by surface S is assumed to be immersed in an incompressible fluid of large volume V. Laplace’s equation (LE) defines the velocity potential cp in V v%$ = 0,

s. ZILLIACUS

166

subject

to the boundary

condition

I>= --fi .v4 and the radiation

on S

~=~~',,PS~=~Q,~.,[-PSBC

= -(p,,,

condition

1

in V

=z

lim Vb = 0, R- I where 11is the fluid particle velocity component along fi; fi is the unit normal of S, positive away from the fluid; V is the gradient operator; R = Ir, - rl; and ru, r are observation, source coordinate radius vectors LE is satisfied by the Green’s function R-‘.

G(r,,r)=

Now

c ‘M,c

1

= 2 L’‘M,rc.

where M, = - (p /2) {[SBC ‘I+ [SBC ‘I’;, fluid added mass matrix; and S is the diagonal matrix of surface areas. First attempt at uppro.vimuting u dingonal fluid muss mutrix Assume the full fluid mass matrix M is known and let

except at r. = r. Since G is an influence function, the solution at r. due to an arbitrary normal velocity source distribution - P(r) on Sis

T

4 (hJ = 0 - ,~

s,= i

Ml,

,=I

where n is the number Now. assume

P(r),

‘I$.,,

of wet surface

Increments

J

dS’ cos y _ ~ V(r) R2 ss

A,(ro) = - 2n B(r,) +

where the first terms on the right result from evaluation at the singularity R = 0 (r,, = r), cos y R?

=

where D, is the term i of the “equivalent” diagonal matrix D. and M,, is the diagonal term i of the full matrix M. It follows that

The above equality

assumption

is a consequence

of the in-

-(R -‘I,,,, M,, < D, < S,

ri, = ii( The surface S considering the form over each formation from

may be divided into subsurfaces. By source distribution - P to be unisubsurface and making the transintegral to matrix form

sdS’

The inequality is deduced from the oscillation of a gas filled cavity of surface area S surrounded by a very large volume of incompressible fluid. Inertia effects of the fluid on S are assumed to be represented by either D or M. If permissible natural vibrational modes are considered, then application of Newton’s second law leads to the inequality.

+ B.

,R

dS’cosy_C R?

FLUIlX3TRUCTURE INTERACTION-FINITE ELEMENT IMPLEMENTATION

’

the analogous column vector p may be eliminated and the analogous potential column vector expressed

The kinetic energy

Application

sc

at Green’s

T is given by

First Theorem

yields

Ve.VOdV.=~sb.~~dS.--~~,~V2~dY’.

Therefore +,.p dS’4

T = ; ss

= ; sJ

4p dS’+.,.

Changes in existing or~erla~~sunrl subroutines An extension of ADINA[S] was developed for computing transient hydrodynamic response of a submerged structure under steep fronted wave excitation according to the standard DAA formulation. It was necessary to modify six ADINA components, as follows: Three input tapes were added to overlay ADINA. The first stores wet surface- -wet node transformation coefficients, the second, the product matrix of inverse fluid mass matrix times projected wet surface areas, the third. projected surface areas. Overlay ADINI which administers input and solution control data was modified by allocating additional permanent memory locations at the begmning of blank common for variables required in hydrodynamic computations. Overlay LOAD which calculates load vectors was modified to the extent that the additional memory locations (defining array starting points) were passed to and through it to various subroutines.

Fluid-structure interaction and ADINA Subroutine LOADEF which handles externally applied loads inside the main time loop was modified similarly. Also, a call to the new hydrodynamic master subroutine DAACAL was added. Subroutines TODPRL and THDPRL which compute nodal point forces due to pressure on 2 and 3 dimensional element faces were altered as in the case of LOAD, above. Also. provisions for calculation of outward unit normals at the wet nodes and certain correspondence tables for wet node-global node transformations were added. New subroutines Fluid behavior was modeled according to the integrated form of the DAA (standard formulation). During program execution, transient interaction pressure p SR is computed and the load vector is modified by the new hydrodynamic subroutines at each time step. Since interaction loading is a function of fluid-structure boundary surface movement, extrapolation to the present time of nodal velocities of this surface is required. The fluid apparent added mass matrix required for DAA is generated externally, operated upon, and read into ADINA from tape. New subroutines required for producing interaction forces are as follows: DAACAL is an executive subroutine which modifies the nodal force vector to account for fluid-structure interaction pressure. DAACAL is called by LOADEF at each time step and in turn calls each of the following but the first (called only once) and the last 5. (1) DAAO is a subroutine which computes geometric relationships between incident wave and structure, and maximum incident wave pressure and decay constant for each wet node of structure. Four types of incident pressure wave are presently allowed: plane step wave, plane wave with exponential time decay, spherical step wave with l/R spatial decay, spherical wave with exponential time and decay.

l/R

(2) DAAIP computes incident wave pressures p’ and fluid particle velocity components r_”normal to structure surface S in contact with fluid for all corner nodes on S (called wet nodes) at each time step. (3) DAAEV computes normal velocity component ,i-‘(t) and @‘(t + At) (extrapolated) for all wet nodes on S at each time step. $‘(r +At) is corrected to compensate for the error @ - eE incurred at the previous time step. (4) DAATP computes transient fluid-structure interaction pressure. Normal effective fluid particle velocity vector c for the wet nodes is transformed to wet surface centroid basis. Interaction impulse vector qSR for the surface centroids is extrapolated ahead one time step. Current interaction pressure vector pSR at centroids of element surfaces on S is calculaied from the integrated form of the DAA. Interaction pressures p SRare transformed from wet surface to wet node basis. (5) WVFRNT (called by DAATP) distributes incident wave pressure p’ and normal particle velocity c’ “equitably” to wet nodes ahead of and behind the

167

incident wave front for those wet surface increments being traversed by the front. The allocation to each wet node is based upon proximity to the front. The purpose of WVFRNT is to provide realistic, smooth loading in the vicinity of the wave front. (6) VECTRA (called by DAATP) transforms node based column vectors to surface based vectors and vice versa. This subroutine is required because the fluid added mass matrix M, and projected surface area matrix A are calculated with respect to wet surfaces. Thus pressures are calculated first for wet surface centroids and then transformed to wet node basis. (7) DAAPSR (called by DAATP) solves the hydrodynamic equation. DAAPSR is called twice at each time t, to obtain, first, the interaction pressure psR(t) from e(t) and qSR(t), and second. pSR(t +At)-from $(t + At) and qSR(t + At), the last two being extrapolated. qs”(t) is corrected at each time t to compensate for-past integration error. (8) MTRMLT (called by DAAPSR) carries out the multiplication

PC’AI (4). One row at a time of the above full, square matrix is read into core. (9) DIAMLT (called by DAAPSR) performs the multiplication

P-‘Al {qj. The entire diagonal matrix above is stored in core. EXAMPLEPROBLEMS.ELASTIC CYLINDRICAL SHELLSUNDERPLANE STEP WAVELOADING “Exact” results-Single and concentric cylinders The “exact” response histories with which ADINA results are compared may be found in Huang[9]. By employing an approach similar to that he had taken in [2], Huang[lO] found the response of an infinitely long submerged cylindrical shell to plane transverse acoustic wave excitation. Subsequently, he extended the analysis to the case of two concentric fluid coupled cylindrical shells [9] Modeling Discretization is shown in Fig. 2 with three dimensional (plane strain) solid and fluid continuum elements. A consistent structure mass matrix was used. Exterior fluid domain was accounted for by integrated form of DAA with full fluid added mass matrix, except where noted. Integration techniques The Newmark integration algorithm (IOPE = 2) was used in all cases. It became apparent that the y(OPVAR( I) = DELTA) = 0.50, combination /?(OPVAR(Z) = ALPHA) = 0.25 worked quite well for this problem. This value of 7 provided no numerical damping. These “standard” Newmark parameters had caused response instability in problems executed before the introduction of subroutine WVFRNT. Results-Response histories of inner cylinder Response histories are depicted in Figs. 3-8 for the near generator (OO),the intermediate (90”) and the far

768

s. ZILLIACUS

00093

lb-sec2/in4

Incident Acoustic Plane wave pressure ,I=pc2

12 Elements/Group 12 Nodes/Element,(shells) 16 Nodes/Element'(fluid)

~j/I,,ne;;;linder

2

226

432

I

0

Steel

II

0

Water

Inner + titer Cylinder II

I

III

26

Steel II

Fig. 2. Modeling scheme for concentric cylindrical shells with 3-D continuum and fluid elements.

generator (180”), where near denotes points nearest the wave source. The figuresrefer to elastic cylindrical plane strain shells excited by a step incident plane wave of pressure p’ = pc2 as shown in Fig. 2. DISCUSSION

OF RESULTS

From Figs. 3 and 4 may be inferred the benefit of smooth loading. Nodal loading for Fig. 3 includes weighting of incident wave pressure and normal particle velocity in the vicinity of the wave front based on proximity of the front to the node in question (subroutine WVFRNT). Nodal loads for Fig. 4 are based on incident wave pressure and normal particle velocity contribution determined solely by whether or not the wave has traversed the wet node in question. This results in sudden changes in loading and associated spurious oscillatory response. Results of Fig. 5 illustrate response to step wave plus hydrostatic preload of equal pressure. Results were obtained in two stages. First, displacement response to hydrostatic loading was obtained. Second, these nodal displacements, along with zero velocity and acceleration, were used as initial conditions for dynamic execution. Radial step (prescribed) pressure and superimposed wave (incident and interaction) pressure loaded the cylinder during the dy-

namic run. Magnitude of the radial step was identical with hydrostatic pressure of the first stage. Wave arrival was delayed 12 time steps to illustrate the stable constant dynamic response. Velocity response after wave arrival is identical with that of Fig. 3 is shifted, while displacement simply builds on the initial field. From Fig. 6 can be inferred the extent to which the thinner outer shell is “transparent” with respect to the incident wave. These results from Huang[9] may be used to compare response of the inner cylinder alone with that of the inner cylinder in the presence of the outer cylinder. It is noted that the time scale for the double cylinder (Problem 2) has been shifted to remove the time lag associated with the space between the inner and outer cylinder. From Fig. 7 a comparison may be made between “exact” results and those from ADINA. Two points need be made about this ‘figure. First, the high frequency oscillation present in the ADINA response corresponds very closely in period to double the space between the concentric cylinders. Therefore, the frequency of this signal has a physical basis, even though the amplitude may be called into question and depends necessarily upon the discretization scheme. The “exact” solution, which was obtained by superposing nodal responses for

2.

-05

15

I 0

a

/

,

I

I

3

I

1

DIMENSIONLESS TIM!3T=(c/a)t

2

1

Fig. 3. Normal velocity response of cylinder.

1

I

s

t

”

m

=:

6

-1

ts

--es

5

--I

1

DIMENSIONLESS TIME

2

T=(c/a)t

3

l

Fig. 4. Normal velocity response of cylinder, subroutine WVFRNT inactive.

0

S

s.

ZILLIACUS

c

n

e/(o‘LP

h

= (O‘.L)fi

Fluid-structure

interaction and ADINA

-

0

1

2

,

3

i

SHIFTED DIMENSIONLESSTIME T'=T-(b-&/a

Fig. 7. Normal velocity response of cylinder.

, 7 (n = number of circumferential half sine waves), does not exhibit the same pronounced high frequency response. It is estimated that the high frequency component corresponds approximately to n = 33. Second, the decay in response at late time may be associated with radiation of sound associated with the high frequency signal. A comparison of outer and inner cylinder response has shown a slight frequency difference which would appear to support this contention. It is not intended that ADINA response should be considered as accurate as the “exact” results of Huang, but only to state that discrepancies can be explained as a natural consequence of finite element discretization. From Fig. 8, response using a diagonal fluid added mass matrix can be compared with that associated with a full matrix. Employment of a diagonal matrix offers the possibility of shorter execution time, important for large problems. In addition, lower natural frequencies of submerged shells may be found in a relatively straightforward manner by assigning the added mass to appropriate nodes. n=O,l,...

Computational

eficiency

Overall and unit execution times are given in Table 1. From the table it can be concluded that:

(1) The penalty for interactive loading in Problem 1 is: (6.56 - 3.75)/3.75 = 75%. (2) The penalty for using a full fluid added mass matrix in lieu of a diagonal one for Problem 2 is: (5.81 - 5.69)/5.69 = 2%. (3) The apparent decrease of UEX between Problems 1 and 2 is: (6.56 - 5.81)/5.81 = 13%. SUMMARY

AND CONCLUSIONS

Hydrodynamic (DAA) interaction loading capacity has been added to ADINA, making it possible to compute the response of a fluid immersed structure to a steep fronted acoustic wave at arbitrary incidence. Some results have been presented for the case of submerged elastic cylindrical shells under plane step wave excitation. The velocity histories for the near and far generators have been shown to agree well with those from an integral transform-modal decomposition continuum solution of Huang. Although the addition of hydrodynamic loading capacity causes an increase in execution time, this penalty does not seem too severe when the nonlinear

s. ?i!JLLIACLJS

0

1

3

2

,

5

SHIFTED DIKENSIONLESS TIMF. T'-T-(b-d/a

Fig. 8. Normal velocity response of cylmder.

Table 1. Execution ttme (EX), unit execution time (UEX) for NSTE time steps and NEQ unconstrained degrees of freedomdynamlc execution Prob. No. ___-

I

2

Fltnd mass mtrx. type _

Type of loading

5 5

NA Full Ftlll

Interactive

I. 8 8

Diagonal Full

Figure

5

WEX = EX/(NSTE*NEQ)

nature of the interaction

formulation

EX

UEX*

isec)

Omioper)

192 192 192

8.64 47.86 62 97

3.75 6.56 6 56

432 432 ~~_.

255.45 261.17

NSTE

NEQ

Interactive

12 38 50

Interacttve Interactive

104 104

Prescribed

5.69 5.81 -...--..-.--. --

miliis~onds/o~ration.

is taken into

account. Considering the discretization and time step size, the results are generally accurate and stable. ADINA results for Problem 1 (Figs. 3-5) and Problem 2 (Figs. 7 and 8) are in sufficiently good agreement to verify adequacy of the modified code for this type of application.

Acknowledgements-The

writer thanks Mr. P. N. Roth of DTNSRDC for his automatic data processing support, without whrch an essential singularity in the writer’s code would not have been localized. He also is most grateful to 1Mr. B. D. Johnson (DTNSRDC) who created the figures and supervised the manuscript preparation. Finally, thanks are due to Dr. H. Huang of NRL who pattently answered the writer’s questions,

Fluid-structure

interaction and ADINA

REFEREES 1. S.

Zilliacus, T. Toridis and T. Giacofci, Analysis of wave excited submerged structures. Proc. ABINA Conf., Mass. Instit. of Technology, 42.5-445 (Aug. 1979). 2. H. Huang, Transient interaction of plane acoustic waves with a spherical elastic shell. J. Acoust. Sot. Am. 45, 66-670 (1969). 3. T. L. Geers, Excitation of an elastic cylind~cal shell by a transient acoustic wave. J. Appl. Mech., ASME. 36, 459-469 ( 1969). 4. T. L. Geers, Residual potential and approximate meth-

ods for three-dimensional fluid-structure interaction problems. J. Acoust. Sot. Am. 49, (Part 2), 1505-1510 (1971). 5. R. D. Mindlin and H. H. Bleich, Response of an elastic cylindrical shell to a transverse, step shock wave. J. A&. Mech., ASME. 20, 189-195 (1953). 6. H. Huang, A qualitative appraisal of the doubly asymp-

totic approximation for transient analysis of submerged

7.

8.

9.

10.

173

structures excited by weak shock waves. Memorandum Rep. 3135, Naval Research Laboratory, Washmgton, D.C. (Sept. 1975). T. L. Geers and J. A. DeRuntz, Transtent response submerged structures. Rep. analysis for LMSC-D313297, Structural Mechanics Laboratory, Lockheed Palo Alto Laboratory (Dec. 1972). K. J. Bathe, ADINA-a finite element program for automatic dynamic incremental nonlinear analysts. Rep. 82448-1, Acoustics and Vibration Laboratory, Mechanical Engineering Department, Massachusetts Institute of Technology, Sep. 1975 (rev. Dec. 1978). H. Huang, Transient response of two fluid-coupled cylindrical elastic shells to an incident pressure pulse. NRL Memorandum Rep. 3821, Naval Res. Lab. (Aug. 1978). H. Huang, An exact analysis of the transient interaction of acoustic plane waves with a cylindrical elasttc shell. J. Appl. Mech., ASME. 37, 1091--1099 (1970).