**An Improved Method for Predicting the Greenfield Stratum Movements Caused by Shield Tunnel Construction**

#### **XinRong Tan 1, Heng Zhang 2,\*, Gang Zhang 2, Yimo Zhu 2,\* and Peng Tu <sup>3</sup>**


Received: 23 September 2019; Accepted: 21 October 2019; Published: 24 October 2019

**Abstract:** Shield tunneling is becoming the preferred construction scheme for metro construction because of its advantages of fast construction speed and small disturbance. However, limited by process defects, the stratum movements induced by the construction of shield tunnels still affects the safety of nearby underground structures and aboveground buildings. Therefore, the reliable prediction of stratum movements is important. Described in this paper is an analysis method of the Greenfield stratum movements (Greenfield is an area of land that has not yet had buildings on it, stratum movements means the movement of various soil layers) caused by shield tunnel construction combining an elastic half-space model of mirror source–sink method with the use of modified analytical method. Based on the theoretical formula in this paper, not only can the curve of surface settlement trough be calculated, but also the three-dimensional displacement field of deep soil can be obtained. By comparing vertical and horizontal contour maps of Greenfield stratum movements, good consistency between theoretical formula results and centrifugal test results are shown. This solves the defects and limitations of existing two-dimensional formulas; furthermore, based on this, it is convenient to evaluate the effect on the other skewed underground structures through the elastic foundation beam and other similar methods; therefore, this paper can provide a wide guidance and service for the design and construction of underground engineering in the future.

**Keywords:** shield tunnel; stratum movements; analytical; mirror source–sink method; centrifuge modelling test

#### **1. Introduction**

Since the 21st century, tunnel and underground engineering in China have made great progress [1–7]. With the continuous development of urbanization in China, the urban population is growing, and urban traffic pressure is also increasing. In order to alleviate the urban traffic pressure, city managers have taken a lot of measures, such as limited number driving, expansion, and transformation of existing roads, but the traffic demand is still far greater than the current urban ground carrying capacity. Metro has the advantages of large passenger capacity, fast speed, and full use of underground space, which is more and more popular in modern society. However, the urban subway has strict requirements for its construction methods. Because of the high density of buildings in cities, the construction of the metro will cause stratum movement and bring adverse effects on ground buildings. Especially when passing through high-rise buildings, the disturbance caused by construction on

stratum must be more strictly controlled. The construction period of Metro is long, and the interruption of traffic should be avoided as far as possible when crossing the road. Usually, subway construction methods include the open-cut method, shallow-buried excavation method, and shield tunneling method. Among them, the shield tunneling method has become the preferred construction scheme for metro construction because of its great advantages of fast construction speed, a high degree of automation, and small disturbance to stratum [8,9].

Although shield tunneling has many advantages mentioned above, and the construction technology has made great progress after many years of development, due to the defects of geological conditions and construction technology, the advance of shield tunneling will inevitably make a disturbance, change the stress state of soil, and cause stratum displacement, as shown in Figure 1. If the disturbance of stratum cannot be clearly understood, it will affect the safety of nearby underground structures and aboveground buildings. For example, the adjacent metro tunnels may cause deformation restrictions beyond the normal operation of the metro, or even lead to train derailment; it may also lead to the breakage of underground pipelines, resulting in a series of problems, such as gas leakage, interruption of urban water use, interruption of communication power system, and so on, affecting the daily life of the city [10–12]. Therefore, a reliable prediction of stratum movements is important.

**Figure 1.** Sources of stratum movements of shield tunnel.

Over the years, previous researchers have done numerous studies on the effects of shield tunnel construction on stratum displacement [13–18]. In the direction of theoretical research, Peck assumed that the curve of land settlement trough satisfies the conditions of Gauss distribution and invariant stratum volume, and deduced Peck's empirical formula based on many field monitoring data [19]. Attewell et al. used Peck's empirical formula to simulate the free displacement of soil at the location of existing tunnels [20–24], and then deduced the displacement of existing tunnels with elastic foundation beams and other models. However, Celestino and Klar [25,26] gradually found that the Peck curve could not accurately describe soil settlement trough in many cases. Voster, O'Reilly et al. [27–30] have proposed fitting curves from different research angles to fit the free displacement field of soils. Litwinszyn [31] put forward the theory of random medium through a sand box test. Liu, Yang, et al. [32–34] introduced this theory into the prediction of stratum displacement caused by the excavation of geotechnical tunnels. Although these studies have the advantages of concise calculation formula and convenient use, they lack a clear theoretical basis of mechanics. They can only be called a mathematical empirical method, not an analytical solution of soil settlement. Because the empirical formula is fitted by a large number of field monitoring data, the parameters of the fitting function are closely related to the actual construction conditions, and lack of a clear theoretical basis, the accuracy of the results predicted by this method is often difficult to meet the requirements, and its

defects are evident. Comparatively, the solving process of the analytic method and modified analytic method is relatively complex, but through more rigorous analysis, as long as the parameters are selected properly, it has better universality. Sagaseta [35] deduced the expression of three-dimensional displacement of stratum by using the mirror source–sink method. This formula can not only calculate the curve of surface settlement trough, but also calculate the displacement and stress of deep soil, so it is widely used [36–43]. However, the movement pattern of actual volume loss is not consistent with the assumed equivalent radial movement pattern in this method and has some deviation. Based on field monitoring data, Loganathan, Jiang, and other researchers [44–51] modified the stratum displacement pattern by exponential function and deduced the prediction formula of surface subsidence under non-equivalent radial movement pattern. However, most of these methods are only suitable for three-dimensional orthogonal plane or two-dimensional cases, which have great defects and need to be improved urgently. Many other researchers concerned the field description, soil classification for engineering purposes [52–54], and the analysis methods of data such as neural networks [55–58]. Saplachidi et al. [59] used measured settlements from the tunneling excavation for the extension of Line 3 of the Athens Metro for the verification and calibration of the empirical formulae. Suwansawat et al. [60] evaluated the potential and the limitations of artificial neural networks (ANN) for predicting surface settlements caused by Earth pressure balance (EPB) shield tunneling and to develop optimal neural network models for this objective. Ahangari et al. [61] created a database from previous research [62–65] and studied the capability of adaptive neuro-fuzzy inference system (ANFIS) and gene expression programming (GEP) methods for settlement prediction. However, those methods belong to the fuzzy solutions without mechanism analysis, training neural networks requires a lot of data from databases or numerical simulation, and there are problems such as difficult convergence or slow convergence speed.

This paper studies the Greenfield stratum movements caused by shield tunnel construction combining an elastic half-space model of mirror source–sink method with the use of modified analytical method. Through theoretical derivation and empirical function modification, the prediction formula of three-dimensional Greenfield stratum movements is improved, and it shows good consistency with the centrifugal test results. The new solution can calculate the three-dimensional stratum displacement accurately, which is unavailable for the previous methods. It can be used in the process of tunnel excavation rather than only for the ultimate settlement. Based on this improved solution, it is convenient to evaluate the effect on the other skewed underground structures through the elastic foundation beam and other similar methods; therefore, this paper can provide a wide guidance and service for the design and construction of underground engineering in the future.

#### **2. Basic Theory and Calculation Method**

#### *2.1. Disturbance Factors*

Usually, there are three factors that cause disturbance of surrounding soil by shield tunneling technology: First, the gap between segment lining and stratum due to the reasons of assembling segment lining, shield overcutting, shield snake-shaped, and shield body adjustment; although synchronous grouting is adopted, the grouting cannot completely eliminate the gap at present [66]; the second is the elastic-plastic deformation of the soil at the shield heading face, which is closely related to the soil characteristics and the additional force on the shield working face. Thirdly is the horizontal friction between the outer surface of the shield and the surrounding soil. Compared with the effect of the gap, the stratum movements caused by the additional force on the shield working face and the frictional force between the outer surface of the shield and the surrounding soil during tunnel construction is very small [67]. Therefore, this paper neglects the effects of the latter two factors when analyzing, and mainly studies the effects of the gap. Sources of stratum movements of shield tunnel are shown in Figure 1 and the coordinate system is set as shown in Figure 2 [68].

**Figure 2.** Surface displacement caused by shield tunnel.

#### *2.2. Gap Parameter*

