Seismic wave propagation numerical simulation method and system suitable for tundra around tunnel

Abstract

The invention discloses a seismic wave propagation numerical simulation method and a seismic wave propagation numerical simulation system suitable for a tunnel peripheral frozen soil layer, and relates to the field of seismic wave simulation propagation. According to the method, a layered model is constructed to describe physical property distribution of the frozen soil layer at different depths, and a viscoelastic model is combined to integrate the layered model into a finite element model of the tunnel peripheral frozen soil layer; dividing the finite element model by adopting an initialized uniform grid, then preliminarily simulating a seismic wave propagation process, and by comparing the spatial change rate of the wave velocity, identifying a region with a remarkably changed wave velocity gradient, and extracting the wave velocity change amplitude and directivity data to determine the boundary and shape of the region; and planning the density distribution of the updated grid in combination with the geological parameters of the gradient change region. According to the method, the boundary and internal characteristics of the wave velocity gradient region are identified, the grid shape and density are planned according to the spatial change rate and distribution characteristics of the gradient, it is ensured that the calculation region can reflect the multi-scale characteristics of the frozen soil layer, and the seismic wave propagation path and the numerical calculation of energy distribution can be optimized.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic wave simulation and propagation, and in particular to a numerical simulation method and system for seismic wave propagation in frozen soil layers surrounding tunnels. Background Art

[0002] During the construction of plateau tunnels, the permafrost surrounding the tunnels is often high-altitude permafrost. High-altitude permafrost is a unique geological environment composed of three phases: ice, water, and soil. Its interior exhibits complex physical properties, such as oriented ice crystals, nonuniform distribution of multiphase media, and dynamic changes in temperature and moisture content. These characteristics cause seismic waves to exhibit significant anisotropy and dispersion during propagation, with particularly prominent nonlinear energy attenuation. Traditional wave propagation models, often based on homogeneous assumptions, struggle to fully reflect the complexity of high-altitude permafrost, especially under the coupled effects of temperature changes, ice-water phase transitions, and changes in moisture content. Existing methods for studying its dynamic response are limited.

[0003] When performing numerical simulations of seismic wave propagation, existing technologies often use uniform meshing for calculations. Although this method can simply discretize the computational domain, when dealing with regions of permafrost with significant anisotropy and dispersion characteristics, the uniform mesh cannot be optimized for areas with significant wave velocity changes. Specifically, for areas of moisture redistribution in such frozen soils, the significant anisotropy and dispersion characteristics cause significant changes in wave velocity gradients and increased energy attenuation. The resolution of traditional uniform meshing is insufficient to capture the rapidly changing characteristics of wave velocity, resulting in deviations in simulation results. Summary of the Invention

[0004] 1) Technical issues solved The present invention provides a numerical simulation method and system for seismic wave propagation in permafrost layers around tunnels, which can capture the wave velocity variation characteristics and energy attenuation laws of high-altitude permafrost layers based on adaptive grid division.

[0005] 2) Technical solution To achieve the above-mentioned object, the present invention provides the following technical solution: a numerical simulation method for seismic wave propagation in the frozen soil layer around a tunnel, comprising: Obtaining geological parameters and environmental parameters of the permafrost layer surrounding the selected tunnel and the tunnel's design magnitude conditions, and setting loading conditions for simulated seismic waves according to the design magnitude conditions, including the frequency, magnitude, and location of the seismic wave source; Based on the environmental parameters of the permafrost layer, a layered model is constructed to describe the distribution of physical properties of the permafrost layer at different depths. The layered range is adjusted in combination with the ground temperature field and the freezing depth. The geological distribution of each layer in the layered model is used as an input parameter of the viscoelastic model, which is integrated into a finite element model of the permafrost layer around the tunnel to simulate the dispersion and energy attenuation characteristics of the permafrost layer at different depths. The region of the finite element model is initially divided using an initialized uniform grid, and the region is discretized into a regular grid shape. According to the set seismic wave loading conditions, the seismic wave propagation process is simulated in a finite element model, the wave velocity distribution in the frozen soil is calculated, and by comparing the spatial change rate of the wave velocity, the area where the wave velocity gradient changes significantly is identified and marked as a wave velocity gradient area. The wave velocity change amplitude and directionality data of the wave velocity gradient area are extracted to determine its boundary and shape; In combination with the shape and gradient direction of the velocity gradient zone, the boundary curve is calculated based on the geometric form and the shape of the updated grid is generated, and the density distribution of the updated grid is planned in combination with the geological parameters of the gradient change area; The finite element model is updated and the grid is divided. The seismic wave propagation is re-simulated according to the set seismic wave loading conditions. The finite difference method is used for numerical calculation, and the waveform data and energy changes in each wave velocity gradient area are analyzed.

[0006] Furthermore, the geological parameters of the permafrost layer obtained include density, ice content, porosity and viscoelastic parameters, and the thermodynamic characteristics of the permafrost layer are analyzed by combining drilling, geophysical detection and temperature monitoring methods; The geometric parameters, support structure type and construction method of the tunnel, including the tunnel cross-sectional shape, lining thickness and support material properties, are obtained to determine the design magnitude condition.

[0007] Furthermore, the environmental parameters include the geothermal field and freezing depth of the permafrost layer. By calculating the temperature distribution of the permafrost layers at different depths and combining the geological parameters, a stratified model of the permafrost layer is constructed.

[0008] Furthermore, based on the geological parameters of the permafrost layer, spatial layers of surface permafrost, high ice content layer and deep permafrost are divided, and the layered model is used to describe the spatial distribution of physical properties of the permafrost layer.

