# The Optimization Algorithm for Large-Scale In Situ Stress Field

Authors: Fei Li; Ningdong Chang
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

In situ stress state is a predominant factor for the design and safe construction of geotechnical engineering. For a real construction site, the amount of calculation using a finite element method for in situ stress field increases dramatically with the increase of the calculation freedom due to large-scale uncertainties. In order to reduce the computing cost without losing the accuracy of the calculation, an optimization algorithm combined with a reduced order model, which is realized by the proper orthogonal decomposition algorithm (POD) for large-scale in situ stress field, is put forward in this paper. The POD algorithm produces a set of orthogonal bases through the extraction of the field variables, combining with the Galerkin finite element method to create a reduced order numerical model. The reduced order model is then calculated with a global optimization algorithm to inversely find the solution for the actual in situ stress field. In order to verify the accuracy and efficiency of the method, two examples are presented to simulate the inverse calculation of the in situ stress field. They showed that the computation time of the POD method could reach 1/10 of the ordinary computation time. Also, the results showed good accuracy with a minimum computational expense, which can provide a reference for inverse calculation of large-scale in situ stress field.


---

## Body

## 1. Introduction

Many rock engineering processes, especially deep mining and slope, are carried out in complex geological bodies (physics field) and geological environments (mechanical field) [1–4]. In situ stress guides and determines the potential occurrence characteristics and degree of a geological disaster. Estimation of the initial in situ stress field accurately in the stage of preplanning and design of construction is to correctly understand and fully grasp the status and regularity of in situ stress for complex geological body, so as to predict and judge the type and scale of major engineering geological disasters [5]. Furthermore, safe technology solutions and disaster prevention and mitigation measures can be made tentatively to weaken and eliminate the potential hazard in the premise [6].At present, the inverse calculation of the in situ stress field is based on two ideas [7]. One is relied on the displacement measured data with consideration of mechanical characteristics of the local rock mass to calculate the quantity and direction of in situ stress [8]; the other is according to the measured stress of certain control points in actual engineering site, creating a numerical model considering topography and geology, lithology, and so on, then using an inverse algorithm to calculate the in situ stress field [9].In previous work [10], a surrogate model accelerated random search algorithm (SMARS) was proposed in the inverse calculation of in situ stress field, which relies on a random search algorithm to keep global search capabilities and surrogate model to accelerate convergence to a final solution. However, the calculation area is large in general geotechnical engineering; and in order to reflect the local landform, the geological model should try to restore the characteristics of actual construction sites; moreover, for the purpose of improving the calculation precision and reflecting the local deformation characteristics, local mesh encryption is done in most cases, which will all result in a huge mesh and slow calculation [11]. Although, the SMARS algorithm can ease the computational expense for the inverse problem of the in situ stress field by reducing the number of forward calculations. The expensive numerical model of geotechnical engineering problem is inevitably calculated a certain number of times.Therefore, an important problem is how to simplify the calculation process and save computing time and memory capacity in the premise of sufficient accuracy. The proper orthogonal decomposition (shortly for POD) is such an efficient dimension reduction method to do a multivariate statistical analysis for data compression, which can give the optimal low dimensional approximation for multidimensional physical processes (that is, optimal series expansion) in the least squares sense [12, 13]. It aims at describing the multidimensional physical process as low dimension approximation with less degree of freedom to express the characteristics of objects, arriving at the purpose of reducing the amount of data storage required for reproducing the physical processes, saving computation time and computation load. The POD method relies entirely on the raw data without any prior assumption. The effect of the POD method is mainly reflected in two aspects: (1) reducing the order: the multidimensional data can be approximated by low dimension data, that is projected to a low dimensional space; (2) feature extraction: the POD method can turn a group of variables to uncorrelated variables by solving orthogonal basis, and as far as possible to keep all information of the original variables [14]. It provides an effective method to describe a physical process. That is expressing the physical process as a sequence of expansion of a time dependent function and a space dependent function, which is optimal in the mean square sense with only a small amount of items.The POD method was first proposed by Loeve (1945) and Karhunen (1946) and has been successfully applied in many different fields [15]. At present, POD is widely used in fluid mechanics, data compression, image processing, signal analysis, pattern recognition, optimization and control of power system, the supersonic jet impact model, thermal processing of food, and the effect of dynamic wind pressure on the building and so on [16]. In recent years, the POD method is also successfully applied to dimension reduction solution and error estimation of the finite difference method and mixed finite element method for the tropical Pacific Ocean Model and finite difference method for non-Navier-Stoke equations [17]. Also, the strength of aeolian sand was predicted by using the PCA-BP neural network model and BP neural network model, which showed the PCA-BP neural network model had high accuracy and good effect [18]. What is more, the POD method can also solve the inverse problem. For the application of order reduction for the finite element method, the POD method is used to create a group of orthogonal bases to replace the polynomial interpolating functions in the Galerkin projection finite element method, reducing the dimension of the partial differential equation, which is formed according to an actual physical and mechanical problem.In the paper, the inverse problem and calculation parameters are described. Then, the formulation of the in situ stress calculation finite element method is presented. Finally, two examples are shown to make the comparison of efficiency and accuracy of the inverse in situ stress calculation problem by the full order model and reduced order model separately.


## 2. Inverse Problem of In Situ Stress Field Calculation

