DGSE method for multi-scale seismic wave field simulation

Abstract

The invention belongs to the technical field of seismic oscillation simulation, and relates to a DGSE method for multi-scale seismic wave field simulation. According to the method, firstly, non-conformal hexahedral grids are adopted to subdivide different areas, then Riemannian transmission is adopted to realize information transfer between the non-conformal grids to establish a global equation, and finally, second-order frog leaping is introduced to be combined with a local time step length to realize integral rapid solving. The method inherits the unit division flexibility of a discontinuous finite element method, and meanwhile, the high-order spectrum unit is adopted, so that the performance of the method in the aspects of solving precision and numerical dispersion is remarkably improved.

Description

Technical Field

[0001] The invention belongs to the technical field of earthquake motion simulation and relates to a DGSE method for multi-scale earthquake wave field simulation. Background Art

[0002] High-precision wave field simulation is the basis for site effects, forward and inversion of earthquake motions, and input of earthquake motions for major projects. At present, numerical simulation methods based on wave equations mainly include finite element method (FEM), boundary element method (BEM), finite difference method (FDM), pseudo-spectral method (PSM), and spectral element method (SEM). Among them, SEM can be regarded as a combination of FEM and spectral method. It uses GLL integration to obtain the diagonal mass matrix, which can greatly reduce the computational energy consumption during wave field simulation.

[0003] However, a prominent weakness of SEM when dealing with multi-scale complex models is that a large number of additional fine grids will be generated when fitting high-velocity contrast stratum interfaces with structured hexahedral grids. In order to ensure accuracy, a very small time step matching the minimum grid size is required, which greatly reduces the solution efficiency of this method. Although some scholars have tried to develop tetrahedral grid SEM, overly complicated processing will greatly reduce its accuracy. Summary of the invention

[0004] Aiming at the multi-scale problem of seismic wave propagation, the present invention provides a DGSE method for multi-scale seismic wave field simulation. The method introduces the discontinuous Galerkin (DG) method in SEM and uses non-conformal technology on the interface to realize the meshing of the media on both sides with different sizes, which can effectively solve the bottleneck problem of SEM in processing multi-scale complex models. Further combined with fast integration algorithms such as local time step, efficient simulation of broadband ground motion in complex sites can be achieved.

[0005] The technical solution adopted by the present invention to solve the technical problem is: a DGSE method for multi-scale seismic wave field simulation, comprising:

[0006] The computational domain is divided into K non-overlapping subdomains Ω according to the medium parameter characteristics. k , non-conformal hexahedrons are used to divide the grids between different domains; then a reversible mapping relationship is established between each subdomain and the standard reference unit, and N is selected in the reference unit. k The GLL collocation points of the first order are transformed back to the global coordinates through the mapping relationship;

[0007] The three-dimensional wave equation is rewritten as a hyperbolic conservation form of the first-order space-time partial derivative, and both sides of the equation are multiplied by the test function. k Integrating gives the weak form of the equation:

[0008] For different subdomainsΩ kSelect the basis functions of corresponding order to expand the unknown field quantity, and triangulate the polygons generated by the non-conformal interface; use Gaussian integral to calculate the area of ​​the triangles respectively;

[0009] Using the Riemann transport condition as the numerical flux on the non-conformal surface, we obtain Ω k Discrete weak form equations in the domain;

[0010] The global equations are obtained by using the continuity of the numerical flux jumping on the boundaries of different subdomains;

[0011] The second-order frog leaping method combined with the local time step is used to realize the integral solution of the computational domain.

[0012] Preferably, when performing integral calculation on the subdomain boundary cells, the displacement field value in the cell is called for the normal direction entering the boundary, while the displacement field value of the adjacent cell is called for the normal direction leaving the cell.

[0013] Preferably, the integral solution process of the computational domain specifically includes:

[0014] (1) Divide the time domain into N overall time steps Δt = T / N. For any current time t n =nΔt(n=2,3,…,N), establish the second-order leapfrog form of the global equation;

[0015] (2) Using the local time step algorithm to solve the integral, the global equation expanded in the second-order leapfrog form is Z-transformed to obtain the integral time step in each subdomain that satisfies the system stability condition;

[0016] (3) Let the overall time step Δt be an integer multiple of the step length of each subdomain, and calculate the displacement field values ​​of the subdomains at their respective sub-times in sequence. Substituting them into the following formula, we can obtain the displacement field values ​​at the GLL nodes in different subdomains at the next time:

[0017]

[0018] in, For the next moment t n+1 Subdomain Ω k The displacement field value within, Δt k Subdomain Ω k The time step within Subdomain Ω k The stiffness matrix inside, Subdomain Ω k The inverse matrix of the internal stiffness matrix, Subdomain Ω k The internal force vector, is the current time t n Subdomain Ω k The displacement field value inside.

[0019] Preferably, a unit stiffness matrix decomposition strategy is adopted to decompose the stiffness matrix into a diagonal matrix and a square matrix with a completely clear sparse structure.

[0020] Compared with the prior art, the present invention has the following beneficial effects: in response to the multi-scale (from crustal scale to geotechnical scale) challenge of wave field simulation involving seismic wave propagation, a DGSE method for multi-scale seismic wave field simulation is proposed. This technology inherits the unit division flexibility of the discontinuous finite element method, and at the same time adopts high-order spectral units, which significantly improves the performance of this method in terms of solution accuracy and numerical dispersion. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] Figure 1 Schematic diagram of the DGSE method flow for multi-scale seismic wave field simulation provided in an embodiment of the present invention

[0022] Figure 2 The calculation area and its relative position to the seismogenic fault in the embodiment of the present invention;

[0023] Figure 3 Snapshots of three components of the seismic wave field at different times in the embodiment of the present invention. DETAILED DESCRIPTION

[0024] In order to facilitate the understanding of the present invention, the present invention is described in more detail below in conjunction with the accompanying drawings and specific embodiments. However, the present invention can be implemented in many different forms and is not limited to the embodiments described in the specification. On the contrary, the purpose of providing these embodiments is to make the understanding of the disclosure of the present invention more thorough and comprehensive.

[0025] In the embodiments of the present invention, DGSE is the abbreviation of Discontinuous Galerkin Spectral Element, which means discontinuous Galerkin spectral unit in Chinese. GLL is the abbreviation of Gauss-Lobatto-Legendre, which means Gauss-Lobatto-Legendre integral in Chinese.

[0026] Example 1 The present invention provides a DGSE method for multi-scale seismic wave field simulation, such as Figure 1 As shown in the figure, firstly, non-conformal hexahedral meshes are used to divide different regions (i.e., regional decomposition of the physical model), then Riemann transmission is used to realize information transmission between non-conformal meshes to establish global equations (establishment of global equations), and finally, second-order frog leapfrogging is introduced in combination with local time steps to realize fast integral solution (integral solution of the computational domain). The specific method steps are as follows:

[0027] (1) Regional decomposition of the physical model

[0028] First, the computational domain is divided into K non-overlapping subdomains Ω according to the medium parameter characteristics. k, non-conformal hexahedrons are used to divide the mesh between different domains; each subdomain is then meshed with a reference unit (i.e., a unit cube unit: (-1,1) 3 ) establishes a reversible mapping relationship (Formula 1), selects N in the reference unit k The GLL configuration points are transformed back to the global coordinates through the mapping relationship.

[0029]

[0030] in, represents the reference unit, φ k is the distance from the reference unit to the subdomain Ω k Through this mapping relationship, the DGSE method can define the basis functions on the local spectral unit and degenerate the large-scale matrix into a small-scale matrix on the local unit by allowing discontinuities between grid units.

[0031] (2) Establishment of global equations

[0032] First, the three-dimensional wave equation is rewritten into the hyperbolic conservation form of the first-order space-time partial derivative (Equation 2), and both sides of the equation are multiplied by the trial function v. In the subdomain Ω k Integrate to obtain the weak form equation (Equation 3); then for different subdomains Ω k Select the basis function of the corresponding order to expand the unknown field quantity (Formula 4), and triangulate the polygons generated by the non-conformal interface and use Gaussian integral to calculate the triangle area respectively; finally, use the Riemann transmission condition as the numerical flux on the non-conformal surface, and use the continuity of the numerical flux jumping on the boundaries of different subdomains to obtain the global equation (Formula 5). The specific process is as follows:

[0033] When performing integral calculations on the subdomain boundary cells, the displacement field value within the cell is called for the normal direction entering the boundary, while the displacement field value of the adjacent cell is called for the normal direction leaving the cell.

[0034] For the displacement-strain form wave equation, it can be rewritten as the hyperbolic conservation form of the first-order space-time partial derivative:

[0035]

[0036] in, is the partial differential symbol, indicating the partial derivative of a variable; t is time; ρ is mass density; U is the displacement vector; F is the force vector; Represents the gradient operator, which is a vector differential operator; C is the tensor form composed of material elastic constants; σ is the stress tensor.

[0037] Take the trial function space v h ={v∈L 1 (Ω):v| Ωi∈P k (Ω i )}, multiply both sides of equation (2) by the test function v, and in the spectral unit Ω k Integrate on and use Gauss's theorem to reduce the requirement for smoothness of the solution. h Defined as: where the function v belongs to L on the entire area Ω 1 Space, L 1 (Ω) represents the set of all functions that are absolutely integrable over the region Ω, that is, the integral of these functions over the entire region Ω is finite; and v| Ωi ∈P k (Ω i ) indicates that in the subdomain Ω i On the polynomial space P, the function v belongs to k (Ω i ), that is, in each subdomain Ω i The above function v can be expressed as a polynomial of degree not exceeding k. This definition ensures that the space of trial functions v h Contains functions that are absolutely integrable over the entire region Ω and can be expressed as polynomials of a certain degree in each subdomain. In this way, by multiplying both sides of equation (2) by the test function v and integrating over the spectral unit, Gauss's theorem can be used to reduce the smoothness requirements of the solution. Thus, we can obtain:

[0038]

[0039] Among them, Γ k Subdomain Ω k The boundary of the volume, dΩ represents the integral of the volume element.

[0040] Assume that the physical quantity U can be represented by an approximate solution represented by the truncated linear superposition of basis functions, and let the approximate solution space be the same as the trial function space. U| Ωk In subdomain Ω k Expand within

[0041]

[0042] Among them, x is the spatial coordinate, t is the time, U i (t) is the basis function φ corresponding to time t i (x), φ i (x) is the basis function.

[0043] Replace the test function in equation (3) with a highly orthogonal GLL basis function and substitute equation (4) into it, integrating all subdomains Ω k (k=1,2,…,K) can obtain the global equation in matrix form considering multiple scales:

[0044]

[0045] Among them, M is the mass matrix and K is the stiffness matrix. The mass matrix M has the property of block diagonal due to the strict orthogonality of the basis function, so it can be filled in parallel in each process for independent inversion operation.

[0046] (3) Integral solution of computational domain

[0047] First, the time domain is divided into N overall time steps. For any current time t n =nΔt (n=2,3,…,N), establish the second-order leapfrog form of the global equation (Equation 6):

[0048] MU n+1 =(2M-Δt 2 K)U n -MU n-1 ; (6)

[0049] Among them, U n+1 For the next moment t n+1 The displacement vector, U n is the current time t n The displacement vector, U n-1 is the previous moment t n-1 is the displacement vector, and Δt is the time step.

[0050] Then, the local time step algorithm is used to solve the integral, and the global equation expanded in the leapfrog format is Z-transformed to obtain the integral time step in each subdomain that meets the system stability conditions. Z-transform is a discrete time system analysis method that simplifies the analysis of system stability by transforming the integral equation in the time domain to the frequency domain. The equation form obtained after Z-transformation helps to determine the stable integral time step in each subdomain.

[0051] Finally, let the overall time step Δt be an integer multiple of the step length of each subdomain (assuming that two adjacent subdomains Ω i and Ω j The internal step length is Δt1 and Δt2), and the displacement field values ​​of the subdomains at their respective sub-times are calculated in sequence. Substituting them into equation (7) can obtain the displacement field values ​​at the GLL nodes in different subdomains at the next time:

[0052]

[0053] in, For the next moment t n+1 Subdomain Ω k The displacement field value within, Δt k Subdomain Ω k The time step within Subdomain Ω k The stiffness matrix inside, Subdomain Ω k The inverse matrix of the internal stiffness matrix, Subdomain Ω k The internal force vector, is the current time t n Subdomain Ω k The displacement field value inside.

[0054] The unit stiffness matrix decomposition strategy is used to decompose K into a diagonal matrix and a square matrix with a completely clear sparse structure, which can effectively reduce the calculation and storage of velocity-stress. The decomposition process of the stiffness matrix K is as follows:

[0055] 1. Decompose K into a diagonal matrix D and a sparse matrix R, that is, K = D + R.

