F8. Characteristic scales

Firstly, to frame the discussion about soil dynamics problems, it is necessary to distinguish the characteristic scales relative to the definition of the physical phenomena studied in these problems. These scales are:

the geological scale,
the site scale,
the geotechnical scale,
the civil engineering structure scale.

The first issue one can encounter is that each problem in soil dynamics involves one or more of these scales. Indeed, it is important to associate in the same model, domains with extremely different dimensions.

Moving from one scale to another corresponds to a specific type of problem, which requires its own solving methods. Taking seismic problems as an example, the transition from geological to site-scale raises the issue of determining local site effects. Moving to the geotechnical scale highlights the problem of characterizing seismic motion at a control point (which may be located at the soil surface or a characteristic depth). Finally, scaling up the structure underlines the problem of soil-structure interactions, i.e. the problem of understanding how the dynamic response of the soil is modified by the presence of the structure and vice versa.

Depending on the problem studied, it is also possible that the scale of the structure interacts directly with the site scale or the geological scale. This is especially true for large structures such as bridges (structure scale ~ site scale) or tunnels (structure scale ~ geological scale).

In this section, the discussion on geotechnical and structural scales (characteristic scales for most civil engineering structures) will be limited to the problems of Dynamic Soil-Structure Interaction (DSSI).

Classification of the methods for DSSI problems

Before proposing a method classification for DSSI problems, it should be noted that the dynamic behaviors of the studied systems are made of two distinct domains:

the structural domain, including the studied structure and possibly an area of soil in the vicinity of this structure. This domain highlights bounded dimensions, and is commonly called the generalized structure domain or superstructure or simply structure,
the soil domain, surrounding this structure, whose dimensions are unbounded.

The two domains are separated by a boundary which is called the interaction horizon. This boundary is fictitious, and its definition is, in a sense, arbitrary. For example, the interaction horizon can be defined as the real interface between the structure and the soil, or it can be placed very far from the structure.

On the other hand, since DSSI problems are like wave propagation problems, the adopted formulation must satisfy a radiation condition. For example, in the case of a load applied on the structure (i.e., the structure is identified as the source of the waves), the radiation condition requires that there can be no waves converging towards the source.

The methods to solve DSSI problems are divided into three main families:

sub-structure methods,
direct methods,
hybrid methods.

Sub-structure methods are the most common for Civil Engineering DSSI studies. The interaction horizon in this case is placed at the physical interface between the soil and the structure. This separation makes it possible to model the "structure" domain using classical finite elements and the "soil" domain by other means, which are better suited for the treatment of unbounded domains, such as the boundary element method (BEM or one of its variants). The radiation condition is rigorously formulated to infinity. The sub-structure method acquires its practical interest if it is combined with the linear hypothesis for the structure and soil domains. In this case, it can be shown, via Kausel's superposition theorem, that the global soil-structure interaction problem can be decomposed into three sub-problems, as shown in the figure below.

Kausel superposition theorem in the substructure method

The kinematic interaction concerns the diffraction of waves in the soil by the presence of the foundation and its substantial loading. The inertial interaction concerns the exchange of energy between the soil and the structure which can be characterized in terms of impedance. These first two problems provide respectively the driving motion and the boundary conditions that can be used in the third sub-problem for the dynamic analysis of the bounded domain. Consequently, the sub-structure method uses finite elements only to model the structure. One can then refer directly to the chapters explaining the use of finite elements in conventional structural dynamics studies. One of the disadvantages of Kausel’s superposition theorem is related to the nature of geomaterials, which have an intrinsic non-linear mechanical behavior. Therefore, to call upon the hypothesis of linear behavior for the "soil" domain, it is necessary to proceed to a calibration of the viscoelastic parameters of the geomaterials according to the anticipated levels of deformations. Defining the effective viscoelastic parameters of the soil layers constitutes the first step to implement the sub-structure method, it allows us to use the superposition theorem in systems where the non-linear behavior remains moderate.

The direct methods are the second large family of approaches for the DSSI and can be applied for systems with strong non-linear behaviors. The interaction horizon in this case is placed very far from the structure, encompassing a large area of the "soil" domain, as shown in the figure below.

The general principle of direct methods