[0009] Furthermore, the viscoelastic model is used to simulate the seismic wave propagation characteristics of each of the spatial layers, wherein the seismic wave velocity at different frequencies is analyzed by Fourier transform, and the energy attenuation characteristics of each spatial layer are characterized in combination with the quality factor. By adjusting the damping parameters, the simulation results can reflect the wave attenuation characteristics of the permafrost layer.

[0010] Furthermore, the layered model provides the physical property distribution of the permafrost layer at different depths, and the viscoelastic model describes the dispersion and energy attenuation characteristics of the permafrost material. The physical property distribution of each layer in the layered model is used as the input parameter of the viscoelastic model and integrated into a finite element model of the permafrost layer around the tunnel, which is used to simulate the propagation characteristics of seismic waves at different depths in the permafrost layer and their energy changes. An absorption boundary is set at the boundary of the calculation domain of the finite element model to prevent the artificial truncation of the calculation domain from causing seismic wave reflection and transmission.

[0011] Furthermore, the loading conditions for simulating seismic waves are set according to the design magnitude conditions. Specifically, a static ground stress field of the frozen soil layer is constructed, and a dynamic seismic load is applied to simulate the stress response and corresponding failure mode under different earthquake intensities. Correspondingly: Based on the geological parameters of the permafrost layer and the geometric parameters of the tunnel, a static geostress field is constructed to reflect the stress distribution of the tunnel and the surrounding soil in its natural state. Apply dynamic seismic loads according to the set seismic wave loading conditions, set the intensity and frequency of dynamic seismic waves according to different seismic intensities, and simulate the vibration of the tunnel structure caused by the earthquake; Set radial, tangential, and axial ultimate load conditions to simulate the bearing capacity of the tunnel lining and surrounding rock in different directions and evaluate the stress response of the structure in different directions.

[0012] Furthermore, the finite difference method is used to simulate the propagation behavior of seismic waves in permafrost. The region is discretized into individual grid points through an initialized uniform grid, and the finite difference formula is used to solve the energy propagation attenuation characteristics of seismic waves between each grid point.

[0013] Furthermore, the velocity gradient area marked specifically defines the portion where the velocity variation amplitude exceeds a set threshold as the velocity gradient area, and extracts the velocity variation amplitude and directionality data of the velocity gradient area to determine the boundary and shape of the velocity gradient area; According to the boundary shape and gradient direction of the velocity gradient zone, the geometric shape of the area is calculated and a grid shape that is consistent with the characteristics of the velocity gradient zone is generated. Combined with the variation law of the elastic modulus and viscosity coefficient in the velocity gradient zone, the density distribution of the updated grid is planned to match the variation law characteristics of the velocity gradient zone. The transition area at the boundary of the velocity gradient zone is smoothed.

[0014] The numerical simulation system for seismic wave propagation in the frozen soil around tunnels includes: The data acquisition unit is configured to: Obtaining geological parameters and environmental parameters of the permafrost layer surrounding the selected tunnel and the tunnel's design magnitude conditions, and setting loading conditions for simulated seismic waves according to the design magnitude conditions, including the frequency, magnitude, and location of the seismic wave source; The modeling unit is configured to: Based on the environmental parameters of the permafrost layer, a layered model is constructed to describe the distribution of physical properties of the permafrost layer at different depths. The layered range is adjusted in combination with the ground temperature field and the freezing depth. The geological distribution of each layer in the layered model is used as an input parameter of the viscoelastic model, which is integrated into a finite element model of the permafrost layer around the tunnel to simulate the dispersion and energy attenuation characteristics of the permafrost layer at different depths. The region of the finite element model is initially divided using an initialized uniform grid, and the region is discretized into a regular grid shape. The update unit is configured as follows: According to the set seismic wave loading conditions, the seismic wave propagation process is simulated in a finite element model, the wave velocity distribution in the frozen soil is calculated, and by comparing the spatial change rate of the wave velocity, the area where the wave velocity gradient changes significantly is identified and marked as a wave velocity gradient area. The wave velocity change amplitude and directionality data of the wave velocity gradient area are extracted to determine its boundary and shape; In combination with the shape and gradient direction of the velocity gradient zone, the boundary curve is calculated based on the geometric form and the shape of the updated grid is generated, and the density distribution of the updated grid is planned in combination with the geological parameters of the gradient change area; The analog unit is configured as: The finite element model is updated and the grid is divided. The seismic wave propagation is re-simulated according to the set seismic wave loading conditions. The finite difference method is used for numerical calculation, and the waveform data and energy changes in each wave velocity gradient area are analyzed.

[0015] 3) Beneficial effects: Compared with the prior art, this invention has the following beneficial effects: The present invention identifies the boundaries and internal characteristics of the velocity gradient zone, and plans the grid shape and density according to the spatial change rate and distribution characteristics of the gradient, ensuring that the calculation area can not only reflect the multi-scale characteristics of the permafrost layer, but also optimize the numerical calculation accuracy of the seismic wave propagation path and energy distribution. It can not only effectively deal with the heterogeneity and temperature coupling of the permafrost layer, but also provide accurate simulation support for the dynamic response of the permafrost layer during tunnel construction. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 A schematic diagram of soil depth-temperature of a high-frozen soil layer in a selected tunnel provided by an embodiment of the present invention; Figure 2 Schematic diagram of the waveform of the numerical simulation method of seismic wave propagation using regular grids in the finite element model; Figure 3 for Figure 2 Schematic diagram of the enlarged waveform in the middle (b) area; Figure 4A schematic flow chart of a numerical simulation method for seismic wave propagation in frozen soil layers around tunnels provided by an embodiment of the present invention; Figure 5 A schematic diagram of a process for evaluating the ultimate stress of a tunnel in a numerical simulation method for seismic wave propagation in frozen soil around a tunnel provided by an embodiment of the present invention; Figure 6 A schematic diagram showing the principle of a numerical simulation method for seismic wave propagation in frozen soil layers around tunnels provided by an embodiment of the present invention; Figure 7 This is a principle block diagram of a numerical simulation system for seismic wave propagation in frozen soil around tunnels provided by an embodiment of the present invention; In the picture: 100, data acquisition unit; 200, modeling unit; 300, updating unit; 400, simulation unit. DETAILED DESCRIPTION

