Seismic Soil Structure Interaction (SSI) Analysis

Authors: Andre Hergersheimer Ruzsicska  @Rushdi Qader and Unbekannter Benutzer (ge92sic) 

Supervised by Sebastian Schopper, M.Sc. and Dr.-Ing. Francesca Taddei

Abstract

SSI studies the changes of stresses, strains and displacements between ground and structure. These changes influence each other simultaneously, soil to structure and structure to soil. This increases the complexity of the system and makes the calculation harder. The purpose of this project is to facilitate a Simulation of all these interactions as near from reality as possible. This project focus on improving two separated functions/codes with distinct features, so in the end they will work as one and simulate the seismic load under the plate.  

Theory

Introduction 

SSI studies the interaction between soil and foundation, applying this theory makes the analysis of dynamically coupled loads between the two substructures possible on the basis of taking into account their dynamic properties. The importance of this theory is that it deals with external and internal loads. The external loads on the investigated plate and ground in this report are:1D-cosine, 1D-seismic-x, 1D-seismic-xy, 2D-sesmic-xy, 2D-Ricker-x, 2D-Ricker-z. Internal loads take place firstly by the internal interaction as a result of transmitted internal loads to the free field substructure i.e. ground. The transmitted stresses that scattered away from the plate propagates as waves and are denoted as radial damping. Secondly, as the plate is stiffer compared to the soil, there will be additional interactions due to the seismic loads because the plate is not able to adapt to the deformations of the ground. The plate reflects the seismic waves and produces a wave field due to the scattering of the waves and the motion will be modified at the interaction nodes between soil and plate.

Numerical model and methodology

Domain decomposition

For finding the dynamic response of the plate that rests on the soil with the determined boundary conditions, the analysis based on domain decomposition is the most suitable method. Firstly, the structure is divided into two substructures to decouple the stiffness matrix at the interaction nodes between soil and plate, then combining the substructures by enforcing the consistence of displacements and equilibrium of resulting stresses at the interface between soil and plate. For the separated free field, as the vertical distance (z) decreases the displacement amplitude decays because of geometric damping and material damping. Here the stages for modeling Soil-Structure-Interaction numerically are illustrated:

1. Developing time-domain analysis using a direct method.
2. Non-linearity of the interface conditions between the two substructures and nonlinear behavior of the plate has to be considered for the analysis.
3. As the soil has high material damping the non-linearity of the viscosity should be considered.
4. For examining the behavior of the structure due to seismic loads a 3D-validation examination has to be performed.
5. Writing all the previously mentioned structural and material parameters into a software which is capable of calculating the displacements and stresses in all the three dimensions(x,y,z).

Soil substructure 

Layered half space (ITM)

The coupled differential equation of the vector field of displacements can be introduced by substituting the relations between stresses and strains and displacements:

$//<![CDATA[ \begin{array}{l}μu^i |_j^i+(λ+μ) u^j |_j^i-ρü^i=0\end{array} //]]>$ 

with i,j = (x,y,z), λ: first lame parameter, μ: second lame parameter, ρ: density and u: displacement. 

The coupled partial differential lame equation is essential for describing the half space which consists of one layer. The adhesion behavior of the layer is characterized by compressional and shear waves P and S respectively: $//<![CDATA[ \begin{array}{l}c_(_p,_l_)=\sqrt{(λ_l+2μ_l)⁄ρ_l }\end{array} //]]>$ and $//<![CDATA[ \begin{array}{l}c_(_s,_l_)=\sqrt{μ_l⁄ρ_l }\end{array} //]]>$(with Ɩ = 1). The lame equation is decoupled by Helmholtz approach and a threefold Fourier transformation is performed to transform the equation into the wavenumber, frequency domain (k_x, k_y, ω):

Fourier transformation

after finding the displacement in the wavenumber frequency domain, a threefold inverse Fourier transformation is performed to obtain the displacements in space domain (x, y, z, t):

