Isogeometric Analysis of Longitudinal Displacement of a Simplified Tunnel Model Based on Elastic Foundation Beam

1  Introduction

Elastic foundation beams are widely utilized in the foundation of industrial, civic, and agricultural building facilities. Analyzing the deformation of elastic foundation beams has always been a research hotspot [1]. The performance of the building is obviously affected by the longitudinal deformation of the beam on the elastic foundation. For example, the tunnel can be viewed as the elastic foundation beam, according to YUKIS Shiba’s longitudinal equivalent continuous model [2]. The design of tunnel concrete lining structure is mostly transverse, with little or no consideration of the influence of longitudinal deformation. However, as more shield tunnels are built and operated, problems such as water seepage, leakage, longitudinal tensile cracking, and excessive differential settlement due to longitudinal deformation are becoming more prevalent [3].

The elastic foundation model mainly includes the Winkler model, semi-infinite elastic foundation model, and two-parameter foundation model proposed by Filonelko-Borodich et al. [4,5]. Furthermore, scholars [6–9] at home and abroad have put forward many theories and calculation methods to calculate the mechanical properties and vertical deformation of beams. In general, the Winkler foundation model corresponds to reality when the compressible soil layer thickness is less than half of the short side of the foundation; the semi-infinite elastic body hypothesis is appropriate for clay foundations; and, while the two-parameter elastic foundation model has been further developed in theory, it is difficult to apply in practice due to the difficulty in determining the parameters.

The numerical simulation based on the finite element approach has been widely employed in the settlement analysis of elastic foundation beams, and several academics [10,11] have conducted extensive research on the subject. For instance, Huang et al. [12] developed a new finite element method for the first time, based on the theory of beams on elastic foundation and the variational principle. Lou et al. [13] developed cubic finite element grids to study super-long beams supported on elastic foundations using a unified dimensionless parameter. Chen [14] suggested a new numerical approach for resolving the beams on elastic foundations problem. The governing differential equations defined on all elements are discretized using differential quadrature (DQ), and the numerical results of the solutions of beams on elastic foundations de-rived using DQEM are reported. Gao et al. [15] recently used the mixed finite element approach of independent discretization for foundation and beam elements to solve the nonlinear contact problem of a large deformation beam on an elastic foundation.

As for other works, Yin [16] expressed their own opinions. An approach was proposed for getting closed-form solutions for a reinforced Timoshenko beam on an elastic basis under any pressure load. Binesh [17] developed the radial point interpolation method, a meshless approach for analyzing two-parameter beams. Eisenberger et al. [18] established the beam element’s exact stiffness matrix on an elastic foundation. Hasan [19] researched crack-induced variations in the eigenfrequencies of the elastic foundation beam. The meshless analysis of beams on an elastic basis is introduced by Dodagoudar et al. [20].

Because the general thin beam plate and shell structure require C1 continuity displacement interpolation function, and the non-uniform rational B-spline [21] (NURBS) has high-order continuity, the IGA approach was employed to compute the settlement of the Winkler foundation beam in this study. Professor Hughes et al. [22] established the IGA method, a new numerical method based on spline theory to unite Computer Aided Design (CAD) and Computer Aided Engineering (CAE). Since the pioneering work of Hughes, IGA has been widely used in a variety of engineering applications [23–28]. Initially, Kiendl et al. [29,30] proposed using the NURBS basis function to perform IGA of multi-piece spliced Kirchhoff-Love elements, ensuring C1 continuity in the interior of each shell element. After that, Benson et al. [31] proposed the IGA method of Ck continuity of arbitrary order. Finally, Li et al. [32] established an IGA method for static analysis of beams and plates based on the third-order shear deformation theory (TSDT).