In order to analyze the effects of the above factors, Rowe et al. [69,70] proposed the concept of gap parameter to reflect the stratum loss caused by shield tunnel construction, transforming the three-dimensional problem to a two-dimensional gap in plane. The composition of the gap takes into account the following factors: Elastic-plastic deformation of the shield at the tunnel excavation face, the over-excavation of the shield machine, the physical gap between the outer skin of shield shell and segment lining, and the construction technology. The formula for calculating the gap parameters *g* under undrained condition is as in Equation (1).

$$\mathbf{g} = \mathbf{G}\_P + \mathbf{L}\_{3D}^\* + \boldsymbol{\omega} \tag{1}$$

where shows that the gap parameters are composed of three parts [44]:

The first part *Gp* is the physical gap, which represents the gap between the outer skin of shield shell and segment lining; its expression is as in Equation (2):

$$G\_P = 2\Lambda + \delta \tag{2}$$

where Δ is the thickness of shield shell and δ is the space formed by segment lining installation.

The second part *U3D* is the elastic-plastic deformation of shield at the tunnel excavation face; its quantitative expressions are as in Equations (3) and (4):

$$
\Delta l\_{3D} = \frac{k}{2} \delta\_y = \frac{k}{2} \frac{\Omega RP\_0}{E} \tag{3}
$$

$$P0 = K\_0'P\_\upsilon' + P\_w - P\_f \tag{4}$$

where *k* is the soil cutting resistance coefficient; Ω is the displacement coefficient; *R* is the tunnel radius; *E* is the elastic modulus; *K'0* is the effective static lateral pressure coefficient; *P'v* is the effective vertical stress at depth of tunnel center; *Pw* is the pore water pressure at depth of tunnel center; and *Pf* is the supporting force provided by tunnel.

*Appl. Sci.* **2019**, *9*, 4522

The third part ω is the shield construction technology coefficient, of which the value is as in Equation (5):

$$\begin{cases} \frac{lI\_i}{R} = 1 - \frac{1}{\left(1 + \frac{2(1 + \mu\_u)C\_u}{E\_u}\right) \exp(\frac{N-1}{2})} \dot{\frac{1}{2}}\\\ N = \frac{\mathcal{V}H - p\_f}{C\_u} \\\ \omega = \min\{0.6G\_{p\_f}, \frac{1}{3}\mathcal{U}\_i\} \end{cases} \tag{5}$$

where *Ui* is the elastoplastic plane strain displacement at the tunnel crown; μ*<sup>u</sup>* is the soil Poisson's ratio under undrained condition; *Cu* is the soil shear strength under undrained condition; *Eu* is the soil elastic modulus under undrained condition; γ is the unit weight of soil; and *H* is the buried depth of tunnel center.

#### *2.3. Loganathan' Solution*

Based on the above formulas, Loganathan proposed to use plane loss rate ε*x,z* to quantify stratum loss with non-equivalent radial movement. The expression is as in Equation (6):

$$\varepsilon\_{\chi z} = \frac{4gR + g^2}{4R^2} \exp\left\{-\left|\frac{1.38x^2}{\left(H + R\right)^2} + \frac{0.69z^2}{H^2}\right|\right\}.\tag{6}$$

Combined with Verruijit's [71,72] analytical formula, he put forward a formula for estimating stratum displacement caused by undrained stratum loss in clays, of which the formulas are as in Equation (7):

$$\begin{cases} S\_z = R^2 \left\{ -\frac{z - H}{\mathbf{x}^2 + (z - H)^2} + (3 - 4\mu) \frac{z + H}{\mathbf{x}^2 + (z + H)^2} - \frac{2z\left[\mathbf{x}^2 - (z + H)^2\right]}{\left[\mathbf{x}^2 + (z + H)^2\right]^2} \right\} \times \varepsilon\_{\mathbf{x}, z} \\\ S\_\mathbf{x} = -R^2 \mathbf{x} \left\{ \frac{1}{\mathbf{x}^2 + (z - H)^2} + \frac{\left(3 - 4\mu\_u\right)}{\mathbf{x}^2 + (z + H)^2} - \frac{4z(z + H)}{\left[\mathbf{x}^2 + (z + H)^2\right]^2} \right\} \times \varepsilon\_{\mathbf{x}, z} \end{cases} \tag{7}$$

Considering that the displacement solution caused by stratum loss derived by Loganathan is a two-dimensional solution, which is often a three-dimensional problem in actual engineering, it is necessary to extend the two-dimensional displacement solution to three-dimensional. In the above deduction process, we can find that the stratum loss parameter ε*x,z* is the key for extending. If we extend the plane stratum loss parameter ε*x,z* to three-dimensional parameter ε*x,z,y*, the displacement solution can be extended to three-dimensional space. Another important method to solve the stratum displacement caused by stratum loss is used here. The concept of stratum loss proposed by Sagaseta proposed the mirror source–sink method and derives the problem of free displacement field by elastic half-space model. Although field monitoring data [40] show that the movement pattern of actual volume loss is not consistent with the assumed equivalent radial movement pattern in this method, its process is a strict mathematical derivation process, so the three-dimensional stratum loss parameter ε*x,z,y* can be extracted from the Sagaseta solution.

#### *2.4. Sagaseta's Solution*

Sagaseta used the mirror source–sink method, assumed that the surface is free, the soil is isotropic, and incompressible. The analysis steps are as follows and in Figure 3:

**Figure 3.** Analysis steps of mirror source–sink method.

First step: The existence of the original surface is ignored, and the half-space problem is transformed into the internal gap problem of full space. Under the action of a source gap, stress field changing, there will be the normal stress σ*<sup>0</sup>* and the shear stress τ*<sup>0</sup>* at the original surface position.

Second step: A sink is set at the mirror position of the source gap above the original surface; there is a volume expansion of the same size at the sink position. The sink expansion will produce a normal stress −σ*<sup>0</sup>* of the same value and opposite direction and a shear stress τ*<sup>0</sup>* of the same value and same direction at the original surface position. Therefore, under the action of source and sink, the normal stress on the original surface can be offset, which does not exist in practice.

Third step: Under the action of source and sink, a shear stress −*2*τ*<sup>0</sup>* of twice the value and opposite direction is supposed at the original surface position and the corresponding displacement field of the surface shear stress below the surface at each point can be obtained.

Final step: The displacements generated by the above three steps at any point below the surface are *S1*, *S2*, and *S3*, respectively. The displacement solution of the practical problem is *S* = *S1* + *S2* + *S3*.

Based on above steps, Sagaseta deduced the expression of surface displacement caused by tunnel construction is as in Equation (8):

⎧

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

$$\begin{aligned} S\_{x0} &= -\frac{v\_{\text{loss}}}{2\pi} \frac{x}{x^2 + H^2} \left| 1 + \frac{y}{(x^2 + y^2 + H^2)\frac{1}{2}} \right| \\ S\_{y0} &= \frac{v\_{\text{loss}}}{2\pi} \frac{1}{(x^2 + y^2 + H^2)\frac{1}{2}} \\ S\_{z0} &= \frac{v\_{\text{loss}}}{2\pi} \frac{H}{x^2 + H^2} \left| 1 + \frac{y}{\frac{(x^2 + y^2 + H^2)\frac{1}{2}}} \right| \\ S\_{x0}(y \to \infty) &= -\frac{v\_{\text{loss}}}{\pi} \frac{\frac{1}{x^2 + H^2}}{x^2 + H^2} \\ S\_{z0}(y \to \infty) &= \frac{v\_{\text{loss}}}{\pi} \frac{H}{x^2 + H^2} \end{aligned} \tag{8}$$

where *Vloss* is stratum loss. Contrastive Equation (8) can be found:

$$\begin{cases} S\_{x0}(y) = S\_{x0}(y \to \infty) \times \frac{1}{2} \times \left| 1 + \cfrac{y}{(x^2 + y^2 + H^2)^{\frac{\pi}{2}}} \right| \\\\ S\_{z0}(y) = S\_{z0}(y \to \infty) \times \frac{1}{2} \times \left| 1 + \cfrac{y}{(x^2 + y^2 + H^2)^{\frac{\pi}{2}}} \right| \\\\ \end{cases} \tag{9}$$

#### *2.5. An Improved Method*

Equation (9) shows that the distribution of stratum displacement caused by shield construction