and thus, illustrated in Fig. (1):

Fig. (1). Decoupling of partial differential equations  

The vibration of the soil depends on two key aspects: 

• The compressional and shear wave velocities which give rise to P and S waves. The sources of those wave velocities are: Strong motions, Ricker waves, impulses.
• Type of boundary, for a homogeneous half space if there are no boundary conditions no wave will be reflected into the surface as there is no reflecting material like bedrock and the waves scatter away

Soil-Plate coupling (ITM-FEM)

The plate substructure which is coupled to the soil at the contact area on the surface of the soil is modelled with a 3D finite element method and discretized in accordance to the nodal locations of the soil that where it coincides. The nodes interact with the soil at 81 interaction nodes and each node has three degrees of freedom.

The used plate is very thin in depth so it can also be treated as a volume element.

Finally, by substituting the required elements into the equation of motion, we obtain a discretized equation of motion:

By enforcing the equilibrium and compatibility conditions, the FEM model will coincide with the ITM model along the discretized points at the interaction nodes and thus we can combine the two substructures, this is achieved by coupling the dynamic stiffness matrix of the soil and plate Kf and Ks respectively, which provides us the total equation for the ITM-FEM model, that is used for calculating the nodal forces and stresses at the soil plate interaction nodes. We can use the contact stresses as the applied load by transforming it into the wavenumber frequency domain at the surface (z = 0), and thus obtaining the coupled displacements and loads.

Fig. (3) Soil-plate system geometry, Nonlinear rocking stiffness of foundations

G. Gazetas(2013)

Given codes 

Now that we have the general theory of our project and its goal, which is to compute the dynamic response of the soil-plate model under the influence of harmonic loads . We use two distinct codes, which already existed in the Structural Mechanics Chair's MATLAB based FEM code, in order to obtain the seismic load between soil and plate

The first Code, called Main, is a free field problem, there is no plate or other kind of structure over the soil. It has a soil modelled using Integral Transform Method (ITM). It will calculate the loads inside the soil under the influence of an external load We(t) which is placed at the bottom of the soil i.e. at z = 300m, it simulates the reaction of all nodes, and then returns their displacement.

Fig. (4) Non-exhaustive schema of the Main

The second Code, called SsiExampleltmfem, is an example based on SoilSuperElement. It’s a another code in Structural Mechanics chair’s Matlab based FEM that will compute the dynamic behaviour of a homogeneous halfspace. SsiExampleltmfem has a soil model, it uses Integral Transform Method (ITM) to calculate the nodes. This soil is coupled with a plate that is modelled with FEM. The stiffness matrix from soil and plate are coupled at the interaction nodes. The original code could calculate the given harmonic loads over the ground or plate, and only the coupling-nodes between soil and plate can receive internal loads. 

Fig. (5) Non-exhaustive schema of the SsiExampleltmfem

Each Code was used separately for other projects for a long time. In order to make them work together, some changes are necessary in one or both codes:

• Discretization
• Size of soil
• Material properties

Then, after adapting our models, we need the displacements at the interaction nodes and the stiffness matrix matched with the frequencies in the Main to be merged to the SsiExampleltmfem code. Finally, we can calculate the seismic load as explained in the theory, by multiplying the stiffness matrix with the displacements.

Adapt the model

The first task is to have the same soil in both codes. For this project, the Main is a code that was designed for seismic analysis, it is then the most adapted to changes in soil size and material properties, that is why it will serve as a reference soil. The SsiExampleltmfem is a code designed for an article and its method is based on the literature, it is less conducive to changes and we will now explain how we proceeded.

First, we had to change the dimensions as well as the properties of the material. The main changes were implemented in two functions that we developed: Material.m and Soil.m. The first one allows us to get the same material properties for both codes, the second one allows us to have both codes in one in a compact form, then we merged the important parameters (i.e. material and structural) in both codes. Below are the functions in more detail. 