[0056] 2. Inverse the diagonal matrix D and calculate D-1.

[0057] 3. Use D-1 and R to perform matrix operations to simplify the calculation process.

[0058] Embodiment 2 The method of the present invention is applied to an actual earthquake wave field simulation case. The Sanhe-Pinggu earthquake that occurred in 1679 was the largest earthquake recorded in the vicinity of Beijing. The inferred magnitude was 7.5, and it was only about 50km away from Beijing. According to historical records, this earthquake caused a large number of casualties and property losses due to its large magnitude, wide range, and long duration of aftershocks. However, due to the lack of earthquake records for many historical earthquakes, especially ancient earthquakes, it is difficult to invert the source rupture. Only a method combining actual geological investigation and numerical simulation can be used to roughly reproduce the seismic response of the surrounding area when the ancient earthquake occurred.

[0059] The method of the present invention is used to simulate the wave field of this historical earthquake. The simulation calculation area is shown in Figure 2 As shown, the model longitude range is 116.0°~117.5°, the latitude range is 39.6°~40.5°, and the calculation model obtained according to the target area is about 130km long from east to west, about 120km wide from north to south, and about 40km long. Figure 3 These are three-component snapshots of the seismic wave field at different times when the Sanhe-Pinggu earthquake occurred. Figure 3 The first line is 4 seconds after the rupture begins, and the interval between two adjacent lines is 5 seconds. It can be seen that the method of the present invention well reproduces the process of seismic waves propagating from the rupture starting point to both ends of the fault and then radiating to the entire model surface. The seismic field forms an obvious strong earthquake concentration area and gradually spreads along the fault trace direction.

Claims

1. A DGSE method for multi-scale seismic wave field simulation, characterized in that: include: The computational domain is divided into K non-overlapping subdomains Ω according to the medium parameter characteristics. k , non-conformal hexahedrons are used to divide the grids between different subdomains; a reversible mapping relationship is established between each subdomain and the standard reference unit, and N is selected in the reference unit. k The GLL collocation points of the first order are transformed back to the global coordinates through the mapping relationship; The three-dimensional wave equation is rewritten as a hyperbolic conservation form of the first-order space-time partial derivative, and both sides of the equation are multiplied by the test function. k Integrating gives the weak form of the equation: For different subdomainsΩ k Select the basis functions of corresponding order to expand the unknown field quantity, and triangulate the polygons generated by the non-conformal interface; use Gaussian integral to calculate the area of ​​the triangles respectively; Using the Riemann transport condition as the numerical flux on the non-conformal surface, we obtain Ω k Discrete weak form equations in the domain; The global equations are obtained by using the continuity of the numerical flux jumping on the boundaries of different subdomains; The second-order frog leaping method combined with the local time step is used to realize the integral solution of the computational domain.

2. The DGSE method for multi-scale seismic wavefield simulation according to claim 1, characterized in that: When performing integral calculations on the subdomain boundary cells, the displacement field value within the cell is called for the normal direction entering the boundary, while the displacement field value of the adjacent cell is called for the normal direction leaving the cell.

3. The DGSE method for multi-scale seismic wavefield simulation according to claim 1, characterized in that: The integral solution process of the computational domain specifically includes: (1) Divide the time domain into N overall time steps Δt = T / N. For any current time t n =nΔt(n=2,3,…,N), establish the second-order leapfrog form of the global equation; (2) Using the local time step algorithm to solve the integral, the global equation expanded in the second-order leapfrog form is Z-transformed to obtain the integral time step in each subdomain that satisfies the system stability condition; (3) Let the overall time step Δt be an integer multiple of the step length of each subdomain, and calculate the displacement field values ​​of the subdomains at their respective sub-times in sequence. Substituting them into the following formula, we can obtain the displacement field values ​​at the GLL nodes in different subdomains at the next time: in, For the next moment t n+1 Subdomain Ω k The displacement field value within, Δt k Subdomain Ω k The time step within Subdomain Ω k The stiffness matrix inside, Subdomain Ω k The inverse matrix of the internal stiffness matrix, Subdomain Ω k The internal force vector, is the current time t n Subdomain Ω k The displacement field value inside.

4. The DGSE method for multi-scale seismic wavefield simulation according to claim 3, characterized in that: The element stiffness matrix decomposition strategy is adopted to decompose the stiffness matrix into a diagonal matrix and a square matrix with a completely clear sparse structure.

Images