PART 1 – THEORETICAL ELEMENTS

$translationBooks

Chapter A. Introduction

Chapter A. Introduction

When performing a FE calculation of a structure, the modeling assumptions as well as the objectives must be defined. Those goals generally consist of solving:

 Among the modeling hypotheses, it is preferable to handle chronologically:

A1. A general formulation for linear elastic calculations

A1. A general formulation for linear elastic calculations

A2. The dimensionality of the model

A2. The dimensionality of the model

A3. The choice of the FE

A3. The choice of the FE

A4. Interaction between the structure and its environment

A4. Interaction between the structure and its environment

A5. Estimation of the quality of the approximated numerical solution

A5. Estimation of the quality of the approximated numerical solution


$translationBooks

A1. General formulation for linear elastic calculations

A1. A general formulation for linear elastic calculations

The two main formulations for linear calculations are displacements-based and stress-based. Most of the software nowadays rely on a displacement formulation because it can easily consider the boundary conditions, unlike a stress-based formulation. Mixed formulations, using displacements and stresses also exist. There are other formulations, but they will not be discussed in this chapter.

In this document, only the displacement formulation (the most commonly used) is introduced: in this context, the basic unknowns are the displacements at the nodes, typically represented by q, from which all the other values are deduced (strains ε, stresses σ  ....)

The displacement at any point of an element, represented by u(x,y,z), are obtained by interpolation between the values at the nodes using interpolation functions (known as shape functions) represented by N(x,y,z), which are low order polynomials (generally of degree 1 to 3):



The strain-displacement relation is obtained by computing the derivative using the differential operator (matrix) D (different according to the situation,) from which the matrix B appears, connecting strains and displacements to the nodes:



To obtain the stresses, it is essential to introduce the constitutive law of the material. Thus, in linear elasticity, the strain-stress relation is written (H representing Hooke’s matrix):


One can observe that all the interesting quantities (displacement at any point, strains, stresses) can be obtained from the displacements at the nodes of the mesh.

Note: as the deformations (and consequently the strains) are obtained by computing the derivative of the displacements at the nodes, there is a loss in accuracy when calculating those values.

For a static calculation in mechanics of structures, identifying the numerical solution of the problem by the Finite Elements Method (FEM) can be summed up in 4 main steps:

Step 1. Determining the elementary matrices and vectors of each element of volume V by applying the following equations (with BT representing the transposed of B):

 elementary stiffness matrix

 elementary vector of equivalent loads 

with fV the volumetric load vector and fS the surface load vector. The vector of equivalent loads allows us to “bring back to the nodes” the loads applied to the elements (in the volume or the surface,) according to the chosen interpolation.

Be careful! Usually, it does not correspond to an even distribution at the nodes of the total load applied on an element. 

Step 2: Determining the stiffness matrix K and the force vector F with respect to the global coordinate system by assembling the elementary quantities and introducing the boundary conditions.

Step 3: Determining the vector of nodal displacements q by solving the linear system (large scale):


The vector q gives the node displacements of the mesh, likewise, the vector F contains the values of the equivalent loads. 


Step 4: Resolution of the quantities of interest such as strains and stresses by post-processing the vector of nodal displacements q.

Going back to the elements, the displacement, the strain, and the stress values at any point can be deducted from the displacements at the nodes by using the previous equations. As opposed to the displacements, the stresses and the strains are calculated element by element (and these quantities are not always continuous from one element to another.)

Note: it is possible to make an analogy between solving, using finite elements, a mechanical problem, and a thermal problem in a linear permanent regime: in this case, the temperatures T at the nodes are obtained by solving a linear system: Λ. T = Φ, with Λ representing the conductivity matrix and Φ the flux vector at the nodes.  

The order of the shape functions N(x,y,z), which is deducted for the solid elements (or continuum mechanics – CM) from the number of nodes in the edge of a finite element (2 nodes: linear element, 3 nodes: quadratic element) has an impact on the degree of strains and stresses evaluated using finite elements. For example, if the differential operator D corresponds to the first-order derivative (which is the case in CM,) the strains and stresses will be constant along with the element for linear elements, and linear along with the element for quadratic elements… As a reminder: in Strength of Materials (SM) to compute the strains we make use of the second-order derivatives (bending of beams, plates….)

Shape functions – Triangular element with 3 nodes (linear functions)

Shape functions – Triangular element with 6 nodes (quadratic functions)

Be careful! It is not possible to mix in the same mesh finite element with shape functions of a different order (for example linear and quadratic elements.) Actually, the continuity of strains cannot be ensured at the interface between elements with shape functions of a different order.

In the software, most finite elements are defined with the same parameters (isoparametric,) meaning that the different finite elements of a mesh are built as a geometric transformation of a reference element (which is itself expressed in unitary coordinate system ξ,η,χ.)


Isoparametric element (quadratic): different transformations (2D)

These geometric transformations, which generally rely on the same nodes used for the interpolation of the displacements (hence “isoparametric”, with the same parameters,) are characterized by their Jacobian matrix J (matrix of the geometric transformation’s partial derivatives of the first order.)  

Isoparametric element (quadratic): local and global coordinate systems (2D)

One of the advantages of this transformation is the possibility of approximating the boundaries of a given domain by using polynomial geometries. However, the main problem is that the calculation of elementary matrices and vectors can become more complex. Their evaluation is made by numerical integration using a quadrature scheme: weighted summation of the values at a certain point located inside the elements (weights ω_i, points, ξ_i) particular to the method (Gauss method, Hammer method….) In this case, the calculation of the elementary stiffness matrix and the elementary vector of equivalent loads are written:

When the matrices (by applying the B matrix) are evaluated this way, it is obvious that one must calculate the strains and stresses at the same points (also as a function of B.) The software provides these quantities at the integration points (Gauss or Hammer points, according to the integration scheme.) Moreover, talking about the values at the nodes of a mesh for strains and stresses does not make sense because there are as many values in a node as the number of elements connected to it. Be careful, there is not necessarily a continuity between the elements! Wanting at all costs a nodal value (of strain or stress) of the mesh implies some data processing: extrapolating from the values at the integration points and the average values of different elements sharing the node, which will, at the very least, smoothen the results (nonetheless, talking about values at the nodes in a particular element is justified.) When the quantities are obtained by analytical integration (in the case of beam elements for instance,) they can be expressed at the nodes, at the center of gravity… the previous observations about the continuity between elements being still valid.

Quadratic elements and integration points (Gauss points)

The way to consider the displacements at the boundary conditions (restraint, support…) might lead to different results depending on the adopted numerical technique.  There are two ways of considering the displacements at the boundaries:

Note: in dynamics, the displacements at any point x,y,z at any moment t is obtained, as in statics, by interpolating the values at the nodes: ux,y,z=Nx,y,z.q(t). These nodes verify the equilibrium condition, given at any time t by:



The mass matrix can be obtained by assembling the elementary mass matrices MV (volumetric mass) and C the damping matrix (viscous). To compute the aforementioned matrix, a system of differential equations needs to be solved using adequate numerical methods (see the chapter on dynamics.)

$translationBooks

A2. The dimensionality of the model

A2. The dimensionality of the model

It is very important to simplify the full-size problem to model the interaction between the structure and its environment. To do so, there are 2 known techniques:

Those 2 techniques help reduce the computational cost. However, the user has to be very careful when adopting those simplifications. For instance, the second technique uses hypotheses that are valid only in a limited domain. It is important to make sure that the adopted simplification stays within the domain of the full-size problem to obtain relevant results.

1) Case of finite elements applied to Strength of Materials

From a finite elements method perspective, the main difference between Continuum Mechanics and Strength of Materials relies on the geometry, which is simplified through additional hypotheses: the initial 3D problem is then transformed into a 2D problem (mean surface for plates and shells) or 1D (mean fiber for beams), but still represented in a 3D space (except for 2D problems, as mentioned above). When the FE method is used to solve a Strength of Materials problem, the finite elements have specific characteristics. It is necessary to assign to those elements geometrical properties (the section, the inertia for beams, the thickness for plates and shells). Moreover, they combine the effects of tension/compression (for beams) or membrane effects (for plates and shells) to the bending effects. The traction/compression and membrane effects are treated as they are treated in Continuum Mechanics, while the bending effects are handled separately.

For the computation of bending effects, simplifications in the geometry lead to a particular definition of deformations, which results in different expressions of the differential operator D according to the situation. Another consequence is due to the definitions of the unknowns at the nodes: while in Continuum Mechanics the unknowns are the displacement components, in Strength of Materials the unknowns corresponding to the rotations are added. This happens because it is no longer possible to evaluate directly the rotations as a consequence of the simplifications of the geometry. The initial choice of treating the problem in the frame of Continuum Mechanics or Strength of Materials implies choosing between different families of finite elements. It is, in principle, not possible to mix finite elements of different families because at the interfaces between elements of different nature the rotations are not transmitted unless they are adapted to take it into account. Moreover, the calculated stresses inside the element are usually generalized stresses (or Strength of Materials efforts: axial, shear, torsion, bending). To obtain the stresses (in Continuum Mechanics) at a given point, it is necessary to provide additional information (position inside the beam’s cross-section, for example).

In this scenario, choosing a finite element brings on an additional difficulty related to the consideration or not of the shear forces (Euler-Bernoulli or Timoshenko beam theories, Love-Kirchoff, or Reissner-Mindlin plate theories). This choice is related to geometric considerations (slenderness of the cross-section or thickness of the plate). Moreover, taking into account the shear forces might lead to numerical errors (shear locking), which renders the use of some finite elements complicated.

The Euler-Bernoulli beam finite element allows to represent exactly the bending moments that vary linearly along with the mean fiber of an element (shape functions of order 3 for the bending displacements and the bending moments obtained by taking the second derivative of the displacements): it is therefore not useful to introduce several elements between 2 nodes subjected to point loads (a dense mesh should be used when the loading is distributed between these 2 points).

Finally, for plate and shell elements, the monotonic convergence is not guaranteed depending on the mesh geometry. This behavior is connected to the formulation of the element itself. This kind of behavior is illustrated in the figure below, showing the evolution of the relative error of the deflections ω and Mx as a function of the number of degrees of freedom (in the case of a rectangular plate supported by its four sides submitted to uniform loading): in green, a non-compliant element (COQ3) for which the results are much less precise than for a compliant element (DKT). Also, for a non-compliant element (COQ3) the convergence is not monotonic (for the bending moments Mx in the figure): more elements leads, paradoxically, to results that can be less precise!

Bending plate: convergence (deflection ω and bending moment Mx) as a function of the number of degrees of freedom

This domain is particularly important to civil engineering, but it presents specific difficulties that will be explored in the next chapters.

2) Two-dimensional calculations

The studied problems are three-dimensional, however, it is faster to calculate two-dimensional problems. In some cases, it is possible to transform a three-dimensional problem into a two-dimensional problem:

Cylinder under axisymmetric pressure

Dam: important thickness of plane strain assumption

Assembly: low thickness/plane stress assumption

3) Taking into account the symmetries

Some problems present symmetries (axis of symmetry, plane of symmetry…) and it can be interesting to use them to make the finite element calculations faster. Although, it is important to notice that solving symmetric problems will provide only symmetric results (especially when computing eigenmodes).

For that matter, it is important to keep in mind that to use those symmetries in the calculations, not only the geometry has to be symmetrical, but also the loading and the boundary conditions as mentioned above. The resolution of this model will not represent the solution of the whole structure unless the symmetry conditions are added to the model. In Continuum Mechanics, one should impose the displacement components perpendicular to the axis/plane of symmetry equal to zero. In Strength of Materials, it is convenient to add no-rotation conditions around the axis/plane of symmetry. Finally, taking into account symmetries might impact the loading: for instance, a punctual load should be applied with half the intensity if applied to an axis of symmetry.


$translationBooks

A3. Choosing the finite elements

A3. Choosing the finite elements

Choosing the element is an important step. The goal is to select the type of element (its shape, the order of the shape functions associated), and the size of the element. The type and size of the elements define the form and the precision of the displacement fields, and consequently the stresses and strains.

Moreover, besides the shape, the element has an aspect. It should be avoided to create decayed or altered elements (flattened or elongated) because it degrades the precision in the resolution of the problem.

Generally, the generated mesh can be structured (a regular division of the elements) or unstructured. It is possible to mix structured and unstructured parts in the same domain depending on the complexity of the geometry.

In the same context, the size of the element chosen depends on the geometry of the structure to model and the loadings it is subjected to. The zones with high-stress variations (gradients) or high stresses (friction or cracking, for example) determine the parts where the mesh needs to be more refined in comparison with the other parts to properly observe the stresses and the strains.

It is important to perform the first simulation with an initial meshing to determine the sensitive zones and then refine the mesh where needed.

The person in charge of modeling must analyze the meshing critically considering the geometry of the structure and the important zones to observe.

To choose correctly the finite elements construct the mesh, it is imperative to think about the type of calculation desired:

Besides the displacement degrees of freedom, the Strength of Material elements have also rotational degrees of freedom (θx, θy, θz), allowing to take into account the non-meshed geometry (cross-section for beam elements, thickness for plates and shells). Moreover, one should ask the question if the shear forces should be taken into account or not (Navier-Bernoulli or Timoshenko beam elements, Love-Kirchoff or Mindlin-Reissner shell elements). In the case of shell elements, as mentioned before, the question about the quality of elements comes to mind (compliant or non-compliant elements). It is particularly difficult to choose between plate and shell finite elements, especially because the literature is scarce in this area. It can be helpful to perform a calculation with known results to test the quality of the available elements.

The different elements that are commonly found in Strength of Materials are described in the table below.


$translationBooks

A4. Interaction between the structure and its environment

A4. Interaction between the structure and its environment

Taking into account the interactions between the structure and its environment: this consideration allows us to reduce the number of discrepancies between the model and the real-life structure. The effect of the environment on the structure is determined from the nodal values such as the nodal displacements and nodal efforts. The first ones are connected to the boundary conditions and the second ones to the external loading [1].

Boundary conditions for a plate with a hole in the center, 2 axes of symmetry, submitted to tension

Afterward, it is essential to eliminate the rigid body movements. A correct finite element model has to restrain free rotations. In 3D, the 6 rigid body movements previously introduced should be avoided (in 2D, there are 3 movements, figure 4).

Possible rigid body movements in 2D: a) possible horizontal translation; b) possible vertical translation; c) rotation around a hinge; d) all rigid body displacements restrained

Once these 2 steps are accomplished, it should be verified that every rigid body movement is correctly restrained and that no rigid body movement was blocked when this mode has already been removed. In the first case, the problem has no solution and in the second case, unexpected strains may appear.

There are also contact forces such as pressure or any other force that needs contact with the structure. They can be surface, linear, or punctual forces. Their application should also be transformed into nodal forces. Special attention should be paid to the transformation of the contact forces when creating a model respecting its validity domain. Using a point load can generate singularities such as stress concentration in the neighborhood of the node where the load was applied. Thus, to avoid this kind of singularity, it is necessary to apply the point load as a surface load in the neighborhood of the node. This technique consists of applying a surface load containing a relatively large area around the application point. Then comes the question of mesh refinement in this area and its influence on the obtained results.

Connection conditions