The in situ stress is mainly composed of gravity and tectonic stress [19]. The gravity can be simulated by assigning density to rock mass accordingly. The formation of the tectonic stress field is due to the movement of the plate, namely extrusion tectonism, shear tectonism in the horizontal direction, and shear tectonism in the vertical direction. So, in order to simulate the tectonic stress field, the corresponding boundary conditions are applied to a numerical model, simplified as horizontal pressure, horizontal shear stress, and vertical shear stress. For general engineering sites, these boundary conditions can be assumed as linear distribution. As shown in Figure 1, the bottom boundary and the right boundary of the model are fixed, and the horizontal pressure on the left boundary is applied in Figure 1(a), the hydrostatic pressure on the left boundary is applied in Figure 1(b), and the vertical shear stress is applied in Figure 1(c).Figure 1
The assumed distribution of boundary load, which contributes to the tectonic stress field.
(a)(b)(c)The core problem of inverse calculation of in situ stress field can be recast as an optimization problem, which is to construct a group of boundary conditions (stress or displacement) applied to the numerical model to generate the calculated stress, then the error between the calculated stress and the measurement stress is minimum. Therefore, the objective function for optimization is as follows:(1)minrl⟶=∑1nfim−fiexpl⟶L2B.All functions belong to the Hilbert spaceL2B=f:∫Bfβ⟶2dβ⟶<∞, which is endowed with the norm fL2B=∫Bfβ⟶2dβ⟶1/2. r is the error between the measurement stress and the stress calculated by the numerical model. fm is the measurement stress, fiexp is the calculated stress. l⟶ is unknown coefficients of the boundary conditions. n is the number of the measurement points. i is the stress components.


## 3. POD Method

Let us consider a set of functionsukk=1m∈L2Ω. u¯ is the average of these functions expressed as u¯=ukx⟶, and the symbol · denotes average. This set of functions can be considered as snapshots of the solution of a field problem at different conditions or time domains. The POD method is to use a group of data sets in n-dimensional spaceuk, which is called snapshot, finding an m-dimensional (m << n) orthogonal subspaceVm=spanΦii=1m and making the mapping error from uk to Φ minimum for any given number of basis function over all possible sets of basic functions such that m≤M [20, 21], namely,(2)Minimizeu¯−Pu¯,ΦL2Ω2=1For simplify, a constraint is applied as follows, whereL2Ω is Hilbert space, ·2 is the L2-norm, P is a projection operator: P:uk⟶Φ, that is,(3)Pu¯=u¯,ΦΦ,ΦΦ,and therefore,(4)u¯−Pu¯L2Ω2=u¯−u¯,ϕϕ,ϕϕ,u¯−u¯,ϕϕ,ϕϕ=u¯2−u¯,ϕ2ϕ,ϕ=u¯2−u¯,ϕ2.Applying a Lagrange multiplier to the minimization problem described above, the equation is transformed into the Lagrange equation, whereλ is Lagrange division.(5)Minimizeu¯2−u¯,Φ2+λΦ2−1⟹Maximizeu¯,Φ2−λΦ2−1.The minimization problem can be turned into a maximization problem. If(6)JΦ=u¯,Φ2−λΦ2−1.Setting the Gateaux derivative to zero, the extreme-value problem turned into the following:(7)ddεJϕ+εψ|ε=0=0⟹ddεu¯,ϕ+εψ2−λϕ+εψ2−1=ddεu¯,ϕ+εu¯,ψ2−λϕ2+2εϕ,ψ+ε2ϕ2−1=ddεu¯,ϕ2+2εu¯,ϕu¯,ψ+ε2u¯,ψ2−λϕ2+2εϕ,ψ+ε2ϕ2−1=2u¯,ϕu¯,ψ+2εu¯,ψ2−λ2ϕ,ψ+2εϕ2=u¯,ϕu¯,ψ−λϕ,ψ=∫uxϕxdx∫uξψξdx−λ∫ϕξψξdx=∫ψξdx∫uxuξϕxdx−λϕξ=0⇒∫uxuξϕxdx=λϕε.Even though the above eigenvalue problem is also computationally expensive as the problem of sizeN×N, the snapshot method proposed by Sirovichd can be used to reduce the problem to a smaller M×M sized problem, that is, getting a set of the orthogonal basis of Vm=spanΦii=1m.


### 3.1. Selection of Basis Functions

POD is used to select an optimalm-dimensional orthogonal subspace Vm=spanϕ1,…,ϕl from a group of snapshots. Each vector ϕii=1,…,min space is called a POD mode. The criterion of selecting the POD basis is the dimension of the selected subspace, which is as small as possible, but containing the maximum information of the data containing maximal energy. The mathematic expression is as follows:(8)Im=∑i=1mλi∑i=1lλi,subjecting to(9)Im≥γ≈1,γ is the percentage of the vector space Vm containing the information in vector space V. As a result, the state variables of the model can be expressed by a linear combination of the POD basis.(10)Vt,x=V¯+∑i=1mαitϕix,αiti=1,…,m represents the corresponding coefficients of each POD basis.


### 3.2. Finite Element Method Based on POD

