Two efficient schemes for dynamic analyses with finite elements

141

Nuclear Engineering and Design 66 ( 1981) 141- 146 North-Holland Publishing Company

TWO EFFICIENT S C H E M E S FOR DYNAMIC A N A L Y S E S WITH FINITE E L E M E N T S Ronald USHIJIMA and Eduardo KAUSEL Dept. of Civil Engineering, M I T . CamhridRe. MA. USA

Received 9 February 1981

This paper proposes two efficient schemes that may be applied to finite element analyses in the frequency domain~ such as those currently used to analyze nuclear power plants under seismic loads. The first scheme consists of a reduction in the number of degrees of freedom necessary to model an axisymmetric system, using what is referred to as the quasi-cylindric approximation (not a plane strain approximation!). The second scheme is the use of Hermite interpolation to decrease the computational effort associated with the determination of the transfer functions. The application of the technique is reviewed for the case of a Waas Lysmer transmitting bounda~ (a similar boundary is used in the computer program FLUSH), and the relative efficiency is discussed.

1. Introduction The special design requirements imposed by certain classes of dynamically loaded structures, such as nuclear power plants, radar tracking stations, offshore towers and very tall buildings, has stimulated significant interest in the development of analytical and numerical tools for the solution of foundation interaction problems. Although significant progress has been made, much remains to be investigated, particularly in the area of seismically loaded embedded structures and structural groups. The need for more research in these areas was also attested to at the recent workshop, "Research Needs and Priorities for Geotechnical Earthquake Engineering Applications," held in Austin, Texas, in June 1977. Recent advances in foundation mechanics problems have taken place along three parallel avenues: (1) development of analytical solutions for relatively simple geometries, such as cylindric and hemispheric (or even ellipsoidal) foundations on elastic halfspaces [2,8,17], as well as arbitrarily shaped surface plates on layered media [4,16]; (2) implementation of computer codes with advanced methodologies [7,10,14]; and (3) development of semianalytical or empirical procedures [ 11,12], etc. The analytical solutions serve several important functions, beyond providing an answer to the particular problem (usually of a highly idealized nature) under consideration. On the one hand, comparison of the solutions available provides insight into the behavior of dynamically loaded foundations; on the other hand,

these solutions consitute valuable yardsticks against which approximate and numerical procedures, as implemented in manu computer codes, can (and should) be tested. However, there is little hope that 'exact' solutions may become available for every imaginable or practicable situation. Hence, it is here, in the domain of complicated foundation geometries and soil properties, that numerical procedures such as the finite element method are most attractive. It is then highly desirable to develop computationally efficient schemes without sacrifice of accuracy. This paper seeks to explore the feasibilities of two of these alternatives. 2. Quasi-cylindric approximation In the axisymmetric formulation, stresses and displacements are expressed in terms of Fourier expansions in the aximuthal direction. For vertical excitation and torsion, only the zeroth order component is needed to define the problem, necessitating two degrees of freedom per node for the vertical case (radial and vertical displacement), and only one degree of freedom per node for the torsional mode (the tangential displacement). Horizontal translation (swaying) and rotation (rocking), on the other hand, involve only the first order component in the expansion, with three degrees of freedom per node: u r = u ri

where

0 0 2 9 - 5 4 9 3 / 8 1 / 0 0 0 0 - 0 0 0 0 / $ 0 2 . 5 0 © 1981 N o r t h - H o l l a n d

cos 0,

u 0 = - u ~ sin 0,

Ur, U o, Uz a r e

u .~= u ! c o.s O ,

(1)

the radial, tangential and vertical

142

R. Ushijima, E. Kausel / l)vnamic ana(lses with l)mtc clement.~

