A DGSE method for multi-scale seismic wavefield simulation

Abstract

The present invention belongs to the field of seismic motion simulation technology and relates to a DGSE method for multi-scale seismic wavefield simulation. This method first uses non-conformal hexahedral meshes to mesh different regions. Riemannian transmission is then used to transfer information between non-conformal meshes to establish global equations. Finally, a second-order frog leapfrog algorithm is introduced in conjunction with a local time step to achieve a fast integral solution. This method inherits the element partitioning flexibility of the discontinuous finite element method (DIFEM) while employing high-order spectral elements, significantly improving its performance in terms of solution accuracy and numerical dispersion.

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 wavefield simulation is fundamental for conducting site effects, forward and inversion seismic motion simulations, and seismic input for major engineering projects. Currently, numerical simulation methods based on wave equations primarily include the finite element method (FEM), boundary element method (BEM), finite difference method (FDM), pseudospectral method (PSM), and spectral element method (SEM). SEM can be considered a combination of FEM and spectral methods. It employs GLL integration to obtain a diagonal mass matrix, significantly reducing computational energy consumption during wavefield simulations.

[0003] However, a significant weakness of SEM when dealing with complex multi-scale models is that fitting high-velocity contrast interfaces with structured hexahedral meshes generates a large number of extraneous fine meshes. To ensure accuracy, a very small time step size must be used to match the minimum mesh size, significantly reducing the method's solution efficiency. While some researchers have attempted to develop tetrahedral mesh SEM, the overly complex processing can significantly reduce its accuracy. Summary of the Invention

[0004] To address the multi-scale problem of seismic wave propagation, this paper proposes a DGSE method for multi-scale seismic wavefield simulation. This method incorporates the discontinuous Galerkin (DG) method into the SEM, employing non-conformal techniques to achieve meshing of the media on either side with different sizes at the interface. This effectively addresses the bottleneck of the SEM in processing complex multi-scale models. Furthermore, combined with fast integration algorithms such as local time steps, this method enables efficient simulation of broadband ground motions in complex sites.