The dynamic equilibrium equation is formulated for the "structure + soil" assembly and involves an external loading vector qf defined along the interaction horizon. The "structure + soil" set is modeled by the FE method and the displacements in the "soil" and "structure" domains are computed simultaneously. The calculations are performed in the time domain. It is then possible to incorporate in the model all the geometrical characteristics of the problem, the material heterogeneities of the ground or the superstructure, and to introduce the constitutive laws necessary for the description of the non-linearities of the system. Despite these important advantages, the direct methods present several issues:

defining the optimal position for the interaction horizon and identifying the external loads applied onto it (especially for seismic problems) are non-trivial tasks, often requiring substantial preparatory work,
in contrast to the sub-structure methods, where the radiation condition is rigorously formulated at infinity, in the direct methods, the radiation condition is roughly formulated at the level of the interaction horizon. This hypothesis allows us to neglect the unbounded domain in the model, but special elements must be implemented along the interaction horizon (which is identified with the boundary of the modeled domain). These elements are supposed to preserve the unbounded character of the modeled domain by ensuring that there will be no spurious wave reflections along the boundary of the modeled domain,
the spatial and temporal discretization (mesh size, time step) depend on the mechanical parameters of the modeled system, and often the mesh size or time step required for a problem can be expensive to compute,
the time integration scheme must satisfy the accuracy and stability requirements of the solution while remaining optimized with respect to the computational cost,
the constitutive laws admitted that describe the non-linear behavior of soils require a calibration to represent the cyclic behavior of soils,
finally, direct methods are still extremely expensive for 3D problems. Therefore, it is necessary to be able to switch from 3D to 2D modeling while preserving the dynamic characteristics (mass/stiffness / modal characteristics) of the studied system.

Finally, hybrid methods are the third large family of methods for DSSI processing. These methods are in-between direct and sub-structure methods. The main idea is to separate the soil into two distinct domains: the first is a field near the structure where all the non-linearities relevant to the DSSI problem are assumed to be developed. The second is the far-field where the soil behavior is not affected by the interaction with the structure. Thus, the far-field can be handled by techniques suitable for linear problems in unbounded domains (typically BEM with a resolution in the frequency domain), while the near field is incorporated in the model of the superstructure and can be handled by a direct method (FE modeling - non-linearity modeling - resolution in the time domain). The general principle of hybrid methods is presented in the figure below.

The general principle of hybrid methods

The issue with hybrid methods is the definition of the boundary between the near field and the far-field (interaction horizon) as well as the formulation of a scheme that restores rigorously and efficiently the coupling between the two fields. Numerically, this operation results in a FEM-BEM coupling scheme, where the BEM component feeds the FEM component with the histories of the interaction forces and impedance functions to be deployed along the interaction horizon.

A special case of modeling using hybrid methods is the development of macro-elements in the nonlinear dynamic behavior of different types of foundations. Macro-elements are connection elements, which are located at the foundation level and have a constitutive law to describe the nonlinear mechanisms developed along the soil-structure interfaces, such as soil sliding, debonding, and plasticizing.

What was said previously suggests that within the approaches to solve DSSI problems, the FE method is mainly used in the families of direct and hybrid methods for near-field modeling.

In the following paragraphs, a methodology is given to subdivide the structure and some issues are discussed regarding the implementation of direct methods for DSSI problems.

More about the sub-structuring method

The first step of the sub-structuring method is setting the viscoelastic parameters of the geomaterials according to the level of deformation. The hysteretic non-linear behavior of soil under cyclic stresses can be approximated by an equivalent shear modulus and viscous damping. Firstly, the curves of the reduction of the shear modulus and the increase of the damping coefficient as a function of the deformation have to be defined according to the results of laboratory tests specific to the project (part of the geotechnical investigation) or, in the absence of a dedicated geotechnical campaign, based on scientific references. Then, a linear 1D analysis of the response of a soil column sample must be conducted using transfer functions which, from the acceleration at the rock and the first set of soil layer properties, gives the displacements, velocities, accelerations, stresses, and shear deformations in each layer. Finally, an iterative procedure is required to ensure that the properties used in the analysis are compatible with the maximum deformations calculated in the layers. The inverse analysis, i.e. the calculation of the rock motion from the surface acceleration, is possible through the same transfer functions and is called deconvolution. The iterative procedure does not always guarantee the uniqueness of the solution in the search for equivalent soil properties in case of large deformations. These analyses can be carried out with dedicated software such as SHAKE 91 for example.

Diffraction and impedance problems can be solved in the frequency domain by adapted calculation tools such as SASSI2010 or MISS3D codes. Alternatively, simplified methods valid under certain conditions can be used.