[0017] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0018] In the description of the present invention, it should be understood that the terms "longitudinal", "transverse", "length", "width", "thickness", "up", "down", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inside", "outside" and the like to indicate orientations or positional relationships based on the orientations or positional relationships shown in the accompanying drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be understood as limiting the present invention.

[0019] In addition, the terms "first", "second", etc., if used, are merely used to distinguish and describe, and should not be understood as indicating or implying relative importance.

[0020] It should be noted that, in the absence of conflict, the features in the embodiments of the present invention may be combined with each other.

[0021] During the excavation of plateau tunnels, the inventors discovered that the permafrost surrounding the tunnels is mostly high-frozen soil. Due to its multiphase structure of ice, water, and soil, as well as the directional arrangement of ice crystals, high-frozen soil exhibits significant anisotropy and dispersion during seismic wave propagation. However, traditional wave propagation models typically assume an isotropic or homogeneous medium and fail to fully consider the microstructure and complex mechanical properties of the permafrost. This simplification can easily lead to computational errors, resulting in significant deviations between simulation results and actual observations, which in turn affects the analysis of the stresses on tunnel structures in seismic environments.

[0022] Temperature fluctuations in high-altitude frozen soil directly affect its physical and mechanical properties, especially factors such as ice-water phase transition, pore water content, and viscoelastic changes. Figures 1 to 3 , Figure 1 This is a diagram showing the correlation between humidity and temperature of the permafrost layer around a tunnel in Jingyangling. It can be seen from the figure that there are multiple areas with large temperature changes in the permafrost layer at different depths. Due to temperature changes, the permafrost layer around the tunnel may have local phase change areas, that is, ice-water conversion occurs in some areas, forming high-humidity soil or local melt areas. Temperature changes in these areas have a significant impact on wave speed, energy attenuation, and wave propagation path. You can refer to this. Figure 2 and Figure 3 This phase transition reduces the wave velocity in the phase transition zone, forming a velocity gradient field that affects the wave propagation path and diffraction effects. Furthermore, the increased moisture content in the melt zone enhances the absorption of seismic waves, leading to increased energy loss. Furthermore, changes in acoustic impedance in some localized phase transition zones can cause seismic waves to refract, reflect, or even scatter, affecting the overall wavefield distribution.

[0023] In order to solve the above-mentioned seismic wave simulation problem of high frozen soil layer in plateau tunnels, combined with Figures 4 to 7 As shown, an embodiment of the present invention proposes a numerical simulation method and system for seismic wave propagation in permafrost layers around tunnels. This method comprehensively considers the dynamic characteristics of permafrost, the laws of seismic wave propagation, and the stress mechanism of tunnels, providing an efficient and accurate means of seismic response analysis for tunnel projects in cold plateau areas, thereby improving the seismic safety of the project.

[0024] Now refer to Figure 4 First, S10 is performed: geological and environmental parameters of the permafrost layer surrounding the selected tunnel, as well as the tunnel's design magnitude conditions. Based on these conditions, the loading conditions for the simulated seismic waves are set, including the frequency, magnitude, and location of the seismic wave source. This section primarily prepares for seismic wave simulation, including the collection of simulation input data and the setting of seismic wave loading conditions.

[0025] Specifically, the geological parameters of the permafrost layer are the key inputs for numerical simulation, usually including density, ice content, porosity and viscoelastic parameters. Among them, density affects the wave velocity and dynamic response characteristics of the permafrost layer, ice content determines the stiffness, elastic modulus and dynamic response behavior of the permafrost, porosity affects the mechanical properties during the freeze-thaw process, and viscoelastic parameters determine the deformation characteristics of the permafrost and the propagation characteristics of seismic waves.

[0026] The environmental parameters of the permafrost layer mainly include temperature field and freezing depth. These parameters reflect the thermodynamic environment in the permafrost layer. Combined with the ground temperature field and freezing depth, they determine the physical properties of the permafrost layer and the location of the freeze-thaw interface, and at the same time set dynamic loading conditions for subsequent seismic wave propagation simulations.

[0027] In some embodiments, these geological parameters can be obtained by coring permafrost samples through drilling, and conducting indoor tests such as triaxial shear tests and thermal properties measurements, and deploying temperature sensors or obtaining geothermal field data through satellite remote sensing to understand how the freezing depth changes over time.

[0028] Regarding the selection of seismic wave loading conditions, including the frequency, magnitude, and location of the seismic wave source, in this embodiment, the design magnitude conditions are determined by obtaining the tunnel's geometric parameters, support structure type, and construction method, including the tunnel cross-sectional shape, lining thickness, and support material properties. The seismic performance level of the tunnel is first assessed based on the tunnel's geometric parameters and support type, and the corresponding seismic design specifications are then matched to it. The relationship between the assessed magnitude and source energy is then used to determine the input wave, that is, the energy of the seismic wave used for simulation.