During the long-term geological evolution, the in situ stress field is constantly influenced by the internal and external environmental effects of the Earth’s crust, which makes the in situ stress field very complex on the temporal and spatial scales. But for a geological age, the in situ stress field can be regarded as a relatively stable field without considering the time factor. A lot of research and engineering practices show that the Earth stress is a function of the following variables [22]:(11)σ=Fx,y,z,E,μ,ΔG,ΔH,T,W,Q.In which,σ is the in situ stress field; x,y,z is the coordinates of the spatial position of the measuring point; E,μ is the elastic modulus and Poisson's ratio of the rock mass; ΔG is gravity factor; ΔH is the geological tectonic factor; T,W is the temperature field and the seepage field; Q is the grey factor, which stands for the other undetermined factors affecting the regional in situ stress field.In the inverse calculation of in situ stress field,x,y,z can be measured in the engineering site; E,μ can be obtained by material test, which can be considered as the known quantity. T,W,Q has little influence on the distribution of the in situ stress field, which can be neglected in general. ΔG, ΔH is the main unknown quantity and can be calculated by assigning the density and gravitational acceleration to the stratum, while ΔH can be modelled by applying corresponding boundary conditions to the numerical model.In general, the governing equation for the calculation of the initial in situ stress field is referred to as follows:(12)∇⋅σx¯+b⟶x⟶=0.The weak form is shown as follows:(13)∫Ωδd∇⋅σx⟶+b⟶x⟶dΩ=0,(14)⇒∫ΓeδdTτdΓe+∫ΩeδdTbdΩe=∫ΩeδεTσdΩe.Setting(15)ux⟶≈∑1mαiϕix⟶,wx⟶≈∑1mdiϕix⟶.Takingd=Nde,ε=Bde,σ=DBde into equation (14), there are the following equations:(16)∫ΩeNδdeTbdΩe=∫ΩeδdeTNTbdΩe,(17)∫ΓeNδdeTτdΓe=∫ΓeδdeTNTτdΓe(18)∫ΩeBδdeTDBdedΩe=∫ΩeδdeTBTDBdedΩe,So, it comes to the following:(19)δdeT∫ΩeNTbdΩe+∫ΓeNTτdΓe−∫ΩeBTDBdedΩe=0,(20)⟹∫ΩeNTbdΩe+∫ΓeNTτdΓe−∫ΩeBTDBdedΩe=0,(21)⟹∫ΩeBTDBdedΩe=∫ΩeNTbdΩe+∫ΓeNTτdΓe.SettingKe=∫ΩeBTDBdΩe,(22)be=∫ΩeNTbdΩe,(23)fe=∫ΓeNTτdΓe,(24)⇒Kede=fe+be,If(25)Kg=∑elementsKe,Pg=∑elementsbe+fe+F,dg=all displacements,(26)⇒Kgdg=Pg.Using POD orthogonal basis, it finally gets the following:(27)ΦTKΦa=pΦ.By using the Snapshot method, the matrix goes fromN × N to M × M. The POD method can greatly reduce the computational effort.The SMARS algorithm combining the POD reduced order numerical model can inversely calculate the in situ stress field in a very short time for any engineering site and give an accurate solution.


### 3.3. The Process of POD Method

The process of generating the POD reduced order model is shown in Figure2. Firstly, different field variables should be obtained by the full order model. Then, eigenvectors and eigenvalues are extracted, and the vector subspace is created afterwards. Finally, the reduced order model is generated using the pod mode.Figure 2
Flow chart of the POD method.


## 3.1. Selection of Basis Functions

POD is used to select an optimalm-dimensional orthogonal subspace Vm=spanϕ1,…,ϕl from a group of snapshots. Each vector ϕii=1,…,min space is called a POD mode. The criterion of selecting the POD basis is the dimension of the selected subspace, which is as small as possible, but containing the maximum information of the data containing maximal energy. The mathematic expression is as follows:(8)Im=∑i=1mλi∑i=1lλi,subjecting to(9)Im≥γ≈1,γ is the percentage of the vector space Vm containing the information in vector space V. As a result, the state variables of the model can be expressed by a linear combination of the POD basis.(10)Vt,x=V¯+∑i=1mαitϕix,αiti=1,…,m represents the corresponding coefficients of each POD basis.


## 3.2. Finite Element Method Based on POD