Material.m 

It’s a function with all input data needed. It will return all necessary and useful soil parameters, but without the material damping. Material damping is important to take in account, it is further coded in both code but differently, we will explain further how we obtain the same final result.

We have made this function because it has the advantage that if we want to change the material properties, we can change the parameters of our function, and it will change for both codes. It was benefitial for saving time and and preserving accuracy in data calculation and stacking . The following flowcharts represent the material function:

Soil.m

The soil function is implemented in the SsiExampleltmfem. First, it will run the Main code, then it will return all important data which is necessary to have the same size of soil and to compute the seismic load. The following flowcharts represent the soil function:

We can now implement the same type of soil, but we must be careful when adapting SsiExampleltmfem to the dimensions of the Main. Indeed, as explained at the beginning, we want to use our project for seismic analysis, the dimensions are then large and we have to keep a certain order of magnitude in order to be able to continue to compare the results with those of the literature. Initially, the floor was very small, 16[m]. To increase its size, we had to keep the number of nodes constant in the x and y directions, so we changed the discretization which is equal to the width of the Main's soil divided by the number of nodes. After changing the dimensions of the soil, we need to change the dimensions of the plate and its load accordingly. To do this, we need to calculate the scale factor between the width of the soil and the plate. Initially we had a 16-meter-wide floor for a 1-meter plate. It follows that we have a scale factor equal to 16 and we have implemented it in the code. In the same way, we found that the load of the plate must be equal to half its size.

Material damping 

Finally, we must now take into account the material damping in our codes. It is implemented in the Main by adding a complex part to the shear modulus G [N/m²] because the soil has a hysteretical damping behaviour:

$//<![CDATA[ \begin{array}{l}G^ ̃=G(1+sign(ω)2iξ)\end{array} //]]>$

with ω the frequency [Hz] and ξ [-] the hysteretic material damping ratio (ζ in SSIExampleltmFem). We can notice that in the code, only negative frequencies are used, so the formula is coded with a "-" sign.

For the SsiExampleltmfem, it is taken into account by adding the complex part to the young's modulus E [N/m²]. The characteristics of an isotropic material must be used: 

$//<![CDATA[ \begin{array}{l}E=2G(1+ν)\end{array} //]]>$

$//<![CDATA[ \begin{array}{l}E^ ̃=E(1+sign(ω)2iζ)\end{array} //]]>$

with ν the Poisson’s ratio [-]

The damping occurs in any system; it is related to the energy that is dissipated in the system and it is constant throughout the material because both of the substructures are homogeneous. Moreover, in a model using ITM, the damping is of great use to «avoid spurious numerical effects within the ITM solution by preventing the superposition of the periodically repeated signal in the wavenumber domain» (1, page 8)

Reshape the displacement of the free-field

As it was mentioned in last chapter in Main-function we obtain the displacement of each node from the soil. The displacement u1_f depends on first dimension of the frequency of load. Then we have the dimensions with coordinates X, Y and Z to inform where on the soil our node is located. The fifth dimension is the degree of freedom, which means a variation between one of the three degrees in continuum. The next figure illustrates the given displacement: 

The following videos show the soil displacement in Main code for 26 different frequencies (left) and 5 different frequencies (right): 

                                                             10Hz_26f.mp4 freq_5_for10Hz.mp4

However, the stiffness matrix in SsiExampleltmfem does not correspond with this type of distribution. The first dimension of the K-matrix will give a displacement in the direction for each node. The second dimension gives the direction of load applied to each node and the third one is the value of the frequency as shown in the following figure:

u1_f_surface

As it was explained in the Theory chapter, in order to be able to calculate the seismic load, we need to multiply the stifness matrix with the displacement. This means a change of dimensions for the displacement u1_f.  To transform the 5D matrix into a vector, we first knew that only the nodes on the surface were important for the stiffness matrix, so we can take only values, where Z=0, reducing one dimension and creating the 4D displacement matrix u1_f_surface. In Matlab it is possible to solve this operation by storing only the corresponding layer for Z or fourth dimension of u1_f. 