[0029] Taking into account the natural frequencies of tunnel support materials and structures, an appropriate frequency range for seismic wave recording can be selected. Furthermore, because frozen soil responds differently to seismic waves of different frequencies, an appropriate frequency range must be selected based on the frequency dispersion characteristics of the frozen soil. Regarding the epicenter location, or the location of the wave source, in some embodiments, the wave source location is set based on the geometry of the tunnel and the distribution of the frozen soil layer, typically at a certain depth below the tunnel.

[0030] After acquiring the data, S20 is performed: a layered model is constructed based on the environmental parameters of the permafrost layer to describe the distribution of physical properties of the permafrost layer at different depths. The stratification range is adjusted in combination with the ground temperature field and the freezing depth. The geological distribution of each layer in the layered model is used as the input parameter of the viscoelastic model, which is integrated into a finite element model of the permafrost layer around the tunnel to simulate the dispersion and energy attenuation characteristics of the permafrost layer at different depths. The region of the finite element model is initially divided using an initialized uniform grid, and the region is discretized into a regular grid shape. It should be understood here that grid division is the basis of finite element analysis in subsequent data model calculations. Its purpose is to discretize the spatial domain of the finite element model into small units. The initialized grid division is to discretize the spatial domain of the finite element model into small units of uniform and regular shape, so that each unit has a regular shape and equal size. In the subsequent calculation of wave propagation, each grid unit serves as the wave propagation medium. In addition, regularized grids help improve computational efficiency.

[0031] High-altitude frozen soil layers have significant heterogeneity, anisotropy and dispersion characteristics. In order to accurately characterize their seismic response, a layered model is needed to describe the physical property distribution of the frozen soil layer.

[0032] In plateau tunnel environments, the permafrost layer can generally be divided into the following areas: Active layer (seasonal frozen soil): The theoretical depth is 0.5m-13m, and the actual observed depth is generally 0.5m-5m. It is greatly affected by seasonal temperature changes. In summer, the active layer is deeper, has a high moisture content, and freeze-thaw cycles are frequent, showing obvious plastic deformation characteristics; Permafrost (primary permafrost): The theoretical depth is approximately 10m-150m, the temperature is below 0°C, and the ice content is high (30%-60%). It is the main impact area for tunnel construction; Deep permafrost: buried more than 150m deep, less affected by temperature changes, with more obvious ice crystal orientation, showing strong anisotropy, and faster energy attenuation during seismic wave propagation.

[0033] Melting layer / transition layer: A melting layer may form in local areas due to thermal disturbance or groundwater activity. The mechanical properties of the permafrost in this area are significantly reduced, which can easily lead to local damage under earthquake action.

[0034] The energy attenuation characteristics of high-altitude frozen soil are primarily determined by the multiphase interaction between ice, water, and soil. Wave propagation exhibits anisotropy and strong dispersion. Therefore, a viscoelastic model is used to describe the dynamic response of frozen soil. Specifically, in some embodiments, the Kelvin-Voigt model can be used to simulate the dynamic response of frozen soil. Its formula is:

[0035] in, is stress, For strain, is the storage modulus of frozen soil, is the viscosity coefficient.

[0036] After completing the physical property modeling of the permafrost layer, a finite element model of the tunnel and supporting structure needs to be established to simulate the dynamic response under earthquake action. Tunnel linings typically use elastic-plastic constitutive relations, such as the Mohr-Coulomb criterion or the Drucker-Prager model, to accurately describe their load characteristics. Furthermore, when performing initial meshing, a regularized mesh is first performed. In some embodiments, triangulation or octree partitioning is used to generate the initial regular mesh.

[0037] In addition, an initial geostress field needs to be applied to simulate the natural stress distribution of the tunnel, and on this basis, dynamic seismic loads are applied to analyze the vibration response of the tunnel structure under earthquake action.

[0038] To avoid the reflection effects of seismic waves at the boundaries of the computational domain, it is necessary to set absorbing boundaries so that the energy of seismic waves entering the boundaries is effectively absorbed. In some embodiments, common absorbing boundary methods can be used, including viscoelastic boundaries and perfectly matched layers. Among them, viscoelastic boundaries dissipate the energy of incident waves by adding high-damping regions at the boundaries. PML technology, on the other hand, constructs special media to gradually attenuate seismic waves after entering the boundary, thereby avoiding the influence of reflected waves. These methods can effectively reduce boundary reflection errors and make numerical simulations more consistent with actual seismic environments.

[0039] It should be considered that since the tunnel structure is located inside the permafrost layer, the propagation of seismic waves not only affects the underground medium but also causes vibrations on the surface. Therefore, free boundary conditions are required to simulate surface fluctuations.

[0040] Specifically, Neumann boundary conditions can be applied at the surface, that is, assuming that the surface shear stress is zero, allowing the surface to vibrate freely. This can more realistically reflect the propagation effect of seismic waves on the surface, and then analyze the impact of surface vibration on the tunnel structure.

[0041] In plateau environments where tunnels pass through high permafrost, the stress characteristics and failure modes under seismic loads are influenced by complex stress states, freeze-thaw effects, and anisotropic material properties. To accurately assess the seismic response of the tunnel lining and surrounding rock, it is necessary to establish a static geostress field and apply dynamic seismic loads to simulate stress conditions under different earthquake intensities and analyze the ultimate stress state and failure modes.

