CN112329311A - Finite difference simulation method and device for seismic wave propagation and computer storage medium - Google Patents
Finite difference simulation method and device for seismic wave propagation and computer storage medium Download PDFInfo
- Publication number
- CN112329311A CN112329311A CN202011245265.XA CN202011245265A CN112329311A CN 112329311 A CN112329311 A CN 112329311A CN 202011245265 A CN202011245265 A CN 202011245265A CN 112329311 A CN112329311 A CN 112329311A
- Authority
- CN
- China
- Prior art keywords
- equation
- rigid part
- temperature
- finite difference
- stress
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 50
- 238000003860 storage Methods 0.000 title claims abstract description 7
- 238000004088 simulation Methods 0.000 title claims description 39
- 238000010521 absorption reaction Methods 0.000 claims abstract description 30
- 238000004364 calculation method Methods 0.000 claims abstract description 22
- 239000007787 solid Substances 0.000 claims abstract description 13
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 8
- 230000035882 stress Effects 0.000 claims description 19
- 238000012545 processing Methods 0.000 claims description 10
- 230000007480 spreading Effects 0.000 claims description 7
- 238000003892 spreading Methods 0.000 claims description 7
- 239000011159 matrix material Substances 0.000 claims description 6
- 230000000694 effects Effects 0.000 claims description 5
- 238000005452 bending Methods 0.000 claims description 3
- 238000004590 computer program Methods 0.000 claims description 3
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 230000006355 external stress Effects 0.000 claims description 3
- 239000004576 sand Substances 0.000 claims description 3
- 238000012360 testing method Methods 0.000 claims description 3
- 230000007704 transition Effects 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 2
- 238000011160 research Methods 0.000 abstract description 6
- 239000011435 rock Substances 0.000 description 9
- 238000005516 engineering process Methods 0.000 description 5
- 238000011065 in-situ storage Methods 0.000 description 4
- 230000005540 biological transmission Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 239000002131 composite material Substances 0.000 description 2
- 230000005489 elastic deformation Effects 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000010494 dissociation reaction Methods 0.000 description 1
- 230000005593 dissociations Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 239000011148 porous material Substances 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention relates to a method and a device for simulating finite difference of seismic wave propagation and a computer storage medium, wherein the method comprises the following steps: a, converting a thermoelastic equation into a temperature-speed-stress equation; b, decomposing a temperature-speed-stress equation into a rigid part and a non-rigid part based on a time splitting algorithm, and obtaining an analytic solution of the rigid part through numerical calculation; c, carrying out numerical solution on the non-rigid part based on a rotary staggered grid finite difference method, and substituting the analytic solution of the rigid part into the non-rigid part to carry out numerical solution; d applying non-splitting convolution to the thermoelastic equation to completely match the layer absorption boundary; and e, implementing a rotation staggered grid finite difference method and a non-splitting convolution complete matching layer absorption boundary on a temperature-velocity-stress equation, and finally obtaining the waveform record of seismic wave propagation in the high-temperature solid. The method can more accurately describe the propagation condition of the seismic waves in the high-temperature environment, and has important significance for the research of the deep structure.
Description
Technical Field
The invention relates to a finite difference simulation method and device for seismic wave propagation and a computer storage medium, belonging to the field of geophysical seismic exploration research.
Background
The forward modeling of the seismic wave field plays an important role in seismic exploration, is not only reflected in the aspect of theoretical research, but also plays a guiding role in production practice, and the increasing maturity of the technology also improves the accuracy of the application of the inversion technology, so that the geophysical researchers pay more and more attention to the inversion technology. With the gradual deepening of the research of the seismic exploration technology, the seismic exploration is gradually changed from the exploration of a shallow normal-temperature oil and gas reservoir to the exploration of a deep high-temperature oil and gas reservoir. However, the conventional seismic wave theory is usually based on normal temperature medium to perform forward modeling, but due to the change of temperature, the physical property and shallow state of a deep reservoir are different, so that in actual deep-ultra deep seismic exploration, the influence of temperature must be considered to obtain a more accurate underground deep structure, and therefore the seismic wave propagation modeling method in the high temperature medium is of great importance.
The density, elastic modulus, wave velocity and other parameters of the rock are seriously influenced by temperature, especially thermal elastic parameters and deep complex crack structures. With the increase of temperature, the traditional elastic wave equation cannot meet the description of a seismic wave field, and the seismic wave propagation characteristics need to be explained by using a medium model and numerical simulation which are more practical, so that theoretical guidance is provided for the research of deep structures. Accordingly, there is a need in the art for a finite difference simulation method and system for seismic wave propagation in high temperature solid media to solve the above problems.
Disclosure of Invention
Aiming at the outstanding problems, the invention provides a seismic wave propagation finite difference simulation method, a device and a computer storage medium, and the method solves the problem that the traditional seismic wave forward model does not consider the influence of temperature, so that the simulation result is inaccurate.
In order to achieve the purpose, the invention adopts the following technical scheme:
the invention provides a finite difference simulation method for seismic wave propagation in a first aspect, which comprises the following steps:
a, converting a thermoelastic equation into a temperature-speed-stress equation;
b, decomposing the temperature-speed-stress equation into a rigid part and a non-rigid part based on a time splitting algorithm, and obtaining an analytic solution of the rigid part through numerical calculation;
c, carrying out numerical solution on the non-rigid part based on a rotation staggered grid finite difference method, and bringing the analytic solution of the rigid part into the non-rigid part for numerical solution;
d, applying non-splitting convolution complete matching layer absorption boundary to the thermoelastic equation to solve the problem of artificial boundary reflection in seismic wave simulation;
and e, implementing a rotation staggered grid finite difference method and a non-splitting convolution complete matching layer absorption boundary on the temperature-velocity-stress equation, and finally obtaining the waveform record of seismic wave propagation in the high-temperature solid.
Preferably, the thermoelastic equation in step a is as follows:
converting the thermoelastic equation into the temperature-velocity-stress equation:
the above formula can be written in a matrix form
Wherein,
v＝[vx vz σxx σzz σxz T Φ]T (1-4)
wherein,
T,isubscripts of (1), i, uj,jThe subscript of (j, j ═ x, z) denotes the first partial derivative;
uj,jisubscripts of (1), ji, ui,jjSubscripts of (1), jj, T,jjThe subscript of (j, j ═ x, z) denotes the second partial derivative; u. ofi，σij(i, j ═ x, z) and vi(i ═ x, z) are displacement, stress, and velocity components, respectively; t is0Is a reference temperature in an unstressed and unstrained state, T is relative to T0Let Φ represent a first derivative of temperature versus time for computational convenience; ρ represents a bulk density; lambda and mu are Lame coefficients; κ represents the thermal conductivity; c represents the specific heat at constant strain; thermal modulus γ ═ (3 λ +2 μ) αTIn which α isTRepresents a coefficient of thermal expansion; τ represents the relaxation time; f. ofi(i ═ x, z) represents a physical component; f. ofij(i, j ═ x, z) represents applied external stress; q represents a heat source; wherein the parameter vi(i＝x,z)、σij(i,j＝x,z)、fij(i,j＝x,z)、Φ、T、ui、uj,jThe points at the upper end of v represent the derivative with respect to time, and the number of points represents the number of derivatives with respect to time.
Preferably, the step b includes the following specific steps:
decomposing the temperature-velocity-stress equation M into the rigid part M based on a time-splitting algorithmsAnd the non-rigid part MrI.e. M ═ Mr+MsThe rigid part M is obtained by numerical calculationsThe analytic solution of (2); for the rigid part MsWherein：Ms37＝Ms47＝-β，
Respectively represent the rigid portions MsThe values of the third row, the seventh column, the fourth row, the seventh column, and the seventh row, the seventh column of the matrix; the remaining component being 0, i.e. the rigid part MsThe equation is:
obtaining the rigid part M through numerical calculation by taking t as n delta t through discrete time variable, wherein delta t is a time step length, n is an integersThe analytical solution of (a) can be expressed as:
where the superscript "+" indicates the solution obtained from the intermediate of the calculation of (n +1) Δ t from n Δ t.
Preferably, the step c includes the following specific steps:
c1 defining all velocity components in one position and all physical parameters involved in the stress, temperature components and the thermo-elastic equation in another position in a rotating staggered grid;
c2 calculating the spatial derivatives of the rotated interleaved mesh along the diagonal of the spatial mesh, rotating the horizontal x and vertical z directions to the diagonal
And
obtaining a spatial derivative along a diagonal direction;
c3, linearly interpolating the spatial derivative calculation result in the step c2 along the direction of the normal coordinate axis to obtain the spatial derivative of the coordinate axis direction under the rectangular coordinate system, and obtaining the discrete form of each physical quantity in the finite difference of the rotation staggered grids;
c4 bending the rigid part MsIs discretized and brought into the non-rigid part MrAnd carrying out numerical solution.
Preferably, in the finite difference simulation method for seismic wave propagation, in the step c2, the relationship between two coordinate axes x and z is:
wherein,
representing the diagonal lengths of the grid elements, Δ x and Δ z represent the grid spacing, so the first spatial derivative in the coordinate axis direction at rectangular coordinates is:
wherein,
and
as a derivative of x and z,
and
presentation pair
And
a derivative of (a);
in step c3, the discrete format of each velocity, stress and temperature component is:
wherein, cnExpressing a difference coefficient, wherein N is an integer, the above formula shows that the precision is 2N, each subscript represents the coordinate of the position of the physical quantity, and the discrete form of each physical quantity in the rotary staggered grid finite difference is obtained by substituting (1-10) into the equation (1-9);
in the step c4, the rigid part M is put into usesUsing the steps c2 and c3 to discretize the non-rigid part M that is brought into the step brA numerical solution is performed, obtaining from n Δ t the value of (n +1) Δ t:
wherein n is an integer.
Preferably, the step d includes the following specific steps:
d1 expanding rectangular coordinate to complex coordinate by introducing complex factor to obtain corresponding differential operator containing expansion function
d2 transformation of spreading function sx；
d3 introducing an auxiliary variable psixTo avoid calculating the convolution by storing the wavefield information at past times;
d4 obtaining auxiliary variables of the components of stress, velocity and temperature through the steps d 1-d 3
The temperature-velocity response is then determined using equation (1-2)The spatial derivative of the force equation in complex coordinates;
d5 absorption parameter for non-splitting perfect matching layer (d)x，χx，αx) The setting is performed.
Preferably, the step d1 of the finite difference seismic wave propagation simulation method includes the following specific steps:
the rectangular coordinates are extended to complex coordinates:
wherein
dx(s) is an attenuation factor, and omega is an angular frequency, and a corresponding differential operator is obtained
Wherein the spreading function sxComprises the following steps:
in the step d2, in order to improve the absorption effect of the large-angle and low-frequency incidence, the spreading function s is modifiedx：
Wherein, χx＞＞1，αx≥0；
In said step d2, an auxiliary variable ψ is introducedxTo avoid storing past time-of-day wavefield information, substituting equations (1-15) into (1-13) yields:
wherein,
is psixThe Fourier transform of (a) the (b),
fourier transform of ω,. phixThe method has the advantages of no dimension,
comprises the following steps:
by using
Iterative solution of auxiliary variable psix：
Wherein,
the superscript n in (1) is an integer;
therefore, instead of the spatial derivative, the wave equation transition is matched exactly to the layer region wave equation as follows:
in said step d4, the auxiliary variables of the components of stress, velocity and temperature are obtained from equations (1-19)
And then, solving the spatial derivative of the temperature-speed-stress equation under the complex coordinate by using an equation (1-2):
in said step d5, for dxThe following settings are made:
wherein,
wherein L (L is more than or equal to 0 and less than or equal to L) is the distance from a calculation point in a PML (perfect matched layer) to the inner boundary of the PML area; l is the thickness of the PML region; m is the order of the polynomial; r is a theoretical reflection coefficient; vmaxTo be the maximum wave velocity in the medium, let:
wherein, χmaxObtained by testing, amax＝f0，f0The dominant frequency of the seismic source.
Preferably, the seismic wave propagation finite difference simulation method is used for obtaining the seismic wave propagation waveform record in the high-temperature solid by implementing a rotation staggered grid finite difference method and a non-splitting convolution complete matching layer absorption boundary on the temperature-velocity-stress equation.
A second aspect of the present invention provides a seismic wave propagation finite difference simulation apparatus, comprising:
the first processing unit is used for converting a thermoelasticity equation into a temperature-speed-stress equation;
the second processing unit is used for decomposing the temperature-speed-stress equation into a rigid part and a non-rigid part based on a time splitting algorithm, and obtaining an analytic solution of the rigid part through numerical calculation;
the third processing unit is used for carrying out numerical solution on the non-rigid part based on a rotary staggered grid finite difference method and bringing the analytic solution of the rigid part into the non-rigid part for carrying out numerical solution;
a fourth processing unit for applying non-splitting convolution perfect matching layer absorption boundary to the thermoelastic equation to solve the problem of artificial boundary reflection in seismic wave simulation
And the fifth processing unit is used for implementing a rotation staggered grid finite difference method and a non-splitting convolution complete matching layer absorption boundary on the temperature-velocity-stress equation to finally obtain the waveform record of seismic wave propagation in the high-temperature solid.
A third aspect of the invention provides a computer storage medium having stored thereon a computer program which, when being executed by a processor, carries out the steps of the above-mentioned method of finite difference simulation of seismic wave propagation.
Due to the adoption of the technical scheme, the invention has the following advantages:
1. the model rock physical parameters input by the method are accurate, wherein the density rho, Lame constants lambda and mu, the thermal conductivity kappa and the thermal expansion coefficient alphaTThe information of six parameters in total, namely the unit initial volume constant strain ratio heat c, needs the in-situ rock of the target area, so that the physical experiment measurement of the hot rock is carried out, and the obtained parameters are more real and reliable.
2. The method utilizes the thermoelastic equation, considers the influence of thermoelastic parameters compared with the traditional elastic wave equation, can more accurately describe the propagation condition of the seismic wave in a high-temperature environment, and has important significance for the research of deep structures.
3. The finite difference method adopted by the invention is simple, convenient and flexible, and has been widely used in the simulation of elastic media, visco-elastic media, pore media and anisotropic media; the staggered grid technology can effectively improve the calculation precision, is simple and easy to operate, and is widely applied to finite difference method simulation.
4. An important problem in seismic wave simulation is the absorption boundary condition, and numerical simulation is the problem of solving the problem in a limited area, which can generate artificial boundary reflection. The basic idea of the recently developed perfect matching layer absorption boundary (CPML) is to use a special medium as the absorption layer that does not exist in reality to cut off the outer domain problem into a bounded problem. The non-split iterative format complete matching layer absorption boundary which is developed in the microwave propagation calculation and does not need explicit convolution calculation improves the absorption effect at large-angle and low-frequency incidence. The method can accurately investigate the influence of high temperature on seismic wave propagation, provides a theoretical basis for the expansion of the deep and ultra-deep oil and gas exploration field, and provides a reasonable explanation for the abnormal phenomenon of the deep observation data of the earth.
Drawings
FIG. 1 is a flow chart of the main steps of a finite difference simulation method for seismic wave propagation in a high-temperature solid medium according to an embodiment of the present invention;
FIG. 2 is a diagram illustrating a rotation interleaved mesh finite difference operator according to an embodiment of the present invention;
FIG. 3 is a wave field snapshot of an embodiment of the present invention in which the homogeneous medium (in the form of high attenuation of thermal waves) is shown in FIGS. 3(a), 3(b), and 3(c) respectively show the stress component σ under heat source conditionszzVelocity component vzAnd a wavefield snapshot of the temperature field;
FIG. 4 is a wave field snapshot of an embodiment of a homogeneous medium (wave-like propagation of thermal waves) according to the present invention, wherein FIGS. 4(a), 4(b), and 4(c) respectively show stress components σ under heat source conditionszzVelocity component vzAnd a wavefield snapshot of the temperature field;
FIG. 5 shows a non-uniform medium (stress component σ under heat source condition) according to the present inventionzz) The wavefield snapshot of the example, wherein the middle layer of the model in fig. 5(a) is of high thermal conductivity, the upper and lower layers are of low thermal conductivity, the middle layer of the model in fig. 5(b) is of low thermal conductivity, and the upper and lower layers are of high thermal conductivity. In the figure, r represents a reflected wave, c represents a converted wave, and T represents a transmitted wave, for example, PtcTtcP is a P wave converted by transmission of a T wave converted from transmission of a P wave;
FIG. 6 shows a non-uniform medium (stress component σ under heat source condition) according to the present inventionzz) Wavefield snapshotting maps of examples wherein the volcanic body of the model in fig. 6(a) is high thermal conductivity and the outer layer is low thermal conductivity, the volcanic body of the model in fig. 6(b) is low thermal conductivity and the outer layer is high thermal conductivity;
FIG. 7 is a wavefield snapshot for an embodiment of the absorption boundary of the present invention, wherein the T-wave propagates to the CPML inner boundary in FIG. 7(a), and to the middle of the CMPL in FIG. 7(b), wherein the boundary thickness is 40 grid intervals;
FIG. 8 is a composite seismic trace of the present invention, where seismic traces passing through the seismic source are selected for observing the wave absorption at the CPML boundary, and FIG. 8(a), (b), and (c) are composite seismic traces with receiver positions 40, 20, and 1 grid interval, respectively, where YRp is the distance from the receiver position to the outer boundary;
FIG. 9 is a specific flow chart of the finite difference simulation device for seismic wave propagation in high-temperature solid medium according to the present invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the technical solutions of the present invention are described clearly and completely below, and it is obvious that the described embodiments are some, not all embodiments of the present invention. All other embodiments, which can be obtained by a person skilled in the art without any inventive step based on the embodiments of the present invention, are within the scope of the present invention.
The heat conduction equation of the thermoelastic material derived from the basic law of thermodynamics, constitutive relation and Helmholtz free energy expression contains a strain rate in addition to a function of a temperature field to be determined. It shows that the temperature field on the object is not only dependent on the heat source and related thermodynamic parameters, but also affected by the elastic deformation strain rate, or the elastic deformation volume strain rate will change the heat transfer on the object to a certain extent, so the heat conduction equation and the thermoelastic motion equation need to be solved in a coupling way.
Finite difference methods have been widely used in elastic, viscoelastic, porous and anisotropic media simulations because of their simplicity and flexibility. Based on a thermoelastic theory and a finite difference method, the invention develops a high-temperature solid medium seismic wave propagation finite difference simulation method, which comprises the following steps:
as shown in fig. 1, the present invention provides a flow chart of main steps of a seismic wave propagation finite difference simulation method, which comprises the following specific steps:
the thermoelastic equation in step a is:
converting the thermoelastic equation into the temperature-velocity-stress equation:
the above formula can be written in a matrix form
Wherein,
v＝[vx vz σxx σzz σxz T Φ]T (1-4)
wherein,
T,isubscripts of (1), i, uj,jThe subscript of (j, j ═ x, z) denotes the first partial derivative;
uj,jisubscripts of (1), ji, ui,jjSubscripts of (1), jj, T,jjThe subscript of (j, j ═ x, z) denotes the second partial derivative; u. ofi，σijAnd vi(i, j ═ x, z) are displacement, stress, and velocity components, respectively; t is0Is a reference temperature in an unstressed and unstrained state, T is relative to T0Let Φ represent a first derivative of temperature versus time for computational convenience; ρ represents a bulk density; lambda and mu are Lame coefficients; κ represents the thermal conductivity; c represents the specific heat at constant strain; thermal modulus γ ═ (3 λ +2 μ) αTIn which α isTRepresents a coefficient of thermal expansion; τ represents the relaxation time; f. ofi(i ═ x, z) represents a physical component; f. ofij(i, j ═ x, z) represents applied external stress; q represents a heat source; wherein the parameter vi(i＝x,z)、σij(i,j＝x,z)、fij(i,j＝x,z)、Φ、T、ui、uj,jThe points at the upper end of v represent the derivative with respect to time, and the number of points represents the number of time derivatives; for convenience of description, equations (1-2) are described using matrices v, s, Mv.
The step b comprises the following specific steps:
decomposing the temperature-velocity-stress equation M into the rigid part M based on a time-splitting algorithmsAnd the non-rigid part MrI.e. M ═ Mr+MsThe rigid part M is obtained by numerical calculationsThe analytic solution of (2); for the rigid part MsWherein: ms37＝Ms47＝-β，
Ms37、Ms47、Ms77Values of a third row, a seventh column, a fourth row, a seventh column and a seventh row, seventh column of the rigid part matrix are respectively represented; the remaining component being 0, i.e. the rigid part MsThe equation is:
obtaining the rigid part M through numerical calculation by taking t as n delta t through discrete time variable, wherein delta t is a time step length, n is an integersThe analytical solution of (a) can be expressed as:
where the superscript denotes the solution obtained from the intermediate of the calculation of (n +1) Δ t from n Δ t.
The step c comprises the following specific steps:
c1 defining all velocity components in one position and all physical parameters involved in the stress, temperature components and the thermo-elastic equation in another position in a rotating staggered grid;
c2 calculating the spatial derivatives of the rotated interleaved mesh along the diagonal of the spatial mesh, rotating the horizontal x and vertical z directions to the diagonal
And
and obtaining the spatial derivative along the diagonal direction, wherein the relation between x and z of two coordinate axes is as follows:
wherein,
representing the diagonal lengths of the grid elements, Δ x and Δ z represent the grid spacing, so the first spatial derivative in the coordinate axis direction at rectangular coordinates is:
wherein,
and
as a derivative of x and z,
and
presentation pair
And
a derivative of (a);
c3, linearly interpolating each speed, stress and temperature component along the direction of the normal coordinate axis to obtain the spatial derivative of the coordinate axis direction under the rectangular coordinate system, and obtaining the discrete form of each physical quantity in the finite difference of the rotation staggered grid, which is as follows:
wherein, cnAnd (2) expressing difference coefficients, wherein N is an integer, the above formula has 2N-order precision, each subscript represents the coordinate of the position of the physical quantity, and the discrete form of each physical quantity in the rotary staggered grid finite difference is obtained by substituting (1-10) into equations (1-9).
c4 bending the rigid part MsUsing steps c2 and c3 to perform the dissociation of (equations (1-7)))Scattering and bringing into step b said non-rigid part MrA numerical solution is performed, obtaining from n Δ t the value of (n +1) Δ t:
in an embodiment of the present invention, the method further includes: and giving rock physical model parameters, position coordinates of each demodulator probe and each shot point in the observation system, grid unit spacing, time step length, simulation area grid points, simulation precision, information of seismic source wavelets and the like. The rock physical model parameters in the embodiment specifically include: density rho, Lame constants lambda and mu, thermal conductivity kappa and thermal expansion coefficient alphaTAnd the constant strain specific heat c per unit initial volume is six parameters. The information of the six parameters needs to be acquired from the in-situ rock of the target area, and the in-situ rock of the target area is subjected to high-temperature in-situ hot rock physical experiment measurement, so that the obtained parameters are more accurate, real and reliable. The information of the seismic source wavelet comprises information of start-stop frequency, center frequency, stop frequency and the like, if the seismic source selects a heat source, only fast longitudinal waves and thermal waves can be observed in the embodiment but transverse waves cannot be observed, and if the selected force source is the seismic source, transverse waves can be observed. As shown in the following table, the values of the constants and some parameters in the examples of the present invention are as follows:
fig. 3 and 4 are wave field snapshots of the stress field, velocity field, and temperature field of a homogeneous model embodiment of the present invention. With the input of specific parameters, the thermal wave may exhibit different modes of propagation: high attenuation modes and wave-like propagation modes. The embodiment of fig. 5 is a non-uniform model. In the traditional elastic wave numerical simulation, only two types of records exist, namely shear wave records and longitudinal wave records, and a new waveform appears when the propagation of seismic waves in a thermal medium is simulated by using an elastic theory, so that a new idea is provided for seismic exploration.
The step d comprises the following specific steps:
d1 expands the rectangular coordinates into complex coordinates by introducing complex factors, taking x as an example (similar in z direction):
wherein
dx(s) is an attenuation factor, and omega is an angular frequency, and a corresponding differential operator is obtained
Wherein the spreading function sxComprises the following steps:
d2 in order to improve the absorption effect of large-angle and low-frequency incidence, the expansion function s is modifiedx：
Wherein, χx＞＞1，αx≥0；
d3 introducing an auxiliary variable psixTo avoid substituting equations (1-15) into (1-13) by storing the wavefield information at past times, we find:
wherein
Is psixThe Fourier transform of (a) the (b),
fourier transform of ω,. phixThe method has the advantages of no dimension,
can be written as:
by using
Iterative solution of auxiliary variable psix：
Wherein,
the superscript n in (1) is an integer;
thus, instead of spatial derivatives, the wave equation transitions can be matched exactly to the layer PML region wave equation as follows:
d4 obtaining auxiliary variables for each component of stress, velocity and temperature from equations (1-19)
And then, solving the spatial derivative of the temperature-speed-stress equation under the complex coordinate by using an equation (1-2):
d5 absorption parameter for non-splitting perfect matching layer (d)x，χx，αx) Is set up as to dxThe following settings are made:
wherein,
wherein L (L is more than or equal to 0 and less than or equal to L) is the distance from a calculation point in the PML to the inner boundary of the PML area; l is the thickness of the PML region; m is the order of the polynomial; r is a theoretical reflection coefficient; vmaxTo be the maximum wave velocity in the medium, let:
wherein, χmaxObtained by testing, amax＝f0，f0The dominant frequency of the seismic source.
And e, implementing rotary staggered grid finite difference and non-splitting convolution complete matching layer absorption boundary on the temperature-speed-stress equation to obtain the waveform record of seismic wave propagation in the high-temperature solid.
FIG. 7 is a graph of the waveform of seismic wave propagation in a high temperature solid obtained by fully matching the layer absorption boundaries by applying a rotating staggered mesh finite difference and non-splitting convolution according to an embodiment of the homogeneous model of the present invention, wherein the thickness of the absorption boundaries is 40 mesh intervals. The absorption effect of the absorption boundary was visually observed in conjunction with the seismic synthetic record of the example (fig. 8).
Finally, it should be noted that: the above examples are only intended to illustrate the technical solution of the present invention, but not to limit it; although the present invention has been described in detail with reference to the foregoing embodiments, it will be understood by those of ordinary skill in the art that: the technical solutions described in the foregoing embodiments may still be modified, or some technical features may be equivalently replaced; and such modifications or substitutions do not depart from the spirit and scope of the corresponding technical solutions of the embodiments of the present invention.
Claims (10)
1. A seismic wave propagation finite difference simulation method is characterized by comprising the following steps:
a, converting a thermoelastic equation into a temperature-speed-stress equation;
b, decomposing the temperature-speed-stress equation into a rigid part and a non-rigid part based on a time splitting algorithm, and obtaining an analytic solution of the rigid part through numerical calculation;
c, carrying out numerical solution on the non-rigid part based on a rotation staggered grid finite difference method, and bringing the analytic solution of the rigid part into the non-rigid part for numerical solution;
d, applying non-splitting convolution complete matching layer absorption boundary to the thermoelastic equation to solve the problem of artificial boundary reflection in seismic wave simulation;
and e, implementing a rotation staggered grid finite difference method and a non-splitting convolution complete matching layer absorption boundary on the temperature-velocity-stress equation, and finally obtaining the waveform record of seismic wave propagation in the high-temperature solid.
2. The finite difference seismic wave propagation simulation method of claim 1, wherein the thermoelastic equation in step a is:
converting the thermoelastic equation into the temperature-velocity-stress equation:
writing equation (1-2) in matrix form
Wherein,
v＝[vx vz σxx σzz σxz T Φ]T (1-4)
wherein,
T,isubscripts of (1), i, uj,jThe subscript of (j, j ═ x, z) denotes the first partial derivative;
uj,jisubscripts of (1), ji, ui,jjSubscripts of (1), jj, T,jjThe subscript of (j, j ═ x, z) denotes the second partial derivative; u. ofi，σij(i, j ═ x, z) and vi(i ═ x, z) are displacement, stress, and velocity components, respectively; t is0Is a reference temperature in an unstressed and unstrained state, T is relative to T0Let Φ represent a first derivative of temperature versus time for computational convenience; ρ represents a bulk density; lambda and mu are Lame coefficients; κ represents the thermal conductivity; c represents the specific heat at constant strain; thermal modulus γ ═ (3 λ +2 μ) αTIn which α isTRepresents a coefficient of thermal expansion; τ represents the relaxation time; f. ofi(i ═ x, z) represents a physical component; f. ofij(i, j ═ x, z) represents applied external stress; q represents a heat source; wherein the parameter vi(i＝x,z)、σij(i,j＝x,z)、fij(i,j＝x,z)、Φ、T、ui、uj,jThe points at the upper end of v represent the derivative with respect to time, and the number of points represents the number of derivatives with respect to time.
3. The finite difference seismic wave propagation simulation method of claim 1, wherein the step b comprises the following specific steps:
decomposing the temperature-velocity-stress equation M into the rigid part M based on a time-splitting algorithmsAnd the non-rigid part MrI.e. M ═ Mr+MsThe rigid part M is obtained by numerical calculationsThe analytic solution of (2); for the rigid part MsWherein: ms37＝Ms47＝-β，
Respectively represent the rigid portions MsThe values of the third row, the seventh column, the fourth row, the seventh column, and the seventh row, the seventh column of the matrix; the remaining component being 0, i.e. the rigid part MsThe equation is:
obtaining the rigid part M through numerical calculation by taking t as n delta t through discrete time variable, wherein delta t is a time step length, n is an integersIs expressed as:
where the superscript "+" indicates the solution obtained from the intermediate of the calculation of (n +1) Δ t from n Δ t.
4. The finite difference seismic wave propagation simulation method of claim 1, wherein the step c comprises the following specific steps:
c1 defining all velocity components in one position and all physical parameters involved in the stress, temperature components and the thermo-elastic equation in another position in a rotating staggered grid;
c2 calculating the spatial derivatives of the rotated interleaved mesh along the diagonal of the spatial mesh, rotating the horizontal x and vertical z directions to the diagonal
And
obtaining a spatial derivative along a diagonal direction;
c3, linearly interpolating the spatial derivative calculation result in the step c2 along the direction of the normal coordinate axis to obtain the spatial derivative of the coordinate axis direction under the rectangular coordinate system, and obtaining the discrete form of each physical quantity in the finite difference of the rotation staggered grids;
c4 bending the rigid part MsIs discretized and brought into the non-rigid part MrAnd carrying out numerical solution.
5. The finite difference seismic wave propagation simulation method of claim 4, wherein in the step c2, the relationship between two coordinate axes x and z is:
wherein,
representing the diagonal lengths of the grid elements, Δ x and Δ z represent the grid spacing, so the first spatial derivative in the coordinate axis direction at rectangular coordinates is:
wherein,
and
as a derivative of x and z,
and
presentation pair
And
a derivative of (a);
in step c3, the discrete format of each velocity, stress and temperature component is:
wherein, cnExpressing a difference coefficient, wherein N is an integer, the above formula shows that the precision is 2N, each subscript represents the coordinate of the position of the physical quantity, and the discrete form of each physical quantity in the rotary staggered grid finite difference is obtained by substituting (1-10) into the equation (1-9);
in the step c4, the rigid part M is put into usesUsing the steps c2 and c3 to discretize the non-rigid part M that is brought into the step brA numerical solution is performed, obtaining from n Δ t the value of (n +1) Δ t:
wherein n is an integer.
6. The finite difference seismic wave propagation simulation method of claim 2, wherein the step d comprises the following specific steps:
d1 expanding rectangular coordinate to complex coordinate by introducing complex factor to obtain corresponding differential operator containing expansion function
d2 transformation of spreading function sx；
d3 introducing an auxiliary variable psixTo avoid calculating the convolution by storing the wavefield information at past times;
d4 obtaining auxiliary variables of the components of stress, velocity and temperature through the steps d 1-d 3
Then, solving the spatial derivative of the temperature-speed-stress equation under complex coordinates by using an equation (1-2);
d5 absorption parameter for non-splitting perfect matching layer (d)x，χx，αx) The setting is performed.
7. The finite difference seismic wave propagation simulation method of claim 6, wherein the step d1 comprises the following specific steps:
the rectangular coordinates are extended to complex coordinates:
wherein
dx(s) is an attenuation factor, and omega is an angular frequency, and a corresponding differential operator is obtained
Wherein the spreading function sxComprises the following steps:
in the step d2, in order to improve the absorption effect of the large-angle and low-frequency incidence, the spreading function s is modifiedx：
Wherein, χx＞＞1，αx≥0；
In said step d2, an auxiliary variable ψ is introducedxSubstituting equations (1-15) into (1-13) yields:
wherein,
is psixThe Fourier transform of (a) the (b),
fourier transform of ω,. phixThe method has the advantages of no dimension,
comprises the following steps:
by using
Iterative solution of auxiliary variable psix：
Wherein,
the superscript n in (1) is an integer;
therefore, instead of the spatial derivative, the wave equation transition is matched exactly to the layer region wave equation as follows:
in said step d4, the auxiliary variables of the components of stress, velocity and temperature are obtained from equations (1-19)
And then, solving the spatial derivative of the temperature-speed-stress equation under the complex coordinate by using an equation (1-2):
in said step d5, for dxThe following settings are made:
wherein,
wherein L (L is more than or equal to 0 and less than or equal to L) is the distance from a calculation point in the PML to the inner boundary of the PML area; l is the thickness of the PML region; m is the order of the polynomial; r is a theoretical reflection coefficient; vmaxTo be the maximum wave velocity in the medium, let:
wherein, χmaxObtained by testing, amax＝f0，f0The dominant frequency of the seismic source.
8. The seismic wave propagation finite difference simulation method of claim 7, wherein the waveform record of seismic wave propagation in a high temperature solid is obtained by performing a rotation-interleaved-mesh finite difference method and convolution-free perfect-matching-layer absorption boundaries on the temperature-velocity-stress equation.
9. A seismic wave propagation finite difference simulation apparatus, comprising:
the first processing unit is used for converting a thermoelasticity equation into a temperature-speed-stress equation;
the second processing unit is used for decomposing the temperature-speed-stress equation into a rigid part and a non-rigid part based on a time splitting algorithm, and obtaining an analytic solution of the rigid part through numerical calculation;
the third processing unit is used for carrying out numerical solution on the non-rigid part based on a rotary staggered grid finite difference method and bringing the analytic solution of the rigid part into the non-rigid part for carrying out numerical solution;
a fourth processing unit for applying non-splitting convolution perfect matching layer absorption boundary to the thermoelastic equation to solve the problem of artificial boundary reflection in seismic wave simulation
And the fifth processing unit is used for implementing a rotation staggered grid finite difference method and a non-splitting convolution complete matching layer absorption boundary on the temperature-velocity-stress equation to finally obtain the waveform record of seismic wave propagation in the high-temperature solid.
10. A computer storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the steps of the seismic wave propagation finite difference simulation method of any of claims 1-8.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011245265.XA CN112329311A (en) | 2020-11-10 | 2020-11-10 | Finite difference simulation method and device for seismic wave propagation and computer storage medium |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011245265.XA CN112329311A (en) | 2020-11-10 | 2020-11-10 | Finite difference simulation method and device for seismic wave propagation and computer storage medium |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112329311A true CN112329311A (en) | 2021-02-05 |
Family
ID=74317873
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011245265.XA Pending CN112329311A (en) | 2020-11-10 | 2020-11-10 | Finite difference simulation method and device for seismic wave propagation and computer storage medium |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112329311A (en) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113281808A (en) * | 2021-04-22 | 2021-08-20 | 南方海洋科学与工程广东省实验室(湛江) | Anti-dispersion seismic wave forward modeling method, system, device and medium |
| CN114442154A (en) * | 2022-04-11 | 2022-05-06 | 中国石油大学(华东) | Method, system and equipment for simulating propagation of thermal physical seismic waves |
| CN115081274A (en) * | 2022-06-10 | 2022-09-20 | 哈尔滨工业大学 | Method for establishing meshless three-dimensional seismic wave field based on recurrent neural network |
| CN115081267A (en) * | 2022-05-18 | 2022-09-20 | 内蒙古农业大学 | Time-space variable step length finite difference seismic wave numerical simulation method based on optimization |
| CN116451347A (en) * | 2023-04-07 | 2023-07-18 | 长安大学 | Seismic wave numerical simulation method and device for high-speed rail mobile seismic source |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20160034612A1 (en) * | 2014-07-30 | 2016-02-04 | Chevron U.S.A. Inc. | Re-ordered Interpolation and Convolution for Faster Staggered-Grid Processing |
| CN105676280A (en) * | 2016-01-21 | 2016-06-15 | 中国矿业大学(北京) | Two-phase medium geological data obtaining method and device based on rotationally staggered grids |
| CN109490947A (en) * | 2018-10-16 | 2019-03-19 | 中国科学院地质与地球物理研究所 | A kind of high-temperature medium seimic wave propagation analogy method |
-
2020
- 2020-11-10 CN CN202011245265.XA patent/CN112329311A/en active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20160034612A1 (en) * | 2014-07-30 | 2016-02-04 | Chevron U.S.A. Inc. | Re-ordered Interpolation and Convolution for Faster Staggered-Grid Processing |
| CN105676280A (en) * | 2016-01-21 | 2016-06-15 | 中国矿业大学(北京) | Two-phase medium geological data obtaining method and device based on rotationally staggered grids |
| CN109490947A (en) * | 2018-10-16 | 2019-03-19 | 中国科学院地质与地球物理研究所 | A kind of high-temperature medium seimic wave propagation analogy method |
Non-Patent Citations (3)
| Title |
|---|
| JOSE M CARCIONE: "Simulation of wave propagation in linear thermoelastic media", 《GEOPHYSICS》 * |
| 张鲁新 等: "不分裂卷积完全匹配层与旋转交错网格有限差分在孔隙弹性介质模拟中的应用", 《地球物理学报》 * |
| 李元燮 等: "基于旋转交错网格伪谱法和 L-S 理论的热弹性介质波场模拟与特征分析", 《中国地球科学联合学术年会 2020》 * |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113281808A (en) * | 2021-04-22 | 2021-08-20 | 南方海洋科学与工程广东省实验室(湛江) | Anti-dispersion seismic wave forward modeling method, system, device and medium |
| CN113281808B (en) * | 2021-04-22 | 2023-10-20 | 南方海洋科学与工程广东省实验室(湛江) | Anti-dispersion seismic wave forward modeling method, system, device and medium |
| CN114442154A (en) * | 2022-04-11 | 2022-05-06 | 中国石油大学(华东) | Method, system and equipment for simulating propagation of thermal physical seismic waves |
| JP7199132B1 (en) | 2022-04-11 | 2023-01-05 | 中国石油大学（華東） | Modified thermophysical property seismic wave propagation simulation method, system, and equipment |
| JP2023155874A (en) * | 2022-04-11 | 2023-10-23 | 中国石油大学（華東） | Variable thermophysical properties seismic wave propagation simulation method, system, and apparatus |
| CN115081267A (en) * | 2022-05-18 | 2022-09-20 | 内蒙古农业大学 | Time-space variable step length finite difference seismic wave numerical simulation method based on optimization |
| CN115081267B (en) * | 2022-05-18 | 2023-08-29 | 内蒙古农业大学 | Optimization-based space-time variable step length finite difference seismic wave numerical simulation method |
| CN115081274A (en) * | 2022-06-10 | 2022-09-20 | 哈尔滨工业大学 | Method for establishing meshless three-dimensional seismic wave field based on recurrent neural network |
| CN115081274B (en) * | 2022-06-10 | 2023-05-30 | 哈尔滨工业大学 | Method for establishing gridless three-dimensional seismic wave field based on cyclic neural network |
| CN116451347A (en) * | 2023-04-07 | 2023-07-18 | 长安大学 | Seismic wave numerical simulation method and device for high-speed rail mobile seismic source |
| CN116451347B (en) * | 2023-04-07 | 2024-04-23 | 长安大学 | Seismic wave numerical simulation method and device for high-speed rail mobile seismic source |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112329311A (en) | Finite difference simulation method and device for seismic wave propagation and computer storage medium | |
| CN104122585B (en) | Seismic forward simulation method based on elastic wave field resolution of vectors and low-rank decomposition | |
| Ford et al. | A nonhydrostatic finite-element model for three-dimensional stratified oceanic flows. Part I: Model formulation | |
| De La Puente et al. | Discontinuous Galerkin methods for wave propagation in poroelastic media | |
| CN105911584B (en) | Implicit staggered-grid finite difference elastic wave numerical simulation method and device | |
| Ladopoulos | Elastodynamics for Non-linear Seismic Wave Motion in Real-Time Expert Seismology | |
| Liu et al. | Higher-order triangular spectral element method with optimized cubature points for seismic wavefield modeling | |
| CN114139335B (en) | Viscous sound wave simulation method based on single relaxation time lattice Boltzmann model | |
| CN109490947B (en) | High-temperature medium seismic wave propagation simulation method | |
| CN110596754A (en) | Three-dimensional TTI medium qP wave and qSV wave field simulation method | |
| CN104965222A (en) | Three-dimensional longitudinal wave impedance full-waveform inversion method and device | |
| Ladopoulos | Non-linear Seismic Wave Motion in Elastodynamics with Application to Real-Time Expert Seismology | |
| Ladopoulos | Non-linear Elastodynamics by Seismic Wave Motion in Real-Time Expert Seismology | |
| CN113221392B (en) | Construction method of non-uniform viscous acoustic wave propagation model in infinite domain | |
| Shevchenko et al. | Boundary and contact conditions of higher order of accuracy for grid-characteristic schemes in acoustic problems | |
| CN110764146B (en) | Space cross-correlation elastic wave reflection waveform inversion method based on acoustic wave operator | |
| Wang et al. | Modeling the effect of multiscale heterogeneities on wave attenuation and velocity dispersion | |
| CN113221409A (en) | Two-dimensional numerical simulation method and device for acoustic waves with coupled finite elements and boundary elements | |
| CN112630829A (en) | Method and system for analyzing tight sandstone elastic wave attenuation property | |
| CN113589362B (en) | Three-dimensional terrestrial coupled wave forward modeling method | |
| CN113917526B (en) | Forward modeling method based on non-split complete matching layer absorption boundary | |
| CN114910957A (en) | Method for constructing high-temperature medium thermal bomb AVO response model and medium | |
| Jabbari et al. | An exact solution for Lord-Shulman generalized coupled thermoporoelasticity in spherical coordinates | |
| Avendaño et al. | Visco-Acoustic modeling in the frequency domain using a mixedgrid finite-difference method and attenuation-dispersion model | |
| Dastour | Optimal Finite Difference Schemes for the Helmholtz Equation with PML |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20210205 |
|
| RJ01 | Rejection of invention patent application after publication |