(stratum loss) along the y-axis driving satisfies the law of <sup>1</sup> <sup>2</sup> <sup>×</sup> ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ <sup>1</sup> <sup>+</sup> *<sup>y</sup>* (*x*<sup>2</sup> + *y*<sup>2</sup> + *H*2) 1 2 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ . Therefore, based

on this law, the Loganathan's solution is extended to the three-dimensional. It is noticed that the coordinate system adopted in this paper is different from that used in the Sagaseta's solution, so the distribution law should be revised, and its expression is as in Equation (10):

$$\kappa\_{\mathcal{Y}} = \frac{1}{2} \times \left| 1 - \frac{y}{\frac{(\mathbf{x}^2 + y^2 + H^2)\frac{1}{2}}{(\mathbf{x}^2 + y^2 + H^2)\frac{1}{2}}} \right|. \tag{10}$$

The plane stratum loss parameter ε*x,z* is extended to three-dimensional parameter ε*x,z,y*:

$$\varepsilon\_{x,y,x} = \varepsilon\_{y} \times \varepsilon\_{x,z} = \left[1 - \frac{y}{\frac{1}{\left(x^2 + y^2 + H^2\right)\tilde{\Delta}}}\right] \times \frac{4gR + g^2}{8R^2} \exp\left\{-\left|\frac{1.38x^2}{\left(H + R\right)^2} + \frac{0.69z^2}{H^2}\right|\right\} \tag{11}$$

ε*x,z* is replaced by ε*x,z,y* in Equation (7) and Equation (12) is the three-dimensional analytical solution of stratum displacement caused by stratum loss. When *y* = 0, *z* = 0, the curves of settlement trough calculated by Equations (7), (8), and (12) are the same, and Equation (8) is the particular solution about the surface settlement value of Equation (12) at the tunnel face.

$$\begin{cases} S\_z = R^2 \left[ -\frac{z - H}{\mathbf{x}^2 + (z - H)^2} + (3 - 4\mu) \frac{z + H}{\mathbf{x}^2 + (z + H)^2} - \frac{2z\left[\mathbf{x}^2 - (z + H)^2\right]}{\left[\mathbf{x}^2 + (z + H)^2\right]^2} \right] \times \varepsilon\_{x, y, z} \\\ S\_x = -R^2 \mathbf{x} \left\{ \frac{1}{\mathbf{x}^2 + (z - H)^2} + \frac{(3 - 4\mu)}{\mathbf{x}^2 + (z + H)^2} - \frac{4z(z + H)}{\left[\mathbf{x}^2 + (z + H)^2\right]^2} \right\} \times \varepsilon\_{x, y, z} \end{cases} \tag{12}$$

For example, let *H* = 13.65 m, μ = 0.4, *R* = 2.325 m, *g* = 0.058 m, then the displacement of soil layers at *z* = 0 m and *z* = 10 m is showed in Figure 4.

**Figure 4.** Settlement trough in different depth stratum (example).

#### **3. Compared with Centrifugal Test Results**

#### *3.1. Material Properties*

British scholar Marshall [73] designed a series of model tests to study the influence of shield tunneling in sandy ground by using the geocentrifugal testing machine of Cambridge University. The design conditions and dimensions of centrifuge model test are shown in Figure 5. The bottom dimension of the test box containing model tunnels and sand is 770 mm × 147.5 mm, and the sand depth after filling is 311 mm; the axis depth of the shield tunnel model is *Zt* = 182 mm, the tunnel diameter is *Dt* = 62 mm, the section dimension of the model tunnel is shown in Figure 6, and the length of the shield tunneling section is *L0* = 147.5 mm.

**Figure 5.** Centrifuge package—tunneling in sandy ground.

**Figure 6.** Centrifuge package schematic where LVDTs (linear variable differential transformers) were used to measure vertical sub-surface soil displacements at the same locations as the lasers. Cameras were used to capture images of the soil body during tests for analysis.

The test soil samples were taken from dry Leighton Buzzard Fraction E silica sand, UK, and its parameters are shown in Table 1. This sand has a typical *D50* of 122 μm [74], a specific gravity of 2.67, and maximum and minimum void ratios of 0.97 and 0.64 [75], respectively. Uniform sand with relative density of 90% was prepared by automatic sand leaker. The elastic modulus *E* of sand soil is 20 MPa and Poisson's ratio μ is 0.4. The whole test process was carried out under the centrifugal acceleration of 75× *g* (g is the acceleration of gravity); namely, the scale ratio of the test size was 75 times. Therefore, this test is equivalent to simulating a shield tunnel with a diameter of 4.65 m (model diameter 62 mm) and a tunnel central burial depth of 13.65 m (model tunnel central burial depth 182 mm) passing through Greenfield in reality. In the process of shield excavation, the simulation of stratum loss is realized by using a cylinder sealed with liquid connected to the tunnel driven by a motor. The stratum loss rate rises at a rate of about 0.3% per minute. The above test parameters are brought into Equation (12) for theoretical calculation of stratum displacement. Considering the clarity and regularity of result image, the test results of medium volume loss of 2.5% are compared.

**Table 1.** Physical parameters of soil samples.


#### *3.2. Comparative Analysis of Results*

The comparative analysis of the results is shown in Figure 7; the vertical displacement contours of 37.5 mm and 30 mm in the central part show good consistency and regularity. In general, the centrifugal test results are nearly straight-line contours, while the theoretical calculation results are arc contours, which also leads to a large difference between 22.5 mm, 15 mm, and 7.5 mm contours, and shows a trend of increasing the difference from the center to the edge.

As for vertical displacement, the numerical value is relatively small, and its regularity is not obvious. This can also be seen from the poor symmetry of the left and right halves of the horizontal displacement in the centrifugal test nephogram. The shape and distribution position of centrifugal contours and theoretical contours are similar, especially the upper 4 mm curve and the central 2 mm curve. However, the difference between theoretical contours and centrifugal contours also shows an increasing trend from the center to the edge, the whole lines shift to the right a little, and the theoretical contour has a wider range. This is consistent with the law that the settlement trough of sand is narrower and deeper than that of clay [76].

**Figure 7.** Contour comparison of Loganathan's solution, theorical calculation of this paper, and centrifugal test results (2.5% volume loss): (**a**) Horizontal displacement contours for Greenfield; (**b**) vertical displacement contours for Greenfield.

Due to sand soil being chosen as the test soil sample in centrifugal test, Sagaseta's solution is applicable to sand soil, but another important part of our theoretical formula, Loganathan's solution, is only applicable to clay soil. Be aware of the exponential function of Loganathan' solution modify the formation loss coefficient by fitting the clay soil data, although the soil parameters substituted in Equation (12) are sand soil's, which cannot eliminate the difference between the Loganathan's ε*x,z* and the properties of sand. This leads to the fact that in the vertical displacement contours, the centrifugal contours are linear and the theoretical contours arc, which is similar to the difference between the two on the slip fracture surface of soil. Generally, clay soil has a complex composition, elastic-plastic deformation occurs easily, deformation has a time-delay property, and many other disturbing factors, all of which are very unfavorable to the test analysis. However, sand soil has a single material, of which properties are easier to control in the test, and the drainage or non-drainage conditions need not be considered. Therefore, the current centrifugal tests are mostly based on sand. We regret that we have not found a suitable clay centrifugal test as a reference. However, considering the inevitable errors in centrifugal test and theoretical calculations, the comparison results are still acceptable.

#### **4. Conclusions**


**Author Contributions:** X.T. wrote the article; H.Z. carried out the calculations and analyzed the calculation results; G.Z. and Y.Z. processed the data; P.T. offered some useful suggestions for the preparation.

**Funding:** The study was supported by the National Natural Science Foundation of China (NSFC) under Grant No. 51508037.

**Acknowledgments:** We also highly appreciate reviewers for useful comments and editors for improving the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **A Precise Prediction of Tunnel Deformation Caused by Circular Foundation Pit Excavation**

**Huasheng Sun 1,\*, Lingwei Wang 1, Shenwei Chen 1,2, Hengwei Deng <sup>1</sup> and Jihua Zhang <sup>1</sup>**


Received: 29 April 2019; Accepted: 30 May 2019; Published: 2 June 2019