[0042] This is because the initial stress state of the tunnel is mainly affected by the geological environment, surrounding rock characteristics and geothermal field distribution, mainly including: self-weight stress (vertical stress caused by the overlying rock layer), horizontal tectonic stress (horizontal stress may be greater than vertical stress due to tectonic movement in plateau areas) and freezing stress (additional stress caused by ice crystal expansion and freezing contraction in the permafrost layer).

[0043] In order to accurately describe the static stress distribution of tunnel surrounding rock, Figure 4 and Figure 5 First, S201 is performed; based on the geological parameters of the frozen soil layer and the geometric parameters of the tunnel, a static geostress field is constructed to reflect the stress distribution of the tunnel and the surrounding soil in a natural state. In some embodiments, a three-dimensional geomechanical model can be used, and its governing equation is: in, are the stress tensor components, is the medium density of the surrounding rock, is the gravitational acceleration component.

[0044] Gravity field loading is used to simulate the self-weight stress, and the density of the frozen soil layer is calculated using the following formula: in, is the vertical self-weight stress, is the thickness of the soil cover, that is, the depth, is the density of the permafrost layer that varies with depth.

[0045] Finally, the horizontal stress is calculated using the tectonic stress coefficient , the formula is: Among them, the horizontal stress coefficient Usually determined according to the regional tectonic background, such as the plateau area , which means that the horizontal stress may be greater than the vertical stress.

[0046] Next, S202 is performed: a dynamic seismic load is applied according to the dynamic loading conditions. The intensity and frequency of the dynamic seismic waves are set according to different earthquake intensities to simulate the vibration of the tunnel structure caused by the earthquake. The dynamic load represents the propagation characteristics of the seismic waves and their real-time effects on the structure, and is generally determined by factors such as the frequency and magnitude of the earthquake source, and the propagation path of the seismic waves. During the simulation, the intensity, frequency, and time history of the dynamic seismic waves will vary depending on the preset earthquake intensities, so different parameters need to be input based on the different intensities of the seismic waves. By setting these dynamic loading conditions, the impact of seismic loads on the tunnel structure and the surrounding frozen soil layer can be effectively simulated.

[0047] Furthermore, the tunnel's ultimate stress conditions are taken into consideration, so S203 is performed: radial, tangential, and axial ultimate stress conditions are set to simulate the bearing capacity of the tunnel lining and surrounding rock in different directions and evaluate the structural response to stresses in these directions. This can also be understood as setting radial, tangential, and axial ultimate stress conditions in the model, i.e., the aforementioned stress tensor directions, gravity stress directions, and horizontal stress directions. These stress conditions represent the bearing capacity of the tunnel lining and surrounding rock in different directions, respectively.

[0048] Among them, radial force reflects the pressure of the surrounding soil on the tunnel structure, tangential force mainly involves the friction between the tunnel lining and the soil, and axial force represents the axial tension or compression of the tunnel structure. By setting these extreme force conditions, the force response of the tunnel structure in different directions can be simulated and the areas most prone to damage can be identified.

[0049] Using the dynamic loads and stress conditions set above, the ultimate stress response and corresponding failure modes of the tunnel lining and surrounding rock are evaluated. This is because tunnel structures can experience different failure modes, such as localized collapse and cracking, under different earthquake intensities. By calculating the stress and deformation responses at different locations in the model, the failure modes of the tunnel lining and surrounding rock under ultimate stress conditions can be determined, providing a basis for subsequent safety assessment and design optimization.

[0050] In summary, it can be understood that the calculation of the static geostress field and the application of loads are carried out on the basis of the constructed finite element model, and the study of its dynamic simulation and failure mode is to optimize the parameters of the finite element model before the subsequent numerical simulation so that it can better fit the actual seismic action in the selected tunnel area.

[0051] Later reference Figure 4 Regarding S30: Based on the previously established finite element model, according to the above-mentioned seismic wave loading conditions, the seismic wave propagation process is simulated in the finite element model, and the wave velocity distribution in the frozen soil is calculated. By comparing the spatial change rate of the wave velocity, the area where the wave velocity gradient changes significantly is identified and marked as the wave velocity gradient area. The wave velocity change amplitude and directionality data of the wave velocity gradient area are extracted to determine its boundary and shape.

[0052] The propagation of seismic waves in the permafrost layer is simulated, and the wave velocity value of each grid cell is calculated. The magnitude, source location and frequency of the seismic waves are applied to the simulation area through the finite element model. According to the geological and temperature parameters of the permafrost layer, the process of seismic wave propagation is simulated to calculate the distribution of wave velocity in the entire simulation area, and a wave velocity field associated with the grid cell is obtained. The obtained wave velocity distribution data provides the basis for subsequent gradient area identification.

[0053] Specifically, by comparing the spatial variation rate of wave velocity, the area where the wave velocity gradient changes significantly is identified, and the area where the wave velocity changes significantly is located, that is, the wave velocity gradient area. It can also be understood as calculating the spatial variation rate of wave velocity of each grid unit, that is, the degree of change of wave velocity relative to position. In some embodiments, the gradient formula is often used, that is, the partial derivative of the wave velocity relative to the position. If the gradient value of a certain area exceeds the set threshold, it is marked as an area with significant change in wave velocity gradient, that is, a wave velocity gradient area, which can reflect the impact of sudden changes in physical properties in the permafrost, such as temperature gradient, ice-water phase change, etc. on the propagation of seismic waves, and is a key area of concern.