There are different types of finite elements, such as volumetric, plates, shells, beams, and bars. Thin plates and shells are elements, which have a thickness smaller than their other 2 dimensions. A plate works only perpendicularly to its plan (3 degrees of freedom (DOF) for a node: 1 translation and 2 rotations), whereas a shell element works along its plan and perpendicularly to it (6 DOF: 3 translations, 3 rotations). The solid elements also have 6 DOF per node. Chapter 2 presents some practical examples.

 

$translationBooks

A5. Estimation of the quality of the approximated numerical solution

A5. Estimation of the quality of the approximated numerical solution

The difference between the exact solution of the problem and the approximate one obtained through the finite element method allows us to assess the quality of the solution: this is called “discretization error”. Since in the general case the exact solution is not known, the idea is to estimate this difference by computing an “error estimator”. The error estimators can be categorized as:

The “error estimator” tools (available in some software, with different computational costs in each one) can be used for two different purposes:

It is important to highlight that certain error estimators, such as the stress smoothing method, only estimate the error, while others, such as the weighted residual method, guaranty the numerical results by computing the range of error.

 

$translationBooks

Chapter B. Structural Dynamics

Chapter B. Dynamics

For many applications such as seismic calculations, collisions, vibrations… it is necessary to consider the dynamic phenomena.  

The dynamic charges applied to a civil engineering structure belong to two different categories:

Concerning seismic motions, if they are theoretically considered as transitory, it is admissible to assimilate them as stationary phenomena while in their strong phase. For the cases in which one tries to model the structure with geometric or material non-linearities, stationary assumptions cannot be considered anymore. 

The means of representing the loading categories can then be distinguished as shown below:

Two big families of analysis can be considered:

The explicit schemes dictate the choice of very small-time steps. Thus, they are the most used to solve problems with small periods (like collision/impact problems.) On the contrary, implicit schemes allow us to use greater time steps and are therefore favorable to study problems occurring on wider time ranges. 

Examples of applications

Applications

Loading representation

Quantities available

Modal

Vibrations analysis

FT

FT

PSD

PSD

Tracking of natural frequencies

ORS

Spectrum extrema of quantities of various interests

Implicit transitory

Seismic Study

Accelerations, velocities, forces, pressures, or displacements as a function of time

Quantities of diverse interests expressed throughout time

Weakening 

Explicit transitory

Fall of on object

Modeling of projectiles in contact, collisions

Quantities of diverse interests expressed throughout time

Plane crash

Accelerations, velocities, forces, pressures, or displacements as a function of time

Once the dynamic problem is discretized in finite elements, the resolution of the equilibrium equation can be written as shown below (cf. chapter 1):

In a modal analysis, the natural frequencies ω_i and the associated modes φ_i are used.

In the time-dependent analysis, the displacements q(t), the velocities q'(t), and the accelerations q"(t) at the nodes are calculated for each time-step t by direct integration of the equilibrium equations. 

The second approach has an advantage as it allows us to handle nonstationary solicitations.

B.1 Analysis based on modal search

B.1 Analysis based on modal search

B.2 Analysis based on a direct temporal integration

B.2 Analysis based on a direct temporal integration

B.3 Considering the damping

B.3 Considering the damping

B.4 Particularities of a seismic analysis

B.4 Particularities of a seismic analysis


$translationBooks

B1. Analysis based on a modal search

B1. Analysis based on a modal search

A reminder of the concept of simple oscillator – Concept of the elastic response spectrum

In the following section, a lot of applied methods return the response of damped oscillators with a degree of freedom often named simple oscillator. 

Its general behavior is briefly reminded here. 

Its behavior is explained starting from the dynamic equation applied to an oscillator with a degree of freedom submitted to a time-dependent harmonic loading and building its transfer function with the following parameters: the ratio between its natural frequency and the natural frequency of the harmonic loading, and its critical damping ratio. 

The dynamic equilibrium equation applied to the mass-spring-damper with only one degree of freedom is reminded below:

It can also be written in its canonical form by dividing all the terms by the mass m such that:

If the harmonic loading is applied  xs=X.ejΩt, one can find the transfer function as an absolute response by establishing the ratio between the input loading and the absolute output response below:

The norm of this question is the amplitude of the transfer function that allows highlighting the dynamic amplification phenomenon illustrated for different values of the critical damping in figure 1. This transfer function’s independent variable is the phase.

Transfer function

A response spectrum is different from a transfer function. The curve giving the acceleration (named spectral acceleration) as a function of the period (or frequency) is called the spectrum of elastic response. The latter corresponds to the maximum absolute acceleration seen by a simple oscillator throughout time as a function of its natural period (or natural frequency) and its critical damping ratio. It dimensions the seismic motion. It is possible to construct an approximate relation between the Fourier spectrum of acceleration and the spectrum of spectral velocities for a damping ratio equal to zero. 

Starting from the dynamic equation of the simple oscillator subjected to an arbitrary acceleration: 

The maximum responses are deduced from its resolution according to time using the method of your choice (Duhamel integration, direct calculation…) to then calculate the relative spectral values of displacement, relative velocity, and absolute acceleration:

Une image contenant texte Description générée automatiquement

The concept of pseudo-displacement, pseudo-velocity, or pseudo-acceleration is often used, which means approximating these quantities, starting from one of them, based on the fact that the damping ratio is low (less than 20%).

Une image contenant texte Description générée automatiquement

A spectrum can be presented as shown in figure 2.

Spectrum

Concept of vibration modes

The structure’s response is based on a sum of harmonic responses, which are part of the domain of solutions q(t) of the equation:

Une image contenant objet, horloge, signe, table Description générée automatiquement

They can therefore be written as follows:

If the Basile hypothesis is used (C is neglected), the algebraic solution becomes:

If the Basile hypothesis is not considered, one is now interested in finding complex modes such that:

and the eigenproblem can be resolved:

Une image contenant objet, horloge Description générée automatiquement

To get there without neglecting C, one had to postulate that:

Une image contenant texte Description générée automatiquement

With the matrices A and B being:

  and    

The vector ϒ is created such that: 

Une image contenant texte Description générée automatiquement

The solutions presented as pairs of combined complex harmonics are calculated as follows:

with:


In any case, real or complex, one must solve an equation in λi2 of the same degree N (2N for complex cases) than the degree of the matrix N (2N for complex cases). The order of this matrix is equal to the discretized system’s number of degrees of freedom (twice as much for complex cases). To solve the mentioned equation one must find the cases in which the determinant of this matrix and its characteristic polynomials are equal to zero.

To simplify the problem, the case of real vibration modes assuming the Basile hypothesis is considered.

For each natural frequency ωi a modal shape vector is associated. The latter will be calculated by fixing one of the components to 1, then solving a system with n-1 parameters. 

Be careful! The modes are defined with a precision level of more or less a multiplying constant. Moreover, they are later normalized using various procedures. The most common for FE software being the normalization according to the mass matrix. This procedure is detailed in the following section. The modes can also be normalized using their greatest modal displacement.   

The vibration modes are orthogonal to the mass matrix, which implies that: 

Une image contenant objet, horloge Description générée automatiquement

When there are no loadings, the modes have no physical meaning. Different methods can then be used to normalize them to make their visualization more comprehensible. Thus, it is common in FE software to:

Une image contenant texte Description générée automatiquement

For the cases in which the software does not normalize this value, it can be denoted by the generalized mass term. It can then be written that for the mode i:

The objective of this document is to explain case studies for which FE methods were applied, so it is important to remember that the vibration modes are normalized according to the mass matrix, and thus the generalized mass is always equal to 1.

The physical value called modal mass is also created, which will allow identifying the structure’s mass quantity that is dragged by a mode in a given direction Δk:

Une image contenant texte Description générée automatiquement

Please note: modal mass is a very different quantity from generalized mass, as the above equation shows.

Use of the modal basis

Time-dependent response by projection onto the modal basis:

A particular solution to the global equation: 

can be broken down on the bases of orthogonal modal shape vectors into N independent problems that are describing the response ri(t) of a simple oscillator:

In the case of earthquakes, the loading f(t) is an inertial loading applied to the entire structure: 

Harmonic Analysis: 

For this type of analysis, the goal is to construct a function describing the response of the structure according to the frequencies at different nodes of the model. This response can be a quantity of arbitrary interest (displacement, velocity, acceleration – absolute or relative -, efforts…) as a function of a variety of input harmonic loadings that can also be described in various ways similar to the output quantity of interest. 

In this case, the imposed loading on an arbitrary set of nodes are of the harmonic type: 

The result of the analysis at all nodes of the structure is a complex transfer function whose norm (amplitude) and argument (phase) are the most generally used. They enable us to identify the resonant frequencies of a structure of equipment and to know the amplitude of its response according to different harmonic loadings. 

Modal basis truncation – General case

Usually, because it is not possible to extract all the modes using calculations (too computationally costly, even numerically impossible in some cases), one must focus on those that are most likely to respond to the loading of the structure. 

If the representative loading frequency saved is fc (Hz), one must extract the modes up to 2fc.

It is also important to make sure the modes inclined to contribute locally are well represented. 

A search for vibration modes must not be accompanied by interpreting the frequency, but by analyzing the modal shapes. 


$translationBooks

B2. Analysis based on a direct temporal integration

B2. Analysis based on a direct temporal integration 

Schemes of integration

Main Principles:

The principle of the different direct integration methods is to divide the studied interval into n smaller intervals of length ∆t = T/n and to verify the equilibrium at the discrete times Ti = i ∆t = i T/n. The differences between the methods (centered differences, Wilson, Newmark…) rely on the hypothesis that is considered about the variation of the kinematic values on the interval ∆t.

Explicit scheme:

If the displacement value can be calculated directly from the value (or values) of the previous step, the method is said explicit. In this case, the equilibrium is considered at the beginning of the increment.

Implicit scheme: 

The implicit methods must consider the equilibrium at the end of the interval, which requires the resolution of a linear system. 

As stated earlier, there exist many schemes of integration. The objective here is not to make an exhaustive description of all of them (see the specialized publications for that). Only two schemes will be presented in this guideline, the centered differences method: an explicit and conditionally stable method, and the Newmark method: unconditionally stable (for linear problems).

A scheme is unconditionally stable if, for any initial condition and time step ∆t, the solution is bounded, particularly when ∆t/T is large. On the contrary, a scheme is conditionally stable if the obtained solution is bounded only if ∆t is smaller than a critical value ∆tcrit.

Precision is a concept that is different from stability and it is extremely important for unconditionally stable schemes. Besides the inevitable rounding errors, the precision has an impact on two sources of approximations: an (artificial) increase of the period and a decrease of the amplitude. The influence of these two phenomena gains importance with the increase of ∆t, but independently. 

Centered Differences Method:
It is a scheme resembling finite differences. It relies on the approximation of the acceleration (Taylor series of the second order):

To obtain the same order of error for the velocity, one uses:

The displacements at t+∆t are obtained considering the equilibrium at t:

Thus, introducing the approximations of the acceleration and the velocity:

Of the form: 

This method requires a starting procedure to calculate q(-∆t) (from the equilibrium at t=0). The damping introduced by the scheme is equal to zero (no decrease in amplitude). 

Newmark method:

The scheme is based on the following approximations of the velocity and the displacement at the end of the intervals:

From the values of both parameters α and δ (ranging between 0 and 1) depends on the precision and the stability of the method: the values (δ = 1/2 and α = 1/4) lead to an unconditionally stable scheme named the Newmark method: it corresponds to considering the mean constant acceleration. Considering the equilibrium at the end of the studied interval (at t+Δt), one obtains: 

with 

and 

the coefficients being defined as follows:

The scheme also requires a starting procedure: the value of q¨(0) is obtained considering the equilibrium at t=0. Alike the centered differences method, the Newmark scheme with basis (δ = 1/2 and α = 1/4) does not introduce numerical damping. The values (α = 0, δ = ½) enable to get back to the centered differences method. 

Choice of the spatial and temporal discretization

Criteria about the element's sizes that satisfy the wavelengths

Different criteria can intervene. They are related to the precision of the expected results and the type of calculations used. 

For the transitory analysis, the general recommendation is to have a range of 8 to 10 elements per wavelength. The stationary waves being composed of a sum of propagation waves, the stationary dynamics analysis (in the case of earthquakes typically) will be using the same type of criteria. 

As a reminder, the general equation relating wavelength λ, frequency f, wave celerity c, and natural frequency ω is:

According to the type of elements used and the type of waves of interest, the following formulations are to be remembered: 

Volumetric elements in an isotropic environment (typical case of the modeling of an elastic soil od with Young’s modulus E, Poisson’s ratio , and density ρ):

Shell isotropic elements:

These formulations can be written in a more general form for anisotropic media with 2 principal directions 1 and 2:

One can find the size required to target a cutoff frequency of 40Hz for shells, plates, and volumetric elements. The size required is given by the numbers below / λ for 10 elements. 

Types of modeling

Target Frequency [Hz]

Thickness [m]

Type of wave

Size coeff. of the element

Shells and blocks – GC

Concrete

40

0.2

Flexural

0.59

40

0.25

Flexion

0.66

40

0.5

Flexion

0.93

40

1

Flexion

1.32

40

1.5

Flexion

1.61

40

-

Shear

6.04

40

-

Compression

6.75

Shells and blocks – Meca

Steel

40

0.01

Flexion

0.16

40

0.02

Flexion

0.22

40

0.1

Flexion

0.50

40

-

Shear

8.02

40

-

Compression

9.58

This criterion is not necessarily sufficient. In the cases where a temporal analysis by transformation to the modal bases, or a spectral response is sought, a sensitivity analysis must be conducted to verify if the refining of the mesh leads to significant variations in the modal participation factors (or the modal mass if this criterion is examined). These elements are described later in the chapter, and their relations are expressed.

The criterion for the time steps

The conditionally stable temporal schemes, like the scheme of centered differences, must satisfy a condition concerning the time step chosen. This condition is often called CFL (Courant-Friedrichs-Levy). 

Example of application for a model of steel truss whose smallest element has a length L:

Let us consider a bar element with two nodes 1 and 2.

The bars have a length L, a section A, a density ρ, and a Young modulus E.

The mass matrix is considered as the total mass uniformly distributed at both ends. 

The matrices K and M are expressed as follows:

The characteristic polynomial is expressed such as:

Which leads to:

Since the propagation velocity of compression waves in a continuous medium is expressed as follows:

One can find that:

Une image contenant objet, horloge Description générée automatiquement

In this model whose mesh is irregular, it is the smallest time step that controls the global time step. 


$translationBooks

B3. Considering damping

B3. Considering damping

The damping of a structure’s vibration is related to dissipation phenomena of various origins that can happen inside of a structure. It can be due to viscous phenomena specific to fluid dampers, which can be found in spring boxes or dampers (cable insulators, Jarret™ type dampers, viscous fluid dampers). It can also be caused by structural dissipation phenomena related to cracks and friction at the interfaces between materials. The energy is consumed in the behavioral hysteresis of the materials. 

In most cases, the damping is considered in the dynamic calculations through a viscous damping matrix C. The components of this matrix can be exact at the local scale in the cases where the damping is completely controlled because the damping equipment used has well-known characteristics. However, most of the time, for civil engineering structures, the damping components are unknown because of the heterogeneity of the structures and the non-viscous characteristic of the materials. 

It is preferable to build the damping matrix using critical damping ratios as defined in the regulation for various materials. Those ratios are defined using a viscous behavior hypothesis proportional to the solicitation velocity. It represents a fraction of the critical damping value that leads to, for a simple oscillator, a return to equilibrium without going through an oscillatory regime.

