1    Introduction

In geotechnical and rock engineering, accurately simulating the behaviour of rock or soil using Finite Element Analysis (FEA) requires careful stress initialisation. The stress state of the rock or soil at the beginning of the analysis greatly influences the subsequent deformation and stress distribution throughout the simulation. Therefore, correctly estimating the in-situ stresses is crucial.

As a first approximation, it can be assumed that the initial vertical stresses, σ_v0, are primarily influenced by the weight of overlying materials, such as soil, rock or water bodies. They are assumed to increase with depth due to the cumulative weight of the overburden material. Conversely, estimating the initial horizontal stresses, σ_h0 can be more intricate, as they depend on the material and the history of its formation.

In the case of rock materials, vertical and horizontal stresses often demonstrate varying magnitudes. The ratio of horizontal stress to vertical stress, known as the coefficient of stress, is subject to variation due to factors such as lithology, tectonic setting, and stress history.

2    Stress Initialisation

In PLAXIS 2D and PLAXIS 3D, users have the flexibility to generate initial stresses through three distinct methods: the K0 procedure, Field stress or the Gravity loading option. These methods offer enhanced control over the in-situ stress field, allowing for more complex problems to be simulated. Below is a description of each of the in-situ stress initialisation methods.

It is important to note that these options are accessible in the Calculation type drop-down menu exclusively in the Initial phase (Figure 1).

Figure 1: Available options for the initial stress generation

2.1   K₀ – procedure

The K₀ – procedure is a calculation method available in PLAXIS for defining initial stresses in a model, considering the loading history of the soil. The procedure involves defining a coefficient, typically denoted as K₀, which represents the ratio of lateral stress to vertical stress. This coefficient characterises the initial stress conditions within the soil or rock mass.

The initial vertical stress, σ_v0, is determined by the weight of the overburden rock, known as overburden pressure:where,

• γ represents the unit weight of the rock material, and
• h denotes the depth

In the presence of porewater pressure, p_w, within the soil or the rock mass, the mechanical behaviour is governed by the effective stress, σ’_v0, instead. This effective stress is related to the total stress by Terzaghi’s principle:

where,

• p_w represents the porewater pressure

The initial effective horizontal stresses, σ’_H and σ’_h (Figure 2) are determined based on the user-defined K_0,H and K_0,h, according to the relationships in the following Equation:

where,

• K_0,H is the coefficient of earth pressure at rest (in horizontal direction 1)
• K_0,h is the coefficient of earth pressure at rest (in horizontal direction 2)

Figure 2: Indication of vertical and horizontal stresses

PLAXIS provides two methods for defining K₀:

• User-defined: You have the flexibility to directly input K₀ values (K_0,H and K_0,h) based on your knowledge of the in-situ conditions.
• Jaky's formula: PLAXIS offers automated computation of K₀ using Jaky's formula, which correlates K₀ with the soil's friction angle (φ): K_0,H = K_0,h = 1 - sin(φ).

Both methods allow the specification of two K₀ values: one for the x-direction and one for the z-direction (PLAXIS 2D) or y-direction (PLAXIS 3D). Following the notation of Eq. 2, these values are denoted as K_0,H and K_0,h, respectively:

Once established, PLAXIS utilises these K₀ values to generate initial stresses throughout the soil model. This process involves:

• Calculation of the vertical total stresses, σ_v, based on the summation of weights above the stress point (referenced by Eq. 1a).
• Calculation of the vertical effective stresses, σ’_v, by subtracting the porewater pressures (referenced by Eq. 1b)
• Calculation of horizontal effective stresses, σ’_H and σ’_h, by multiplying the K₀ values with the vertical effective stresses. This ensures a specific relationship between horizontal and vertical stresses, reflecting the at-rest condition of the soil (as described in Eq. 2).

Note:

• In a K₀ – procedure the stress field equilibrium is not checked at the end of the initial phase.

2.2   Field Stress

Field stress allows for setting up a homogeneous initial stress state in the model, with eventually the consideration of a rotation of principal stresses. This might be relevant for applications in deep soil or rock layers the geological formation of which has caused a rotation of principal stresses (shearing).

In addition to the selection of Field stress as Calculation type for the initial phase, the user needs to define the magnitude of the three principal stresses σ₁, σ₂, σ₃, as well as their respective orientation (Figure 3).

Note that in PLAXIS, field stress can be defined globally for the entire model or locally for individual rock clusters. This latter option is particularly useful for scenarios involving rock layers with disturbance zones, where the weaker material may have different initial stress states compared to the surrounding rock mass.

Figure 3: Field stress in Model explorer

Notes:

• The main interest of the field stress option is its simplicity for defining uniform state of stress without consideration of the overburden.
• The calculation type Field stress does not consider the increase of stresses with depth due to gravity. Therefore, the ΣMweight multiplier is always set to zero when using the Field stress option.
• The Field stress option may generate shear stresses along model boundaries in case principal directions are not aligned with the global axis directions, which cannot be resisted by the default boundary conditions. Therefore, the user should in that case manually set all global model boundaries (BoundaryXMin, BoundaryXMax, BoundaryYMin, BoundaryYMax, BoundaryZMin, BoundaryZMax) to fully fixed in the Deformations section of the Model explorer.