**Abstract:** In comparison with tetragonal retaining structures, circular retaining structures have an advantage in terms of controlling the deformation caused by foundation excavation, and are a reasonable choice in engineering practice. Many results have been obtained regarding the effect of tetragonal excavation on the deformation of an adjacent tunnel. Nevertheless, a sufficient understanding of the circular excavation's effect on the deformation of an adjacent tunnel is currently lacking. Therefore, this study focused on the problem of precise predicting tunnel deformation below a circular excavation. A numerical model was established to calculate the tunnel deformation caused by the circular excavation. An advanced nonlinear constitutive model, known as a hypoplasticity model, which can capture path-dependent and strain-dependent soil stiffness even at small strains, was adopted. The models and their associated parameters were calibrated by centrifuge test results reported in the literature. The deformation mechanism was revealed, and the calculated results were compared with those obtained with a square excavation and the same excavation amount. The differences between the deformations caused by these two types of excavation shapes were analyzed. It was found that under equal excavation area conditions, the excavation-induced deformations of the metro tunnel below a circular excavation were approximately 1.18–1.22 times greater than those below a square excavation. The maximum tunnel tensile bending strain caused by the circular excavation was 32% smaller than that caused by the square excavation. By comparing with the measured results, it is proved that the proposed numerical method can provide effective reference for engineers to analyze soil-structure problems.

**Keywords:** advanced model; precise prediction; circular foundation pit; tunnel deformation

#### **1. Introduction**

With urban development, a large number of excavation works have been conducted. Stress relief caused by basement excavation will inevitably affect surrounding subway tunnels. Many research results have been obtained with regard to the effect of tetragonal excavation on the deformation of an adjacent metro tunnel. The main contents of this study are as follows: geometry of the excavation [1,2]; relative position of tunnel and excavation [3–7]; excavation conditions and reinforcement method [8–11]; tunnel lining stiffness [4] and diameter of tunnel lining [1]; different soil constitutive models [12]; soil relative density and retaining wall stiffness [13]; influence zone [14] and calculation methods [15–17].

However, existing studies have mainly focused on analyzing the tunnel response caused by the unloading of the tetragonal basement excavation. With the rapid development of urban construction, the development and utilization of underground space is rapidly increasing and more and more foundation basements are constructed as circular foundation basements [18–21]. Unlike the traditional

tetragonal rectangular basement, a circular basement can make full use of the circular arch effect, such that the force of the structure becomes more reasonable.

Existing circular basement studies have mainly focused on the soil pressure around the basement and the retaining structure's deformation caused by the excavation [22–30]. However, few studies have focused on the deformation of surrounding foundations and buildings, which is caused by the unloading of the circular excavation. Considering the many differences between circular and tetragonal basements, it is theoretically and practically significant to evaluate the effect of circular basement excavation on the deformation of the surrounding structure to reduce damages.

The objective of this study was to investigate the contribution of circular basement construction to the deformation of a nearby existing tunnel, and to reveal the tunnel deformation mechanism. Three-dimensional numerical analyses using advanced constitute model combined with accurate numerical simulation techniques were conducted based on the reported centrifuge model test. The effect of the excavation shape and the tunnel and soil responses were analyzed with the expectation that the results may provide a useful reference for engineering practice.

#### **2. Introduction of a Centrifuge Model Test**

A centrifuge model test, which models the interaction between a new basement excavation and an existing tunnel was carried, out as shown in Figure 1 [7]. The centrifuge acceleration of the test was 60 g. The length, width, and depth of the excavation was 300 × 300 × 150 mm, which is equal to 18 × 18 × 9 m in the prototype, respectively. The Young's modulus of the model diaphragm wall was 70 GPa, and was equal to that of the model tunnel. The model tunnel had a diameter of 100 mm, length of 1,200 mm, and thickness of 3 mm, which is equal to 6 m, 72 m, and 0.18 m in the prototype, respectively. The details of the proposed centrifuge test have been reported by Ng et al. (2013) [7]. It is also noted that the following measured results in this paper are all reported by Ng et al. (2013) [7].

**Figure 1.** Layout of the proposed centrifuge model test [7].

#### **3. Finite Element Model**

#### *3.1. Finite Element Mesh and Boundary Condition*

In this study, the finite element software ABAQUS [31] was used. A three-dimensional finite element mesh was built as shown in Figure 2. The finite element mesh in this model had a length of 1200 mm, width of 990 mm, and depth of 750 mm. The square basement excavation (on plan) was the same as that in the centrifuge model test, with a length and width of 300 × 300 mm (equal to 18 × 18 m in the prototype). Preliminary calculations checked the influence of the mesh density (13,540–86,590 elements in a structured or unstructured mesh, respectively) and of the element type (triangle or quadrilateral elements with linear interpolation). The influence of these factors on the results has been shown to be small. However, for precise and quick convergence, 86,950 elements in a structured mesh with quadrilateral element type was adopted. Based on the principle of the excavation amount being equal to the square excavation amount recorded in the test, a circular excavation was conducted with a

diameter of 339 mm (equal to 20.3 m in the prototype). The basement excavation depth was 150 mm, which is equal to 9 m in the prototype. The embedment depth of the diaphragm wall was 225 mm, which is equal to 13.5 m in the prototype. The sand and the diaphragm wall were simulated by an eight-node brick element, and the tunnel was simulated by a four-node shell element. To consider the friction boundary condition, pin supports were adopted both at the vertical sides and at the base of the model. Interface elements were adopted at the soil-diaphragm wall and at the soil-tunnel interfaces. Each interface element used is described by zero-thickness slip element, assigned with the Coulomb friction law. The friction coefficient (μ) and limiting relative displacement (γlim) at which slippage occurs are controlled by two input parameters for each slip element. The interface friction coefficient, μ, is derived from μ = tanδ, where δ is the interface friction angle, which is taken as 20◦ (i.e., 2/3 of the critical friction angle of soil). The limiting displacement of 5 mm is assumed to achieve full mobilization of the interface friction. It is noted that the initial compaction due to air pressure was not taken into account.

**Figure 2.** Finite element model (in model scale, unit: mm).

#### *3.2. Constitutive Model and Model Parameters*

In the references paper [12], the ability of three different soil constitutive models to predict a tunnel's response to basement excavation has been evaluated. It is found that the hypoplastic model yielded the best predictions of tunnel heave among the adopted models. The nonlinearity of soil can be captured by a hypoplastic (HP) constitutive model. Various HP models have been developed in the 1990s (Kolymbas 1991; Gudehus 1996; Von Wolffersdorff 1996; Wu et al. 1996) [32–35] as well as recently (Mašín 2012 2014, 2013) [36–38]. The model proposed by Von Wolffersdorff (1996) [34] was adopted to describe the behavior of Toyoura sand. This model was incorporated into the software package ABAQUS using open-source implementation which can be freely downloaded from the web (Gudehus et al. 2008) [39]. The model specifies eight material parameters (ϕc, *h*s, *n*, *e*d0, *e*c0, *e*i0, α and β). Niemunis and Herle (1997) [40] improved the model for predictions of small-strain stiffness and the recent stress history, leading to five additional parameters (*m*T, *m*R, *R*, β<sup>r</sup> and χ). See Table 1 and the literature mentioned above for their physical meaning.


**Table 1.** Soil parameters adopted in the hypoplastic model.

Note: (1) Herle and Gudehus, 1999 [41]. (2) Obtained by fitting test results from Maeda and Miura (1999) [42] and Iwasaki et al. (1978) [43].

Six parameters of Toyoura sand (ϕc, *h*s, *n*, *e*d0, *e*c0 and *e*i0) were obtained from Herle and Gudehus (1999) [41], while the triaxial test results reported by Maeda and Miura (1999) [42] were used to calibrate the parameters of α and β. Five parameters (*m*T, *m*R, *R*, β<sup>r</sup> and χ) of the intergranular strain can be calibrated from the stiffness degradation curve of Toyoura sand (Iwasaki et al. 1978) [43], as shown in Figure 3. The void ratio of soil was considered as a state variable in the HP model. The parameter selection list is presented in Table 1. The parameter selection procedure and values have been reported by Ng et al. (2015) [12]. The tunnel and diaphragm wall can be considered as an elastic material. Their elastic modulus and Poisson's ratio were 70 GPa and 0.2, respectively.

#### *3.3. Numerical Modelling Procedures*

The gradual unloading method was adopted to simulate the excavation. Figure 3 shows the flow chart of specific simulation process and as follows.


**Figure 3.** Flow chart of the numerical modelling procedures.

#### **4. Result Interpretation**

#### *4.1. Tunnel Heave in Longitudinal Direction*