For the cases in which the damping components are known, the matrix C is constructed explicitly, and it is necessary to use a suitable temporal calculation method to respect its characteristics:

The analysis of rapid dynamic phenomena (explosion, impact) requires specific considerations. 

It is not necessary to consider damping for local resistance verifications or vibrations over short periods, not long enough to “activate” harmonic responses of the structure. It is then not necessary for those cases to consider damping. It is interesting to note that this kind of calculations, realized with dynamic algorithms applying explicit time-dependent integration schemes (LS-Dyna, Radioss, PAM-Crash, Europlexus…), do not use the classical approach, which consists of inversing the stiffness matrix. On the contrary, for vibration analysis induced on a period following the impact duration, it is necessary to consider and characterize the damping. 

Note: for rapid dynamic calculations (milliseconds), the damping component is often neglected because its order of magnitude is much smaller than the one of inertial local stiffness terms during the period of analysis

For the more “classic” cases, the damping matrix is generally constructed from the critical damping ratios defined for materials and a linear combination of the M and K matrices that were previously calculated by the FE computational model. This process enables us to considerably simplify the stages of calculation. This approach is related to the “Basile” hypothesis for modal analysis. This hypothesis was formulated to identify the structural damping (mainly in aeronautics), and is explained as follows:

“Even when modal damping coupling is present, the modal equations of motion are dynamically decoupled, for structures with a low damping ratio, if the separation of the modes in frequencies is satisfied”.

Be careful! Considering the Basile hypothesis leads to an underestimation of the dynamic response elsewhere than at the edge of the materials or at the supports, which are subjected to a high damping ratio compared to a time-dependent calculation by means of direct integration, or using more complex methods such as the complex modes method. This justifies the need for direct integration calculations of the motion equations, or complex modes. 

Most of the time, the Basile hypothesis is considered and the matrix is constructed starting from the Rayleigh method: 

The terms α and β must be introduced. To be pragmatic, the following steps must be considered:

Une image contenant texte Description générée automatiquement

Une image contenant texte Description générée automatiquement

f1 is a lower bound of the first significant vibration mode of the structure, fn, and the first mode encountered after the cut-off of the elastic response spectrum characterizing the applied temporal signal. The value of the reduced damping is calculated for all values of ω, applying the following formula: 

Be careful! Since the damping is largely overestimated after the 2 frequencies of interest, it is necessary to master these two values. If fn is chosen adequately, the overestimation of the damping has no consequences because rigid body motions are found on these ranges of frequencies, which will be excited in the response, and the modes with frequencies greater than the critical frequency are insensitive to damping. However, one must be certain that there are no local vibration modes with frequencies worth considering.

The impact of an error in the choice of the first frequency can have the worst consequences if a significant mode was forgotten because its response will then be neglected.

It is important to note that in between the two frequencies, the damping is underestimated. The evaluation of the response of the structure is therefore conservative.  


$translationBooks

B4. Specificities of the seismic analysis

B4. Specificities of the seismic analysis

Spectral response – Specific case of earthquakes

The principle of the method is, for a given seismic direction, to construct the maximum responses from the loading spectrum at all points, mode by mode, then, to accumulate all of them using various methods. 

When the seismic response is obtained for a given direction, the seismic directions are combined to get the global response. 

Let us be in a case where the modal vectors are normalized according to the mass matrix, which is the case for most of the FE computational models.

One can write the expression for a degree of freedom in a canonical form:

The term pik is called the weighting coefficient. It is an important concept in the case of seismic analysis because it restores the physical contribution of excitation in a given direction k. 

It is determined for the cases in which the modal shape vectors are normalized according to the mass matrix:

Moreover, for a seismic motion q¨s(t) in a given direction Deltak, the mean solutions greater than the maximum values of r¨ik(t), r°ik(t), and rik(t) in terms of pseudo-accelerations, pseudo-velocities, and displacements via the respective spectra of this motion, are known:

This response is constructed for a modal damping ξi that can be evaluated as a proportion of the modal strain energy (see §1.2.1 of Génie parasismique Tome 3, J. Betbeder-Matibet).

with:

Then, from the simple oscillator response, one can construct the response in the structure for each mode, for the acceleration, the velocity, and the displacement respectively:

These modal values must then be combined for a given loading direction to obtain the global response of the structure. There exist several cumulation methods (from now on, ak,p, for the sake of example, denotes the component p of the ai,k vector, such as φi,p with φi) as shown below:

with Pij denoting the quadratic coupling coefficient of the modes i and j, such that:

It is important to note that, in this case, it concerns a cumulation of algebraic terms, but the signs are all positive, placing it under the square root don’t cause any problem. 

Let us consider the modal damping values ξi and ξj equal to a variable such that x = {0.2, 4, 7, 20 et and 30%}, one can graph, for the natural frequency ratios ωi/ωj, the curve (figure 3) that highlights the cumulation of the neighboring modes. The cumulation coefficient decreases rapidly when the frequency ratio increases, and even more when the damping ratio is low. 

Graphic of the coupling coefficient for different damping ratios varying as a function of the natural frequency ratios.

This last observation shows that the TQC method is more appropriate than the SRSS because the latter does not consider the cumulation of modes close to one another. 

On the contrary, there are no advantages of using the TQC method as a comparison to the SRSS method when the nodes are far from one another. 

Be careful! The cumulated modal responses must not be used to calculate other quantities of interest. For instance, the efforts deduced from the total quadratic cumulation method cannot be calculated using: : this calculation is incorrect in comparison to the evaluation by the cumulation of the modal efforts. 

One must, therefore, for this special case, calculate as follows:

1.
2. For each term p of the fk vector, it can be stated that:

It is also important to keep in mind that criteria allowing to evaluate the state of cracks based on the invariants of the stress tensor, and thus, of the principal stresses, (cf. for instance the criterion proposed in the annex LL of the EN1992-2) cannot be considered after using a quadratic cumulation of any kind.

The spectral responses of each mode, after being cumulated for one seismic direction, must be combined to obtain the global response. In this case, it is called spatial cumulation. 

The spatial cumulations can be handled using different methods:

Truncation of modal bases – Seismic case

For the specific case of earthquakes, and normalizing the modes according to the mass matrix:

It corresponds to the mass driven by the mode i in the direction k. This mass is therefore related to a specific direction. For a mode i, there will be 3 modal masses as a function of the different loading directions (if we are in a 3D space, 2 if we are in 2D, and only one for a 1D problem). Figure 4 highlights a 2D example for a frame. 

The modal displacements multiplied by the weighting coefficient represent the deformation of the structure whose product at all points with the response of a simple oscillator gives the time-dependent response of a given model. This concept is detailed in the following chapter. 

Mode 1 and displacement in the ey direction                                          Mode 1 and displacement in the ex direction

Illustration of the mass displacements for the same nodes in 2 different directions

As explained earlier if the model contains N degrees of freedom, there will be at most N natural modes of vibration.

However, since the search for natural frequencies and vibration modes is conducted numerically, not all the modes are extracted. In theory, one must have:

In practice, the algorithm will stop running after a threshold value of the frequency fixed by the user is reached. Then, the user must verify that there are enough modes to restore the percentage of mass required by the engineering rules followed. Usually, the criterion is to restore 90% of the mass. 

Be careful! If at the threshold frequency, the percentage of mass targeted is not reached, it might be necessary to include a pseudo-mode. 

Conversely, if a significant percentage of mass is reached at a low frequency, way below the threshold frequency, it is necessary to incorporate a pseudo-mode or to complete the modal base. Indeed, the participating mass of a floor for a local flexural mode inside a sizable structure represents a very small percentage of the total mass. 

One should also be careful concerning the antisymmetric modes with a modal mass that can be equal to zero because the masses displaced around one axis of the structure will balance one another. 

Pseudo-mode or static correction

As formulated before, the selection criterion of the natural modes concerning the cumulation of the mass is equal to 90% for a frequency of the same order of magnitude as the seismic response’s threshold frequency i.e. 40Hz maximum. 

Whenever this value cannot be reached, the lowest percentage withheld will be completed with an additional mass associated with a “pseudo-mode” of vibration. 

Let us remember that for a seismic direction k, there is an accumulated response on n modes such that:

To be rigorous, if there are N degrees of freedom, the formulation becomes:

As seen previously, the modes after the threshold frequency are rigid body motion modes. The structure reacts in phase with the seismic loading that it is subjected to with a relative displacement equal to zero. 

To complete the modal base, one can construct the pseudo-modes considering that the total response is the sum of the “dynamic” response, taking into account the modal base of the n first modes withheld, and the term proportional to the seismic acceleration of the support q̈s(t). Let us remember the following equation:

with Pk the displacement of the structure in the direction k subjected to the static loading equivalent to the mass of the accelerated structure at the acceleration of the last mode n extracted from the natural frequency:

Since the relative acceleration response for the structure is equal to zero after the threshold frequency, one can simply write:

Reformulating the previous equation as a function of , it will be easier to link the pseudo-accelerations spectrum, considering the frame of spectrum analysis:

Thus, the following formulation can also be considered:

After evaluating these vectors, they must be cumulated.

Indeed, the temporal solution for a seismic direction is given using the summation of the components on the modal base:

However, this combination cannot be applied to the maximum values that were just calculated previously. 

 

$translationBooks

Chapter C. Static non-linear calculations

Chapter C. Nonlinear Static calculations

Nonlinear calculations are generally time-demanding: a finite element software solves a nonlinear problem as a series of linear problems (incremental process). To prepare this type of calculation the user has to make decisions such as defining the increments, choosing the algorithm, etc. Therefore, the process requires a certain amount of experience from the user.

Most of the nonlinear computational models can, in principle, run until one or multiple conditions are satisfied, or in the case of contact problems, until a specified minimum number of connections (or supports) disappear. Therefore, before starting a nonlinear calculation, it is important to make a preliminary estimation, which will enable us to know when to stop the calculations.

C.1 Nonlinear mechanical problems

C.1 Nonlinear mechanical problems

C.2 Why performing nonlinear calculations?

C.2 Why performing nonlinear calculations?

C.3 Implementation

C.3 Implementation

C.4 Convergence problems? Symptoms and solutions

C.4 Convergence problems? Symptoms and solutions


$translationBooks

C1. Mechanical non-linear problems

C1. Mechanical non-linear problems

  1. Description of possible non-linearities

The non-linearity of a mechanical problem comes from the fact that the coefficients of the equilibrium equation depend on the displacements of the solid itself in equilibrium. In other words, the equilibrium equation is generally implicit.

There are several non-linearity categories for mechanical static problems:

All the above non-linearities can be coupled if the algorithm allows it, but the resolution of the problem becomes more complex.

  1. Principle of resolution of a non-linear problem: Newton method

When solving a finite element problem, one looks for the displacement field u, such that the inner forces Lint are equal to the external forces Lext:

, which is a non-linear problem as a function of u.

Generally, to solve the non-linear static problem, an incremental algorithm is used. For that matter, the problem is parameterized in terms of t (with t representing a pseudo-time, unlike the t parameter used in dynamics). This parameter is used to index the successive load-steps applied on the structure. More precisely, it consists of searching for the equilibrium state corresponding to the successive load-steps F1, F2, …

This separation leads to solving a series of quasi-linear problems as shown in the figure below and to determine the state of the structure at the time-step t (displacements, strains, stresses) knowing the solution at the state t-1. The greater the number of load-steps, the better the precision.

time

 increment.    

Principle of parametrisation in function of t


At each increment ti the discrete problem is Ki x qi = Fi where qi is the unknown displacement vector under the applied imposed loading Fi. While in the linear case seen in chapter 1 the K matrix was explicit, when the problem is non-linear, Ki is a matrix with its terms depending implicitly on the value of qi. So, qi cannot be determined directly by computing the inverse of the matrix K.

The most used method to solve this non-linear equation is to use a Newton-type algorithm. The idea is to build a good approximation for the equation’s solution

considering its first-order Taylor expansion


One must start from an initial point (close enough to the solution) and then compute by iterations

 

At each iteration, one should evaluate the residual vector F(qk) until it exceeds (in absolute value) the value arbitrarily close to zero. This convergence criterion must be chosen with care, respecting the standard used by the calculation code (see section 3.3 for more details).

Note: With the Newton method, at each iteration, one should compute the tangent matrix at the considered point:

The computational cost of this matrix can be time-consuming. If using this matrix allows having a quadratic convergence (so, in fewer iterations), it is not essential to use this matrix. Other strategies can be adopted to estimate this matrix, namely the quasi-Newton methods. It is conceivable to use the tangent matrix without updating it at each iteration, but also to use the elastic matrix (figure b) or the secant matrix in the case of a damage model. An illustration of the successive iterations according to the used matrix is shown below.

Illustration of the Newton or quasi-Newton method (elastic matrix)

In general, using the tangent matrix allows a faster convergence (in fewer iterations) but the alternatives might be more effective or more robust according to the situation.

As the method is iterative, the process should be stopped when the stop criterion is reached, in other words, when it is verified that a given value (or several values) becomes negligible. The global algorithm can be written as follows:


by defining the increments, i indexing the Newton iterations and ε being a positive value close to zero.

Note: The Newton algorithm is used to solve the equilibrium at each time step. It can also be used to find the stresses in each Gauss point (at all iterations of the Newton problem on the global scale) when the constitutive law requires it.


$translationBooks

C2. Why performing non-linear calculations

C2. Why performing non-linear calculations

As in many areas of physics, it is only reasonable to undertake a non-linear calculation after having a good idea of the "end of the story". In other words, how the structure will evolve until it becomes unstable. Below are presented the good and more questionable reasons for undergoing non-linear calculations.

  1. Good reasons

These procedures are reasonable, because:

  1. Questionable reasons

Those motivations often end up leading to non-convergence or “zero-pivot” problems. When the calculation is interrupted at a given stage of the iterative process, the deformation of the structure (given by the deformation of the mesh, or the strain field) is misleading:


$translationBooks

C3. Implementation

C3. Implementation

The correct method consists of always proceeding step by step so that non-linearities are not introduced all at once:

Every step should be verified properly, and it is necessary to make sure that the solution converges when the mesh and the time-step are refined. 

Finally, it should be clear that in non-linear cases, the problem to be solved may have no unique solution (instability, bifurcation, etc.).

  1. Choosing the discretization

The discretization must be chosen focusing on the fact that it must be consistent with the loadings that one wishes to apply, especially if those are not monotonic.

The more important the non-linearities are, the smaller the load-increments must be. Also, if the behavior of the model is complex, the integral of the constitutive law is not completely implicit: in that case, it should be verified that the increments are sufficiently small for the model to converge.

According to the standards, the management of the time-steps can be automated to a certain extend. This allows us to reduce the time-steps if there is a convergence problem, or even to increase it if the non-linearity is weak and the convergence is fast. For the sake of diligence, it is advised to verify that the convergence is reached by performing a second calculation with smaller load increments.

  1. Choosing the boundary conditions 

It should be clear that in the case of softening-types constitutive laws or if a limit load exists, the loading with imposed forces might become not admissible, as shown in the figure below.

In the first case, the closer the imposed forces get to the limit load, the harder it will be for the model to converge until it becomes impossible. In the second case, it will be impossible to exceed the effort corresponding to the first peak and obtain a solution after that.


Examples of response with limit load or loss of stiffness

There are other cases, where there is no unique solution for a given effort or displacement. It happens, for example, in the case of buckling of thin shells or when the problem has bifurcations, for instance, dissipative branch and elastic branch (see figure below).