During the long-term geological evolution, the in situ stress field is constantly influenced by the internal and external environmental effects of the Earth’s crust, which makes the in situ stress field very complex on the temporal and spatial scales. But for a geological age, the in situ stress field can be regarded as a relatively stable field without considering the time factor. A lot of research and engineering practices show that the Earth stress is a function of the following variables [22]:(11)σ=Fx,y,z,E,μ,ΔG,ΔH,T,W,Q.In which,σ is the in situ stress field; x,y,z is the coordinates of the spatial position of the measuring point; E,μ is the elastic modulus and Poisson's ratio of the rock mass; ΔG is gravity factor; ΔH is the geological tectonic factor; T,W is the temperature field and the seepage field; Q is the grey factor, which stands for the other undetermined factors affecting the regional in situ stress field.In the inverse calculation of in situ stress field,x,y,z can be measured in the engineering site; E,μ can be obtained by material test, which can be considered as the known quantity. T,W,Q has little influence on the distribution of the in situ stress field, which can be neglected in general. ΔG, ΔH is the main unknown quantity and can be calculated by assigning the density and gravitational acceleration to the stratum, while ΔH can be modelled by applying corresponding boundary conditions to the numerical model.In general, the governing equation for the calculation of the initial in situ stress field is referred to as follows:(12)∇⋅σx¯+b⟶x⟶=0.The weak form is shown as follows:(13)∫Ωδd∇⋅σx⟶+b⟶x⟶dΩ=0,(14)⇒∫ΓeδdTτdΓe+∫ΩeδdTbdΩe=∫ΩeδεTσdΩe.Setting(15)ux⟶≈∑1mαiϕix⟶,wx⟶≈∑1mdiϕix⟶.Takingd=Nde,ε=Bde,σ=DBde into equation (14), there are the following equations:(16)∫ΩeNδdeTbdΩe=∫ΩeδdeTNTbdΩe,(17)∫ΓeNδdeTτdΓe=∫ΓeδdeTNTτdΓe(18)∫ΩeBδdeTDBdedΩe=∫ΩeδdeTBTDBdedΩe,So, it comes to the following:(19)δdeT∫ΩeNTbdΩe+∫ΓeNTτdΓe−∫ΩeBTDBdedΩe=0,(20)⟹∫ΩeNTbdΩe+∫ΓeNTτdΓe−∫ΩeBTDBdedΩe=0,(21)⟹∫ΩeBTDBdedΩe=∫ΩeNTbdΩe+∫ΓeNTτdΓe.SettingKe=∫ΩeBTDBdΩe,(22)be=∫ΩeNTbdΩe,(23)fe=∫ΓeNTτdΓe,(24)⇒Kede=fe+be,If(25)Kg=∑elementsKe,Pg=∑elementsbe+fe+F,dg=all displacements,(26)⇒Kgdg=Pg.Using POD orthogonal basis, it finally gets the following:(27)ΦTKΦa=pΦ.By using the Snapshot method, the matrix goes fromN × N to M × M. The POD method can greatly reduce the computational effort.The SMARS algorithm combining the POD reduced order numerical model can inversely calculate the in situ stress field in a very short time for any engineering site and give an accurate solution.


## 3.3. The Process of POD Method

The process of generating the POD reduced order model is shown in Figure2. Firstly, different field variables should be obtained by the full order model. Then, eigenvectors and eigenvalues are extracted, and the vector subspace is created afterwards. Finally, the reduced order model is generated using the pod mode.Figure 2
Flow chart of the POD method.


## 4. Case Study by Employing the POD Method

POD algorithm is realized by the software FEniCS (free software for automated solution of differential equations).


### 4.1. A 2D Simple Slope

A simple 2D model, which is 500 m × 400 m, is created to estimate the POD algorithm. The model is simply shown as in Figure3. The gravity and the horizontal pressure of 3 MPa are applied to the model. The elastic modulus is 28 GPa, Poisson’s ratio is 0.23, and density is 2800 kg/m3.Figure 3
Schematic diagram of the calculation model.


#### 4.1.1. Accuracy Test of Reduced Order Model

The mesh of the model is 500 × 400, and each small rectangle is divided into two triangular meshes. Firstly, the bottom boundary in the vertical direction and the left boundary in the horizontal direction are fixed. A gravity load and a uniform pressure 3 MPa on the right boundary are applied to the model. The model is calculated by the full order model and reduced model, respectively. The results are shown in Table1 of both the full order model and reduced order model.Table 1
The displacement in given node got by reduced order model and full order model.
Node numberHorizontal displacementsVertical displacementsAbaqus (m)FEniCS (m)ErrorAbaqus (m)FEniCS (m)Error50233−4.1500 × 10−3−4.1452 × 10−30.1157%−3.0400 × 10−2−3.0428 × 10−20.0921%52449−9.8100 × 10−3−9.8088 × 10−30.0122%−3.0700 × 10−2−3.0692 × 10−20.0261%103867−6.8000 × 10−3−6.8035 × 10−30.0515%−5.2800 × 10−2−5.2831 × 10−20.0587%107093−1.6800 × 10−2−1.6768 × 10−20.1904%−5.1100 × 10−2−5.1134 × 10−20.0665%155967−8.6300 × 10−3−8.6329 × 10−30.0336%−6.4300 × 10−2−6.4272 × 10−20.0435%164215−2.4500 × 10−2−2.4516 × 10−20.0653%−6.1100 × 10−2−6.1136 × 10−20.0589%200542−8.7200 × 10−3−8.7153 × 10−30.0062%−6.6100 × 10−2−6.6118 × 10−20.0272%200901−3.9800 × 10−2−3.9811 × 10−20.0276%−5.4000 × 10−2−5.4042 × 10−20.0778%It can be seen from the table that the results obtained by the full order model using the commercial software ABAQUS and by reduced order model based on FEniCS seem very close. The relative error is small. Therefore, it can be proved the reduced order model could give a correct solution. What is more, it cost 44.4 s by full order model while it costs 3.9 s by reduced order model, which is 8.78% of the full order model, greatly improving the computational efficiency.


#### 4.1.2. Comparison of Full Order Model and Reduced Order Model