Figure 4 shows the measured tunnel heave computed by numerical analyses for both the square and the circular basement excavations. It can be seen that the computed tunnel heave was generally in agreement with the measured magnitude and distribution values. The tunnel heave in the existing tunnel was caused by stress relief, as shown in Figure 5. With the increase of the basement excavation depth, the tunnel heave increased. After the basement excavation, the measured tunnel heave reached 0.07% *H*e (*H*e is the basement depth) at a distance of 0.4 *H*e from the basement center, and reduced to zero at a distance of 2.4 *H*<sup>e</sup> from the basement center. The computed maximum value of the tunnel heave was 0.084% *H*<sup>e</sup> and 0.099% *H*<sup>e</sup> for the square and circular basement excavations, respectively, with a similar trend distribution. By comparing the computed tunnel heave caused by a different basement shape, it was found that the maximum tunnel heave caused by the unloading of the circular excavation was 1.18 times greater than that of the square excavation with the same excavation amount. According to the provisions of the Land Transport Authority [44], the allowable maximum tunnel heave value is 15 mm (i.e., 0.170% *H*e, as shown in Figure 4). The computed tunnel heave ranged in the standard allowable value and also in the measured value. From the above analysis, it was determined that the tunnel heave caused by the circular basement excavation was 18% larger than that caused by the square basement excavation. Additionally, if a circular basement is adopted, and assuming that the tunnel heave caused by the square basement is in a critical state, it is essential to take corresponding reinforcement measures.

**Figure 4.** Measured and computed tunnel heave.

**Figure 5.** Vertical stress changes in soil element around tunnel crown.

#### *4.2. Stress Distribution around the Tunnel Lining*

To explain the reason for the occurrence of tunnel heave, Figure 5 shows the stress distribution in the soil element around the tunnel crown. The stress increases are shown as positive values, while the stress decreases are shown as negative values. The basement excavation with two different shapes predicted almost similar vertical stress patterns for the soil elements around the tunnel crown. A great decrease in the vertical stress beneath the basement was observed after the excavation. Additionally, it was observed that the vertical soil stress was uniform beneath the basement, and increased greatly just beneath the diaphragm wall. The vertical stress increased by 64 kPa and 75 kPa for the square and circular basement, respectively. The largest stress change for the soil element was caused by the circular basement around the tunnel crown and was 17% larger than that of the square basement. Thus, the tunnel heave caused by the circular basement excavation was larger than that caused by the square basement excavation, as shown in Figure 4. At a distance 0.6 *H*<sup>e</sup> from the diaphragm wall location, the stress change in the soil element was less than 20 kPa. Based on the Building department of the

government of Hong Kong Special Administrative Region, HKSAR [45], the stress change in this case did not exceed the proposed limit value (i.e., ±20 kPa).

#### *4.3. Tunnel Diameter Change*