displacements; u,,I U~ u~1 are the amplitudes which are functions of r and z only (the index '1' referring to the first order component); and r, z, 0 are the cylindric coordinates (radial, vertical and azimuthal, respectively). For the points of the plate, as well as the soil in immediate contact with the plate, the radial and tangential amplitudes are identical, i.e., ul~ = u~ for points in contact with the plate. This fact follows from the rigidity condition in the plate. Dissimilar values would result in an ovoidal deformation of a previously circular ring, which is not possible within the rigid plate. Points along the vertical axis of symmet~' must also satisfy this condition, since the~, cannot be associated with two different displacements without rupturing the axis. If the soil is supported by a rigid bedrock, then this interface too satisfies the conditions of equalityof tangential and radial components of displacement, Notice also that most classical analyses of axisymmetric superstructures (for instance, nuclear containment shells) use this condition implicitly when the structure is modeled as a cantilever shear or bending beam. The in-plane (ovoidal) distortion is thus neglected. On the other hand, experience gained with a large number of analyses with a sophisticated axisymmetric finite element computer code incorporating transmitting boundaries [5] has shown that for points in the soil underneath and away from the foundation, the condition stated above is approximately satisfied, with

discrepancies of less than 10% (even at points relatively removed from the foundation). This suggests equating a priori the tangential component to the radial component, and performing a reduction in the mnnbcr of degrees of freedom from three to two: this pwcedurc may be labeled "quasicylindric approximation.' While the resulting finite element model has only two degrees of freedom per node, it is no[ a plane-strain model, even though it compares favorably (in terms of efficiency) wit codes based on a plane-strain formulation. As will be shown, this scheme results in onlv a slight loss of accuracy, while realizing dramatic savings in computer storage and execution time. A decrease in the number of degrees of freedom by a factor of 2/'3 reduces the size of the stiffness matrices by (2/3) ~ ,0.44 (less than 1/2!), and the number of operations in the Gaussian reduction by (2/3) ~ = 0.3. To explore the accuracy of the quasi-cylindric approximation, a number of actual cylindric foundations embedded in an elastic stratum were analyzed with the approximate and rigorous models (fig. 1). Table 1 shows a comparison of the static stiffnesses as computed by programs based on the two models. As expected, the approximate model yields slightly larger stiffnesscs than the true cylindric model, because of the displacement constraints which have been added: however, the results are essentially identical. Figs. 2 and 3, on the other hand, show the dynamic

For all cases: R=

1

G=

I

~

SWAYING

ROCK I NG

'V'= 1/3

/ ////I////

///////////

/ ///////

E=O

.a-

//ll/J~ i

E=~2

ii/

Fig. 1. Models used for analyses.

/

////////

E=I

////11/7:

E=