Then, a simple assumed 2D in situ stress field is generated and calculated using the full order model and reduced order model separately combining the SMARS algorithm. The time consuming for full order model is 9936.4557 s, and the coefficients for boundary loads are 2.9466 for uniform pressure, 0.1555 for triangular pressure, −0.0273 for shear load, 10.0182 for gravity; for reduced order model, it cost 1366.9811 s, 3.0572 for uniform pressure, 0.081625 for triangular pressure, −0.028068for shear load, −9.9622for gravity. The horizontal displacement, vertical displacement, and stress distribution of the two methods are plotted, respectively.It is apparent in Figures4–6 and Table 2 that the distribution of the stress field is basically the same as the optimization results using the two methods. The maximum and minimum displacements are very close. The errors of the calculated horizontal stress and vertical stress with the measured stress are small, but the error of shear stress is relatively large. In a word, the precision using the reduced order model to do the optimization can be guaranteed.Figure 4
The distribution nephogram of horizontal displacement (m).
(a)(b)Figure 5
The distribution nephogram of vertical displacement (m).
(a)(b)Figure 6
The distribution nephogram of Mises stress (Pa).
(a)(b)Table 2
The sample stress got by reduced order model and full order model.
Element numberHorizontal stressVertical stressShear stressSample stress (Pa)Error of PODError of FOMSample stress (Pa)Error of PODError of FOMSample stress(Pa)Error of PODError of FOM50233−3.83 × 1060.522%0.783%−1.01 × 10700−3.21 × 10600.312%52449−3.65 × 1060.274%0.822%−1.00 × 1070.200%0−3.14 × 10600.318%103867−3.01 × 10601.66%−7.72 × 1060.130%0.130%−2.47 × 10600.405%107093−3.50 × 1060.571%0.857%−8.51 × 1060.470%0−2.76 × 10600.362%155967−2.99 × 1060.334%1.67%−6.09 × 1060.164%0.164%−2.09 × 10600.478%164215−3.14 × 1060.637%0.955%−6.84 × 1060.439%0−2.30 × 1060.435%0200542−2.93 × 1060.683%1.02%−5.68 × 1060.352%0−1.98 × 10600.505%200901−2.98 × 1060.671%0.671%−5.17 × 1060.193%0.193%−1.87 × 1060.535%0.535%


### 4.2. A Slope with Multiple Stratums

In order to apply the POD method to the general geotechnical engineering problem, a slope is created with three stratums as shown in Figure7, the mechanical parameters of the soil layer are as shown in Table 3.Figure 7
A slope with multiple stratums (unit: m).Table 3
The mechanical parameters of the soil layer.
Material numberDensity (kg/m3)Young’s modulus (GPa)Poisson’s ratio1253017.40.142264041.30.133241028.60.17In the traditional inverse calculation of in situ stress field, the distribution of the boundary condition which is related to the tectonic stress field should be assumed at first, because the distribution of the boundary load has the corresponding relation with the stress field, that is when the boundary load is applied to the area to be inspected, all the stress distribution information of the regional stress field can be obtained. However, it is generally assumed that the boundary condition obeys a simple linear distribution. For a region with complex geological conditions, the inverse calculation of the in situ stress field cannot be simply assumed that the boundary conditions are linearly distributed or conform to a certain form so that a multinomial function distribution of the boundary load can be constructed to create an arbitrary complex in situ stress field distribution.The boundary load expressed by polynomials is of great value for practical engineering. It is proved that the calculation accuracy of the in situ stress field by the uniform distribution of boundary load can meet the needs of engineering design and construction for simple geological conditions. Therefore, a three-degree polynomial distribution can fit the general engineering site accurately; that is, the boundary load distributions are subordinated to the function as follows:(28)y=ax3+bx2+cx+d.ABAQUS is used to realize applying the boundary conditions to the numerical model, and six points are selected as “measured points,” and the stress components are exported in Table4.Table 4
“The measured data” of the in situ stress.
Element numberSxx (Pa)Syy (Pa)Sxy (Pa)50136−2.54E + 06−3.09E + 063.49E + 0551467−1.50E + 06−2.12E + 06−1.60E + 058221−3.57E + 06−8.14E + 06−5.55E + 057250−2.14E + 06−3.33E + 06−5.43E + 0529248−2.85E + 06−7.25E + 063.37E + 0427870−2.49E + 06−5.45E + 062.74E + 04Table5 is the in situ stress in the “measured points” using a full order model with ABAQUS combining a surrogate model accelerated random search algorithm. Table 6 is the in situ stress in the “measured points” using a reduced order model. Based on the data in Tables 5 and 6, the calculation error by full order model and reduced order model is analysed and a histogram is plotted in Figure 8. No matter whether the full order model or the reduced order model is used to do the inverse calculation, the calculated stress is very close to “the measured stress.” At some individual measuring points, the stress error in the vertical direction is zero, and the error of all points is less than 1.25%.Table 5
The calculated stress in measured points by full order model.
Element numberSxx (Pa)ErrorSyy (Pa)ErrorSxy (Pa)Error50136−2.54E + 060%−3.09E + 060%3.49E + 050.00%51467−1.50E + 060%−2.12E + 060%−1.60E + 050.00%8221−3.57E + 060%−8.14E + 060%−5.54E + 050.18%7250−2.14E + 060%−3.33E + 060%−5.42E + 050.18%29248−2.85E + 060%−7.26E + 060.14%3.39E + 040.59%27870−2.49E + 060%−5.45E + 060%2.76E + 040.73%Table 6
The calculated stress in measured points by the reduced order model.
Element numberSxx (Pa)ErrorSyy (Pa)ErrorSxy (Pa)Error50136−2.56E + 060.79%−3.09E + 060.00%3.48E + 050.29%51467−1.51E + 060.67%−2.12E + 060.00%−1.62E + 051.25%8221−3.60E + 060.84%−8.14E + 060.00%−5.56E + 050.36%7250−2.16E + 060.93%−3.33E + 060.00%−5.45E + 050.55%29248−2.87E + 060.70%−7.25E + 060.14%3.39E + 040.00%27870−2.51E + 060.80%−5.44E + 060.18%2.73E + 041.09%Figure 8
The comparison between the full order model and the reduced order model.
(a)(b)(c)(d)In addition, the inverse calculation of the in situ stress field by full order model using ABAQUS is 9554.386 s, while that of the reduced order model is 931.173 s, which is 9.746% of the time of the full order model. The inverse calculation of the in situ stress field uses a reduced order model method greatly improving the efficiency of calculation on the premise of high accuracy.