The study of kinematic interaction is of interest in the case of deep foundations but is neglected in the case of shallow foundations. For a seismic analysis, the simplified method consists of the imposition of free-field soil deformations on the foundation piles. This stress, of kinematic origin, must be added to the stresses resulting from the inertial interaction for the dimensioning of the piles.

Using a dedicated software (e.g. SASSI2010, MISS3D), the study of the inertial interaction leads to the calculation of the stiffness and damping at the monitoring point of the foundation as a function of frequency. For the calculation of impedance functions, the foundation is by default considered a rigid substructure. For a shallow foundation, the impedance terms relative to displacements and rotations are used, unlike for deep foundations, and the coupling terms between the displacements and rotational stiffnesses must also be considered. The impedance values, retained for the structural problem of dynamic analysis, are set based on the eigenmodes of soil-structure interaction with an iterative procedure explained below. A first modal analysis is conducted with a set of test impedance values. The next set of impedance values is calculated from the frequencies of the eigenmodes, identified as those of ground-structure interaction for each degree of freedom of the (supposedly rigid) foundation. The second modal analysis is used to correct the frequencies of the selected eigenmodes and so on. The procedure stops when the convergence of the eigenfrequencies is reached.

Simplified methods for shallow foundations generally give formulae for the direct calculation of the impedance. The best-known methods are those of Newmark-Rosenblueth, Deleuze, Eurocode 7, SETRA Guide, Veletsos, and Gazetas.

The main condition for the validity of the simplified methods is the assumption of homogeneous, single-layer soil, except for the Gazetas method that allows the consideration of a double-layer soil for the cases of foundations of approximate circular shape. Deleuze's method considers the frequency of the fundamental vibration mode of the structure. The methods of Veletsos and Gazetas consider the embedding of the foundation in the soil.

For deep foundations, the Gazetas method and the Eurocode 8-5 method give formulas applicable in the following three soil conditions:

soil whose Young's modulus varies linearly with depth,
soil whose Young's modulus varies with the square root of the depth,
soil whose Young's modulus remains constant with respect to depth.

Another simplified method for deep foundations is Winkler's method, which requires the modeling of a beam resting on a system of springs and dampers with properties based on those of each layer and representing the frontal pressures and frictional forces.

After solving the dynamic analysis problem, the Winkler method is also used to calculate the stresses along the piles for their design. This method makes it possible to model an elasto-plastic behavior for each soil layer and to incorporate a simplified analysis of the kinematic interaction by applying the soil displacements imposed on the pile.

Defining the mesh size and time step

The choice of mesh size is a very important step in soil dynamics modeling because the discretization influences the characteristics of the waves that can be modeled with a given mesh. The usual criterion to model a wave propagating in the soil with acceptable accuracy is to have in the mesh, 10 elements per wavelength. If we denote h the mesh size and λ the wavelength, the criterion to be imposed on the precision of discretization is written as follows:

This criterion can also be formulated as a function of the minimum propagation speed V_min that needs to be considered in the model and the maximum characteristic frequency f_max of the waves to be represented in the mesh. The criterion for the mesh size is then written as:

For seismic problems, the quantity Vmin usually corresponds to the minimum value of the effective velocities (compatible with the seismic level) of propagation of S waves for the modeled soil layers. The quantity fmax corresponds to the cut-off frequency of the dimensioning spectrum (~30 Hz) or a maximum frequency of interest for the ground-structure interaction (~15 Hz).

Apart from the mesh size h, the choice of a time step Δt is also conditioned by the mechanical characteristics of the modeled medium. Especially, by analogy with the wavelength, it is necessary to be able to describe the smallest period that must be considered in the analyses by at least 10-time increments (the precision criterion is identical in the spatial and temporal discretization). The criterion about the time step can be formulated as follows:

where V_max is the velocity corresponding to the maximum value of the (effective) propagation velocities of the S waves and h corresponds to the mesh size.

Note that for an applied problem where V_min = 150 m/s, V_max = 500 m/s and f_max = 15 Hz, the obtained values would be h ≤ 1m and Δt ≤ 0.002s. One must pay attention to non-linear dynamic problems because the characteristics of interest to define the mesh size (in particular V_min) evolves with time (softening of the soil by passing in elastoplastic regime).

Cyclic behavior of soils - Damping