R. Ushifima, E. Kausel / Dynamic analyses with finite elements Table 1 (Poisson's ratio = 0.33)

l.O

Stiffness component

Stratum depth H/R

Embedment E/R

Quasicylindric (approx.)

Cylindric ('exact')

Swaying

2

0 1/2 1 3/2

6.14 11.46 16.82 25.57

6.09 11.42 16.79 25.55

0 1/2 I 3/2

5,70 10.08 13.75 17.73

5.65 10.04 13.72 17.69

0 1/2 1 3/2

4.49 9.44 18.45 35.35

4.35 9.22 18.16 34.84

0 1 3/2

4.34 8.74 16.23 28.38

4.22 8.57 16.04 28.17

0 1/2 1 3/2

-0.22 1.80 5.85 12.59

-0.31 1.71 5.76 12.51

0

-0.29

-0.37

oa ts.i 0.6

-

3

Rocking

2

3

l/2

Coupling

2

3

1/2

1.43

1.36

1 3,/2

4.71 9.47

4.67 9,47

stiffness coefficients (i.e., the dynamic stiffnesses divided by the static stiffness) for some typical cases. As before, the agreement between the fully cylindric and quasi-cylindric models is excellent. Frequently the computer codes used for soilstructure interaction analyses are based on a twodimensional formulation [10]. However, plane strain idealizations do not, in general, provide a fully satisfactory representation of an axismmetric structure. The quasi-cylindric approximation, on the other hand, has the same number of degrees of freedom per node as a plane strain model, while providing a vastly better representation of the soil and structure. Thus, the scheme proposed here is an attractive alternative to the less accurate plane-strain formulations currently in use, and to the more costly fully cylindric models. 3. Hermite interpolation

It is well known that finite element formulations in the frequency domain are obtained by determining the

143

~

QUASI-CYLINDRIC

k/ J

cY,,,oRic

04

ec ~

0.2

0.0

o~o

G,

d.2

;.3

G.4

&

&

;.7

().6

().7

FREQUENCY (NZ)

i.o

o.B

aJ o_- o6

3 ~ o.4 ',

t~ 0 2

o.o

0.0

0.1

0.2

03

0.4

FREQUENCY

0.5

( HZ )

Fig. 2. Comparison of the cylindric and quasi-cylindric models H / R = 2, E / R = 0, Poisson = 1/3, damping = 0.05. Normalized swaying stiffness. smooth transfer functions at a few hundred frequencies, which are ten interpolated to match the frequency spacings in the Fourier spectrum of the excitation. Since the computation of every point requires the solution of a system of linear equations, it is desirable to minimize the number of sample points without sacrificing the resolution of the transfer functions, particularl in the neighborhood of resonant peaks. Such a goal can be achieved with Hermite interpolation; that is, finding not only the (complex) amplitudes of the transfer functions, but the slopes (derivatives with respect to frequency) as well. This permits essentially halving the sampling rate as well as the number of operations. The procedure, which was brought to the authors' attention by Waas [15], is as follows: Let M, C, K be the global mass, damping and stiffness matrices, and U, F the displacement and load-

R. Ushijima, E. Kausel/ Dynamic ana(rses wtthjinite e/emc,lI~

144

parenthesis 08

Kd =-w2M + iwC + K

O6

is the dynamic stiffness matrix. The vector U is found in the classical scheme by Gaussian elimination. Formally,

QUASI-CYLINDRIC

(3)

U= Kd IF.

(4)

04

Taking derivatives of the above equation with respect to oa (and since M, C, K, F are independent of o0 one obtains

02"

O0

, OI

0.0

(~ 2

014

v 03

FREQUENCY

I 0,5

1 06

( - 2 ~ o M + i t ) U + K~U'= o,

, 07

(5)

where

(HZ)

U' = dU/dto.

(6)

I0

It follows then that

K o U ' = ( 2 t o M - iC)U,

08 }z uJ 06

QUASI-CYLINDRIC

L) o 0.4 ,=I

_~

L> 02

O0

r'

OI

0.0

02

0.3

/ CYLINDRIC r

r

04

0.5

FREQUENCY

I

06

07

(HZ)

Fig. 3. Comparison of the cylindric and quasi-cylindric models H/R=2, E/R=O, Poisson= 1/3, damping=0.05. Normalized rocking stiffness.

ing vectors. For a solution in the frequency domain, the dynamic equilibrium equations are

(-to:M + itoC + K ) U = F,

(2)

with ~0= frequency of excitation. The expression in

//•v'2 "4 FI X1

x

x 2

Fig. 4. Hermite interpolation.

,

(7)

which has the same form as the original equation, with only a modified loading vector. It is well known that most of the operations in the Gaussian elimination are related to the triangularization of K d, not to the back substitution. Therefore, determination of U' requires only minor additional effort after U has been determined. In the above equation, it has been assumed that K, the stiffness matrix, is independent of frequency. When working with macroelements or transmitting boundaries, on the other hand, this is not the case, and additional terms are needed. (Transmitting boundaries are special elements placed at ihe edges of the finite element grid to prevent wave reflections or echoes.) To examine this point, consider first the case of a Waas-Lysmer transmitting boundary in the two-dimensional space (plane strain). Following Waas [14], the (symmetric) b o u n d a ~ stiffness matrix R is given by an expression of the form

R=iAXKX

~+D,

(8)

with i = ~ 1 ;A,D = matrices assembled with the material and geometric properties of the layered soil (A is symmetric, while D is not); X = matrix whose columns are the natural modes of wave propagation in the soil (obtained from a quadratic eigenvalue problem); and K is a diagonal matrix with the wavenumbers associated with the propagation modes. Details of the formulation are not important in this context; the interested reader may find them in the work of Waas [14]. The quadratic eigenvalue problem from which the modes X and wavenumbers K are obtained is of the form (see Waas [14]):

AXK'~-+ i ( D r - D ) X K + ( G

~:M)X=(),

(9)

R. Ushijima, E. Kausel / Dynamic analyses with finite elements with G , M being again symmetric matrices assembled with material and geometric properties of the soil; (the matrix M above is not the same as the mass matrix of the finite element systems in eq. (2)); 0 is the null matrix, and co is the frequency of excitation. Multiplying by i X t from the right, and factoring out the term X K X - ~, one obtains (iAXKX-'

+ D - D T ) X K X - ' + i(G - ~2M) = 0,(10)

but i A X K X -~ + D - DT = R - DT = R T - DT = R T D T = ( R - D) T (since the boundary matrix can be shown to be symmetric). Also, i X K X - I = A t( R -- D). Hence, (R-D)TA-'(R-D)-(G-o~2M)=O.

(ll)

Taking derivatives with respect to frequency, one obtains R'A

I(R-D)4-(R-D)TA-1R'=-2wM,

(12)

or defining q = XKX

'

(13)

eq. (12) may be written as R'Q 4- QTR' = 2i~M.

(14)

This matrix equation is linear in R' for a given frequency. Since the product X = X K X - ~ is available from the generation of the boundary matrix, in principle, the calculation of the elements of the symmetric matrix R' can be performed solving a system of ½n( n 41) equations, with n being the number of rows (or columns) of R. Such a direct scheme is not attractive, because of the considerable expense involved in solving a system of equations so large. However, some computational advantages can be achieved taking into consideration the structure of the matrix Q. Multiplying from the left by K -~/2xT, and from the right by XK I/2, o n e obtains K

I/2XTR,XKI/2 + K I / 2 X T R , X K =2i~0K

I/2XTMXK

I/2

I/2.

(15)

Defining the elements of the symmetric matrices x T R ' x and X v M X by ri/ and m,i respectively, one obtains from the above equation [ k ~t/2

+:

l~,l

m,,

]r,/=2io~(kik:),/2