[0054] Finally, the finite difference method is used to simulate the propagation of seismic waves in the permafrost layer. This method discretizes the calculation area into grid points and uses differential formulas to solve the propagation process of seismic waves, effectively capturing the dispersion effect and energy attenuation characteristics of the permafrost layer. During this process, the physical parameters, freeze-thaw state, and structural complexity of the permafrost layer are precisely set through attributes in the model. During the numerical simulation phase, that is, after calculating the static stress field and simulating the application of loads in the constructed finite element model, S40 is performed: based on the shape and gradient direction of the velocity gradient zone, the boundary curve is calculated based on the geometric form and the shape of the updated grid is generated. The density distribution of the updated grid is planned based on the geological parameters of the gradient change area. After the updated grid, S50 is performed: the finite element model is partitioned into an updated grid, and seismic wave propagation is re-simulated according to the set seismic wave loading conditions. The finite difference method is used for numerical calculations, and the waveform data and energy changes of each velocity gradient zone are analyzed.

[0055] Based on the above, it can be understood that the initial model meshing stage uses an initial uniform regular grid. Although it ensures computational efficiency during the calculation process, it cannot capture enough information, especially for areas with significant fluctuation gradient changes such as the above. It cannot capture enough detailed information. Therefore, it is necessary to dynamically adjust the shape and density of the grid to ensure detailed analysis of high gradient areas. Combined with the above-mentioned extraction of the velocity variation amplitude and directionality data of the velocity gradient area, the purpose is to determine the boundary and shape of the gradient area, thereby generating the preliminary geometric shape of the velocity gradient area, providing a basis for subsequent mesh re-division.

[0056] Specifically, directional data on velocity variations within a velocity gradient region is extracted. This data is typically expressed as a gradient direction vector, where the gradient direction indicates the dominant direction of the velocity variation. The amplitude of the velocity variation is used to determine the boundary location, with regions with smaller amplitudes considered the boundaries of the velocity gradient region. Therefore, it can be understood that the boundaries of the velocity gradient region can be delineated by combining gradient amplitude and directional data.

[0057] The boundary curve is then calculated based on the geometry of the velocity gradient region and the updated mesh shape is generated. This is to optimize the meshing to match the physical properties of the velocity gradient region, thereby improving simulation accuracy and optimizing computational efficiency. In some embodiments, a contour algorithm can be used to first extract the boundary of the gradient region. The extracted boundary curve can be represented by a set of discrete points.

[0058] The boundary curve extracted by the boundary points of the velocity gradient zone needs to extract the inflection points and curvature distribution characteristics of the boundary and use them as a reference for updating the grid division. Among them, for the extracted high curvature area, the velocity changes drastically in this area, so it is necessary to increase the density of the updated grid and make the shape of the updated grid irregular. Here, a four-sided grid or a polygonal grid above can be used to adapt to the curved shape of the curve. In addition, the minimum step size of each grid in the updated grid set in the high curvature area should be smaller than the step size of the regular grid initially set. In the low curvature area, the velocity changes smoothly, and a low-density irregular grid is used. The minimum step size of each grid in this part of the updated grid is also smaller than the step size of the regular grid initially set.

[0059] In summary, it can be understood that the updated grid density distribution should be planned based on the geological parameters of the gradient variation area. Different mesh densities should be set according to the characteristics of the velocity gradient area. Based on the geological parameters of the velocity gradient area, such as the freezing depth and temperature gradient mentioned above, higher mesh densities should be set in areas with larger gradient variations to capture more details. In areas with smaller gradient variations, the mesh density can be appropriately reduced to reduce the computational effort.

[0060] By obtaining the optimized distribution of grid density through logical operations and dynamically adjusting the grid density, we can focus on key areas while maintaining the rationality and efficiency of the overall calculation.

[0061] Furthermore, after updating the grid, some embodiments also consider comparing simulation results with existing permafrost wave propagation patterns, modifying permafrost parameters, and updating the grid's boundary conditions or density. The calculations are then iterated until the error converges to a preset range. This step is intended to gradually converge the deviation between the simulation results and theoretical or experimental data to a preset error range, ultimately ensuring the accuracy and reliability of the simulation.

[0062] Specifically, the simulated seismic wave propagation characteristics are first compared and analyzed with existing permafrost wave propagation patterns, focusing on wave velocity, dispersion characteristics, and energy attenuation. If there are significant deviations, the physical properties of the permafrost layer, such as the density, elastic modulus, and damping coefficient mentioned above, are adjusted to match the dynamic characteristics of actual permafrost.

[0063] Through continuous iterative calculations, the deviation between the simulation results and theoretical or experimental data gradually converges to a preset error range, ultimately ensuring the accuracy and reliability of the simulation and providing a reliable calculation basis for tunnel seismic analysis.

[0064] In summary, by identifying the boundaries and internal characteristics of the velocity gradient zone and planning the grid shape and density according to the spatial change rate and distribution characteristics of the gradient, we ensure that the calculation area can not only reflect the multi-scale characteristics of the permafrost, but also optimize the numerical calculation accuracy of the seismic wave propagation path and energy distribution. This can not only effectively cope with the heterogeneity and temperature coupling of the permafrost, but also provide accurate simulation support for the dynamic response of the permafrost during tunnel construction.

[0065] In some embodiments of the present invention, a simulation system is also proposed to realize the above-mentioned numerical simulation method of seismic wave propagation in the frozen soil layer around the tunnel, which can be referred to herein. Figure 7 , including a data acquisition unit 100, a modeling unit 200, an updating unit 300, and a simulation unit 400.