Response (a) with multiple branches (b) non-monotonic loading and displacement

If the model requires imposed forces, the problem can be solved by continuation or arc length methods (see Code_Aster [R5.03.80] documentation for instance).

For an algorithm to work effectively, it is then necessary to conceive its model and to prepare its data, so that one looks for a unique solution by introducing multiple restrictions, such as:

  1. Evaluation of the convergence of the global non-linear algorithm

The stopping criterion for the Newton algorithm must be considered carefully. It is advised to verify what limit criterion is programmed by default in the software. For example, certain algorithms test the convergence using the norm of the residual (in that case, the criterion depends on the unit used for the forces in the calculation, the precision is not the same if the required residual is 1N or 1MN). Other software use a criterion where the norm of the residual relates to the norm of the loading (in this case, the convergence does not depend on the chosen physical unit). Special precautions must be taken when different elements are mixed (solid, structural elements) or when compound elements are used.

In general, one must be careful when the convergence criterion is modified from the recommended value to reach certain results. 


$translationBooks

C4. Convergence issues? Symptoms and solutions

C4. Convergence issues? Symptoms and solutions

As opposed to linear elastic calculations, non-linear calculations can perform iterations indefinitely, never verifying the stopping criterion, but they can also be interrupted abruptly at a given load increment (after the first one), with the warning “zero-pivot”.  

Let us examine the causes of this type of errors.

  1. Non-convergence of iterations

This phenomenon is almost always a direct consequence of the inadequacy of the used iterative method: the problem to be solved certainly does not satisfy the mathematic conditions, which ensure convergence (non-contracting operator, non-positive operator, non-convex operator, insufficient restrictions in the interval in which one looks for the solution, etc.).

This phenomenon may be a consequence of accumulated round off numerical errors that can be attributed to the algorithm itself; however, it does not happen often.

The evaluation of the residuals over the iterations may provide information about the nature of the problem:

After a given state, the residuals increase abruptly to a certain value and then start decreasing again, etc. This behavior appears when the algorithm finds multiple possible solutions: it converges to one of them, then “bounces” towards another one. In other words, there are multiple possible equilibrium configurations!

  1. Abrupt stop

The abrupt stop of the algorithm in a non-linear calculation is often reported as a “zero-pivot” warning when solving a linear system. This expression refers to the “pivot” in the Gaussian elimination method, which is the algorithm used (with small variations) to solve linear systems. Most parts of iterative algorithms replace the problem with a linear system.

This situation corresponds to the structure becoming unstable; in elastoplasticity or damage, it means that the structure lost its cohesion. For large displacements or contact problems, it reflects the formation of a mechanism, especially, the possibility of rigid body movement: all interface nodes slide, and the restraints are not enough to reach equilibrium.

It may also happen that the algorithm stops for reasons not connected to the state of the structure whatsoever: disk space overflow or programming bug. The first one is easy to solve.

  1. When nothing works

If a non-linear calculation does not converge, and one cannot find in the explanations above the reason, it is recommended to ask oneself the following questions:

 

$translationBooks

Chapter D. Civil Engineering

Chapter D. Civil Engineering

D.1 Civil engineering materials

D.1 Civil engineering materials

D.2 The categories of structural elements

D.2 The categories of structural elements

D.3 The steps of construction

D.3 The steps of construction


$translationBooks

D1. Civil engineering materials

D1. Civil engineering materials

The objective of this chapter is to highlight the peculiarities of civil engineering materials regarding FE modeling. Among those, there are the particularities of the civil engineering materials themselves (concrete, timber, steel), and their use within the structure. 

  1. Concrete

Concrete is a compound material composed of cement paste, and aggregates of various sizes. The cement paste that contains water and additives acts as a binding agent (glue) between all the aggregates. Even though concrete is a heterogeneous material, it revolutionized the world of construction because of its characteristics such as its versatility, adaptability (offering a variety of aspects: form, color, texture), mechanical resistivity, fire resistivity, acoustic insulation, and diversity (high-performance concrete, precast concrete, etc.). However, concrete is subjected to degradation processes. The latter can be due to environmental factors such as humidity, rain, cold, freeze-thaw, heat, wind, drought, aggressions from deicing rock salt, alkali-reactions, internal or external sulfate-reactions, corrosion, mechanical solicitations due to collisions, and load increase. Depending on the degradation phenomenon, the deterioration of the structural element is progressive and initiates as soon as the concrete tensile stresses are applied. Degradation leads to the emergence of small cracks and localized deformations. The increase of the deformation leads to the decrease of the material’s strength, which can then lead to softening behavior and the collapse of the structure. 

Constitutive law – In the FE calculations, generally, only some aspects of the concrete behavior are considered because of its complexity. To express the most relevant constitutive law, it is necessary to well understand the characteristics of concrete. Thus, it is essential to know the behavior of the cement paste, as well as the behavior of the aggregates, both responsible for the stiffness properties of the material. The paste is the cause of the mechanical strength of the concrete and, because of its water content, is also the cause of the deficiencies such as the increase of porosity, the decrease in mechanical strength, the delayed effects (shrinkage, creep), and the transfer of aggressive agents. 

Concrete is considered as a “quasi-fragile” material because of its complex behavior. Its main characteristics are highlighted in the literature. The most significant one is its fragile behavior in tension and more ductile in compression and when subjected to unilateral forces in a tension/compression cycle. Other properties must be considered such as the shrinkage and creep according to the environment and the loading applied to the concrete. 

In the literature, the two main approaches to model the behavior of concrete are:

However, in the literature, the two most used model categories are the ones based on mechanical damage and the ones considering explicit discontinuities. 

Although the elastic component is not exactly linear ([Baron] p.276; [Sargin 1971]), the calculations are generally conducted using the isotropic linear elastic assumption (Hooke’s law) to represent the elastic phase of the concrete behavior. The values of Young’s modulus and Poisson’s ratio are computed by characterization of the concrete. In the eventuality that other mechanical characterizations of the concrete cannot be conducted, some properties such as the compressive strength can be deducted using model laws that can be found in model-codes such as Eurocode 2 for instance. 

Note – Design offices use the model-codes except for study cases when they conduct detailed analyses on a specific structural element on which particular mechanical phenomena occur: cracks, heave, or local heterogeneity. Indeed, there exist finite elements of “cracked element” type, damage models, viscoplastic laws, or even poro-mechanical coupling laws. However, to conduct these studies, data from on-site or sample analyses of the material must be available to orient the engineer towards adequate constitutive laws, for which, coefficients and moduli will have to be determined anyway.  

Elastoplastic model – For the sake of simplicity, concrete is often modeled as an isotropic homogeneous elastoplastic material, which is acceptable until the concrete experiences cracking. 

Regarding an elastoplastic law, concrete cracking is not directly modeled: the cracking zones are characterized by important anelastic deformations (>1-2%) and a fixed constraint state between tensile ft and compressive strength fc of concrete. 

As shown in chapter 3, an elastoplastic constitutive law for concrete is composed of an elastic law and a plastic flow rule associated with a “criterion of plasticity”.

However, the post-cracking behavior can be approximated using a stress-strain curve going further than the tensile strength. The latter can cover the effects of tensile softening (representing the necessary work to open the crack) and tensile hardening (concrete contribution in between the cracks, stresses due to reinforcement adherence).

If the decreasing behavior of the material leads to a global decreasing behavior, one must pay attention to the eventual local effects: the size of the finite elements will limit the size of the anelastic zones and the solution will depend on the used mesh. Different numerical techniques enable us to solve or limit this issue. 

Damage – A damage law is a law that considers one of the main macroscopic effects of concrete cracking: the loss of material stiffness. The fundamental idea is to renounce following the eventual cracks (their apparition and propagation) and consider that the concrete of a given structure deteriorates by multiplication of cracks in the damaged areas. This type of law allows us to describe the decrease in material stiffness undergoing small cracks. This stiffness loss is measured by an internal variable called damage, denoted D, that evolves from 0 (undamaged material) to 1 (totally damaged material). This variable is generally a scalar value. 

To best represent the behavior of concrete, the damage laws consider post-peak softening behavior. This enables determining the stresses as follows:

with  

The advantage of this method is to consider the concrete as a “continuous medium”, for which the FE methods are perfectly adapted. 

In the case of concrete behavior, the two main damage modeling families are the anisotropic and isotropic ones. The isotropy characterizes the invariance of the physical properties of concrete regardless of the direction. On the other hand, the anisotropy depends on the direction. An anisotropic law presents different responses due to loading according to its orientation. 

One of the most used damage models in the industrial and research world is the Mazars model [Mazars, 1984]. It is certainly the first damage model for concrete that works accurately. 

The main difficulties raised by the damage models are:

Several regularization classes exist including the non-local regularization and the regularization by cracking energy [Hillerborg, 1976] (that solves the problem only partially). Among those non-local methods, it is possible to cite the integral methods [Pijaudier et al., 1987], [Giry et al., 2011] or gradient methods (deformation gradient or internal variable gradient [De Borst et al., 1992], [Peerlings et al., 1996], [Nedjar, 2005], [Lorentz, 2017]. These methods require the use of relatively thin meshes, which renders the analysis computationally expensive. 

Delayed effects of the stress relocation – When analyzing the behavior of a concrete structure several weeks after pouring, but especially for the long term, it is necessary to consider the delayed effects such as shrinkage and creep. 

These phenomena specific to concrete can, in principle, be modeled adopting a visco-elasto-plastic constitutive law (Bingham law), or sometimes only visco-elastic (“scientific” creep, [Eymard]): this approach is generally applied by research laboratories to analyze tests on materials. However, in the case of a refined model considering the delayed effects, it is necessary to incorporate phenomena such as drying and cracking in addition to shrinkage and creep. Indeed, all these phenomena interact with one another, so it is important to model them numerically. 

The deformation caused by the shrinkage is induced by the drying of the concrete due to the environmental effects. The shrinkage leads to a differential deformation, meaning that the stresses are more significant where the concrete dries faster. This phenomenon causes tensile stresses at the surface, which then leads to cracks and compressive stresses inside the body. 

Concerning creep, its deformation is generally separated into two deformations, one is due to creep itself, and the other one to shrinkage from the drying process. The latter can be explained by the strong dependence between relative humidity and creep. 

In Eurocode 2, it is possible to determine concrete deformations due to the delayed effects (without external loading). For this purpose, one must compute the deformation due to the endogenous shrinkage (caused by the internal humidity) and the desiccation deformation (caused by the drying process and the size of the structural element). According to the §3.1.9 of the EN1992-1-1, the creep deformation under compressive stresses σc can be written as:

with Ec the tangent modulus (equal to 1.05 Ecm) and φ the creep coefficient. 

The Eurocode 2 (EN 1992-1-1) describes a simplified calculation method of φ (∞, t0). Annex B gives a more complete method enabling the estimation of φ(t, t0) as well as the shrinkage evolution (see the annex B of EN 1992-2). Moreover, it is important to remember that the calculation of the delayed deformation is relative to the type of cement paste. 

  1. Structural steel

The behavior of steel is much simpler than the behavior of concrete for several reasons: it is an isotropic material with identical strength and moduli in tension and compression. Moreover, it is subjected to industrial quality control to ensure its homogeneity. 

Elastic models – Even though steel behaves essentially as anisotropic elastoplastic material (Eurocode 3, part 1.1, §5.4.3), the typical models for steel structures, or compound containing steel, adopt a linear elastic behavior. One must then check that the elastic constraint threshold fy was not reached. In the case of beam or column type elements, the codes allow, if the sections have a size that guarantees a sufficient local ductility, to exceed in the analysis the elastic strength and consider reaction moments based on the plastic distribution of the stresses. 

Elastoplasticity and strain hardening – The elastoplastic theory was developed from the study of steel alloys, mainly to predict the rolling and forging stresses ([Hill], [Nadai]). For steel, the usual model is composed of Hooke’s law for the elastic deformation and the plasticity criterion denoted “Von Mises criterion”, as explained in the Eurocode 3 (part 2, §7.3; part 1.5 §10 and part 1.7 §5.2.3.2): for this criterion, only the elastic constraint threshold fy of steel must be provided. 

Considering the FE calculations, another issue is the strain hardening of the material, meaning the hardening in the purely plastic phase. This aspect is notably explained in the annex C6, part 1.5 of the Eurocode 3. Since it is hard to verify the model, it is preferable to conduct a preliminary computational model with no strain hardening: indeed, one can then verify the quality of the obtained results in the plastic zones by observing the eigenvalues of the Von Mises stresses (the plastic zones must be about monochrome). Nevertheless, it is important to note that the use of an elastic perfectly plastic law with no strain hardening can lead to convergence issues of the linear analysis. Indeed, the plastic zones have zero tangent modulus and, therefore, no stiffness. Considering strain hardening can stabilize the numerical resolution. 

In the case of strain hardening, until a stress fy+X is reached, the possible loading/unloading phases are elastic, with a modulus equal to the initial one. Thus, fy+X becomes the new elastic yield stress of the material. Moreover, the behavior of the steel is close to the kinematic strain hardening model. As a first approach, one can consider a tensile strain hardening up to a certain stress value (fy+X) that will decrease the compressive elastic yield strength to fy-X, and vice versa. 

The most common steel model is the isotropic strain hardening, which depends on the accumulated plastic deformation: 

where  and fy an increasing function.

In this case, there is no distinction between tension and compression, so a strain hardening in one direction leads to an increase of the compressive elastic yield stress. This hypothesis is accurate only if the considered loading is monotone and not cyclic. 

In most numerical computational models, one can retrieve the eigenvalues of the cumulated plastic deformation. Thus, it is possible to verify that the Von Mises constraint is distributed according to this deformation in the plastic zones. 

In the specific case of stainless steel, the hypothesis of a linear behavior up the limit of elasticity fy corresponding to plastic deformation of 0.2% is not respected, so it is recommended by the Eurocode 3-1-4, § 4.2. to consider the effects of the nonlinear behavior on the force-displacement relationship in the deflection calculations. The material law is of Ramberg-Osgood type and recommendations to implement this model can be found in annex C of this section of the Eurocode. 

  1. Pre-stressing steel

Regarding their behavior, steel cables (braces, tension cables) and steel for pre-stressing differ from other structural steel only by an elastic limit about three times greater (from 1680 to 2140MPa in France), which is justified by the serviceability solicitations, generally greater than 1000MPa. One can then, firstly, model the pre-stressing reinforcement as linear elastic elements or elastoplastic of Von Mises type. Steel elements for pre-stressing are all hot-rolled, so the linear isotropic strain hardening law in the plastic phase should be taken into account. The strain hardening coefficient h can be calculated using the characteristics guaranteed by the manufacturer:

with Rm the constraint to failure, fp0,1 the usual elastic limit, E the Young modulus, and Agt the elongation at the point of failure. 

A direct consequence of this high level of stress is the initiation of the relaxation mechanism. Relaxation of steel is a non-elastic delayed effect that depends essentially on the time spent since the loading was applied. It leads to a progressive decrease of the constraint for bars and cables subjected to a constant elongation. This mechanism initiates at ambient temperature only for elongations exceeding about 60% of the elastic limit (≈1000MPa). Relaxation increases slightly with temperature. 

Alike concrete creep, it is possible to model the relaxation using a linear visco-elastic constitutive law. However, this approach is mostly used in research work. It is not common to model steel for pre-stressing: generally, it is accounted for by introducing a distributed force within the pre-stressed concrete. Nevertheless, for detailed analysis, in which the interaction between the steel cables and the surrounding medium is of interest, relaxation can be considered using an incremental time-dependent calculation with the relaxation losses introduced as initial stress conditions. 

  1. Passive steel

Normal mechanical actions of passive steels – For FE calculations, passive steel elements are frequently modeled by linear 1D elements of “bar” type. According to the historical reinforced concrete methods, the reinforcement is considered to work only in tension/compression along its own axis, which is exactly what FE models enable. Connecting those bar elements to the nodes in the solid elements representing the concrete makes it possible to avoid mechanisms specific to “bar-only” assemblies.

In most computational software, the bar elements are considered by default as elements with a linear elastic behavior. Considering the 1D characteristic of these elements, the law is of the type N = E A u/L, with u the axial displacement, and N the normal stress at the nodes. It must be verified afterward that the axial force N/A stays within ±fy.

Friction – What was said in the above paragraphs is valid only if there is a continuous friction force acting between the steel reinforcement and the concrete. 

  1. Timber

Timber is a material that possesses some inherent peculiarities: it is not homogeneous (this is true at various scales), it is not isotropic, it does not show symmetrical behaviors, and it can be subjected to ductile or fragile failures according to the solicitation and its orientation. Timber is sensitive to humidity, which has a direct impact on its dimensional characteristics, stiffness, and strength. The duration of the loading has an important impact on the strength and the deformation of the wooden elements. The variations in humidity can amplify the deformations (mechano-sorption). Relying on the previous observation, one can imagine a method to implement the FE analysis in this particular case. 

Homogenization – The hypothesis must be validated for the Representative Elementary Volume (REV), or at least, the volume of the smallest finite element. Knowing that according to the different timber species the growth ring can exceed 1cm, it becomes difficult to assume this homogenization hypothesis close to the connections (peaks, pins, bolts…) with a diameter of the same order of magnitude, or smaller. 

The presence of nodes is rarely considered when modeling structures or structural elements. 

Orthotropy –Timber has a structure and characteristics that are a function of three directions, the lengthwise directions (the axis of the trunk, wood fibers), the radial directions, and the tangential directions perpendicular to the lengthwise direction. The last two constitute a plan, often named the cross-sectional plane of the beams, on which the growth rings can be visible. The orthotropy coordinate system is therefore quasi-cylindrical. On the contrary, the coordinate system in which the elements are inscribed is cartesian. The representation of this orthotropy, when considered, is in most cases limited by an isotropic transverse hypothesis (radial and tangential axes having identical characteristics) in a cartesian frame. The slope of fibers associated with the presence of nodes is generally not considered in the FE calculations. Rather, it is modeled in the assembly calculation especially for the assemblies by contact.  

The flexibility matrix (inverse of the stiffness matrix) can be defined as follows for an orthotropic plan:

Note: the orthotropy and symmetry hypotheses of Sij reduce the number of independent terms from 36 (most general 3D case) to 9 terms. 

Timber – orthotropy hypothesis

Elastic modeling – One must simply know the isotropic transverse matrix of behavior for a 2D or 3D model. The moduli between the lengthwise/longitudinal, radial, and tangential directions can present ratios of the order of 20. They evolve as a function of the loading period (creep), timber humidity when constructing, and during its lifetime (environment). The representation of the FE method will depend on the pertinence of the parameters considered. 

Plastic modeling, failure criterion – For an isotropic transverse non-symmetric material, the Hill, Tsai… criterion can be used knowing that it will be necessary to describe the perpendicular tensile and shear fragile failures. The significant strength variability makes it harder to fix the criteria parameters. The sequence of growth rings or assembled layers with different mechanical characteristics can render a homogeneous model difficult to fix in terms of strength.  Indeed, the strength and stiffness are highly correlated for timbers and “systemic” effects appear quickly in terms of element strength. Consequently, there are strength limits in axial tension, axial compression, bending….

Finally, timber structures are particularly sensitive to the behavior of their assembly. Thus, the elements do not need to be modeled in detail, only the assembly zones do. Nevertheless, problems can arise because of the contact zones, the multitude of materials, some of them entering plasticity, the fragile failure of others, and the homogenization limits reminded earlier. 

It is clear that the modeling effort is related to the scale of the investigation or the advancement of the project. Provided the above singularities are considered, the modeling of timber inside of a structure can be conducted like any other material. 

Modeling of the delayed effects and the interaction with the hydric/water phenomena – Timber is a hygroscopic material sensitive to the relative variation of humidity of the air. Moreover, it is subjected to creep under loading. If one wants to model these phenomena, a rheological viscoelastic model in a variable environment and respecting the thermodynamic principles can be adopted. The generalized viscoelastic behavior of Kelvin-Voigt type, characterized by rheologic parameters and dependent on the level and history of humidity can be associated with the Ranta-Maunus non-viscoelastic behavior to characterize the shrinkage-swelling and the mechano-sorption. 


$translationBooks

D2. Different categories of structural elements 

D2. Different categories of structural elements 

Here are presented the peculiarities of these elements regarding FE calculations. 

  1. Reinforced concrete elements

Considering the delayed phenomena – Usually, the project engineer is interested in the redistribution of stresses and the delayed deformations that accompany the aging of the materials. 

To achieve this, the static calculation by time increments is particularly well adapted and largely sufficient. The shrinkage and the creep of concrete for instance depend only on the time that has elapsed since the pouring phase. Using the incremental calculation, they will be introduced as imposed deformations at each node of the mesh. It is possible to compute beforehand, for each time increment the shrinkage and creep related deformation map: the FE software incorporates the initial state seeking an equilibrium of the elastoplastic material. 

Nevertheless, attention should be paid to the interaction between the delayed effects and the construction phases, as it will be explained later.

Consideration of cracking - When concrete cracks in tension, cracks develop towards the nearest face, as well as along the steel-concrete interface ([Goto], figure below).







Internal cracking of the reinforced concrete

When considering a tensioned connecting rod of cracked reinforced concrete as a homogenized continuous medium, the behavior law N(ε) (with N the normal tensile stress) follows the trend of the figure below.









Model of the tension stiffening effect

The Model Code 1990 of the CEB-FIP expresses an analytical formulation. This relationship can be used by modeling the section of cracked concrete as a bar element whose stiffness would be updated in increments (depending on the value of normal stress or elongation, since the relationship is invertible). It can also be used as a law for the behavior of a fiber of a multi-fiber element.

Cast Concrete - For cast concrete areas (e.g. by passive reinforcement), it is possible to take into account the residual post-crush strength, as indicated in EN1992-1, using Sargin's law:

with εc1 the peak value such that

with

For

The EN1992-1 §3.1.9 also suggests an increment of strength and deformation of concrete when it is subjected to a confinement stress of σ2=σ3 such that: 

Une image contenant texte Description générée automatiquement

However, this law is only one-dimensional and is only valid under monotonic loading. Therefore, it can only be used with the finite element method when modeling concrete confined by a bar element (so the modulus would be incrementally modified), or as a constitutive law of a fiber of a multi-fiber beam element.

  1. Pre-stressed Concrete Elements in Pre-tensioning

Prestressing by pre-tensioning, which is characteristic of industrially prefabricated concrete products (beams, pre-slabs, honeycomb slabs, girder, etc.), consists of tensioning cables (wires, strands or bars) on manufacturing benches, then pouring the concrete before slackening the cables when a minimum concrete strength is reached (called concrete sag resistance). The amount tension in the cable (which must not exceed the maximum prestressing force allowed by design codes), the number of cables, and the strength of the concrete are adjusted according to the loads that the prestressed floor or element must resist.

When the cables are loosened, instantaneous losses, of the order of 8% for a prestressed pre-slab for example, must be taken into account (losses due to the entry of anchors, to the relaxation of the prestressing reinforcement during the period between the tensioning of the cables and the transfer of the constraint, to the elastic shortening of the concrete under the compressive stress imposed by the constraint) when designing in the provisional phase of construction with a prestressing transmission length to be considered starting from the end of the precast concrete element.

In the longer term, delayed losses due to shrinkage of the concrete, relaxation of the steel, or creep of the concrete, which will finally reach a total loss of 20% for a prestressed pre-slab, must also be considered.

In many FE computational software, prestressing can be incorporated in beam-type finite elements representing the cables connected to volumetric finite elements representing the concrete. Depending on the study area, it may be necessary to consider the transmission length of the prestresses in the cables, for example in the end zones. The distribution of the actual prestressing force (considering the instantaneous or delayed losses as a function of the moment in the service life of the product where the calculation/verification is to be performed) is then variable along this transfer length. A linear distribution of this prestressing force is allowed in most computational software and remains safe for dimensioning purposes compared to the more realistic parabolic distribution.

For reasons of complexity and need (designs limited to elastic deformations), the modeling of these elements is generally carried out with linear assumptions (constitutive laws of concrete and linear elastic steel, perfect contact between concrete and steel...). However, for refined studies, non-linear hypotheses can be followed according to the need such as damage type constitutive laws for concrete, elasto-plastic type constitutive laws for steel, and the introduction of steel-concrete interface elements.

  1. Pre-stressed Concrete Elements in Post-tensioning

Elastic losses - Post-tensioning of the prestressing cables is accompanied by instantaneous losses: friction, anchor recoil, and loss due to elastic elongation.

Delayed effects: creep, shrinkage, and relaxation - Delayed effects are considered using an incremental calculation. Creep and shrinkage can be introduced as incremental volume deformations given at each node of the mesh.

  1. Steel Structures

Choosing the type of analysis

Due to their high slenderness, steel structures are very deformable. As a result, the traditional assumption of reaching stress equilibrium in the initial configuration is not always valid. It is then necessary to establish the internal stress distribution in the deformed configuration.

The sensitivity to these non-linear effects, sometimes called second-order effects, is expressed through the critical multiplier αcr: the load multiplier leading to the Eulerian instability of the structure. In the current version of Eurocode 3:

In the last two cases, if the global geometrical irregularities (or imprecisions) and the element irregularities influence the global behavior, they must be considered.

Element irregularities include:

The latter can be represented by an equivalent geometric irregularity, whose value can be found in the current standards.

These irregularities must be included in any analysis of a structural member incorporating non-linear effects.

In principle, structural or Strength of materials models (bars, beams, plates, and shells) are ideally suited for structural steel calculations.

Wireframe structures

Elastic analysis

For modeling purposes, it is necessary to analyze precisely:

Plastic analysis

When the elementary ductility and/or the assembly ductility is ensured, it is possible to carry out structural analyses incorporating plasticity. Different methods can be used (EC3, §5.4.3):

Torsion

Steel elements are rarely massive, they are composed of thin layers to form profiles, open or closed. Especially in the first case, the response to torsional stresses is both in uniform torsion, known as Saint Venant torsion, and in non-uniform torsion, leading to the buckling of the section. The latter phenomenon is usually not considered in commercial software, even though it can have a significant effect on the response of structures. In this case, two solutions are possible:

Decomposition of torsion into uniform and non-uniform parts: simplified biflexion hypothesis

2D or 3D elements 

The linear elastic analysis of two- or three-dimensional steel elements does not present any specific problems, so the general rules apply.

On the other hand, a non-linear analysis may be required to study instability phenomena. Indeed, structures composed of steel plates, whether flat or curved, are subject to buckling phenomena.

In the case of flat elements, called plates, buckling is a relatively stable phenomenon: the initiation of buckling of the plate does not lead to failure, the maximum load is reached after buckling. This is called post-critical behavior.

In the case of curved elements (shells), instability leads to the immediate and often brutal failure of the structure. From a numerical point of view, in the non-linear static analysis, this leads to a rapid decrease of the load beyond the maximum.

In both cases, the maximum load is strongly dependent on the initial deformation applied, in amplitude and shape. The amplitude is fixed by standards. The shape is usually chosen affine to the first mode of instability. However, this choice is not necessarily the most detrimental. It is advised to supplement it with local modes when the structure presents panels of strongly different dimensions.

For example, in the case of an orthotropic deck, the deformation affine to the global buckling mode must be supplemented by deformations affine to the buckling of the under panels.

In the case of shell-type structures, the problem is even more critical. It is advised, once a first calculation has been performed according to the above assumptions, to adopt in a second calculation a shape of initial deformation affine to the deformation obtained at failure.

  1. Compound structures: steel-concrete

The alternative of a compound steel-concrete construction is sometimes preferred for certain types of industrial buildings and bridges with small to medium spans (midspan < 100 m). The combination of these two materials by making them "work" in their strength domains (concrete in compression and steel in tension) can result in a strong and lightweight design. To achieve this result, the connection between these two materials must be correctly designed. A distinction is made between:

Overall structural analyses are usually carried out in the elastic domain by homogenizing the section, or by representing the two materials separately. This second way of proceeding can lead to difficulties in the processing of the results. Indeed, it requires recalculating the stresses on the globalized section to apply the current standards. 

When more detailed analyses are needed, given the diversity of materials involved in this type of combination both geometrically and in terms of non-linear behaviors, 3D finite element models are necessary to conduct local analyses, which include the processing of the various interfaces (using contact finite elements for example). For studies at the scale of the structure or the elements, relatively high-performance 2D models have been developed over the last decade, such as those based on a fiber cross-section cut-out (fiber model) to enable the section stiffness to be estimated by numerical integration.

The assembly of composite bridge girders is another relatively complex detail. Regardless of the assembly solution chosen (butt-strap joint, joist trimmer, or diaphragm), 3D models using solid finite elements are preferred to simplified 2D models, but they are more computationally expensive.

Cracking - Composite girders are usually made of a steel section that is connected to a reinforced concrete floor or deck by means of connectors. As a result, in areas undergoing a positive moment, where the slab is compressed, strength and stiffness are particularly consequential, whereas in areas undergoing a negative moment, cracking results in significantly lower mechanical characteristics. This cannot be neglected in the overall structural analysis. Different levels of modeling are allowed in the codes:

Connection

Except in the case of non-linear analysis, it is not useful to model the connection. Current standards allow us to consider the effect of a partial connection on the strength of the elements.

For modeling purposes, beam-type finite elements (2D or 3D) are generally used to model the point connection. If a distributed connection hypothesis is considered, appropriate models are available in the literature.

Collaborative widths

Steel beams are connected to particularly wide members, and shear dragging can result in a non-uniform stress distribution across the width of the slab.

In wireframe models, this phenomenon is usually taken into account by modeling a slab of reduced width, at constant stress, since it is required for beam models.

Strictly speaking, since shear dragging is related to the transmission of shear forces through the connection to the slab, it is therefore dependent on the shape of the moment diagram. The collaborating width should vary from one combination to another. However, the standards allow a single width to be adopted for all calculations.

The variability of collaborative widths will be considered if the slab is modeled using surface type elements.

Delayed effects and shrinkage

Since the delayed effects of concrete influence the stiffnesses and stress distribution throughout the structure, they must be taken into account. It should be noted that the Eurocode adjusts the value of the steel-concrete equivalence coefficient according to the type of loading. The shrinkage solicitation induces a self-balanced stress distribution throughout the section; thus it should also be taken into account.

Non-linear analyses

Non-linear computational models of composite structures assume the material and geometric irregularities hypotheses used when modeling concrete and steel materials separately.

As mentioned earlier, the connection must be modeled with its own stiffness. However, it should be noted that its numerical processing remains difficult. Usually, it is assumed that the concrete and steel parts cannot be separated transversally and that only longitudinal sliding is possible. The finite element formulation of such a connection can lead to the "locking" phenomenon when the sliding is blocked by the looseness of the transverse stiffness projection preventing the uplift when the equilibrium is reached in the deformed configuration. Therefore, when this type of element is used, it is advised to verify the consistency of the reaction forces in the connection with the tensile/compressive forces of the steel elements.

  1. Braces, tension cables, and suspension cables

Introduction to the calculations - Some cable elements can be modeled as bars of the equivalent section. This is the case for vertical tension cables, or prestressing cables, guided in their duct. The curvature of these cables at equilibrium is practically independent of their linear mass. In other cases, the cables tend to be bent by their own weight: the vertical direction is a specific direction that plays a role in the stiffness of these elements.

When, at the scale of the cable's span, its curvature becomes important, it strongly influences the stresses and strains it can transmit to the rest of the structure. As an example, for braces, the stiffness of the cable depends mainly on its deflection and elongation. Deflection affects the stiffness of the cable; thus, the tension/deflection relationship of the cable is inherently non-linear.

The cable elements present in most codes are based on the [Gimsing] model, which assumes that the cables deform as a parabolic curve, an assumption that is valid if the deflection/span ratio is less than 1/12. This model decouples the elastic elongation from the bending of the rope.

The catenary effect - When the deflection of a rope is greater than 1/12th of the span, it is no longer possible to decouple the elastic elongation from the bending because the tension lever arm becomes the dominant factor in the expression of the bending moment. Several computational models suggest catenary elements to account for this geometrical effect. These elements connect two nodes of the mesh so that a single element is sufficient to model the wire rope.

Nodal displacements of a catenary element

For this type of element, the nodal unknowns are the vertical displacement (in the direction of gravity) and the horizontal displacement (see figure 4). The associated reactions are the variation of the horizontal component of the cable tension, and the vertical reaction at the anchors. Since the relationship between nodal displacements and these reactions depend on the cable tension, finding the equilibrium becomes a non-linear problem, although the structure is globally elastic.


$translationBooks

D3. The different construction phases 

D3. The different construction phases 

Studying the phases of construction has two objectives:

There are various effects related to assembling. They are due to:

Taking these effects into account can be relatively complex and, in the most difficult cases, it may be essential to use software that can model the evolution of the structure step by step.

However, it is often possible to proceed by superimposing various linear analyses.

Effect of the construction phases in the case of the keying of a cantilevered bridge span

The main issues are related to the effects of the delayed deformation of concrete. Indeed, how can one evaluate the final state midspan moment in the above example? In the case of a single keying point, it is possible to use the method known as the "coefficients method" (Figure 6). This approach is based on the following arguments:

Final state = (E(t0, t1)/E(t0, t∞)) x Final state not keyed + (1-(E(t0, t1)/E(t0, t∞))) x Final state without phasing.

with t0 the duration of load application, t1 the duration of the keying, t∞ the time considered for the final state, and E(t0, t1) the concrete modulus for obtaining the deformation of the concrete at time t1 for a load applied at t0.

This method, in the case of a single keying point, outputs the exact theoretical final state. However, it is difficult to extend it to the case of multiple keying points and can lead to absurd results.

The unitary case used for the coefficients method

It is best to externalize the effects of changes in the static state as follows (Figure 7) :

This technique can be extended more easily to cases where the static configuration of the structures is modified many times.

Unitary cases used for the superposition method

It is worth remembering that, under the effects of delayed deformation, concrete reacts with an apparent deformation modulus called relaxation modulus, which is lower than the corresponding creep modulus. If the classical ratio between the moduli of steel and concrete, including creep effects, is about 18, in the case of an imposed displacement, this ratio increases up to 24. This tends to make the adjustments by support elevation and actuators less efficient when the latter leads to imposing a deformation on the structure.

From the point of view of FE simulations, the large number of intermediate states that need to be processed multiplies the risk of errors. The verifications must concern:

In construction studies, one must also remember that the actual creep of concrete can deviate strongly from the theoretical formulations. The model must therefore be constructed such that it is easily adaptable to restore the deformations appearing in the first phases, and thus improve the prediction of the following ones.

 

$translationBooks

Chapter E. Typical post-treatment in Civil Engineering

Chapter E. Typical post-treatment in Civil Engineering

Introduction

In civil engineering there is a great variety of quantities of interest because of the numerous subjects that can be studied:

Moreover, the quantities of interest may be directly or indirectly accessible when performing  FE calculations. They can be obtained from primary results later processed to provide a given quantity of interest. The latter is then compared to some criterion or used as data for the next steps (reinforcement area, rebar orientation, …).

The Eurocodes enables us to know the quantities that should be analyzed according to the various situations mentioned above. However, the way to obtain these values is not established in the codes and standards. This chapter highlights the maneuvers to access those quantities and the mistakes to avoid.

E.1 Generalities

E.1 Generalities

E.2 The quantities of interest in dynamics

E.2 Quantities of interest in dynamics

E.3 The particularities of reinforced concrete

E.3 Particularities of the reinforced concrete


$translationBooks

E1. Generalities

E1. Generalities

The different post-treated quantities may change significantly according to the objective of the finite element analysis.

Usually the studied quantities in the pre-dimensioning phase (a sketch or preliminary draft) are different from those searched in the pre-design and design phases. Other quantities of interest are needed when analyzing existing structures to quantify structural risks.

 It is tempting to consider:

It is important to keep in mind that not all quantities of interest are accessible for every type of element and constitutive models. The type of element determines the nature of the degrees of freedom: the nodes and moments on solid elements are not accessible. Analogously, a model simulating plasticity cannot provide damage results. These details may seem trivial, but they are often forgotten by the users, so it is important to mention it. 

  1. Stresses and strains

As exposed in the “Generalities” chapter, the stress and strain fields are not obtained directly when solving a mechanical problem with a FE solver. They are obtained from the displacements by computing the derivative of the interpolated displacement fields. Using the same notation as in the “Generalities” chapter:

In linear-elasticity, the strain-stress relation can be written as (H being Hooke’s matrix): 

In this kind of procedure, the strain and stress fields evolve according to the interpolation functions of degree inferior to one of the displacement functions.

For non-linear analyses, the stresses and strains cannot be obtained by simply multiplying Hooke’s matrix by the derivative of the nodal displacements. They are based on the values of the element at the Gauss points, and the shape functions are used to extrapolate the results to the rest of the element and nodes. This procedure allows the stresses to evolve according to the same degree as the interpolation functions and not their derivatives. 

As an example, on an element with n nodes, npg Gauss points located at the coordinates ξnpg with shape functions Ni, we compute the least-square minimization between the interpolated field evaluated from the nodal values searched and the known Gaussian values. 

Building elementary nodal values starting from the values at the Gauss points in 1D

So, either it is a matter of minimizing the function:

Or for each node i among the n nodes of the element:

Which can be expressed in matrix form. The matrices are finally computed for all the isoparametric elements of reference:

And so, it is possible to directly obtain the nodal stresses:

Thus, it is possible to have several nodal values of strains or stresses for a given node, shared by multiple elements. The other nodal values can be deduced from those values.

Discontinuity of the elementary nodal values in 1D

Several methods can be used to recover the continuity of strains or stresses (if the material is the same in the 2 adjacent elements) to have a single nodal value of strain (or stress). Those methods consist of smoothening the stress and strain fields by using certain algorithms on the totality of the structure (by least-square minimization). It can also be done by computing the average of the values on a given node shared by multiple elements as shown in the figure below. 

Example of the calculation of an averaged nodal value starting from elementary nodal values in 1D

For convenience one might prefer obtaining smooth elementary values. However, when we look for precision, it is convenient to give priority to the values obtained at the Gauss points.

  1. Going from continuum mechanics to structural mechanics stresses

The automatic methods to compute reinforcement rely on algorithms that take into account combined axial forces and bending moments based on plate and shell torsors with 8 components, or deviated combined axial forces and bending moments based on beam torsors with 6 components.  

Some structures may need, for various reasons, to be modeled using solid finite elements (dam, structures such as prestressed tanks: liquefied natural gas tanks, retaining structures). To avoid computations based on 8-component torsors, it is necessary to rebuild the Structural Mechanics internal forces starting from the stress fields along the segment as shown in the figure below. 

Segment along which the stresses are used as a reference to rebuilding a shell type torsor

We consider the stress tensor expressed along the segment shown in the previous figure:

To be accurate, for thin plates satisfying the Kirchhoff-Love hypothesis, we consider the following torsor with 10 components that will be associated with the stresses according to the following equivalence principle:

This continuous integration is discrete along the segment and can typically be computed from the stress values expressed at the Gauss points or the nodes.

Multiple integration techniques can be used (Trapezoidal rule, Gauss, Lobatto, etc …)

A continuous expression will be kept for the sake of simplicity, even if it is not rigorous in the discrete context of the values obtained by the FE method.  

Taking into consideration the Cauchy theorem equivalent to the shear stresses at the adjacent faces, this torsor is reduced to 8 components:

Considering that this plate is an isotropic slab implies that the contribution of membrane efforts is neglected.

Then, it becomes a 5 component torsor (considering that the notation M_x represents the bending moment activating the rebars in the direction e_x), which are the internal mechanical actions in a plate element:


(1) Shear efforts outside of the plane (z-direction) with x as the normal direction

(2) Shear efforts outside of the plane (z-direction) with y as the normal direction


(3) Bending moment activating the rebars in the x-direction


(4) Bending moment activating the rebars in the y-direction


(5) Torsion moment on the cross-section of the slab with x or y as the normal direction


  1. Section method: elementary reduction elements (EF) to beam-type structural reduction elements

Depending on the stress state it might be convenient to build a model of the whole structure with plate elements including the zones where these are not appropriate such as posts and lintels.

In these zones, one must restore a beam-type torsor to take into consideration more accurately the local behavior, in particular, to estimate the required reinforcement.

A widely used technique to solve this problem consists of “cutting” the structure fictively and estimating the efforts at the cross-section. 

This “cut” must be chosen wisely to respect the Euler-Bernoulli hypothesis stating that cross-sections must remain plane. It is usually the case for posts and lintels.

Example of cross-sections in a lintel or post passing or not through nodes

We then build in the local coordinate system of the cross-section the beam-type reduction elements substituting the shell-type reduction elements. The latter can be the efforts at the nodes or evaluated at any point of the element coinciding with the section line and analyzed utilizing elementary form functions.

For the example illustrated above, it can be written for the cutting line, where local and global coordinate systems coincide, that: 


(1) Axial forces at the cross-section of the associated beam, i.e. the “cut plane”

(2) Shear forces on the plane

(3) Shear forces off the plane


(4) Bending moment on the plane


(5) Bending moment off the plane

(6) Torsion moment 


$translationBooks

E2. Quantities in Dynamics

E2. Quantities in Dynamics

  1. Time-dependent analyses

Post-processing of time-dependent quantities does not present any major difficulty as long as one is not interested in targeting a precise instant or in characterizing a single representative value at all times of the analysis.

For cases where the loading is controlled and not random, a single analysis may be enough. If it is preferred to extract a scalar quantity representative at all times of the analysis, the choice of the standard is the responsibility of the engineer who analyses the results, and it must be justified. A statistical characterization may also be relevant (typically expressed as a percentage).

If the engineer is interested in extracting a set of quantities constituting a vector, the problem of the concomitance of the normalized quantities (absolute, maximum, minimum values, etc.) arises. The time t at which a bending moment is maximal in a point for a torsor does not necessarily coincide with the point of maximal shear stress.

When the loading has a random character, it is important to multiply the number of calculation cases by integrating the random character of the loading. As an example, if we consider N cases of calculations and we construct for a quantity of interest gi(t) the corresponding dimensioning value Gi for a case i (it can be the maximum positive value and also the maximum negative value to then incorporate them into an approach of actions concomitance). We consider:

With λ(N) calculated from the Student variable for N samples (calculation cases) for a 95% confidence level (one-tailed) such that:

The table below provides estimations of the mean for various calculation cases.

  1. Seismic spectral analysis

The seismic spectral analyses provide representative quantities of the average of an extremum or period.

The purpose of this paragraph is to draw the reader's attention towards important precautions that must be taken to avoid making significant errors in the calculation of the results coming from a quadratic accumulation.

As mentioned in the dynamic chapter, it is possible to quadratically cumulate simply or completely the spectral responses of each mode. The result of these cumulations is a positive value.

It is important to consider that any operation that one wishes to perform on a quantity of interest that is the subject of a quadratic cumulation (simple or complete) of the spectral responses of each mode must be performed before the cumulation.

Let us take a simple illustration with the difference between the displacement of two nodes A and B cumulated simply on modes 1 and 2 (an illustration with a CQC combination would just be more troublesome to write using cross-terms):

Strictly speaking, the difference evoked as an example can be written:

Therefore, it is obvious that estimating the difference afterward is largely erroneous:

The first value is always positive, but it integrates algebraic differences of quantities in the same mode. The second expression might lead to a drastic underestimation of an extremum response. It can also be demonstrated by using the Cauchy-Schwartz inequality that the first estimate is the upper bound of the absolute value of the erroneous calculation presented in the second formula.

A more complex illustration can also be used with the estimation of a von Mises stress based on the main stresses for 2 modes :


An ‘a posteriori’ estimation of the accumulation of the von Mises stresses is largely erroneous:

Example of a specific application: estimation of the opening of a joint between two buildings

In seismic engineering, it is necessary to estimate the opening of a joint between two buildings under seismic conditions to ensure that there is no risk of collision. The accepted practice is to independently calculate the extreme values of displacement of the envelope of each building based on a quadratic accumulation and then to evaluate the maximum opening of the joint by calculating the difference between these two positive values. Thus, an unfavorable phase opposition of the maximum responses considered.


$translationBooks

E3. The specific case of reinforced concrete

E3. The specific case of reinforced concrete

  1. Overview of the common methods for designing the reinforcement of plate elements

There are mainly 3 methodologies applicable to plates that can be consulted using the following references:

  1. Example of the Capra-Maury method

A set of facets is defined, centered at the point of computation of the finite element code. This can be a node, a Gauss point, or any point where forces are interpolated.

Its normal vector rotates from this point in the plane tangent to the mean sheet. The facet is marked by the angle θ between the OX axis and the normal vector of the element in its local coordinate system (see figure 2.1-a). The angle θ is discretized regularly from -90º to +90º (here with a step of 5º). The Ox and Oy axes are the axes of the reinforcement layers.

Reference facet parallel to the beam’s cross-section in equilibrium subjected to combined bending and axial force

For each of these facets, the bending moment (M), the membrane force (N), and the shear force (V) are applied as a function of the effort tensors using the following equations:

By using combined bending calculations, it is possible to determine the reinforcement area in upper and lower layers AS(θ) and AI(θ) required in the θ direction to reach equilibrium in the chosen reinforcement concrete section.

The compressive strengths in the θ direction of the two reinforcement layers can be evaluated using the following expressions:

where fyd is the maximum allowable stress of the steel (identical in both directions).
The strength is ensured if the resistive force is greater than the applied force, which is expressed as:

Thus, by considering an orthonormal coordinate system with AXS on the abscissa and AYS on the ordinate, one can finally solve for the upper and lower reinforcement:

  for all angles θ

  for all angles θ

and

and minimum.

The resistance inequalities define for each θ a half-space limited by a straight line with negative slopes that reflects a validity domain, as shown in the following figure.

Resistance domain for a facet θ

By browsing through all the values of θ, one can obtain the validity domain shown in the following figure, delimited by the line ABCD.

Resistance domain for all facets

For each point P in the validity domain, the total area of reinforcement can be obtained by projecting the point P in Q onto the first bisector. The distance OQ then represents the value with AS = AXS + AYS.

Therefore, it can be noticed that the optimum reinforcement corresponds to one of the 36 points (considering the facet rotation step chosen if one facet is taken every 10 degrees for example) of the boundary (illustrated by the 4 points of the line ABCD ...) whose projection on the first bisector is closest to the origin of the axes. Researching this point can be done using "dichotomy" type methods.

  1. Connecting rod and tie-rod methods based on a finite element calculation result

In the presence of structural elements subjected to important point loads or presenting abrupt modifications of their section and geometry, the classical plane section methods of analysis are not satisfactory. These locations are generally reinforced using good practice rules based on experience or empirical guidelines. The strut and tie method (STM) is a rational design procedure for complex structural elements. The procedure is based on mechanics but is simple enough to be easily applied in design.

In general, the STM involves the idealization of a complex structural element into a simple structure capable of representing the path of stresses within the element.

The truss model consists of struts representing the compression fields in the concrete, ties representing the elastic steel reinforcement, and nodes representing the localized areas where the components interconnect or the areas where the elastic steel is anchored in the concrete.  The struts and ties can carry only uniaxial forces. This mechanism must be stable and maintain the equilibrium with the applied loads.

The failure of the structure is dictated by the failure of one or more ties or by excessive compressive stresses within the struts or nodes. Ideally, only the first mode of failure should occur.

Example of D region and system modeling to obtain stress fields and then loads at the ST truss supports

This method is applied to the so-called D regions. To characterize these regions, the distribution of deformations over the depth of the element is considered to be non-linear.

Consequently, the underlying hypotheses about the cross-section procedure are not valid. The principle of equivalence cannot be applied here.

According to Saint Venant's principle, an elastic stress analysis indicates that stresses due to axial forces and bending moments come near a linear distribution at a distance approximately equal to the depth of the element, h, away from the discontinuities.

In other words, a non-linear stress distribution exists in the depth of a member from the location where the discontinuity is introduced. Then, it can be stated that the D regions are assumed to extend to a distance h from the applied load and support reactions. In general, a region of a structural member is assumed to be dominated by a non-linear profile, or a D region, when the ratio extent/depth (a / h) is smaller than 2 or 2.5. The shear extent (a) is defined as the distance between the applied load and the nearest support in simple elements.

The approach consisting in defining a strut-and-tie truss can be summarized in Figure 10.

The general principle of design by STM

The verification methods are expressed in Eurocode 2 part 1-1.

The approach presented in Figure 10 is difficult to implement and depends strongly on the engineer who implements it.

More and more approaches are being developed to automate this procedure. In France, the algorithm with code_aster elements integrates a CALC_BT operator that renders the procedure semi-automatic based on:

Reference example on the left - Solution obtained automatically by the CALC_BT operator of code_aster

This method requires a high level of experience and knowledge of the field.

 

$translationBooks

Chapter F – Geotechnical calculations

Chapter F – Geotechnical calculations 

Introduction

Conventional methods for the design and dimensioning of geotechnical structures mainly aim to analyze the resistance to failure of an isolated structure. These analytical or semi-analytical methods only consider very simple geometries and provide little or no information on the ground deformations surrounding the structures.

The increasingly intensive use of urban spaces and soils, occupied by a variety of structures close to each other, makes it necessary to control the interactions between structures. The designer of a structure must justify that the displacements induced by its construction remain below the threshold set by the client. Traditional methods do not meet this requirement, which explains the increasingly frequent use of numerical modeling and software dedicated to geotechnics that are adapted to design offices. Specifically, digital modeling is used in two different situations:

Depending on the case, models with various levels of complexity can be used, considering the uncertainty of the natural soil's behavior and their spatial variability. Consequently, it becomes difficult to choose a constitutive law and determine its parameters. The user often has to choose between a robust model whose operation is well understood, but which does not account for all the complexity of soil behavior and a model that is potentially more accurate to describe the actual behavior of soils, but which includes many parameters whose roles are sometimes difficult to identify and quantify. 

The philosophy behind the design rules promotes the first approach, to use simple and robust models and ensure the safety of the design by inputting appropriate factors and carrying out parametric studies to estimate the sensitivity of the results with respect to the said parameters. However, this approach may lead to an oversimplification of the problem and can lead to overly conservative and unnecessarily costly designs.

The other approach is to use constitutive models that attempt to better represent different aspects of soil behavior. A very large number of models have been proposed, but their practical use remains difficult. Furthermore, the idea of having one "universal" model that describes all phenomena is an illusion and should be abandoned. It is necessary to define the limits of the model chosen: a perfectly plastic elastic model or one considering isotropic strain-hardening, even with a complex load surface and strain-hardening law, cannot predict a progressive accumulation of deformation in a soil mass.

Therefore, the use of finite elements requires the user to have a critical perspective on the soil constitutive model. Though, this is not the only aspect to consider: the result of the calculations can depend essentially on the 3D geometry of the structure and the geometry of the soil layers. It is sometimes (not always) justified to simplify the representation of the soil behavior and to focus on the representation of the construction process. In any case, the results of numerical simulations should be carefully verified.

Finally, it is clear that geotechnical structures are made of soil, in contact with soil, or totally or partially buried; thus, numerical modeling must represent the mechanical interaction between the soil and the structure.

This chapter is organized as follows: it summarizes the main features of numerical calculations in geotechnics before briefly presenting the verification principles adopted in Eurocode 7. Once these elements are presented, recommendations are suggested for good practice using FE calculations in geotechnics. It also deals with the features of dynamic calculations.

F.1 Geometrical aspects

F.1 Geometric aspects

F.2 Material non-linearities

F.2 Material non-linearities

F.3 Soil-structure interactions

F.3 Soil-structure interactions

F.4 Hydraulic effects

F.4 Hydraulic effects

F.5 Uncertainties and recommendations

F.5 Uncertainties and recommendations

F.6 Normative aspects: Principles of Eurocode 7

F.6 Normative aspects: Principles of Eurocode 7

F.7 Dynamic modeling

F.7 Dynamic modeling

F.8 Characteristic scales

F.8 Characteristic scales

Conclusions and references

This chapter presents an overview of the issues raised specifically by FE calculations in geotechnics, in statics or dynamics. Specific precautions should be taken when performing the calculations concerning the meshing, boundary conditions, phasing, the choice of constitutive models, and the setting of parameters.

With a certain amount of practice and the constant concern to control the obtained results, a wide range of problems and structures can be dealt with. However, there are still situations that are difficult and/or costly to deal with numerically.

Firstly, difficulties may arise due to the ambiguous mathematical structure of the studied problems. In dynamics, choosing a time discretization coherent with all the input data remains difficult. In statics, managing several non-linearities of different natures (unilateral contact, damage) may lead to numerical issues hard to overcome.

One should also remember that the physics of certain phenomena remain poorly understood and described, which is reflected in the numerical analysis. It is the case for various problems such as the initiation and evolution of landslides, implementation conditions of tasks such as the compaction of backfill materials, the initial state of stresses, water content, and other parameters linked to the history of the structural elements, etc.


$translationBooks

F1. Geometric aspects

F1. Geometric aspects 

One of the advantages of the FE method is the possibility of describing the exact geometry of the structures even during the various construction phases. CAD-like pre-processing tools make it easy to generate very complex geometries. One of the features of geotechnical structures is that, generally, the earth mass that constitutes or surrounds the structure must be considered within the mesh.

  1. Boundaries of the studied domain

The first challenge is to identify the boundaries of the domain considered for the analysis. For a geotechnical structure, the horizontal and the lower boundaries of the studied domain are rarely determined with precision: the extent of the domain is then bounded by vertical planes whose position is generally fixed by using empirical rules.

Using the plane strain assumption, for example, the position at which the lower limit of the mesh is set has a direct impact on the calculated settlement for a strip footing or above a tunnel. This relation is clear in the case of a homogeneous linear elastic soil mass. It can be reduced by taking into account elastic moduli, which increases with depth. Nevertheless, it is likely to induce a significant error in the calculated displacements. The ideal case is when a rigid bedrock has been found at a given depth, which assumes that the research has been carried out to a sufficient depth. To model a tunnel, for instance, it is advised to collect soil sampling way beyond the depth of the shaft, which is generally not the case in real projects.

In the lateral directions, considering a domain that is too small can also significantly modify the response of the numerical model. Fixed displacements lead to an overestimation of the stiffness of the solid, and "smooth contact" type conditions lead to an overestimation of the displacements. The choice of the mesh dimensions adapted to a structure remains a largely open problem even if some authors have suggested practical rules. However, those rules should not be taken as absolute prescriptions (see Mestat and Prat, 1999).

Consequently, choosing the size of the domain considered for the mesh is an important step in the modeling of geotechnical structures, even for relatively simple static analyses. In the case of dynamic calculations, the question of the size of the meshing domain raises specific issues and is a crucial part of the modeling strategy. This will be detailed in section 8.

  1. Soil heterogeneities

In some regions, such as London or Frankfurt, geology allows the soil in the vicinity of the structure to be considered homogeneous (in the sense that its mechanical and hydraulic behaviors can be represented by a single model). However, in other contexts, particularly in the Paris region, it is common for the studied area to include several layers of soil of very different natures and characteristics (particularly mechanical). Therefore, the elaboration of a computational model begins, as for traditional methods, with a detailed study of the soil layers in the area of interest. The goal is not to reproduce the exact geometry of the geological layers (which may be locally thin), but to define geotechnically homogeneous sets.

  1. Discontinuities

An important particularity of geotechnical calculations is the presence, within the solid, of failure planes that existed before the construction of the structure or the implementation of the studied load. They produce a discontinuity of displacement between the blocks located on either side of the failure plane. The FE method is well adapted to search for continuous displacement fields, and the consideration of this type of discontinuities requires the implementation of special techniques (specific elements are generally used) or even the use of another calculation method such as the discrete element method.

Behaviors such as landslides are extremely complicated to predict. Nowadays, it is extremely difficult to predict the occurrence and development of failure planes. Usually, the engineer is forced to consider an existing plane whose position has been identified with more or less precision using the adapted equipment (for example, inclinometers that follow the deformation of a slope).

Solid rocks often contain a large number of fractures with directions that are seemingly parallel to one or two main directions (local diaclases). The fracture distribution is diffuse, and random, or at least impossible to characterize completely at the scale of the solid. If the global behavior of the solid is of interest, it can be modeled as a continuous medium employing a constitutive model, which incorporates, at the calculation scale, the effects of the discontinuities: homogenization methods are generally used.

One may also have to consider large discontinuities in a solid bedrock (fractures), which can be treated as the failure planes of landslides.

  1. An "open" material system and construction techniques to be modeled

As stated in the introductory chapter, a peculiarity of civil engineering calculations is the need to consider construction phases such as clearing or backfilling, implementing a concrete mass, tensioning cables, etc. Taking these construction phases into account within the framework of the FE method can be a lengthy and complicated process: the method consists of reducing the problem to the calculation of a stiffness matrix, a vector of nodal forces, then solving the system obtained by considering the boundary conditions. To model the construction phases, one must perform a series of calculations taking into account the change in stiffness of some elements, disappearance or modification of some supports, changes in the point of application of the loads, etc. 

Therefore, the issue is to propose simulation techniques that can take into account a large number of constructive arrangements within the relatively narrow framework of the FE method. It is up to the user to decide whether the modeling tools proposed by the software correctly reflect the phenomena involved.

In the case of tunnels dug by TBM, the stresses applied onto the ground during the various digging phases are complex as the soil is closing in on the ring formed by the segments as the TBM advances. Pile driving is another example, which is difficult to simplify to the FE method framework.

 

$translationBooks

F2. Material non-linearities 

F2. Material non-linearities 

In geotechnics, it is very rare to be able to study a structure exclusively using linear behavior assumptions (except for certain dynamic analyses). It may be useful to first perform a linear calculation, to check that the geometry and boundary conditions are correct and get a preliminary idea of the deformation that the loading may cause. Nevertheless, the results should be put in perspective as they can be completely wrong in some cases: for an excavation in front of a cast wall for example, the kinematics calculated with linear behavior are clearly different from the one observed.

  1. Constitutive laws

Even when limited to cases that are perfectly saturated or perfectly dry, the behavior of soils is complex. In practice, elastoplastic models are most often used, which give a force-displacement relationship that is non-linear but independent of time. The effects of creep and viscosity can be considered for particular applications when required - and if the corresponding parameters can be obtained experimentally - but the use of such models remains limited.

Among the elastoplastic models, perfectly plastic linear elastic models are still most often used (see the survey cited by Gilleron, 2016). The use of non-linear elastic models associated with one or more plastic collapsible mechanisms is gradually becoming widespread, particularly under the influence of software developers: in some cases, they give results that are much more representative of reality (for example, for excavation in front of a molded wall), but the influence of each parameter of the models is not necessarily well controlled. Generally, the choice of a constitutive model for soils must consider the objectives set for the calculation, the type of structure (and the type of stress to which the soil will be subjected), the level of precision of the available reconnaissance and laboratory tests.

  1. Initial state

For non-linear models, the stiffness of the material depends on the initial state of the stresses in the studied soil mass (and possibly other parameters, such as strain-hardening or damage). Therefore, the determination of the initial stresses has a decisive influence on the results. Unfortunately, the initial stresses are generally evaluated in a very simple way: they are either assimilated to a "geostatic" stress field (for a solid with a horizontal surface), or obtained by applying gravity to the whole mesh from a state of zero stress. These processes are relatively poor compared to the complexity of the rheological models used for soils. However, they remain unavoidable in practice, for lack of a better way of estimating initial stresses in the soil.

 

$translationBooks

F3. Soil-structure interactions

F3. Soil-structure interactions

Geotechnical structures often combine layers of soil with metal or concrete structures, which are generally much stiffer than soil. The interaction may be limited to a few points where the structure rests on the ground, or it may be continuous over a significant surface, such as the top surface of a tunnel, or retaining wall.

The interaction will be treated more or less precisely depending on the case.

In the case of tunnels, for example, it is very common to consider a perfect bond between the ground and the vault because of the construction method: with traditional (sequential) methods, the retaining wall is made of concrete sprayed directly on the surface of the ground uncovered by the excavation, which in principle ensures good continuity of movement. During tunneling performed by a TBM, efforts are made to ensure good force transmission between the ground and the segments by carrying out backfilling injections into the space between the arch and the ground. However, in older tunnels, the pathologies observed (or the tests carried out in situ) may suggest that contact is locally lost between the vault and the ground. For instance, it can be due to water seepage that may have washed out the ground: the modeling must then describe more precisely the contact conditions between the ground and the vault.

Covered trench tunnels present a different problem, as the ground around the structure is backfilled. The modeling of this operation may require explicit consideration, using specific elements, of the interface between the vault and the ground.

It is also common to introduce explicit modeling of the interface between the soil and the structure for retaining structures, when backfilling behind a wall (the phenomenon of soil sliding at the interface with the structure being similar to that involved in covered trenches), or when excavating in front of a cast wall, for example, the mass of soil supported may slide and present a vertical displacement greater than the height of the wall.

The question of modeling an interface between soils and structures must be considered on a case-by-case basis. One can introduce contact or interface elements specifically intended to represent the mechanical interaction between the two, but these elements introduce new parameters, which can be difficult to identify (such as the normal and tangential stiffnesses of the interface). This modeling approach presents a risk: the interface elements tend to control the behavior of the structure and to blur the role of the soil behavior, giving the impression that the response of the structure hardly depends on the soil anymore.

Reinforced structures

In many cases, the soil is reinforced by inclusions with very high stiffness and strength characteristics. These inclusions are discreetly distributed in the soil and very slender: piles, micro piles, tie rods, wall reinforcement in reinforced soil. This geometrical particularity poses various difficulties. First of all, strictly speaking, a row of piles is not equivalent to a continuous wall, and the use of plane strain assumption is not justified. In practice, one would adopt, for the planar calculations of the wall, mechanical characteristics "equivalent" to those of a row of piles, using assumptions that can be more or less difficult to justify. The same applies to the parameters of the mechanical interface between the soil and the piles/wall. The difficulty is the same if the wall is represented by surface elements or by linear beam-type elements.

To overcome this difficulty, 3D modeling can be used. However, because of the dimensions of the cross-section of the inclusions, whenever there are more than a few units, it becomes impossible to represent in the mesh the real geometry of the inclusions: for a reinforced earth wall, with reinforcements of 5 mm x 45 mm section, at a rate of 4 to 6 reinforcements per 0.75 m x 0.75 m, and for a volume whose dimensions are of the order of ten meters, the number of nodes of a mesh that would respect the real geometry of the inclusions and that would give an acceptable discretization exceeds the current calculation capacities. We can therefore propose to represent the inclusions by 1D elements (with or without considering bending effects). This approach is questionable from a theoretical point of view because the introduction of a linear density of force exerted by the inclusion in a 3D medium is not compatible with the classical representation of internal forces by a stress tensor. It can be used, however, one must be careful in the interpretation of the results, at least with respect to the stresses in the vicinity of the inclusions.

An alternative solution is to adopt homogenization-type approaches to take into account the influence of inclusions on the overall behavior of the structure. More or less complex models have been developed and implemented in some software.

Whatever the choices made (calculation in plane strains or in 3D condition, discretization of inclusions - by linear or non-linear elements - or homogenized approach), it is necessary to represent the mechanical interaction between the pile and the soil that occurs at the contact between the soil and the sidewall of the pile, and also between the soil and the footing of the pile. Modeling the mechanical interaction at the footing of the pile is particularly difficult to master. Modeling a single pile by volume elements, possibly with interface elements with the surrounding soil, gives results that depend on the mesh and the constitutive model used for the soil. It is necessary to use a model that reproduces the soil failure in compression if one is interested in the pile failure. For instance, models such as Mohr-Coulomb or Drucker Prager are not suitable in this context.

In some 2D or 3D models where inclusions are represented by bar or beam elements, a fictitious end (e.g. a horizontal beam element perpendicular to the pile) is associated with the inclusion, in an attempt to better represent the interaction at the footing of the pile: performing studies to determine the sensitivity of the results to the dimension of the added elements seems appropriate to verify the relevance of this approach.

Finally, other modeling techniques are available, which propose to explicitly integrate an interaction model for lateral friction and another one for peak interaction via ad hoc elements.

Without going more into detail, let us highlight the fact that the user is free to choose between numerical simulation techniques and models, which have a direct influence on the obtained results.

 

$translationBooks

F4. Hydraulic effects

F4. Hydraulic effects

  1. Hydromechanical coupling

Another peculiarity of geotechnical calculations concerns the role of water in soils. When a mechanical load is rapidly applied to a saturated layer of soil, instantaneous deformation, and pressurization of the fluid in the vicinity of the applied load occurs. Depending on the hydraulic boundary conditions, the gradient of the hydraulic load causes the fluid to be set in motion, leading to pressure redistribution and delayed deformation of the soil.

Therefore, a problem of hydromechanical coupling must be dealt with. A solid theoretical framework was established by Biot (1941) and developed by Coussy (1991). In terms of numerical resolution, the coupled problem is much more difficult to deal with than a classical problem, for several reasons:

Complete processing of the hydromechanical coupling is rarely performed. One generally tries to limit oneself to a decoupled approach to the problem, in which one calculates the evolution of the pressure field while neglecting the deformations of the solid. However, this decoupling has the consequence of strongly underestimating the duration over which the pressure redistribution occurs.

  1. Unsaturated soils

The treatment of unsaturated soils complicates even more the problem as the transition from near-saturated to unsaturated zones introduces additional unknowns, non-linearities, and parameters. Again, the complete treatment of unsaturated soils remains rare. It is preferred to propose simplified solutions, neglecting partially saturated zones for example.

Defining the initial state in the case of unsaturated soils is extremely complex because of the lack of experimental methods to characterize it in situ.

 

$translationBooks

F5. Uncertainties and recommendations

F5. Uncertainties and recommendations 

  1. Uncertainties

There are many sources of uncertainty in geotechnics. The first is the relative lack of knowledge of the geometry of the soil layers that make up the solid being studied. In addition to the geological knowledge about the area where the structure stands, the main sources of information are the boreholes drilled on the project site. The extrapolation between boreholes, especially if they are far apart and not located at the exact location of the future structure (which may be modified after reconnaissance, for example in the case of tunnels), does not necessarily give an accurate representation of local variations in layer thickness.

In urban areas, the presence of heterogeneities (cellars, wells, foundations of previous structures) is often difficult to detect.

The other source of uncertainty already mentioned concerns the initial state of stresses (and possibly pore pressures) in the soil mass. It can have a major influence on the results of the calculation: this is particularly clear in the case of tunnels, where the loads considered depend on the initial stresses.

Finally, the choice of constitutive models and the determination of the parameters of these models introduce a significant uncertainty on the representativeness of the calculations: if the constitutive model does not capture a phenomenon that controls the behavior of the structure, the result may be qualitatively and quantitatively very far from reality.

  1. Recommendations

In general, the user must be aware of the objectives of the calculation he is undertaking: the approach is different depending on whether one is trying to justify dimensions or to evaluate the influence of certain constructive provisions (the number and position of the struts, for example).

One must also be aware of the modeling choices on which the calculation is based (even if these choices are sometimes partly imposed by the software). One must be able to identify the phenomena to be taken into account, which leads to the choice of a quasi-static or dynamic analysis, the consideration or not of the hydromechanical coupling, etc.

One must choose between a 2D or 3D calculation. Three-dimensional calculations remain rare for the moment because of the time required to prepare the calculations. However, for some problems, it is clear that two-dimensional calculations can only give a poor indication of the behavior of the studied structure regardless of how long the engineer took in determining the soil parameters. For instance, the study of the stability of the tunnel face cannot be considered outside a three-dimensional context. The same is true for the study of bolt reinforcements of the tunnel face. The development of pre-processors specific to each application should simplify the use of 3D calculations and improve the representativeness of many finite element analyses.

In geotechnics, special focus should be given to the choice of soil parameters: it could be the subject of a whole book. Most advanced constitutive models do not come without a detailed and robust parameter identification procedure, mainly because the model equations cannot be solved even for a simple problem such as triaxial compression. One must therefore calibrate the parameters so that the modeling of triaxial tests, for instance, outputs results in satisfactory agreement with the test results. A trial and error procedure is used, and since the agreement obtained is evaluated subjectively (because we can choose to better reproduce one part or another of the experimental curves), it is not guaranteed that two users will obtain the same parameter values from the same tests. Thus, all constitutive models do not have the same qualities. Some have many parameters, each of which influences a particular aspect of the soil response (but which does not necessarily appear in the available test results). Other models, on the contrary, have a relatively small number of parameters, but each of them can simultaneously modify several aspects of soil deflection. It makes the calibration much more complex. 

The last recommendation that must be kept in mind is to check the calculations as much as possible. There are not yet any general tools to measure the quality of a calculation: research work aims to provide error estimators, but their use in geotechnics remains rare. It is therefore necessary to look in-depth at the results: some inconsistencies are sometimes easy to detect. If there are any doubts, it is useful to have the results checked by someone else. In any case, it is highly recommended to carry out parametric studies to get an idea of the influence of certain factors, in particular soil parameters, if it is not certain that their influence on the results is moderate and that their value has been determined with acceptable accuracy.

 

$translationBooks

F6. Normative aspects: Principles of the Eurocode 7

F6. Normative aspects: Principles of the Eurocode 7

Numerical modeling always had a particular relationship with calculation standards, which essentially focus on the verification of ultimate limit states by calculating safety coefficients or the equilibrium of forces integrating partial coefficients. Indeed, numerical modeling provides first and foremost values of displacements and deformations and is therefore a tool very well suited to the verification of serviceability limit states (SLS) for which the partial coefficients are equal to 1 Its use for the verification of ultimate limit states (ULS) therefore initially appeared to be limited. From now on, particularly with procedures aiming to reduce the shear properties of soils (for example, the "c-phi reduction" procedure), it is easy to calculate a safety coefficient. It is also possible, through the procedures suggested by certain calculation standards, in particular Eurocode 7, to consider the results of numerical modeling both for the verification of the ultimate limit states and for that of the serviceability limit states.

However, Eurocode 7 in its current version is not very clear about how the analysis and processing of the results using numerical modeling should be conducted. Indeed, the three calculation approaches proposed by Eurocode 7, which allow applying partial coefficients on the actions, the effects of the actions (bending moments and shear forces in a retaining wall, axial force in a pile, etc.), the properties of the soils (c and or cU), and the geotechnical resistances, have been designed to be used with limit equilibrium methods. The verifications related to the application of the Eurocodes are mainly transcribed in the form of comparisons between actions or effects of actions and resistances. Nevertheless, several publications (Potts and Zdravkovic, 2012, Tschuchnigg et al., 2015, etc.) and reports of different working groups convened for the development of the second generation of the Eurocodes make it possible to outline verification procedures. First of all, it is important to underline that the weighting "at the source" of the soil or rock properties, i.e. the possibility of reducing the cohesion and the angle of friction before carrying out the calculation, is not allowed because it leads to results that cannot be interpreted. Different procedures can however be used and are summarized in the table below.

Synthesis of the different types of ULS verifications

Type of ULS verifications 

1 - Structural ULS

2 - Geotechnical (and structural) ULS

3 - Geotechnical ULS

4 - Structural (and geotechnical) ULS

Type of procedures

Multiplication of the effects of the actions by 1.35

Reduction of the shear properties of the soils

Estimation of the resistance centered around a specific geotechnical structure (pile, tie rod, etc.)

Increase of the loads applied to the geotechnical structure

Comments 

Must be combined with the verification of geotechnical ULS

Must be combined with the verification of structural ULS

Must be combined with the verification of structural ULS

A priori, this approach is sufficient on its own

The approach that tends to be used is to carry out a calculation with characteristic values for both the loads and the properties of the soil to reach the first state of equilibrium. Concerning this state of equilibrium, two types of verifications must be carried out. The first verification, namely type 1, is relative to the structural ULS and consists of multiplying the effects of the calculated actions by 1.35 (i.e. the coefficient usually considered in the Eurocodes for unfavorable permanent loads). The second verification, type 2, is related to geotechnical ULS and consists of progressively reducing the shear resistance properties of the modeled soil to reveal a failure mechanism. The reduction factor applied can be considered as a safety factor, but its interpretation can be subjected to a discussion. Indeed, during the reduction of shear properties, different elements may interact such as strain-hardening functions or flow rules. In the case of calculations associating volumetric elements and structural elements such as "bar" or "beam" type elements, this verification can also be used to check structural ULS. The displacements accumulated during the procedure of reduction of the shear resistance properties generate forces in the structural elements whose interpretation is a matter of debate among the engineering researcher’s community. Therefore, the user must be extremely careful about the way to implement this type of procedure. Figure 7.1 shows the sequence of type 1 and 3 verifications in the case of a phased calculation.

Geotechnical ULS can also be specifically apprehended by considering the resistance centered around the geotechnical structure to be dimensioned, a pile, an anchor, etc. and comparing it to the mobilized resistance provided by the standards, for instance, the pressuremeter norms for the calculation of the bearing capacity of piles: this is the type 3 procedure and it must be associated with the type 1 for structural ULS.

In certain configurations, specifically for piles or footings, it is also possible to increase the actions applied to the geotechnical structure to be dimensioned to directly obtain the maximum applicable force and to deduce the safety coefficient: this is the object of the type 4 approach.

More generally, the methods used to justify geotechnical ULS must allow either to identify the failure mechanism and verify that there is a sufficient margin to trigger the mechanism (this is the object of types 2 and 4 verifications) or to compare the mobilized resistance with a mobilizable resistance that can be calculated elsewhere (this is the purpose of type 3 verifications).

Sequencing of the ULS verifications in a phased computation

When the failure mechanism is associated with the action of water special precautions must be taken. Indeed, the procedures described above do not enable us to correctly compute the safety coefficient if the failure mechanism results exclusively from the action of water. Firstly, it is necessary to define whether the failure mechanism is related to a mechanical equilibrium defect (for example, water rising too high behind a retaining wall), or to a hydraulic problem (for example, a hydraulic gradient being too high). In the case of a mechanical equilibrium defect, it is necessary to carry out a parametric study relative to the level of the water table to estimate its effect by implementing the ULS verifications presented in the previous table. In the case of a hydraulic problem, it is possible to compare the calculated gradient with a critical gradient calculated elsewhere and to verify that the effective stresses remain positive at any point of the solid (because the condition on the critical gradient is more conservative than the condition on the effective stress).

 

$translationBooks

F7. Modeling in Dynamics

F7. Modeling in Dynamics

 

Because of the continuous and unconfined nature of geomaterials, the treatment of dynamic problems in "soil" media differs from the treatment of dynamic problems associated with classical civil engineering structures. While for most of the latter, mechanical modeling through the assembly of masses and discrete sources of stiffness seems sufficient to understand their dynamic behavior, geomaterials must be handled as continuous and unconfined media. Their dynamic response must be studied as a wave propagation mechanical problem.

Concerning the numerical methods used in soil dynamics, it should first be noted that modeling wave propagation requires to use of a range of modeling and calculation methods, which differ from conventional geotechnical calculations using the FE method. Although the FE method has an important place in the inventory of methods available to the engineer for soil dynamics problems, it is often coupled with, or replaced by, other numerical computation methods that are better suited to model wave propagation in unbounded environments.

Additionally, the equations that describe dynamic problems involve a new variable (apart from spatial variables): time. The computational algorithms are formulated in the time domain, and alternatively, in the frequency domain. The methods used to solve soil dynamics problems are implemented in specific software that requires the user to have a very good understanding of the underlined theory, and specific skills to determine the parameters, the calibration of numerical models and constitutive laws, the optimization of calculations, etc. 

Among the common soil dynamics problems, the following ones can be cited:

 

$translationBooks

F8. Characteristic scales

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 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).

  1. 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:

  1. 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, 

  2. 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 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:

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.

  1. 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:

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.

  1. 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 Vmin that needs to be considered in the model and the maximum characteristic frequency fmax 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 Vmax 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 Vmin = 150 m/s, Vmax = 500 m/s and fmax = 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 Vmin) evolves with time (softening of the soil by passing in elastoplastic regime).

  1. 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:

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.

  1. 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:

  1. 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.