Both FEA and IGA employ the isoparametric concept, which means that the same basis is used for geometry and analysis. One obstacle, however, is that NURBS is not as straightforward as Lagrange polynomials. In order to solve this problem, many scholars have used IGA based on the Bézier extraction method which is an isometric analysis element structure similar to FEM and simplifies the implementation of isometric analysis in the finite element environment to solve crack [33] and plate buckling analysis [34] problems. Moreover, Nguyen [35] gave the Bézier extraction process of NURBS and T-splines in detail and the design of the IGA program based on Bézier extraction. Lai et al. [36] used the three-dimensional IGA method based on Bézier extraction to calculate the problem of a thick-walled cylinder and hollow thick-walled sphere. The results show that the convergence and accuracy of IGA based on Bézier extraction and traditional IGA are the same. Borden et al. [37,38] introduced an isogeometric finite element data format that extracts NURBS and T-splines using Bézier and can be immediately incorporated into the finite element software. Furthermore, based on the Bézier extraction operator, Scott et al. [39] proposed the theory of constructing arbitrary complex structures with hierarchical splines, making the surface of complex bodies smoother.

Although IGA has the advantages of high geometric accuracy, high continuity, and high precision, it is different from the C0 continuous shape function of the conventional FEM, and its programming is complicated. In this paper, we introduce a method combining IGA and Bezier extraction operators to study the vertical deformation of tunnels. The method is to decompose the NURBS function into a linear combination of Bernstein polynomials, thereby realizing the decomposition of NURBS elements into C0 continuous Bézier elements, which are similar to Lagrangian elements, making it easy to integrate into existing FEM frameworks and ensure the same accuracy. In addition, applying it to the simplified tunnel model based on the elastic foundation beam can effectively exert its advantages of higher accuracy than the FEM, thereby providing a theoretical basis for the longitudinal design of the tunnel.

2  Isogeometric Analysis

2.1 B-Spline and NURBS Basis Functions

The B-spline basis function is composed of a non-subtractive real sequence of node values, called a node vector Ξ={ς1,ς2,…,ςi,ςi+1,…,ςn+pmax+1}(ςi≤ςi+1), where ςi is the node. pmax is the highest polynomial order, and the number of basis functions is represented by n. After determining the B-spline basis function’s polynomial order, the node vector and the B-spline basis function group correspond one to one. Given the node vector, the B-spline basis function [40] is as follows:

Ni,0={1 if ςi≤ς<ςi+10 otherwise, for p=0,(1)

Ni,p(ς)=ς−ςiςi+p−ςiNi,p−1(ς)+ςi+p+1−ςςi+p+1−ςi+1Ni+1,p−1(ς), for p≥1.(2)

The B-spline curve can be evaluated by basis function, and control point coordinates set {Pi} as

C(ς)=∑i=1nNi,p(ς)Pi.(3)

The B-spline curve can be extended to the B-spline surface. The node vectors of the two directions are Ξ={ς1,ς2,…,ςn+p+1} and H={τ1,τ2,…,τm+q+1}, and the control point is Pi,j, where m and n are unary basis functions in both directions, and the B-spline surface equation can be expressed as

S(ς,τ)=∑i=1n∑j=1mMj,q(τ)Ni,p(ς)Pi,jPi,j∈R2.(4)

The quadratic NURBS basis functions and curves are shown in Figs. 1a and 1b, respectively. According to the definition method of the B-spline basis function and the introduction of weight, the definition of the NURBS basis function is as follows:

Rip(ς)=Ni,p(ς)ωiW(ς)=Ni,p(ς)ωi∑i=1nNi,p(ς)ωi.(5)

The NURBS curve is defined as:

C(ς)=∑i=1nRi,p(ς)Pi.(6)

The NURBS surface can be obtained from the tensor product of ς and τ in two coordinate directions:

S(ς,τ)=∑i=1n∑j=1mRi,jp,q(ς,τ)Pi,j.(7)

Ri,jp,q(ς,τ) is the basic function of the NURBS surface:

Ri,jp,q(ς,τ)=Ni,p(ς)Mj,q(τ)ωi,j∑i=1n∑j=1mNi,p(ς)Mj,q(τ)ωi,j.(8)

Figure 1: Quadratic NURBS basis functions and curves for a knot vector

2.2 Bézier Extraction Method

A sequence of NURBS basis functions is decomposed into linear combinations of Bernstein polynomials using the Bézier extraction procedure. Thus, the NURBS element is decomposed into a C0 continuous Bézier element. The Bernstein polynomial is defined as [41]

Bi,p(ς)=12[(1−ς)Bi,p−1(ς)+(1+ς)Bi−1,p−1(ς)].(9)

where B1,0(ς)≡1.

The expression of the Bézier curve is as follows:

C(ς)=∑i=1p+1Bi,p(ς)Pi=PTB(ς).(10)

The node ς∗∈[ςm,ςm+1) is inserted into the node vector Ξ to form a new node vector Ξ={ς1,ς2,…,ςm,ς∗,ςm+1,…,ςn+p+1}. After the node vector is updated, the control points and weights previously matched with the node vector must also be updated. The weights of the control points and their coordinate update equations are obtained as follows:

{wi1=w1i=0,1,…,m−pwi1=(1−αi)wi−1+αiwii=m−p+1,m−p+2,…mwi1=wi−1i=m+1,m+2,…,n+1(11)

{Pi1=Pii=0,1,…,m−pPi1=(αiwiPi+(1−αi)wi−1Pi−1)/wi1i=m−p+1,m−p+2,…mPi1=Pi−1i=m+1,m+2,…,n+1(12)

αi={1i=0,1,…,m−pς∗−ςςi+p−ςii=m−p+1,m−p+2,…m0i=m+1,m+2,…,n+1.(13)

The shape of the B-spline curve is the same as that of the Bézier curve if the existing nodes are inserted into the original B-spline’s node vector and the degree of repetition is equivalent to the curve’s order. At this time, the continuity of the curve and the continuity between the elements do not change [42].

Bézier decomposition is a node embedding operation. After getting the expression of αi, the Bézier extraction operator Cj is deduced from the new node {ς1¯,ς2¯,…ςj¯,…,ξ¯k}. αji represents the i th α after embedding the j th node vector, where i=1, 2,…, n+j. The following matrix could define the operator:

Cj=[α11−α2000α21−α300⋱αn+j+11−αn+j].(14)

The change of control point after node embedding is as follows:

P¯j+1=(Cj)TP¯j,P¯1=P.(15)

Let the final control point P¯k+1=Pc, where k is the number of node embeddings.

CT=(Ck)T(Ck−1)T⋯(C1)T(16)

Pc=CTP(17)

According to the B-spline curve Eq. (3), the geometric parameters of the Bézier curve after embedding the nodes and the original B-spline curve are unchanged. Therefore, the NURBS curve is shown below:

C(ς)=PTN(ς).(18)

We can deduce the link between the B-spline basis function and the Bernstein polynomial from the preceding equation.

N(ς)=CB(ς).(19)

According to Eq. (19), C is only related to node vectors and embedded new nodes, but not to control points or basis functions, so this extraction operator can also be used in NURBS.

For the denominator of the NURBS basis function, let W(ς)=∑i=1nNi,p(ς)wi and correlate with W(ς) as follows:

W(ς)=∑i=1nNi,p(ς)wi=wTN(ς)=(CTw)TB(ς)=(wb)TB(ς)=Wb(ς).(20)

where wb=CTw, wb (w is the weight of NURBS) is the weight of Bézier. Therefore, the basis function equation of NURBS using the Bézier extraction operator becomes

R(ς)=1Wb(ς)WN(ς)=1Wb(ς)WCB(ς).(21)

where W is the diagonal matrix of weight. So the relationship between Bézier and NURBS control points is as follows:

Pc=(Wc)−1CTWP.(22)

2.3 Numerical Tests

A circular beam which is a cantilever beam is used as a numerical example to demonstrate the validity of both the classic IGA and the IGA based on the Bézier extraction method. At the free end, the beam is subjected to the specified displacement u0=−0.01. Fig. 2 displays its geometry, boundary conditions, and material properties. The material is linearly elastic and under plane stress. Zienkiewicz et al. [43] provided an exact result for the strain energy of this circular beam.

U=1π(ln⁡2−0.6)≈0.029649668442377(23)

Figure 2: The geometry of the circular beam with material properties, boundary conditions and end shear

Fig. 3 displays the NURBS and Bézier element meshes that were used to model the circular beam. The number of tangential elements was chosen to be twice the number of radial elements, both polynomials of order 2, and the IGA polynomial of order 1 was not considered, as this geometry cannot be modeled with a 1st-order NURBS surface. For the rigor of IGA vs. finite element comparison, the number of global degrees of freedom is made as close as possible to that of FEM, while still keeping the number of tangential elements double the number of radial elements.

Figure 3: Circular beam, meshes of IGA, NURBS and Bézier physical mesh with control points

Table 1 shows the strain energy results to 14 decimal places for different grids and methods. The results of the three methods were validated by Zienkiewicz et al. [43]. As expected, for approximately the same number of global degrees of freedom and the same order of elements, the IGA shows strain energy that is closer to the exact solution than the FEM. The results of the FEM using the Q4 element are the furthest, while the results using the traditional IGA and the IGA method based on the Bézier extraction are very close.

3  Isogeometric Equation for Winkler Foundation Beam

3.1 Basic Equation

The deformation of the elastic foundation beam and the soil is consistent under load. As per the elastic foundation’s local deformation law, as shown in Fig. 4, the expressions of foundation reaction Pf and deflection ω of foundation beam are as follows:

Pf=kω.(24)

Figure 4: Winkler foundation model force diagram

k is the foundation reaction coefficient. After considering the strain energy, the total potential energy of the beam is as follows:

I=12∫EI(d2ωdx2)2dx−q(x)ωdx−∑i=1nfpiω(xi)−∑j=1mMjϑ(xj)+12Pfωdx.(25)

The first four terms on the right side of the equation are the strain energy of the beam, the distributed load potential energy, the concentrated load potential energy, and the concentrated moment load potential energy. Where EI is the bending stiffness, q(x) is the distributed load, fp is the concentrated load, and M is the concentrated moment. The fifth item Id=12∫Pfωdx is the contribution of foundation soil deformation energy to the total potential energy of a structural system.

After substituting Eq. (24) into Id, it is as follows:

Id=12∫lPfωdx=12∫lkω2dx.(26)

The form of unit superposition is as follows:

Id=∑eIde=∑e12∫ekω2dx.(27)

Suppose the element displacement mode of the beam is as follows:

ω=Nd.(28)

where N and d are NURBS element shape functions and nodal displacement vectors, respectively.

Simultaneous Eqs. (27) and (28):

Id=∑e12∫edTNTkNddx=∑dT(12∫eNTkNdx)d.(29)

By substituting the total potential energy Eq. (25), the additional term of foundation stiffness Kde in the element stiffness matrix can be obtained by taking the extreme value.

Kde=∫eNTkNdx.(30)

The additional term of foundation stiffness can be calculated by Eq. (30). The total stiffness matrix Ke¯ of the beam element on the elastic foundation is formed by the superposition of the Ke and the Kde.

Ke¯=Ke+Kde.(31)

where Ke is the foundation stiffness matrix. Ke=∬ABTDBdxdyt is structural element stiffness matrix and B is the strain matrix formed by the derivative of NURBS basis function. When p=q=1, Bi is as follows:

Bi=[Ri,x00Ri,yRi,yRi,x](i=1,2,3,4).(32)

where Ri,x and Ri,y are the partial derivatives of NURBS basis functions R to x and y, respectively. D is a matrix of material constants.

The expressions of nodal force and nodal displacement of beams on elastic foundation under total stiffness are as follows:

F=K¯δ.(33)

where δ is the whole node displacement array, K¯ is the total stiffness matrix of the beam element and F is the whole node load array.

3.2 Numerical Integration of Equivalent Control Point Forces

The equilibrium differential equation of the element is as follows:

Fe=Ke¯δe.(34)

The right end of the equation is the equivalent internal force. The stiffness matrix Ke¯ of the control point is multiplied by the array of displacement δe. The left end is the equivalent external force, including the following concentrated force, surface force, and physical strength (only these three cases are considered in this paper). The stiffness matrix K¯ and external force F can be obtained by numerical integration. In order to simplify the force analysis of the element, the load of the element is moved to the control point according to the principle of static equivalence.

(1)   Concentrated force: assuming that there is a concentrated load fp=[fpx,fpy]T at any control point c, the equivalent external force is as follows:

Fie=[FixeFiye]=(Ri)cTfpt.(35)

(2)   Surface force: there is a surface force q=[qxqy]T at a certain element boundary. The equivalent external force is as follows:

Fie=[FixeFiye]=∫τRiTqtds.(36)

(3)   Body force: f=[fxfy]T, the equivalent external force is as follows:

Fie=[FixeFive]=∫−11∫−11RiTfdxdyt.(37)

4  Verification

In this section, the bottom of the foundation is fully constrained, and the beam is coupled with the foundation. Therefore, this model does not consider the separation of the beam and the foundation. The structural parameters of beams on elastic foundations are shown in Table 2.

Fig. 5 shows the geometrical conditions of the elastic foundation beam and three different boundary conditions. Where μ is Poisson’s ratio and ql=382kN/m is the uniformly distributed load. A set of polynomial degrees p(p=2,3,4) with fine meshes is used. When p=2, the physical mesh and Bézier element mesh of the beam on the elastic foundation are shown in Fig. 6. The beam is subjected to uniformly distributed loads, and the three boundary conditions are as follows: clamped at both ends (C-C), simply supported at both ends (S-S), and one end clamped and one end simply supported (C-S). Fig. 7 depicts the displacement convergence of the beam’s centre point under an evenly distributed load and compared to the reference solution [44]. It is clear that the IGA method is capable of solving the elastic foundation beam problem. The polynomial order affects the precision of convergence.

Figure 5: Boundary conditions of elastic foundation beams under uniformly distributed loads

Figure 6: Elastic foundation beams meshing (8 × 8 meshing): (a) the Bézier physical mesh; (b) the Bézier control points

Figure 7: The beam’s displacement convergence of the middle point under an evenly distributed load and different boundary conditions

5  Application and Discussion

5.1 Equivalent Elastic Foundation Beam of Tunnel

This part investigates the deformation behaviour of the shield tunnel longitudinal structure, which serves as a reference basis for longitudinal design, using the equivalent continuous model and the theory of beam on elastic foundation. The actual tunnel structure is a tubular structure formed by bolted segments. The moment of inertia and bending stiffness of the section should be computed according to the actual section and material of the tunnel structure in order to simplify it to an elastic beam. The model assumes that tunnel materials are equally distributed in the transverse direction, and tunnel stiffness and structural features are the same as the simplified model in the longitudinal direction. The expression of equivalent elastic bending stiffness is as follows:

(EI)eq=cos3⁡φcos⁡φ+(φ+π2)⋅sin⁡φEtIt.(38)

The following equation determines the position of the neutral axis:

φ+ctgφ=π(12+mkblrEtAt).(39)

where kb=EkAklk is the joint bolt’s average rigidity per unit length, Ek is the bolt’s modulus of elasticity, and Ak is the bolt’s section area. The bolt length is lk, and the number of bolts is m. lr is the ring’s width; Et is the tunnel section’s modulus of elasticity; At is the tunnel’s section area, and It=π64[Dt4−(Dt−ts)4] is the tunnel’s vertical inertia moment. Dt is the tunnel’s diameter; ts is the shield tunnel’s segment thickness; φ is the corresponding angle of the position of the neutral axis in Fig. 8.

Figure 8: Schematic diagram of neutral axis position

5.2 Example of Tunnel Settlement

The two ends of a subway tunnel are connected with the station, and the soil at the front and rear of the entrance and exit section is reinforced by Φ 800 mm cement jet grouting pile. The reinforcement range is 7 m on both sides, as shown in Fig. 9. According to engineering experience in the Shanghai area, the reinforced foundation reaction coefficient is K1=5×107N/m3, and the uniform foundation reaction coefficient is K2=5×106N/m3. The tunnel’s submerged depth is 10.3 meters, and the soil layer’s average bulk density is 18×103N/m3. The self-weight of the tunnel is 51×103N per extension meter. Terzaghi earth pressure theory calculates the average linear load to be 1.2×106N/m (excluding ground overloading). The tunnel has a 6200 mm exterior diameter and a 5500 mm interior diameter. There are 17 M30, 8.8-grade bolts spread irregularly across the circumference. Table 3 shows the fundamental information about the segments and bolts. The elastic equivalent bending stiffness (EI)eq=6.68×1010N⋅m2 is computed using the equivalent continuous beam model, with the angle of the neutral axis φ=0.9635.

Figure 9: Interval tunnel calculation model

The length of the interval tunnel is L=42m. The foundation reaction coefficient is divided into three types. One is the homogeneous soil layer k1, the other is the homogeneous soil layer k2, and the last is shown in Fig. 9.

(1)   The foundation reaction coefficient is k1.

(2)   The foundation reaction coefficient is k2.

(3)   The foundation reaction coefficient is k2 in the middle and k1 at both ends.

The 10 × 10, 20 × 20, and 30 × 30 meshing of quadratic, cubic, and quartic Bezier elements is utilized in the computation to depict the displacement history of each control point of the tunnel model. The working condition combinations of three different foundation reaction coefficients are depicted in Figs. 10–12. Fig. 10 shows the comparison of the results of calculating the vertical settlement of the tunnel by the initial-parameter method [45] and the IGA method based on Bézier extraction under the different grid and basis function orders when the foundation reaction coefficient is k1. Figs. 11 and 12 are similar to Fig. 10. The foundation reaction coefficient in Fig. 11 is k2, while the foundation reaction coefficient of the tunnel in Fig. 12 is k2 in the middle and k1 at both ends. As shown in Figs. 10a and 10b, when the foundation reaction coefficient is k1 and the order of the basic function of the IGA is 2 and 3, there is a certain error between the calculation results and the results of the initial-parameter method. Compared to the corresponding graphs in Figs. 11 and 12, the error is relatively small.

Figure 10: Interval tunnel loaded by a uniform pressure: displacement of each point and convergence to the initial parameter solution in the first situation

Figure 11: Interval tunnel loaded by a uniform pressure: displacement of each point and convergence to the initial parameter solution in the second situation

Figure 12: Interval tunnel loaded by a uniform pressure: displacement of each point and convergence to the initial parameter solution in the third situation

When the foundation reaction coefficient is k1, the deformation of the foundation is the smallest, when the foundation reaction coefficient is k2, the deformation of the foundation is the largest, and when the foundation reaction coefficient is the combination of k1and k2, the distortion of the foundation is in the middle. Furthermore, when the calculated results are compared to the results of the initial-parameter method, the results suggest that utilizing a fine grid and a high-order polynomial calculation approach can result in higher precision data. Simultaneously, the data obtained by IGA and FEM are presented in Tables 4–9 to show the displacement of the tunnel model’s central point more clearly. As expected, the displacement given by IGA is closer to the exact solution than that shown by finite element analysis for roughly the same number of global degrees of freedom and the same order of elements. The higher-order element reduces the error. The finite element analysis with the Q4 element deviates the most from precise results, whereas the finite element analysis with the Q9 element is closer to the exact solution of the finest mesh.

6  Conclusion

The Bézier C0 element is used to analyze the vertical displacement of beams on an elastic basis in this work.

(1)   The investigation of the convergence and accuracy of vertical deformation of beams on elastic foundations under simple loads is established to validate the efficiency of this method.

(2)   After obtaining the stiffness of the homogeneous cylinder from the equivalent continuous calculation model, the shield tunnel is simplified to a uniform continuous beam with equivalent stiffness. The simulation results are in good agreement with the initial-parameter solution, and the accuracy is higher than that of the FEM solution, providing a foundation for longitudinal design.

(3)   The foundation reaction coefficient has a major influence on tunnel settlement, as shown by the comparison of the three different workings.

Acknowledgement: The authors gratefully acknowledge the support from the National Natural Science Foundation of China (52079128).

Funding Statement: The paper is supported by the National Natural Science Foundation of China (52079128).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.