[0066] Specifically, data acquisition unit 100 collects geological parameters (density, elastic modulus, ice content, etc.) and environmental factors (ground temperature field, freezing depth) of the permafrost layer surrounding the tunnel, as well as tunnel structural parameters (geometry, support type) to determine the tunnel's design magnitude conditions. Based on these design magnitude conditions, the loading conditions for the simulated seismic waves are set, including the frequency, magnitude, and location of the seismic wave source. This data supports modeling of the permafrost's physical properties and provides basic parameters for subsequent calculations.

[0067] The modeling unit 200 constructs the distribution of physical properties of the permafrost layer based on the layered model to ensure that the simulation calculation can reflect the non-uniform characteristics of the permafrost layer. The viscoelastic model is used to characterize the dispersion and energy attenuation characteristics of the permafrost, providing more reasonable material parameters for seismic wave propagation, and constructing a finite element model of the tunnel and support structure. The calculation accuracy is optimized by combining the refined grid division strategy. In the process of finite element modeling, an initialized uniform grid is used to initially divide the area of the finite element model, and the area is discretized into a regular grid shape. The surface fluctuations are simulated by setting free boundaries, and absorbing boundaries are used to prevent seismic wave reflections, thereby improving the stability and reliability of the simulation. In addition, the modeling unit 200 can also generate seismic wave fields under different earthquake intensities according to the regional seismic motion characteristics, and apply dynamic loads. By setting radial, tangential, and axial stress conditions, it is used to analyze the stress state of the tunnel structure and surrounding rock under different working conditions.

[0068] The update calculation unit 300 uses the finite difference method to calculate the propagation behavior of seismic waves in the permafrost layer, simulating its wave velocity, energy attenuation, and scattering characteristics. Based on the calculated wave velocity distribution within the permafrost, the unit compares the spatial rate of change of the wave velocity to identify areas with significant changes in the wave velocity gradient. These areas are labeled as wave velocity gradient zones. The wave velocity variation amplitude and directionality data of the wave velocity gradient zones are extracted to determine their boundaries and shapes. Based on the boundary shape and gradient direction of the wave velocity gradient zone, the unit calculates the geometry of the region and generates a grid shape that matches the characteristics of the wave velocity gradient zone. Combined with the variation patterns of the elastic modulus and viscosity coefficient within the wave velocity gradient zone from the data acquisition unit 100, the unit plans the density distribution of the updated grid to match the variation patterns of the wave velocity gradient zone. Furthermore, the update unit 300 smoothes the transition regions at the boundaries of the wave velocity gradient zones to ensure the continuity of the updated grid division.

[0069] Finally, the simulation unit 400, combined with the updated mesh generated by the update unit, re-meshes the finite element model. It then re-simulates seismic wave propagation based on the configured seismic wave loading conditions, performs numerical calculations using the finite difference method, and analyzes waveform data and energy variations in each velocity gradient region. Through parameter adjustments and multiple rounds of iterative calculations, the simulation error is ensured to converge to a preset range, enhancing the reliability of the calculation results.

[0070] In summary, the simulation system achieves comprehensive simulation and evaluation of seismic wave propagation in the permafrost layer around the tunnel through a process of data acquisition, modeling, calculation, and analysis optimization, thereby improving the scientificity and accuracy of the tunnel's seismic design.

[0071] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. The scope of patent protection of the present invention shall be based on the claims. Any equivalent structural changes made using the description and drawings of the present invention shall be included in the scope of protection of the present invention.

Claims

1. A numerical simulation method for seismic wave propagation in permafrost around tunnels, characterized by: include: Obtaining geological parameters and environmental parameters of the permafrost layer surrounding the selected tunnel and the tunnel's design magnitude conditions, and setting loading conditions for simulated seismic waves according to the design magnitude conditions, including the frequency, magnitude, and location of the seismic wave source; Based on the environmental parameters of the permafrost layer, a layered model is constructed to describe the distribution of physical properties of the permafrost layer at different depths. The layered range is adjusted in combination with the ground temperature field and the freezing depth. The geological distribution of each layer in the layered model is used as an input parameter of the viscoelastic model, which is integrated into a finite element model of the permafrost layer around the tunnel to simulate the dispersion and energy attenuation characteristics of the permafrost layer at different depths. The region of the finite element model is initially divided using an initialized uniform grid, and the region is discretized into a regular grid shape. According to the set seismic wave loading conditions, the seismic wave propagation process is simulated in a finite element model, the wave velocity distribution in the frozen soil is calculated, and by comparing the spatial change rate of the wave velocity, the area where the wave velocity gradient changes significantly is identified and marked as a wave velocity gradient area. The wave velocity change amplitude and directionality data of the wave velocity gradient area are extracted to determine its boundary and shape; In combination with the shape and gradient direction of the velocity gradient zone, the boundary curve is calculated based on the geometric form and the shape of the updated grid is generated, and the density distribution of the updated grid is planned in combination with the geological parameters of the gradient change area; The finite element model is updated and the grid is divided. The seismic wave propagation is re-simulated according to the set seismic wave loading conditions. The finite difference method is used for numerical calculation, and the waveform data and energy changes in each wave velocity gradient area are analyzed.

2. The numerical simulation method for seismic wave propagation in permafrost around tunnels according to claim 1 is characterized in that: The geological parameters of the permafrost layer obtained include density, ice content, porosity and viscoelastic parameters, and the thermodynamic characteristics of the permafrost layer are analyzed by combining drilling, geophysical detection and temperature monitoring methods; The geometric parameters, support structure type and construction method of the tunnel, including the tunnel cross-sectional shape, lining thickness and support material properties, are obtained to determine the design magnitude condition.