## 4.1. A 2D Simple Slope

A simple 2D model, which is 500 m × 400 m, is created to estimate the POD algorithm. The model is simply shown as in Figure3. The gravity and the horizontal pressure of 3 MPa are applied to the model. The elastic modulus is 28 GPa, Poisson’s ratio is 0.23, and density is 2800 kg/m3.Figure 3
Schematic diagram of the calculation model.


### 4.1.1. Accuracy Test of Reduced Order Model

The mesh of the model is 500 × 400, and each small rectangle is divided into two triangular meshes. Firstly, the bottom boundary in the vertical direction and the left boundary in the horizontal direction are fixed. A gravity load and a uniform pressure 3 MPa on the right boundary are applied to the model. The model is calculated by the full order model and reduced model, respectively. The results are shown in Table1 of both the full order model and reduced order model.Table 1
The displacement in given node got by reduced order model and full order model.
Node numberHorizontal displacementsVertical displacementsAbaqus (m)FEniCS (m)ErrorAbaqus (m)FEniCS (m)Error50233−4.1500 × 10−3−4.1452 × 10−30.1157%−3.0400 × 10−2−3.0428 × 10−20.0921%52449−9.8100 × 10−3−9.8088 × 10−30.0122%−3.0700 × 10−2−3.0692 × 10−20.0261%103867−6.8000 × 10−3−6.8035 × 10−30.0515%−5.2800 × 10−2−5.2831 × 10−20.0587%107093−1.6800 × 10−2−1.6768 × 10−20.1904%−5.1100 × 10−2−5.1134 × 10−20.0665%155967−8.6300 × 10−3−8.6329 × 10−30.0336%−6.4300 × 10−2−6.4272 × 10−20.0435%164215−2.4500 × 10−2−2.4516 × 10−20.0653%−6.1100 × 10−2−6.1136 × 10−20.0589%200542−8.7200 × 10−3−8.7153 × 10−30.0062%−6.6100 × 10−2−6.6118 × 10−20.0272%200901−3.9800 × 10−2−3.9811 × 10−20.0276%−5.4000 × 10−2−5.4042 × 10−20.0778%It can be seen from the table that the results obtained by the full order model using the commercial software ABAQUS and by reduced order model based on FEniCS seem very close. The relative error is small. Therefore, it can be proved the reduced order model could give a correct solution. What is more, it cost 44.4 s by full order model while it costs 3.9 s by reduced order model, which is 8.78% of the full order model, greatly improving the computational efficiency.


### 4.1.2. Comparison of Full Order Model and Reduced Order Model

Then, a simple assumed 2D in situ stress field is generated and calculated using the full order model and reduced order model separately combining the SMARS algorithm. The time consuming for full order model is 9936.4557 s, and the coefficients for boundary loads are 2.9466 for uniform pressure, 0.1555 for triangular pressure, −0.0273 for shear load, 10.0182 for gravity; for reduced order model, it cost 1366.9811 s, 3.0572 for uniform pressure, 0.081625 for triangular pressure, −0.028068for shear load, −9.9622for gravity. The horizontal displacement, vertical displacement, and stress distribution of the two methods are plotted, respectively.It is apparent in Figures4–6 and Table 2 that the distribution of the stress field is basically the same as the optimization results using the two methods. The maximum and minimum displacements are very close. The errors of the calculated horizontal stress and vertical stress with the measured stress are small, but the error of shear stress is relatively large. In a word, the precision using the reduced order model to do the optimization can be guaranteed.Figure 4
The distribution nephogram of horizontal displacement (m).
(a)(b)Figure 5
The distribution nephogram of vertical displacement (m).
(a)(b)Figure 6
The distribution nephogram of Mises stress (Pa).
(a)(b)Table 2
The sample stress got by reduced order model and full order model.
Element numberHorizontal stressVertical stressShear stressSample stress (Pa)Error of PODError of FOMSample stress (Pa)Error of PODError of FOMSample stress(Pa)Error of PODError of FOM50233−3.83 × 1060.522%0.783%−1.01 × 10700−3.21 × 10600.312%52449−3.65 × 1060.274%0.822%−1.00 × 1070.200%0−3.14 × 10600.318%103867−3.01 × 10601.66%−7.72 × 1060.130%0.130%−2.47 × 10600.405%107093−3.50 × 1060.571%0.857%−8.51 × 1060.470%0−2.76 × 10600.362%155967−2.99 × 1060.334%1.67%−6.09 × 1060.164%0.164%−2.09 × 10600.478%164215−3.14 × 1060.637%0.955%−6.84 × 1060.439%0−2.30 × 1060.435%0200542−2.93 × 1060.683%1.02%−5.68 × 1060.352%0−1.98 × 10600.505%200901−2.98 × 1060.671%0.671%−5.17 × 1060.193%0.193%−1.87 × 1060.535%0.535%