There is a growing interest in the implementation of direct methods for DSSI problems concerning the possibility of rigorously modeling the non-linear behavior of geomaterials, at least within the domain encompassed by the interaction horizon, which is modeled with finite elements. The designer can accurately model the geometry of the ground layers and assign to each formation the most suitable constitutive law to reproduce the anticipated cyclic behavior.

Calibration of non-linear models - The first step in this modeling procedure concerns the calibration of numerical parameters involved in defining the chosen constitutive laws so that they accurately reproduce the cyclic behavior of the soils. The latter is characterized, in most Civil Engineering applications, from the degradation curve of the shear modulus G as a function of the distortion γ (G~γ curve). It should be remembered that the G~γ curve varies according to the type of soil (sand/clay/silt/gravel), the effective vertical stress, the plasticity index, etc. The procedure for the calibration of non-linear models includes the following steps:

organizing the mesh and showing the areas governed by the same G~γ curve,
defining the elastic parameters (low distortion parameters) for each zone,
adopting an adequate constitutive law for each zone (fracture surface, strain hardening law, flow law, behavior in the deviation plane and the hydrostatic axis, etc.),
defining the G~γ curve characterizing each zone as a function -among others- of its mean effective stress,
calibrating the numerical parameters that characterize the non-linear constitutive law to reproduce the G~γ curve of each zone.

The calibration procedure results in a table that provides for each zone, the type of constitutive law adopted, and the set of numerical parameters necessary to define it.

Drained and undrained behaviors – When a water table is present in the soil domain, it is important to define whether the response of the geomaterials during dynamic loading will be governed by drained or undrained behaviors. This definition depends on the nature of the soils (cohesive materials; pulverulent materials) but also on the duration of the studied phenomenon. For phenomena of limited duration (of the order of a few seconds, such as for earthquakes) it is often difficult to prescribe whether the soil will respond in drained or undrained conditions. For these situations, it is advised to carry out two separate calculations, both under the hypotheses of drained and undrained behavior, and to define the overall response as the envelope of these two calculations.

Damping - When modeling the behavior of soil layers via an elastoplastic law, the anticipated damping of the material is obtained directly via the dissipation deduced from the constitutive law. However, this definition of damping differs from its conventional characterization (usually accompanying the definitions of the G-γ curves) via a β-γ curve, where β represents the equivalent viscous damping rate. This difference comes from the fact that the type of damping is not the same in both cases: on the one hand there is the viscous damping β, on the other hand, there is the hysteretic damping directly derived from the elastoplastic law. Thus, it is inaccurate to calibrate the non-linear model according to the β~γ curve. However, it is important to verify afterward that the damping rates generated by the calibration of the elastoplastic model are not far from the values of β which come from the theoretical β-γ curves. In addition to the hysteretic damping, it is necessary to define an "elastic" damping, i.e. a basic damping rate to consider as soon as the load is applied. Defining this damping component (ranging between 1% - 3%) is useful for the proper conditioning of the dynamic equilibrium equation and can be introduced via Rayleigh damping. Finally, it is possible to introduce a third damping component, namely numerical damping, which results from choosing a dissipative integration scheme. This third component should be used with care and avoided for common applications.

Liquefaction - Special calibration of soil models is required in the context of liquefaction risk studies. In this case, a two-phase model is required for potentially liquefiable areas. As a result, it is possible to model the increase in pore pressure that occurs in these low permeability materials during earthquakes, a phenomenon that can lead to liquefaction. The calibration, in this case, consists of showing that a liquefied state is obtained on a soil sample that was consolidated in the state of stress of the potentially liquefiable zone when it is solicited by the theoretical earthquake leading to liquefaction. The most relevant parameters for calibration with respect to liquefaction are those governing the response of the soil loaded along the hydrostatic axis (volumetric parameters), such as the angle of expansion or the rule for defining plastic volumetric deformations, etc.

Boundary conditions: Absorbent boundaries - Floating models

Dimensions of the area to mesh - Another complexity concerning the implementation of a direct approach, concerns the "dimensioning" of the mesh and the definition of boundary conditions. First, the dimensions of the mesh must be sufficient to describe the phenomena involved. In the vertical direction, the mesh is usually extended to the bedrock level, where a linear behavior can be assumed. In the horizontal direction, the mesh is extended at a characteristic distance equal to 4 to 5B on either side of the structure, where B is the characteristic dimension of the structure. It follows that the minimum dimension of a mesh for the implementation of a direct DSSI method is 10 times the characteristic dimension of the structure. It is understood that this criterion, combined with the criteria for mesh size, presented above, can make the cost of a DSSI calculation by the direct method prohibitive.