Dlil

Hence R'=iooX

r(~ij}X-'

(18)

While the above expression is more efficient than the direct solution of eq. (14), it still demands a significant number of calculations. Specifically, forming the prod u c t X T M X = m i : requires somewhat more than ½n 3 operations ( M is tridiagonal); the modified elements mij : m o / l ( k i 4- kj) are formed with about ½rt 2 operations, while the product x - T ( ~ , j ) X I implies some n 3 operations (the inverse, X t, or at least the triangular form of X - t , is available from the generation of the boundary matrix). Thus, computation of the derivative of the boundary matrix with respect to frequency requires on the order of 2n 3 operations; this is somewhat less than four times the effort required to form the boundary matrix itself. Nevertheless, this expense is relatively small compared to the operations involved in solving the system of equations for the finite element grid plus the boundary. For example, if the finite element grid consisted of the boundary endc-I columns of finite elements, the total number of degrees of freedom would be nc, with bandwidth b = n. The solution of this system would require some N = n 3 ( 3 c - 2)/6 + n2c/2 operations. If the finite element mesh consisted of some five columns (a minimum typical number), then c = 6 and N,~ 3n 3, to which the ½n 3 operations to form the boundary matrix must be added. Thus, at least for the plane-strain case, forming the derivative of the boundary matrix and using Hermite interpolation is computationally attractive. In the case of cylindric coordinates, on the other hand, the boundary matrix has a more complex structure. While it is also assembled with the same propagation modes as the plane-strain case (see [5, 6]), the individual components are weighted with cylindric functions (Hankel functions), which have the wavenumbers as argument. Even though it is possible to determine the derivatives with respect to frequency of each mode, wavenumber and cylindric function, the implied procedure is too cumbersome to make the procedure attractive. It is possible that a simpler scheme to determine R' for the cylindric case may exist, but it has so far eluded the authors. Nevertheless, Hermite interpolation may still be used for cylindric formulations if consistent transmitting boundaries are not employed.

