CC BY
14
ГИАБ. Горный информационно-аналитический бюллетень / MIAB. Mining Informational and Analytical Bulletin, 2021;(4-1):81—90 ОРИГИНАЛЬНАЯ СТАТЬЯ / ORIGINAL PAPER
УДК 699.841+004.94 DOI: 10.25018/0236_1493_2021_41_0_81
ОЦЕНКА ВЛИЯНИЯ ЗЕМЛЕТРЯСЕНИЙ НА ПОДЗЕМНЫЕ СООРУЖЕНИЯ НА ОСНОВЕ РАСЧЕТОВ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ
Ян Прушка1, Вероника Павельцова1
1 Чешский технический университет в Праге, Чешская Республика
Аннотация: Статья посвящена оценке влияния землетрясения на подземные сооружения на основе расчетов методом конечных элементов. Сейсмические воздействия на подземные сооружения в целом могут быть исследованы либо с помощью квазистатического анализа, либо с помощью динамических расчетов. Динамические расчеты основаны на непосредственном применении полученных ускорений к реальной расчетной модели туннеля. При квазистатическом анализе на основе метода конечных элементов обычно предполагают, что вертикальное распространение волн давления и сдвига можно описать в рамках одномерной модели, в которой все переменные не зависят от времени и определяются исключительно одной вертикальной координатой. В данной работе был выполнен расчет поведения подземных сооружений, подверженных землетрясениям, с использованием программного обеспечения GEO FEM, которое реализует динамический анализ с использованием комбинации статических граничных условий, а также граничных условий свободного поля вдоль вертикальных границ. При этом считается, что на нижней границе волны полностью поглощаются. Этот более сложный подход к динамическим расчетам может подойти для инженерной практики, поскольку позволяет учитывать эффекты взаимодействия грунта и конструкции.
Ключевые слова: землетрясение, свободное поле, подземное сооружение, метод конечных элементов, сейсмическая нагрузка, динамический анализ, тоннель, квазистатический метод.
Для цитирования: Ян Прушка, Вероника Павельцова Оценка влияния землетрясений на подземные сооружения на основе расчетов методом конечных элементов // Горный информационно-аналитический бюллетень. — 2021. — № 4-1. — С. 81—90. DOI: 10.25018/0236_1493_2021_41_0_81.
Evaluation of the effect of earthquake on underground structures using
the finite element method
Jan Pruska1, Veronika Pavelcové1
1 Czech Technical University in Prague, Czech Republic
Abstract: The paper assesses the impact of earthquakes on underground structures by means of the finite element method. Seismic effects on underground structures can be investigated using either quasi-static analysis or dynamic calculations. Dynamic calculations involve using the observed accelerations in a real design model of a tunnel. In quasi-static finite element analysis, it is usually assumed that the vertical propagation of pressure and shear waves can be
© Ян Прушка, Вероника Павельцова. 2021
described by a one-dimensional model whose variables are all independent of time, determined exclusively by the vertical coordinate. The paper analyzes the behavior of underground structures during earthquakes, which is simulated by GEO FEM software; the software performs dynamic analysis using a combination of static boundary conditions and free-field boundary conditions along vertical boundaries. It is assumed that waves are completely absorbed at the lower boundary. This more sophisticated approach to dynamic analysis can be useful in engineering practice because it takes into account the effects of soil-structure interaction. Key words: earthquake, free-field, underground structure, finite element method, seismic load, dynamic analysis, tunnel, pseudo static.
For citation: Jan Pruska, Veronika Pavelcova Evaluation of the effect of earthquake on underground structures using the finite element method. MIAB. Mining Inf. Anal. Bull. 2021;(4-1):81—90. [In Russ]. DOI: 10.25018/0236_1493_2021_41_0_81.
Introduction
The response of underground structures on the seismic load can be solved in several ways:
An analytical solution commonly based on the free-field shear deformation method [1-7]. This approach assumes that the deformation of the structure should conform to the deformation of the soil in the free-field under the design earthquakes. It is the most conservative condition and this method has only a limited application. On the other hand, it provides a first-order estimation of the value of the stresses in lining and thus evaluating the accuracy of more complex calculations.
Pseudo-static finite element calculation provides solution for any shape of underground structure in the nonhomogeneous layered rock/soil mass with non-linear material behaviour. The basic idea of this approach is that the maximum shear deformation in free-field condition of the layer represents the maximum dynamic earthquake stress in this layer and can be input to the numerical model as boundary condition. The dynamic loading is further calculated by the finite element method as static problem.
A fully dynamic analysis of is builds upon the application of d'Alembert principle that requires the sum of all forces, including inertia forces, acting on a structure be zero. Application of this
principle to the model with many degrees of freedom yields a set of differential equations to be integrated in time. Unlike traditional static analyses, the results of dynamic analysis of earthquake effects is greatly affected by application of proper boundary conditions. The seismic waves follow traditional laws of mechanics. Thus, if the incoming wave hits the interface of two materials, part of the wave may proceed through such interface and part might be reflected back. In most cases, the geological profile cannot be simplified to a homogenous domain. Therefore, it is necessary to describe the wave behavior at the interface by employing proper boundary/interface conditions, which is not an easy task.
Free-field method
The basic principle of free-field analysis is based on the determination of deformation of original earth body free of the underground structure associated with the maximum acceleration of shaking soil reduced with respect to the depth of expected underground structure, e.g. tunnel. Such deformations, strains, are then directly applied to the structure, see Fig. 1. Shear strain y of the analyzed body caused by seismic loading is then given by:
Deformation of structure Fig. 1. Deformation of underground structure in free-field method [1]
where gV — vertical stress, G — shear modulus, PGA — peak ground acceleration at the terrain level, g — gravitational acceleration, R — reduction coefficient accounting for depth.
The horizontal displacement of the structure crest without considering the influence of soil then takes the form:
= YH,
(2)
where y — shear strain developed in the analyzed domain due to earthquake, H — structure height.
Soil-structure interaction can be taken into account by introducing the soil-structure stiffness ration as:
f = g w
= k h '
(3)
where KS — structure stiffness, W — structure width.
When addressing structures with low overlying it is necessary to unify results provided by the free-field method and by traditional solutions adopting the influence of earth pressures caused by earthquake according to, e. g. Mononobe — Okabe [8]. One option is to introduce the flexibility ratio R:
R =
4(1 -vm ) F
3 - 4v + F'
(4)
where vm — Poisson's ratio of the soil, F — soil-structure stiffness ratio.
The horizontal displacement at the structure crest for such conditions then follows from:
AS = R ■ Adiff ■
(5)
In case of layered subsoil one may employ the relations pertinent to the homogeneous environment by adopting an effective shear modulus G
hom'
Two
options are available. Exploiting the Voigt assumption of constant strain in all layers the effective shear modulus Ghom is determined as an arithmetic weighted average of shear moduli of individual layers. The average shear strain is determined analogically, i.e. as a weighted average of shear strains developed in individual layers. Satisfying the assumption of piecewise constant shear strain in individual layers requires the values of strains of individual layers to comply with stiffness ratio of individual soils to satisfy stress continuity across the interface of the two materials. Assumption of constant stresses through all layers then yields the Reuss estimate of the shear stiffness as a weighted average of compliances:
Ghom ( h G1 + h G2 '
(6)
where h1 — thickness of layer 1, h2 — thickness of layer 2, G1 — shear modulus of layer 1, G2 — shear modulus of layer 2, h — sum of thicknesses of individual layers.
Pseudo-static analysis using fem
This approach allows for considering a nonhomogeneous layered subsoil and a locally nonlinear response of soil. The maximum shear strain in individual layers follows from the free-field analysis. These are introduced in the form of displacements prescribed along vertical boundaries of the analyzed domain. In case of layered subsoil we consider piecewise constant variation of strains. The analysis is then carried out with the help of finite element method similarly to standard static tasks. See [9] for details.
Fully dynamic analysis using fem
The numerical solution of dynamic load effects is based on the D'Alembert principle, which states that the sum of all forces acting on the body is equal to zero. Applying this principle to a problem with many degrees of freedom with domain loaded by vertically propagating waves we obtain upon discretization and with reference to Figure 2 the system of second order differential equations in the form:
MuR (t) + CÙ r (t)9 + KuR (t) = -Mu0(t) + CÙ10(t) \y=o,
(7)
cases where the geometrical domain has to be truncated in the vertical direction within a given layer. The load is then introduced in terms of the incoming wave u0(t) = uI0(t) , recall Fig. 2.
Providing the domain can be truncated at the interface between compliant and infinitely rigid layers, the absorbing boundary condition can be replaced by the so-called fixed boundary. Because u0(t) then corresponds to the total motion we have uR(x,y = 0,t) = 0 and the second term on the right-hand side of Eq. (7) drops out. This boundary condition should, however, be used with caution as it may greatly overestimate the strains developed in the domain. Given this, only the model with absorbing boundary will be considered henceforth. This will apply also to the calculation of free-field strains adopted in the pseudo-static analysis [10].
As seen in Fig. 2 the absorbing boundary is imagined in the form of viscous dampers where the pressure and shear wave velocities cs and cp are provided by:
cp =
; Cs =
(9)
where M, C and K are the mass, damping and stiffness matrices, respectively.
Equation (7) assumes that the total displacement vector u(x,t) can be split into the displacement u0(t) which is constant within the whole domain, thus independent of the position x, and the relative displacement uR(x,t) such that:
u( x, t) = Ùo(t) + ùr ( x, t). (8)
This particular form of Eq. (7) assumes the so called absorbing or quiet boundary conditions applied at the bottom boundary of the domain. This condition allows for absorbing the downward travelling wave, so it is not reflected back to the domain. Such boundary is typically introduced in
where Eoed, G, p are the oedometric modulus, shear modulus and density of the soil, respectively. For particular derivation of these boundary conditions we refer the interested reader to Zienkiewicz et al. (1986).
Apart from the bottom boundary we should also accord our attention to the boundary conditions applied along the lateral edges. As pointed out in Zienkiewicz et al. [11] such boundary conditions should account for the disturbance from the free-field conditions caused by the presence of structure, e.g. tunnel. This is schematically shown in Figure 2, suggesting that only the difference (uR - uR) between the actual uR and free-field motion uR should be absorbed. The free-
field motion in particular is derived from an independent one-dimensional free-field column analysis assuming the same loading conditions as in case of the two-dimensional analysis [12]. To ensure that the solution of the 1D free-field column analysis is properly reproduced by the 2D analysis in case of no structure, it is necessary to supplement these boundary conditions by the prescribed shear stress provided by the 1D column analysis, see [13], for additional details. Thus, in case of no material damping Eq. (7) becomes:
MuR(t) + C +KuR(t)
R (t)
BBuc
(t)
y=o
+CLBÙr
(t)|
-Mu I0(t) + CBB ' '
u„
+CLBu FF'
-RAt)\
(t)|y=0
0 +R,(t)lx=
(10)
The damping matrices CBB and CLB reflect the effect of artificial viscous dampers applied along the absorbing boundaries (the superscripts BB and LB identify the bottom and lateral boundaries). Finally, the vector RT stores the nodal forces associated with the prescribed boundary shear stress.
Case study
The analyzed structure is a circular track tunnel in Baku, the capital city of
Azerbaijan. The section is 1 500 m long with two directional curves of 1 200 m and 1 600 m radius. The track tunnel crosses the height difference of 25,65 m and has a slope of 1,90 %. Tunnel excavation was carried out using a TBM with a diameter of boring head of 6 m. The tunnel lining consists of reinforced concrete segments with a thickness of 0,3 m. The geological conditions in this area are: under the thin layer of sediments with sand and sandstones deposits there is a subsoil formed by clay layers (medium and low plasticity) and a rock mass of limestone.
The computational model was simplified considering only clay (top) and limestone (bottom) layers — see Fig. 3. In numerical analysis the limestone was assumed to behave linearly elastic while the response of the clayey layer is governed by the Mohr-Coulomb model. The adopted material parameters are listed in Tables 1 and 2. As evident from Table 2 the MC model considered the modulus of elasticity varying with depth. Upon unloading during excavation step the material stiffness is governed by the unloading/reloading modulus Eur. In case of dynamic analysis, soils having in the original static analysis the elastic modulus E
cpp(ûR-ù/F)\ csp(vR-vRFF)^i
[¿cpp(ùR-ùRFF)
csp(vR - v,0FF)
fcpp(ùR-û/F) 8
3
\^csp(vR-vRFF)^ w
csp(vR-vI0 )
^xy
Ur-ÙrFF)
[¿cpp(ùR-ûRFF)
csp(vR - vwFF)
u (x,t) = uI0 (t) + UR (x,t) Fig. 2. Boundary and loading conditions for both free-field and 2D analysis
Fig. 3. Geometry and finite element mesh
Table 1
Material parameters of limestone — linear-elastic material model
Parameter symbol value unit
Unit weight Y 21,70 kN/m3
Poisson number V 0,30 -
Elastic modulus E 630,00 MPa
Biot number a 1,00 -
Coefficient of earth Ko 0,60 -
pressure at rest
Unit weight of satu- y sat 22,00 kN/m3
rated soil
Table 2
Material parameters of clay — Mohr-Coulomb material model
Parameter symbol value unit
Unit weight Y 19,60 kN/m3
Poisson number V 0,40 -
Elastic modulus on the surface E 5,00 MPa
Increase in elastic modulus with depth Kd 0,51 MPa/m
Excavation modulus E,,r 80,00 MPa
Parameter symbol value unit
Dynamic modulus Edyn 80,00 MPa
Biot number a 1,00 -
Coefficient of earth K> 0,70 -
pressure at rest
Angle of internal Vef 25,00 o
friction
Cohesion cef 10,00 kPa
Angle of dilatancy V 0,00 O
Unit weight of Ysat 20,60 kN/m3
saturated soil
evolving with depth should be replaced with the soil having constant properties in the given layer. Also, the value of elastic modulus must be increased to ensure that the velocity of a propagating wave corresponds to that observed in situ, i. e. the dynamic modulus of elasticity Edyn must be considered. To meet these requirements, the value of Eur, set equal to the dynamic modulus Edyn, was assigned to the whole clayey layer.
Fig. 4 plots the acceleration and velocity profiles are derived from the real earthquake recorded near the area of the tunnel construction. The record was modified so that it corresponds to the requirements of the Azerbaijani design standards and specifications for seismic design of building structures. Design acceleration was set to 0,25g and time integration step to 0,01 s. Thus modified record was adopted as dynamic loading to the computational model.
The calculation is performed in GEO5 FEM program using fully dynamic analysis. The analysis itself has four subsequent computational phases. The first three computational phases solve the static part of the analysis:
1. Computation of the initial geostatic stress.
2. Tunnel excavation and activation of 30 % of the load on the unreinforced excavation.
3. Installation of reinforced concrete segment lining and activation of the remaining 70 % of the load.
While the last, fourth, phase solves the dynamic part of the analysis. The computational model is loaded by the earthquake and additional stress contribution caused by the seismic load in being computed in this phase.
The distribution of maximum values of internal forces around the perimeter of tunnel lining comparing the results at the end of excavation (third computational phase) with values of the initial forces achieved during earthquake (fourth computational phase) appears in Figure 5 — 7. The initial forces pertinent to the end of excavation plotted in grey are superposed by absolute maximal values achieved during earthquake in red.
The fully dynamic 2D analysis was carried in the light of Eq. (10). The results of the analysis were compared with other methods to suggest an optimal boundary conditions (at least in the light of bending moments,) and as such to promote more simple methods at the preliminary stage of design. Comparing the results clearly shows that using fixed boundary conditions is inappropriate as it provides unrealistically large values of internal forces. The reason is that the rock mass can hardly be considered as infinitely stiff in comparison with the stiffness of a
Fig. 5. Maximal normal forces in the tunnel lining
Fig. 6. Maximal shear force in the tunnel lining
tp [rad]
Fig. 7. Maximal bending moments in the tunnel lining 88
soil layer. In general, it can be stated that employing absorbing boundary conditions on the bottom boundary appears more suitable and the fixed boundary conditions should be used with caution [14, 15]
Conclusions
In general, two groups of approaches are available to analyze underground structures subjected to earthquake — the pseudo-static analysis and fully dynamic analysis using FEM. In the application of fully dynamic analysis in GEO5 FEM program there must be
REFERENCES
accorded an attention to implementation of various types of boundary conditions. While along vertical boundaries the combination of static and free-field boundary condition is the only option, the bottom boundary allows for choosing either fixed or absorbing boundary conditions. In this regard, the absorbing boundary condition is given the preference. Although more complex, this approach to dynamic analysis is suitable for engineering practice and allows us to account for mutual soil-structure interaction effects.
1. Wang J. Seismic design of tunnels. A simple state-of-the-art design approach. Monograph 7, Parsons, Brinckerhoff Quade and Douglas Inc., 1993.
2. Hashash Y. M. A. Seismic design and analysis of underground structures. Tunneling and undergrounding Space Technology, 2001, Vol. 16, pp. 247—293. D0l:10.1016/S0886 — 7798(01)00051-7.
3. Zienkiewicz O. C., Bicanic N., Shen F. Q. Generalized smith boundary — a transmitting boundary for dynamic computation. Institute for Numerical Methods in Engineering, University College of Swansea, 1986, Vol. 207, pp. 15—26.
4. ITA AITES home page, available at http://www.ita-aites.org/en, 2020.
5. Power M. S., Rosid D., Kaneshiro J. Strawman: screening, evaluation, and retrofit design of tunnels. Report Draft, volume Ill. New York: National Center for Earthquake Engineering Research, 1996, 354 p.
6. Kravcov A. N., Svoboda P., Pospichal V., Morozov D. V., Ivanov P. N. Limit depth of rock mine shafts for underground shelters. Key Engineering Material, 2017, Vol. 755, pp. 198—201.
7. Kravtsov A., Svoboda P., Pospichal V., Zdebsky J. Experimental Studies on Process of Transition of Explosion to Deflagration due to Methane Gas Explosion in Underground Structures. International Conference on Military Technologies 2015, Rio de Janeiro. IEEE Institute of Electrical and Electronics Engineers Inc., 2015, pp. 125 — 133.
8. Towhata I. Geotechnical Earthquake Engineering (Springer Series in Geomechanics and Geoengineering). Springer Berlin Heidelberg, 2008, 684 p.
9. Poklopova T., Pavelcova V., Janda T., Sejnoha M. Evaluation of real underground structure subjected to earthquake — pseudostatic analysis. Acta Polytechnica CTU Proceedings, 2018, Vol. 15, pp. 91—93. D0I:10.14311/APP.2018.15.0088.
10. Pavelcova V., Poklopova T., Janda T., Sejnoha M. The influence of boundary conditions on the response of underground structures subjected to earthquake. Acta Polytechnica CTU Proceedings, 2018, Vol. 26, pp.64—79. DOI: 10.14311/APP.2020.26.0064.
11. Zienkiewicz O. C., Bicanic N., Shen F. Q. Generalized smith boundary — a transmitting boundary for dynamic computation. Institute for Numerical Methods in Engineering, University College of Swansea, 1986, Vol. 207, pp. 109 — 138.
12. Kucera D. Analysis of geotechnical structures subjected to earthquake. Diploma thesis. CTU in Prague, 2017, 15 p.
13. Kravtsov A., Svoboda P. Experimental Studies of the Blast Pressure due to an Explosion in the Tunnel. Proceedings from the Fourth International Symposium on Tunnel Safety and Security. Stockholm: Royal Institute of Technology, 2014. pp. 281 — 288.
14. Pruska J., Pavelcova V., Poklopova T., Janda T., Sejnoha M. The Response of Underground Structures on the Seismic Loadings. Proceedings of the WTC 2019 ITA-AITES World Tunnel Congress (WTC 2019). Naples, Italy, 2019, 10 p.
15. Sejnoha M., Janda T., Pruska J., Broucek M. Metoda konecnych prvku v geomechanice: Teoretické zâklady a Inzenyrské aplikace. Prague, 2015, 314 p. ü^re
ИНФОРМАЦИЯ ОБ АВТОРАХ
Ян Прушка1 — доцент;
Вероника Павельцова1 — РЬЮ^^еп^
1 Чешский технический университет в Праге, факультет гражданского строительства, кафедра Геотехники, Такурова 7, 166 29 Прага 6 — Дейвице, Чешская Республика.
INFORMATION ABOUT THE AUTHORS
Jan Pruska1, Assistant professor, Pruska@fsv.cvut.cz; Veronika Pavelcova1, PhD-student;
1 Czech Technical University in Prague, Faculty of Civil Engineering, Department of Geo-technics, Thakurova 7, 166 29 Prague 6 - Dejvice, Czech Republic.
Получена редакцией 20.01.2021; получена после рецензии 17.02.2021; принята к печати 10.03.2021. Received by the editors 20.01.2021; received after the review 17.02.2021; accepted for printing 10.03.2021.