3. The numerical simulation method for seismic wave propagation in permafrost around tunnels according to claim 1, characterized in that: The environmental parameters include the ground temperature field and freezing depth of the permafrost layer. By calculating the temperature distribution of permafrost layers at different depths and combining the geological parameters, a stratified model of the permafrost layer is constructed.

4. The numerical simulation method for seismic wave propagation in permafrost around tunnels according to claim 3 is characterized in that: Based on the geological parameters of the permafrost layer, the spatial layers of surface permafrost, high ice content layer and deep permafrost are divided, and the layered model is used to describe the spatial distribution of physical properties of the permafrost layer.

5. The numerical simulation method for seismic wave propagation in permafrost around tunnels according to claim 1, characterized in that: The viscoelastic model is used to simulate the seismic wave propagation characteristics of each of the spatial layers, wherein the seismic wave velocity at different frequencies is analyzed by Fourier transform, and the energy attenuation characteristics of each spatial layer are characterized by combining the quality factor. By adjusting the damping parameters, the simulation results can reflect the wave attenuation characteristics of the permafrost layer.

6. The numerical simulation method for seismic wave propagation in permafrost around tunnels according to claim 5, characterized in that: The layered model provides the physical property distribution of the permafrost layer at different depths, and the viscoelastic model describes the dispersion and energy attenuation characteristics of the permafrost material. The physical property distribution of each layer in the layered model is used as the input parameter of the viscoelastic model and integrated into a finite element model of the permafrost layer around the tunnel, which is used to simulate the propagation characteristics of seismic waves at different depths in the permafrost layer and their energy changes. An absorption boundary is set at the boundary of the computational domain of the finite element model to prevent the reflection and transmission of seismic waves caused by the artificial truncation of the computational domain.

7. The numerical simulation method for seismic wave propagation in permafrost around tunnels according to claim 1, characterized in that: The loading conditions for simulating seismic waves are set according to the design magnitude conditions. Specifically, a static ground stress field of the frozen soil layer is constructed, and dynamic seismic loads are applied to simulate the stress response and corresponding failure modes under different earthquake intensities. Correspondingly: Based on the geological parameters of the permafrost layer and the geometric parameters of the tunnel, a static geostress field is constructed to reflect the stress distribution of the tunnel and the surrounding soil in its natural state. Apply dynamic seismic loads according to the set seismic wave loading conditions, set the intensity and frequency of dynamic seismic waves according to different seismic intensities, and simulate the vibration of the tunnel structure caused by the earthquake; Set radial, tangential, and axial ultimate load conditions to simulate the bearing capacity of the tunnel lining and surrounding rock in different directions and evaluate the stress response of the structure in different directions.

8. The numerical simulation method for seismic wave propagation in permafrost around tunnels according to claim 6, characterized in that: The finite difference method is used to simulate the propagation behavior of seismic waves in frozen soil. The region is discretized into grid points through an initialized uniform grid, and the finite difference formula is used to solve the energy propagation attenuation characteristics of seismic waves between grid points.

9. The numerical simulation method for seismic wave propagation in permafrost around a tunnel according to any one of claims 1 to 8, characterized in that: Specifically, the velocity gradient area marked is defined as the portion where the velocity variation amplitude exceeds a set threshold value, and the velocity variation amplitude and directionality data of the velocity gradient area are extracted to determine the boundary and shape of the velocity gradient area; According to the boundary shape and gradient direction of the velocity gradient zone, the geometric shape of the area is calculated and a grid shape that is consistent with the characteristics of the velocity gradient zone is generated. Combined with the variation law of the elastic modulus and viscosity coefficient in the velocity gradient zone, the density distribution of the updated grid is planned to match the variation law characteristics of the velocity gradient zone. The transition area at the boundary of the velocity gradient zone is smoothed.

10. A numerical simulation system for seismic wave propagation in permafrost around tunnels, characterized by: include: The data acquisition unit is configured to: Obtaining geological parameters and environmental parameters of the permafrost layer surrounding the selected tunnel and the tunnel's design magnitude conditions, and setting loading conditions for simulated seismic waves according to the design magnitude conditions, including the frequency, magnitude, and location of the seismic wave source; The modeling unit is configured to: Based on the environmental parameters of the permafrost layer, a layered model is constructed to describe the distribution of physical properties of the permafrost layer at different depths. The layered range is adjusted in combination with the ground temperature field and the freezing depth. The geological distribution of each layer in the layered model is used as an input parameter of the viscoelastic model, which is integrated into a finite element model of the permafrost layer around the tunnel to simulate the dispersion and energy attenuation characteristics of the permafrost layer at different depths. The region of the finite element model is initially divided using an initialized uniform grid, and the region is discretized into a regular grid shape. The update unit is configured as follows: According to the set seismic wave loading conditions, the seismic wave propagation process is simulated in a finite element model, the wave velocity distribution in the frozen soil is calculated, and by comparing the spatial change rate of the wave velocity, the area where the wave velocity gradient changes significantly is identified and marked as a wave velocity gradient area. The wave velocity change amplitude and directionality data of the wave velocity gradient area are extracted to determine its boundary and shape; Combined with the shape and gradient direction of the velocity gradient zone, the boundary curve is calculated based on the geometric morphology and the shape of the updated grid is generated. The density distribution of the updated grid is planned based on the geological parameters of the gradient change area. The simulation unit is configured to: The finite element model is updated and the grid is divided. The seismic wave propagation is re-simulated according to the set seismic wave loading conditions. The finite difference method is used for numerical calculation, and the waveform data and energy changes in each wave velocity gradient area are analyzed.

Images