[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 are transformed back to the global coordinates through the mapping relationship;

[0007] Rewrite the three-dimensional wave equation into a hyperbolic conservation form of the first-order space-time partial derivative, and multiply both sides of the equation by a test function. k Integrating this gives the weak form of the equation:

[0008] For different subdomains Ω kSelect the basis functions of the 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 ​​each triangle;

[0009] Using the Riemann transport condition as the numerical flux on the non-conformal surface, we get Ω 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 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.

[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 subjected to Z-transformation to obtain the integral time step in each subdomain that satisfies the system stability conditions;

[0016] (3) Let the overall time step Δt be an integer multiple of the step length of each subdomain. Calculate the displacement field values ​​of the subdomains at their respective sub-times in sequence. Substituting the values ​​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 is the subdomain Ω k The time step within is the subdomain Ω k The stiffness matrix inside, is the subdomain Ω k The inverse matrix of the internal stiffness matrix, is the subdomain Ω k The internal force vector, is the current time t n Subdomain Ω k The displacement field value within.

[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 existing technology, the present invention has the following beneficial effects: in response to the multi-scale (from crustal scale to geotechnical scale) challenges 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 the three components of the seismic wave field at different times in the embodiment of the present invention. DETAILED DESCRIPTION

[0024] To facilitate understanding of the present invention, the present invention will be described in more detail below with reference to 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 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 spectrum 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, first, non-conformal hexahedral meshes are used to divide different regions (i.e., regional decomposition of the physical model), then Riemannian transmission is used to realize information transmission between non-conformal meshes to establish the global equation (establishment of the global equation), and finally, second-order frog leaping is introduced in combination with local time steps to achieve fast integral solution (integral solution of the computational domain). The specific 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 the reference unit (i.e., the 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 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, rewrite the three-dimensional wave equation into the hyperbolic conservation form of the first-order space-time partial derivative (Equation 2), and multiply both sides of the equation by the test function v. In the subdomain Ω k Integrate to obtain the weak form equation (Equation 3); then for different subdomains Ω k The corresponding order basis functions are selected to expand the unknown field (Equation 4), and the polygons generated by the non-conformal interface are triangulated and the areas of the triangles are calculated using Gaussian integrals. Finally, the Riemann transmission condition is used as the numerical flux on the non-conformal surface, and the continuity of the numerical flux jumping on the boundaries of different subdomains is used to obtain the global equation (Equation 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 the elastic constants of the material; σ is the stress tensor.

[0037] Take the trial function space v h ={v∈L 1 (Ω):v| Ωi∈P k (Ω i )}, multiply the left and right sides of equation (2) by the test function v, and in the spectrum 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 ) represents 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 test functions v h Contains functions that are absolutely integrable over the entire region Ω and can be expressed as polynomials of a specific degree in each subdomain. In this way, by multiplying both sides of Equation (2) by the test function v and integrating over the spectral elements, Gauss's theorem can be used to reduce the smoothness requirements of the solution. Thus, we can obtain:

[0038]

[0039] Among them, Γ k is the subdomain Ω k The boundary of , 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, let the approximate solution space be the same as the trial function space, and U| Ωk In the 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 The expansion coefficient on (x), φ i (x) is the basis function.

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

[0044]

[0045] Where M is the mass matrix and K is the stiffness matrix. The mass matrix M has the property of being block diagonal due to the strict orthogonality of the basis functions, and can therefore be populated in parallel in each process for independent inversion operations.

[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), and 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] The local time-step algorithm is then used to solve the integral. The global equations expanded using the leapfrog scheme are then subjected to a Z-transform to determine the integration time steps within each subdomain that satisfy the system stability conditions. The Z-transform is a discrete-time system analysis method that simplifies the analysis of system stability by transforming the integral equations in the time domain to the frequency domain. The resulting equation form helps determine the stable integration time step within 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 formula (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 is the subdomain Ω k The time step within is the subdomain Ω k The stiffness matrix inside, is the subdomain Ω k The inverse matrix of the internal stiffness matrix, is the subdomain Ω k The internal force vector, is the current time t n Subdomain Ω k The displacement field value within.

[0054] The element 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 amount and storage amount 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. Invert 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] Example 2: The method of the present invention is applied to an actual earthquake wavefield simulation case. The Sanhe-Pinggu earthquake that occurred in 1679 is the largest earthquake recorded in the vicinity of Beijing, with an estimated magnitude of 7.5. It is only about 50 km 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 seismic records for many historical earthquakes, especially paleoearthquakes, it is difficult to invert the source rupture. Only a method combining actual geological surveys and numerical simulations can be used to roughly reproduce the ground motion response of the surrounding area at the time of the paleoearthquake.

[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°, and the latitude range is 39.6°~40.5°. According to the target area, the calculated model is about 130km long from east to west, about 120km wide from north to south, and about 40km long. Figure 3 These are snapshots of the three components of the seismic wave field at different times when the Sanhe-Pinggu earthquake occurred. Figure 3 The first row is 4 seconds after the rupture begins, and the interval between each row is 5 seconds. As can be seen, the method of the present invention effectively reproduces the process of seismic waves propagating from the rupture initiation point to both ends of the fault, and then radiating across the entire surface of the model. The seismic field forms a clear concentrated area of ​​strong earthquakes, which gradually spreads along the fault trace.

Claims

1. A DGSE method for multi-scale seismic wavefield simulation, characterized in that: include: The computational domain is divided into non-overlapping K subdomain Ω 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 a N k The GLL collocation points are transformed back to the global coordinates through the mapping relationship; Rewrite the three-dimensional wave equation into a hyperbolic conservation form of the first-order space-time partial derivative, and multiply both sides of the equation by a test function. k Integrating this gives the weak form of the equation: For different subdomains Ω k Select the basis functions of the 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 ​​each triangle; Using the Riemann transport condition as the numerical flux on the non-conformal surface, we get Ω k Discrete weak form equations in the domain; The global equation is obtained by utilizing the continuity of the numerical flux jumping on the boundaries of different subdomains, and the integral solution of the computational domain is realized by using second-order frog leaping combined with local time steps; specifically: (1) Divide the time domain into N overall time step , for any current moment 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 subjected to Z-transformation to obtain the integral time step in each subdomain that satisfies the system stability conditions; (3) Let the overall time step be is an integer multiple of the step size of each subdomain. The displacement field values ​​of the subdomains at their respective sub-times are calculated in sequence. Substituting it into the following formula 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 is the subdomain Ω k The time step within is the subdomain Ω k The stiffness matrix inside, is the subdomain Ω k The inverse matrix of the internal stiffness matrix, is the subdomain Ω k The internal force vector, For the current moment t n Subdomain Ω k The displacement field value within.

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 element stiffness matrix decomposition strategy is used to decompose the stiffness matrix into a diagonal matrix and a square matrix with a completely clear sparse structure.

Images