(16)

Acknowledgement

and from here, ri: = i o~ ½(/q + k,) -

145

(17)

The research reported in this paper was made possible by National Science Foundation Grant No. ENG 79-08080.

R. Ushijima, K Kausel / l)vnamic analyses with finite element~

146

Appendix A. Hermite interpolation formula Consider a function y = f ( x ) a n d its derivative y ' = which are b o t h sampled (or are known) at intervals A x. It is desired to determine interpolated values within the interval Ax. It can be shown that the appropriate interpolation formula is

df/dx,

y = (1 -- t 2 ) ( l + 2 t ) y 1 + (3 -- 2t)t2y2

+ axt(t-

l)[(t-

l)y{ + ty~],

in which t = (x - x,)/(x~

- x,) = (x - x, )/a~

and Yl, Yz, Y{, Y~ are the values of the function a n d its derivative at the beginning and end of the interval, as s h o w n in fig. A. 1. In particular, for the m i d p o i n t t = 1/2, and

Y,/2 =½(Y, +Y2) +½ A x ( y [ - - y ~ ) .

References [1] A.H. Ang and N.M. Newmark, Development of a transmitting boundary for numerical wave motion calculations, Report to Defense Atomic Support Agency, Contract DASA-01-0040, Washington, D.C. ( 1971). [2] R.J. Apsel and J.E. Luco, Torsional response of rigid embedded foundation, J. Engrg. Mech. Div. ASCE, (Dec. 1976) 957-970. [3] P. Bettess, Infinite elements, Internat. J. Numer. Meth. Engrg. II (1977) 53-64. [4] G.C. Gazetas, Dynamic stiffness functions of strip and

rectangular footings on layered ~oil, S.M. Thesis, Dept. of Civil Engrg. Mass. Inst. Technol., Cambridge, MA. 1975. [51 E. Kausel, Forced vibrations of circular foundalions on layered media. Research Report R74-1 I. Dept. of Civil Engrg., Mass. Inst. of Technol., Cambridge, MA, 1974. {6] E. Kausel and J.M. Roessct, Scmianalytic hypcrelement for layered strata, J. Engrg. Mech. Div., AS(T'. (August 1977) 569-588. [7] R.L. Kuhlemeyer, Vertical vibrations of footings embedded in layered media, Ph.D. Thesis, Univ. of California, Berkeley,CA, 1969. [8] J.E. Luco, Torsion of a rigid cylinder embedded in an elastic halfspace, J. Appl. Mech. 43 (E3) (Sept. 1976) 419 423. J. Lysmer and R.L. Kuhlemeyer. Finite dynamic model for infinite media. J. Engrg. Mech. Div.. ASCE (Aug. 1969). [10] J. Lysmer et al., FLUSH A computer program for approximate 3-D analysis of soil-.,,tructure interaction problems, Report No. EERC 75-30, Earthquake Engineering Research Center, University of California. Berkeley, CA, Nov. 1975. [I 1] J.W. Meek and A.S. Veletsos, Simple models for foundations in lateral and rocking motion. 5th World Conference on Earthquake Engineering, Rome, 1973, paper No. 331. [12] M. Novak, Vibrations of embedded footings and structures. ASCE National Structural Engineering Meeting. April 1973. San Francisco, CA, Meeting preprint 2029. [13] W.D. Smith, A non-reflecting boundaD for wave propagation problems, J. Comput. Phys. 15 (1974) 492. 5(12. [14] G. Waas, Linear two-dimensional analysis of soil dynamics problems in semi-infinite layer media, Ph.D. Thesis, Univ. of California, Berkeley, CA, 1972. [15] G. Waas, Personal communication. [16] H.L. Wong and J.E. Luco, Dynamic response of rigid foundations of arbitrary shape, Earthq. Engrg. Struct. Dyn.(1976) 579 587. [17] H.L. wong and M.D. Trifunac, Interaction of a shear wall with the soil for incident planet SH waves: elliptical rigid foundations. Bull. Seism. Soc. Amer. 64 (1974) 1825-- 1842.