- Email: [email protected]
An efficient numerical algorithm for transient analysis of exterior scalar wave propagation in a homogeneous layer
. IB
JQzJ
Computer methods in applied mechanics and engineering
mm
mm
Comput.
EISEVIER
Methods
Appl. Mech. Engrg.
167 ( 1998) 261-273
An efficient numerical algorithm for transient analysis of exterior scalar wave propagation in a homogeneous layer Mu-thy N. Guddati”“, “Tercu Imtitutr
,jbr Cornputationctl hDepartmerzt
cud Applied
of Ci\?l Engineering,
Received
John L. Tassoulas b3*2
Mathematics,
The UniversiQ
The University
18 November
qf
cf Tern... ut Austin,
Texas ut Austin.
1996; revised
I8 August
Austin.
Austin.
TX 7X7/2- 1076, USA
TX 787/2- 1076. USA
1997
Abstract In this paper, a computationally efficient alternative for the analysis of transient wave propagation in an infinite layer is presented. The new method is based on a simple cell centered characteristic finite difference technique. This approach can be used in place of the existing methods based on closed form Green’s functions. as it would save considerable computational time and storage. Numerical examples are presented to illustrate the performance of the method. Some guidelines are suggested to further increase the efficiency of the proposed method. 0 1998 Elsevier Science S.A. All rights reserved.
1. Introduction Solution of a scalar wave equation in an infinite layer is generally encountered in fluid-structure interaction and soil-structure interaction problems. A typical example is the analysis of dam-reservoir interaction. The reservoir can be idealized as an infinite fluid layer fixed at the base. Another problem of similar type is the simulation of an elastic layer under out-of-plane motion. In this case the elastic vector wave equation reduces to a scalar wave equation and the procedure is essentially the same as the dam-reservoir interaction analysis. The elastic layer problem is considered as the model problem for the purposes of our discussion. Various techniques are used for the analysis of these types of problems [5,12]. The boundary element method is used for the analysis of the dam-water interaction problem by Tsai and Lee [lo]. They model the infinite surface and the interaction surface with boundary elements. The disadvantage of this approach is that the problem size is considerably large and requires the evaluation of singular integrals. An alternative to the boundary element method is the use of a Dirichlet-to-Neumann (DtN) map to replace the infinite layer [5]. Local DtN maps are sometimes used, but they are known to be crude and are good only for specific problems [3,8]. A more accurate way of modeling the effect of the infinite layer would be the use of a DtN map which is global in both time and space. Separation of variables is used for the derivation of the Green’s function for this system. The Green’s function consists of Bessel functions which are expensive to compute. The resulting Green’s function is used in conjunction with convolution to do transient wave propagation analysis of these types of systems. This approach demands high computational cost and the storage of the Green’s function. Semidiscretization techniques, which are equivalent to truncation of analytic Green’s functions. are often used to obtain a discrete Green’s function. Instead of Bessel functions, consistent
* Correspondmg author ’ Graduate student. ’ Professor. 00457825/98/$19.00 0 1998 Elsevier Science S.A. All rights reserved PII: SOO4S-7825(98)00123-6
262
M.N. Guddati. J.L. Tassoulas
/ Comput. Methods Appl. Mech. Engrg.
167 (1998) 261-273
transmitting boundaries can be used to calculate the dynamic stiffness matrix (see [9]) and Fourier transformation can be used to arrive at the impulse response function (Green’s function in time domain). These methods are extremely expensive as they involve either Fourier transformation or the computation of Bessel functions and their integrals. Tsai et al. [ 1 I] use modal decomposition to get the discrete Green’s function and the technique they propose is more efficient than the other existing methods. They, however, do not utilize the full decoupling which will reduce the computation time by an order of magnitude. Furthermore, in their approach, it is necessary to compute the integrals of Bessel functions. Another disadvantage of this method is that there is considerable loss of efficiency when the time step size is not constant. In this paper, an efficient alternative to the above methods is proposed. This method, based again on separation of variables and dimensional reduction technique, is considerably simpler to implement and extremely efficient. Semidiscrete equations in space-time are solved based on a characteristic finite difference technique and a symmetric positive definite DtN map is derived. There is no change in efficiency when the time steps are of variable size. The presentation of the rest of the paper is as follows: Section 2 states the model problem and the governing differential equations. In Section 3, dimensional reduction based on separation of variables technique is discussed. The derivation of the DtN map is carried out in Section 4 by means of the solution of the resulting system of hyperbolic equations. In Section 5, the computational cost of the current method is compared with the existing methods. Section 6 contains numerical examples which illustrate the performance of the method. We conclude the paper by discussing possible other applications and limitations of the proposed method in Section 7.
2. Problem
statement
Although the problem that is being solved here has broader motivation, we give a particular example of shear wave propagation in an infinite elastic layer. We have an elastic soil layer resting on a rigid rock and the loading is such that there is only out-of-plane deformation. The aim is to simulate the interaction between the layer and the structure embedded in or resting on it. The deformation is considered invariant with the out-of-plane position, i.e. the stress state in the soil is antiplane shear. The layer is assumed to be homogeneous. The general governing differential equation in the layer is the vector wave equation [I]: v* (C : Vu) - pu,, =f ,
(1)
where u is the three-dimensional displacement field, C is the elasticity tensor andf is the body force. With the restrictions described in the previous paragraph, the vector wave equation reduces to a scalar one: GAu-pu,,=f,
(2)
where u is the out-of-plane displacement field, G is the shear modulus andf is the corresponding source term resulting from the body force. The following procedure is used to analyze the structure (Fig. 1): a finite region of the soil bounded by straight vertical edges is modeled by finite elements. The rest of the infinite layer on either side of this region is simulated using a Dirichlet-to-Neumann (DtN) map. The aim of this work is to get an appropriate time dependent DtN map for general transient analysis of this system. The problem can be stated as follows: find the relation DtN : u,. -+s, in which u,. : f X I + !‘H is the displacement on the boundary r at time t E I= [0, T,,,] and s : r X I + % is the traction on the boundary. The applicable differential equation and the imposed boundary and initial conditions on the layer are given by Ghu-pu,,=f
forxE0,
tEZ,
where R = !h’ X (-H, 0) is the domain of the soil layer. Homogeneous Dirichlet condition bottom of the layer (z = -H) and Neumann condition on the top of the layer (z = 0): u(x, z, t) = 0
and
at z = -H
(3) is applied on the
(4)
M.N.
Guddati,
J.L. Tassoulas
I Comput.
Methods
Appl.
Mech.
Engrg.
167 (1998)
261-27.1
263
Soil Layer Rock
(a) Physical Problem
DtN map
DtN map
Bounded Soil (FJN Fixed Base (b) Mathematical Model Fig. I. General
u:(x, z, t) = so(x, t) The soil is assumed
media.
(5)
to start from rest: at t = 0 .
reduction
(6)
of the problem
The separation of variables technique field by the following approximation: 2 = rt
for unbounded
at z = 0.
u(x, 2, t) = u,(x, z, t) = 0
3. Dimensional
method of analysis
4,(Z)U,(X, f) = QT(z)U(4
is used to reduce the problem dimension.
t) .
We discretize the displacement
(7)
In the above form of approximation, a(z) can be viewed as the basis function vector in the ; direction and U(X, t) is the nodal degrees of freedom vector. We apply the standard Galerkin method using this approximation. The test function is of the same form as the solution: uh = ,; The Galerkin
&(Z)U,(X, r) = @T(z)V(X, t) . approximation
of the variational
u”(G AU” - p$, -f)
Integrating
dz = 0
by parts with respect to z, we get
(8) form for Eq. (3) is (9)
264
M.N. Gudduti, J.L. Tussoulas
ad au”
cl
I
-G Substituting
a2uh
--dz+G
a2uh
v”--dz
u”Tdz-p
az az
-H
I Comput. Methods Appl. Mech. Engrg. 167 (1998) 261-27.3
-
uhf dz=O.
at2
the approximation
(10)
for uh and u”, we get 0
-GVT
@$j;
dz U + GVT
I
c&PDTdz U.,,
-H
0
I
-pvT
QWT
dz U,, -VT
-H
which reduces to the following
I
” wdz+GVT@\;=,s,=O
vVVEMN
If
(11)
form:
-GU+AU,,-L*U,,-F=Q, c2
(12)
where c =J/G/p is the wave speed. In the above equation, G and A are constant N X N matrices and U and F are vector functions of x and t. G, A and F are given by
(13)
(14) and F(x, t) =
(15)
Note that the degree of freedom corresponding to the base is eliminated and the matrices G and A turn out to be symmetric and positive definite. Our modified problem is to find a map between the displacement vector U, and the traction vector consistent to the approximation made in displacement, i.e. if S is the traction vector, we need to satisfy the work equivalence:
(16) where sI, is the approximation relation
of the traction
s. As s = -Gu,
we have s,, = - GTGIJ,
and we will have a
0
I
VTS = -VT
PN
@GaTdzU,
This implies that the traction S = -AU,
b’VEaN.
vector is nothing
(17) but
.
(18)
We now present the discrete version of our original problem: Find the map DtN : U, -+ S where U, : I + sHN is the displacement vector at x = 0 and S : I + SN is the corresponding traction vector. We satisfy the following differential equations and initial conditions: AUxx-LAU,,-GU-F=O
forxE%+,
tEZ
(19)
C2
U=U,=Q
forxE%+,
t=O.
(20)
M.N. Gudduti, J.L. Tussoulas I Comput. Methods Appl. Mec,h. Engrg. 167 (1998) 261-273
26.5
Derivation of the DtN map In this section, the hyperbolic nature of the system of equations derived from the semidiscretization procedure used to arrive at an efficient DtN map. We first decouple the system of equations by spectral decomposition of with respect to A. The system can be written as QQ’~.,-~QQ’c!,-QAQ’u-F=U
(21)
where Q is the modal matrix and A is the diagonal Premultiplying this equation with Q-’ gives
matrix containing
the eigenvalues
o,,--L q,,- /@-g =0 c2
of G with respect to A.
(22)
in which i?=Q’U and @=Q-‘F. As an aside, we note here that the above system of equations are uncoupled second-order hyperbolic equations. This fact can be used effectively to increase both time and storage efficiencies when using discrete Green’s functions. If we use the approach based on Bessel functions as in [ 111, we can store the values of cinstead of U and increase the computational efficiency by an order of magnitude. We, however, do not utilize this modification as it is not necessary to resort to the use of Bessel functions in the proposed approach. We note that the uncoupled system of equations (22) have the same characteristic speed. This gives us the advantage of transforming all the equations into standard canonical form using the same transformation. For all the equations, we substitute .Y+ ct t=2 The canonical
and
x - ct 7 = 2 .
(23)
form will then look like
ir,,-n&F=0 Our objective is to compute rule, we have
(24) l?, at (x, t) = (0, 7’,,) or, equivalently,
at (5,~) = ([,,, v,,) (see Fig. 2). By the chain
(25)
‘t
_.-C
* q=-cT,/2 11= -CT n-, I2
‘-“1;--i_.g; Fig. 2. Finite difference
grid for updating
0.
266
M.N. Guddati, J.L. Tassoulas / Comput. Methods Appl. Mech. Engrg. 167 (1998) 261-273
Integrating
Eq. (24) along the characteristic
On account of Eq. (25)
5 = &, (with respect to v), we find
we obtain
It is important to note that 0 and therefore fit can be computed along the &characteristic at v= 0 provided that the source term (8) and the initial conditions on fi (at t = 0, x 2 0) are known. This follows from the fact that the domain of dependence [2,4] of any point on the v = 0 characteristic is contained in the line {X2 0, t = 0). The Cauchy data (6, 6,) on this line, together with the source term (P) uniquely determine the value of 6 on the characteristic 7 = 0 (see [2,4]). This in turn determines fit at r] = 0. Thus, the computation to arrive at o& &, 0) is carried out, in many cases explicitly, prior to the initiation of the transient analysis. For a medium starting from rest (0 = tit = 0) with vanishing source term (P = 0), it is clear that 0 = c6 = 0 along 7 = 0. In this case, we can replace Eq. (27) by
We now divide the space-time domain into cells which have the edges parallel to the characteristics (Fig. 2). fi is assumed to be bilinear in each cell. We then require the values of 6 at the grid points only. We use the following notation for the purposes of presentation: 0: is the value of e at (5, q) = ([,, Q) = (cT,/2, -cTk /2). The interval between T,, ~, and T, is considered as the nth interval. We now describe the procedure for obtaining the values of IJ?on the v characteristic corresponding to the current time T,,, i.e. 5 = [,,. The set of active cells in this step are shown in Fig. 2. The ordering of the analysis is also shownk in the same, figAu;e: we stA$rom cell (n, 1) and proceed till (n, n - 1). For each cell (n, k), we need to obtain U,,, given U,,_, , U,j_, and U,, . We use cell centered finite difference technique for this purpose. We satisfy the differential equation at the center of the cell:
(29) As we have A&,, = c AT, 12 and Aqk = -c AT, 12 (we are sweeping following form:
in negative
77 direction),
we arrive at the
o:(i)=a;@:_,(i) +ii;-‘(i), - o;::(i) + p;Fyq, where 4
-- A(i) c’ AT,, AT, 4 cu,=
4
+ A(i) ’
c2bTn AT,
4
(31)
and
P;=
41 +4) c2 AT,, AT,
’
(32)
4
In the above equations, 6(i) and P(i) stand for the ith element of 8 and P, respectively. A(i) is the ith diagonal element of the matrix A. Note here that the coefficient cy, is always less than unity for any mode. It is expected that this fact is crucial for the stability of the method. A proof of stability is not available at this time. Numerical experiments indicate that the method is stable, as long as 4 > 1. A general proof of stability is currently under investigation.
M.N. Gudduti, .I.L. TassoulaJ I Comput. Methods Appl. Mech. Eqrg.
167 (1998) 2hlk-77.1
267
Thus, with l? known on the previous characteristic and at v= 0, we can get 6 on the new characteristic using Eq. (30). Once 6 is known on this characteristic, we use standard numerical integration to compute the integral in Eq. (27). The most appealing integration procedure would be the one which is consistent with the approximation of 0 over the mesh. This suggests the piecewise linear interpolation on the characteristic and we thus use trapezoidal rule for numerical integration. It should be noted that higher order difference schemes or finite element schemes could also be used for the updating, and consistent integration schemes could be used for the integration. Using the trapezoidal rule, Eq. (27) can be replaced by the following summation
(33) which can be rewritten
as
(34) We premultiply
the above equation
by -GQ
to get
(35) since A = QQ’ and G = QAQ’. It can be seen that the left-hand side of this equation is nothing but the traction vector that we are interested in. We thus have the following relation between the traction and displacement at x = 0: S=KU+CZ+H
at
(36)
3
where c AT K= 4G:
(37)
C=+A,
(38)
and c(AT,
+
AT,, ,)
(AO;;+$;)+yF,s
4
CAT
no
-A-$. I
au Cl (39)
Notice that both the coefficient matrices K and M are symmetric and positive definite. Thus, for conventional time-stepping schemes, we get a symmetric positive definite coefficient matrix from the boundary which has definite advantage for linear solution purposes. Another advantage, when the discretization is local, is that the DtN map in Eq. (36) preserves the sparsity of the interior coefficient matrices. This is due to the fact that the matrices G and A are sparse when the interpolation function @J is sparse (see Eqs. (13) and (14)). Before we conclude this section, we mention an alternative way of obtaining the DtN map. In contrast to Eq. (25), we can use fi,, to compute o,r: (40) l?, at (x, t) = (0, T,,) can be obtained
locally by a backward
difference
approximation:
M.N. Guddati, J.L. Tassoulas I Comput. Methods Appl. Mech. Engrg. 167 (1998) 261-273
268
(41)
Higher-order the resulting
difference approximations could be used as well. Skipping relationship between the traction and displacement:
all the intermediate
details, we present
(42)
where (43)
and Hz-
&Qc::-,
.
(45)
Note that this form does not involve the cost of computing the integral. However, it has several disadvantages. The most important one is the possible lack of stability: the matrix c is negative definite and when augmented with the interior coefficient matrices, it can lead to an unstable analysis procedure because of the loss in positive definiteness. Furthermore, it can be shown that the computation of traction from Eqs. (36)-(39) i.e. by means of integration along the characteristic, is second-order accurate. On the other hand, first-order difference approximations, like the one in Eq. (4) limit the accuracy of the overall scheme. Higher-order difference approximations can circumvent this difficulty, but require greater smoothness of the field variable. The boundary condition in Eqs. (36)-(39) does not suffer from these drawbacks and we choose to use it in spite of the need to evaluate the integrals. To summarize this section, we have derived a DtN map which relates the displacement vector on the boundary with the corresponding traction vector. The relationship is global in space, but local whenever the interpolation functions used for the semidiscretization are local. However, it is pseudo-local in time as it uses only the displacement vector and the time derivative in a direct sense. The past history is stored in a different form with the displacement vector stored along the characteristics. It should be noted that, unlike methods with higher order derivatives, this method is expected to be stable as it relies only on the first time derivative and integration. We now go ahead with the presentation of the computational cost after some comments on implementation issues.
5. Computational cost From the previous section it is clear that the values of t?? on the current characteristic, and the matrices G, A and Q must be stored. For a general setting (i.e. assuming the matrices G and A are dense), the initial required storage would be 3N2 where N is the number of degrees of freedom on the boundary. When piecewise polynomials are used for the purposes of semidiscretization, the matrices G and A are banded and the storage would be comparatively less. The storage required for the history (6 on the characteristic) at kth time step is given by k X N. Also, it is noted that the full vector 6 need not be stored if we are analyzing for loading which would excite only few global modes, i.e. we will just need k XG, where N is the number of modes required to be stored. This implies that the maximum storage required for the proposed algorithm is 3N2 + nG, where n is the maximum number of time steps the structure is analyzed for. Analyzing the computational time: for initial processing, we need to get the spectral decomposition of G with respect to A which would require O(N’) operations. Also, for kth time step, we need to do k cell updates of cost 3% operations each, and integration which requires k?? operations. Also, it involves the multiplication of Q with the integral of the transformed vector 6 on the characteristic. In summary. we would require an initial cost of O(N’) operations and for kth time step, we would require 4kz + NE.
M.N. Gudduti, J.L. Tassoulas / Cornput. Methods Appl. Mcch. Engrg. 167 (lYY8) 26/-27.1
169
If a convolution-based approach were used, we would need to store the Green’s function. This means that, if we like to store the Green’s function without resorting to interpolation, we would require space for nN’ numbers. Alternatively, we can store the Green’s function at fewer time points and interpolate in-between. The latter approach gives the flexibility of variable time steps, but increases the computation time considerably. In the former case of storing all the Green’s functions, we need the computation of Bessel functions at different time points in addition to the need for the spectral decomposition required in our approach. It is well known that the computation of Bessel functions is expensive and it is tricky to control the truncation error. For convolution purposes, the computational time required at kth time step is kN’ for the former and kN +kN’ for the latter. Here, k is the number of stored Green’s functions that influence the time until kth time step. The additional kN term comes from the need for interpolation because of the non-matching time grid. If we use the method proposed by Tsai et al., and if we use variable timestep size, we need to compute the integral of the Bessel function for N times for each step and the interval of integration keeps growing. Thus, the computational expense of computing the effect of past history is proportional to k%. The proportionality constant, however, is very high because of the computation of the Bessel functions. Also, the coefficient matrix (sr$+~~ nzatvir) is obtained by matrix multiplication for each time step and would require Ni operations. Neglecting the lower-order terms, we now compare the computational time required for each time step: for original Green’s function approach we need kN’ operations. For interpolated Green’s function approach we need kN’, with k growing with k. Modified Tsai’s method would require CkE +%N’. where C is the cost associated with Bessel function computation (certainly, C >> 4). The proposed approach would require 4kN + Nz. It is clear that the proposed method would in general be an order of magnitude faster than even the interpolated Green’s function approach. In fact, we can even do an analogue to the interpolated Green’s function approach: we can coarsen the characteristic cell granularity away from x = 0 which would give similar advantages. It is important to note here that we-can also truncate the number of modes depending on the expected type of loading. This would imply that N could be considerably smaller than N. We also mention the advantage of the proposed method in preserving the sparsity of the finite element system matrix of the interior, when the DtN map is coupled with it (in cases where the semidiscretization is performed using local approximations such as piecewise polynomials). The storage required in the proposed approach is minimal: only the history is stored. Whereas, in the conventional Green’s function approach, we need to store N X N matrices at different times which may be enormous. Another advantage is that we do not need to know the time range of analysis beforehand. In conclusion, the proposed method has distinct advantage over the other existing methods in terms of efficiency in both storage and computation.
6. Numerical
examples
In this section, the performance of the proposed method is illustrated with the use of some numerical examples. Two example problems are solved: one with single degree of freedom and another with two degrees of freedom. Both steady state (under time harmonic load) and transient analyses are performed to compare the results in both frequency and time domains. The frequency domain comparisons are made because the exact solutions in frequency domain are simple to compute. The time domain computations are performed to illustrate the effectiveness of the method when used with variable time step size. 6.1. Axial rod on elastic foundation First, the method is used to analyze an axial rod on elastic foundation (Fig. 3). This is generally considered as the model problem for out-of-plane shear problems as it is similar to a layer discretized with just one degree of freedom. The impulse response function is well known to be a Bessel function [2] and the frequency domain stiffness is a simple function [ 121: ( J1-w?). The method is used to perform steady state analysis of the system with the time step size At = min(T, T,)/6, where T is the time period of excitation and T, is the cut-off time period of the system. The results are compared with the exact frequency domain solution in Fig. 4. It can be seen that the method performs extremely well for all practical purposes. The smoothing effect at the points of discontinuity can be attributed to the leakage effect of sampling (discretization in time).
M.N. Guddati, J.L. Tassoulas I Comput. Methods Appl. Mech. Engrg. 167 (1998) 261-273
270
f
Fig. 3. One-dimensional tion.
test problem: axial bar on elastic founda-
Fig. 4. Comparison quency domain.
of the solution
with exact solution
in fre-
The method is then used for transient analysis of the same system under a square sine pulse (Fig. 5). The time stepping for the analysis is chosen in such a way that the there is little error in approximating the loading function and the upper bound is chosen considering the cut-off frequency. For this particular problem, the first 2 time-units are discretized using 10 steps of 0.2 and the remaining 18 time-units are discretized using 18 steps of unit size. The first time step size is governed by the prevention of aliasing, and the second by the cut-off frequency of the system. The results are compared in Fig. 6 with the results obtained by Fourier transforming the load and applying the frequency domain stiffness matrix. It is evident from the figure that the error is barely apparent. This illustrates the flexibility and accuracy of the method for general transient analysis. 6.2. Two DOF layer on a rigid rock Although it is clear that the extension to multiple degrees of freedom case is fairly straight forward, just for the sake of completeness, a two-degree-of-freedom model of a layer in antiplane shear is analyzed using this method. The shape functions defining the mode of deformation are assumed to be quadratic (Fig. 7). We first performed the steady state analysis and the results are compared in Figs. 8- 11 with the results of the consistent
2
1.5-A
-0.21 0
I 2
4
Fig. 5. Displacement analyses.
6
8
loading
10 time
function
12
14
16
18
20
-2
0
2
4
6
8
10
12
14
15
18
20
time
used
in the
transient
Fig. 6. Comparison domain.
of the solution
with exact solution
in time
271
M.N. Guddati, J.L. Tassoulas I Comput. Methods Appl. Mech. Engrg. 167 (1998) 261-273
Fig 7. Two DOF test problems: a rigid rock.
Out of plane motion of a layer on
Fig. 8. Comparison of the solution with the hyper element in frequency domain: Force at node I due to displacement at node 1.
boundaries [9]. The results compare well except at the cut-off frequencies. This is again attributed to the leakage effect. Transient analysis is performed on this system with the same displacement loading function as in the axial rod case (Fig. 5) with variable time steps. The time steps are chosen following similar guidelines as discussed in previous case: 10 steps of size 0.2 and 45 subsequent steps of size 0.4. The results from the consistent boundaries are used in conjunction with FFT to get the results for this transient analysis. The time resolution and frequency resolution are chosen fine enough to avoid any considerable approximation errors. The displacement is first applied at node 1 and the reacting force histories at both the nodes are compared with their counterparts from consistent boundary analysis (Figs. 12 and 13). Same comparison is performed with displacement being applied at node 2 (Figs. 14 and 15). Again, we can observe the superior performance of the method for transient analysis.
0
0.5
1
IS
2
2.5
3
3.5
4
4.5
5
Fig. 9. Comparison of the solution with the hyper element in frequency domain: Force at node 1 due to displacement at node 2.
0
0.5
1
1.5
2
2.5 fWuency
3
3.5
4
4.5
5
Fig. IO. Comparison of the solutions with the hyper element in frequency domain: Force at node 2 due to displacement at node 1.
212
M.N. Guddati, J.L. Tassoulas I Comput. Methods Appl. Mech. Engrg. 167 (1998) 261-273
DisplacementApplied at N&e 5.22: Fona at Node 2 due to Dlsdacantat Nodo 2
1
0.5 I‘ 0.4 -
_
0.3r
B fi
0.2.
K IL 0.1 I
O-
::i
2
4
6
8
10 time
12
14
16
1s
20
Fig. 12. Comparison of the solutions with the hyper element in time domain: Force at node 1 due to displacement at node I.
Fig. Il. Comparison of the solutions with the hyper element in frequency domain: Force at node 2 due to displacement at node 2.
Displacement Applied
0
at Node 1
Dispbement Applied at Node 2
-0.4
-0.5 I\/
-0.81 0
I
I 2
4
Fig. 13. Comparison time domain: Force
5
8
10 time
12
14
18
16
0
20
2
I 4
6
8
10
12
14
16
18
DisplacementApplied at Node 2
9
1.5.
I% p
0.5
8 9 O-
-0.5
-I’0 Fig. 15. Comparison
of the solutions
2
4
6
20
Fig. 14. Comparison of the solutions with the hyper element in time domain: Force at node 2 due to displacement at node I.
of the solutions with the hyper element in at node 1 due to displacement at node 2.
2-
w,
8
10 time
12
14
with the hyper element in time domain:
16
I8
Force at node 2 due to displacement
at node 2
M.N. Guddati, J.L. Tassoulas I Comput. Methods Appl. Mech. En~r8. 167 (199X) -761-27.1
7. Summary
273
and conclusions
In this paper, we propose a simple, but very efficient method of analysis of exterior scalar wave propagation problems in layered solids/fluids. This method can also be viewed as an efficient way of mesh extension as we are essentially dealing with a growing mesh in space-time. The method is simple to use and does not involve any functions that are complicated to compute numerically. It relies mainly on the semidiscretization idea in the vertical direction and the use of canonical transformation into characteristic coordinates. It is shown that the computational cost for this method is an order of magnitude less than the corresponding global methods. The proposed approach requires minimal storage. Also, if the semidiscretization is local, the resulting DtN map would be local in space and pseudo-local in time and thus preserve the sparsity of the interior coefficient matrix. The performance of the proposed method is illustrated with both steady state and transient analysis of single and multi-degree-of-freedom systems. The proposed method can be extended to cylindrical and spherical boundaries. but requires some modifications based on radial scaling, upwinding and exact integration (see [6] for details). The proposed method, however, can not be readily extended to the case of general elastic wave propagation. Also, it cannot be extended to scalar wave propagation in layered or inhomogeneous media. The reason for this limitation is that these problems have different characteristic speeds for different modes. At this time, the authors feel that it is not possible to extend this method efficiently. For these types of problems, a relatively less efficient technique based on a space-time finite element method is developed in [7]. The efficiency of the method can further be increased by a posteriori analysis of the gradient of the displacement vector and corresponding coarsening of the mesh. However, the coarsening should be performed with care to ensure that accuracy is not lost. The strategies for mesh-coarsening will be studied at a later time. Another venue of improvement of the efficiency is mesh truncation. As this method is essentially a mesh extension technique, it might be possible to truncate the mesh and use an appropriate treatment of the truncated boundary. The authors feel that the use of local boundary conditions in conjunction with mesh truncation can curtail the growth of the problem with time, but still retain good accuracy. This would be specifically advantageous in cylindrical and spherical geometries where the wave amplitude has a decay tendency. It is hoped that this method used with local boundary conditions would provide simple, efficient and accurate tool for the analysis of exterior wave propagation in homogeneous media.
References 111 ]2] [3] [4] (51 [6] [7] [8] ]9]
[IO] ] I 11
[ I21
A. Bedford and D.S. Drumheller, Introduction to Elastic Wave Propagation (John Wiley and Sons Ltd., Chichester. 1994). R. Courant and D. Hilbert, Methods of Mathematical Physics, Vol. II (Interscience Publishers, New York, 1966). B. Engquist and A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Math. Comput. 3 1 ( 1977) 629-65 I. P.R. Garabedian, Partial Differential Equations (John Wiley and Sons, New York, 1964). D. Givoli, Numerical Methods for Problems in Infinite Domains (Elsevier Science Publishers B.V., 1992). M.N. Guddati and J.L. Tassoulas, Characteristics methods for transient analysis of wave propagation m unbounded media, Comput. Methods Appl. Mech. Engrg. 164 (1998) 187-206. M.N. Guddati and J.L. Tassoulas, Space-time finite elements for the analysis of transient wave propagation in unbounded layered media, Int. J. Solids Struct., in press. E. Kausel, Local transmitting boundaries, ASCE J. Engrg. Mech. 114 (1988) IOIl- 1027. J.L. Tassoulas, Elements for the numerical analysis of wave propagation in layered media, Research Report R81-2. Dept. of Civil Engg., Order No. 689, Mass. Inst. of Tech., Cambridge, MA, 1981. C.S. Tsai and G.C. Lee, Arch dam-fluid interactions: by FEM-BEM and substructure concept, Int. J. Numer. Methods Engrg. 24 (1987) 236772388. C.S. Tsai, G.C. Lee and R.L. Ketter, A semi-analytical method for time-domain analyses of dam-reservoir interactions. Int. J. Numer. Methods Engrg. 24 (1990) 913-933. J.P. Wolf, Soil-Structure-Interaction in Time Domain (Prentice-Hall, Englewood Chffs, NJ. 1988).
JQzJ
Computer methods in applied mechanics and engineering
mm
mm
Comput.
EISEVIER
Methods
Appl. Mech. Engrg.
167 ( 1998) 261-273
An efficient numerical algorithm for transient analysis of exterior scalar wave propagation in a homogeneous layer Mu-thy N. Guddati”“, “Tercu Imtitutr
,jbr Cornputationctl hDepartmerzt
cud Applied
of Ci\?l Engineering,
Received
John L. Tassoulas b3*2
Mathematics,
The UniversiQ
The University
18 November
qf
cf Tern... ut Austin,
Texas ut Austin.
1996; revised
I8 August
Austin.
Austin.
TX 7X7/2- 1076, USA
TX 787/2- 1076. USA
1997
Abstract In this paper, a computationally efficient alternative for the analysis of transient wave propagation in an infinite layer is presented. The new method is based on a simple cell centered characteristic finite difference technique. This approach can be used in place of the existing methods based on closed form Green’s functions. as it would save considerable computational time and storage. Numerical examples are presented to illustrate the performance of the method. Some guidelines are suggested to further increase the efficiency of the proposed method. 0 1998 Elsevier Science S.A. All rights reserved.
1. Introduction Solution of a scalar wave equation in an infinite layer is generally encountered in fluid-structure interaction and soil-structure interaction problems. A typical example is the analysis of dam-reservoir interaction. The reservoir can be idealized as an infinite fluid layer fixed at the base. Another problem of similar type is the simulation of an elastic layer under out-of-plane motion. In this case the elastic vector wave equation reduces to a scalar wave equation and the procedure is essentially the same as the dam-reservoir interaction analysis. The elastic layer problem is considered as the model problem for the purposes of our discussion. Various techniques are used for the analysis of these types of problems [5,12]. The boundary element method is used for the analysis of the dam-water interaction problem by Tsai and Lee [lo]. They model the infinite surface and the interaction surface with boundary elements. The disadvantage of this approach is that the problem size is considerably large and requires the evaluation of singular integrals. An alternative to the boundary element method is the use of a Dirichlet-to-Neumann (DtN) map to replace the infinite layer [5]. Local DtN maps are sometimes used, but they are known to be crude and are good only for specific problems [3,8]. A more accurate way of modeling the effect of the infinite layer would be the use of a DtN map which is global in both time and space. Separation of variables is used for the derivation of the Green’s function for this system. The Green’s function consists of Bessel functions which are expensive to compute. The resulting Green’s function is used in conjunction with convolution to do transient wave propagation analysis of these types of systems. This approach demands high computational cost and the storage of the Green’s function. Semidiscretization techniques, which are equivalent to truncation of analytic Green’s functions. are often used to obtain a discrete Green’s function. Instead of Bessel functions, consistent
* Correspondmg author ’ Graduate student. ’ Professor. 00457825/98/$19.00 0 1998 Elsevier Science S.A. All rights reserved PII: SOO4S-7825(98)00123-6
262
M.N. Guddati. J.L. Tassoulas
/ Comput. Methods Appl. Mech. Engrg.
167 (1998) 261-273
transmitting boundaries can be used to calculate the dynamic stiffness matrix (see [9]) and Fourier transformation can be used to arrive at the impulse response function (Green’s function in time domain). These methods are extremely expensive as they involve either Fourier transformation or the computation of Bessel functions and their integrals. Tsai et al. [ 1 I] use modal decomposition to get the discrete Green’s function and the technique they propose is more efficient than the other existing methods. They, however, do not utilize the full decoupling which will reduce the computation time by an order of magnitude. Furthermore, in their approach, it is necessary to compute the integrals of Bessel functions. Another disadvantage of this method is that there is considerable loss of efficiency when the time step size is not constant. In this paper, an efficient alternative to the above methods is proposed. This method, based again on separation of variables and dimensional reduction technique, is considerably simpler to implement and extremely efficient. Semidiscrete equations in space-time are solved based on a characteristic finite difference technique and a symmetric positive definite DtN map is derived. There is no change in efficiency when the time steps are of variable size. The presentation of the rest of the paper is as follows: Section 2 states the model problem and the governing differential equations. In Section 3, dimensional reduction based on separation of variables technique is discussed. The derivation of the DtN map is carried out in Section 4 by means of the solution of the resulting system of hyperbolic equations. In Section 5, the computational cost of the current method is compared with the existing methods. Section 6 contains numerical examples which illustrate the performance of the method. We conclude the paper by discussing possible other applications and limitations of the proposed method in Section 7.
2. Problem
statement
Although the problem that is being solved here has broader motivation, we give a particular example of shear wave propagation in an infinite elastic layer. We have an elastic soil layer resting on a rigid rock and the loading is such that there is only out-of-plane deformation. The aim is to simulate the interaction between the layer and the structure embedded in or resting on it. The deformation is considered invariant with the out-of-plane position, i.e. the stress state in the soil is antiplane shear. The layer is assumed to be homogeneous. The general governing differential equation in the layer is the vector wave equation [I]: v* (C : Vu) - pu,, =f ,
(1)
where u is the three-dimensional displacement field, C is the elasticity tensor andf is the body force. With the restrictions described in the previous paragraph, the vector wave equation reduces to a scalar one: GAu-pu,,=f,
(2)
where u is the out-of-plane displacement field, G is the shear modulus andf is the corresponding source term resulting from the body force. The following procedure is used to analyze the structure (Fig. 1): a finite region of the soil bounded by straight vertical edges is modeled by finite elements. The rest of the infinite layer on either side of this region is simulated using a Dirichlet-to-Neumann (DtN) map. The aim of this work is to get an appropriate time dependent DtN map for general transient analysis of this system. The problem can be stated as follows: find the relation DtN : u,. -+s, in which u,. : f X I + !‘H is the displacement on the boundary r at time t E I= [0, T,,,] and s : r X I + % is the traction on the boundary. The applicable differential equation and the imposed boundary and initial conditions on the layer are given by Ghu-pu,,=f
forxE0,
tEZ,
where R = !h’ X (-H, 0) is the domain of the soil layer. Homogeneous Dirichlet condition bottom of the layer (z = -H) and Neumann condition on the top of the layer (z = 0): u(x, z, t) = 0
and
at z = -H
(3) is applied on the
(4)
M.N.
Guddati,
J.L. Tassoulas
I Comput.
Methods
Appl.
Mech.
Engrg.
167 (1998)
261-27.1
263
Soil Layer Rock
(a) Physical Problem
DtN map
DtN map
Bounded Soil (FJN Fixed Base (b) Mathematical Model Fig. I. General
u:(x, z, t) = so(x, t) The soil is assumed
media.
(5)
to start from rest: at t = 0 .
reduction
(6)
of the problem
The separation of variables technique field by the following approximation: 2 = rt
for unbounded
at z = 0.
u(x, 2, t) = u,(x, z, t) = 0
3. Dimensional
method of analysis
4,(Z)U,(X, f) = QT(z)U(4
is used to reduce the problem dimension.
t) .
We discretize the displacement
(7)
In the above form of approximation, a(z) can be viewed as the basis function vector in the ; direction and U(X, t) is the nodal degrees of freedom vector. We apply the standard Galerkin method using this approximation. The test function is of the same form as the solution: uh = ,; The Galerkin
&(Z)U,(X, r) = @T(z)V(X, t) . approximation
of the variational
u”(G AU” - p$, -f)
Integrating
dz = 0
by parts with respect to z, we get
(8) form for Eq. (3) is (9)
264
M.N. Gudduti, J.L. Tussoulas
ad au”
cl
I
-G Substituting
a2uh
--dz+G
a2uh
v”--dz
u”Tdz-p
az az
-H
I Comput. Methods Appl. Mech. Engrg. 167 (1998) 261-27.3
-
uhf dz=O.
at2
the approximation
(10)
for uh and u”, we get 0
-GVT
@$j;
dz U + GVT
I
c&PDTdz U.,,
-H
0
I
-pvT
QWT
dz U,, -VT
-H
which reduces to the following
I
” wdz+GVT@\;=,s,=O
vVVEMN
If
(11)
form:
-GU+AU,,-L*U,,-F=Q, c2
(12)
where c =J/G/p is the wave speed. In the above equation, G and A are constant N X N matrices and U and F are vector functions of x and t. G, A and F are given by
(13)
(14) and F(x, t) =
(15)
Note that the degree of freedom corresponding to the base is eliminated and the matrices G and A turn out to be symmetric and positive definite. Our modified problem is to find a map between the displacement vector U, and the traction vector consistent to the approximation made in displacement, i.e. if S is the traction vector, we need to satisfy the work equivalence:
(16) where sI, is the approximation relation
of the traction
s. As s = -Gu,
we have s,, = - GTGIJ,
and we will have a
0
I
VTS = -VT
PN
@GaTdzU,
This implies that the traction S = -AU,
b’VEaN.
vector is nothing
(17) but
.
(18)
We now present the discrete version of our original problem: Find the map DtN : U, -+ S where U, : I + sHN is the displacement vector at x = 0 and S : I + SN is the corresponding traction vector. We satisfy the following differential equations and initial conditions: AUxx-LAU,,-GU-F=O
forxE%+,
tEZ
(19)
C2
U=U,=Q
forxE%+,
t=O.
(20)
M.N. Gudduti, J.L. Tussoulas I Comput. Methods Appl. Mec,h. Engrg. 167 (1998) 261-273
26.5
Derivation of the DtN map In this section, the hyperbolic nature of the system of equations derived from the semidiscretization procedure used to arrive at an efficient DtN map. We first decouple the system of equations by spectral decomposition of with respect to A. The system can be written as QQ’~.,-~QQ’c!,-QAQ’u-F=U
(21)
where Q is the modal matrix and A is the diagonal Premultiplying this equation with Q-’ gives
matrix containing
the eigenvalues
o,,--L q,,- /@-g =0 c2
of G with respect to A.
(22)
in which i?=Q’U and @=Q-‘F. As an aside, we note here that the above system of equations are uncoupled second-order hyperbolic equations. This fact can be used effectively to increase both time and storage efficiencies when using discrete Green’s functions. If we use the approach based on Bessel functions as in [ 111, we can store the values of cinstead of U and increase the computational efficiency by an order of magnitude. We, however, do not utilize this modification as it is not necessary to resort to the use of Bessel functions in the proposed approach. We note that the uncoupled system of equations (22) have the same characteristic speed. This gives us the advantage of transforming all the equations into standard canonical form using the same transformation. For all the equations, we substitute .Y+ ct t=2 The canonical
and
x - ct 7 = 2 .
(23)
form will then look like
ir,,-n&F=0 Our objective is to compute rule, we have
(24) l?, at (x, t) = (0, 7’,,) or, equivalently,
at (5,~) = ([,,, v,,) (see Fig. 2). By the chain
(25)
‘t
_.-C
* q=-cT,/2 11= -CT n-, I2
‘-“1;--i_.g; Fig. 2. Finite difference
grid for updating
0.
266
M.N. Guddati, J.L. Tassoulas / Comput. Methods Appl. Mech. Engrg. 167 (1998) 261-273
Integrating
Eq. (24) along the characteristic
On account of Eq. (25)
5 = &, (with respect to v), we find
we obtain
It is important to note that 0 and therefore fit can be computed along the &characteristic at v= 0 provided that the source term (8) and the initial conditions on fi (at t = 0, x 2 0) are known. This follows from the fact that the domain of dependence [2,4] of any point on the v = 0 characteristic is contained in the line {X2 0, t = 0). The Cauchy data (6, 6,) on this line, together with the source term (P) uniquely determine the value of 6 on the characteristic 7 = 0 (see [2,4]). This in turn determines fit at r] = 0. Thus, the computation to arrive at o& &, 0) is carried out, in many cases explicitly, prior to the initiation of the transient analysis. For a medium starting from rest (0 = tit = 0) with vanishing source term (P = 0), it is clear that 0 = c6 = 0 along 7 = 0. In this case, we can replace Eq. (27) by
We now divide the space-time domain into cells which have the edges parallel to the characteristics (Fig. 2). fi is assumed to be bilinear in each cell. We then require the values of 6 at the grid points only. We use the following notation for the purposes of presentation: 0: is the value of e at (5, q) = ([,, Q) = (cT,/2, -cTk /2). The interval between T,, ~, and T, is considered as the nth interval. We now describe the procedure for obtaining the values of IJ?on the v characteristic corresponding to the current time T,,, i.e. 5 = [,,. The set of active cells in this step are shown in Fig. 2. The ordering of the analysis is also shownk in the same, figAu;e: we stA$rom cell (n, 1) and proceed till (n, n - 1). For each cell (n, k), we need to obtain U,,, given U,,_, , U,j_, and U,, . We use cell centered finite difference technique for this purpose. We satisfy the differential equation at the center of the cell:
(29) As we have A&,, = c AT, 12 and Aqk = -c AT, 12 (we are sweeping following form:
in negative
77 direction),
we arrive at the
o:(i)=a;@:_,(i) +ii;-‘(i), - o;::(i) + p;Fyq, where 4
-- A(i) c’ AT,, AT, 4 cu,=
4
+ A(i) ’
c2bTn AT,
4
(31)
and
P;=
41 +4) c2 AT,, AT,
’
(32)
4
In the above equations, 6(i) and P(i) stand for the ith element of 8 and P, respectively. A(i) is the ith diagonal element of the matrix A. Note here that the coefficient cy, is always less than unity for any mode. It is expected that this fact is crucial for the stability of the method. A proof of stability is not available at this time. Numerical experiments indicate that the method is stable, as long as 4 > 1. A general proof of stability is currently under investigation.
M.N. Gudduti, .I.L. TassoulaJ I Comput. Methods Appl. Mech. Eqrg.
167 (1998) 2hlk-77.1
267
Thus, with l? known on the previous characteristic and at v= 0, we can get 6 on the new characteristic using Eq. (30). Once 6 is known on this characteristic, we use standard numerical integration to compute the integral in Eq. (27). The most appealing integration procedure would be the one which is consistent with the approximation of 0 over the mesh. This suggests the piecewise linear interpolation on the characteristic and we thus use trapezoidal rule for numerical integration. It should be noted that higher order difference schemes or finite element schemes could also be used for the updating, and consistent integration schemes could be used for the integration. Using the trapezoidal rule, Eq. (27) can be replaced by the following summation
(33) which can be rewritten
as
(34) We premultiply
the above equation
by -GQ
to get
(35) since A = QQ’ and G = QAQ’. It can be seen that the left-hand side of this equation is nothing but the traction vector that we are interested in. We thus have the following relation between the traction and displacement at x = 0: S=KU+CZ+H
at
(36)
3
where c AT K= 4G:
(37)
C=+A,
(38)
and c(AT,
+
AT,, ,)
(AO;;+$;)+yF,s
4
CAT
no
-A-$. I
au Cl (39)
Notice that both the coefficient matrices K and M are symmetric and positive definite. Thus, for conventional time-stepping schemes, we get a symmetric positive definite coefficient matrix from the boundary which has definite advantage for linear solution purposes. Another advantage, when the discretization is local, is that the DtN map in Eq. (36) preserves the sparsity of the interior coefficient matrices. This is due to the fact that the matrices G and A are sparse when the interpolation function @J is sparse (see Eqs. (13) and (14)). Before we conclude this section, we mention an alternative way of obtaining the DtN map. In contrast to Eq. (25), we can use fi,, to compute o,r: (40) l?, at (x, t) = (0, T,,) can be obtained
locally by a backward
difference
approximation:
M.N. Guddati, J.L. Tassoulas I Comput. Methods Appl. Mech. Engrg. 167 (1998) 261-273
268
(41)
Higher-order the resulting
difference approximations could be used as well. Skipping relationship between the traction and displacement:
all the intermediate
details, we present
(42)
where (43)
and Hz-
&Qc::-,
.
(45)
Note that this form does not involve the cost of computing the integral. However, it has several disadvantages. The most important one is the possible lack of stability: the matrix c is negative definite and when augmented with the interior coefficient matrices, it can lead to an unstable analysis procedure because of the loss in positive definiteness. Furthermore, it can be shown that the computation of traction from Eqs. (36)-(39) i.e. by means of integration along the characteristic, is second-order accurate. On the other hand, first-order difference approximations, like the one in Eq. (4) limit the accuracy of the overall scheme. Higher-order difference approximations can circumvent this difficulty, but require greater smoothness of the field variable. The boundary condition in Eqs. (36)-(39) does not suffer from these drawbacks and we choose to use it in spite of the need to evaluate the integrals. To summarize this section, we have derived a DtN map which relates the displacement vector on the boundary with the corresponding traction vector. The relationship is global in space, but local whenever the interpolation functions used for the semidiscretization are local. However, it is pseudo-local in time as it uses only the displacement vector and the time derivative in a direct sense. The past history is stored in a different form with the displacement vector stored along the characteristics. It should be noted that, unlike methods with higher order derivatives, this method is expected to be stable as it relies only on the first time derivative and integration. We now go ahead with the presentation of the computational cost after some comments on implementation issues.
5. Computational cost From the previous section it is clear that the values of t?? on the current characteristic, and the matrices G, A and Q must be stored. For a general setting (i.e. assuming the matrices G and A are dense), the initial required storage would be 3N2 where N is the number of degrees of freedom on the boundary. When piecewise polynomials are used for the purposes of semidiscretization, the matrices G and A are banded and the storage would be comparatively less. The storage required for the history (6 on the characteristic) at kth time step is given by k X N. Also, it is noted that the full vector 6 need not be stored if we are analyzing for loading which would excite only few global modes, i.e. we will just need k XG, where N is the number of modes required to be stored. This implies that the maximum storage required for the proposed algorithm is 3N2 + nG, where n is the maximum number of time steps the structure is analyzed for. Analyzing the computational time: for initial processing, we need to get the spectral decomposition of G with respect to A which would require O(N’) operations. Also, for kth time step, we need to do k cell updates of cost 3% operations each, and integration which requires k?? operations. Also, it involves the multiplication of Q with the integral of the transformed vector 6 on the characteristic. In summary. we would require an initial cost of O(N’) operations and for kth time step, we would require 4kz + NE.
M.N. Gudduti, J.L. Tassoulas / Cornput. Methods Appl. Mcch. Engrg. 167 (lYY8) 26/-27.1
169
If a convolution-based approach were used, we would need to store the Green’s function. This means that, if we like to store the Green’s function without resorting to interpolation, we would require space for nN’ numbers. Alternatively, we can store the Green’s function at fewer time points and interpolate in-between. The latter approach gives the flexibility of variable time steps, but increases the computation time considerably. In the former case of storing all the Green’s functions, we need the computation of Bessel functions at different time points in addition to the need for the spectral decomposition required in our approach. It is well known that the computation of Bessel functions is expensive and it is tricky to control the truncation error. For convolution purposes, the computational time required at kth time step is kN’ for the former and kN +kN’ for the latter. Here, k is the number of stored Green’s functions that influence the time until kth time step. The additional kN term comes from the need for interpolation because of the non-matching time grid. If we use the method proposed by Tsai et al., and if we use variable timestep size, we need to compute the integral of the Bessel function for N times for each step and the interval of integration keeps growing. Thus, the computational expense of computing the effect of past history is proportional to k%. The proportionality constant, however, is very high because of the computation of the Bessel functions. Also, the coefficient matrix (sr$+~~ nzatvir) is obtained by matrix multiplication for each time step and would require Ni operations. Neglecting the lower-order terms, we now compare the computational time required for each time step: for original Green’s function approach we need kN’ operations. For interpolated Green’s function approach we need kN’, with k growing with k. Modified Tsai’s method would require CkE +%N’. where C is the cost associated with Bessel function computation (certainly, C >> 4). The proposed approach would require 4kN + Nz. It is clear that the proposed method would in general be an order of magnitude faster than even the interpolated Green’s function approach. In fact, we can even do an analogue to the interpolated Green’s function approach: we can coarsen the characteristic cell granularity away from x = 0 which would give similar advantages. It is important to note here that we-can also truncate the number of modes depending on the expected type of loading. This would imply that N could be considerably smaller than N. We also mention the advantage of the proposed method in preserving the sparsity of the finite element system matrix of the interior, when the DtN map is coupled with it (in cases where the semidiscretization is performed using local approximations such as piecewise polynomials). The storage required in the proposed approach is minimal: only the history is stored. Whereas, in the conventional Green’s function approach, we need to store N X N matrices at different times which may be enormous. Another advantage is that we do not need to know the time range of analysis beforehand. In conclusion, the proposed method has distinct advantage over the other existing methods in terms of efficiency in both storage and computation.
6. Numerical
examples
In this section, the performance of the proposed method is illustrated with the use of some numerical examples. Two example problems are solved: one with single degree of freedom and another with two degrees of freedom. Both steady state (under time harmonic load) and transient analyses are performed to compare the results in both frequency and time domains. The frequency domain comparisons are made because the exact solutions in frequency domain are simple to compute. The time domain computations are performed to illustrate the effectiveness of the method when used with variable time step size. 6.1. Axial rod on elastic foundation First, the method is used to analyze an axial rod on elastic foundation (Fig. 3). This is generally considered as the model problem for out-of-plane shear problems as it is similar to a layer discretized with just one degree of freedom. The impulse response function is well known to be a Bessel function [2] and the frequency domain stiffness is a simple function [ 121: ( J1-w?). The method is used to perform steady state analysis of the system with the time step size At = min(T, T,)/6, where T is the time period of excitation and T, is the cut-off time period of the system. The results are compared with the exact frequency domain solution in Fig. 4. It can be seen that the method performs extremely well for all practical purposes. The smoothing effect at the points of discontinuity can be attributed to the leakage effect of sampling (discretization in time).
M.N. Guddati, J.L. Tassoulas I Comput. Methods Appl. Mech. Engrg. 167 (1998) 261-273
270
f
Fig. 3. One-dimensional tion.
test problem: axial bar on elastic founda-
Fig. 4. Comparison quency domain.
of the solution
with exact solution
in fre-
The method is then used for transient analysis of the same system under a square sine pulse (Fig. 5). The time stepping for the analysis is chosen in such a way that the there is little error in approximating the loading function and the upper bound is chosen considering the cut-off frequency. For this particular problem, the first 2 time-units are discretized using 10 steps of 0.2 and the remaining 18 time-units are discretized using 18 steps of unit size. The first time step size is governed by the prevention of aliasing, and the second by the cut-off frequency of the system. The results are compared in Fig. 6 with the results obtained by Fourier transforming the load and applying the frequency domain stiffness matrix. It is evident from the figure that the error is barely apparent. This illustrates the flexibility and accuracy of the method for general transient analysis. 6.2. Two DOF layer on a rigid rock Although it is clear that the extension to multiple degrees of freedom case is fairly straight forward, just for the sake of completeness, a two-degree-of-freedom model of a layer in antiplane shear is analyzed using this method. The shape functions defining the mode of deformation are assumed to be quadratic (Fig. 7). We first performed the steady state analysis and the results are compared in Figs. 8- 11 with the results of the consistent
2
1.5-A
-0.21 0
I 2
4
Fig. 5. Displacement analyses.
6
8
loading
10 time
function
12
14
16
18
20
-2
0
2
4
6
8
10
12
14
15
18
20
time
used
in the
transient
Fig. 6. Comparison domain.
of the solution
with exact solution
in time
271
M.N. Guddati, J.L. Tassoulas I Comput. Methods Appl. Mech. Engrg. 167 (1998) 261-273
Fig 7. Two DOF test problems: a rigid rock.
Out of plane motion of a layer on
Fig. 8. Comparison of the solution with the hyper element in frequency domain: Force at node I due to displacement at node 1.
boundaries [9]. The results compare well except at the cut-off frequencies. This is again attributed to the leakage effect. Transient analysis is performed on this system with the same displacement loading function as in the axial rod case (Fig. 5) with variable time steps. The time steps are chosen following similar guidelines as discussed in previous case: 10 steps of size 0.2 and 45 subsequent steps of size 0.4. The results from the consistent boundaries are used in conjunction with FFT to get the results for this transient analysis. The time resolution and frequency resolution are chosen fine enough to avoid any considerable approximation errors. The displacement is first applied at node 1 and the reacting force histories at both the nodes are compared with their counterparts from consistent boundary analysis (Figs. 12 and 13). Same comparison is performed with displacement being applied at node 2 (Figs. 14 and 15). Again, we can observe the superior performance of the method for transient analysis.
0
0.5
1
IS
2
2.5
3
3.5
4
4.5
5
Fig. 9. Comparison of the solution with the hyper element in frequency domain: Force at node 1 due to displacement at node 2.
0
0.5
1
1.5
2
2.5 fWuency
3
3.5
4
4.5
5
Fig. IO. Comparison of the solutions with the hyper element in frequency domain: Force at node 2 due to displacement at node 1.
212
M.N. Guddati, J.L. Tassoulas I Comput. Methods Appl. Mech. Engrg. 167 (1998) 261-273
DisplacementApplied at N&e 5.22: Fona at Node 2 due to Dlsdacantat Nodo 2
1
0.5 I‘ 0.4 -
_
0.3r
B fi
0.2.
K IL 0.1 I
O-
::i
2
4
6
8
10 time
12
14
16
1s
20
Fig. 12. Comparison of the solutions with the hyper element in time domain: Force at node 1 due to displacement at node I.
Fig. Il. Comparison of the solutions with the hyper element in frequency domain: Force at node 2 due to displacement at node 2.
Displacement Applied
0
at Node 1
Dispbement Applied at Node 2
-0.4
-0.5 I\/
-0.81 0
I
I 2
4
Fig. 13. Comparison time domain: Force
5
8
10 time
12
14
18
16
0
20
2
I 4
6
8
10
12
14
16
18
DisplacementApplied at Node 2
9
1.5.
I% p
0.5
8 9 O-
-0.5
-I’0 Fig. 15. Comparison
of the solutions
2
4
6
20
Fig. 14. Comparison of the solutions with the hyper element in time domain: Force at node 2 due to displacement at node I.
of the solutions with the hyper element in at node 1 due to displacement at node 2.
2-
w,
8
10 time
12
14
with the hyper element in time domain:
16
I8
Force at node 2 due to displacement
at node 2
M.N. Guddati, J.L. Tassoulas I Comput. Methods Appl. Mech. En~r8. 167 (199X) -761-27.1
7. Summary
273
and conclusions
In this paper, we propose a simple, but very efficient method of analysis of exterior scalar wave propagation problems in layered solids/fluids. This method can also be viewed as an efficient way of mesh extension as we are essentially dealing with a growing mesh in space-time. The method is simple to use and does not involve any functions that are complicated to compute numerically. It relies mainly on the semidiscretization idea in the vertical direction and the use of canonical transformation into characteristic coordinates. It is shown that the computational cost for this method is an order of magnitude less than the corresponding global methods. The proposed approach requires minimal storage. Also, if the semidiscretization is local, the resulting DtN map would be local in space and pseudo-local in time and thus preserve the sparsity of the interior coefficient matrix. The performance of the proposed method is illustrated with both steady state and transient analysis of single and multi-degree-of-freedom systems. The proposed method can be extended to cylindrical and spherical boundaries. but requires some modifications based on radial scaling, upwinding and exact integration (see [6] for details). The proposed method, however, can not be readily extended to the case of general elastic wave propagation. Also, it cannot be extended to scalar wave propagation in layered or inhomogeneous media. The reason for this limitation is that these problems have different characteristic speeds for different modes. At this time, the authors feel that it is not possible to extend this method efficiently. For these types of problems, a relatively less efficient technique based on a space-time finite element method is developed in [7]. The efficiency of the method can further be increased by a posteriori analysis of the gradient of the displacement vector and corresponding coarsening of the mesh. However, the coarsening should be performed with care to ensure that accuracy is not lost. The strategies for mesh-coarsening will be studied at a later time. Another venue of improvement of the efficiency is mesh truncation. As this method is essentially a mesh extension technique, it might be possible to truncate the mesh and use an appropriate treatment of the truncated boundary. The authors feel that the use of local boundary conditions in conjunction with mesh truncation can curtail the growth of the problem with time, but still retain good accuracy. This would be specifically advantageous in cylindrical and spherical geometries where the wave amplitude has a decay tendency. It is hoped that this method used with local boundary conditions would provide simple, efficient and accurate tool for the analysis of exterior wave propagation in homogeneous media.
References 111 ]2] [3] [4] (51 [6] [7] [8] ]9]
[IO] ] I 11
[ I21
A. Bedford and D.S. Drumheller, Introduction to Elastic Wave Propagation (John Wiley and Sons Ltd., Chichester. 1994). R. Courant and D. Hilbert, Methods of Mathematical Physics, Vol. II (Interscience Publishers, New York, 1966). B. Engquist and A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Math. Comput. 3 1 ( 1977) 629-65 I. P.R. Garabedian, Partial Differential Equations (John Wiley and Sons, New York, 1964). D. Givoli, Numerical Methods for Problems in Infinite Domains (Elsevier Science Publishers B.V., 1992). M.N. Guddati and J.L. Tassoulas, Characteristics methods for transient analysis of wave propagation m unbounded media, Comput. Methods Appl. Mech. Engrg. 164 (1998) 187-206. M.N. Guddati and J.L. Tassoulas, Space-time finite elements for the analysis of transient wave propagation in unbounded layered media, Int. J. Solids Struct., in press. E. Kausel, Local transmitting boundaries, ASCE J. Engrg. Mech. 114 (1988) IOIl- 1027. J.L. Tassoulas, Elements for the numerical analysis of wave propagation in layered media, Research Report R81-2. Dept. of Civil Engg., Order No. 689, Mass. Inst. of Tech., Cambridge, MA, 1981. C.S. Tsai and G.C. Lee, Arch dam-fluid interactions: by FEM-BEM and substructure concept, Int. J. Numer. Methods Engrg. 24 (1987) 236772388. C.S. Tsai, G.C. Lee and R.L. Ketter, A semi-analytical method for time-domain analyses of dam-reservoir interactions. Int. J. Numer. Methods Engrg. 24 (1990) 913-933. J.P. Wolf, Soil-Structure-Interaction in Time Domain (Prentice-Hall, Englewood Chffs, NJ. 1988).