## 4.1.1. Accuracy Test of Reduced Order Model

The mesh of the model is 500 × 400, and each small rectangle is divided into two triangular meshes. Firstly, the bottom boundary in the vertical direction and the left boundary in the horizontal direction are fixed. A gravity load and a uniform pressure 3 MPa on the right boundary are applied to the model. The model is calculated by the full order model and reduced model, respectively. The results are shown in Table1 of both the full order model and reduced order model.Table 1
The displacement in given node got by reduced order model and full order model.
Node numberHorizontal displacementsVertical displacementsAbaqus (m)FEniCS (m)ErrorAbaqus (m)FEniCS (m)Error50233−4.1500 × 10−3−4.1452 × 10−30.1157%−3.0400 × 10−2−3.0428 × 10−20.0921%52449−9.8100 × 10−3−9.8088 × 10−30.0122%−3.0700 × 10−2−3.0692 × 10−20.0261%103867−6.8000 × 10−3−6.8035 × 10−30.0515%−5.2800 × 10−2−5.2831 × 10−20.0587%107093−1.6800 × 10−2−1.6768 × 10−20.1904%−5.1100 × 10−2−5.1134 × 10−20.0665%155967−8.6300 × 10−3−8.6329 × 10−30.0336%−6.4300 × 10−2−6.4272 × 10−20.0435%164215−2.4500 × 10−2−2.4516 × 10−20.0653%−6.1100 × 10−2−6.1136 × 10−20.0589%200542−8.7200 × 10−3−8.7153 × 10−30.0062%−6.6100 × 10−2−6.6118 × 10−20.0272%200901−3.9800 × 10−2−3.9811 × 10−20.0276%−5.4000 × 10−2−5.4042 × 10−20.0778%It can be seen from the table that the results obtained by the full order model using the commercial software ABAQUS and by reduced order model based on FEniCS seem very close. The relative error is small. Therefore, it can be proved the reduced order model could give a correct solution. What is more, it cost 44.4 s by full order model while it costs 3.9 s by reduced order model, which is 8.78% of the full order model, greatly improving the computational efficiency.


## 4.1.2. Comparison of Full Order Model and Reduced Order Model

Then, a simple assumed 2D in situ stress field is generated and calculated using the full order model and reduced order model separately combining the SMARS algorithm. The time consuming for full order model is 9936.4557 s, and the coefficients for boundary loads are 2.9466 for uniform pressure, 0.1555 for triangular pressure, −0.0273 for shear load, 10.0182 for gravity; for reduced order model, it cost 1366.9811 s, 3.0572 for uniform pressure, 0.081625 for triangular pressure, −0.028068for shear load, −9.9622for gravity. The horizontal displacement, vertical displacement, and stress distribution of the two methods are plotted, respectively.It is apparent in Figures4–6 and Table 2 that the distribution of the stress field is basically the same as the optimization results using the two methods. The maximum and minimum displacements are very close. The errors of the calculated horizontal stress and vertical stress with the measured stress are small, but the error of shear stress is relatively large. In a word, the precision using the reduced order model to do the optimization can be guaranteed.Figure 4
The distribution nephogram of horizontal displacement (m).
(a)(b)Figure 5
The distribution nephogram of vertical displacement (m).
(a)(b)Figure 6
The distribution nephogram of Mises stress (Pa).
(a)(b)Table 2
The sample stress got by reduced order model and full order model.
Element numberHorizontal stressVertical stressShear stressSample stress (Pa)Error of PODError of FOMSample stress (Pa)Error of PODError of FOMSample stress(Pa)Error of PODError of FOM50233−3.83 × 1060.522%0.783%−1.01 × 10700−3.21 × 10600.312%52449−3.65 × 1060.274%0.822%−1.00 × 1070.200%0−3.14 × 10600.318%103867−3.01 × 10601.66%−7.72 × 1060.130%0.130%−2.47 × 10600.405%107093−3.50 × 1060.571%0.857%−8.51 × 1060.470%0−2.76 × 10600.362%155967−2.99 × 1060.334%1.67%−6.09 × 1060.164%0.164%−2.09 × 10600.478%164215−3.14 × 1060.637%0.955%−6.84 × 1060.439%0−2.30 × 1060.435%0200542−2.93 × 1060.683%1.02%−5.68 × 1060.352%0−1.98 × 10600.505%200901−2.98 × 1060.671%0.671%−5.17 × 1060.193%0.193%−1.87 × 1060.535%0.535%


## 4.2. A Slope with Multiple Stratums