Boundary Conditions - For problems where external stresses come from the structure, defining the boundary conditions consists of introducing absorbent boundary elements along the contour of the mesh whose role is to annihilate the parasitic reflections of the waves that reach the edges of the mesh. The absorbent boundary elements commonly used in practice are standard viscous dampers that are only accurate (as absorbent boundaries) for a given direction of propagation and wave frequency. To overcome these limitations, several improved absorbent boundary formulations are available. Nevertheless, their implementation remains relatively cumbersome and these formulations are still little used in current studies.

For seismic problems, the role of absorbent boundaries is to reduce the parasitic reflections and to introduce the incident seismic motion along the contour of the domain modeled in finite elements. To achieve this second objective, it is often necessary to determine the motion at the bedrock level (since seismic standards generally provide the dimensioning earthquake at the ground surface): this motion is used as the incident motion at the base of the model. Additionally, the motion (under free-field conditions) on a soil column (from the bedrock to the surface) must be determined, which will be used to establish boundary conditions along the sidewalls of the model. Furthermore, the mesh needs to be sufficiently large in the horizontal direction because the free field movement at its lateral boundaries must be expressed. For asymmetric meshes, the calculation of the soil column in the free field must be conducted separately for the columns on the left and the right side of the mesh.

It must always be verified that the calculated motion at a node close to the lateral boundaries of the model is identical or very close to the free field motion. This verification ensures that the mesh is correctly dimensioned for the interaction problem under study.

Phasing - As with all geotechnical problems, it is also necessary for dynamic problems to establish the true state of stress in the geomaterials before introducing dynamic loading. Usually, any DSSI studies are initiated by a gravity loading phase. At the end of this phase, the system displacements are set to zero and the stress states are kept throughout the entire mesh. At the end of the gravity phase, it is possible to introduce dynamic loading into the system.

Floating Models - Still for seismic problems, it is important to adopt a model that allows simultaneous injection of the horizontal and vertical components of the seismic motion. To consider the vertical component, it is often necessary to use the so-called "floating model" method, which consists of the following steps:

blocking the nodes at the base of the model and initiating under self-weight,
identifying the reactions of the nodes blocked under self-weight,
continuing the calculations under self-weight by vertically releasing the nodes at the base and simultaneously applying the gravity forces and the reactions calculated in the previous step: the model is not fixed at any point in the vertical direction (hence the name "floating model"), the equilibrium is ensured by setting the equality of the reaction forces and the gravity loads,
finally, after the gravity load, the earthquake is introduced. Moreover, having released the vertical nodes also introduces the vertical component of the seismic motion.

Choice of the integration scheme

The implementation of the direct method for the DSSI is completed with the definition of the time integration scheme of the dynamic equilibrium equation. The principles that apply to general non-linear dynamic problems are also applicable in the case of DSSI problems. As these are non-linear problems, the time integration scheme is coupled with a method for integrating non-linear constitutive laws, typically the Newton-Raphson method or the modified Newton-Raphson method, etc. The time integration scheme is then used as a basis for the integration of the non-linear constitutive laws.

Generally, problems of very short duration (impacts, explosions, etc.) are handled by explicit schemes. These schemes are very efficient if the mass matrix can be triangulated and are conditionally stable (stability criteria often lead to very small time-steps). For seismic problems, implicit schemes are used, the most used being the Newmark integration scheme (which is unconditionally stable in a linear regime for a given set of integration parameters).

The choice of a dissipative scheme (leading to numerical damping) can also be proposed to facilitate the convergence of the solution for highly non-linear problems.

In any case, the choice of an integration scheme, specific to a given problem, with the definition of all the associated parameters (tolerances, number of iterations, strategy to manage the divergence, etc.) requires a very good knowledge of the theory and also a significant experience in the typology of the problems treated (seismic, impacts, fast dynamics, vibration problems, etc.). As a general principle, it is advised to start the analyses with non-dissipative schemes and with relatively strict tolerances (of the order of 10-4). At the end of each analysis, the file containing the detailed results should be carefully examined to verify that the residuals of each increment do not exceed the tolerances and that convergence is obtained on all the meshes of the model.