Sorted.m

After that, we need to sort the nodes under the plate from the remaining soil nodes which do not interact with the plate model. This will be the first part of the Sorted.m code, it takes the coordinates (X and Y) of soil and its displacement u1_f_surface (U) from Main function,  and thus the coordinates of the plate (plate) from SsiExampleltmfem. With this the program will create variables to respectively store the length, first and last node of plate and how many x-points there are in Main function. As we will see, the combination of x- and y-values for the plate works in such a way that the first node gives the two first values and the last node gives the last two values. Knowing the length (l) we can store these values in variables (beginn_x, beginn_y, ende_x and ende_y). In our project the plate is square and combines nine x-values with nine y-values, which gives the following coordinates to the nodes (index 2 and 3):

Then, with a loop the corresponding indexes for start and end of our plate in the Main function (index_beginn_x and index_end_x) will be searched for in the X-coordinates (lx). For reasons of time, the Sorted.m code was created simutaneously with the adaptation of soil size in Main function. The result was that the nodes in Main function were not always the same in SsiExampleltmfem function, so again with a loop we will check how many there are on plate. As represented in the last table, all x-values on the plate are given in the first series to generate the nodes on the plate, so we can in this loop increase the index for checking the number by 1 per cycle. The moment we finish the first series, the values will repeat, in order to save time and memory, we create a break when the first number repeats. The last loop will ensure  that the corresponding indexes from Main function are stored for the x-values on the plate. Although the plate in our example comes with a squared shape, we have developed a general code for rectangular plates. This means that the code should be able to calculate how many y-values there are on the plate separate from x-values. Since the node table for our plate shows that the y-values are generated in an outer part of a loop, we need to multiply in the second loop the current index of the number of x-values and increase by 1 per cycle, again when the code checks for a repetition of the number it will break the loop. The next flowchart represents that whole part of Sorted.m code:

Then, the second part of the Sorted.m code will generate a new displacement vector with three index types. The first should give the combination of coordinates with degrees of freedom. The second index indicates that the new u1 is a column vector. The last index is the frequency number. In our example:

• Plate's node: 81
• Degrees of freedom: 3
• First dimension of new u1: 243 (81*3)
• Second dimension of new u1: 1
• Third dimension of new u1: 5

In this way, the displacement u1_f_surf_sorted matches with the stiffness matrix. In the next step of the Sorted.m code the length of stored indexes (lp_x_end and lp_y_end) will be created. And generates the number of degree of freedom per node (nDof) and frequency (omega). These values will define our new vector with already known sizes and zeros values, which will save memory in MATLAB and will not allow MATLAB to risk changing the number type when the value given for displacement is zero. It will then start a series of loops to combine the coordinates with the degrees of freedom by frequency and store the corresponding displacement value of u1_f_surface in u1_f_surf_sorted, as illustrated in the next flowchart:

Generating the seismic load

Onlypos.m

Since the code for stiffness matrix K_soil_found_f and K_soil_found_pos_f were given by the Structural Mechanics chair, it will not be delved into in this project how the calculation works. But it is important to note some differences between the Main and SsiExampleltmfem code with respect to frequency input. The Main code can generate a large number of frequencies (negative and positive) always with frequency zero as center value and the SsiExampleltmfem could only receive by direct encoding positive non-zero frequencies. To use the same frequencies in both codes, we store the frequencies from Main code in f_main through Soil.m and with a new function Onlypos.m we take only the positive frequencies and the frequency zero as shown below:

To be able to calculate the stiffness matrix in SsiExampleltmfem for frequency = 0 , it was necessary to implement in the same way as it was done for Main code a if-case for the frequency loop, when frequency = 0 the code will change the value to 0,001 in the loop. In our example with 5 frequencies (-10, -5, 0, 5, 10) the frequency steps for the stiffness are:

-K_soil_found_f

• 1° frequency = 0
• 2° frequency = -5
• 3° frequency = -10

-K_soil_found_pos_f

• 1° frequency = 0
• 2° frequency = 5
• 3° frequency = 10

The frequency steps for the displacement vector are:

-u1_f_surf_sorted

• 1° frequency = -10
• 2° frequency = -5
• 3° frequency = 0
• 4° frequency = 5
• 5° frequency = 10

Seismicload.m

The last step to combine the stiffness matrix and the new displacement vector is to sort the values with the correct frequencies. To do this we create the function Seismicload.m with K_soil_found_f, K_soil_found_pos_f and u1_f_surf_sorted, first the code will check how many frequencies are in K Matrix (n_freq). Similarly to the Sorted.m code, we will generate a vector S with already known sizes (n_dofs and n_freq_vector) and values of zeros. The first loop will store the values out of multiplication of stiffness with displacement for negative frequencies, for this the loop goes with a step of -1. Keeping the index of the last loop (j) we now store in S the multiplication for positive frequencies. It is not necessary to repeat the calculation with frequency = 0 in second loop. The next flowchart illustrates this process:

Final results

To facilitate our work in visualization means and running time, we use in this project only five frequencies (-10, -5, 0, 5 and 10), as well as focusing on the simplest load case, a constant cosines excitation in the  x-direction for all bottom nodes. With the new functions presented in this report the soil parameters for both codes are:

• B_y = 960 [ m ]
• B_x = 960 [ m ]
• H = 300 [ m ]
• Layer = 1

The plate dimension are:

• B_y = 60 [ m ]
• B_x = 60 [ m ]
• H = 0,3 [ m ]

For the chosen load case the maximum displacement and seismic load were found in the corner of the plate:

• U = 2,5*e-5 [ m ]
• S = 2,3*e4 [ N ]

To give a general idea of the absolute seismic load in the plate by frequency, the next graphic plots the result for the middle node:

Conclusion

The given MATLAB codes have been successfully merged into one code file. The SsiExampleltmfem code can now call the soil parameters, size, dimension, frequency as well as displacement of soil from the Main code and implement it as its own variables. With this it is possible to calculate the seismic load using the displacement generated in the free field and the stiffness matrix at the coupling nodes. Despite having a limited number of frequencies, the results obtained for displacement and seismic load are plausible with those known in the literature. In the final results it is possible to observe the symmetry property between the negative and positive frequencies, with only a peak at one of frequencies, which could indicate an approximation with the natural frequency of the plate. The new functions are designed to adapt and accept different rectangular plate shapes and to expand or reduce the soil dimension if necessary taking into account the permitted level of deformations.

Further Research 

For future research, the program can vary the quantity of layers and their properties, as well as working with a large frequency spectrum for different types of loads is possible. In this way, the SsiExampleltmfem is closer to representing the deformations as a result of a seismic load on a structure.

Now that both codes are merged and the flow of variables is easier, it is possible to further improve the new created functions, so that they could save more runtime and memory for more extensive calculations. 

Reference 

[1] J. Freisinger, “Dynamic response of three-dimensional rigid and flexible foundations on layered soils with local inhomogeneities”, February 2022. [Online]. Available: https://doi.org/10.1016/j.soildyn.2021.107007 [Accessed 18.02.23]

[2] L. Chen, “Numerical models for the analysis of soil, structure and their interaction”, January 2016. [Online]. [Accessed 17.02.23]

[3] N. Tsai, G. Housner, “Calculation of surface motions of a layered half-space”, October 1970. [Online]. Available: https://authors.library.caltech.edu/49716/1/1625.full.pdf [Accessed 17.02.23]

[4] Braja M. Das, “Advanced Soil Mechanics”, Fifth Edition