2.3   Gravity loading

Gravity loading is a type of Plastic calculation, in which initial stresses are generated based on the volumetric weight of the soil. If Gravity loading is adopted, then the initial stresses are set up by applying the soil self-weight in the first calculation phase starting from a zero state of stresses. This is achieved by setting ΣMweight = 1.0.

Gravity loading provides limited control over the initial stress ratio (K₀) due to its dependence on the chosen constitutive model and its parameters. For example, in the Mohr-Coulomb model (an elastic-perfectly plastic soil model), K₀ is highly sensitive to the assumed Poisson's ratio (ν). In one-dimensional compression, elastic calculations yield:

This formula offers a simplified picture, assuming the immediate application of gravity to an elastic material with no lateral movement. In reality, complex factors like locked-in stresses from tectonic shifts or overburden removal, render this scenario rare. Therefore, selecting appropriate Poisson's ratio values is crucial to achieve the desired initial stress ratio.

As Poisson's ratio must be lower than 0.5, it is not straightforward to generate K₀ values larger than 1 using Gravity loading. As such, Gravity loading provides limited control over the initial stress ratio, especially for rock material where K₀ values larger than 1 are often desired and typically the initial horizontal effective stresses are not equal, due to the material anisotropy. If K₀ values larger than 1 are desired, it becomes almost necessary to use the K₀ – procedure (see Section 2.1).

3    Considerations for stresses in rock materials

3.1   Simulating large overburden pressures

Large overburden pressures in deep rock structures, such as mines or tunnels, can be challenging for engineers. One way to address this issue is to simulate the overburden at a smaller scale in the project. For example, you can create a thin layer to represent the omitted overburden and assign a unit weight estimated from the overburden. This can help to compensate for the missing soil weight and avoid generating non-realistic stresses.

 

Figure 4: Simulating overburden pressures using a virtual layer

As shown in Figure 4, the example considers that the model is at a depth of 130.0m (h_real) and the unit weight is γ_real = 18.0 kN/m³. To simulate the overburden pressure of this soil layer, while reducing the size of the model geometry, we can introduce a virtual layer of 1.0m (h_virtual) on top of the existing model and assign a unit weight estimated from the overburden, which will be:

3.2   Simulating locked-in stresses

Another challenge is the consideration of the specific in-situ stress conditions, which can vary significantly from one geological formation to another. This can lead to initial locked stress in certain rocky layers that must be accounted for.

For example, in areas with active tectonic processes, such as fault zones or regions near plate boundaries, understanding the locked-in stress components is essential for assessing geological stability. These stress components can influence fault behaviour, rock deformation, and the likelihood of seismic activity. Also, in civil and geotechnical engineering projects, particularly those involving underground structures like tunnels or mines, the consideration of locked-in stresses is vital for ensuring structural integrity and safety. These stresses can exert significant pressure on surrounding rock or soil, affecting excavation, support system design, and overall project stability. Therefore, in both geological and engineering contexts, proper consideration of locked-in stress component (σ_tect) is crucial.

To quantify the initial stress conditions within geological formations, we employ the following equations, denoted as (6a) and (6b), which delineate the distribution of stress. These equations provide valuable insights into the relationship among vertical stress (σ′_v), horizontal stresses (σ_H' and σ’_h), and the locked-in stress component (σ_tect). Equation 6a illustrates that horizontal stresses primarily comprise two components: one originating from gravitational stress and the other from tectonic stress.

When considering a stress ratio, β, between the major and minor horizontal effective stresses, Equation 6a can be expressed as:

Although, there is no straightforward way to apply this stress component of tectonic origin in PLAXIS, we can employ the following strategy:

1. Define the initial stresses according to what was discussed in the Chapter of Stress Initialisation, without taking into account the stress σ_tect.
2. Introduce σ_tect by applying a prescribed strain leading to the desired level of tectonic stress. Namely:

           or

Notes:

• The application of the prescribed strain should be done in combination with a “dummy” linear elastic material. This material has zero Poisson ratio (ν = 0) and a Young’s modulus equal to the stiffness of the rock (E = E_rock).
• In case of multiple layers of rock, the same technique can be employed with the addition of frictionless interfaces between the layers where the prescribed strain is applied avoiding the generation of spurious shear stress around their contact surfaces due to discontinuous prescribed strain field in each rock layer

3.2.1  Solved Example

The aim of this example is to illustrate the modelling technique described in the previous Chapter. For more information, please refer to the "Downloads" section below.

Downloads

• 
•  PLAXIS 2D 2024.1 - Initial Stress Definition with Consideration of Locked-in Stress [command log file]
•  PLAXIS 3D 2024.1 - Initial Stress Definition with Consideration of Locked-in Stress  [command log file]