In order to apply the POD method to the general geotechnical engineering problem, a slope is created with three stratums as shown in Figure7, the mechanical parameters of the soil layer are as shown in Table 3.Figure 7
A slope with multiple stratums (unit: m).Table 3
The mechanical parameters of the soil layer.
Material numberDensity (kg/m3)Young’s modulus (GPa)Poisson’s ratio1253017.40.142264041.30.133241028.60.17In the traditional inverse calculation of in situ stress field, the distribution of the boundary condition which is related to the tectonic stress field should be assumed at first, because the distribution of the boundary load has the corresponding relation with the stress field, that is when the boundary load is applied to the area to be inspected, all the stress distribution information of the regional stress field can be obtained. However, it is generally assumed that the boundary condition obeys a simple linear distribution. For a region with complex geological conditions, the inverse calculation of the in situ stress field cannot be simply assumed that the boundary conditions are linearly distributed or conform to a certain form so that a multinomial function distribution of the boundary load can be constructed to create an arbitrary complex in situ stress field distribution.The boundary load expressed by polynomials is of great value for practical engineering. It is proved that the calculation accuracy of the in situ stress field by the uniform distribution of boundary load can meet the needs of engineering design and construction for simple geological conditions. Therefore, a three-degree polynomial distribution can fit the general engineering site accurately; that is, the boundary load distributions are subordinated to the function as follows:(28)y=ax3+bx2+cx+d.ABAQUS is used to realize applying the boundary conditions to the numerical model, and six points are selected as “measured points,” and the stress components are exported in Table4.Table 4
“The measured data” of the in situ stress.
Element numberSxx (Pa)Syy (Pa)Sxy (Pa)50136−2.54E + 06−3.09E + 063.49E + 0551467−1.50E + 06−2.12E + 06−1.60E + 058221−3.57E + 06−8.14E + 06−5.55E + 057250−2.14E + 06−3.33E + 06−5.43E + 0529248−2.85E + 06−7.25E + 063.37E + 0427870−2.49E + 06−5.45E + 062.74E + 04Table5 is the in situ stress in the “measured points” using a full order model with ABAQUS combining a surrogate model accelerated random search algorithm. Table 6 is the in situ stress in the “measured points” using a reduced order model. Based on the data in Tables 5 and 6, the calculation error by full order model and reduced order model is analysed and a histogram is plotted in Figure 8. No matter whether the full order model or the reduced order model is used to do the inverse calculation, the calculated stress is very close to “the measured stress.” At some individual measuring points, the stress error in the vertical direction is zero, and the error of all points is less than 1.25%.Table 5
The calculated stress in measured points by full order model.
Element numberSxx (Pa)ErrorSyy (Pa)ErrorSxy (Pa)Error50136−2.54E + 060%−3.09E + 060%3.49E + 050.00%51467−1.50E + 060%−2.12E + 060%−1.60E + 050.00%8221−3.57E + 060%−8.14E + 060%−5.54E + 050.18%7250−2.14E + 060%−3.33E + 060%−5.42E + 050.18%29248−2.85E + 060%−7.26E + 060.14%3.39E + 040.59%27870−2.49E + 060%−5.45E + 060%2.76E + 040.73%Table 6
The calculated stress in measured points by the reduced order model.
Element numberSxx (Pa)ErrorSyy (Pa)ErrorSxy (Pa)Error50136−2.56E + 060.79%−3.09E + 060.00%3.48E + 050.29%51467−1.51E + 060.67%−2.12E + 060.00%−1.62E + 051.25%8221−3.60E + 060.84%−8.14E + 060.00%−5.56E + 050.36%7250−2.16E + 060.93%−3.33E + 060.00%−5.45E + 050.55%29248−2.87E + 060.70%−7.25E + 060.14%3.39E + 040.00%27870−2.51E + 060.80%−5.44E + 060.18%2.73E + 041.09%Figure 8
The comparison between the full order model and the reduced order model.
(a)(b)(c)(d)In addition, the inverse calculation of the in situ stress field by full order model using ABAQUS is 9554.386 s, while that of the reduced order model is 931.173 s, which is 9.746% of the time of the full order model. The inverse calculation of the in situ stress field uses a reduced order model method greatly improving the efficiency of calculation on the premise of high accuracy.


## 5. Conclusions

The proper orthogonal decomposition method is used to reduce the dimension of the numerical model, and the optimization is achieved in the sense of capturing the maximum energy of the system so as to reduce the amount of calculation and save the calculation time and the effect of CPU load.(1)
The paper describes the basic principle of the proper orthogonal decomposition method. It is derived from a mathematical point of view for solving the basis function. Then, the selection criteria of the basis function are introduced. Firstly, all eigenvalues and eigenvectors should be reordered once they are calculated. Secondly, several eigenvalues and corresponding eigenvectors are kept on condition that the number of eigenvectors is minimum but contains most of the system energy.(2)
Two examples are carried out to make the comparison. The reduced order model is treated by FEniCS which is a software to solve differential equations, while the full order model is using ABAQUS. The two inverse calculations are done combining a SMARS strategy. For the first example, it turned out the relative error range between the full order model and the reduced order model is 0.0122%–0.1904%, but time consuming is just 8.78% of the full order model by ABAQUS. For the second example, the inverse calculation results by either the full order model or the reduced order model are very close to the true solution. Most errors by full order model are 0. In addition, the vertical stress of four points is exactly the same as the “measured stress” by the reduced order model, and the maximal error is 1.25% which is relatively small. The time consumption is just 9.746% of the full order model.(3)
It can be seen that the calculation accuracy and efficiency of the reduced order model are very high. When the freedom of the model is higher, the advantage of the reduced order model will be more obvious.



---

*Source: 5539898-2021-03-01.xml*