American Journal of Interdisciplinary Research and Development ISSN Online: 2771-8948 Website: www.ajird.journalspark.org Volume 13, Feb., 2023 124 | P a g e NUMERICAL INVESTIGATION OF INTERACTION OF UNDERGROUND STRUCTURE (TUNNELS) WITH ELASTIC WAVES PROPAGATED IN GROUND MEDIUM Gaynazarov Sultan National University of Uzbekistan named after Mirzo Ulugbek, dr. phys. - math., Professor of the Faculty «Applied Mathematics and Intellectual Technologies» gaynazarovsm@mail.ru Rakhimdjanova Nodira National University of Uzbekistan named after Mirzo Ulugbek, 2nd year master student of the Faculty of «Applied Mathematics and Intellectual Technologies» nodiraziyo1990@gmail.ru ABSTRACT In this article considers a dynamic method for calculating the analysis of an underground structure (tunnels) interacting with the soil environment. On the basis of numerical simulation, a diffraction pattern for an underground structure (tunnels) was constructed. Keywords: finite element method, Newmark method, Bate method, stiffness matrix, damping matrix, mass matrix. Introduction Seismic waves - vibrations of rocks in the Earth, resulting from natural (earthquakes) or artificial processes of their excitation. At present, the study of seismic waves generated by earthquakes is necessary for understanding the nature of earthquakes and their prediction. Seismic waves caused by earthquakes or explosions elastic waves propagating in the body of the Earth whether. Of the body waves, the primary, or P, wave has the higher speed of propagation and so reaches a seismic recording station faster than the secondary, or S, wave. P waves, also called compressional or longitudinal waves, give the transmitting medium whether liquid, solid, or gas a back and forth motion in the direction of the path of propagation, thus stretching or compressing the medium as the wave passes any one point in a manner similar to that of sound waves in air. Seismic action is a special term, which in the practice of calculating structures for seismic resistance means the oscillatory movement of the soil during an earthquake, which creates kinematic excitation of vibrations of building structures. https://www.britannica.com/science/longitudinal-wave https://www.britannica.com/science/wave-water https://www.merriam-webster.com/dictionary/propagation https://www.britannica.com/science/secondary-wave https://www.britannica.com/science/longitudinal-wave https://www.britannica.com/science/motion-mechanics https://www.britannica.com/science/sound-physics American Journal of Interdisciplinary Research and Development ISSN Online: 2771-8948 Website: www.ajird.journalspark.org Volume 13, Feb., 2023 125 | P a g e Literature Review O. K. Zenkevich, The finite element method in engineering science. The monograph is devoted to the presentation of the fundamentals of the finite element method – one of the most effective modern methods for the numerical solution of engineering, physical and mathematical problems using computers. L. Segerlind, Application of the finite element method. The monograph reflects the work of many researchers. The order in which the material is arranged depends on the results of the author's experience. K. Bathe, Finite element procedures. Finite element procedures are now an important and frequently indispensable part of engineering analysis and design. Finite element computer programs are now widely used in practically all branches of engineering for the analysis of structures, solids, and fluids. Formulation of the problem. The movement of an underground structure (tunnel) and its surrounding elastic medium (soil) is caused by the propagation of a longitudinal seismic wave in the medium in the area 𝛺 = {(𝑥, 𝑦): 𝑥 ∈ [0, 𝐿], 𝑦 ∈ [0, 𝐻]}, and with external borders 𝛤𝑙 = {(𝑥, 𝑦): 𝑥 = 0, 𝑦 ∈ [0, 𝐻]}, 𝛤𝑟 = {(𝑥, 𝑦): 𝑥 = 𝐿, 𝑦 ∈ [0, 𝐻]}, 𝛤𝑡 = {(𝑥, 𝑦): 𝑥 ∈ [0, 𝐿], 𝑦 = 𝐻}, 𝛤𝑏 = {(𝑥, 𝑦): 𝑥 ∈ [0, 𝐿], 𝑦 = 0}, as well as with the boundaries of an underground structure (tunnel) of a rectangular area 𝛤𝑙𝑖𝑛 = {(𝑥, 𝑦): 𝑥 = 𝐿1, 𝑦 ∈ [0, 𝐻1 + 𝐻2]}, 𝛤𝑟𝑖𝑛 = {(𝑥, 𝑦): 𝑥 = [𝐿1 + 𝐿2], 𝑦 ∈ [0, 𝐻1 + 𝐻2]}, 𝛤𝑡𝑖𝑛 = {(𝑥, 𝑦): 𝑥 ∈ [0, 𝐿1 + 𝐿2], 𝑦 = [𝐻1 + 𝐻2]}, 𝛤𝑏𝑖𝑛 = {(𝑥, 𝑦): 𝑥 ∈ [0, 𝐿1 + 𝐿2], 𝑦 = [0, 𝐻1]}, The mathematical model in a strict formulation (differential equation) has the form: 𝜌�̈� + 𝑐�̇� − 𝑑𝑖𝑣(𝜎) = 𝑟 in area Ω (1) 𝜎𝑥𝑥 = 𝑐𝑝𝜌�̇�, 𝜎𝑦𝑦 = 𝑐𝑠𝜌�̇�, 𝜎𝑥𝑦 = 0 at the borders Γt and Γb (2) 𝜎𝑥𝑥 = 2𝑐𝑝𝜌𝑓(𝑡)𝑛𝑥, 𝜎𝑦𝑦 = 2𝑐𝑝𝜌𝑓(𝑡)𝑛𝑦, 𝜎𝑥𝑦 = 0 at the border Γl (3) 𝜎𝑥𝑥 = 𝑐𝑝𝜌�̇�, 𝜎𝑦𝑦 = 𝑐𝑠𝜌�̇�, 𝜎𝑥𝑦 = 0 at the border Γr (4) 𝑢 = 0, �̇� = 0 при t=0 (5) Given the Koshi relation 휀𝑥𝑥 = 𝜕𝑢𝑥/𝜕𝑥; 휀𝑦𝑦 = 𝜕𝑢𝑦/𝜕𝑦; 휀𝑥𝑦 = 𝜕𝑢𝑥/𝜕𝑦 + 𝜕𝑢𝑦/𝜕𝑥 (6) and Hooke's law 𝜎 = 𝐵휀, 𝐵 = [ 𝜆 + 2𝜇 𝜆 0 𝜆 𝜆 + 2𝜇 0 0 0 𝜇 ] (7) 𝜆 and 𝜇 Lame coefficients. American Journal of Interdisciplinary Research and Development ISSN Online: 2771-8948 Website: www.ajird.journalspark.org Volume 13, Feb., 2023 126 | P a g e In this problem, a longitudinal wave falls on the left edge at a given angle, with an initial impulse. Conditions for the transparency of borders have been established on all sides. Since the problem is solved by the finite element method, which involves setting a weak form (variational form) of the equilibrium equations, we will use the Lagrange variational principle, as well as the d'Alembert principle. As a result, the variational equilibrium equation for an isotropic body has the form: ∫ 𝜆𝑑𝑖𝑣(𝑢)𝑑𝑖𝑣(𝛿𝑢) + 2𝜇휀𝑇휀 − (𝑟 − 𝛺 𝜌�̈� − 𝑐�̇�)𝛿𝑢𝑑𝛺 − ∫ 𝑃 𝛤 𝛿𝑢𝑑𝛤 = 0 (8) where 𝑟 – volume load vector, 𝑃 – surface load vector. Applying the method of partial discretization by finite elements i.e. movements in the element, representing in the form: 𝑢𝑒 = ∑ 𝑢𝑖(𝑡)𝜑(𝑥, 𝑦)𝑖 𝑛 𝑖=1 (9) where 𝑛 – the number of degrees of freedom in the element 𝑒, 𝑢𝑖(𝑡) – desired nodal displacements depending on time, 𝜑(𝑥, 𝑦)𝑖 – form functions on the element. We obtain a system of linear ordinary equations: 𝑀�̈� + 𝐶�̇� + 𝐾𝑈 = 𝑅 (10) With initial conditions 𝑈 = 0, �̇� = 0 at t=0 (11) Equations (11) obtained from the consideration of static equilibrium at time t can be written as: 𝐹𝑖(𝑡) + 𝐹𝑑(𝑡) + 𝐹𝑒(𝑡) = 𝑅(𝑡) (12) To solve systems of equations (12) with initial conditions (11), the implicit difference methods of Newmark and Bate are used. Step-by-step solution using Newmark integration method. A. Initial calculations: 1. Form stiffness matrix 𝐾, mass matrix 𝑀 and damping constant 𝐶. 2. Initialize 𝑈0, �̇�0 and �̈�0. 3. Select time step ∆𝑡, and parameters 𝛼 and 𝛿 and calculate integration constants: 𝛿 ≥ 0,50; 𝛼 ≥ 0,25(0,5 + 𝛿)2; 𝑎0 = 1 𝛼∆𝑡2 ; 𝑎1 = 𝛿 𝛼∆𝑡 ; 𝑎2 = 1 𝛼∆𝑡 ; 𝑎3 = 1 2𝛼 − 1 ; 𝑎4 = 𝛿 𝛼 − 1 ; 𝑎5 = ∆𝑡 2 ( 𝛿 𝛼 − 2) ; 𝑎6 = ∆𝑡(1 − 𝛿) ; 𝑎7 = 𝛿∆𝑡. 4. Form effective stiffness matrix �̂�: �̂� = 𝐾 + 𝑎0𝑀 + 𝑎1𝐶. 5. Triangularize �̂�: �̂� = 𝐿𝐷𝐿𝑡. American Journal of Interdisciplinary Research and Development ISSN Online: 2771-8948 Website: www.ajird.journalspark.org Volume 13, Feb., 2023 127 | P a g e B. For each time step: 1. Calculate effective loads at time 𝑡 + ∆𝑡 : �̂�𝑡+∆𝑡 = 𝑅𝑡+∆𝑡 + 𝑀(𝑎0𝑈𝑡 + 𝑎2�̇�𝑡 + 𝑎3�̈�𝑡) + 𝐶(𝑎1𝑈𝑡 + 𝑎4�̇�𝑡 + 𝑎5�̈�𝑡). 2. Solve for displacements at time 𝑡 + ∆𝑡 : 𝐿𝐷𝐿𝑡𝑈𝑡+∆𝑡 = �̂�𝑡+∆𝑡. 3. Calculate accelerations and velocities at time 𝑡 + ∆𝑡 : �̈�𝑡+∆𝑡 = 𝑎0(𝑈𝑡+∆𝑡 − 𝑈𝑡) − 𝑎2�̇�𝑡 − 𝑎3�̈�𝑡; �̇�𝑡+∆𝑡 = �̇�𝑡 + 𝑎6�̈�𝑡 + 𝑎7�̈�𝑡+∆𝑡 . Step-by-step solution using the Bathe integration method. A. Initial calculations: 1. Form stiffness matrix 𝐾, mass matrix 𝑀 and damping constant 𝐶. 2. Initialize 𝑈0, �̇�0 and �̈�0. 3. Select time step ∆𝑡 and calculate integration constants: 𝑎0 = 16 ∆𝑡2 ; 𝑎1 = 4 ∆𝑡 ; 𝑎2 = 9 ∆𝑡2 ; 𝑎3 = 3 ∆𝑡 ; 𝑎4 = 2𝑎1 ; 𝑎5 = 12 ∆𝑡2 ; 𝑎6 = − 3 ∆𝑡2 ; 𝑎7 = − 1 ∆𝑡 . 4. Form effective stiffness matrices 𝐾1̂ 𝑎𝑛𝑑 𝐾2̂ : �̂�1 = 𝐾 + 𝑎0𝑀 + 𝑎1𝐶, �̂�2 = 𝐾 + 𝑎2𝑀 + 𝑎3𝐶. 5. Triangularize 𝐾1̂ 𝑎𝑛𝑑 𝐾2̂: �̂�1 = 𝐿1𝐷1𝐿1 𝑡, �̂�2 = 𝐿2𝐷2𝐿2 𝑡. B. For each time step: First sub-step: 1. Calculate effective loads at time 𝑡 + ∆𝑡 2 : �̂� 𝑡+ ∆𝑡 2 = 𝑅 𝑡+ ∆𝑡 2 + 𝑀(𝑎0𝑈𝑡 + 𝑎4�̇�𝑡 + �̈�𝑡) + 𝐶(𝑎1𝑈𝑡 + �̇�𝑡). 2. Solve for displacements at time 𝑡 + ∆𝑡 2 : 𝐿1𝐷1𝐿1 𝑡𝑈 𝑡+ ∆𝑡 2 = �̂� 𝑡+ ∆𝑡 2 . 3. Calculate accelerations and velocities at time 𝑡 + ∆𝑡 2 : �̈� 𝑡+ ∆𝑡 2 = 𝑎1 (�̇� 𝑡+ ∆𝑡 2 − �̇�𝑡) − �̈�𝑡; �̇� 𝑡+ ∆𝑡 2 = 𝑎1 (𝑈 𝑡+ ∆𝑡 2 − 𝑈𝑡) − �̇�𝑡 . Second sub-step: 1. Calculate effective loads at time 𝑡 + ∆𝑡 : �̂�𝑡+∆𝑡 = 𝑅𝑡+∆𝑡 + 𝑀 (𝑎5𝑈 𝑡+ ∆𝑡 2 + 𝑎6𝑈𝑡 + 𝑎1�̇� 𝑡+ ∆𝑡 2 + 𝑎7�̇�𝑡) + 𝐶(𝑎1𝑈𝑡+∆𝑡/2 + 𝑎7𝑈𝑡). 2. Solve for displacements at time 𝑡 + ∆𝑡 : American Journal of Interdisciplinary Research and Development ISSN Online: 2771-8948 Website: www.ajird.journalspark.org Volume 13, Feb., 2023 128 | P a g e 𝐿2𝐷2𝐿2 𝑡𝑈𝑡+∆𝑡 = �̂�𝑡+∆𝑡. 3. Calculate accelerations and velocities at time 𝑡 + ∆𝑡 : �̈�𝑡+∆𝑡 = −𝑎7�̇�𝑡 − 𝑎1�̇� 𝑡+ ∆𝑡 2 + 𝑎3�̇�𝑡+∆𝑡 , �̇�𝑡+∆𝑡 = −𝑎7𝑈𝑡 − 𝑎1𝑈 𝑡+ ∆𝑡 2 + 𝑎3𝑈𝑡+∆𝑡 . Analysis and results. The study was conducted on the following problem: Meaning L=400m, H=100m, Young's modulus E=1.e8Pa, Poisson's ratio ν=0.3, ))2-(1)+/((1E)),1(2/(  =+= E - Lame coefficients, density ρ=1800Pa, coefficient of the longitudinal and transverse waves of the medium were calculated by the formula  /)2(,/ +== sp cc . A triangular element with six nodes, the second degree of accuracy, was used as finite elements. A study was carried out on the influence of the number of elements on the accuracy of the solution, the following table shows the value of the displacements in the center of the two-dimensional region at the point (200,50), the calculation was carried out by the Newmark method, with a step dt=0.001sec at a point in time t=1sec: It can be seen from the table that with an increase in the partition, the accuracy of the result increases and when the sides are divided into 100/25 we have 4 correct digits. Next, a comparative analysis of the calculations by the Newmark method and the Bathe method was carried out. The following table shows the offset values ux at the point (200,50) at a point in time t=1 sec. with the same physical parameters as above. It can be seen from the table that at a step 0.0001 Newmark's method yields an exact solution, while Bathe's method yields an exact solution with a step 0.001. Thus, the Bathe method produces a more accurate result with a larger step. Next, consider the oscillation of the midpoint of a two-dimensional body with the above physical and geometric characteristics. The figure shows the movement of a point for 4 seconds. Amount of elements\Bias 212(20/5) 858(40/10) 1896(60/15) 5306(100/25) ux -0.00728348 -0.00739348 -0.00742446 -0.00743622 uy 0.00130321 0.00126898 0.00129406 0.00129469 Step\ Method 0.05 0.01 0.005 0.001 0.0001 Newmark -0.006348 -0.007408 -0.007421 -0.007436 -0.007424 Bathe -0.006994 -0.007417 -0.007421 -0.007424 -0.007424 American Journal of Interdisciplinary Research and Development ISSN Online: 2771-8948 Website: www.ajird.journalspark.org Volume 13, Feb., 2023 129 | P a g e The dynamics of wave propagation along the line y=50 is shown in the following figure Conclusion Studies have shown that for dynamic calculation it is not necessary to thicken the mesh, the use of the Bate method is preferable to the use of the Newmark method, the use of the boundary transparency condition correctly simulates the passage of the wave. References 1. Сегерлинд Л. Применение метода конечных элементов – Москва.: «МИР», 1979. – 392с. [Segerlind, L. 1979. Application of the finite element method. Moscow. Mir. 392p.] (in Russian). 2. Зенкевич, О. Метод конечных элементов в технике. – Москва: «МИР», 1975. – 541с. [Zienkiewicz, O. 1975. The finite element method in engineering science. Moscow. Mir. 541p.] (in Russian). 3. K. Bathe, 2014. Finite element procedures. Second edition. Watertown. 1043p.