Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Three applications of the reciprocal theorem in soil-structure interaction
(USC Thesis Other)
Three applications of the reciprocal theorem in soil-structure interaction
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
THREE APPLICATIONS OF THE RECIPROCAL THEOREM IN SOIL-STRUCTURE INTERACTION by Kirsten McKay A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (CIVIL ENGINEERING) May 2009 Copyright 2009 Kirsten McKay Acknowledgements First and foremost I want to thank my Heavenly Father for giving me everything I have and am. I am forever indebted to my advisor, H.L. Wong, for providing me with the opportunity, motivation, and means to come to this point. His contribution cannot be quantified. I wish to thank Professor J.E. Luco and Dr. R.J. Apsel for the use of their Green’s Function program, without it, this research cannot be completed. I am also deeply grateful to my parents for their unceasing love and support in this and every endeavor I pursue. ii Table of Contents Acknowledgements........................... ii List of Figures............................. v Abstract ...............................vii Chapter 1 – Introduction ......................... 1 Chapter 2 – Applications of the Reciprocal Theorem ............. 8 2.1 – Application One, A Body Force Formulation ............ 12 2.2 – Application Two, A Direct Boundary Integral Equation Method ..... 14 2.3 – Application Three, A Radiation Boundary for Discrete Methods ..... 17 Chapter3–A Body Force Formulation ................... 23 3.1 – Numerical Procedure for Integration of Singularities ......... 28 3.2 – Application of Green’s Functions Obtained in Cylindrical Coordinates . . 29 3.3 – Convergence Rate of of Distributed Force Method .......... 32 3.3.1 – Solution with a Point-Displacement Compatibility Condition ..... 33 3.3.2 – Solution with an Average-Displacement Compatibility Condition . . . 37 3.3.3 – Convergence of the Dynamic Solution .............. 40 3.4 – Comparison with Published Results ................ 44 3.5 – Results of a Square Embedded Foundation ............. 53 3.6 – Response of a Massless Embedded Foundation to Incident Waves .... 57 Chapter 4 – Application Two, the Boundary Integral Equation Method ...... 62 4.1 – Numerical Procedure for the Integration of Singularities ........ 66 4.2 – Application of Green’s Functions Obtained in Cylindrical Coordinates . . 67 4.3–ATwo-Dimensional SH-Wave Example .............. 72 4.4 – Clues to Numerical Instability Problem............... 76 4.4.1 – Application of l’Hospital’s Rule ................. 78 Chapter 5 – Application Three, a Radiation Boundary ............. 80 5.1 – Three-Dimensional Finite Element Formulation ........... 82 5.2 – Creation of a Radiation Boundary ................. 84 5.3 – A Modal Solution ....................... 89 5.3.1 – Richardson Extrapolation for the Missing Submatrix ........ 93 Chapter 6 – Conclusion ......................... 95 References .............................. 99 iii AppendixA–A USGS Damage Report of the Olive View Hospital Building . . . 103 Appendix B – Example of the Reciprocal Theorem ............. 106 B.1 Displacements and Stresses Generated by a Vertical Point Load .... 107 B.1.1 – Displacements ...................... 108 B.2 Displacements and Stresses Generated by a Point Load in theX-Direction 109 B.2.1 – Displacements ...................... 111 B.3 Displacements and Stresses Generated by a Point Load in theY -Direction 111 B.3.1 – Displacements ...................... 113 AppendixC–V olume Integral of Static Singularities ............ 114 AppendixD–Wave Motion in a Semi-Infinite Medium ........... 118 D.1 – Two-Dimensional Plane-Strain Problem ............. 118 D.2 – Incident P-Wave Solution ................... 119 D.2.1 – Wavefields for P-wave Incidence ............... 124 D.3 – Incident SV-Wave Problem .................. 125 D.3.1 – Wavefields for SV-wave Incidence............... 129 D.4 – Incident SH-Wave Solution .................. 131 D.4.1 – Wavefields for SH-wave Incidence .............. 132 D.5 – Rayleigh Wave Solution ................... 132 Appendix E – Three-Dimensional Isoparametric Elements .......... 135 iv List of Figures Figure 2.1 – Loading Conditions of a Solid Body. ............... 8 Figure 2.2 – Body Force Distribution Inside an Embedded Foundation. ...... 13 Figure 2.3 – A Schematic of the Direct Boundary Integral Equation Method. .... 15 Figure 2.4 – Formulation of a Radiation Boundary. .............. 18 Figure 3.1 – Distributed Forces Inside an Embedded Foundation.......... 24 Figure 3.2 – Six Degree of Freedoms for a Rigid Embedded Foundation....... 26 Figure 3.3 – The Definition of the Cylindrical Coordinate System. ........ 31 Figure 3.4 – A Schematic for aRigid Embedded Foundation. .......... 34 Figure 3.5 – Performance of the Distributed Force Method for Various Embedment Ratios h/b. .............. 36 Figure 3.6 – Convergence of Impedance Functions for h/b=1 and h/b = 2. ..... 38 Figure 3.7 – Convergence Rate of Average Displacement Compatability versus Point-Displacement Compatability. ............. 39 Figure 3.8 – Convergence Trend of Average-Displacement Method......... 41 Figure 3.9 – An Approximate Grid for a Cylindrical Foundation. ......... 44 Figure 3.10 – Vertical and Torsional Impedances for Various Poisson’s Ratios. . . . 45 Figure 3.11 – Horizontal, Rocking and Coupling Impedances for Various Poisson’s Ratios. ............... 47 Figure 3.12 – A Closer Observation of the Horizontal Impedance. ........ 48 Figure 3.13 – Comparisons ofK vv andK tt to Published Results. ........ 49 Figure 3.14 – Comparisons ofK hh ,K mm andK mh to Published Results. ..... 51 Figure 3.15 –K vv andK tt for a Square Embedded Foundation. ......... 53 Figure 3.16 –K hh ,K mm andK mh for a Square Embedded Foundation....... 54 v Figure 3.17 – Horizontal, Vertical and Rocking Input Motion Due to Incident SH-Waves. ................ 56 Figure 3.18 – Horizontal, Vertical and Rocking Input Motion Due to Incident P-Waves................... 58 Figure 3.19 – Horizontal, Vertical and Rocking Input Motion Due to Incident SV-Waves. ................ 60 Figure 4.1 – Subdivision of the Rigid Embedded Foundation Surface. ....... 63 Figure 4.2 – Definition of the Cylindrical Coordinate System. .......... 68 Figure 4.3–ATwo-Dimensional Circular Cylindrical Foundation. ........ 74 Figure 4.4 – Antiplane Impedance Function as a Function of Frequency. ...... 75 Figure 5.1 – A Three-Dimensional Finite Element Model. ........... 81 Figure 5.2 – Nodal Definition for Radiation Boundary. ............. 82 Figure 5.3 – Simplified Modal Analysis Model. ............... 90 Figure 5.4–APartially Restrainted Finite Element Model. ........... 94 Figure B.1 – VerticalQ Z Point Force Configuration. ............ 107 Figure B.2 – HorizontalQ X Point Force Configuration. ........... 110 Figure B.3 – HorizontalQ Y Point Force Configuration. ........... 112 Figure D.1 – Incident and reflected P-wave. ................ 119 Figure D.2 – Incident and reflected SV-wave. ............... 126 Figure D.3 – Incident and reflected SH-wave. ............... 131 Figure E.1 – Node Numbering Scheme for a Three-Dimensional Solid Element. . . 136 vi Abstract Using the Reciprocal Theorem and the principle of superposition, this dissertation examines three procedures for obtaining the impedance functions and driving force vectors necessary for a soil-structure interaction analysis. The first procedure employs a volume integral of body forces and the displacement Green’s Functions to obtain the frequency dependent force-displacement relationship (impedance functions) for a rigid embedded foundation and the generalized force vector required to hold the rigid foundation fixed while subjected to incident waves (driving force). Although a formulation based on a volume integral requires more computation than the others based on surface integrals, this procedure proves to be numerically stable and is an excellent tool for linear soil-structure interaction analyses. The second procedure is based on integrating the surface displacements with the traction Green’s Functions and surface tractions with the displacement Green’s Function, over the surface of a rigid foundation. This well-known formulation has numerical instability problems at the resonant frequencies of the medium interior to the foundation surface. This dissertation explains the reason behind the instability using an analytical solution of the Hilbert-Schmidt method and recommends an algorithm which uses the l’Hospitale Rule to obtain the solution at those critical frequencies. The development herein for this method is not meant to be a practical solution, but an academic exercise to explain a long-standing, puzzling problem. The third and final procedure is to design a radiation boundary for a three-dimensional finite element grid using the surface displacements and surface tractions at the discrete artificial boundary, along with the Green’s Functions, to calculate outgoing wave motion at nodes immediately outside of the artificial boundary, thereby eliminating unwanted reflection of the outgoing waves. A procedure is recommended to obtain required data from a commercially available software program. A modal method can be employed to maximize the efficiency of the analysis. vii Chapter One Introduction Most earthquake engineering researchers agree that a soil-structure interaction (SSI) analysis is essential for designing structures with a heavy, voluminous foundation (Choi et al, 2001, Gueguen et al, 2002). One example of this type of structure is a nuclear power plant which has a massive underlying containment unit which prevents leakage of nuclear byproducts. When soil-structure interaction first began to be studied, little emphasis was attached to an SSI analysis for high-rise buildings because they were viewed as light- weight structures. The earthquake response of these types of structures has generally been analyzed without considering SSI effects. This view is changing because damage assessment reports of buildings in recent, destructive earthquakes such as the 1985 Mexico City, 1989 Loma Prieta, and 1994 Northridge earthquakes have shown that SSI effects are important. These tremblors left many buildings with significant damage, but structural engineers were impressed by the fact that the damage was not more severe. Duan (1999), USGS (1996) and others have suggested that SSI is the reason for this surprising result. At the same time, however, the 1995 devastation of the Kobe earthquake was amplified by SSI effects (Mylonakis and Gazetas, 2000). This opens a great debate about what variables control SSI and further augments the need to verify a method for including soil-structure interaction in the analysis of all buildings (Jennings and Bielak, 1973). In the past, engineers have considered the exclusion of a soil-structure interaction analysis to be a conservative practice. But as more and more evidence of the importance of soil-structure interaction is collected, many researchers and practicing engineers are requesting that a thorough soil-structure interaction analysis be part of municipal building codes (e.g. Stewart et al, 2004). There have been previous movements towards including SSI effects in building codes. The first of these was the 1978 ATC-3 (Applied Technology 1 Council) code. Refinements have been made over the years. In 1996 the ATC-40 code made some provisions for a flexible supporting structure. A year later, the National Earthquake Hazards Reduction Program (NEHRP) adopted a simplified sub-structuring technique as a recommendation for soil-structure interaction analysis in everyday design. Most recently, the European Building Code (Eurocode 7) has endeavored to include a variety of recommendations for SSI in design (Frank, 2006). Due to the lack of validation and verification of many SSI methods, most building codes have not mandated this type of analysis (Celebi and Crouse, 2001). A myriad of methods has been proposed to solve the soil-structure interaction problem. To simplify the problem, linear analysis techniques have been developed. One of the most commonly used approaches is a substructuring method which allows the problem to be analyzed in parts (Kintzer et al, 1982). The system is decomposed into superstructure and foundation sub-systems, thus allowing the dynamic response of each component to be investigated separately. The analysis of the foundation system can be reduced to the calculation of frequency dependent soil-springs and soil-dampers (more generally known as impedance functions) and driving forces from incident waves. The impedance functions include the effects of the mass, stiffness, and energy dissipation properties of the underlying soil as well as radiation damping. The kinematic interaction of the foundation with the incident waves is implemented in the form of a driving force vector. This method is also known as the Continuum Mechanics approach. Determining the soil impedance functions and the driving forces from incident waves is a rather complex process, especially when the geometric configuration of the foundation and the material properties of the soil medium are complicated. Several methods have been introduced. The Boundary Integral Equation Method has been used by Apsel (1979), Ohsaki (1973), and Tassoulas (1989), to name a few. Because of the difficult task of dealing 2 with the singularities of the Green Functions, many investigators pull the source points of the Green’s Functions away from the foundation surface, which results in the Indirect Boundary Integral Equation method. The basic premise of this method is to reduce a three-dimensional problem to a two-dimensional problem by considering only values on the boundary of the foundation surface. The Reciprocal Theorem (also known as Betti’s Law or Maxwell’s Reciprocal Theorem) is the basis for this formulation. Based on a theory of solid mechanics, the Reciprocal Theorem states that in a linear elastic solid, the work done by a set of forces acting through the displacements due to a second set of forces is equal to the work done by this second set of forces acting through the displacements of the first set of forces (Fung, 1965). With the ever-increasing amount and low-cost of computer power now available, the Finite Element modeling method is being used increasingly often in structural analysis. Unfortunately, there is a large difference between the finite element method applied to a structural model and that of a soil medium. The Finite Element Model of the structure can be highly refined because the material behavior of building elements is well understood. On the other hand, the standard Finite Element Method has difficulties modeling the semi-infinite space on which the structure rests. An artificial boundary is necessary to limit the total number of degrees of freedom. Unfortunately, this same boundary creates the problem of trapped wave energy. The outgoing waves reflect off the artificial boundary instead of being radiated through the large, underlying medium (Kintzer et al, 1982). This creates standing wave patterns, or resonance, which is exactly the opposite of what happens in a real soil medium. In the real, physical world, the large spatial dimensions of the soil medium allow the wave energy to damp out. This phenomena, known as radiation damping, is a major factor in reducing damage to the structure from powerful earthquake wave energy. 3 A popular area of research in earthquake engineering is nonlinear soil analysis in soil- structure interaction, which can be handled by the more ambitious Finite Element models. According to Duan (1999), the kinematic interaction between the embedded foundation and the incident waves has a greater effect on the structural response than the soil conditions. If this is correct, soil nonlinearity is then of secondary importance in the analysis. Nevertheless, a nonlinear analysis will help promote the understanding of the limitations of linear analyses. This thesis will focus on the linear realm of soil-structure interaction. A number of experimental investigations by Japanese researchers have been most helpful in determining the effect of SSI on an actual structure (Fujimori et al, 1992; Akino et al, 1996; Mizuhata et al, 1988; Iguchi et al, 1988; Urao et al, 1988; Watakabe et al, 1992; Imamura et al, 1992). These were performed in-situ, and are therefore, models which cannot be duplicated in a laboratory. This type of model most closely represents an actual building configuration, and therefore provides a standard by which analytical results can be compared. Urao et al (1988) used the Thin Layer Method (TLM) to compare their experimental results with an analytical model and were successful with some data correlation. Kurimoto et al (1992) developed a complex soil spring for irregularly embedded foundations by combining bottom and side soil springs. Parmlee and Kudder (1973) used a discrete parameter model solution to investigate embedment effects of tall buildings in a flexible foundation medium with a rigid embedded portion. Using a cone model, Takewaki et al (2003) developed a fast and practical method for evaluating SSI effects in embedded structures. This model is based on the assumption that the force transmitting mechanism of a foundation on the ground subjected to dynamic disturbances can be represented approximately by a cone “chopped” (truncated) by the foundation. Aviles and Perez-Rocha (1998) modeled soil- structure interaction effects by using free-field ground motion as the foundation input motion and then modified the fundamental period and associated damping of the structure. The 4 Finite Element Method was compared with soil springs in an embedment investigation by Seed and Idriss (1973) for massive embedded structures. Spryakos and Xu (2004) implemented a hybrid BEM (Boundary Element Model)-FEM model and performed several parametric studies within the soil-structure interaction problem. In the theoretical realm, Lysmer (1970) developed a lumped mass model for a layered, elastic half-space subjected to Rayleigh waves. Sandi (1974) solved the problem for an elastic medium with geometric and physical linearity. Gueguen et al (2002) studied how contamination of the free-field vibration occurs when several buildings are located within a small area. The formulation within this dissertation will follow the continuum mechanics approach. Three different applications of the Reciprocal Theorem will be made, each aiming at the same goal: to obtain impedance functions and driving force vectors. Although there have been previous efforts made using methods derived from the same principle, revisiting this approach in the 21st century may prove more fruitful for several reasons. First, the computational cost has diminished to nearly zero because of the availability of personal computers. Another significant advance in technology is the abundance of memory and disk storage, making next-to-impossible problems of the past now plausible. With the need for a soil-structure interaction analysis increased by recent damage assessments, validation of methods used in practical applications has become imperative. The first application of the Reciprocal Theorem will be a volume approach, which will be used to solve the displacement compatibility problem by adjusting the body forces placed inside the volume carved out by the boundary of an embedded foundation. The drawback to this method is similar to that of the Finite Element Method, namely, a three- dimensional formulation is used, requiring significant computational resources. However, this is a much smaller task than trying to model the infinite soil medium surrounding the foundation. This method also provides an alternative to the Boundary Integral Equation 5 Method by requiring only the displacement Green’s Function. In contrast, the Boundary Integral Equation Method requires both displacement and stress Green’s Functions. The second application of the Reciprocal Theorem is the direct formulation of the Boundary Integral Equation Method (BIEM). This formulation is mathematically elegant, but has many, well-known numerical problems. Rizzo (1989), Liu and Chen (1999), Burton and Miller (1971), and Pyl et al (2003) are among those who have tried to resolve the problem of fictitious eigenfrequencies of the interior medium which appear when the Boundary Integral Equation Method is used. This thesis will attempt to solve the problem with a new scheme. It will be shown in Chapter Four that the Traction and the Displacement Green’s Function matrices are Normal Matrices. This means that the two matrices have the same eigenvalues and eigenvectors. Also, while investigating the mathematical formulation of the Hilbert-Schmidt Method for solving an acoustic problem (Pao and Mow, 1973), it is clear that when using the BIEM, both the traction and the displacement Green’s Function matrices are controlled by a frequency-dependent function which causes both matrix determinants to approach zero at the same frequencies. It is interesting to note that these frequencies are the same as the resonant frequencies of the interior medium carved out by the boundary of the rigid foundation. In the mathematical formulation of the Hilbert-Schmidt method, this function, which appears with both matrices, can be eliminated by simply cancelling the function from both sides of the equation. But when a numerical method is used, this function cannot be removed. Experiments will be made to remove this function without disturbing the underlying mathematics. At the resonant frequencies of the interior medium, an eigenvector can be obtained for that zero eigenvalue and the zeroes in both matrices can be replaced by the ratio of the slopes of the functions, a direct application of l’Hospital’s Rule. 6 The third application of the Reciprocal Theorem is to use the integral representation as a basis to replace the artificial boundary of the Finite Element Method with a radiation boundary. In discrete formulations such as the Finite Difference and the Finite Element Methods, the number of nodes must be controlled due to limited computer memory. Since the calculation of the response of a particular node depends on those of nodes in its proximity, the truncation at the artificial boundary completely changes the nature of the model. It has been proposed by Chien (2007) to use the Reciprocal Theorem to generate values for one additional layer of nodes immediately outside the artificial boundary. The calculation of the values at the exterior nodes is based on the complete integration of displacements and tractions on the artificial boundary using the Green’s Functions of the problem. Because the integration involves all the nodes on the artificial boundary, it is not an efficient method. But considering the importance of the soil-structure interaction problem, this method is attractive because it combines the flexibility of the Finite Element Method with the unique, large soil medium representation made possible by continuum methods. In all three applications, the Green’s Function developed by Apsel (1979), Luco and Apsel (1983a, 1983b) will be used. Without this significant development, the present research would not be feasible. 7 Chapter Two Applications of the Reciprocal Theorem To solve the elastic wave propagation problem in a semi-infinite configuration, it is important to write an integral representation to relate the displacements, stresses and body forces near an embedded foundation. This integral representation theorem can be derived by using the Betti-Rayleigh relationship (Fung, 1965) to a body Ω as shown in Fig. 2.1. Figure 2.1 – Loading Conditions of a Solid Body. The bodyΩ, bounded by a surfaceS, is subjected to certain loading conditions, with ~ t being the traction vector applied at the surface and ~ f being the body force per unit volume applied within Ω. The displacement vector resulting from this loading condition is ~ u. Consider the action of another loading condition with ~ t and ~ f , they will create a new displacement vector ~ u . The reciprocal relationship is written as Z S ~ u T ~ tdS+ Z Ω ~ u T ~ fdV = Z S ~ u T ~ t dS + Z Ω ~ u T ~ f dV : (2:1) 8 Ω S : ~ f ^ e n ^ e n ~ u(~ r) ~ r p ~ t(~ r) ~ u(~ r) The above equation indicates that in a linear viscoelastic solid, the work done by a set of forces acting through the corresponding displacement produced by a second set of forces is equal to the work done by the second set of forces acting through the corresponding displacement produced by the first set of forces. The inner products of all four terms in Eq. (2.1) are scalars. This formulation is based on principles of physics and solid mechanics rather than a mathematical hypothesis. Since the inner product is a scalar, the inner-products on the right hand side can be reordered. The equation can then be written as Z S ~ u T ~ tdS+ Z Ω ~ u T ~ fdV = Z S ~ t T ~ udS+ Z Ω ~ f T ~ udV : (2:2) To make Eq. (2.2) useful for general applications, let the loading conditions represented by “*” be those generated by unit point loads in orthogonal directions. More specifically, consider 3 different cases where the body force ~ f is a point load applied at location~ r p = [x p ;y p ;z p ] T : CASE 1: Body force in thex-direction at point~ r p : ~ f x = h (~ r −~ r p ) ; 0 ; 0 i T (2:3a) CASE 2: Body force in they-direction at point~ r p : ~ f y = h 0 ;(~r−~r p ) ; 0 i T (2:3b) CASE 3: Body force in thez-direction at point~ r p : ~ f z = h 0 ; 0 ;(~r−~r p ) i T (2:3c) The resulting displacement vectors ~ u at~ r=[x;y;z] T for the above loading conditions are 9 CASE 1: ~ u x = h U xx (~ r p j~ r) ;U xy (~ r p j~ r) ;U xz (~ r p j~ r) i T (2:4a) CASE 2: ~ u y = h U yx (~ r p j~ r) ;U yy (~ r p j~ r) ;U yz (~ r p j~ r) i T (2:4b) CASE 3: ~ u z = h U zx (~ r p j~ r) ;U zy (~ r p j~ r) ;U zz (~ r p j~ r) i T (2:4c) in whichU ij is the displacement in ther j -direction caused by a point load in ther i -direction. The column vectors in Eq. (2.4) represent the displacement Green’s functions. The index notation forr i is defined asr 1 =x,r 2 =y andr 3 =z. Following similar steps, the traction vectors ~ t at~ r resulting from the same loading conditions are CASE 1: ~ t x = h T xx (~ r p j~ r) ;T xy (~ r p j~ r) ;T xz (~ r p j~ r) i T (2:5a) CASE 2: ~ t y = h T yx (~ r p j~ r) ;T yy (~ r p j~ r) ;T yz (~ r p j~ r) i T (2:5b) CASE 3: ~ t z = h T zx (~ r p j~ r) ;T zy (~ r p j~ r) ;T zz (~ r p j~ r) i T (2:5c) in whichT ij is the traction developed at~ r in ther j -direction by a point load at~ r p in the r i -direction. The traction vectors in Eq. (2.5) can be calculated by performing the matrix product of the Green’s Function stress tensor and the direction cosines of the outward normal vector ^ e n on surfaceS: 2 4 T x T y T z 3 5 = 2 4 xx xy xz yx yy yz zx zy zz 3 5 2 4 n x n y n z 3 5 (2:6) 10 in whichn x ,n y andn z are the direction cosines of ^ e n . By substituting the 3 sets of loading conditions into Eq. (2.2), a matrix equation of the following form results, Z S U(~ r p j~ r) ~ t(~ r)dS+ Z Ω U(~ r p j~ r) ~ f(~ r)dV = Z S T(~ r p j~ r) ~ u(~ r)dS+~ u(~ r p ); (2:7) in which U(~ r p j~ r) = 2 6 6 6 4 U xx (~ r p j~ r) U xy (~ r p j~ r) U xz (~ r p j~ r) U yx (~ r p j~ r) U yy (~ r p j~ r) U yz (~ r p j~ r) U zx (~ r p j~ r) U zy (~ r p j~ r) U zz (~ r p j~ r) 3 7 7 7 5 ; (2:8) and T(~ r p j~ r) = 2 6 6 6 4 T xx (~ r p j~ r) T xy (~ r p j~ r) T xz (~ r p j~ r) T yx (~ r p j~ r) T yy (~ r p j~ r) T yz (~ r p j~ r) T zx (~ r p j~ r) T zy (~ r p j~ r) T zz (~ r p j~ r) 3 7 7 7 5 ; (2:9) and the vector~ u(~ r p ) on the right-hand side of Eq. (2.7) is the result of the integral Z Ω 2 4 (~ r −~ r p )0 0 0 (~r−~r p )0 00 (~r−~r p ) 3 5 2 4 u 1 (~r) u 2 (~r) u 3 (~r) 3 5 dV : (2:10) The typical application of the integral representation of Eq. (2.7) is to generate the wave motion at the location~ r using the displacements and tractions at the surfaceS: ~ u(~ r p )= Z S U(~ r p j~ r) ~ t(~ r)dS − Z S T(~ r p j~ r) ~ u(~ r)dS : (2:11) This formulation is generally preferred because an integral evaluated over a surface areaS is one dimension smaller than an integral evaluated over a volumeΩ. In this proposal, three different methods will be explored, including one which employs the body force integral. 11 2.1 – Application One, A Body Force Formulation As an application of the integral representation in Eq. (2.7), it is proposed here to eliminate the existence of the surfaceS and formulate the embedded foundation problem using the volume integral of body forces only. The Reciprocal Theorem can therefore be written as Z Ω U(~ r p j~ r) ~ f(~ r)dV =~ u(~ r p ) ; (2:12) and the displacement at any point in the medium~ r p can be calculated as an integral of the displacement Green’s Function and the body forces distributed inside a volume Ω.Asan application to find the impedance of an embedded foundation, consider the schematic shown in Fig. 2.2. Let the volume Ω be that occupied by the foundation and ~ f(~ r) the distribution of unknown body forces inside the volume which would constrain the displacement in the volume to be rigid body motion, with 3 translational components, 2 rocking modes, and 1 torsional mode. The driving force vector induced by incident wave motion can also be obtained by determining the body force distribution which would hold the foundation mass fixed while subjected to the incoming wave motion. The numerical solution of Eq. (2.12) can be accomplished by subdividingΩ into small subregions and assuming a zeroth order approximation of a constant body force within that subregion. Prescribe the displacement~ u(~ r) to match the motion of a rigid foundation and Eq. (2.12) can be discretized into a matrix equation of the form ~ f =~u; (2:13) in which ~ f contains the values of the constant body forces inside the subregions and ~ u represents the prescribed displacements. The elements of the matrix [] are defined as the 12 Figure 2.2 – Body Force Distributions Inside an Embedded Foundation. integral of the displacement Green’s Function inside the subregionv j : ij = Z V j U(~ r i j~ r) dV : (2:14) The disadvantage of this method is that many subregions might be needed to obtain an accurate solution that satisfies the rigid body motion displacement constraints within Ω. Since the volume integral is three-dimensional, the number of subregions can become large very quickly. The advantage is that the formulation is simple and only the displacement Green’s Function is needed. The Boundary Integral Equation Method would also require the stress Green’s Function. After the body force vector ~ f is determined by the inversion of Eq. (2.13), the translational impedance can be obtained as K = Z Ω ~ fdV +! 2 Z Ω dV ; (2:15) 13 x y z V j V i in which is the mass density of the soil and the second term is the inertial force of the mass contained inΩ for an embedded foundation oscillating at the frequency!. Eq. (2.15) is used because the conventional definition of the impedance function does not include the mass of the rigid foundation. The application of this formulation appears to be straightforward and the convergence of the results should be excellent and stable. The only concern is that the second term of Eq. (2.15) may become very large at high frequencies, requiring an alternate formulation for that case. Also, the singularity of the dynamic Green’s Function needs to be treated with care. Fortunately, it is well known that the singularity of the dynamic Green’s Function is the same as that of the static Green’s Function. In fact, the singularity of the Green’s Function is the same in an infinite space as that in a semi-infinite space. It is proposed that the integration of the diagonal terms of Eq. (2.14) be performed in two parts: (i) Integrate the static Green’s Functions analytically since they are simple algebraic expressions. (ii) Integrate the difference of the dynamic and static Green’s Function numerically using Gaussian Quadrature. The numerical integral in part (ii) contains no singularity as it has been removed through subtraction. 2.2 – Application Two, A Direct Boundary Integral Equation Method This application assumes there are no distributed body forces in the wave propagation medium. Therefore, the Reciprocal Theorem as given in Eq. (2.7) can be simplified as shown earlier by Eq. (2.11) and becomes ~ u(~ r p )= Z S U(~ r p j~ r) ~ t(~ r)dS − Z S T(~ r p j~ r) ~ u(~ r)dS ; (2:16) 14 Figure 2.3 – A Schematic of the Direct Boundary Integral Equation Method. in which the volume integral over Ω has been eliminated. The wave propagation medium is that exterior toS and, therefore, the outward normal vector ^ e n is directed away from the exterior medium as shown in Fig. 2.1. The Green’s Function matrices, [U] and [T], have singularities if j~ r −~ r p j =0. This is the reason why many publications using this continuum approach have avoided putting the source point~ r p on the surfaceS. For several years, those methods were called Sources methods because the sources were placed inside the surfaceS and the displacement compatibility conditions at the boundary were satisfied using a least squares approximation. This method is now identified as the Indirect Boundary Integral Equation Method. If the source point ~ r p is moved onto S, then the integral over S would involve a point in whichR = j~ r −~ r p j=0. For the [U] matrix, the Green’s Function entries have singularities of the form: 1=R,x=R 2 ,y=R 2 andxy=R 3 . When cylindrical coordinates are used,x = Rcos andy = Rsin , and all singularities are of the form 1=R with some 15 ~ r p ~ t(~ r);~ u(~ r) x z y azimuthal dependencies. Therefore, all the singularities of the matrix[U] are integrable over S. The singularities of the matrix [T], however, are one order higher at 1=R 2 and are not integrable over the surfaceS. If the principal value integral is used for the term involving [T], then the limit when~ r p approachesS is the integral equation, 1 2 ~ u(~ r p )= Z S U(~ r p j~ r) ~ t(~ r)dS − Z S T(~ r p j~ r) ~ u(~ r)dS ; ~ r p onS; (2:17) in which the second integral on the right hand side is the principal value integral. Chapter Six will discuss a scheme which subtracts out the static singularities and performs those integrals in closed form and then integrates the difference of the two functions numerically. This procedure is only required for the subregion ofS where~ r p and~ r coincide. Unusual as it may seem, the numerical evaluation of the principal value integral of[T] is simple. Although the singularities are not integrable, they are odd functions. Therefore, the resulting numerical value of the integral is equal to that of the even part, or the nonsingular part, of the function. Using Gaussian Quadrature with an even number of sample points per axis of integration, i.e., an algorithm which does not require the center point, the point where the singularity resides, convergence can be obtained using as few as two points per integration direction. The numerical solution of Eq. (2.17) can be accomplished by subdividingS into small subregions and assuming, as a zeroth order approximation, the unknown traction to be constant within that subregion. Now prescribe the displacement~ u(~ r) to match the motion of a rigid foundation and Eq. (2.17) can be discretized into a matrix equation: U ~ t = 1 2 [I]+[T] ~u; (2:18) in which ~ t contains the values of the unknown tractions at the subregions of S and ~ u represents the prescribed displacements. The block elements of the matrices [U] and [T] 16 are defined as the integral of the displacement and traction Green’s Functions, respectively, over the subregion surfaceS j : U ij = Z S j [U(~ r i j~ r)]dS ; (2:19) and T ij = Z S j [T(~ r i j~ r)]dS : (2:20) After the numerical solution of Eq. (2.18) is obtained, the translational impedance functions can be computed as K = X j ~ t Z S j dS : (2:21) The calculation of rotational impedances is similar to the calculation of moments and will be covered in detail in Chapter Six. 2.3 – Application Three, A Radiation Boundary for Discrete Methods For discrete formulations such as the finite difference and the finite element methods, the number of nodes must be limited due to memory constraints. Thus, an artificial boundary is created to curb the size of the model. For many applications outside of the earth sciences, this is not a severe constraint. But for wave propagation models, it completely changes the nature and behavior of the model. This modeling misrepresentation is created by the fact that the calculation of values at a particular node depends on nodal values in its proximity. The truncation at the artificial boundary reflects the outgoing waves instead of allowing them to pass through, and thereby creates unwanted standing wave patterns. It has been proposed by Chien (2007) to use the Reciprocal Theorem to generate values for one additional layer of nodes immediately outside of the artificial boundary in order to eliminate the artificial boundary. The calculation of the values at the exterior nodes is based on the integration of the entire artificial boundary using the appropriate Green’s Functions. The integration of 17 the displacements and tractions provided by the discrete model at the artificial boundary is also necessary. Since the integration involves all the nodes on the artificial boundary, the banded matrix property of the discrete method will be destroyed, rendering this an inefficient method. But when the problem is important enough, efficiency is of secondary importance. As a tradeoff, the size of the model of the interior medium can be reduced significantly. This method is worth exploring for the embedded foundation problem because it combines the power of discrete and continuum methods. Figure 2.4 – Formulation of a Radiation Boundary. Shown in Fig. 2.4 is a schematic of the proposed radiation boundary for a discrete model. Although the illustration is two-dimensional, there is no limitation for this method. It can be expanded to handle three-dimensional models in a semi-infinite space. The nodes are identified by the subscripts E, A, S and I. The subset A contains the nodes on the artificial boundary. E contains the nodes exterior to the artificial boundary. S contains the interior nodes adjacent toA, andI represents all the remaining interior nodes, which constitute the most important portion of the analysis. The subset S was extracted from 18 Subset A Subset S Subset E Subset I subsetI for the purpose of calculating stresses and tractions at the artificial boundary. The displacement values atE will be determined from the displacements and tractions at the artificial boundary using the Green’s Functions of the problem and Eq. (2.11). The indices of the nodes are ordered starting with the nodes in subsetE are numbered first, followed by subsetsA,S andI. The formulation of the discrete model can be written in a matrix form as [Z]~ u = ~ f; (2:22) in which[Z] is a frequency dependent impedance matrix of the discrete model which includes stiffness, damping, and inertial effects of the soil condition near the embedded foundation. The forcing function, ~ f, is included to provide excitation from the interior of the model. The type of excitation can be concentrated or distributed loads or spherical waves generated as reflections of incident waves from irregular wave scattering boundaries. According to the pre-assigned nodal orders, the above matrix equation can be partitioned as follows: 2 6 6 6 6 6 6 6 4 [Z EE ][Z EA ] [0] [0] [Z AE ][Z AA ][Z AS ] [0] [0] [Z SA ][Z SS ][Z SI ] [0] [0] [Z IS ][Z II ] 3 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 4 ~ u E ~ u A ~ u S ~ u I 3 7 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 4 ~ 0 ~ 0 ~ 0 ~ F I 3 7 7 7 7 7 7 7 5 ; (2:23) in which ~ F I is the forcing function applied at the interior part of the model. It is assumed that no forces are applied at or beyond the artificial boundary. For a model with a fixed artificial boundary, Eq. (2.23) would reduce to 2 4 [Z SS ][Z SI ] [Z IS ][Z II ] 3 5 2 4 ~ u S ~ u I 3 5 = 2 4 ~ 0 ~ F I 3 5 : (2:24) 19 This is not a useful assumption for harmonic wave problems as it does not possess any characteristics of a wave propagation model. Unwanted reflected waves from the artificial boundary exist regardless of how far the boundary is placed away from the essential locations of the model. The solution of Eq. (2.24) has many resonant frequencies because the reflected waves from the artificial boundary interfere constructively with the outgoing waves. To remove the problems caused by an artificial boundary, it is proposed to calculate the displacement,~ u E , as defined in Eq. (2.23) using the Green’s Functions of the problem and a representation theorem of the form, ~ u(~ r E )= Z S h U(~ r E j~ r A ) i ~ t A (~ r A ) − h T(~ r E j~ r A ) i ~ u A (~ r A ) dS ; (2:25) in which~ u A and ~ t A are the displacement and traction vectors at the artificial boundary. The integration is to be performed over a closed boundary,S, in this case, the artificial boundary. The Green’s function matrices, [U(~ r E j~ r A )] and [T(~ r E j~ r A )], contain the displacement and traction, respectively, observed at~ r A , generated by a point source located at~ r E . The above formulation allows the displacements,~ u E , to be calculated without information from nodes outside of the domain ofE, therefore, the waves can propagate outward based on the motion of the interior nodes. For numerical calculation, approximate the integrals in Eq. (2.25) by numerical quadratures of the form, Z S f(~ r A )g(~ r A )dS = X j w j f(~ r j )g(~ r j ) ; (2:26) in which w i are weighting factors and ~ r j are the sample points selected for numerical integration. With this discretization scheme in place, the right hand side of Eq. (2.26) can be written as a matrix product and it can be replaced by the matrix equation ~ u E = −[T zz ]~ u A +[U zz ] ~ t A ; (2:27) 20 in which the matrices[T zz ] and[U zz ] are not necessarily dimensionally square because the number of nodes in subsetE is not likely to be equal to that of subsetA. The vector~ u A in Eq. (2.27) appears as an unknown in Eq. (2.23), therefore, the first term on the right hand side of Eq. (2.27) can be implemented simply by substitution. The traction vector ~ t A , on the other hand, must be computed from the stress distributions near the artificial boundary and is therefore dependent on nodes from SubsetS. In matrix form, the traction vector ~ t A can be expressed in terms of~ u A and~ u S as ~ t A =[D A ]~ u A +[D S ]~ u S ; (2:28) in which[D A ] and[D S ] are sparse matrices containing mostly rows with one nonzero value except for those associated with the corner nodes on the artificial boundary. Eq. (2.27) can now be written as ~ u E =(−[T zz ]+[U zz ][D A ])~ u A +[U zz ][D S ]~ u S ; (2:29) with the aid of Eq. (2.28). The matrix products,[U zz ][D A ] and[U zz ][D S ], are best performed by the strategic use of indices by combining scalar multiples of the rows of [U zz ] since the [D] matrices are not far from being identity matrices. With Eq. (2.29) available, Eq. (2.23) can now be condensed because~ u E can be replaced by~ u A and~ u S . Eliminating the first row of the partitioned matrix equation, Eq. (2.23) can be re-ordered as 2 6 6 6 4 [Z AA ][Z AS ] [0] [Z SA ][Z SS ][Z SI ] [0] [Z IS ][Z II ] 3 7 7 7 5 2 6 6 6 4 ~ u A ~ u S ~ u I 3 7 7 7 5 = 2 6 6 6 4 −[Z AE ]~ u E ~ 0 ~ F I 3 7 7 7 5 : (2:30) 21 The substitution of Eq. (2.29) for~ u E on the right hand side of Eq. (2.30) and subsequent rearrangement of unknown variables yields a complex matrix equation for the finite difference model as 2 6 6 6 4 [Z AA +Z AE (−T zz +U zz D A )] [Z AS +Z AE U zz D S ] [0] [Z SA ][Z SS ][Z SI ] [0] [Z IS ][Z II ] 3 7 7 7 5 2 6 6 6 4 ~ u A ~ u S ~ u I 3 7 7 7 5 = 2 6 6 6 4 ~ 0 ~ 0 ~ F I 3 7 7 7 5 : (2:31) The above equation will provide the mechanism for outgoing waves to propagate outward without reflection from the radiation boundary. The important difference between Eq. (2.31) and Eq. (2.24) is that the new formulation is complex and, therefore, includes radiation damping. Resonant behavior caused by the artificial boundary will be eliminated. 22 Chapter Three A Distributed Force Formulation As an application of the integral representation in Eq. (2.7), the need for surface S is eliminated and and the rigid embedded foundation problem is formulated using only a volume integral of distributed forces. The Reciprocal Theorem can be simplified to Z Ω U(~ r p j~ r) ~ f(~ r)dV =~ u(~ r p ) : (3:1) In the above expression, the displacement vector at any point in the medium,~ r p , can be calculated as an integral of the displacement Green’s Function,U, and the distributed forces inside volume Ω. The major advantage of this formulation is that only the displacement Green’s Functions are needed and they are more easily obtained than the traction Green’s Functions because tractions are sensitive to the detailed geometry of the scatterer’s surface. The disadvantage of this formulation is the increased computational effort due to the additional dimension in the volume integral over a surface integral. However, the “1/R” type of singularity in the displacement Green’s Functions is efficiently removed by integration over 3 independent variables. To find the impedance of a rigid embedded foundation, let the volume Ω of Eq. (3.1) be that occupied by an embedded foundation and ~ f(~ r) be the vector of unknown distributed forces inside the foundation volume which would constrain the displacement within the volume to move as a rigid body. The six modes of rigid body motion include three translational components, two rocking components and one torsional component. In a separate formulation to obtain the driving force vector induced by incident wave motion, 23 Eq. (3.1) can also be used to determine the body force distribution that would hold the foundation mass fixed while being subjected to the incoming wave motion. The body force discussed here is different from what is typically implied in other applications. It is most commonly associated with gravitational or electromagnetic forces on the material within a volume. In the present application, it is used as a vehicle to prescribe rigid body motion for the foundation volume, to simulate a rigid embedded foundation. In the remainder of this chapter, the terms body forces and distributed forces will be used interchangeably. Figure 3.1 – Distributed Forces Inside an Embedded Foundation. A numerical solution of Eq. (3.1) can be accomplished by first subdividing Ω intoN subregions,V i ;i=1;2;:::;N. Then, Eq. (3.1) can be written as N X j=1 Z V j [U(~ r p j~ r)] ~ f(~ r)dV =~ u(~ r p ) : (3:2) 24 x y z V j V i The next step is to approximate ~ f(~ r) as a constant vector, ~ f j , insideV j . Eq. (3.2) becomes N X j=1 Z V j [U(~ r p j~ r)]dV ! ~ f j =~ u(~ r p ) : (3:3) Assign the source point of the Green’s function, ~ r p ,to be~r i , the centroid of one of the subregions,V i , then Eq. (3.3) can be expressed as a 3 3 matrix equation N X j=1 [ ij ] ~ f j =~ u i ;i=1;2;:::;N; (3:4) in which [ ij ]= Z V j [U(~ r p j~ r)]dV ; (3:5) will be referred to as the influence matrix which relates the displacement vector at~ r i , the centroid of subregionV i , to the uniform body force ~ f j in subregionV j . Now replace the summation in Eq. (3.4) by a matrix representation of order 3N 3N expressed as ~ f =~u; (3:6) in which ~ f and~ u contains the body force vector and the displacement vectors at the centroids of allN subregions, respectively. To adapt Eq. (3.6) to the calculation of the impedance matrix of an embedded foundation, let the displacement vector, ~ u, be defined by rigid body motion as defined for each subregionV i as ~ u i =[ i ] ~ U 0 ; (3:7) in which ~ U 0 is a generalized displacement describing the six-degree-of-freedom foundation motion as ~ U 0 = 2 6 6 6 6 6 4 U x U y U z b x b y b z 3 7 7 7 7 7 5 ; (3:8) 25 and [ i ] is the rigid-body kinematic matrix defined as [ i ]= 2 4 100 0 z i =b −y i =b 010 −z i =b 0 x i =b 001 y i =b −x i =b 0 3 5 : (3:9) The parameterb in Eq. (3.8) and Eq. (3.9) has a unit of length related to the lateral dimension of the foundation so that all entries in the generalized displacement,~ u, will have units of length. Figure 3.2 – The Six Degrees of Freedom for a Rigid Embedded Foundation. The solution of the matrix equation, Eq. (3.7), can be obtained numerically as f r = −1 ~ U 0 ; (3:10) 26 x z y U x U y U z x (x i ;y i ;z i ) y z in which [f r ] is a 3N 6 matrix containing the 6 body force distributions which satisfy the displacement compatibility conditions specified by the 3N 6 matrix []. The matrix [ i ] in Eq. (3.9) is a subset of[] in Eq. (3.10); the latter contains all the centroids of theN subregions. The final step of the formulation is to accumulate the body forces into three forces and three moments using the transpose of the rigid body kinematic matrix as ~ F 0 = T f r = T −1 ~ U 0 ; (3:11) in which ~ F 0 = 2 6 6 6 6 6 4 F x F y F z M x =b M y =b M z =b 3 7 7 7 7 7 5 ; (3:12) is the generalized force, including three forces and three moments, normalized to units of force. On the right hand side of Eq. (3.11), the matrix product, T −1 , is normally the definition of the66 impedance matrix, however, the body forces f r have the additional task of moving the mass within the volume. The impedance matrix for a massless embedded foundation must now be defined as K = T −1 +! 2 m ; (3:13) in which m contains the mass of the soil within each subregionV i . The inertia term in Eq. (3.13), in detailed form is 27 T m =(3:14) N X i=1 m i 2 6 6 6 6 6 4 10 0 0 z i =b −y i =b 01 0 −z i =b 0 x i =b 00 1 y i =b −x i =b 0 0 −z i =b y i =b (y 2 i +z 2 i )=b 2 −x i y i =b 2 −x i z i =b 2 z i =b 0 −x i =b −x i y i =b 2 (x 2 i +z 2 i )=b 2 −y i z i =b 2 −y i =b x i =b 0 −x i z i =b 2 −y i z i =b 2 (x 2 i +y 2 i )=b 2 3 7 7 7 7 7 5 ; containing mass and normalized moment of inertia terms. 3.1 – Numerical Procedure for the Integration of Singularities The integral in Eq. (3.5) can be evaluated accurately by numerical quadratures except for the diagonal elements ii because the source point,~ r p , and the observation point,~ r, are located in the same subregionV i . The singularities of the displacement Green’s Functions are of the type 1=R. They are easily integrable if the integrals are evaluated analytically, but it is not possible to evaluate them numerically because the denominator becomes zero when the source and observation points are at the same location. It is comforting to recognize that for this difficult numerical problem there is a relatively simple solution. It is well known in the theory of solid mechanics that the form of the singularities for all Green’s Functions are the same, regardless if the problem is dynamic or static, an infinite space or a semi-infinite space configuration. It is proposed to evaluate ii of Eq. (3.5) in two parts as [ ii ]= Z V i [U(~ r p j~ r)]dV = Z V i [U(~ r p j~ r)] −[U (~ r p j~ r)] dV + Z V i [U (~ r p j~ r)]dV ; (3:15) 28 in which [U (~ r p j~ r)] is any appropriate Green’s Function matrix which can be expressed in closed form. For this particular application, the static point load solution (Love, 1906) in an infinite space will be used. With the procedure introduced in Eq. (3.15), the first integral can be evaluated numerically since it is now void of singularities and the second integral can be evaluated analytically. The static point load solution by Love (1906) can be rearranged for this application in the form, U (~ r p j~ r) = (1+γ 2 ) 8R 2 4 100 010 001 3 5 (3:16) + (1 −γ 2 ) 8R 3 2 4 (x −x p ) 2 (x −x p )(y −y p )(x−x p )(z −z p ) (x −x p )(y −y p )(y−y p ) 2 (y−y p )(z −z p ) (x −x p )(z −z p )(y−y p )(z −z p )(z−z p ) 2 3 5 ; in whichR = p (x −x p ) 2 +(y−y p ) 2 +(z−z p ) 2 ,γ = = is the ratio of the shear wave velocity to the compressional wave velocity. The constant is the shear modulus. The integral ofU as shown in Eq. (3.16) can be performed analytically as Z V i [U (~ r p j~ r)]dV = (1+γ 2 ) 8 2 4 S 0 00 0 S 0 0 00 S 0 3 5 + (1 −γ 2 ) 8 2 4 S xx S xy S xz S xy S yy S yz S xz S yz S zz 3 5 ; (3:17) in which the expressions forS, in closed form, are given in Appendix C. 3.2 – Application of Green’s Functions Obtained in Cylindrical Coordinates Most Green’s Functions derived in classical wave propagation research are given in cylindrical coordinates because the configuration of a point load adapts well to that system. For time harmonic Green’s Functions, the solutions are expressed in Hankel Transforms and numerical evaluation is tedious. As an example, the Green’s Function computer program created by Luco and Apsel (1979a, 1979b) required a number of years to develop. 29 For the present method, only the displacement Green’s Function are needed. The applicable functions are f rz and f zz for the vertical point load and f rr , f r and f zr for the horizontal point load. The functions,f, are all nonsingular. The Green’s Functions are obtained by dividing by the factor 1=(r). Consider now the coordinate system illustrated in Fig. 3.3. Note that the polar coordinate system is defined with the z-axis pointing downward, which is the typical convention for classical geophysical problems. The origin of the polar coordinate system is defined at the source point, ~ r p , therefore, the observation point, ~ r, is located at the coordinates, (r; ;z), in which r = j~ r −~ r p j= q (x −x p ) 2 +(y−y p ) 2 (3:18) and =arg(~ r −~ r p )=tan −1 y −y p x −x p : (3:19) Thez-dependency of the Green’s Functions is included in the functionsf. Using the polar coordinate system, all horizontal point forces can be represented by P r ( 0 ) because the reference angle can be varied to match any orientation. Therefore, the general displacement-force relationship can be written in a form of a3 2 matrix as 8 < : u r (r; ) u (r; ) u z (r; ) 9 = ; = 1 r 2 4 f rr cos( − 0 ) f rz f r sin( − 0 )0 f zr cos( − 0 ) f zz 3 5 P r ( 0 ) P z : (3:20) To obtain a displacement-force relationship in the Cartesian coordinate system, the following mapping of the forces and the displacements may be used: 8 < : P x (~ r p ) P y (~ r p ) P z (~ r p ) 9 = ; = 8 < : P r (0) P r (=2) P z 9 = ; ; (3:21) 30 Figure 3.3 – The Definition of the Cylindrical Coordinate System. and 8 < : u x (~ r) u y (~ r) u z (~ r) 9 = ; = 2 4 cos −sin 0 sin cos 0 00 1 3 5 8 < : u r (r; ) u (r; ) u z (r; ) 9 = ; : (3:22) The first, Eq. (3.21), can be used with Eq. (3.20) to relate the displacements in polar coordinates to the forces in Cartesian Coordinates as 8 < : u r (r; ) u (r; ) u z (r; ) 9 = ; = 1 r 2 4 f rr cos f rr sin f rz f r sin −f r cos 0 f zr cos f zr sin f zz 3 5 8 < : P x (~ r p ) P y (~ r p ) P z (~ r p ) 9 = ; ; (3:23) Now apply the transformation (3.22) to both sides of Eq. (3.23). The result is a displacement- force relationship in Cartesian coordinates written as 8 < : u x (~ r) u y (~ r) u z (~ r) 9 = ; = 1 r G 8 < : P x (~ r p ) P y (~ r p ) P z (~ r p ) 9 = ; ; (3:24) 31 r x y z 0 P r P z P y P x in which G = 2 6 4 f rr cos 2 −f r sin 2 (f rr +f r )sin cos f rz cos (f rr +f r )sin cos f rr sin 2 −f r cos 2 f rz sin f zr cos −f zr sin f zz 3 7 5 : (3:25) Using Cartesian coordinates, the factors, sin and cos , can be evaluated simply as sin = y −y p r ; (3:26) and cos = x −x p r : (3:27) If the elements of matrix [G] are defined with subscripts as G = 2 4 G xx G xy G xz G yx G yy G yz G zx G zy G zz 3 5 ; (3:26) then the matrix [U(~ r p j~ r)] in Eq. (3.5) can be formed as U(~ r p j~ r) = 1 r 2 4 G xx G yx G zx G xy G yy G zy G xz G yz G zz 3 5 : (3:27) A column in[G] represents the displacements caused by a point load in a particular direction, following the order x, y, z, respectively. In the matrix [U(~ r p j~ r)], however, the same displacement components are stored as a row. The reason for this transposition can be observed in Eq. (2.2) in which the Green’s Functions entries, i.e., ~ u T are written as a row vector in preparation for an inner product. 3.3 – Convergence Rates of the Distributed Force Method To test the rate of convergence for the method proposed in the previous section, the static Green’s Functions provided by Mindlin (1936) will be used. After some 32 modifications, using the parameterγ ==, the ratio between the shear wave velocity and the compressional wave vecocity, the five static displacement Green’s Functions written in cylindrical coordinates can be expressed as f rr = r(1 −γ 2 ) 8 1+γ 2 1−γ 2 1 R 1 + r 2 R 3 2 + 1 R 2 + r 2 R 3 1 + 2cz R 3 2 1 − 3r 2 R 2 2 + 2γ 2 (1 −γ 4 )(R 2 +z +c) 1 − r 2 R 2 (R 2 +z +c) ; (3:28) f r = r 3 (1 −γ 2 ) 8 1 R 3 1 + 1+γ 2 (1 −γ 2 )R 3 2 − 6cz R 5 2 − 2γ 2 (1 −γ 4 )R 2 (R 2 +z +c) 2 −f rr ; (3:29) f zr = r 2 (1 −γ 2 ) 8 z −c R 3 1 + 1+γ 2 1−γ 2 (z−c) R 3 2 − 6cz(z +c) R 5 2 + 2γ 2 (1 −γ 4 )R 2 (R 2 +z +c) ; (3:30) f rz = r 2 (1 −γ 2 ) 8 z −c R 3 1 + 1+γ 4 (1 −γ 2 )R 3 2 − 2γ 2 (1 −γ 4 ) 1 R 2 (R 2 +z +c) + 6cz(z +c) r 5 2 ; (3:31) and f zz = r(1 −γ 2 ) 8 z −c R 3 1 + 1+γ 4 (1 −γ 2 ) 2 R 2 + (z −c) 2 R 3 1 + 1+γ 2 1−γ 2 (z+c) 2 −2cz r 3 2 + 6cz(z +c) 2 r 5 2 : (3:32) The above equations can now be used along with Eq. (3.25) to create the Green’s Function matrix necessary for the numerical solution. 3.3.1 – Solution with a Point-Displacement Compatibilty Condition To develop a model for an embedded foundation, consider the special case of a square foundation of width2b and depthh. The dimension2b has been used as the width in many previous publications because early results for embedded foundation were obtained for a 33 cylindrical foundation. To make an easy comparison of the two results, it is convenient to use the half-width of the rectangular foundation. The Poisson’s Ratio,, is set to 1/3 for the results in section 3.3; the corresponding ratio of body wave velocities,γ, is therefore ==1=2. Figure 3.4 – A Schematic for a Rigid Embedded Foundation. Cubic subregions are used for the numerical solution. The first task is to prescribe the displacement compatibility conditions at the centroid of each subregion. In subsection 3.3.2, the displacement compatibility will be prescribed with the average displacement throughout the subregion. Since the subregions are cubical in shape, the available combinations ofh=b are not totally flexible. However, the simplest subregion geometry is used so the rate of convergence can be tested without involving too many parameters. Theh=b ratios included in this section are 0.25, 0.50, 0.75, 1.00, 1.50, and 2.00. The number of subregions used is the product ofN x ,N y andN z , the number of subregions in thex,y andz direction, respectively. For example, ifh=b is 0.25, the possible combinations for that model are(N x ;N y ;N z ) equal to 34 h 2b (8,8,1), (16,16,2), (24,24,3) or (32,32,4). There are no intermediate combinations possible unless non-cubical subregions are used. There are twoh=b ratios that are amenable to a more refined grid of combinations,h=b equal to 1 and 2. These will be used as the basis for a more detailed convergence test. Shown in Fig. 3.5 are the convergence tests for embedded ratios h=b equal to 0.25, 0.50 and 1.50. These are the more restrictive models because of their aspect ratios. For the cases withh=b equal to 0.25 and 1.50, only 4 data points are available for each graph. The horizontal axis of each graph in the figure is the value ofN x . In parts (a), (b) and (c), the five impedance functions: K hh , K vv ,K mm , K tt , K hm , are plotted. They are dimensionless values defined as K hh = F h bU h ;K vv = F z bU z ;K mm = M h b 3 h ; K tt = M z b 3 z ;K hm = F h b 2 h ; in which is a reference shear modulus, h is a horizontal direction, either x or y, v is the vertical or z direction, m is a moment about one of the horizontal axes, t is the moment about the z-axis, or the torsional component. The convergence trend exists for all 5 components, but it is clear that the translational components converge faster than the rotational components. This is to be expected. The radial distance from the center of rotation amplifies the inaccuracy of the higher stresses near the surface of the foundation. The stress distribution near the center of the rigid foundation contributes very little to the rotational components. The vertical scales were expanded to show the variations clearly, but the last two data points on each graph have a percentage difference of less than 5%. As mentioned previously, the number of subdivision combinations possible using a cubical subregion is the largest for the special cases ofh=b equal to 1 and 2. Forh=b=1, 35 Figure 3.5 – Performance of the Distributed Force Method for Various Embeddment Ratiosh=b. 36 8 162432 6 7 8 8 162432 6 12 18 8 162432 1 2 4 8 12 16 20 24 7 8 9 4 8 12 16 20 24 7 9 11 13 15 17 19 21 4 8 12 16 20 24 2 3 4 4 8 12 16 9 10 11 12 13 4 8 12 16 25 30 35 40 4 8 12 16 10 11 12 13 (a) h/b=0.25 (b) h/b=0.25 (c) h/b=0.25 (d) h/b=0.50 (e) h/b=0.50 (f) h/b=0.50 (g) h/b=1.50 (h) h/b=1.50 (i) h/b=1.50 K hh K vv K mm K tt K hm K hh K vv K mm K tt K hm K hh K vv K mm K tt K hm the combinations can be (N x ;N y ;N z )=(N x ;N x ;N x =2) for any even numberN x .For h=b=2, the combinations can be (N x ;N y ;N z )=(N x ;N x ;N x ) for all positive integer values ofN x . Fig. 3.6 shown the results forh=b=1 in parts (a) and (b) forN x values from 2 to 20 in increments of 2. The results forh=b=2 are shown in parts (c) and (d) of the same figure asN x varies between 1 to 16 in increments of 1. The large number of cases analyzed shows smooth convergence trends for both ratios and, therefore, demonstrates the stability of this method. The large number of subregions can be taxing on a computer. Some of the largest cases took over 10 hours of CPU time on a medium speed personal computer. But this is an academic exercise to provide validation of the method. Practical engineering applications would not require such precision because the determination of the soil properties is not an exact science. 3.3.2 – Solution with Average-Displacement Compatibilty Condition The results in section 3.3.1 were obtained by prescribing the rigid body displacement at the centroid of each subregion. The displacements away from the center points, however, may deviate from the rigid body condition. When this method is employed in the frequency domain, and when the excitation frequency is high, the deviation can be even larger. For numerical solutions with a smaller number of prescribed subregions, the compatibility conditions throughout the volume Ω are not strictly enforced. One improvement that has been made to other rigid foundation analyses is to apply displacement compatibility to the average displacement within the volume of a subregion. This would ensure that the displacement within the subregion, on the average, is equal to the prescribed condition. 37 Figure 3.6 – Convergence of Impedance Functions forh=b=1 andh=b=2. 38 2468 101214161820 4 6 8 10 12 2468 101214161820 8 10 12 14 16 18 20 22 24 26 28 30 1 4 7 10 13 16 6 8 10 12 14 16 18 20 1 4 7 10 13 16 0 10 20 30 40 50 60 70 (a) h/b=1 (b) h/b=1 (c) h/b=2 (d) h/b=2 K hh K vv K hm K mm K tt K hh K vv K hm K mm K tt 3.7 – Convergence Rate of Average-Displacement Compatibilty verus Point-Displacement Compatibility. 39 1 4 7 10 13 16 10 11 12 13 14 15 1 4 7 10 13 16 8 9 10 11 12 13 1 4 7 10 13 16 0 10 20 30 40 50 60 70 1 4 7 10 13 16 0 10 20 30 40 50 60 (a) K hh (b) K vv (c) K mm (d) K tt and K hm Center-Point Average K hm K tt Fig. 3.7 shows the rate of convergence for the two compatibility conditions. The results for the average-displacement scheme are plotted with a dashed line while those for the point- displacement scheme are plotted with a solid line. Since this method underestimate the exact theoretical solution, it is clear that the average-displacement scheme converges consistenly faster than the point-displacement result. Observing the plotted results suggests that fewer subregions can be used if the average-displacement compatibility scheme is applied. This will significantly reduce the matrix dimensions and, therefore, the amount of required computer memory when the matrix is inverted. The overall computer analysis time may not decrease, however, because integration must now be performed for the distributed forces as well as for the displacements, six dimensions of integration. But an alternative formulation always increases the understanding of the numerical algorithm, so the exercise is not without value. 3.3.3 – Convergence of the Dynamic Solution With the experience gained in the previous two subsections, the average displacement compatibility condition will be the condition of choice throughout the remainder of Chapter Three. For the higher frequencies, there can be several wavelengths within the model’s boundaries, therefore, an average displacement approach would guarantee better results. To present the dynamic time-harmonic results, it is convenient to define a dimensionless frequency, a 0 = !b=, in which ! is the circular frequency, b is the half width of the foundation and is the reference shear wave velocity. Since the width of the foundation is 2b, the value ofa 0 = represents a condition when the shear wave wavelength is exactly equal to the width of the foundation. The dynamic results will be presented for the range ofa 0 from 0 to 6 in the remaining figures of this chapter. 40 Figure 3.8 – Convergence Trend of Average-Displacement Method. 41 5 10 15 20 25 5 10 15 10 15 20 25 0 2 4 6 8 10 12 0123456 0 2 4 6 8 10 12 a 0 0123456 3 4 5 6 7 a 0 k hh c hh k mm c mm k mh c mh h/b=0.75 ν=1/3 32x32x12 24x24x9 16x16x6 8x8x3 From the practical point of view, the range ofa 0 described in the previous paragraph can be converted to real world unit in the following manner. For example, if the shear wave velocity of a relatively firm half space is 500 m/sec, and if the total width of a square foundation is 50 m, then a 0 =2 represents a frequency of 20 cycles/sec, or a circular frequency of 40 rad/sec. The reference point for the impedance functions is the origin of the coordinate system, defined at the center and at the top of the embedded square foundation. This will be the convention used, except for the results for a cylindrical foundation in Section 3.4 because the results will be compared to those provided by Apsel (1979) and his results are referred to the bottom of the cylindrical foundation. Shown in Fig. 3.8 is the horizontal, rocking and coupled horizontal-rocking components of the impedance function for a square foundation embedded in a semi-infinite space with an embedment ratio of h=b=0:75. The Poisson’s Ratio of the viscoelastic medium is again 1/3; lightly damped values of the medium are measure usingQ factors,Q = 100 for compressional waves, andQ =50 for shear waves. The top row of two figures in Fig. 3.8 shows the convergence of results for the horizontal impedance,K hh . The left figure shows the soil spring value as defined as the real part of the impedance, i.e., k hh = Re(K hh ); whereas the right figure shows the soil damper defined as the imaginary part of the impedance divided by the dimensionless frequency, i.e., c hh = Im(K hh )=a 0 . Four different subdivisions of the square foundation with embedment ratio ofh=b=0:75 is considered. The smallest model has number of subregions in thex,y andz direction as 8, 8, and 3, respectively. The more precise models have factors of 2, 3 and 4 times larger in each dimension and, therefore, 8, 27 and 64 times as many subregions in the foundation volume. Since the computer time for solving a set of simultaneous equations is of the orderN 3 , the most precise model would take more than a quarter million times longer than the simplest case. The computer time devoted for the most precise case is about 42 4 hours per frequency on a personal computer and the least precise case is about 10 seconds per frequency. Noting that the scales of the plots are magnified to show the differences, the rate of convergence is excellent except for the least precise case. The computer time, although quite long, is not unreasonable because there are many hours of idle time available on personal computers at the present time. It is not unreasonable to spend several hours of computer time to obtain good results for a seismic analysis. More evidence for the convergent trends are illustrated in the Rocking and Coupling impedance functions in the four lower figures. Uniformly, the imaginary parts converge with a higher degree of accuracy. But the real parts have some struggles in the higher frequency end. It is important to point out that at high frequencies, the imaginary parts are an order of magnitude larger than the real part and therefore the real parts lose one digit of accuracy due to the domination by the imaginary parts. It is also important to note that the less accurate results for the soil springs will not affect the accuracy of a soil-structure interaction analysis at high frequencies because the inertia and damping terms would dominate. This effect can be illustrated quite clearly by the forces of a single degree of freedom oscillator in the frequency domain, −! 2 m+i!c+k, where the inertia term is largest at high frequencies. The convergence test was conducted using uniform size cubical subregions for integration, therefore, it is not the most efficient way of attaining accuracy quickly. Too many small subregions were spent deep into the interior of embedded foundation where the smaller stresses do not contribute significantly to the final results. Future work with variable size subregions could enhance this method by placing smaller subregions near the foundation’s embedded surfaces and larger subregions in the interior. It is clear that this would be especially effective for rotational impedances as the near surface stresses, with larger moment arms, contribute more strongly to the moment calculations. 43 Over 70% of the computer time is devoted to the solution of the complex simultaneous equations after the matrix is assembled. It is hopeful that recent development of multiple core processors will decrease the computational time significantly. Research has shown that a Dual-Core or Quad-Core processor, with the aide of an optimized BLAS (Basic Linear Algegra System), can increase the speed of calculation by a factor of 2 to 5. At the time of this writing, the marketing of single core processors is effectively discontinued. Figure 3.9 – An Approximate Grid for a Cylindrical Foundation. 3.4 – Comparison with Published Results There are not many published impedances for embedded foundation to compared. One set of data which is ideal for comparison is the impedance functions for embedded circular cylindrical foundation prepared by Apsel (1979). Those results were calculated using the same set of Green’s Functions but using an indirect boundary integral equation method. 44 The present method is formulated over the volume of the foundation model while Apsel’s formulation placed sources inside the foundation surface and impose the displacement compatibility conditions at the surface approximately by minimizing the energy of a scalar potential. As mentioned in the previous section, Apsel’s results have the reference point set at the bottom of the cylindrical foundation, therefore, all results in Section 3.4 will be transformed to the bottom of the foundation. Figure 3.10 – Vertical and Torsional Impedances for Various Poisson’s Ratios. 45 -5 0 5 10 8 9 10 11 12 13 14 15 16 0123456 10 15 20 a 0 0123456 0 2 4 6 8 a 0 k vv c vv k tt c tt ν=0.25 ν=0.33 ν=0.45 Using cubical subvolumes, a circular cylinder can only be approximated using the present method. But a large number of subregions, as many as 12,992 for the embedment ratioh=b=1, the approximation of a circle is excellent as shown in Fig. 3.9. Since there are three component of displacements in each subregion, the number of degrees of freedom is nearly 40,000. A large computation effort was performed to ensure the results are as accurate as possible using a personal computer. Apsel’s impedance results were computed for five components: horizontal, vertical, rocking, torsion and a coupling term between horizontal and rocking. The publised results were obtained for the Poisson’s Ratio of 0.25. At this point in the analysis, it is important to point out how the Poisson’s Ratio,, does have a large effect on the impedances. This parameter, , affects the ratio of shear wave velocity to compressional wave velocity,γ = =, therefore, it is a critical parameter to characterize a soil medium. The results to be used for this comparison use the embedment ratio ofh=b=0:75. As shown in the lower figures of Fig. 3.10, the torsional impedance is not affected at all by the Poisson’s Ratio because only shear waves are generated by torsional motion of a cylindrical foundation. The rough edges of the approximated circle as shown in Fig. 3.9, however, may have caused some small inaccuracies. In the upper figures of Fig. 3.9, the vertical impedance is shown. For an embedded foundation, the vertical motion generates compressional waves from the bottom of the foundation and shear waves from its vertical side boundaries. Therefore,K vv is highly dependent of. For the higher value of =0:45, which represents soil media with a large moisture content, the soil spring actually becomes negative in the high frequency range. The examination of the impedance functions and their dependence on Poisson’s Ratio continues in Fig. 3.11, where the horizontal, rocking and coupling terms are shown. Although these impedances are frequency dependent, there appears a surprising jump on the horizontal components at a 0 1:9 when =0:25. Knowing some other methods 46 Figure 3.11 – Horizontal, Rocking and Coupling Impedances for Various Poisson’s Ratios. 47 0 2 4 6 8 10 12 0 5 10 15 0 5 10 15 0 1 2 3 4 5 6 0123456 -2 0 2 4 a 0 0123456 2 3 4 a 0 k hh c hh k mm c mm k mh c mh ν=0.25 ν=0.33 ν=0.45 Figure 3.12 – A Closer Observation of the Horizontal Impedance. based on Green’s Functions have had spurious frequencies in their solutions, this jump of approximately 10% relative to the magnitude of the impedance has to be monitored carefully. It appears this characteristic only appear when the soil medium has a Poisson’s Ratio of 0.25 or less. Typically, approximate Poisson’s Ratios for soil representation are more near 0.50 than 0.25 because an alluvial valley has high moisture content. The value=0:25 is a typical value used for steel and aluminum. After a detailed analysis of the horizontal impedance with a very fine increment in the frequency axis, it is concluded that the change in value is not a characteristic of the present method, but rather, a characteristic of the soil properties or of the Green’s Function calculations. In Fig. 3.12, a closer observation of the horizontal impedance with a reduced frequency axis and a refined frequency spacing is presented. The vertical axis was also magnified to show the details of the change in values. It is clear that only the curve for 48 0123 7 8 9 10 11 12 13 14 a 0 0123 8 9 10 11 12 13 a 0 k hh c hh ν=0.25 ν=0.33 ν=0.45 Figure 3.13 – Comparisons ofK vv andK tt to Published Results. =0:25 has the jump and that it is not a random numerical instability, but rather a characteristic of the underling medium because the jump is actually a smooth pulse in the frequency domain. Shown in the lower figures of Fig. 3.13 is the comparison for the torsional impedance, K tt , for embedment ratios ofh=b=0:25, 0:5, and 1. The results of the present method are shown as solid lines while Apsel’s results are plotted with dashed lines. The left side 49 0 5 10 15 0 5 10 15 0123456 5 10 15 20 25 a 0 0123456 0 5 10 a 0 k vv c vv k tt c tt Current Method Apsel of the figure contains the real part of the impedance or the soil spring,k tt . The right side of the figure contains the imaginery part of the impedance divided by the dimensionless frequency, or the soil damper,c tt . The soil damper from the impedance function includes both the material damping and the effects of radiation damping; the latter is a result from the loss of energy due to spreading of outgoing waves. Since the figures are quite busy, there are no labels for theh=b ratios, it would be correct to assume that the highest curves represent the results of h=b=1, the middle curves represent h=b=0:5 and the lowest curves show the results forh=b=0:25. The torsional results are remarkably similar, much better than the othercomponents to be shown later. One explanation is that for a cylindrical foundation model, the torional component generates only shear waves, therefore, there is no mode conversion between shear and compressional waves. The torsional results are the simplest to obtain as compared to the other components. A closed-form solution for a hemispherical embedded foundation was prepared by Luco (1976). Shown in the upper figures of Fig. 3.13 is a comparison for the vertical impedance, K vv , for three embedment ratios. It clear that the results do not match as well as the case of the torsional component. The largest deviation is in the high frequency range, when dimensionless frequency, a 0 , is larger than 3. Physically, when the value ofa 0 = , the wavelength of the shear wave is equal the width of the foundation, 2b. One immediate concern is the deviations at high frequency. The explanation is not entirely clear because other results generated for a Poisson’s Ratio of 0.45, the real part decreases sharply for high frequency. A general rule about high frequency results is that the imaginery part dominates at high frequency to the order ofk+i!c, therefore, the loss of accuracy for the real part is a common phenomenon. 50 Figure 3.14 – Comparisons ofK hh ,K mm andK mh to Published Results. 51 4 6 8 10 12 0 5 10 15 0 5 10 15 0 2 4 6 8 10 0123456 -1 0 1 2 3 4 a 0 0123456 0 2 4 6 8 10 a 0 k hh c hh k mm c mm k mh c mh Current Method Apsel Another complicated component, the horizontal impedanceK hh is shown in the top figures of Fig. 3.14. For a surface foundation, the horizontal impedance produces only shear waves and the results are predictable like the torsional case. For embedded foundations, however, the leading vertical surface of the foundation generates compressional waves while the bottom surface generates shear waves; opposite to the vertical component. One particular phenomenon of note is that the present results show a small jump neara 0 1:7. This is not present in Apsel’s results and it is also not present in Day’s (1977) results which Apsel benchmarked against. At first, a numerical instability is assumed, but a refined analysis with a very fine frequency increment reveals that it is perhaps a physical characteristic of the soil medium. Apsel’s results and the present results were calculated with a Poisson’s Ratio of 1/4 and that translates to a ratio of shear wave to compressional wave velocities of γ=1= p 3. Results generated for Poisson’s Ratios of 1/3 (γ=1=2) and 0.45 (γ=0:3015) did not result in a jump. It is supposed that for a small value ofγ the compressional waves have a smaller wavelength and that could have contributed to that jump. As an experiment, a low value of 0.20 for the Poisson’s Ratio also resulted in a jump at a similar value ofa 0 . But realistically, the Poisson Ratio for a soil medium, especially for agricultural areas, is much nearer 0.5 than 0.0, therefore, the speed ratio,γ, should be smaller than p 3. Shown in the middle figures of Fig. 3.14 is the Rocking component of the impedance function,K mm . Again, there is no label for the results of different embedment ratios,h=b. The reader could follow the results by knowing that the higher values are from the larger h=b ratios; they are 1.0, 0.5 and 0.25. In bottom figures of Fig. 3.14 the coupling impedance between the rocking and the horizontal components,K mh , is displayed. Like most other components, the real part has deviations at high frequency while the imaginary parts match 52 Figure 3.15 –K vv andK tt for a Square Embedded Foundation. nearly identically. From previous figures, it has been shown that the high frequency range of the real part does not have the same difficulty when the Poisson Ratio is great than 1=4. 3.5 – Results of a Square Embedded Foundation The cylindrical embedded foundation results of Apsel were calculated with axisymmetric sources, known as ring loads. The vertical and torsional components were 53 0 5 10 15 5 10 15 20 0123456 10 15 20 25 a 0 0123456 0 5 10 15 a 0 k vv c vv k tt c tt h/b=0.75 h/b=0.50 h/b=0.25 Figure 3.16 –K hh ,K mm andK mh for a Square Embedded Foundation. 54 0 5 10 15 0 5 10 15 0 5 10 15 20 0 2 4 6 8 10 12 0123456 0 2 4 6 8 a 0 0123456 0 2 4 6 8 a 0 k hh c hh k mm c mm k mh c mh h/b=0.75 h/b=0.50 h/b=0.25 analyzed with zeroth order harmonic for the vertical and azimuthal component of ring forces, respectively. The horizontal, rocking and coupling components were calculated with the first order harmonic for the horizontal and vertical ring forces, respectively. With these axisymmetric properties, the formulation becomes two-dimensional in nature and very accurate results could be obtained. The ring load method, however, cannot be extended to arbitrary shaped foundation. As a simple application of the present method, which is formulated for an arbitrarily shaped foundation, the case of a square foundation with various embedment ratio is presented. The results for thex andy direction are the same except for a sign in the coupling terms, therefore, a more concise presentation can be made. For a rectangular foundation with a different aspect ratio, the two horizontal and the two rocking component would be completely different. Furthermore, a totally nonsymmetrical foundation has a full 6 6 impedance matrix as results and the matrix is frequency dependent. Since the impedance matrix is symmetric, there is a maximum of 21 components for an isolated foundation. The impedance functionsK vv andK tt are shown in Fig. 3.15. In each figure, three ratios of embedment, h=b =0.75, 0.50 and 0.25 are shown. All the results were calculated using a Poisson’s Ratio of 1/3. The numerical solution is stable throughout the frequency domain and the results show that the present method does not have a numerical instability problem. It is amenable for practical engineering application. The main drawback is that a large number of subregion is required to produce excellent results, therefore, the computational effort can be quite tedious. Additional results for the horizontal, rocking and coupling impedances are shown in Fig. 3.16. 55 Figure 3.17 – Horizontal, Vertical and Rocking Input Motion Due to Incident SH-Waves. 56 -2 -1 0 1 2 -2 -1 0 1 2 -.5 .0 .5 -.5 .0 .5 0123456 -1 0 1 a 0 0123456 -1 0 1 a 0 Re U x Im U x Re R y Im R y Re T z Im T z 90 o 60 o 30 o 0 o 3.6 -– Response of a Massless Embedded Foundation to Incident Waves An important component of the soil-structure interaction equation of motion is the driving force vector from the incoming waves. Physically, it is defined as the force required the hold the foundation fixed while subjected to a prescribed incident wave. The driving force vector can also be replaced by the input motion vector and they are related through the impedance matrix. The input motion vector represents, physically, the response of a massless embedded foundation to incident waves. The way the foundation respond to a wave is identified as kinematic interaction. This effect is considered to be an important parameter in a seismic analysis because it reduces the horizontal motion of the foundation as compared to the free field motion. Later results show that while the horzontal component may be reduced for the higher frequencies, the torsional component becomes prominent for oblique incidence. The Modified Ohsaki’s Method has flexibility, the only difference between the impedance calculation and the driving force calculation is the boundary condition. The only requirement is to set displacement condition of the new problem as the negative of the incident wave field within the foundation’s volume. Effectively, this condition sets the displacements within the foundation to zero, i.e., a fixed foundation. The summation of the resulting body forces yields the driving forces. Since the foundation is motionless, there is no need to correct for the intertia forces generated by the movement of the foundation as in the previous case of impedance calculation. To present results which are easier to understand visually, the driving force vectors are converted to input motion vectors by multiplying it by the inverse of the impedance matrix. Three different sets of results are given. The response of the massless foundation toincident body waves of the type SH, P and SV were considered. The free field motion of these body 57 Figure 3.18 – Horizontal, Vertical and Rocking Input Motion Due to Incident P-Waves. 58 -2 -1 0 1 2 -1 0 1 -1 0 1 2 3 -2 -1 0 1 2 0123456 -1 0 1 a 0 0123456 -1 0 1 a 0 Re U x Im U x Re U z Im U z Re R y Im R y 30 o 60 o 90 o waves, plus that of the Rayleigh Surface Wave, are given in Appendix D. A minor change of coordinate axes is needed to use the results in Appendix D, the derivation in the appendix was given with they-axis defined downward and thez-axis into the plane. In the present system, thez-axis is downward and they-axis is outward from the plane. The response of the square massless foundation to a SH-wave with an embedment ratio ofh=b=0:75 is presented in Fig. 3.17. The incident wave travels in they-direction, therefore, the particle motion of the wave is in thex-direction. The horizontal shear wave causes no vertical component but the twisting motion generates a torsional response. Four different incidence angles, , were used. They are 0 , 30 , 60 and 90 , measured with respect to the y-axis. The last angle, 90 , is the vertical incident angle and no torsional component will result from that. Shown in the upper figures of Fig. 3.17 is the horizontal response to SH-waves,U x . The limit to the zero frequeny is2 for an incident wave amplitude of1 because the reflected wave interfere constructively with the incident wave. The real and imaginary parts alternate, an indication there is a phase factor in its response. The rocking component is largest when =90 while the torsional component is largest for=0 . At their peaks, the rocking and twisting angles, multiplied by the half width of the foundation, is about1=4 of the maximum horizontal response. Therefore, the rotational components contribute significantly to the response of the foundation and therefore, the superstructure. There is no mode conversion for the SH-component, the free field motion is of pure shear. Shown in Fig. 3.18 is the response of a massless foundation to an incident P-wave. The wave travels in thex-direction with its particle motion in thez and thex-directions. Except for the vertical incidence of90 , the incident P-wave causes a mode conversion and a reflected SV-wave results. Three angles of incidence are considered, they are 30 , 60 59 Figure 3.19 – Horizontal, Vertical and Rocking Input Motion Due to Incident SV-Waves. 60 -1 0 1 2 3 4 -1 0 1 2 3 -1 0 1 -1 0 1 0123456 -1 0 1 a 0 0123456 -1.5 -1.0 -.5 .0 .5 a 0 Re U x Im U x Re U z Im U z Re R y Im R y 60 o 70 o 80 o and90 . measured with respect to thex-axis. The horizontal and vertical response to these incident P-waves is given in a table in Appendix D. For the case of the vertical incident wave, there is no induced rocking. The wavelength of the incident P-wave is twice as long as that of the incident S-waves, therefore, the kinematic interaction is less prominent. Shown in Fig. 3.19 is the response of a massless foundatio to an incident SV-wave. The wave travels in thex-direction with its particle motion in thez andx-directions. The SV-wave has a very limited range of angles which are admissable, therefore, the angles used are60 ,70 and80 . The vertical incident angle of90 was not used because it would be the same as that of the SH-wave. Similar to the P-wave case, only the horizontal, vertical and rocking components are excited. The free surface amplitudes of these incident SV-waves are listed in Appendix D. In summary, the input motion vectors, along with the impedance matrix, completes the substructural information required for a soil-structure interaction analysis. The response of a structure to various incident waves can be calculated first in the frequency domain and then transformed into the time domain using Fourier Transformation. 61 Chapter Four Application Two, the Boundary Integral Equation Method As shown in Eq. (2.17), a boundary integral equation can be written using the reciprocal theorem. By eliminating the body force term, only the tractions and the displacements at the surface of the wave scatterer need to be analyzed. The limit of the source point~ r p onto surfaceS leaves an integral equation in the form 1 2 ~ u(~ r p )= Z S U(~ r p j~ r) ~ t(~ r)dS − Z S T(~ r p j~ r) ~ u(~ r)dS ; ~ r p onS: (4:1) In the above equation, the displacement vector at any point on the surface, ~ r p , can be calculated as a surface integral of the displacements and tractions at the same surface multiplied by the traction and displacement Green’s Functions, respectively. Eq. (4.1) is a vector formulation of a boundary value problem which yields a unique solution to the vector wave equation. After the boundary values are determined, Eq. (2.16) can be used to calculate the wave solution at other locations away from the surfaceS. The major advantage of this formulation is that only the surface of the wave scatterer is involved, thus reducing a three-dimensional problem to two dimensions. Another advantage is the elegance of the formulation which has attracted competent researchers from mathematics, physics, and engineering to devote significant effort to its study over the past 30 years. The disadvantage of this formulation is the instability of the numerical solution. Spurious resonant frequencies of the corresponding interior problem appear in the solution of the exterior problem, rendering the solution useless for many practical applications. An attempt will be made in this chapter to clearly understand the numerical instability based 62 on previous solutions and devise an improved workable solution to eliminate the interior resonant frequencies. As an application to find the impedance of a rigid, embedded foundation, let the surfaceS of Eq. (4.1) be an embedded, rigid foundation and ~ t(~ r) be the unknown traction which would constrain the displacement vector~ u(~ r) on the same surface to move as a rigid body. The six modes of rigid body motion include three translational, two rocking, and one torsional modes. In a separate formulation to obtain the driving force vector induced by incident wave motion, Eq. (4.1) can be used to determine the traction distribution that would hold the foundation surface fixed while being subjected to the incoming wave motion. Figure 4.1 – Subdivision of the Rigid Embedded Foundation Surface. 63 ~ r p ~ t(~ r);~ u(~ r) x z y A numerical solution of Eq. (4.1) can be accomplished by first subdividingS intoN subregions,S i ;i=1;2;:::;N. Then, Eq. (4.1) can be expressed as 1 2 ~ u(~ r p )= N X j=1 Z S j U(~ r p j~ r) ~ t(~ r)dS − N X j=1 Z S j T(~ r p j~ r) ~ u(~ r)dS ; (4:2) with ~ r p placed on S. The next step is to approximate ~ t(~ r) as a constant vector, ~ t j and ~ u(~ r) as a constant vector,~ u j ,onS j , as a zeroth order approximation. It has been shown in previous formulations of this type (Wong 1975) that the higher order approximations often do not improve the accuracy because the traction distribution has extremely high values at the edges. Therefore, a constant value which represents the average traction distribution in the subregions provides a good approximation. With the above assumptions, Eq. (4.2) can be simplified to N X j=1 Z S j [U(~ r p j~ r)]dS ! ~ t j = 1 2 ~ u(~ r p )+ N X j=1 Z S j [T(~ r p j~ r)]dS ! ~ u j : (4:3) Assign the source point of the Green’s Function,~ r p , to be equal to~ r i , the centroid of the subregions,S i , then Eq. (4.3) can be expressed as a 3 3 matrix equation N X j=1 [U ij ] ~ t j = 1 2 ~ u i + N X j=1 [T ij ]~ u j ;i=1;2;:::;N; (4:4) in which [U ij ]= Z S j [U(~ r p j~ r)]dS ; (4:5) and [T ij ]= Z S j [T(~ r p j~ r)]dS ; (4:6) will be referred to as the influence matrices. Replace the summation in Eq. (4.4) by a matrix equation of order 3N 3N and it can be expressed as U ~ t = 1 2 I + T ~u; (4:7) 64 in which~ u contains the displacement vectors at the centroids of theN subregions. To adapt Eq. (4.7) to the calculation of the impedance matrix of an embedded foundation, let the displacement vector, ~ u, be defined by rigid body motion as defined for each subregionS i as ~ u i =[ i ] ~ U 0 ; (4:8) in which ~ U 0 is the six-degree-of-freedom foundation motion defined as ~ U 0 = 2 6 6 6 6 6 4 U x U y U z b x b y b z 3 7 7 7 7 7 5 ; (4:9) and [ i ] is the rigid-body kinematic matrix defined in this manner: [ i ]= 2 4 100 0 z i =b −y i =b 010 −z i =b 0 x i =b 001 y i =b −x i =b 0 3 5 : (4:10) The parameterb in Eq. (4.9) and Eq. (4.10) has a unit of length related to the size of the foundation so that all entries in the generalized displacement vector,~ u, have units of length. The solution of the matrix equation, Eq. (4.7), can be obtained numerically as t r = U −1 1 2 I + T ~ U 0 ; (4:11) in which [t r ] is a 3N 6 matrix containing the 6 columns of traction distributions which satisfy the displacement compatibility conditions specified by the 3N 6 matrix []. The matrix [ i ] in Eq. (4.10) is a subset of [] in Eq. (4.11). The latter contains the rigid body motion of all the centroids of theN subregions. 65 The final step of the formulation is to accumulate the effects of the surface tractions into three forces and three moments using the transpose of the rigid body kinematic matrix: ~ F 0 = T t r = T U −1 1 2 I + T ~ U 0 (4:12) in which ~ F 0 = 2 6 6 6 6 6 4 F x F y F z M x =b M y =b M z =b 3 7 7 7 7 7 5 ; (4:13) is the generalized force, including three forces and three moments, normalized to the unit of force. On the right hand side of Eq. (4.12), the matrix product, K = T U −1 1 2 I + T (4:14) is the definition of the 6 6 impedance matrix. 4.1 – Numerical Procedure for the Integration of Singularities The integrals in Eqs. (4.5) and (4.6) can be evaluated accurately by numerical quadrature except for the diagonal elementsU ii andT ii , because the source point,~ r p , and the observation point,~ r, are located in the same surface subregionS i . The singularities of the displacement Green’s Functions are of the order 1=R and they are easily integrable analytically, but it is an illegal operation to numerically divide by zero. The singularities of the traction Green’s Functions are of the order 1=R 2 and they are not integrable even analytically. For this formulation, only the principal value integral is required. As shown in Chapter Three, the form of the singularities for all Green’s Functions are the same regardless of whether the problem is dynamic or static, or whether the model has 66 an infinite space configuration or a semi-infinite space configuration. Using this knowledge, U ii of Eq. (4.5) can be evaluated in two parts: [U ii ]= Z S i [U(~ r p j~ r)]dS = Z S i [U(~ r p j~ r)] −[U (~ r p j~ r)] dS + Z S i [U (~ r p j~ r)]dS ; (4:15) in which [U (~ r p j~ r)] is any appropriate Green’s Function matrix which can be expressed in closed form. For this particular application, the static point load solution in an infinite space will be used (Love, 1906). With the procedure introduced in Eq. (4.15), the first integral can now be evaluated numerically since it is void of singularities and the second integral has a reasonably simple closed form solution. The static point load solution by Love (1906) was shown in Eq. (3.16) The integral of U is evaluated over a surface in this application as opposed to over a volume, like Application One. Therefore, depending on the orientation of the surface normal atS i , the static integral varies. 4.2 – Application of Green’s Functions Obtained in Cylindrical Coordinates As shown in Chapter Three, there is a need to adapt the Green’s Functions developed in the cylindrical coordinate system to the Cartesian coordinate system used by the embedded foundation problem. For the integral equation method, the displacement Green’s Function f rz andf zz for the vertical point load andf rr ,f r andf zr for the horizontal point load are needed, as was the case in Chapter Three. Additionally, the stress Green’s Functions, g rr ,g ,g zz andg rz for the vertical point load andg rr ,g ,g zz ,g r ,g rz andg z for the horizontal point load are also required. The functionsf ii are all regular functions and they 67 Figure 4.2 – Definition of the Cylindrical Coordinate System. are to be normalized by the factor 1=(r) while the functionsg ii are to be normalized by the factor 1=(r 2 ). Now consider the coordinate systems illustrated in Figure 4.2. Again, thez-axis of both systems is defined pointing downward. The origin of the polar coordinate system is defined at the source point, ~ r p , therefore, the observation point, ~ r, is located at the coordinates, (r; ;z), in which r = j~ r −~ r p j = q (x −x p ) 2 +(y−y p ) 2 ; (4:16) and = arg(~ r −~ r p ) = tan −1 y −y p x −x p : (4:17) Thez-dependency of the Green’s Functions is included in the functionsf andg. 68 r x y z 0 P r P z P y P x Using the polar coordinate system, all horizontal point forces can be represented by P r ( 0 ) because the reference angle can be varied to match any orientation. Therefore, the general displacement-force relationship can be written in a form of a 3 2 matrix: 8 < : u r (r; ) u (r; ) u z (r; ) 9 = ; = 1 r 2 4 f rr cos( − 0 ) f rz f r sin( − 0 )0 f zr cos( − 0 ) f zz 3 5 P r ( 0 ) P z : (4:18) To obtain a displacement-force relationship in the Cartesian coordinate system, the following mapping of the forces and the displacements may be used: 8 < : P x (~ r p ) P y (~ r p ) P z (~ r p ) 9 = ; = 8 < : P r (0) P r (=2) P z 9 = ; ; (4:19) and 8 < : u x (~ r) u y (~ r) u z (~ r) 9 = ; = 2 4 cos −sin 0 sin cos 0 00 1 3 5 8 < : u r (r; ) u (r; ) u z (r; ) 9 = ; : (4:20) During the first step of the transformation, the relationship in Eq. (4.19) can be substituted into Eq. (4.18), then expanded by one column to relate the displacements in polar coordinates to the forces in Cartesian coordinates: 8 < : u r (r; ) u (r; ) u z (r; ) 9 = ; = 1 r 2 4 f rr cos f rr sin f rz f r sin −f r cos 0 f zr cos f zr sin f zz 3 5 8 < : P x (~ r p ) P y (~ r p ) P z (~ r p ) 9 = ; ; (4:21) Now apply the transformation (4.20) to both sides of Eq. (4.21), the result is a displacement- force relationship in Cartesian coordinates written as 8 < : u x (~ r) u y (~ r) u z (~ r) 9 = ; = 1 r G 8 < : P x (~ r p ) P y (~ r p ) P z (~ r p ) 9 = ; ; (4:22) in which G = 2 6 4 f rr cos 2 −f r sin 2 (f rr +f r )sin cos f rz cos (f rr +f r )sin cos f rr sin 2 −f r cos 2 f rz sin f zr cos −f zr sin f zz 3 7 5 : (4:23) 69 Using Cartesian coordinates, the factors, sin and cos , can be evaluated simply as sin = y −y p r ; (4:24) and cos = x −x p r : (4:25) If the elements of matrix [G] are defined with the following subscripts, G = 2 4 G xx G xy G xz G yx G yy G yz G zx G zy G zz 3 5 ; (4:26) then the matrix [U(~ r p j~ r)] in Eq. (4.5) can be formed as U(~ r p j~ r) = 1 r 2 4 G xx G yx G zx G xy G yy G zy G xz G yz G zz 3 5 : (4:27) A column in[G] represents the displacements caused by a point load in a certain direction, in the orderx,y, andz, respectively. In the matrix[U(~ r p j~ r)], however, the same displacement components are stored as a row. The reason for this transposition can be observed in Eq. (2.2) in which the Green’s Functions entry, i.e., ~ u T , is written as a row vector in preparation for an inner product. The stresses caused by a vertical point load can be written in the form, (P z ) = 1 r 2 2 4 g rr 0 g rz 0 g 0 g rz 0 g zz 3 5 ; (4:28) while the stresses caused by a horizontal point load in the 0 direction, i.e.,P r ( 0 ), are (P r ( 0 )) = 1 r 2 2 4 g rr cos( − 0 ) g r sin( − 0 ) g rz cos( − 0 ) g r sin( − 0 ) g cos( − 0 ) g z sin( − 0 ) g rz cos( − 0 ) g z sin( − 0 ) g zz cos( − 0 ) 3 5 : (4:29) 70 Using Eq. (4.29), the stresses created by a point load in thex-direction can be obtained by setting 0 to 0, in the form, (P x ) = 1 r 2 2 4 g rr cos g r sin g rz cos g r sin g cos g z sin g rz cos g z sin g zz cos 3 5 ; (4:30) and the stresses created by a point load in they-direction can be obtained by setting 0 to =2, in the form, (P y ) = 1 r 2 2 4 g rr sin −g r cos g rz sin −g r cos g sin −g z cos g rz sin −g z cos g zz sin 3 5 : (4:31) To calculate the tractions required for the matrix [T] in Eq. (4.6), the stress tensors in Eqs. (4.28), (4.30) and (4.31) must be transformed in the xyz through two orthogonal transformations of Eq.(4.20). Numerically, it is more efficient to transform vectors than matrices, therefore, it is suggested that the surface normal vector,~ e n , be first transformed to the cylindrical coordinate system by the equation 8 < : n r (r; ) n (r; ) n z (r; ) 9 = ; = 2 4 cos sin 0 −sin cos 0 00 1 3 5 8 < : n x (~r) n y (~r) n z (~r) 9 = ; ; (4:32) then the normalized traction vectors at the observation point (r; ;z) subjected toP x ,P y andP z , are, respectively, 8 < : t rx t x t zx 9 = ; = 2 4 g rr 0 g rz 0 g 0 g rz 0 g zz 3 5 8 < : n r (r; ) n (r; ) n z (r; ) 9 = ; ; (4:33) 8 < : t ry t y t zy 9 = ; = 2 4 g rr cos g r sin g rz cos g r sin g cos g z sin g rz cos g z sin g zz cos 3 5 8 < : n r (r; ) n (r; ) n z (r; ) 9 = ; ; (4:34) and 8 < : t rz t z t zz 9 = ; = 2 4 g rr sin −g r cos g rz sin −g r cos g sin −g z cos g rz sin −g z cos g zz sin 3 5 8 < : n r (r; ) n (r; ) n z (r; ) 9 = ; : (4:35) 71 The next step is to convert the traction vectors in cylindrical coordinates to Cartesian coordinates using the transformation, 2 4 H xx H xy H xz H yx H yy H yz H zx H zy H zz 3 5 = 2 4 cos −sin 0 sin cos 0 00 1 3 5 2 4 t rx t ry t rz t x t y t z t zx t zy t zz 3 5 ; (4:36) and finally, the matrix [T(~ r p j~ r)] in Eq. (4.6) can be formed as T(~ r p j~ r) = 1 r 2 2 4 H xx H yx H zx H xy H yy H zy H xz H yz H zz 3 5 : (4:37) 4.3–ATwo-Dimensional SH-Wave Example The simplest case in elastodynamic wave propagation is the two-dimensional SH- wave problem. It has just one antiplane component and the SH-wave equation is directly analogous to the scalar acoustic wave equation. The Green’s Function for this problem is easily expressed in closed form for time harmonic waves of the forme i!t . In this particular case, the displacement Green’s Function matrix has only one element, [U(~ r p j~ r)] =U zz (x p ;y p jx;y)= i 4 H (2) 0 (kR 1 )+H (2) 0 (kR 2 ) ; (4:38) in which is the shear modulus, H (2) 0 is the Hankel Function of the second kind and of zeroth order, and the distances,R 1 andR 2 , are calculated as R 1 = q (x p −x) 2 +(y p −y) 2 ; and R 2 = q (x p −x) 2 +(y p +y) 2 ; respectively. The traction Green’s Function matrix also has only one element, it can be computed as [T(~ r p j~ r)] =T zz (x p ;y p jx;y)= zx n x + zy n y ; (4:39) 72 in whichn x andn y are the two-dimensional direction cosines of the outward normal atS and the stresses, zx = @U zz @x = ik 4 H (2) 1 (kR 1 ) @R 1 @x +H (2) 1 (kR 2 ) @R 2 @x ; (4:40a) zy = @U zz @y = ik 4 H (2) 1 (kR 1 ) @R 1 @y +H (2) 1 (kR 2 ) @R 2 @y ; (4:40b) in whichH (2) 1 is the Hankel Function of the second kind and of first order, and the partial derivatives of the distances can be calculated simply as @R 1 @x = x p −x R 1 ; @R 2 @x = x p −x R 2 ; @R 1 @y = y p −y R 1 ; @R 2 @y = y p +y R 2 : The information above can now be used in Eq. (4.4) and the vector at each node is now a scalar, as there is only one antiplane component to be considered. The surface integrals in Eq. (4.5) and Eq. (4.6) can be reduced to line integrals for a two-dimensional geometry. The diagonal elements, U ii of Eq. (4.5) and T ii of Eq. (4.6), must be treated with care. The displacement Green’s Function in Eq. (4.38) has a log(R) singularity as~ r approaches ~ r p . The scheme presented in Eq. (4.15) can be applied without difficulty because log(R) is easily integrable analytically. The traction Green’s Function as shown in Eq. (4.39) and Eq. (4.40) has a 1=R singularity and is not integrable even analytically. However, only the principal value integral is required and the integrand is extremely well-behaved, except when~ r is exactly at~ r p . The limit is taken there and the result is -1/2 for a smooth surface, which explains why the integral equation in Eq. (4.1) has a factor of 1/2. The other 1/2 was removed by the 1=R singularity of the traction Green’s Function. Shown in Fig. 4.3 is an SH-wave model of a circular cylindrical foundation. It is the most basic of all models, but it will be used here to illustrate a fundamental problem which exists for boundary integral equation formulations. The radius of the cylinder isa and the 73 Figure 4.3–ATwo-Dimensional Circular Cylindrical Foundation. results will be presented using the dimensionless frequencya 0 =!a=, in which! is the frequency in rad/sec and is the SH-wave velocity. Shown in Fig. 4.4 are four plots of the antiplane component of impedance as a function of the dimensionless frequencya 0 . The four plots represent the models created with the number of subregions, N, set to 10, 20, 40, and 80. The real part of the impedance is plotted with a solid line while the imaginary part, divided bya 0 , is plotted with a dashed line. Physically, the real part is likened to a spring and the imaginary part, after it is divided by frequency, is likened to a viscous dashpot. The results for the case whenN =10 are totally unstable, but the backbones of the curves do show some correct trends of the classical solution. The results for the case when N =20 are much improved, but the instability is still too great for practical application. Further improvement can be made by continuing to increase the number of subregions to 40 and then to 80, but no matter how largeN might be, there are locations of instability. The 74 a y x Figure 4.4 – Antiplane Impedance Function as a Function of Frequency. for Various Embeddment Ratiosh=b. 75 0 1 2 3 4 02468 10 0 1 2 3 4 Dimensionless Frequency a 0 02468 10 Dimensionless Frequency a 0 Real (K) Imag (K) / a 0 Real (K) Imag (K) / a 0 N=10 N=20 N=40 N=80 case whenN = 160 is not shown here, but it can be reported that the numerical instability still exists. The frequencies at which this numerical instabilty occurs are the resonant frequencies of the interior problem. The determinant of the [U] matrix in Eq. (4.7) is zero at those frequencies. It is puzzling why the frequencies of the interior model appear in the solution of the exterior model. Many researchers have been frustrated by this phenomenon for an extended period of time. 4.4 – Clues to Numerical Instability Problem For the past 30 years, many researchers have recognized the instability problem of the boundary integral equation formulation. Most recognized that the troubled frequencies are those of the interior problem with the same scattering surface. In this section, an explanation will be given based on the Hilbert-Schmidt classical solution for SH-wave scattering. This method formulates the integral equation solution by expanding the Green’s Function as a linear combination of orthogonal functions which are specific to the scatterer’s surface. The actual procedure is similar to the method of separation of variables. A recap of this method was given in Pao and Mow (1971) for an SH-wave example in an infinite space. The Green’s Function for an infinite space, somewhat similar to that expressed in Eqn. (4.38), is G(~ r p j~ r)=− i 4 1 X m=0 cosm( − 0 ) 8 < : J m (ka)H (2) m (kr);r>r 0 ; J m (kr)H (2) m (ka);r<r 0 ; ; (4:41) in whicha is the radius of the cylindrical wave scattering surface and is the parameter which describes the surfaceS. As shown in Eq. (4.41), the Green’s Functions for the interior problem and the exterior problem are very similar. Only the direction of their normal vectors differ, one pointing out, the other pointing in. 76 Using the same orthogonal function, the incident wave field can be expressed as u i (~ r p )=S 1 X m=0 m J m (kr)cosm ; (4:42) in which 0 =1 and m =2 form 6=0. Express the solution for the surface traction as a linear combination of orthogonal functions with unknown coefficientsB n , in this manner, t(~ r)= 1 X n=0 B n cosn 0 : (4:43) With all the prepared expansions, let dS = ad 0 , and perform the integration over the surfaceS: u i (~ r p )=− i 4 Z 2 0 1 X n=0 B n cosn 0 !" X m J m (ka)H (2) m (kr)cosm( − 0 ) # ad 0 = − ia 4 X n B n J n (ka)H (2) n (kr)cosn ; r onS (4:44) Replacingu i (~ r) by its series approximation, we have the algebraic equation S 1 X m=0 m J m (ka)cosn = − ia 4 1 X n=0 B n J n (ka)H (2) n (ka)cosn ; (4:45) The orthogonality of the eigenfunctions allows each coefficient to be determined uniquely and independently, B n = 4iS n aH (2) n (ka) : (4:46) The simplicity of the last step allows many problems in applied mathematics to be solved routinely. The termJ m is cancelled from both sides of the equation after the inner product is performed to separate the terms. It is mathematically straightforward to cancel a term from each side of the equation. Numerically, however, theJ m function goes to zero at certain frequencies, more specifically, the resonant frequencies of the interior problem. If one of 77 the eigenvalues of the left-hand-side matrix [U] becomes zero, the determinant of [U] also becomes zero. 4.4.1 – Application of l’Hospital’s Rule One of the major discoveries in this dissertation is the fact that the matrices of each side of Eq. (4.7) are Normal Matrices, i.e., U 1 2 I + T = 1 2 I + T U : (4:47) This fact has been confirmed numerically. This shouldn’t come as a surprise, though, because a matrix product is a collection of inner products. The inner product of displacement and traction forces, regardless of the order in which it is performed, is equivalent to energy. The most useful property of Normal Matrices is that they share the same eigenvalues and eigenvectors. Applying this knowledge, numerical schemes can be designed based on an eigenvalue solution. Since both matrices are frequency dependent, it is not as simple as the eigenvalue problem for constant matrices. However, it is reasonable to assume that an eigenvector is changing slowly within a narrow band of frequencies. In fact, the eigenvector for a circular geometry is not a function of frequency, though it is expected that a more general geometry would yield frequency dependent eigenvectors. If this assumption is acceptable, let ~ be the eigenvector of interest, and its product with the matrices would yield ~ T [U(!)] ~ =f(!) ; (4:48) ~ T 1 2 [I]+[T(!)] ~ =g(!) ; (4:49) 78 and at !, bothf( !) andg( !) are zero. If the matrices are available for two other frequencies, ! 1 and! 2 , one on each side of !, then an approximation of the limit can be taken using l’Hospital’s rule as q = g(! 2 ) −g(! 1 ) f(! 2 ) −f(! 1 ) ; (4:49) and the matrix equation in Eq. (4.7) can be repaired by U( !) + ~ (c) ~ T ~ t= 1 2 I + T( !) + ~ (cq) ~ T ~u; (4:50) in whichc is an arbitrary constant. The repair factor will change the determinant of[U( !)] to a finite non-zero value. Since the same change is made to the right-hand-side, the ratio for the contribution of that particular eigensolution is correct according to Eq. (4.49). This method is a proposal for future work in a more vigorous fashion. Most of the effort of this dissertation was devoted to the development of Application One, the Modified Ohsaki method. The contribution made in this chapter could provide a way for the boundary integral equation method to eventually be stabilized. 79 Chapter Five Application Three, a Radiation Boundary The equation of motion for a three-dimensional elastodynamic problem can be expressed in terms of horizontal displacements, u and v, and vertical displacement, w, as @ 2 u @t 2 =(+) @ @x + @ 2 u @x 2 + @ 2 u @y 2 + @ 2 u @z 2 ; (5:1a) @ 2 v @t 2 =(+) @ @y + @ 2 v @x 2 + @ 2 v @y 2 + @ 2 v @z 2 ; (5:1b) @ 2 w @t 2 =(+) @ @z + @ 2 w @x 2 + @ 2 w @y 2 + @ 2 w @z 2 ; (5:1c) in which and are the Lame’ constants and is the three-dimensional dilation defined as = @u @x + @v @y + @w @z : (5:2) The Lame’ constants are related to the Young’s modulus,E, and the Poisson’s ratio,, E = (3+2) + ; (5:3) and = 2(+) ; (5:4) respectively. The strong coupling between the three equations for the three-dimensional problem in Eqs. (5.1) makes the implementation of a finite difference model extremely difficult. The impracticality stems from the coupling of the displacement components in 80 Figure 5.1 – A Three-Dimensional Finite Element Model. the boundary conditions for the tractions because the stresses are defined in this manner, xx =+2 @u @x ; (5:5a) yy =+2 @v @y ; (5:5b) zz =+2 @w @z ; (5:5c) xy = @u @y + @v @x ; (5:5d) yz = @v @z + @w @y ; (5:5e) zx = @w @x + @u @z ; (5:5f) Therefore, even a simple stress-free boundary condition is difficult to formulate. For most cases, other than the two-dimensional SH-wave models or acoustic wave problems, it is much simpler to use a finite element model. 81 5.1 – Three-Dimensional Finite Element Formulation Consider a three-dimensional embedded foundation model in a semi-infinite space as illustrated in Fig. 5.1. The model is constructed using isoparametric block elements as described in Appendix G (Zienkiewicz, 1971). The 8-node element stiffness matrix and element mass matrix have a size of 24 24, with 3 displacements per node. Using a grid generator to prepare the global and local node indexing scheme, an undamped time-harmonic model can be represented by an equation of motion of the form, −! 2 h M i + h K i ~ u = ~ f; (5:6) in which[M] is the global mass matrix,[K] is the global stiffness matrix and ~ f is the nodal force vector. The simplicity of the finite element formulation is obvious as a stress-free boundary at the spherical cavity can be implemented by merely setting the nodal forces at the free boundary nodes to zero, as shown in Fig. 5.1. Figure 5.2 – Nodal Definition for Radiation Boundary. 82 Subset A Subset S Subset E Subset I For most applications, the nodal numbering scheme is optimized in order to minimize the bandwidth of[K]. To implement the radiation boundary, however, it is better to number the nodes by subsets. Start with subsetE for the exterior nodes, subsetA for the artificial boundaries nodes, subsetS for an additional layer of nodes for the purpose of calculating stresses and tractions, and subsetI for all other interior nodes of the model. The organization of the nodes is shown in Fig. 5.2. Eq. (5.6) can be rewritten in the form, [Z]~ u = ~ f; (5:7) in which[Z] is a complex, frequency-dependent impedance matrix which includes stiffness, damping, and inertial effects of the soil medium. The forcing function, ~ f, can be used to simulate internal forces or to set boundary conditions. For incident wave scattering problems, the force vector can be set opposite to the incident wave stresses to create a stress- free boundary. For the latter case, the solution of Eq. (5.7) would then be the scattered wave field. The sum of the incident wave field and the scattered wave field would constitute the total solution. According to the pre-assigned nodal order, Eq. (5.7) can be partitioned as follows: 2 6 6 6 6 6 6 6 4 [Z EE ][Z EA ] [0] [0] [Z AE ][Z AA ][Z AS ] [0] [0] [Z SA ][Z SS ][Z SI ] [0] [0] [Z IS ][Z II ] 3 7 7 7 7 7 7 7 5 2 6 6 6 6 6 6 6 4 ~ u E ~ u A ~ u S ~ u I 3 7 7 7 7 7 7 7 5 = 2 6 6 6 6 6 6 6 4 ~ 0 ~ 0 ~ 0 ~ F I 3 7 7 7 7 7 7 7 5 ; (5:8) in which ~ F I is the forcing function applied to the interior portion of the model. It is assumed that no forces are applied at the artificial boundaries, or exterior to them. 83 5.2 – Creation of a Radiation Boundary Each node in the finite element model is coupled to adjacent nodes through the stiffness matrix and sometimes through the consistent mass matrix or damping matrix. For example, one node in subsetA is related to nodes from ubsetsE,A, andS. Other interior nodes have similar surrounding nodes. The nodes in subsetE, however, have missing adjacent nodes which would allow propagation of outgoing waves to the far-field. To remove the problems caused by the truncation of a model by an artificial boundary, it is proposed to calculate the displacement,~ u E , as defined in Eq. (2.25), using the Green’s functions of a semi-infinite medium and a representation theorem in the form, ~ u(~ r E )= Z S h U(~ r E j~ r A ) i ~ t(~ r A ) − h T(~ r E j~ r A ) i ~ u(~ r A ) dS ; (5:9) in which ~ u(~ r A ) and ~ t(~ r A ) are the displacement and traction vectors at the artificial boundaries. The integration is to be performed over a radiating boundary, S, and in this case, the artificial boundaries. The Green’s Function matrices,[U(~ r E j~ r A )] and[T(~ r E j~ r A )], contain the displacement and traction values, respectively, at observer~ r A generated by a point source located at~ r E . Eq. (5.9) was derived based on the Betti-Rayleigh reciprocal theorem; it allows the displacement at a node in subsetE to be computed by integrating the effects of the displacement and traction values at a radiative surface,S, based on Huygen’s Principle. The detailed derivation of Eq. (5.9) is given in Appendix A. The above formulation allows the displacements,~ u E , to be calculated without nodes exterior to model. This allows the waves to propagate outward, based on the motion of the interior nodes. To calculate the traction at a given location, an outward normal to the exterior medium is defined as a unit vector toward the interior nodes from the artificial boundary. The reason for this convention is that for the outgoing waves, the exterior medium is the propagating 84 medium, the finite element model of interest is actually outside of the propagating medium. The relationship between the traction vector, ~ t, and the stress tensor is ~ t = 2 4 t x t y t z 3 5 = 2 4 xx xy xz yx yy yz zx zy zz 3 5 2 4 n x n y n z 3 5 ; (5:10) in whichn x ,n y andn z are the direction cosines for the outer normal vector, ^ e n . For three- dimensional problems, Eq. (5.9) is quite tedious, it has nowx,y andz components. The derivation of the Green’s Function matrix of a semi-infinite medium for point loads in three directions was given by Luco and Apsel (1979). In this particular case, the displacement Green’s Function matrix is [U(~ r E j~ r A )] = 2 4 U 11 (~ r E j~ r A ) U 12 (~ r E j~ r A ) U 13 (~ r E j~ r A ) U 21 (~ r E j~ r A ) U 22 (~ r E j~ r A ) U 23 (~ r E j~ r A ) U 31 (~ r E j~ r A ) U 32 (~ r E j~ r A ) U 33 (~ r E j~ r A ) 3 5 ; (5:11) and the traction Green’s Function matrix is [T(~ r E j~ r A )] = 2 4 T 11 (~ r E j~ r A ) T 12 (~ r E j~ r A ) T 13 (~ r E j~ r A ) T 21 (~ r E j~ r A ) T 22 (~ r E j~ r A ) T 23 (~ r E j~ r A ) T 31 (~ r E j~ r A ) T 32 (~ r E j~ r A ) T 33 (~ r E j~ r A ) 3 5 ; (5:12) in which the matrix elements,T ij , can be calculated as 2 4 T 11 T 12 T 13 3 5 = 2 4 xx X xy X xz X yx X yy X yz X zx X zy X zz X 3 5 2 4 n x n y n z 3 5 ; (5:13) 2 4 T 21 T 22 T 23 3 5 = 2 4 xx Y xy Y xz Y yx Y yy Y yz Y zx Y zy Y zz Y 3 5 2 4 n x n y n z 3 5 ; (5:14) and 2 4 T 31 T 32 T 33 3 5 = 2 4 xx Z xy Z xz Z yx Z yy Z yz Z zx Z zy Z zz Z 3 5 2 4 n x n y n z 3 5 : (5:15) In Eq. (5.13) through (5.15),n x ,n y andn z are the direction cosines of the outer normal vector of the boundary surface at the point of integration and the stress matrices are those generated by forces in theX,Y , andZ directions, respectively. 85 For numerical calculation, approximate the integrals in Eq. (5.9) by numerical quadratures of the form, Z S f(~ r A )g(~ r A )dS = X j w j f(~ r j )g(~ r j ) ; (5:16) in whichw i are weighting factors for integration and~ r j contains the sample points selected for numerical integration. For most applications in this chapter, an equally spaced two- dimensional Trapezoidal Rule is sufficiently accurate for rectangular artificial boundaries. For a cubical model, there are six, two-dimensional surfaces. The values of the weights are w i =h 2 =4 at the corner of a rectangular integration surface,S, andw i =h 2 =2 at the edges andw i =h 2 for all interior points of the rectangular area. The value ofh is the grid spacing in the direction of integration, assumed to be the same in both lateral directions, without any loss of generality. With this discretization scheme in place, the right hand side of Eq. (5.9) can be written as matrix products: ~ u E = −[T]~ u A +[U] ~ t A ; (5:17) in which the matrices[T] and[U] are generally not square in dimension because the number of nodes in subsetE is usually not equal to that of subsetA. The vector~ u A in Eq. (5.17) is one of the unknown vectors in matrix equation Eq. (5.8), therefore, the first term on the right hand side of Eq. (5.17) can be simply implemented by substitution. The traction vector ~ t A , on the other hand, has to be computed from the stress distributions with the formula ~ t(~ r A )= 2 4 t x t y t z 3 5 = 2 4 xx (~ r A ) xy (~ r A ) xz (~ r A ) xy (~ r A ) yy (~ r A ) yz (~ r A ) xz (~ r A ) yz (~ r A ) zz (~ r A ) 3 5 2 4 n x n y n z 3 5 ; (5:18) 86 and that requires the stresses at the nodes in subset A and the direction cosines of the boundary at the same corresponding locations,~ r A . To obtain a second order approximation for the stresses, similar to that provided by the Central Difference Formula in Chapter 3, it is proposed to use the values of the average stresses of the eight surrounding block elements. The concept of a stress average was explained in Chapter 4 by Eq. (4.18). For an isoparametric block element, the stresses can be calculated as ~ = 2 6 6 6 6 6 4 xx yy zz xy yz zx 3 7 7 7 7 7 5 =[D][L][N]~u; (5:19) in which [D]= E(1 −) (1+)(1 −2) 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 1 (1 −) (1 −) 000 (1 −) 1 (1 −) 000 (1 −) (1 −) 1 000 000 (1 −2) 2(1 −) 00 000 0 (1 −2) 2(1 −) 0 000 0 0 (1 −2) 2(1 −) 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ; (5:20) 87 h L i = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 @ @x 00 0 @ @y 0 00 @ @z @ @y @ @x 0 0 @ @z @ @y @ @z 0 @ @x 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ; (5:21) [N]= h [N 1 ][N 2 ][N 3 ][N 4 ][N 5 ][N 6 ][N 7 ][N 8 ] i ; (5:22) [N i ]= 2 4 N i (x;y;z)0 0 0 N i (x;y;z)0 00 N i (x;y) 3 5 ; (5:23) and the shape functions, as defined in Appendix G in terms of the natural coordinates and , are N i (;;)= 1 8 (1+ i )(1+ i )(1+ i ) : (5:24) The range of the natural coordinates is defined as, −1 1, −1 1 and −1 1. In matrix form, the traction vector ~ t A can be represented in terms of~ u E ,~ u A and~ u S as ~ t A =[D E ]~ u E +[D A ]~ u A +[D S ]~ u S ; (5:25) in which the matrices, [D E ], [D A ] and [D S ] contain weighting factors that include the derivative computation, geometrical information, and material constants of the elements near the artificial boundary. Substitution of Eq. (5.25) into Eq. (5.17) yields a new expression for~ u E : ~ u E = −[T]~ u A +[U]([D E ]~ u E +[D A ]~ u A +[D S ]~ u S ) : (5:26) 88 It is clear that ~ u E appears on both sides of Eq. (5.26). The expression can be written in explicit form by introducing [V E ]=[I]−[U][D E ] ; (5:27) ~ u E is expressed explicitly as ~ u E =[V E ] −1 (−[T]+[U][D A ])~ u A +[V E ] −1 [U][D S ]~ u S : (5:28) The complex matrix equation for a finite element model, equipped with a radiation boundary with the unknown variable ~ u E eliminated, now has a new left hand side given by 2 6 6 6 4 [Z AA + Z AE V −1 E (−T + UD A )] [Z AS + Z AE V −1 E UD S ] [0] [Z SA ][Z SS ][Z SI ] [0] [Z IS ][Z II ] 3 7 7 7 5 2 6 6 6 4 ~ u A ~ u S ~ u I 3 7 7 7 5 = 2 6 6 6 4 ~ 0 ~ 0 ~ F I 3 7 7 7 5 (5:29) The above equation is ready for application to three-dimensional models in a semi-infinite medium. Only the boundary conditions pertinent to the model need to be implemented. 5.3 – A Modal Solution The formulation of the previous section requires a great deal of coding to create a finite element model with an exact radiation boundary. With many outstanding finite element programs available on the market, it is preferable to take advantage of these to obtain the necessary solution components instead of duplicating what has already been done. In this section, a new formulation which makes use of extractable information from commercially available software will be implemented. Without any loss of generality, the scheme was tested using the SAP2000 program. All components necessary for the radiation boundary implementation were obtained successfully and straighforwardly. Some technical 89 tricks were developed to extract key information which is usually eliminated during the matrix formation process. For example, all the degrees of freedom of the constraint nodes are automatically eliminated. That information is vital for the formulation of a radiation boundary. The information for these constraint degrees of freedom was extracted using the scheme described in Section 5.2.1. Figure 5.3 – Simplified Modal Analysis Model. Consider the simpler model as shown in Fig. 5.3. As compared to the earlier model in Fig. 5.2, the nodal points in subsetS are eliminated and the stresses will be calculated using modal stress information provided by the finite element program. Numbering the nodes in the sequence, subsetE,A, and thenI, the equation of motion can be written in a partition matrix form: 2 6 6 6 4 [Z EE ][Z EA ] [0] [Z AE ][Z AA ][Z AI ] [0] [Z IA ][Z II ] 3 7 7 7 5 2 6 6 6 4 ~ u E ~ u A ~ u I 3 7 7 7 5 = 2 6 6 6 4 ~ 0 ~ 0 ~ F I 3 7 7 7 5 ; (5:30) 90 SubsetE SubsetI SubsetA The submatrices[Z] include all the effects of the mass, damping, and stiffness of the interior model. The submatrices in Eq. (5.30) are not easy to extract from a standard finite element software program. As restraints are placed on a certain node, the corresponding degrees of freedom are condensed out during the matrix assembly process. More specifically, the outer part of the model would be fixed or partially restrained so that the submatrices[Z EE ], [Z EA ]=[Z AE ] T are eliminated. A simple scheme based on Richardson Extrapolation will be introduced later to recover these submatrices. Condense Eq. (5.30) so that the contribution if~ u E is moved to the right side of the equation. Eq. (5.30) can then be re-written in this form: 2 4 [Z AA ][Z AI ] [Z IA ][Z II ] 3 5 2 4 ~ u A ~ u I 3 5 = 2 4 −[Z AE ]~ u E ~ F I 3 5 : (5:31) The matrix on the left side of the above equation is the same as that of a model which has all the nodes in subsetE fixed. Now perform a modal analysis of this fixed boundary model and generate a sufficient number of modal frequencies and mode shapes for later application. The matrix on the left side of Eq. (5.31) can be optimized for a small bandwidth by the finite element program. Renumbering of the results will therefore be necessary in the later stages. Define the modal matrix, [] as a matrix which contains the mode shapes as columns so that the displacements at the artificial boundary,A, and in the interior of the model,I, can be expressed as linear combinations of the mode shapes: ~ u A ~ u I = []~; (5:32) in which~ contains the modal ordinates. The properties of the modal matrix, [] are such that [] T 2 4 [M AA ] [0] [0] [M II ] 3 5 []=[I] ; (5:33) 91 and [] T 2 4 [K AA ][K AI ] [K IA ][M II ] 3 5 []=[diag(! 2 i )] : (5:34) Substituting Eq. (5.32) into Eq. (5.31) and premultiplying the equation by [] T yields the following equation, [] T 2 4 [K AA ][K AI ] [K IA ][M II ] 3 5 []~ =[D(!)]~ = [] T 2 4 −[Z AE ]~ u E ~ F I 3 5 ; (5:35) in which the dynamic matrix, [D(!)] is defined as [D(!)] = diag(−! 2 +! 2 i +2i!! i i ) : (5:36) In the above equation,! i contains the modal frequencies of the fixed boundary model. The submatrix, [Z AE ], on the right side of Eq. (5.35), can be simplified by assuming that the material damping between nodes in subsetsE andA is negligible. Also, recognizing that SAP2000 uses a lumped mass model, we obtain [Z AE ]=[K AE ], with the latter matrix being the stiffness matrix between nodes of those two subsets. The displacement vector, ~ u E , is the missing displacement that a finite element model cannot produce. This vector is to be calculated using the reciprocal theorem, as stated in Eq. (5.17), with the Green’s Functions of a semi-infinite medium. Let [ A ] contain the mode shape ordinates at the selected nodes of subset A. The displacement vector,~ u A , can then be expressed as [ A ]~ . Define [ S ] as the mode shape ordinates at the selected nodes surrounding subset A. This will provide enough stress information for the calculation of the traction vector, ~ t, through a matrix [N],as ~ t = [N][ S ]~ . Substitution of these two new definitions into Eq. (5.17) yields the vector~ u E in terms of the modal ordinates: ~ u E =(−[T][ A ]+[U][N][ S ])~: (5:37) 92 The right side of Eq. (5.35) can now be separated: [] T 2 4 −[Z AE ]~ u E ~ F I 3 5 = [] T 2 4 ~ 0 ~ F I 3 5 −[ A ] T [K AE ](−[T][ A ]+[U][N][ S ])~ Gathering the term related to ~ to the left side, Eq. (5.35) now has the form [R]~ = [] T 2 4 ~ 0 ~ F I 3 5 ; (5:38) in which [R]=[D(!)]+[ A ] T [K AE ](−[T][ A ]+[U][N][ S ]) : (5:39) The formulation of [R], being a full, complex matrix, guarantees that resonant behavior of the interior medium will not occur as radiation damping terms are added to the interior medium’s dynamic matrix,[D(!)]. The size of matrix[R] is equal to the number of modes used, therefore, it is not unreasonably large. Furthermore, the formation of the mode shape matrices is performed just once and calculation over the frequency domain should not require a great computational effort. The application of the above method to soil-structure interaction is based on the type of forces used in ~ F I and the displacement produced for ~ u I . Many other interesting boundary value problems can be solved when complications from the artificial boundary are eliminated. 5.3.1 – Richardson Extrapolation for the Missing Submatrix One important entry in the matrix[R] in Eq. (5.39) is the submatrix[K AE ]. In a modal analysis model, the nodes in subset E are fixed, therefore, [K AE ] will not be present in the output information. To extract this important submatrix, the nodes in subset E must be 93 Figure 5.4–APartially Restrainted Finite Element Model. free. To ensure the stability of the model, several springs with arbitrary values ofk are attached at various points of the model, as shown in Fig. 5.4. SAP2000 and many other finite element programs have this capability. For the proposed algorithm, the value of is not necessarily small, but the stiffness matrix must be formed twice, once with the spring value ofk and the second time with the spring value of2k. Define the resulting stiffness matrix for the first case to be [K()] and for the second case to be [K(2)]. Then Richardson Extrapolation will generate the matrix without the attached springs, [K(0)]: [K(0)] = 2[K()] −[K(2)] : (5:40) The above formula will efficiently generate the required matrix [K AE ]. At the locations where there are no springs attached, the stiffness matrix element is calculated as 2k −k = k. At the special locations where springs are attached, the stiffness matrix elements are calculated as2(k+k s )−(k+2k s )=k, in whichk s is the arbitrary value of the spring. This algorithm for obtaining [K AE ] completes the extraction of required information for a finite element analysis with an exact radiation boundary. 94 Subset E Subset A k k k k k k Chapter Six Conclusion Three applications of the Reciprocal Theorem for three-dimensional soil-structure interaction analysis of rigid, embedded foundations have been presented. The major requirement for all three applications is that the problem must be linear, otherwise, the principle of superposition cannot be applied. Method One, a volume approach, is a modification of a method introduced by Ohsaki for static results. It uses the body forces within the volume of the foundation to satisfy the displacement compatibility conditions. This method has the advantage of requiring only the displacement Green’s Functions, and therefore, less computer code development. It has proven to be the most viable method, of the three studied, for routine engineering applications. Method One generates stable results and (a) can analyze arbitrarily shaped foundations, (b) is able to handle layered viscoelastic media, (c) can analyze the response of a massless foundation to all types of incident waves for three translations and three rotations, and (d) has the potential to analyze multiple foundation configurations using an iterative algorithm without severely increasing computer memory requirements. The drawback to Method One, the Modified Ohsaki Method, is that it is a volume formulation, therefore, it will require significant computational effort if the foundation is deeply embedded. Future improvement of this method should include non-cubical subregions (only cubical elements were used in this dissertation) to allow for more flexibility when modeling an oddly shaped foundation. Another consideration is that the stress concentration occurs near the rigid foundation’s surface. A scheme could be implemented 95 using smaller subregions near the surface and larger subregions within the interior of the foundation, thereby significantly reducing the number of degrees of freedom for the foundation. Over 70% of the computation time is devoted to the solution of complex, simultaneous equations. Computer time can be greatly diminished by reducing the number of subregions. Through the use of the Boundary Integral Equation created by the reciprocal relationship, Method Two was formulated using integrals over only the surface of the embedded foundation model. Therefore, the spatial dimension of the numerical method is reduced by one compared to Method One. This is a very attractive approach, as a significant reduction of unknowns can be achieved for detailed models. The main problem with this approach, as indicated by many researchers, is that the resonant frequencies of the interior portion of the foundation volume affect the numerical solution of the exterior problem. A simple explanation was given in Chapter Four, based on a two-dimensional Hilbert-Schmidt solution, for why the interior frequencies are involved in an exterior formulation. An algorithm to combat the instability was developed by recognizing that the surface integrals of tractions and displacements, written in matrix form, are normal matrices. As a byproduct of this discovery, it is easy to prove that both matrices share the same eigenvalues and eigenvectors. A scheme was suggested in Chapter Four that substitutes the division of two small numbers by a ratio of slopes, an application of l’Hospital’s Rule for the eigenvalue under consideration. From a practical perspective, Method Two was not explored as an alternative for routine engineering applications in this dissertation. Instead, as an academic contribution, some recommendations were made for improving its stability. As mentioned previously, only surface integrals are required, leading to a computational reduction. This advantage 96 over Method One, however, is lessened because Method Two requires both traction and displacement Green’s Function matrices and is more mathematically difficult. The third application of the Reciprocal Theorem uses the same integral representation as the Boundary Integral Equation Method except the displacement is calculated at points away from the surface of the wave scatterer. Additionally, the surface for this formulation is not the foundation model’s surface, but rather, a surface which encloses the soil medium model around the foundation site. The wave motion inside this surface is modelled by a discrete method such as the Finite Element method, and the integral representation, based on the Reciprocal Theorem, is used to radiate the outgoing waves into the surrounding semi- infinite medium. Chien (2007) had some success using this approach for the simpler SH wave problems and for scattering problems in an infinite medium. The present application is designed for a three-dimensional, semi-infinite medium with more difficult Green’s Functions, not expressible in closed form. To tap into the capabilities of existing Finite Element programs, a scheme was developed using the concepts of Richardson Extrapolation to extract pertinent matrix information for the implementation of Method Three. It has been shown that the normal modes and their accompanying mode frequencies, modal displacements, and modal stresses can be routinely generated as output. The extraction of the stiffness matrix for the nodes at the artificial boundary, however is not simple because those degrees of freedom are eliminated during the assembly process as restraints are placed on those nodes to satisfy the boundary conditions. It was proposed that an adequate number of springs be attached to an otherwise unrestrained structural model and to then output the corresponding stiffness matrix. Then repeat the process by doubling the value of the springs and output this second stiffness matrix. Using these two matrices, an extrapolation can be made from double- valued springs to single-valued springs to a spring-free model. After this simple process, 97 the influence of the springs is totally removed. This algorithm is applicable regardless of the values used for the springs. In conclusion, the three applications of the Reciprocal Theorem have resulted in (a) the development of a practical method for calculating impedance functions and driving forces needed for a soil-structure interaction analysis, (b) the recommendation of a scheme to improve the stability of the Boundary Integral Method, and (c) a method to extract information from commercially available finite element programs for the implementation of an exact radiation boundary. All three methods provide a basis for future development. 98 References Abramowitz, M.A., and I.A. Stegun (1970).“Handbook of Mathematical Functions,” Dover 0-486-61272-4,New York. Akino, K., Ohtsuka, Y ., Fukuoka, A., Ishida, K (1996). “Experimental studies on Embedment Effects On Dynamic Soil-Structure Interaction,” 11th World Conference on Earthquake Engineering, Paper No. 59. Apsel, R.J. (1979). “Dynamic Green’s Functions for Layered Media and Applications to Boundary-Value Problems,” Ph.D. Dissertation, University of California, San Diego, La Jolla, California. Brinkgreve, R.B.J. and Swolfs W.M. (2008). “Possibilities and Limitations of the Finite Element Method for Geotechnical Applications,” 4th Conference on Advances and Applications of GiD, Ibiza, Spain. Burton, A.J. and Miller, G.F. (1971). “The application of integral equation methods to the numerical solution of some exterior boundary-value problems, Proceedings of the Royal Society of London V ol. 323, 201-210. Celebi, M. and Crouse, C.B. (2001) “Recommendations for Soil Structure Interaction (SSI) Instrumentation,” COSMOS Workshop on Structural Instrumentation, Emeryville, California. Chien, A.Y . (2007). “A Complete Time-Harmonic Radiation Boundary for Discrete Elastodynamic Models,” Ph.D. Dissertation, University of Southern California, Los Angeles, California. Choi,J.-S., Yun, C.-B., Kim, J.-M. (2001). “Earthquake response analysis of the Hualien soil-structure interaction system based on updated soil prperties using forced vibration test data,” Earthquake Eng. and Str. Dyn. V ol. 30, 1-26. Day, S.M. (1977). “Finite Element Analysis of Seismic Scattering Problems,” Ph.D. Thesis, University of California, San Diego, Ca. Duan, X. (1999). “A Soil Structure-Interaction Analysis of Tall Buildings,” Ph.D. Dissertation, University of Southern California, Los Angeles, California. Ewers, J., Hitzschke, U., Krutzik, N.J., Papandreou, D., and Schutz, W. (1993). “Time Versus Frequency Domain Analysis of Nuclear Power Plant Building Structures,” 12th International Conference on Structural Mechanics in Reactor Technology, Stuttgart, Germany, 307-312. Ewing, W.M., W.S. Jardetsky, and F. Press (1957).“Elastic Waves in Layered Media,” McGraw-Hill Book Company, Inc., New York. Frank, R. (2006). ”Some aspects of soil-structure interaction according to Eurocode 7 ’Geotechnical design’,” Journal of Civil Engineering of the University of the Minho, Number 25, http://www.civil.uminho.pt/cec/revista/Num25/n 25 pag 5-16.pdf. 99 Fujimori, T., Tsundoa, T., Izumi, M., and Akino, K. (1992). “Partial Embedment Effects on Soil-Structure Interaction,” Proc. 10th World Conf. on Earthquake Eng., V ol. 3, 1713-1718. Fung, Y .C. (1965). “Foundation of Solid Mechanics,” Prentice-Hall, ISBN-013329920. Gueguen, P., Bard, P.-Y ., Chavez-Garcia, F.J. (2002). “Site-City Seismic Interaction in Mexico City-Like Environments: An Analytical Study,” Bulletin of the Seismological Society of America, V ol. 92, No.2, 794-811. Iguchi, M. ,Akin, K., Jido, J.I., Kawamura, S., Ishikawa, Y ., Nakata, M. (1988). “Large- Scale Model Tests on Soil-Reactor Building Interaction Part I: Forced Vibration Tests,” Proc. Of Ninth World Conf. on Earthquake Eng., V ol. 3, 697-702. Inukai, T., Imazawa, T., Kusakabe, K., Yamamoto, M. (1992). “Dynamic behavior of embedded structure on hard rock site,” Proc. 10th World Conf. on Earthquake Eng., V ol. 3, 1695-1699. Imamura, A., Watanabe, T., Ishizaki, M., and Motosaka, M. (1992). “Seismic response characteristics of embedded structures considering cross interaction,” Proc. 10th World Conf. on Earthquake Eng., V ol. 3, 1719-1724. Jennings, P.C. and Bielak, J. (1973). “Dynamics of Building-Soil Interaction,” Bulletin of the Seismological Society of America, V ol. 63, No.1, 9-48. Konno, T., Motohashi, S., Izumi, M., and Iizuka, S. (1993). “Study on Vertical Seismic Response Model of BWR-Type Reactor Building,” 12th International Conference on Structural Mechanics in Reactor Technology, Stuttgart, Germany, 217-222. Kurimoto, O., Seki, T., and Omote, Y . (1992). “Dynamic characteristics of irregularly embedded foundation,” Proc. 10th World Conf. on Earthquake Eng., V ol. 3, 1725-1730. Lamb, H. (1903). “On the Propagation of Tremors Over the Surface of an Elastic Solid,” Rhil. Transaction of Royal Rhil. Soc. Of Lond., 1-4-2. Liu, T. and Chen, S. (1999). “A new form of the hypersingular boundary integral equation for 3-D acoustics and its implementaion withC O boundary elements,” Computer Methods in Applied Mechanical Engineering, V ol. 173, 375-386. Love, A.E.H. (1926). “A Treatise on the Mathematical Theory of Elasticity,” ISBN 0-486- 60174-9, Dover Publications, New York. Luco, J.E. (1976). “Torsional response of structures for SH waves: The case of hemispherical foundations,’ Bulletin of Seismological Society of America, V ol. 66, No. 1, 109-123. Luco, J.E. and R.J. Apsel (1983). “On the Green’s Functions for a Layered Half-Space, Part I,” Bulletin of Seismological Society of America, V ol. 73, No. 4, 907-929. Luco, J.E. and R.J. Apsel (1983). “On the Green’s Functions for a Layered Half-Space, Part II,” Bulletin of Seismological Society of America, V ol. 73, No. 4, 931-951. 100 Lysmer, J. (1970). “Lumped mass method for Rayleigh waves,” Bulletin of the Seismological Socity of America, V ol. 60, No. 1, 89-104. Meek, J.W. and Wolf, J.P. (1994). “Cone Models for Embedded Foundation,” J. of Geotechnical Engineering, V ol. 120, No. 1, 60-79. Mindlin, R. (1936). “Force at a Point in the Interior of a Semi-Infinite Solid,” Phys.7, 195–202. Mizuhata, K., Kusakabe, K., and Shirase, Y . (1988). “Study on Dynamic Characteristics of Embedded Mass and its Surrounding Ground,” Proc. Of Ninth World Conf. on Earthquake Eng., V ol. 3, 679-685. Mylonakis, G. and Gazetas, G. (2000). “Seismic Soil-Structure Interaction: Beneficial or Detrimental?” Journal of Earthquake Engineering, V ol. 4, No. 3, pp. 277-301. Ohsaki, Y . (1973). “On Movements of a Rigid Body in Semi-Infinite Elastic Medium,” Proc. of the Japan Earthquake Engineering Symposium, Tokyo, 245-252. Pao, Y .H. and C.C. Mow (1973). “Diffraction of Elastic Waves and Dynamic Stress Concentration,” Crane-Russak, ISBN-9780844801551. Parmlee, R.A. and Kudder, R.J. (1973). “Seismic Soil-Structure Interaction of Embedded Buildings,” Proc. Fifth World Conf. on Earthquake Eng., Rome, Italy, 1941-1950. Pyl, L., Clouteau, D., Degrande, G. (2003). “The Soil Impedance of Embedded Structures in the Higher Frequency Range, Computed With the Boundary Element Method,” 16th ASCE Engineering Mechanics Conference, Seattle, Washington. Rizzo, F.J., (1989). “The Boundary Element Method: Some Early History–A Personal View,” Boundary Element Methods in Structural Analysis, Ed. D.E. Besksos, ASCE, 1-16. Sandi, H. (1960). “A Theoretical Investigation of the Interaction Between Ground and Structure During Earthquakes,” Proc. Second World Conf. on Earthquake Eng., 1327- 1343. Seed, H.B. and Idriss, I.M. (1973). “Soil-Structure Interaction of Massive Embedded Structures During Earthquakes,” Proc. Fifth World Conf. on Earthquake Eng., 1881-1890. Shippy, D.J. (2003). “Early Development of the BEM at the University of Kentucky,” Electronic Journal of Boundary Elements, V ol. 1, 26-33. Spyrakos, C.C., and Xu, C. (2004). “Dynamic analysis of flexible massive strip-foundations embedded in layered soils by hybrid BEM-FEM,” Computers and Structures, V ol. 82, 2541- 2550. Stewart, J.P, Comartin, C., and Moehl, J.P. (2004). “Implementation of Soil-Structure Interaction Models in Performance Based Design Procedures,” Proc. 3rd UJNR Workshop on Soil-Structure Interaction, Contributions 4. 101 Takewaki,I.,Takeda,N.,Uetani, K. (2003). “Fast,practical evaluation of soil-structure interaction of embedded structures,” Soil Dyn. and Earthquake Eng., V ol. 23, 13-20. Tassoulas, J.L. (1989). “Dynamic Soil-Structure Interaction,” Boundary Element Methods in Structural Analysis, Ed. D.E. Besksos, ASCE, 273-308. Trifunac, M.D. (2005). “Power Design Method,” Proc. Earthquake Engineering in the 21st Century, Skopje-Ohrid, Macedonia. Urao, K., Masuda, K., Kitamura, E., Saski, F., Ueno, K., Miyamoto, Y ., and Moroi, T. (1988). “Force Vibration Test and its Analytical Study for Embedded Foundation Supported by Pile Group,” Proc. Of Ninth World Conf. on Earthquake Eng., V ol. 3, 673-678. Watakabe, M., Matsumoto, H., Ariizumi, K., Fukahori, Y ., Shikama, Y ., Yamanouchi, K., Kuniyoshi, H. (1992). “Earthquake observation of deeply embedded building structure,” Proc. 10th World Conf. on Earthquake Eng., V ol. 3, 1831-1833. Wong, H.L. (1975). “Dynamic Soil-Structure Interaction,” Ph.D. Dissertation, California Institute of Technology, Pasadena, California. Wong, H.L. and Trifunac, M.D. (1988). “A Comparison of Soil-Structure Interaction Calculations with Results of Full-Scale Forced Vibration Tests,” Soil Dynamics and Earthquake Enigneering, V ol.7, 22-31. http://www.engin.brown.edu/courses/En222/Notes/potenergy/potenergy.htm http://pubs.usgs.gov/of/1996/ofr-96-0263/damagebe.htm 102 Appendix A USGS Damage Report of the Olive View Hospital Building This is an excerpt from the 1996 USGS damage assessment report for the 1994 Northridge Earthquake. The information provided is relevant to the material presented in this thesis. The following paragraphs review the design and construction of the Olive View Hospital as well as the damage it sustained in the 1994 Northridge Earthquake. The values of the peak displacements and accelerations are given, together with an explanation for the damage that was done. Olive View Hospital (OVH) Damage Report The OVH building was designed in 1976 to withstand increased levels of seismic forces based on the disastrous fate of its predecessor. The structural system for resisting lateral forces is a mixed design of concrete and steel shear walls. The foundation consists of spread footings and concrete slab-on-grade for the ground floor. The ground floor and second floor are typically built of concrete shear walls 25 centimeters thick that extend along several column lines. At the third level, the plan of the building changes to a cross shape, making a 4-story tower with steel shear walls surrounding the perimeter. USGS engineers examined the building performance by using data from both the Northridge and Whittier Narrows earthquakes for comparisons. The building was designed for two levels of performance. The first level, 0.52g, represented accelerations at which the building would not be badly damaged. The second level, 0.69g, represented accelerations at which the building would survive, perhaps with major damage but without catastrophic failure. 103 Data from sensors in the OVH building show that the building escaped the impact of long-period (>1 second) pulses generated by the Northridge quake. The OVH records indicated that very large peak accelerations at the ground level (0.91g) and the roof level (2.31g) were accommodated by the building without structural damage. Analyses of the data indicate that the structure was in resonance at frequencies between 2.5 and 3.3 hertz. These frequencies are also within the site-response frequencies of 2-3 hertz, calculated from examining the geological materials upon which the structure was built. The effective structural frequencies derived from the Northridge and Whittier earthquakes data are different and exhibit variations attributable to nonlinear effects. One effect is soil- structure interaction which was seen to be more pronounced during the Northridge event. The building possibly experienced rocking at 2.5 hertz in the north-south direction during that event, and there is the possibility that radiation damping (wherein the building dissipates energy into the surrounding soil) contributed to reducing that response. Damping ratios for the building were 10%-15% (north-south) and 5%-10% (east-west) for the Northridge effects, and 1%-4% (north-south) and 5%-8% (east west) for the Whittier effects. These nonlinear effects that tended to reduce the shaking of the building during large ground motions are consistent with the different damping ratios observed during the Northridge and Whittier Narrows events. There was also nonlinear behavior due to minor structural damage during the Northridge earthquake. It is also likely that the cruciform wings responded with a different frequency than that of the overall building. The Olive View Hospital building was conceived and designed as a very strong and stiff structure, particularly in response to the disastrous performance of the original OVH building during the 1971 San Fernando earthquake. However, the resulting design placed the fundamental frequency of the building within the frequency range of the site (2-3 hertz), thereby producing conditions for resonance. This case study indicates that determining the 104 site frequencies needs to be emphasized in developing design response spectra. Despite the site resonances, the performance of the OVH building was considered to be a great success during the Northridge event. The hospital sustained no structural damage under very strong shaking (greater than 2g) at the roof level. 105 Appendix B Example of the Reciprocal Theorem The reciprocal theorem is based on a physical principle of work-energy, it is therefore not a mathematical principle. It may be difficult to grasp how it works, therefore a simple example will be shown here to verify, numerically at least, that it does indeed hold true. The solution of the elastic wave equation in an infinite space will be used as the basis for the verification. Shown in later sections in this appendix will be the formulas for displacements generated by point loads inx,y andz directions. Let the two locations be represented by~ r=(x;y;z) coordinates as ~ r 1 = 2 4 2 2 0:5 3 5 and ~ r 2 = 2 4 1 7 4 3 5 : and the point forces ~ P =(P x ;P y ;P z ) be given as at the respected locations as ~ P 1 = 2 4 1 −3 −3 3 5 and ~ P 2 = 2 4 0 5 −4 3 5 : The resulting displacement vector at the opposite location are computed as ~ u 1 (~ r 2 )= 2 4 0:0396+0:0052i −0:1734−0:0192i 0:1324−0:0164i 3 5 and ~ u 2 (~ r 1 )= 2 4 −0:0114−0:0008i 0:1178+0:0204i −0:0090−0:0105i 3 5 : The two set of inner products are the same as shown by ~ P T 2 ~ u 1 (~ r 1 )=−0:33768−0:030367i; and ~ P T 1 ~ u 2 (~ r 2 )=−0:33768−0:030367i: 106 B.1 Displacements and Stresses Generated by a Vertical Point Load Shown in Fig. B.1 is a concentrated loadQ in the positivez-direction. The solution of the three-dimensional wave equation can be written in terms of two potential functions and as shown in Lamb’s paper (1904) as = −Q 4k 2 @ @z e −ihr r ; (B:1) = Q 4k 2 e −ikr r ; (B:2) in whichr = p x 2 +y 2 +z 2 ,x =x o −x s ,y =y o −y s andz =z 0 −z s are the relative position of the observation point with respect to the source point in thex,y andz-directions, respectively. Figure B.1 – VerticalQ Z Point Force Configuration. 107 x y Q Z z The argument of , hr = !r=, is a dimensionless frequency normalized by the compressional wave velocity, implying that is a potential for compressional waves. , on the other hand, is the shear wave potential becausekr = !r= is normalized by the shear wave velocity. The exponential function with a negative argument is used in this derivation because it represents an outgoing wave as r!1 when associated with the harmonic time factore i!t . For a unit applied loadQ in the vertical (Z) direction, let the amplitudeQ=1. Also, define for convenience the parameters γ k = 1 4k 2 ; (B:3) " h = e −ihr r ; (B:4) and " k = e −ikr r ; (B:5) so that equations (B.1) and (B.2) can be written as =−γ k @" h @z and =γ k " k : B.1.1 – Displacements Using the potentials and, the displacements for a vertical load in thex,y, andz- directions, respectively, can be expressed as, U 31 = @ @x + @ 2 @x@z (B:6) 108 = γ k @ 2 @x@z " k −" h U 32 = @ @y + @ 2 @y@z (B:7) = γ k @ 2 @y@z " k −" h and U 33 = @ @z + @ 2 @z 2 +k 2 (B:8) = γ k @ 2 @z 2 " k −" h +k 2 " k in which, the derivatives of" k are represented by @ 2 " k @x@z = xz r 3 3 r 1 r +ik −k 2 e −ikr ; (B:9a) @ 2 " k @y@z = yz r 3 3 r 1 r +ik −k 2 e −ikr ; (B:9b) and @ 2 " k @z 2 = 1 r 3 (3z 2 −r 2 ) r 1 r +ik −k 2 z 2 e −ikr ; (B:9c) The derivatives of" h are in the same form except the subscript and variablek should be replaced byh. B.2 Displacements and Stresses Generated by a Point Load in theX-Direction Shown in Fig. B.2 is an illustration of the horizontal point load in thex-direction using the(x;y;z) coordinate system. Also in the same figure is the(x 0 ;y 0 ;z 0 ) coordinate system. 109 The (x 0 ;y 0 ;z 0 ) system is rotated from the (x;y;z) system by 90 about the y-axis. The orthogonal transformation between these systems can be expressed as [Q(xyz − x 0 y 0 z 0 )] = 2 4 001 100 010 3 5 : To obtain the results for the horizontal point loadQ X , two steps are required: Figure B.2 – HorizontalQ X Point Force Configuration. (1) Calculate the displacements and stresses in the(x 0 ;y 0 ;z 0 ) system using the expressions presented in Section B.1. This is done becauseQ X is in thez 0 -direction, formerly the vertical direction. The values for the prime coordinates can be obtained from the position vector of the present configuration as 2 4 x 0 y 0 z 0 3 5 =[Q] T 2 4 x y z 3 5 = 2 4 010 001 100 3 5 2 4 x y z 3 5 = 2 4 y z x 3 5 : 110 x,z y,x z,y Q X (2) Use the calculated displacements in the (x 0 ;y 0 ;z 0 ) system and transform them to the (x;y;z) system using 2 4 u x u y u z 3 5 =[Q] 2 4 u x 0 u y 0 u z 0 3 5 = 2 4 001 100 010 3 5 2 4 u x 0 u y 0 u z 0 3 5 = 2 4 u z 0 u x 0 u y 0 3 5 : B.2.1 – Displacements If the displacements in Section B.1 can be written in functional form as 2 4 U 31 (x 0 ;y 0 ;z 0 ) U 32 (x 0 ;y 0 ;z 0 ) U 33 (x 0 ;y 0 ;z 0 ) 3 5 then the displacements caused by a horizontal point load in thex-direction can be written as 2 4 U 11 (x;y;z) U 12 (x;y;z) U 13 (x;y;z) 3 5 = 2 4 U 33 (y;z;x) U 31 (y;z;x) U 32 (y;z;x) 3 5 (B:10) in which the coordinatesx,y andz are the only parameters displayed because the material properties of the viscoelastic medium remain unchanged. B.3 Displacements and Stresses Generated by a Point Load in theY -Direction Shown in Fig. B.3 is an illustration of the horizontal point load in they-direction using the(x;y;z) coordinate system. Also in the same figure is the(x 0 ;y 0 ;z 0 ) coordinate system. The (x 0 ;y 0 ;z 0 ) system is rotated from the (x;y;z) system by 90 about thex-axis. The orthogonal transformation between these systems can be expressed as [Q(xyz − x 0 y 0 z 0 )] = 2 4 010 001 100 3 5 : 111 Figure B.3 – HorizontalQ Y Point Force Configuration. To obtain the results for the horizontal point loadQ Y , the same type of procedures as those used in Section B.2 can be applied. First calculate the displacements and stresses in the (x 0 ;y 0 ;z 0 ) system using the expressions presented in Section B.1. This is done becauseQ Y is in thez 0 -direction, formerly the vertical direction. The values for the prime coordinates can be obtained from the position vector of the present configuration as 2 4 x 0 y 0 z 0 3 5 =[Q] T 2 4 x y z 3 5 = 2 4 001 100 010 3 5 2 4 x y z 3 5 = 2 4 z x y 3 5 : Next use the calculated displacements in the (x 0 ;y 0 ;z 0 ) system and transform them to the (x;y;z) system using 2 4 u x u y u z 3 5 =[Q] 2 4 u x 0 u y 0 u z 0 3 5 = 2 4 010 001 100 3 5 2 4 u x 0 u y 0 u z 0 3 5 = 2 4 u y 0 u z 0 u x 0 3 5 : 112 x,y y,z Q Y z,x B.3.1 – Displacements In terms of the displacements in Section B.1, the displacements caused by a horizontal point load in they-direction can be written as 2 4 U 21 (x;y;z) U 22 (x;y;z) U 23 (x;y;z) 3 5 = 2 4 U 32 (z;x;y) U 33 (z;x;y) U 31 (z;x;y) 3 5 (B:11) 113 Appendix C V olume Integrals of Static Singularities The Green’s Functions for all three-dimensional configurations have the same type of singularities, regardless if the problem is dynamic or static, full space or semi-infinite space. This appendix will provide the exact integrals for the static full-space Green’s Functions (Love, 1906) over a rectangular volume. These integrals will be used to remove the singularities from numerical integration in Chapter 3. For an isotropic infinite solid, the displacements caused by a vertical loadP (Love, 1906) are 2 4 u x u y u z 3 5 = 1 8 0 @ (1+γ 2 ) 2 4 0 0 P=r 3 5 +(1−γ 2 ) P(z−z 0 ) r 2 2 4 (x−x 0 )=r (y−y 0 )=r (z−z 0 )=r 3 5 1 A : (C:1) Similarly, the displacements caused by a point load in thex-direction are 2 4 u x u y u z 3 5 = 1 8 0 @ (1+γ 2 ) 2 4 P=r 0 0 3 5 +(1−γ 2 ) P(x−x 0 ) r 2 2 4 (x−x 0 )=r (y−y 0 )=r (z−z 0 )=r 3 5 1 A ; (C:2) and the displacements caused by a point load in they-direction are 2 4 u x u y u z 3 5 = 1 8 0 @ (1+γ 2 ) 2 4 0 P=r 0 3 5 +(1−γ 2 ) P(y−y 0 ) r 2 2 4 (x−x 0 )=r (y−y 0 )=r (z−z 0 )=r 3 5 1 A : (C:3) To account for all types of singularities, the integration of functions of the type: 1=r, (x−x 0 ) 2 =r 3 , (y−y 0 ) 2 =r 3 , (z−z 0 ) 2 =r 3 , (x−x 0 )(y−y 0 )=r 3 , (x−x 0 )(z−z 0 )=r 3 and (y−y 0 )(z−z 0 )=r 3 , in which r = p (x−x 0 ) 2 +(y−y 0 ) 2 +(z−z 0 ) 2 , will be performed. 114 Volume Integrals of singularities Consider first the volume integral of1=r, designatedS 0 in Chapter 3. The more difficult form S 0 = Z a −a Z a −a Z a −a dxdydz p (x−x 0 ) 2 +(y−y 0 ) 2 +(z−z 0 ) 2 ; (C:4) can be simplified using a change of variables: X 1 =−a−x 0 ; X 2 =a−x 0 ; Y 1 =−a−y 0 ; Y 2 =a−y 0 ; Z 1 =−a−z 0 ; Z 2 =a−z 0 ; Eq. (C.3) can then be expressed as S 0 = Z Z 2 Z 1 Z Y 2 Y 1 Z X 2 X 1 dxdydz p x 2 +y 2 +z 2 ; (C:5) a more manageable form. The analytical integral obtained by Mathematica generated 48 terms, but by defining the functions, q 1 (d;e;f)=− f 2 2 tan −1 de f p d 2 +e 2 +f 2 ! ; (C:6) q 2 (d;e;f 1 ;f 2 )=de log f 2 p d 2 +e 2 +f 2 2 f 1 p d 2 +e 2 +f 2 1 ; (C:7) S 0 can be expressed as S 0 = 2 X i=1 2 X j=1 2 X k=1 (−1) i+j+k [q 1 (Y j ;Z k ;X i )+q 1 (X i ;Z k ;Y j )+q 1 (X i ;Y j ;Z k )] + 2 X i=1 2 X j=1 (−1) i+j [q 2 (X i ;Y j ;Z 1 ;Z 2 )+q 2 (Y i ;Z j ;X 1 ;X 2 )+q 2 (X i ;Z j ;Y 1 ;Y 2 )] : (C:8) Continuing as before, the integralS zz has a similar form: S zz = Z a −a Z a −a Z a −a (z−z 0 ) 2 dxdydz (x−x 0 ) 2 +(y−y 0 ) 2 +(z−z 0 ) 2 3=2 ; (C:9) 115 and its solution can be expressed as S zz = 2 X i=1 2 X j=1 2 X k=1 (−1) i+j+k [q 1 (Y j ;Z k ;X i )+q 1 (X i ;Z k ;Y j )−q 1 (X i ;Y j ;Z k )] + 2 X i=1 2 X j=1 (−1) i+j q 2 (X i ;Y j ;Z 1 ;Z 2 ) : (C:10) Using transformation matrices, other integrals can be obtained by rearranging the independent variables as S xx (x;y;z)=S zz (z;x;y) ; (C:11) S yy (x;y;z)=S zz (y;z;x) : (C:12) The off-diagonal terms have a simpler form, S xz = Z a −a Z a −a Z a −a (x−x 0 )(z−z 0 )dxdydz (x−x 0 ) 2 +(y−y 0 ) 2 +(z−z 0 ) 2 3=2 ; (C:13) with an algebraic solution of the form S xz =− 1 2 2 X i=1 2 X j=1 2 X k=1 (−1) i+j+k y j q X 2 i +Y 2 j +Z 2 k +(X 2 i +Z 2 k )log(Y j + q X 2 i +Y 2 j +Z 2 k : (C:14) Similarly, S yz = Z a −a Z a −a Z a −a (y−y 0 )(z−z 0 )dxdydz (x−x 0 ) 2 +(y−y 0 ) 2 +(z−z 0 ) 2 3=2 ; (C:15) S xy = Z a −a Z a −a Z a −a (x−x 0 )(y−y 0 )dxdydz (x−x 0 ) 2 +(y−y 0 ) 2 +(z−z 0 ) 2 3=2 ; (C:16) 116 have solutions of the form, S yz =− 1 2 2 X i=1 2 X j=1 2 X k=1 (−1) i+j+k X i q X 2 i +Y 2 j +Z 2 k +(Y 2 j +Z 2 k )log(X i + q X 2 i +Y 2 j +Z 2 k ; (C:17) and S xy =− 1 2 2 X i=1 2 X j=1 2 X k=1 (−1) i+j+k Z k q X 2 i +Y 2 j +Z 2 k +(X 2 i +Y 2 j )log(Z k + q X 2 i +Y 2 j +Z 2 k : (C:18) 117 Appendix D Wave Motion in a Semi-Infinite Medium D.1 – Two-Dimensional Plane-Strain Problem A general three-dimensional wave field can be represented by the superpostion of plane waves, therefore, it is a useful exercise to prepare an expression for several fundamental waves. The displacements in elastic wave propagation can be expressed in terms of derivatives of two potential functions, and ~ . The compressional wave potential,, is a scalar but the shear wave potential, ~ , is a vector. For the two-dimensional plane strain wave problem under consideration, the z-component of the vector, z , is used and will be denoted as the scalar symbol herein. The in-plane horizontal displacement, u, and the vertical displacement,v, can be written as (Ewing et al, 1957), u = @ @x + @ @y ; (D:1a) and v = @ @y − @ @x ; (D:1b) respectively. Defining the LaPlacian operator in two dimensions as r 2 = @ 2 @x 2 + @ 2 @y 2 ; (D:2) the two-dimensional stresses are then xx =r 2 +2 @ 2 @x 2 +2 @ 2 @x@y (D:3a) =(+2)r 2 −2 @ 2 @y 2 +2 @ 2 @x@y ; 118 yy =r 2 +2 @ 2 @y 2 −2 @ 2 @x@y (D:3b) =(+2)r 2 −2 @ 2 @x 2 −2 @ 2 @y@x ; and xy = yx =2 @ 2 @x@y + @ 2 @y 2 − @ 2 @x 2 : (D:3c) D.2 – Incident P-Wave Solution Let i be the incident P-wave potential function, r be the reflected P-wave potential function and r be the reflected SV-wave potential function. An additional shear wave reflection, represented by r , is also necessary to satisfy both normal and shear stress boundary conditions at the free surface as shown in Fig. D.1. Figure D.1 – Incident and reflected P-wave. 119 i r r e e f x y Define the potentials as i =Ae −ih(xcos e−ysin e) ; (D:4a) r =Be −ih(xcos e+ysin e) ; (D:4b) r =Ce −ik(xcos f+ysin f) ; (D:4c) withh =!= andk =!= defined as the wave numbers. The second partial derivatives of the potential functions can be written in the form @ 2 i @x 2 =(−ihcose) 2 i ; (D:5a) @ 2 r @x 2 =(−ihcose) 2 r ; (D:5b) @ 2 r @x 2 =(−ikcosf) 2 r ; (D:5c) @ 2 i @y 2 =(ihsine) 2 i ; (D:5d) @ 2 r @y 2 =(−ihsine) 2 r ; (D:5e) @ 2 r @y 2 =(−iksinf) 2 r ; (D:5f) @ 2 i @x@y =(−ihcose)(ihsine) i ; (D:5g) @ 2 r @x@y =(−ihcose)(−ihsine) r ; (D:5h) @ 2 r @x@y =(−ikcosf)(−iksinf) r : (D:5i) Now apply the boundary conditions, yy j y=0 =0 and yx j y=0 =0. The first boundary condition yields, yy j y=0 =(−h 2 cos 2 e)[A+B]e −ihxcos e +(+2)(−h 2 sin 2 e)[A+B]e −ihxcos e −2(−k 2 sinf cosf)Ce −ikxcos f ; (D:6) 120 and for the second boundary condition, yy j y=0 =0, to be true over the entire range ofx, the following requirement, −ihxcose =−ikxcosf; (D:7) is necessary. The above relationship defines the angle,f, based on the incident angle,e. The first boundary condition yields the following expression for the unknown coefficients, [(+2)(−h 2 )+2h 2 cos 2 e]A (D:8) +[(+2)(−h 2 )+2h 2 cos 2 e]B+[2k 2 sinf cosf]C=0: The second boundary condition yields the relationship, yx j y=0 =2(h 2 cosesine)[A−B]e −ihxcos e +(−k 2 sin 2 f)Ce −ikxcos f −(−k 2 cos 2 f)Ce −ikxcos f ; (D:9) and with the relationship,−ihxcose =−ikxcosf, the second simultaneous equation for the unknown coefficients is [2h 2 cosesine]A−[2h 2 cosesine]B +[k 2 (cos 2 f−sin 2 f)]C=0: (D:10) The simultaneous equations for unknownsB=A andC=A can be written in matrix form as (+2)(−h 2 )+2h 2 cos 2 e 2k 2 sinf cosf 2h 2 cosesine −k 2 (cos 2 f−sin 2 f) B=A C=A = (+2)h 2 −2h 2 cos 2 e 2h 2 cosesine : (D:11) Define now the important material constant, γ 2 = h 2 k 2 = ! 2 = 2 ! 2 = 2 = 2 2 = 2 ; (D:12) 121 which can be related to the Lame constant as γ 2 = 2 2 = = (+2)= = +2 : (D:13) With these new definitions, the matrix equation simplifies to −γ 2 +2γ 4 cos 2 e 2γ 2 sinf cosf 2γ 2 cosesine 1−2cos 2 f B=A C=A = γ 2 −2γ 4 cos 2 e 2γ 2 cosesine : (D:14) The determinant of the matrix on the left side is =(−γ 2 +2γ 4 cos 2 e)(1−2cos 2 f)−4γ 4 sinf cosf cosesine =−γ 2 (1−2γ 2 cos 2 e)(1−2cos 2 f)−4γ 4 sinf cosf cosesine : (D:15) More simplification can be accomplished by recognizing thathcose = kcosf, using the relationship, cosf = h k cose =γcose; (D:16) The determinant can be written in a better form: =−γ 2 (1−2γ 2 cos 2 e) 2 −4γ 4 sinf cosf cosesine: (D:17) Using Cramer’s Rule for the matrix equation, the solution for the unknown coefficients is B A = 1 γ 2 (1−2γ 2 cos 2 e) 2 2γ 2 sinf cosf 2γ 2 cosesine 1−2cos 2 f = 1 γ 2 (1−2γ 2 cos 2 e) 2 −4γ 4 sinf cosf cosesine ; (D:18) and C A = 1 −γ 2 (1−2γ 2 cos 2 e) 2 γ 2 (1−2γ 2 cos 2 e) 2 2γ 2 cosesine2γ 2 cosesine = 1 −4γ 4 (1−2γ 2 cos 2 e)cosesine : (D:19) 122 The above expressions are amenable to numerical computation, but a more symmetrical expression can be obtained using the expression cosf = γcose. Manipulation by trigonometric identities yields, sin 2 f =1−cos 2 f =1−γ 2 cos 2 e; (D:20) tan 2 f = sin 2 f cos 2 f = 1−γ 2 cos 2 e γ 2 cos 2 e = sin 2 e+cos 2 e−γ 2 cos 2 e γ 2 cos 2 e = 1 γ 2 tan 2 e+ (1−γ 2 ) γ 2 ; (D:21) and sinf cosf = tanf cos 2 f =γ 2 tanf cos 2 e: (D:22) With the angle e given as the input parameter, the first step is to determine tanf using Eq. (D.21) and then calculate the coefficientsB=A andC=A using the expressions B A = −(1−2γ 2 cos 2 e) 2 +4γ 4 tanf sinecos 3 e (1−2γ 2 cos 2 e) 2 +4γ 4 tanf sinecos 3 e ; (D:23) and C A = 4γ 2 (1−2γ 2 cos 2 e)cosesine (1−2γ 2 cos 2 e) 2 +4γ 4 tanf sinecos 3 e : (D:24) The value of A can be determined from the amplitude of the incident P-wave,jPj,of a compressional displacement pulse. The horizontal and vertical displacements of the incident pulse are u i = @ i @x =(−ihcose) i ; (D:25a) and v i = @ i @y =(ihsine) i ; (D:25b) respectively. Now compute the magnitude,jPj, from the displacements, jPj = p (u i ) 2 +(v i ) 2 = p i 2 h 2 cos 2 eA 2 +i 2 h 2 sin 2 eA 2 =ihA ; (D:26) 123 and in turn the coefficientA can be specified as A = jPj ih : (D:27) D.2.1 – Wavefields for Incident P-Waves With the unknown coefficients determined, the horizontal displacement for allx and fory 0,is u= @ @x + @ @y =−ihcose i + r −iksinf r : (D:28) Using the relationships kγ = ! = ! =h; (D:29) and ksinf =ktanf cosf = tanfkγcose =htanf cose; (D:30) the horizontal displacement within the semi-infinite medium can be written as u(x;y)=−ihcose i (x;y)+ r (x;y)+tanf r (x;y) : (D:31) and the horizontal displacement at the free surface is simply, ju(x;0)j =P cose 1+ B A + C A tanf : (D:32) The vertical displacement for allx andy 0 is, v = @ @y − @ @x =ihsine i − r −(−ikcosf) r : (D:33) Using the relationship between the angles, kcosf =kγcose = ! cose =hcose; (D:34) 124 then v(x;y)=ihsine i (x;y)− r (x;y) +ihcose r (x;y) : (D:35) The vertical displacement at the free surface of the semi-infinite medium is jv(x;0)j =P sine 1− B A + C A cose : (D:36) Since the above expressions are strongly influenced by the ratio of shear to compressional wave speeds,γ, it is convenient to relateγ to the Poisson’s Ratio, in this manner: γ 2 = (1−2) 2(1−) : (D:37) Using=1=3,γ 2 =1=4, the normalized amplitudes for a unit, incident P-wave are eu ff v ff 30 1.39 1.12 60 0.96 1.74 90 0.00 2.00 D.3 – Incident SV-Wave Problem Let i be the incident SV-wave potential function with a yet to be defined amplitude D, r be the reflected SV-wave potential function and r be the reflected P-wave potential function written in the form, i =De −ik(xcos f−ysin f) ; (D:38a) r =Ee −ik(xcos f+ysin f) ; (D:38b) and 125 Figure D.2 – Incident and reflected SV-wave. r =Fe −ih(xcos e+ysin e) ; (D:38c) in whichh =!= andk =!= are the wave numbers. The unknown coefficients,E and F are to be determined later using the boundary conditions. Similar to the incident P-wave problem, r is necessary to satisfy both boundary conditions at the free surface. The second partial derivatives of the potential functions can be listed as @ 2 i @x 2 =(−ikcosf) 2 i ; (D:39a) @ 2 r @x 2 =(−ikcosf) 2 r ; (D:39b) @ 2 r @x 2 =(−ihcose) 2 r ; (D:39c) @ 2 i @y 2 =(iksinf) 2 i ; (D:39d) @ 2 r @y 2 =(−iksinf) 2 r ; (D:39e) 126 y r i r e f f x @ 2 r @y 2 =(−ihsine) 2 r ; (D:39f) @ 2 i @x@y =(−ikcosf)(iksinf) i ; (D:39g) @ 2 r @x@y =(−ikcosf)(−iksinf) r ; (D:39h) @ 2 r @x@y =(−ihcose)(−ihsine) r : (D:39i) Now apply the boundary conditions, yy j y=0 =0 and yx j y=0 =0. The first boundary condition yields, yy j y=0 =(−h 2 cos 2 e)Fe −ihxcos e +(+2)(−h 2 sin 2 e)Fe −ihxcos e −2(−k 2 sinf cosfD−k 2 sinf cosfE)e −ikxcos f ; (D:40) For the second boundary condition, yy j y=0 =0, to be true over the entire range ofx, the following requirement, −ihxcose =−ikxcosf; (D:41) is necessary. The above relationship defines e based on f, the incident angle. The first boundary condition also yields the following expression for the unknown coefficients, −h 2 (+2sin 2 e)F−[2k 2 sinf cosf](D−E)=0 : (D:42) The second boundary condition yields the relationship, yx j y=0 =2(−h 2 sinecose)Fe −ihxcos e +(−k 2 sin 2 f)[D +E]e −ikxcos f −(−k 2 cos 2 f)[D +E]e −ikxcos f ; (D:43) With the relationship,−ihxcose =−ikxcosf, the second simultaneous equation for the unknown coefficients is −2h 2 sinecoseF+k 2 (cos 2 f−sin 2 f)[D +E]=0 : (D:44) 127 The simultaneous equations for unknownsE=D andF=D can be written in matrix form as 2k 2 sinf cosf −h 2 (+2sin 2 e) k 2 (sin 2 f−cos 2 f)2h 2 sinecose E=D F=D = 2k 2 sinf cosf k 2 (cos 2 f−sin 2 f) : (D:45) Now define the important material constant, γ 2 = h 2 k 2 = ! 2 = 2 ! 2 = 2 = 2 2 = 2 ; (D:46) which can be related to the Lame constant as γ 2 = 2 2 = = (+2)= = +2 : (D:47) With the new definitions, the matrix equation simplifies to 2sinf cosf −1+2γ 2 cos 2 e sin 2 f−cos 2 f 2γ 2 sinecose E=D F=D = 2sinf cosf (cos 2 f−sin 2 f) : (D:48) The determinant of the matrix on the left side of the equation is =(1−2γ 2 cos 2 e)(1−2cos 2 f)+4γ 2 sinf cosf sinecose: (D:49) More simplification can be accomplished by recognizing thathcose =kcosf, and using the relationship, cosf = h k cose =γcose; (D:50) the determinant can be written in a better form: =(1−2γ 2 cos 2 e) 2 +4γ 2 sinf cosf sinecose: (D:51) Using Cramer’s Rule for the matrix equation, the solution for the unknown coefficients is, E D = 1 2sinf cosf −1+2γ 2 cos 2 e (cos 2 f−sin 2 f)2γ 2 sinecose = −(1−2γ 2 cos 2 e) 2 +4γ 2 sinf cosf sinecose (1−2γ 2 cos 2 e) 2 +4γ 2 sinf cosf sinecose ; (D:52) 128 and F D = 1 2sinf cosf 2sinf cosf sin 2 f−cos 2 f (cos 2 f−sin 2 f) = −4sinf cosf(1−2γ 2 cos 2 e) (1−2γ 2 cos 2 e) 2 +4γ 2 sinf cosf sinecose : (D:53) The value ofD can be determined from the amplitude of the incident SV-wave,jSj,of a shear displacement pulse. The horizontal and vertical displacements of the incident pulse are u i = @ i @y =(iksinf) i ; (D:54a) and v i =− @ i @x =(ikcosf) i ; (D:54b) respectively. Now compute the magnitude,jSj, from the displacements, jSj = p (u i ) 2 +(v i ) 2 = q i 2 k 2 sin 2 fD 2 +i 2 k 2 cos 2 fD 2 =ikD ; (D:55) or D = jSj ik : (D:56) D.3.1 – Wavefields for Incident SV-waves With the unknown coefficients determined, the horizontal displacement for allx and fory 0,is u= @ @x + @ @y =−(ihcose) r +iksinf i − r : (D:57) Using the relationship kγ = ! = ! =h; (D:58) 129 and ksinf =ktanf cosf = tanfkγcose =htanf cose; (D:59) the horizontal displacement within the semi-infinite medium can be written as u(x;y)=ik (−γcose) r (x;y)+sinf i (x;y)− r (x;y) ; (D:60) and the horizontal displacement at the free surface is simply, ju(x;0)j =S 1− E D sinf− F D (γcose) : (D:61) The vertical displacement for allx andy 0,is v= @ @y − @ @x =−(ihsine) r +(ikcosf) i + r : (D:62) Using the relationship kγ =h; (D:63) then v(x;y)=ik (−γsine) r (x;y)+cosf i (x;y)+ r (x;y) : (D:64) The vertical displacement at the free surface of the semi-infinite medium is jv(x;0)j =S 1+ E D cosf− F D γsine : (D:65) Using=1=3 orγ 2 =1=4, the normalized amplitudes for a unit, incident SV-wave are eu ff v ff 60 3.46 0.00 70 1.93 0.63 75 1.94 0.50 80 1.97 0.34 90 2.00 0.00 130 D.4 – Incident SH-Wave Solution Unlike the plane strain wave problems, the SH, or anti-plane wave solution does not require potential functions, its solution satisfies the scalar wave equation and its harmonic wave solution satisfies the scalar Helmholtz equation, r 2 u z +k 2 u z =0 ; (D:66) in whichk =!= is the wave number for shear waves. Figure D.3 – Incident and reflected SH-wave. There is no mode conversion in this case, therefore, the compressional wave potential is absent. The incident SH-wave problem can be written easily by defining the incident wave,u i z , and the reflected wave,u r z ,as u i z =Se −ik(xcos f−ysin f) ; (D:66a) and 131 i r f f x y u r z =Se −ik(xcos f+ysin f) ; (D:66b) respectively. S is the amplitude of the incident wave pulse. It has been shown that a reflected wave with the same amplitude,S, would satisfy the shear stress boundary condition at the free surface as yz = @u z @y y=0 =0 ; (D:67) in whichu z =u i z +u r z is the total displacement in the semi-infinite medium. D.4.1 – Wavefields for Incident SH-Waves The anti-plane displacement,u z , in the semi-infinite medium can be written as u z (x;y)=Se −ikxcos f cos(kysinf) : (D:68) This is the summation of the incident and reflected waves as shown in Eq. (D.66). This shear wave only solution is not dependent on Poisson’s Ratio,. The displacement amplitude at the free surface is a constant, ju z (x;0)j=2S: (D:69) D.5 – Rayleigh Wave Solution Unlike the body waves described in the previous sections, a semi-infinite medium is amenable to a surface wave and its larger amplitudes are confined near the free surface. If a harmonic wave function has the form =Ae −i ! (xcos e−ysin e) ; (D:70) it has an apparent velocityc on the free surface (y=0). Consider the factor: ! xcose = !x =cose = !x c ; (D:71) 132 the apparent velocity is defined asc ==cose for this particular case. Ife! 90 ,c!1 because all particles on the surface appears to move in the same direction, therefore, it appears on the surface as if the wave is moving with an infinite speed. Consider the case whenc<< by defining the compressional potential and shear potential as =Ae −i ! c (x p c 2 2 −1 y) ; (D:72a) =Be −i ! c (x q c 2 2 −1 y) ; (D:72b) respectively. For the condition,c<<: q c 2 2 −1 is complex. The second partial derivatives of the potential functions can be listed as @ 2 @x 2 =(−i ! c ) 2 ; (D:73a) @ 2 @x 2 =(−i ! c ) 2 ; (D:73b) @ 2 @y 2 = i ! c r c 2 2 −1 ! 2 ; (D:73c) @ 2 @y 2 = i ! c s c 2 2 −1 !2 ; (D:73d) @ 2 @x@y = −i ! c i ! c r c 2 2 −1 ! ; (D:73e) @ 2 @x@y = −i ! c i ! c r c 2 2 −1 ! ; (D:73f) The first boundary condition yields, yy j y=0 =(+2)r 2 −2 @ 2 @x 2 −2 @ 2 @y@x =(+2) − ! 2 c 2 − ! 2 c 2 c 2 2 −1 Ae −i!x=c +2 ! 2 c 2 Ae −i!x=c −2 " ! 2 c 2 s c 2 2 −1 # e −i!x=c =0 133 Remove the harmonic time factor and divide the equation by (!=c) 2 and (+2), the equation simplies to c 2 2 −2γ 2 A 2γ 2 s c 2 2 −1 ! B=0 (D:75) Similarly, the second boundary condition, yx j y=0 =0=2 @ 2 @y@x + @ 2 @y 2 − @ 2 @x 2 (D:76) can be simplified to 2 r c 2 2 −1 ! A+ c 2 2 −2 B=0 (D:77) Using the relationship c 2 2 1 γ 2 = c 2 2 2 2 = c 2 2 The first equation can be written as 2− c 2 2 A2 s c 2 2 −1 ! B=0 (D:78) The simultaneous equations can be written as a matrix equation as 2 6 6 6 4 2− c 2 2 2 s c 2 2 −1 2 r c 2 2 −1 2− c 2 2 3 7 7 7 5 8 < : A B 9 = ; = 8 < : 0 0 9 = ; The only nontrivial solution which exists is when the determinant of the above equation is zero, i.e., 2− c 2 2 2 +4 r c 2 2 −1 ! s c 2 2 −1 ! =0 A simpler transcendental equation can be obtained by defining a parameter,q = c= and using the fact thatc 2 = 2 =γ 2 q 2 , that q 6 −8q 4 +(24−16γ 2 )q 2 +16(γ 2 −1)=0 134 Appendix E Three-Dimensional Isoparametric Element Using the concepts of an isoparametric element, the shape function defined in natural coordinates,, and, are N i (;)= 1 8 (1+ i )(1+ i )(1+ i ) ; (E:1) in which i ; i ; i are the natural coordinates at the eight nodes defined in the order: ( 1 ; 1 ; 1 )=(−1;−1;−1) ; (E:2a) ( 2 ; 2 ; 2 ) = (+1;−1;−1) ; (E:2b) ( 3 ; 3 ; 3 ) = (+1;+1;−1) ; (E:2c) ( 4 ; 4 ; 4 )=(−1;+1;−1) : (E:2d) ( 5 ; 5 ; 5 )=(−1;−1;+1) ; (E:2a) ( 6 ; 6 ; 6 ) = (+1;−1;+1) ; (E:2b) ( 7 ; 7 ; 7 ) = (+1;+1;+1) ; (E:2c) ( 8 ; 8 ; 8 )=(−1;+1;+1) : (E:2d) The physical coordinates,x andy, within the element can be determined as x = 8 X i=1 N i (;;)x i ; (E:3) y = 8 X i=1 N i (;;)y i ; (E:4) z = 8 X i=1 N i (;;)z i ; (E:5) 135 Figure E.1 – Node Numbering Scheme for a Three-Dimensional Solid Element. in which (x i ;y i ;z i ) are the physical coordiantes of thei-th node, the limits of the natural coordinates are−1 1,−1 1 and−1 1. To calculate strains and stresses, the derivatives of the shape functions are required. Use the chain rule to relate the derivatives in the natural coordinates and the physical coordinates as 2 6 6 6 6 6 6 4 @N i @ @N i @ @N i @ 3 7 7 7 7 7 7 5 = h J i 2 6 6 6 6 6 6 4 @N i @x @N i @y @N i @z 3 7 7 7 7 7 7 5 ; (E:6) 136 x y z 1 2 3 4 8 7 5 6 in which [J] is the three-dimensional Jacobian matrix defined as h J i = 2 6 6 6 6 6 6 4 @x @ @y @ @z @ @x @ @y @ @z @ @x @ @y @ @z @ 3 7 7 7 7 7 7 5 : (E:7) Inversely, 2 6 6 6 6 6 6 4 @N i @x @N i @y @N i @z 3 7 7 7 7 7 7 5 = h J i −1 2 6 6 6 6 6 6 4 @N i @ @N i @ @N i @ 3 7 7 7 7 7 7 5 : (E:8) The partial derivatives in the above expressions can be computed as @N i @ = 1 8 i (1+ i )(1+ i ) ; (E:9a) @N i @ = 1 8 i (1+ i )(1+ i ) ; (E:9b) @N i @ = 1 8 i (1+ i )(1+ i ) ; (E:9c) and @x @ = 1 8 8 X i=1 x i i (1+ i )(1+ i ) ; (E:10a) @x @ = 1 8 8 X i=1 x i i (1+ i )(1+ i ) ; (E:10b) @x @ = 1 8 8 X i=1 x i i (1+ i )(1+ i ) ; (E:10c) @y @ = 1 8 8 X i=1 y i i (1+ i )(1+ i ) ; (E:10d) @y @ = 1 8 8 X i=1 y i i (1+ i )(1+ i ) : (E:10e) 137 @y @ = 1 8 8 X i=1 y i i (1+ i )(1+ i ) : (E:10f) @z @ = 1 8 8 X i=1 z i i (1+ i )(1+ i ) ; (E:10g) @z @ = 1 8 8 X i=1 z i i (1+ i )(1+ i ) : (E:10h) @z @ = 1 8 8 X i=1 z i i (1+ i )(1+ i ) : (E:10i) To calculate the three-dimensional strains, define the strain vector, ~ " and the derivative matrix, [L],as ~ "= 2 6 6 6 6 6 4 " xx " yy " zz γ xy γ yz γ zx 3 7 7 7 7 7 5 = h L i 2 4 u x (x;y;z) u y (x;y;z) u z (x;y;z) 3 5 = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 @ @x 00 0 @ @y 0 00 @ @z @ @y @ @x 0 0 @ @z @ @y @ @z 0 @ @x 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 2 4 u x (x;y;z) u y (x;y;z) u z (x;y;z) 3 5 : (E:11) Let the displacements, u x (x;y;z), u y (x;y;z) and u z (x;y;z), be represented as linear combinations of the shape function,N i (x;y), and nodal displacements,u xi ,u yi andu zi , respectively as u x (x;y;z)= 8 X i=1 N i (x;y;z)u xi ; (E:12a) u y (x;y;z)= 8 X i=1 N i (x;y;z)u yi ; (E:12b) u z (x;y;z)= 8 X i=1 N i (x;y;z)u zi ; (E:12c) 138 or in matrix form, 2 4 u x (x;y;z) u y (x;y;z) u z (x;y;z) 3 5 =[N]~u; (E:13) in which [N]= h [N 1 ][N 2 ][N 3 ][N 4 ][N 5 ][N 6 ][N 7 ][N 8 ] i ; (E:14) [N i ]= 2 4 N i (x;y;z)0 0 0 N i (x;y;z)0 00 N i (x;y;z) 3 5 ; (E:15) ~ u = 2 6 6 6 6 6 6 6 6 6 4 ~ u 1 ~ u 2 ~ u 3 ~ u 4 ~ u 5 ~ u 6 ~ u 7 ~ u 8 3 7 7 7 7 7 7 7 7 7 5 : (E:16) and ~ u i = 2 4 u xi u yi u zi 3 5 : (E:17) The constitutive properties applicable for stresses in relation to strains can be represented by defining the material constant matrix, [D]= E(1−) (1+)(1−2) 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 1 (1−) (1−) 000 (1−) 1 (1−) 000 (1−) (1−) 1 000 000 (1−2) 2(1−) 00 000 0 (1−2) 2(1−) 0 000 0 0 (1−2) 2(1−) 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ; (E:18) 139 in whichE is the Young’s modulus and is the Poisson’s Ratio. The shear stress vector ~ = 2 6 6 6 6 6 4 xx yy zz xy yz zx 3 7 7 7 7 7 5 =[D][L][N]~ u=[D][B]~u; (E:19) in which [B]=[L][N] is 624 matrix. The 2424 element stiffness matrix, [K], can now be evaluated as [K]= Z V [B] T [D][B]dV ; (E:20) using numerical integration in the natural coordinate system with the help of the Jacobian transformation to convert, dxdydz =Jddd : (E:21) The 2424 consistent mass matrix can also be evaluated with the shape functions as [M]= Z V [N] T [N]dV ; (E:22) in which is the mass density. 140
Abstract (if available)
Abstract
Using the Reciprocal Theorem and the principle of superposition, this dissertation examines three procedures for obtaining the impedance functions and driving force vectors necessary for a soil-structure interaction analysis. The first procedure employs a volume integral of body forces and the displacement Green's Functions to obtain the frequency dependent force-displacement relationship (impedance functions) for a rigid embedded foundation and the generalized force vector required to hold the rigid foundation fixed while subjected to incidentwaves (driving force). Although a formulation based on a volume integral requires more computation than the others based on surface integrals, this procedure proves to be numerically stable and is an excellent tool for linear soil-structure interaction analyses. The second procedure is based on integrating the surface displacements with the traction Green's Functions and surface tractions with the displacement Green's Function, over the surface of a rigid foundation. This well-knownformulation has numerical instability problems at the resonant frequencies of the medium interior to the foundation surface. This dissertation explains the reason behind the instability using an analytical solution of the Hilbert-Schmidt method and recommends an algorithm which uses the l'Hospitale Rule to obtain the solution at those critical frequencies. The development herein for this method is not meant to be a practical solution, but an academic exercise to explain a long-standing, puzzling problem. The third and final procedure is to design a radiation boundary for a three-dimensional finite element grid using the surface displacements and surface tractions at the discrete artificial boundary, along with the Green's Functions, to calculate outgoing wave motion at nodes immediately outside of the artificial boundary, thereby eliminating unwanted reflection of the outgoing waves. A procedure is recommended to obtain required data from a commerci
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Soil structure interaction in poroelastic soils
PDF
The role of rigid foundation assumption in two-dimensional soil-structure interaction (SSI)
PDF
Diffraction of SH-waves by surface or sub-surface topographies with application to soil-structure interaction on shallow foundations
PDF
Train induced vibrations in poroelastic media
PDF
In-plane soil structure interaction excited by plane P/SV waves
PDF
Wave method for structural system identification and health monitoring of buildings based on layered shear beam model
PDF
The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems
PDF
Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations
PDF
Modeling and dynamic analysis of coupled structure-moving subsystem problem
PDF
Modeling and simulation of circulating tumor cells in flow. Part I: Low-dimensional deformation models for circulating tumor cells in flow; Part II: Procoagulant circulating tumor cells in flow
PDF
Novel multi-stage and CTLS-based model updating methods and real-time neural network-based semiactive model predictive control algorithms
PDF
Effects of the thermo-mechanical history on the linear shear viscoelastic properties of uncrosslinked elastomers
Asset Metadata
Creator
McKay, Kirsten
(author)
Core Title
Three applications of the reciprocal theorem in soil-structure interaction
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Civil Engineering (Earthquake Engineering)
Publication Date
03/06/2009
Defense Date
01/16/2009
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
embedded foundations,foundation impedance,Green's functions,OAI-PMH Harvest,reciprocal theorem,soil-structure interaction
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Wong, Hung Leung (
committee chair
), Lee, Vincent W. (
committee member
), Moore, James Elliott, II (
committee member
)
Creator Email
kirsten.mckay@wgint.com,kmckay@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m2004
Unique identifier
UC173640
Identifier
etd-McKay-2628 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-214820 (legacy record id),usctheses-m2004 (legacy record id)
Legacy Identifier
etd-McKay-2628.pdf
Dmrecord
214820
Document Type
Dissertation
Rights
McKay, Kirsten
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
embedded foundations
foundation impedance
Green's functions
reciprocal theorem
soil-structure interaction