Figure 6 shows the tunnel diameter (*D*) change varying with *H*ec/*C* (defined as the unloading ratio), where *H*ec is the current excavation depth and *C* is the tunnel cover depth. The positive value represents the elongation, while the negative value represents compression. The tunnel lining caused by both basement excavation types was elongated vertically and compressed horizontally. It is clear that with the increase of the basement excavation depth, the elongation (Δ*DV*) and compression (Δ*DH*) of the tunnel lining increased. The square and circular basement resulted in a maximum change in the tunnel diameter by 0.09% *D* and 0.11% *D*, respectively. The maximum tunnel diameter change caused by the circular basement excavation was 22% larger than that caused by the square basement.

To better understand the tunnel diameter change for the two different types of basement excavation, the earth pressure change around the tunnel is shown in Figure 7, where it can be seen that the maximum earth pressure change occurred in the soil around the tunnel crown, while a smaller change occurred at the springlines, and little change occurred at the invert. Because the earth pressure reduction was larger at the crown than at the springlines, the tunnel lining was elongated vertically and compressed horizontally as shown in Figure 6. By comparing the two types of basement excavation, the earth pressure change at the tunnel crown and invert was almost the same because of the same excavation amount assumed by the model. However, owing to the effect of arching, the earth pressure change for the other soil elements of the square basement was smaller than that in the circular basement excavation. Thus, the tunnel diameter change caused by the circular basement excavation was larger than that caused by the square basement, as shown in Figure 6. Based on the recommendations of the British Tunnelling Society [46], the maximum tunnel lining change ((Δ*DV* + Δ*DH*)/*D*) should not exceed 2%. Thus, it was determined that the maximum deformation of the tunnel (i.e., 0.16% *D*) was smaller than the proposed allowable limit.

**Figure 6.** Comparison of tunnel diameter change.

**Figure 7.** Earth pressure change around the tunnel lining (absolute value).

#### *4.4. Bending Strain at the Tunnel Lining*

The bending strains computed at the tunnel lining in the transverse direction were compared with the measured values, as shown in Figure 8a. The positive value represents the tensile strain, and the negative value represents the compressive strain. In comparison with the measured results, the computed results for both basement types exhibited similar trends. Moreover, the symmetrical bending strain distribution with regard to the tunnel centerline could be determined. The compressive strain occurred at the tunnel springlines, while the tensile strain occurred at the other tunnel elements. In the transverse direction, the maximum tensile bending strain was 91 με and 108 με for the square and circular basement excavation, respectively. There was only a small change in the transverse direction for these two excavations types, not only in the magnitude, but also in the distribution. Thus it was concluded that, in the transverse direction, the basement shape had little influence on the tunnel's bending strain. Based on the American Concrete Institute [35], the maximum tensile strain was 150 με. The computed bending strains did not exceed the given allowable limit nor the measured value, assuming that zero bending strain occurred before the basement construction.

Figure 8b shows the bending strains computed after the excavation in the longitudinal direction in comparison with the measured values. The positive sign denotes the tensile strain, while the negative sign denotes the compressive strain occurring at the tunnel lining crown. Similar trends were observed in the measured and computed results. With regard to the basement center, it was found that the distribution of the bending strain was symmetrical. For the square and circular basement, the maximum tensile strain computed by the numerical method was 86 με and 65 με, respectively. The computed maximum tunnel lining bending strain in the longitudinal direction, which was caused by the circular excavation, was 32% smaller than that caused by the square excavation. Based on ACI [47], cracks may occur in the tunnel lining assuming that the tensile strain exceeds 85 με before the circular basement excavation.

**E Figure 8.** Additional bending strains in tunnel lining along the (**a**) transverse and (**b**) longitudinal direction.

#### **5. Conclusions**

This paper studies the tunnel deformation caused by circular excavations; the problem is tackled numerically by means an advanced nonlinear hypoplasticity constitutive model, which is able to capture the path-dependent and strain-dependent soil stiffness, even at small strains. The numerical analysis is accompanied by an experimental calibration of the parameters based on centrifuge tests. The effect of the excavation shape and the tunnel and soil responses were analyzed. The precise prediction method is useful for practical applications and can represent a key reference the design of

tunnel structures, accounting for the soil-structural interaction. The comparisons of the predicted results can be regarded as the verification of different design analyses, carried out by practicing engineers. Based on the comparisons between measured and computed results, the following conclusions may be drawn:


**Data Availability Statement:** The computed data used to support the findings of this study are included within the article. Previously reported measured data were used to support this study and are available at [https://doi.org/10.1139/cgj-2012-0423, https://doi.org/10.1139/cgj-2014-0361]. These prior studies (and datasets) are cited at relevant places within the text as references [7,12].

**Author Contributions:** Conceptualization, H.S.; methodology, H.S.; software, H.S. and S.C.; validation, H.S. and L.W.; formal analysis, H.S. and L.W.; investigation, H.S. and H.D.; resources, H.S. and S.C.; data curation, H.S. and H.D.; writing—original draft preparation, H.S. and L.W.; writing—review and editing, H.S. and J.Z.; visualization, H.S.; supervision, H.S.; project administration, H.S.; funding acquisition, H.S.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 51708245, the Natural Science Foundation of Jiangsu Province, grant number BK20160426, and the Natural Science Foundation for Colleges and Universities in Jiangsu Province (17KJB620002, 17KJB130003, 17KJA560001).

**Acknowledgments:** We thank the anonymous reviewers for their valuable comments and suggestions. We thank Liwen Bianji, Edanz Editing China (www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript.

**Conflicts of Interest:** The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Stochastic Natural Frequency Analysis of Composite Structures Based on Micro-Scale and Meso-Scale Uncertainty**

#### **Shufeng Zhang \* and Xun Chen**

Science and Technology on Integrated Logistics Support Laboratory, National University of Defense Technology, Changsha 410073, China

**\*** Correspondence: sfzhang@nudt.edu.cn; Tel.: +86-731-8457-3396

Received: 6 June 2019; Accepted: 25 June 2019; Published: 27 June 2019

**Abstract:** Composite structure often shows undesirably significant uncertainty in its mechanical properties, which may consequently result into large stochastic variation of its natural frequency. This study provides stochastic natural frequency analysis of typical composite structures based on micro-scale (constituent-scale) and meso-scale (ply-scale) uncertainty. Uncertainty propagation across micro-scale and meso-scale is investigated. Response surface method (RSM) based on finite element modeling is employed to obtain approximate natural frequency of structures with complex shape or boundary conditions, and mean value and standard deviation of natural frequency of composite plate and cylindrical shell are derived. Differences in natural frequency statistics of composite plates and cylindrical shells derived by considering uncertainty at different scales are quantified and discussed. Significant statistical correlation between ply elastic properties and ply density is observed, and the statistical correlation is demonstrated to lay great influence on the statistics of structure natural frequency.

**Keywords:** composite; stochastic; natural frequency; uncertainty

#### **1. Introduction**

Fibre reinforced composite material and structures are extensively used in diverse engineering applications such as aircraft, rockets, ships, wind turbines, buildings, etc. Due to complexity in material microstructure and manufacture process, composite structure often shows undesirably large uncertainty in its mechanical performance, such as structure deflection [1,2], load capacity [3,4], buckling behaviour [5,6] and dynamic response [7]. Therefore, structure design in stochastic or probabilistic approaches is especially important for composite structures [8].

Since engineering composite structures are commonly required to have natural frequencies away from the operating value to avoid structure resonance, stochastic natural frequency of composite structures has received increased attention in recent years. Oh and Librescu [9] adopted stochastic finite element method (SFEM) to investigate stochastic natural frequency of composite cantilever beams. Dey et al. [10,11] derived natural frequency statistics of composite plates by an RS-HDMR approach and later used a surrogate model. Chakraborty et al. [12] investigated stochastic free vibration behaviour of laminated composite plates using polynomial correlated function expansion and also response surface method (RSM). Tawfik et al. [13] presented a neural network-based second order reliability method to obtain probabilistic natural frequency of composite plates. Yin et al. [14] investigated stochastic variability of natural frequencies of laminated plates by modal stability procedure probabilistic approach. Singh et al. [15] provided a free vibration analysis of composite cylindrical panels with random material properties by first-order perturbation technique, and similar study is shown in [16] but for composite conical shell structure. It is noticed that these current studies are based on meso-scale

(ply-scale) uncertainty such as uncertainty in ply elastic properties, ply orientation, ply density and thickness, and these meso-scale random variables are considered to be independent with one another. It is known that uncertainty of ply properties is essentially caused by stochastic variation of micro-scale (constituent-scale) parameters such as properties of fibre/matrix, fibre volume ratio, micro-structure, void, etc. [17]. Although probabilistic analysis based on meso-scale uncertainty takes less computation effort, uncertainty consideration started from micro-scale would provide insight into propagation of uncertainty across different scales [18] and also possible coupling between meso-scale random variables [19–24]. Still, it seems that little work has been reported about achieving comprehensive understanding on stochastic natural frequency of composite structures by considering uncertainty at different scales. Naskar et al. [25] showed difference in natural frequency statistics of composite plates derived from different scale uncertainty, but the uncertainties at different scales are somehow not connected. Therefore, it would be valuable to conduct comprehensive comparison on achieving natural frequency statistics of composite structures from meso-scale and micro-scale uncertainty, and any difference on the results would provide good indication or enlightenment to achieve complete description over uncertainty at different scales.

This study provides natural frequency analysis of typical composite structures based on micro-scale and meso-scale uncertainty. Meso-scale statistics are derived from micro-scale uncertainty by a combination of Monte-Carlo simulation and micromechanical model, and especially statistical correlation between meso-scale random variables is discussed. Response surface method (RSM) based on finite element modeling is employed to obtain approximate natural frequency of structures with complex shape or boundary condition. Statistics of natural frequencies of laminate plate and cylindrical shell are derived by a combination of Monte-Carlo simulation and the RSM. Statistics of structure natural frequency derived from different uncertainty considerations are compared and discussed.

#### **2. Material Statistics**

#### *2.1. Micro-Scale Statistics*

Statistics of mechanical and physical properties of common E-glass fibre and epoxy are shown in Table 1. Random variables considered in micro-scale mainly include *Ef* (Young's modulus of fibre), ν*<sup>f</sup>* (Poisson's ratio of fibre), ρ*<sup>f</sup>* (density of fibre), *Em* (Young's modulus of matrix), ν*<sup>m</sup>* (Poisson's ratio of matrix), ρ*<sup>m</sup>* (density of matrix) and *Vf* (fibre volume ratio). It has been experimentally shown that E-glass fibre tends to be mechanically isotropic [26]. Experimental observation on the statistics of fibre and matrix properties is generally rarely reported. In Table 1, the coefficient of variation (standard deviation/mean value) is assumed to be 5% if a reference is not given. All micro-scale properties shown in Table 1 are assumed to follow normal distribution in this study.


**Table 1.** Statistics of properties of fibre and matrix.

#### *2.2. Meso-Scale Statistics*

Statistics of meso-scale properties of unidirectional (UD) E-glass/epoxy are derived from micro-scale properties by a combination of Monte-Carlo simulation and bridging micromechanical model [28]. For E-glass/epoxy composite, the bridging micromechanical model is expressed as:

$$E\_1 = E\_f V\_f + E\_m V\_m \tag{1}$$

$$E2 = \frac{(V\_f + V\_m a\_{11})(V\_f + V\_m a\_{22})}{(V\_f + V\_m a\_{11})(V\_f S\_{f22} + a\_{22} V\_m S\_{m22}) + V\_f V\_m (S\_{m12} - S\_{f12})a\_{12}} \tag{2}$$

$$\nu\_{12} = \nu\_f V\_f + \nu\_m V\_m \tag{3}$$

$$\mathbf{G}\_{12} = \frac{(\mathbf{G}\_f + \mathbf{G}\_m) + V\_f(\mathbf{G}\_f - \mathbf{G}\_m)}{(\mathbf{G}\_f + \mathbf{G}\_m) - V\_f(\mathbf{G}\_f - \mathbf{G}\_m)} \tag{4}$$

where *a*11, *a*22, *a*12, *Sf*11, ... , *Sm*<sup>12</sup> are functions of the elastic properties of fibre and matrix as given in Huang [28]; *Vm* denotes matrix volume ratio (*Vm* = 1 − *Vf*). *Gm*(*Gf*) is the shear modulus of matrix (fibre) which is dependent on *Em* and ν*<sup>m</sup>* (*Ef* and ν*f*) due to isotropy of matrix (fibre). *E*<sup>1</sup> denotes normal modulus in the fibre direction; *E*<sup>2</sup> denotes the normal modulus in the transverse direction; ν<sup>12</sup> denotes the in-plane major Poisson's ratio; *G*<sup>12</sup> denotes the in-plane shear modulus. In a composite material coordinate, 1 denotes fibre direction and 2 denotes transverse direction (perpendicular to fibre direction). Lamina density is derived by the rule of mixture which is written as:

$$
\rho = \rho\_f V\_f + \rho\_m V\_m \tag{5}
$$

In Monte-Carlo method, the accuracy (ζ) on the estimation of a parameter *P* is expressed as [29]:

$$
\zeta \approx \frac{s/\sqrt{N}}{\hat{P}} \tag{6}
$$

where *P*ˆ is estimation of *P*, *N* is sampling number, and *s* is standard deviation of *P* which is estimated as:

$$s = \sqrt{\frac{1}{N} \sum\_{k=1}^{N} g(\mathbf{X}\_k)^2 - \left(\frac{1}{N} \sum\_{k=1}^{N} g(\mathbf{X}\_k)\right)^2} \tag{7}$$

where X*<sup>k</sup>* is the *k*-th sample of random variables, and *g*(X*k*) is a function on X*<sup>k</sup>* such as Equations (1)–(5). ζ represents the relative difference between the mathematical expectation of a group of samples and the true mean value of corresponding random variable. Equation (6) indicates that a large enough number of samples would provide very accurate estimation on random parameters in Monte-Carlo simulation. In this work, 1 <sup>×</sup> 105 group of samples were drawn for the micro-scale random variables listed in Table 1 and they were substituted into Equations (1)–(5) to derive statistics of meso-scale properties. According to Equation (6), the sampling number of 1 <sup>×</sup> 105 provides estimation accuracy on meso-scale random variables better than 0.02%. The derived mean values, coefficient of variation (CoV) and linear correlation coefficients of UD E-glass/epoxy are shown in Table 2. The probabilistic distribution of derived ply mechanical properties was determined by Kolmogorov–Smirnov and Anderson–Darling tests, and it was shown that *E*11, *E*22, ν12, *G*<sup>12</sup> and ρ all follow normal distribution with a confidence level at 95%.


**Table 2.** The statistics of UD E-glass/Epoxy composite.

Table 2 shows that significant statistical correlation exists between ply mechanical properties (*E*1, *E*2, *G*12, υ12), which agrees with observation in [19,20,22]. A new observation is that statistical correlation between ply density (ρ) and ply elastic properties (*E*1, *E*2, *G*12) was also significant, and this correlation seems to be rarely reported or noticed in previous studies. Actually, experimental data in Huang [28] provided a clear positive correlation between *Vf* (fibre volume ratio) and *E*1, *E*2, *G*12. Since it is known that ply density is positively dependent on the *Vf*, the theoretically predicted positive correlation between ply density and *E*1, *E*2, *G*<sup>12</sup> shown Table 2 is reasonable. Statistical correlation associated with ply density may affect stochastic vibration of composite structures as structure vibration is known to be mass dependent. Studies in later sections will demonstrate the effect of statistical correlation between *E*1, *E*2, *G*<sup>12</sup> and ρ on the stochastic free vibration of composite structures.

#### **3. Analysis of Stochastic Natural Frequency**

#### *3.1. Laminated Plate*

In this work, stochastic natural frequency of several different configurations of laminated plates subjected to different boundary conditions as listed in Table 3 was investigated, including cross-ply laminate, angle-ply laminate, quasi-isotropic laminate and laminate of random ply orientations. For each laminate configuration, statistics of first order natural frequency were derived in three scenarios: (a) considering micro-scale uncertainty of constituent material properties as listed in Table 1, termed in short as 'micro-scale uncertainty'; (b) considering meso-scale uncertainty as listed in Table 2, and their statistical correlation was also considered, termed in short as 'correlated meso-scale uncertainty'; (c) considering meso-scale uncertainty of ply properties as listed in Table 2, but the statistical correlation between ply mechanical properties and also the density was neglected, termed in short as 'independent meso-scale uncertainty'. The flowchart of calculation process of stochastic natural frequency by considering the three different scenarios of uncertainty is shown in Figure 1. It is important to notice that most current studies (as introduced in Section 1) were based on 'independent meso-scale uncertainty', and it is necessary to discuss whether 'independent meso-scale uncertainty' provides reasonable estimation on the statistics of resonance frequency of composite structures.



**(a) Micro-scale uncertainty**

#### **(b) Correlated meso-scale uncertainty**

**(c) Independent meso-scale uncertainty**

**Figure 1.** Flowchart of calculation of natural frequency statistics by different uncertainty consideration.

The statistics of first order natural frequency were derived by a combination of Monte-Carlo simulation and FE modeling. Four-node linear-elastic shell element (Shell181 in Ansys 14.5) is used to model the laminate. A mesh convergence study was conducted, which showed that a mesh size of 0.01 × 0.01 m provided a convergent solution on the first order natural frequency. Considering heavy computation cost of direct Monte-Carlo simulation on FE modeling, response surface method (RSM) with quadratic expression was employed to derive an approximate surrogate model to map random variables to the first order natural frequency, as shown in Equation (8):

$$\log(\mathbf{X}) = a + b\_i \sum\_{i=1}^{n} x\_i + c\_{ij} \sum\_{i=1}^{n} x\_i \sum\_{i=1}^{n} x\_j \tag{8}$$

where *g* represents the first order natural frequency, *xi* represents micro-scale or meso-scale random variable shown in Table 1 or Table 2, and *a*, *b* and *c* are regression terms. In this work, the regression terms were determined by central composite design method [31]. The central composite design was composed of one central point, two axis points per input variable and factorial points at the corners of the hypercube. Twenty-seven design points were employed for the meso-scale uncertainty scenario and 79 design points were employed for the micro-scale uncertainty scenario. In comparison to other experimental design methods such as the Box–Behnken method, the central composite design method

could provide high quality predictions over the entire design space, and it is more flexible considering the number of input random variables [32].

For all laminate configurations shown in Table 3, 1 <sup>×</sup> 105 group samples of random variables were drawn to conduct stochastic analysis based on the surrogate model derived by the RSM, and corresponding relative error on estimation of the mean value of first order natural frequency was less than 0.01% (calculated by Equation (6)). For the scenario of 'correlated meso-scale uncertainty', Gibbs sampling was employed to draw samples of meso-scale random variables to account for their correlation [21], by the procedure as follows:

	- 1. *x* (*t*+1) <sup>1</sup> <sup>∼</sup> *<sup>p</sup>*(*x*1|*x<sup>t</sup>* <sup>2</sup>, *xt* <sup>3</sup>, ··· *<sup>x</sup><sup>t</sup> n*) 2. *xt*<sup>+</sup><sup>1</sup> <sup>2</sup> <sup>∼</sup> *<sup>p</sup>*(*x*2|*xt*+<sup>1</sup> <sup>1</sup> , *xt* <sup>3</sup>, ··· *<sup>x</sup><sup>t</sup> n*) 3. ········· 4. *xt*<sup>+</sup><sup>1</sup> *<sup>j</sup>* <sup>∼</sup> *<sup>p</sup>*(*xj*|*xt*+<sup>1</sup> <sup>1</sup> , ··· , *<sup>x</sup>t*+<sup>1</sup> *<sup>j</sup>*−<sup>1</sup> , *<sup>x</sup><sup>t</sup> <sup>j</sup>*+<sup>1</sup> ··· *<sup>x</sup><sup>t</sup>* 5. ·········
	- 6. *xt*<sup>+</sup><sup>1</sup> *<sup>n</sup>* <sup>∼</sup> *<sup>p</sup>*(*xn*|*xt*+<sup>1</sup> <sup>1</sup> , *xt*<sup>+</sup><sup>1</sup> <sup>2</sup> , ··· *xt*<sup>+</sup><sup>1</sup> *<sup>n</sup>*−1)

where *t* denotes the number of sample, and *p*() denotes probability density function.

*n*)

Table 4 lists mean value and standard deviation of derived first order natural frequency of different laminate plates. It is seen that mean values of the first order natural frequency derived from the three different consideration of uncertainty were almost identical, with a difference less than 0.01%. The standard deviations of the first order natural frequency derived from 'micro-scale uncertainty' and 'correlated meso-scale uncertainty' were also very similar, with a difference less than 2%. This small difference could be introduced by the process of probabilistic density function fitting of meso-scale random variables. However, significant difference was seen between the standard deviation derived by 'independent meso-scale uncertainty' and 'correlated meso-scale uncertainty' or 'micro-scale uncertainty', where the standard deviation derived by 'independent meso-scale uncertainty' was around 10% larger.


**Table 4.** Statistics of first order natural frequencies of laminate plate (unit in Hz).

\* Difference from its counterpart value of 'Correlated Meso-Scale'.

Figure 2 shows histogram and probability distribution of the first order natural frequency of Laminate 1 derived from 'independent meso-scale uncertainty' and 'correlated meso-scale uncertainty'. From the histogram shown in Figure 2a, it is shown that the 'correlated meso-scale uncertainty' provides more samples at the middle region and less samples at the tail region. By Kolmogorov–Smirnov and Anderson–Darling tests, the probability distribution of the first order natural frequency tended to follow normal distribution. From the fitted probability density function by normal distribution as shown

in Figure 2b, it is seen that there is notable difference in the tail region between the two probability density functions. The cumulative probability of the first order natural frequency is shown in Figure 3, which shows that the 'independent meso-scale uncertainty' provides a notable overestimation on the cumulative probability. This demonstrates a necessity to account for statistical correlation between meso-scale uncertainty to derive accurate stochastic natural frequency.

**Figure 2.** Probabilistic distribution of first order frequency of Laminate 1 by considering correlated and independent meso-scale random variables: (**a**) histogram; (**b**) fitted by normal distribution.

**Figure 3.** Cumulative probability of first order frequency of Laminate 1, derived by considering correlated and independent meso-scale random variables.

To achieve better understanding on the contribution of different uncertainty to the stochastic variation of the natural frequency, linear correlation coefficients between the first order natural frequency of Laminate 1 and meso-scale random variables were derived, as shown in Figure 4. The linear correlation coefficient reflected a stochastic dependence or sensitivity of stochastic variation of the natural frequency over input random variables. It is seen that neglecting correlation between meso-scale uncertainty suppressed the stochastic dependence on *E*<sup>22</sup> and *G*12, which was similar to observation in [1] where stochastic deflection of a laminate structure was investigated. An interesting observation is on the dependence over laminate density. If the statistical correlation between the meso-scale input random variables was ignored, the stochastic dependence of the first order natural

frequency over ply density was negative. It seems to be under expectation as it is known that larger density results into smaller natural frequency if plate stiffness and boundary conditions remain unchanged. However, from a multi-scale viewpoint, the ply density is actually positively correlated with ply modulus, and a larger ply density also indicates a larger stiffness. Hence, for the case of 'correlated meso-scale uncertainty', a positive correlation was seen between the natural frequency and the ply density. This highlights the importance to consider statistical correlation between meso-scale mechanical properties to achieve comprehensive understanding on the stochastic free vibration of laminate plates.

**Figure 4.** Linear correlation coefficient between first order natural frequency of laminate plate and meso-scale random variables.

For Laminates 2, 3 and 4 listed in Table 3, similar characteristics of the standard deviation, cumulative probability and sensitivity of the first order natural frequency were obtained as compared to Laminate 1, considering the discrepancy between the three different uncertainty considerations.

#### *3.2. Laminated Cylindrical Shell*

Laminated cylindrical shell structures are widely used in airspace, aerospace and energy engineering fields, with applications such as the dual launch system (SYLDA) and interstage skirt structure (ISS) of Ariane 5 launcher [33]. Especially, circular cut-outs are made on the laminated cylindrical shell of the SYLDA. In this section, stochastic natural frequency of a laminated cylindrical shell with circular cut-outs is investigated. The geometry of the laminated cylindrical shell is shown in Figure 5a, which was based on a scaled model of the SYLDA of Ariane 5 launcher where the shell diameter is equal to the shell length and the cut-out diameter is around 0.1–0.2 times of shell length. The laminate was set as eight-layered E-glass fibre/epoxy material, and the ply orientation was (±θ)4. A ply with an orientation at 0 indicates that fibre is aligned in the cylinder length direction. Ply thickness was selected as 0.45 mm. Both ends of the laminated cylindrical shell were clamped. FE modeling was employed to obtain its first order natural frequency, as shown in Figure 5. The laminate was modeled by four-node linear-elastic shell element (Shell181 in Ansys 14.5), and a mesh size of 5 × 5 mm was found to be fine enough to provide a convergent solution. FE solution on the first order vibration mode of (±45)4 laminate cylindrical shell is shown in Figure 5b.

Similar to the analysis of laminate plates in last section, a combination of RSM and Monte-Carlo simulation is adopted to obtain the statistics of the first order natural frequency, by considering 'micro-scale uncertainty', 'independent meso-scale uncertainty' and 'correlated meso-scale uncertainty'. The derived standard deviation of the first order natural frequency is shown in Figure 6. It is seen that the difference on the standard deviation between the three difference uncertainty consideration was fairly small for θ at 0◦, 30◦, 45◦ and 60◦. However, for the situation of θ = 90◦, the standard

deviation derived from 'independent meso-scale uncertainty' was 10.4% larger than that of 'correlated meso-scale uncertainty' or 'micro-scale uncertainty'. This difference agrees with that of laminate plates shown in Table 4, which again demonstrates that ignoring the statistical correlation between meso-scale uncertainty probably leads to overestimation of the standard deviation of structure natural frequency.

**Figure 5.** (**a**) Geometrical dimension and mesh of laminate cylindrical shell; (**b**) first order vibration model of (±45)4 laminate.

**Figure 6.** Standard deviation of first order natural frequency of (±θ)4 laminate cylindrical shell derived by considering different uncertainty.

#### **4. Discussion**

Fibre reinforced composite is essentially a type of multi-scale material. Uncertainty consideration from different scales would provide more comprehensive understanding on accurate uncertainty characterization. Previous work has highlighted possible existence of significant correlation between ply mechanical properties [19,22]. This correlation would affect reliability of composite structures considering structure failure defined by structure deformation, fracture or buckling [20]. It is known that natural frequency of composite structures not only depends on ply elastic properties but also laminate density. Therefore, as a new observation on the statistical correlation associated to ply density, it would be interesting to see how much the correlation contributes to the statistics of structure natural frequency. Here, a further comparison is conducted on deriving the statistics of structure natural frequency by considering two scenarios: (1) considering only correlation between ply mechanical properties (*E*11, *E*22, *G*12, ν12); (2) considering correlation between ply mechanical properties and also density (*E*11, *E*22, *G*12, ν12, ρ). The derived standard deviations of the first order natural frequency of Laminate 1 (see Table 3) and (90)8 laminate cylindrical shell are shown in Table 5. It is interesting to observe that neglecting correlation associated with ply density results into 30–50% overestimation of

the standard deviation, while neglecting all correlation between meso-scale uncertainty results into around 10% overestimation.


**Table 5.** Standard deviation of first order natural frequency considering different uncertainty.

\* Difference from its counterpart value of 'meso-scale correlation (*E*11, *E*22, *G*12, ν12, ρ)'.

Overestimation of standard deviation by 30–50% would lead to very large error in corresponding cumulative probability estimation. Therefore, it is important to account for the statistical correlation associated with ply density to obtain accurate stochastic free vibration behaviour of laminated structures.

A recent study [25] reports that larger standard deviation of natural frequency was obtained for the case of 'micro-scale uncertainty' than that of the 'independent meso-scale uncertainty', which seems to contradict with our observations in Table 4 and Figure 5. Actually, the comparison in [25] was conducted by assuming identical coefficient of variation (CoV) for both meso-scale and micro-scale uncertainties, for example the CoV of all of meso-scale and micro-scale random variables were assumed to be 5%. However, Tables 1 and 2 tell that micro-scale uncertainty with CoV at 5% actually results in CoV of meso-scale uncertainty at around 10%. The uncertainty amplification from micro-scale to meso-scale was also addressed in [18]. The uncertainty amplification explains why larger standard deviation was obtained for the case of 'micro-scale uncertainty' in [25]. Meanwhile, statistical correlation between meso-scale uncertainty was not considered and discussed in [25].

#### **5. Conclusions**

The present study provides statistical natural frequency of composites structures derived from micro-scale and meso-scale uncertainty, for the circumstance composite structures are composed of plies with unidirectional fibre reinforcement. Meso-scale uncertainty is obtained from micro-scale uncertainty by a combination of Monte-Carlo simulation and bridging micromechanical model. It is observed that significant statistical correlation may not only exist between ply elastic properties but also associated with ply density. The statistical correlation between meso-scale physical properties is commonly ignored in current study on stochastic free vibration of composite structures. Comparison is conducted on natural frequency statistics of composite plates and cylindrical shell considering 'micro-scale uncertainty', 'correlated meso-scale uncertainty' and 'independent meso-scale uncertainty'. Several conclusions are drawn as follows:


It is very important to achieve accurate statistics of structure natural frequency in engineering design, so that potential structure resonance could be well avoided. Therefore, this work further suggests establishing standard experimental approaches to obtain statistical correlation between ply elastic properties and ply density.

**Author Contributions:** Conceptualization, S.Z. and X.C.; methodology, S.Z.; software, S.Z.; formal analysis, S.Z.; investigation, S.Z. and X.C.; writing—original draft preparation, S.Z.; writing—review and editing, X.C.; funding acquisition, S.Z.

**Funding:** This research was funded by National Natural Science Foundation of China, grant number 51405501.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
