CN114065585A - Three-dimensional electrical source numerical simulation method based on coulomb specification - Google Patents
Three-dimensional electrical source numerical simulation method based on coulomb specification Download PDFInfo
- Publication number
- CN114065585A CN114065585A CN202111386036.4A CN202111386036A CN114065585A CN 114065585 A CN114065585 A CN 114065585A CN 202111386036 A CN202111386036 A CN 202111386036A CN 114065585 A CN114065585 A CN 114065585A
- Authority
- CN
- China
- Prior art keywords
- space
- field
- domain
- electric field
- background
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 62
- 238000004088 simulation Methods 0.000 title claims abstract description 26
- 230000005684 electric field Effects 0.000 claims abstract description 87
- 238000004364 calculation method Methods 0.000 claims abstract description 46
- 230000008569 process Effects 0.000 claims abstract description 7
- 230000002159 abnormal effect Effects 0.000 claims description 8
- 239000011541 reaction mixture Substances 0.000 claims description 6
- 238000011160 research Methods 0.000 claims description 5
- 238000013461 design Methods 0.000 claims description 4
- 230000035699 permeability Effects 0.000 claims description 3
- 230000005672 electromagnetic field Effects 0.000 abstract description 12
- 238000005516 engineering process Methods 0.000 abstract description 4
- 230000008901 benefit Effects 0.000 abstract description 3
- 238000012937 correction Methods 0.000 abstract description 2
- 239000000243 solution Substances 0.000 description 17
- 238000004422 calculation algorithm Methods 0.000 description 9
- 230000006872 improvement Effects 0.000 description 7
- 230000006870 function Effects 0.000 description 5
- 238000005070 sampling Methods 0.000 description 3
- 230000007547 defect Effects 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000008094 contradictory effect Effects 0.000 description 1
- 230000032798 delamination Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000003345 natural gas Substances 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 239000012088 reference solution Substances 0.000 description 1
- 238000012795 verification 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
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)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a coulombic-specification-based three-dimensional electrical source numerical simulation method, which is characterized in that when background field calculation is carried out, a layer coefficient in bessel integral is calculated in advance, Fourier transform is adopted when a secondary electromagnetic field is solved, a three-dimensional partial differential equation is converted into a one-dimensional ordinary differential equation under wave number, and finally a tightening operator is adopted to correct an electric field in the solving process. The method is based on a secondary field method, divergence correction is not needed, by utilizing the advantage, the Maxwell equation set is converted into a control equation of coulomb scalar bit vector position by utilizing coulomb specifications, a fast method is adopted when a background field is calculated, and a Fourier transform technology is utilized when a secondary field is calculated, so that the calculation efficiency of numerical simulation of the electromagnetic field of the three-dimensional electric source is greatly improved. The method is applied to the field of numerical simulation.
Description
Technical Field
The invention relates to the technical field of numerical simulation, in particular to a three-dimensional electrical source numerical simulation method based on coulomb specifications.
Background
The artificial source electromagnetic method can overcome the defect of weak magnetotelluric source signals, has strong anti-interference capability and is widely applied to the fields of mine resources, natural gas, geothermal and hydrological environmental engineering exploration and the like. The analytic solution of the electromagnetic field of the electrical source in the complex medium is difficult to obtain, so the three-dimensional forward modeling of the frequency domain is mainly numerical simulation, and the currently commonly used three-dimensional forward modeling method is mainly a finite difference method (Streich r.2009.3d fine-difference frequency-domain modeling of controlled-source electromagnetic data: direct solution and optimization for high access. geomatics, 74 (5)):
F95-F105.), finite element method (liu glu, lie country, korea 2017. two-dimensional adaptive vector finite element forward modeling of controlled source electromagnetic field, geophysical press, 60 (12): 4874 and 4886,), finite volume method (penronghua, huxiangyun, hanbo et al 2016. frequency domain controlled source three-dimensional forward computation based on mimicry finite volume method geophysical report, 59 (10): 3927-: 1549-1562). In recent years, the efficiency of three-dimensional forward modeling is still a research hotspot, and the forward modeling technology of the three-dimensional electromagnetic method of the high-efficiency and high-precision electric source still has a great development space.
Disclosure of Invention
Aiming at the defects in the prior art, the invention provides a coulomb-specification-based three-dimensional electrical source numerical simulation method, which solves the problems of large calculation amount, high storage requirement, low calculation efficiency and the like in the analysis process of an electrical source electromagnetic field in a complex medium in the prior art.
In order to achieve the above object, the present invention provides a coulomb specification-based three-dimensional electrical source numerical simulation method, which comprises the following steps:
step 4, performing two-dimensional Fourier transform on the control equation in the horizontal direction to obtain an ordinary differential equation of the space-wave number domain vector bit scalar position, weighting the ordinary differential equation under each wave number by using a Galerkin method by adopting Gaussian offset wave numbers, and solving by using finite elements to obtain the space-wave number domain secondary field vector bit scalar position under each wave number;
step 6, adding the background electric field and the spatial domain secondary electric field to obtain a spatial domain total electric field, and judging whether the spatial domain total electric field meets a given iterative convergence condition:
if the sum of the space-wavenumber domain secondary magnetic field and the background magnetic field is met, obtaining a space-wavenumber domain secondary magnetic field based on a space-wavenumber domain secondary field vector bit scalar bit, obtaining a space domain secondary magnetic field after two-dimensional inverse Fourier transform in the horizontal direction, obtaining a space domain total magnetic field by adding the background magnetic field, and then outputting a space domain total electric field and a space domain total magnetic field;
if not, modifying the total field of the spatial domain by adopting a tightening operator, multiplying the total field by the abnormal conductivity to obtain new scattering current, updating the scattering current of the right-end term in the control equation, and repeating the step 4-6.
As a further improvement of the above scheme, in step 2, bessel integral terms in the background electric field and the background magnetic field are obtained based on a fast filter coefficient method, which specifically comprises the following steps:
the bessel integral is of the form:
setting a background field with N layers of media, wherein the background field is a background electric field or a background magnetic field; c in formula (1)lAnd DlCoefficients in each layered medium for the background field;
mlis an exponential coefficient, klIs a quantity related to the frequency and the conductivity, permittivity and permeability of the respective layer, l ═ 1,2,3 … N; j. the design is a squareiIs an order bessel function;
r is the distance in the horizontal direction, (x, y) is the horizontal coordinate of the calculation point, (x)s,ys) Is a horizontal coordinate of the source point; m is an integral variable; z is the coordinate of the calculation point in the z direction;
when the background field is layered into N layers, the coefficient C of each layer is firstly calculatedl、DlAnd when calculating the bessel integral, multiplying the coefficient corresponding to each node into the integral expression.
As a further improvement of the above scheme, in step 3, the control equation satisfied by the secondary field vector bit scalar bit is:
in the formula
Is the vector bit of the secondary field in the x, y, z directionssIs a quadratic field scalar bit;
kbis the background field wave number, i is an imaginary unit, ω is an angular frequency, ω is 2 π f, f is the calculated frequency, μ0Is a vacuum magnetic conductivity;
as background admittance ratio, σbFor background conductivity,. epsilon.is the dielectric constant;
scattering current in x, y, z direction, σaFor exceptional conductivity, Ex, Ey, Ez are the total electric field in x, y, z directions.
As a further improvement of the above scheme, in step 3, the acquisition process of the space-wavenumber domain secondary field vector scalar bit is as follows:
performing horizontal two-dimensional Fourier transform on the formula (2) as follows:
in the formula (3), the sign on the variable represents a variable in a space-wave number domain, the sign on the variable represents a variable in a space domain, kx and ky represent wave numbers corresponding to x and y directions in Fourier transform;
adopting Gauss offset wave number, weighting ordinary differential equation under each wave number by Galerkin method, and forming finite element method equation set as follows:
in the formula (4), the symbol on the variable represents a variable in the space-wavenumber domain, and the symbol-free represents a variable in the space domain, NeIs a subdivision in the z direction, NiFor Galerkin's weight function, there is a relation Nz＝2×Ne+1，kpIs the wave number, nzIs a normal vector; and (3) loading the boundary condition of the upper boundary with the boundary condition of the upper traveling wave and loading the boundary condition of the lower traveling wave on the lower boundary, and solving an equation set (4) to obtain the vector scalar bit of the space-wavenumber domain secondary field vector.
As a further improvement of the above solution, in step 5, the relation between the space-wave number domain secondary field vector scalar digit and the electric field in the space-wave number domain is:
in the formula (5), the reaction mixture is,
is a secondary electric field in the x direction of the space-wave number domain,
is a space-wave number domain secondary electric field in the y direction,
is a secondary electric field in the z direction of the space-wave number domain.
As a further improvement of the above solution, in step 6, the relation between the space-wavenumber domain secondary field vector scalar bit and the space-wavenumber domain secondary magnetic field is:
in the formula (6), the reaction mixture is,
is a space-wave number domain x-direction secondary magnetic field,
is a space-wave number domain secondary magnetic field in the y direction,
is a secondary magnetic field in the z direction of the space-wave number domain.
As a further improvement of the above scheme, in step 6, the format of the tightening operator to modify the total electric field is:
e on the left side of the equation in the formula (7)(n)(xi,yj,zk) For the corrected total electric field of the space domain, the right side of the equation is the last corrected total electric field E of the space domain(n-1)(xi,yj,zk) And a new space domain total electric field E obtained by the forward calculation(n)(xi,yj,zk)。
As a further improvement of the above scheme, in step 6, the iterative convergence condition of the total electric field in the spatial domain is:
in the formula, if the above formula is satisfied, the iterative convergence condition is satisfied, otherwise, the iterative convergence condition is not satisfied; err is error, E(n)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) The nth solution equation set of the node of (A) is added with the background electric field to obtain a new total electric field of the space domain, E(n-1)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) When n is 1, E is the total electric field of the space domain obtained by the last calculation of the node (a)(0)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) The background electric field of the node.
The three-dimensional electric source numerical simulation method based on the coulomb specification is based on a secondary field method, divergence correction is not needed, the advantage is utilized, Maxwell equation sets are converted into control equations of coulomb scalar bit vector positions by utilizing the coulomb specification, a rapid filter coefficient method is adopted when a background field is calculated, and a Fourier transform technology is utilized when a secondary field is calculated, so that the calculation efficiency of three-dimensional electric source electromagnetic field numerical simulation is greatly improved.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, it is obvious that the drawings in the following description are only some embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to the structures shown in the drawings without creative efforts.
FIG. 1 is a schematic flow chart of a coulombic-specification-based three-dimensional electrical source numerical simulation method according to an embodiment of the invention;
FIG. 2 is a schematic illustration of an exemplary investigation region in an embodiment of the present invention, wherein (a) is a top view and (b) is a side view;
FIG. 3 is a schematic diagram showing the comparison between the electric field and the magnetic field of the present embodiment and the INTEM3D software with the measurement line y being 0m, wherein (a) is the electric field ExComponent (b) is the electric field EyComponent (c) is the electric field EzComponent (d) is a magnetic field HxComponent (e) is the magnetic field HyComponent (f) is the magnetic field HzAnd (4) components.
The implementation, functional features and advantages of the objects of the present invention will be further explained with reference to the accompanying drawings.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
It should be noted that all the directional indicators (such as up, down, left, right, front, and rear … …) in the embodiment of the present invention are only used to explain the relative position relationship between the components, the movement situation, etc. in a specific posture (as shown in the drawing), and if the specific posture is changed, the directional indicator is changed accordingly.
In addition, the descriptions related to "first", "second", etc. in the present invention are only for descriptive purposes and are not to be construed as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Thus, a feature defined as "first" or "second" may explicitly or implicitly include at least one such feature. In the description of the present invention, "a plurality" means at least two, e.g., two, three, etc., unless specifically limited otherwise.
In the present invention, unless otherwise expressly stated or limited, the terms "connected," "secured," and the like are to be construed broadly, and for example, "secured" may be a fixed connection, a removable connection, or an integral part; the connection can be mechanical connection, electrical connection, physical connection or wireless communication connection; they may be directly connected or indirectly connected through intervening media, or they may be connected internally or in any other suitable relationship, unless expressly stated otherwise. The specific meanings of the above terms in the present invention can be understood by those skilled in the art according to specific situations.
In addition, the technical solutions in the embodiments of the present invention may be combined with each other, but it must be based on the realization of those skilled in the art, and when the technical solutions are contradictory or cannot be realized, such a combination of technical solutions should not be considered to exist, and is not within the protection scope of the present invention.
The embodiment discloses a coulombic-specification-based three-dimensional electrical source numerical simulation method, which is used for calculating a layer coefficient in bessel integral in advance when background field calculation is carried out, so that the solving of the coefficient is reduced, and the calculation efficiency is improved. And Fourier transform is adopted when the secondary electromagnetic field is solved, the three-dimensional partial differential equation is converted into a one-dimensional ordinary differential equation under the wave number, the storage dimension of the unknown quantity is reduced from three dimensions to one dimension, the calculation quantity and the storage requirement are greatly reduced, and the calculation efficiency is higher. And finally, correcting the electric field in the solving process by adopting a tightening operator, wherein the algorithm is stable and convergent and is suitable for any model with the conductivity not being 0. After the background field calculation and the secondary field calculation are processed, the three-dimensional electrical source electromagnetic field is rapidly and accurately numerically simulated, a new technology is provided for the forward modeling of the large-scale electrical source high-efficiency high-precision numerical simulation, and the method has an important invention value for three-dimensional high-efficiency fine inversion. Referring to fig. 1, the method includes the following steps 1-6.
The analytic solution of the background electromagnetic field has a bessel integral term, and the bessel integral term in the background electric field and the background magnetic field is obtained based on a fast filter coefficient method in the embodiment, and the bessel integral is in the form of:
setting a background field with N layers of media, wherein the background field is a background electric field or a background magnetic field; c in formula (1)lAnd DlCoefficients in each layered medium for the background field;
mlis an exponential coefficient, klIs a quantity related to the frequency and the conductivity, permittivity and permeability of the respective layer, l ═ 1,2,3 … N; j. the design is a squareiIs an order bessel function;
r is the distance in the horizontal direction, (x, y) is the horizontal coordinate of the calculation point, (x)s,ys) Is a horizontal coordinate of the source point; m is an integral variable; z is the calculated point vertical z coordinate. When the bessel integral of the formula is calculated using a filtering algorithm, m is related to r, which is found in detail in the literature (Chua John 2014 fast Hankel transform and its application in forward calculations 29(03): 1384-.
Because z and r of each node are different, when the traditional filter coefficient method is used for calculating the integral, each node calculates C oncelAnd Dl. When the number of nodes in the x, y and z directions is Nx、Ny、NzThe number of times of calculation of the layer coefficient is Nx×Ny×Nz. This embodiment takes into account the layering property of the background field, and when the background field is layered into N layers, the coefficient C of each layer is first calculatedl、DlWhen calculating the bessel integral, the coefficient corresponding to each node is multiplied into the integral expression, at this time, the layer coefficient only needs to calculate Nx×NyTimes.2 times. Three-dimensional electromagnetic numerical simulation, delamination of rarely more than 10 layers is generally considered. And when the three-dimensional forward is played, the number of nodes calculated in the vertical direction is generally not less than 10. At the moment, the rapid filter coefficient algorithm can reduce the calculation times of the coefficient, thereby improving the calculation efficiency.
Finally, calculating a background field to obtain a background electric field of each subdivision node
Background magnetic field
Wherein (x)i,yj,zk) Denotes a subdivision node coordinate of (i, j, k), i being 1,2, … Nx，j＝1,2,…,Ny，k＝1,2,…,Nz。
And 3, converting the Maxwell equation set into an equation about vector bits and scalar bits by adopting a coulomb specification, obtaining a control equation satisfied by the scalar bits of the vector bits of the secondary field by using a secondary field method, and taking the product of the background electric field and the background magnetic field as the scattering current of a right-end item in the control equation. Wherein the control equation satisfied by the secondary field vector bit scalar bit is:
in the formula
Is the vector bit of the secondary field in the x, y, z directionssIs a quadratic field scalar bit;
kbis the background field wave number, i is an imaginary unit, ω is an angular frequency, ω is 2 π f, f is the calculated frequency, μ0Is a vacuum magnetic conductivity;
as background admittance ratio, σbFor background conductivity,. epsilon.is the dielectric constant;
scattering current in x, y, z direction, σaFor exceptional conductivity, Ex, Ey, Ez are the total electric field in x, y, z directions.
And 4, performing two-dimensional Fourier transform on the control equation in the horizontal direction to obtain an ordinary differential equation of the space-wave number domain vector bit scalar position, weighting the ordinary differential equation under each wave number by using a Galerkin method by adopting Gaussian offset wave numbers, and solving by using finite elements to obtain the space-wave number domain secondary field vector bit scalar position under each wave number. The acquisition process of the space-wavenumber domain secondary field vector bit scalar bit comprises the following steps:
performing horizontal two-dimensional Fourier transform on the formula (2) as follows:
in the formula (3), the symbol on the variable represents a variable in the space-wave number domain, the symbol on the variable represents a variable in the space domain, and kx and ky represent wave numbers corresponding to x and y directions in fourier transform.
This example uses a gaussian offset wavenumber with 4 gaussian points. The wave number is selected according to the number of sampling points and the sampling interval, and the standard wave number is selected according to the following rule: when the number of samples is M and the sampling interval is Deltax, the fundamental wave Deltak is 2 pi/(M x Deltax), the wave number kpThe number of the selected M is also M.
When M is an even number, the selection rule is as follows:
when M is an odd number, the wave number kpThe selection is as follows:
the gaussian offset wavenumber is chosen as follows:
kp＝(p+tg)·Δk,g＝1,2,3,4
4 Gaussian points tgAnd its corresponding coefficient AgRespectively as follows:
t1＝-0.8611363115940526，A1＝0.3478548451374538；
t2＝-0.3399810435848563，A2＝0.6521451548625461；
t3＝0.3399810435848563，A3＝0.6521451548625461；
t4＝0.8611363115940526，A4＝0.3478548451374538。
the number of wave numbers is 4M.
Weighting ordinary differential equations under each wave number by using a Galerkin method, wherein a finite element method equation system is formed by:
in the formula (4), the symbol on the variable represents a variable in the space-wavenumber domain, and the symbol-free represents a variable in the space domain, NeIs a subdivision in the z direction, NiFor the Galerkin weight function, the present embodiment uses a quadratic interpolation shape function. Existence of a relation Nz＝2×Ne+1，kpIs the wave number, nzIs a normal vector; and (3) loading the boundary condition of the upper boundary with the boundary condition of the upper traveling wave and loading the boundary condition of the lower traveling wave on the lower boundary, and solving an equation set (4) to obtain the vector scalar bit of the space-wavenumber domain secondary field vector.
And 5, obtaining a space-wavenumber domain secondary electric field based on a relational expression between a space-wavenumber domain secondary field vector bit scalar digit and the electric field in the space-wavenumber domain, and obtaining the space domain secondary electric field by utilizing horizontal two-dimensional inverse Fourier transform. The relation between the space-wavenumber domain secondary field vector bit scalar position and the electric field in the space-wavenumber domain is as follows:
in the formula (5), the reaction mixture is,
is a secondary electric field in the x direction of the space-wave number domain,
is a space-wave number domain secondary electric field in the y direction,
the space-wave number domain z-direction secondary electric field can be obtained by using the formula (5)And obtaining a space-wave number domain secondary electric field.
And 6, adding the background electric field and the spatial domain secondary electric field to obtain a spatial domain total electric field, and judging whether the spatial domain total electric field meets a given iterative convergence condition. The iterative convergence condition of the total electric field of the spatial domain is as follows:
in the formula, if the above formula is satisfied, the iterative convergence condition is satisfied, otherwise, the iterative convergence condition is not satisfied; err is error, E(n)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) The nth solution equation set of the node of (A) is added with the background electric field to obtain a new total electric field of the space domain, E(n-1)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) When n is 1, E is the total electric field of the space domain obtained by the last calculation of the node (a)(0)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) The background electric field of the node.
If the space-wavenumber domain secondary magnetic field is satisfied, obtaining a space-wavenumber domain secondary magnetic field based on a space-wavenumber domain secondary field vector scalar digit, obtaining a space domain secondary magnetic field after two-dimensional inverse Fourier transform in the horizontal direction, adding a background magnetic field to obtain a space domain total magnetic field, and outputting a space domain total electric field and a space domain total magnetic field, wherein the relation between the space-wavenumber domain secondary field vector scalar digit and the space-wavenumber domain secondary magnetic field is as follows:
in the formula (6), the reaction mixture is,
is a space-wave number domain x-direction secondary magnetic field,
is a space-wave number domain secondary magnetic field in the y direction,
is a secondary magnetic field in the z direction of the space-wave number domain.
If not, modifying the total field of the spatial domain by adopting a tightening operator, multiplying the total field by the abnormal conductivity to obtain new scattering current, updating the scattering current of a right end item in the control equation, and repeating the step 4-6, wherein the format of the total field modified by the tightening operator is as follows:
e on the left side of the equation in the formula (7)(n)(xi,yj,zk) For the corrected total electric field of the space domain, the right side of the equation is the last corrected total electric field E of the space domain(n-1)(xi,yj,zk) And a new space domain total electric field E obtained by the forward calculation(n)(xi,yj,zk) Corrected total electric field introduction
And solving a new scattering current, and then carrying out a new round of solution.
In the following, the correctness and efficiency of the coulomb specification-based three-dimensional electrical source numerical simulation provided in this embodiment are verified by combining with a design model example, and the tested computer is an intel (r) core (tm) i7-6700HQ CPU with a main frequency of 2.60GHz and a memory of 16GB and 64-bit win10 system.
(1) Correctness verification
The XOY plane projection of the model is shown in FIG. 2, the background is a uniform half-space medium, the upper half-space is air, and the air conductivity σ is0＝10-12S/m, background conductivity of lower half space σbThe method comprises the following steps of (0.01S/m), the frequency is 1Hz, the calculation range is-1000 m in the x direction, the calculation range is-1000 m in the y direction, the calculation range in the z direction is 0-1000 m, the number of nodes of a subdivision grid is 101 multiplied by 101, the three directions are uniformly subdivided, and subdivision intervals in the x direction are dividedThe subdivision intervals delta y in the delta x and y directions are both 20m, and the subdivision interval delta z in the z direction is 10 m. The range of the abnormal body is-100 to 100m in the x direction, -200 to 200m in the y direction, and 400 to 800m in the z direction, and the conductivity sigma of the abnormal body is 0.1S/m. The electrical source was 200m long along the x-direction, the source had a center coordinate of (0, -2000,0.01), and the current was 10A. The correctness of the method is verified by taking the calculation result of three-dimensional forward modeling software INTEM3D based on the integral equation method developed by the university of Utah in the United states as a reference.
Fig. 3 is a comparison between the method of the present embodiment and the software INTEM3D for calculating the electric field and the magnetic field with the measurement line y being 0m, and the numerical solution of the present invention is the numerical solution of the method of the present embodiment; the INTEM3D reference solution is a numerical solution for the software INTEM3D, and as can be seen in FIG. 3, the two numerical solutions match well, with a zero value point and EzThe deviation is large because the numerical value is small, but the precision meets the requirement. In conclusion, the calculation error of the electromagnetic field is small, the correctness of the method of the embodiment is verified, and the method of the embodiment has high calculation precision.
(2) Primary and secondary field computational efficiency
When the electromagnetic fields of 101 × 101 sections are calculated by respectively using a traditional filter coefficient algorithm and a rapid filter coefficient algorithm, the calculation time of the traditional algorithm is 70.3s, and the occupied memory is 102.5 MB; the calculation time of the fast algorithm is 39.8s, and the occupied memory is 120.1 MB; under the condition of almost no memory, the calculation time is shortened by 2/5, so that the calculation speed of the rapid algorithm is high, and the calculation efficiency is high.
For the second field calculation, the electric source electromagnetic field value simulation calculation of the embodiment achieves the expected convergence with the relative residual error of 10-4Then, 12 iterations are needed, and forward time of 101 × 101 nodes is 100.35s, and memory 2034MB is occupied. And the calculation time of the background field is added to be 39.8s, and the forward modeling consumes 140s in total, so that the algorithm has high calculation speed, small occupied memory and high calculation efficiency.
The above description is only a preferred embodiment of the present invention, and is not intended to limit the scope of the present invention, and all modifications and equivalents of the present invention, which are made by the contents of the present specification and the accompanying drawings, or directly/indirectly applied to other related technical fields, are included in the scope of the present invention.
Claims (8)
1. A three-dimensional electrical source numerical simulation method based on coulomb specifications is characterized by comprising the following steps:
step 1, carrying out hexahedral subdivision on a research area along x, y and z directions, wherein uniform subdivision is carried out in the x and y directions, uniform subdivision is carried out in the z direction, or non-uniform subdivision is carried out in the z direction, and the number of subdivision nodes in the x, y and z directions is N respectivelyx、Ny、Nz；
Step 2, acquiring parameters of the earth electric model, including background layer parameters, calculation frequency and abnormal conductivity parameters, and acquiring a background electric field and a background magnetic field based on the calculation frequency, the source position and the given background layer parameters;
step 3, converting the Maxwell equation set into an equation about vector bits and scalar bits by adopting a coulomb specification, obtaining a control equation satisfied by the scalar bits of the vector bits of the secondary field by using a secondary field method, and taking the product of the background electric field and the background magnetic field as the scattering current of a right-end item in the control equation;
step 4, performing two-dimensional Fourier transform on the control equation in the horizontal direction to obtain an ordinary differential equation of the space-wave number domain vector bit scalar position, weighting the ordinary differential equation under each wave number by using a Galerkin method by adopting Gaussian offset wave numbers, and solving by using finite elements to obtain the space-wave number domain secondary field vector bit scalar position under each wave number;
step 5, obtaining a space-wavenumber domain secondary electric field based on a relational expression between a space-wavenumber domain secondary field vector bit scalar digit and an electric field in a space-wavenumber domain, and obtaining the space domain secondary electric field by utilizing horizontal two-dimensional inverse Fourier transform;
step 6, adding the background electric field and the spatial domain secondary electric field to obtain a spatial domain total electric field, and judging whether the spatial domain total electric field meets a given iterative convergence condition:
if the sum of the space-wavenumber domain secondary magnetic field and the background magnetic field is met, obtaining a space-wavenumber domain secondary magnetic field based on a space-wavenumber domain secondary field vector bit scalar bit, obtaining a space domain secondary magnetic field after two-dimensional inverse Fourier transform in the horizontal direction, obtaining a space domain total magnetic field by adding the background magnetic field, and then outputting a space domain total electric field and a space domain total magnetic field;
if not, modifying the total field of the spatial domain by adopting a tightening operator, multiplying the total field by the abnormal conductivity to obtain new scattering current, updating the scattering current of the right-end term in the control equation, and repeating the step 4-6.
2. The coulombic-specification-based three-dimensional electrical source numerical simulation method of claim 1, wherein in the step 2, the bessel integral terms in the background electric field and the background magnetic field are obtained based on a fast filter coefficient method, and the specific process is as follows:
the bessel integral is of the form:
setting a background field with N layers of media, wherein the background field is a background electric field or a background magnetic field; c in formula (1)lAnd DlCoefficients in each layered medium for the background field;
mlis an exponential coefficient, klIs a quantity related to the frequency and the conductivity, permittivity and permeability of the respective layer, l ═ 1,2,3 … N; j. the design is a squareiIs an order bessel function;
r is the distance in the horizontal direction, (x, y) is the coordinate of the calculation point in the horizontal direction, (x)s,ys) Is the horizontal direction coordinate of the source point; m is an integral variable; z is the coordinate of the calculation point in the vertical direction;
when the background field is layered into N layers, the coefficient C of each layer is firstly calculatedl、DlAnd when calculating the bessel integral, multiplying the coefficient corresponding to each node into the integral expression.
3. The coulombic-specification-based three-dimensional electrical source numerical simulation method of claim 1, wherein in the step 3, the quadratic field vector scalar bit satisfies the control equation:
in the formula
Is the vector bit of the secondary field in the x, y, z directionssIs a quadratic field scalar bit;
kbis the background field wave number, i is an imaginary unit, ω is an angular frequency, ω is 2 π f, f is the calculated frequency, μ0Is a vacuum magnetic conductivity;
as background admittance ratio, σbFor background conductivity,. epsilon.is the dielectric constant;
scattering current in x, y, z direction, σaFor exceptional conductivity, Ex, Ey, Ez are the total electric field in x, y, z directions.
4. The coulombic-specification-based three-dimensional electrical source numerical simulation method of claim 3, wherein in the step 3, the acquisition process of the space-wave number domain secondary field vector scalar digits is as follows:
performing horizontal two-dimensional Fourier transform on the formula (2) as follows:
in the formula (3), the sign on the variable represents a variable in a space-wave number domain, the sign on the variable represents a variable in a space domain, kx and ky represent wave numbers corresponding to x and y directions in Fourier transform;
adopting Gauss offset wave number, weighting ordinary differential equation under each wave number by Galerkin method, and forming finite element method equation set as follows:
in the formula (4), the symbol on the variable represents a variable in the space-wavenumber domain, and the symbol-free represents a variable in the space domain, NeIs a subdivision in the z direction, NiFor Galerkin's weight function, there is a relation Nz＝2×Ne+1，kpIs the wave number, nzIs a normal vector; and (3) loading the boundary condition of the upper boundary with the boundary condition of the upper traveling wave and loading the boundary condition of the lower traveling wave on the lower boundary, and solving an equation set (4) to obtain the vector scalar bit of the space-wavenumber domain secondary field vector.
5. The coulombic-specification-based three-dimensional electrical source numerical simulation method of claim 4, wherein in the step 5, the relationship between the electric field and the secondary field vector scalar potential in the space-wave number domain is as follows:
in the formula (5), the reaction mixture is,
is a secondary electric field in the x direction of the space-wave number domain,
is a space-wave number domain secondary electric field in the y direction,
is a secondary electric field in the z direction of the space-wave number domain.
6. The coulombic-specification-based three-dimensional electrical source numerical simulation method of claim 4, wherein in step 6, the relationship between the space-wavenumber-domain secondary field vector scalar bit and the space-wavenumber-domain secondary magnetic field is:
in the formula (6), the reaction mixture is,
is a space-wave number domain x-direction secondary magnetic field,
is a space-wave number domain secondary magnetic field in the y direction,
is a secondary magnetic field in the z direction of the space-wave number domain.
7. The coulombic-specification-based three-dimensional electrical source numerical simulation method of claim 4, wherein in the step 6, the format of the total electric field corrected by the tightening operator is as follows:
e on the left side of the equation in the formula (7)(n)(xi,yj,zk) For the corrected total electric field of the space domain, the right side of the equation is the last corrected total electric field E of the space domain(n-1)(xi,yj,zk) And a new space domain total electric field E obtained by the forward calculation(n)(xi,yj,zk)。
8. The coulombic-specification-based three-dimensional electrical source numerical simulation method according to any one of claims 1 to 7, wherein in the step 6, the iterative convergence condition of the total electric field in the spatial domain is as follows:
in the formula, if the above formula is satisfied, the iterative convergence condition is satisfied, otherwise, the iterative convergence condition is not satisfied; err is error, E(n)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) The nth solution equation set of the node of (A) is added with the background electric field to obtain a new total electric field of the space domain, E(n-1)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) When n is 1, E is the total electric field of the space domain obtained by the last calculation of the node (a)(0)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) The background electric field of the node.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111386036.4A CN114065585B (en) | 2021-11-22 | 2021-11-22 | Three-dimensional electrical source numerical simulation method based on coulomb specification |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111386036.4A CN114065585B (en) | 2021-11-22 | 2021-11-22 | Three-dimensional electrical source numerical simulation method based on coulomb specification |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114065585A true CN114065585A (en) | 2022-02-18 |
| CN114065585B CN114065585B (en) | 2024-05-10 |
Family
ID=80278884
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111386036.4A Active CN114065585B (en) | 2021-11-22 | 2021-11-22 | Three-dimensional electrical source numerical simulation method based on coulomb specification |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114065585B (en) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114297905A (en) * | 2022-03-10 | 2022-04-08 | 中南大学 | Quick numerical simulation method of two-dimensional earth electromagnetic field |
| CN114781307A (en) * | 2022-06-17 | 2022-07-22 | 北京智芯仿真科技有限公司 | Non-uniform sampling method and device for integrated circuit Hankel transform filter |
| CN115017782A (en) * | 2022-08-08 | 2022-09-06 | 中南大学 | Three-dimensional natural source electromagnetic field calculation method considering medium anisotropy |
| CN116244877A (en) * | 2022-09-05 | 2023-06-09 | 中南大学 | Three-dimensional magnetic field numerical simulation method and system based on 3D AS-FT |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2021068527A1 (en) * | 2019-10-12 | 2021-04-15 | 中南大学 | Numerical simulation method for footprint-guided high-efficiency airborne electromagnetic survey |
| CN113051775A (en) * | 2021-04-13 | 2021-06-29 | 中南大学 | Horizontal directional drilling track optimization method based on improved radial movement algorithm |
| CN113051779A (en) * | 2021-05-31 | 2021-06-29 | 中南大学 | Numerical simulation method of three-dimensional direct-current resistivity method |
-
2021
- 2021-11-22 CN CN202111386036.4A patent/CN114065585B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2021068527A1 (en) * | 2019-10-12 | 2021-04-15 | 中南大学 | Numerical simulation method for footprint-guided high-efficiency airborne electromagnetic survey |
| CN113051775A (en) * | 2021-04-13 | 2021-06-29 | 中南大学 | Horizontal directional drilling track optimization method based on improved radial movement algorithm |
| CN113051779A (en) * | 2021-05-31 | 2021-06-29 | 中南大学 | Numerical simulation method of three-dimensional direct-current resistivity method |
Non-Patent Citations (2)
| Title |
|---|
| YUGUO LI等: "Finite element modelling of marine controlled-source electromagnetic responses in two-dimensional dipping anisotropic conductivity structures", IEEE, 31 May 2011 (2011-05-31), pages 622 - 636 * |
| 姚大为;王大勇;王刚;朱威;李永博;张振宇;: "CSAMT 2.5维有限元数值模拟", 物探化探计算技术, no. 03, 15 May 2016 (2016-05-15), pages 297 - 306 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114297905A (en) * | 2022-03-10 | 2022-04-08 | 中南大学 | Quick numerical simulation method of two-dimensional earth electromagnetic field |
| CN114297905B (en) * | 2022-03-10 | 2022-06-03 | 中南大学 | Quick numerical simulation method of two-dimensional earth electromagnetic field |
| CN114781307A (en) * | 2022-06-17 | 2022-07-22 | 北京智芯仿真科技有限公司 | Non-uniform sampling method and device for integrated circuit Hankel transform filter |
| CN114781307B (en) * | 2022-06-17 | 2022-08-23 | 北京智芯仿真科技有限公司 | Non-uniform sampling method and device for integrated circuit Hankel transform filter |
| CN115017782A (en) * | 2022-08-08 | 2022-09-06 | 中南大学 | Three-dimensional natural source electromagnetic field calculation method considering medium anisotropy |
| CN116244877A (en) * | 2022-09-05 | 2023-06-09 | 中南大学 | Three-dimensional magnetic field numerical simulation method and system based on 3D AS-FT |
| CN116244877B (en) * | 2022-09-05 | 2023-11-14 | 中南大学 | Three-dimensional magnetic field numerical simulation method and system based on 3D Fourier transform |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114065585B (en) | 2024-05-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114065585A (en) | Three-dimensional electrical source numerical simulation method based on coulomb specification | |
| CN106980736B (en) | A kind of ocean controllable source electromagnetic method finite element forward modeling method of anisotropic medium | |
| CN114065586B (en) | Three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method | |
| CN113158527B (en) | Method for calculating frequency domain electromagnetic field based on implicit FVFD | |
| CN111079278B (en) | Processing method for three-dimensional time domain hybridization discontinuous Galerkin method with additional electromagnetic source item | |
| CN112347687B (en) | Self-adaptive degree-of-freedom electromagnetic-temperature multi-physical-field coupling analysis method | |
| CN113792445B (en) | Three-dimensional magnetotelluric numerical simulation method based on integral equation method | |
| Christophe et al. | An implicit hybridized discontinuous Galerkin method for the 3D time-domain Maxwell equations | |
| CN105717547A (en) | Anisotropy medium magnetotelluric meshless value simulation method | |
| CN113051779B (en) | Numerical simulation method of three-dimensional direct-current resistivity method | |
| CN108197389A (en) | Quick, the high resolution numerical simulation method in two-dimentional ferromagnetic magnetic field | |
| CN110852025A (en) | Three-dimensional electromagnetic slow diffusion numerical simulation method based on hyperconvergence interpolation approximation | |
| CN115292973B (en) | Arbitrarily sampled space wave number domain three-dimensional magnetic field numerical simulation method and system | |
| Liu et al. | Three-dimensional magnetotellurics modeling using edgebased finite-element unstructured meshes | |
| CN111638556A (en) | Magnetotelluric forward modeling method and device based on geospatial solution strategy and storage medium | |
| CN114036745A (en) | Anisotropic medium three-dimensional magnetotelluric forward modeling method and system | |
| CN103018729A (en) | Method for calculating radar scattering cross section of metal cylindrical calibration body | |
| CN113868863A (en) | Statistical electromagnetic DGTD (differential gradient time division) calculation method for uncertain structure target | |
| CN111339688B (en) | Method for solving rocket simulation model time domain equation based on big data parallel algorithm | |
| Sarakorn | 2-D magnetotelluric modeling using finite element method incorporating unstructured quadrilateral elements | |
| Li et al. | A modified dual-level fast multipole boundary element method for large-scale three-dimensional potential problems | |
| Garza et al. | High-order Chebyshev-based Nyström methods for electromagnetics | |
| CN115017782B (en) | Three-dimensional natural source electromagnetic field calculation method considering medium anisotropy | |
| CN114970289A (en) | Three-dimensional magnetotelluric anisotropy forward modeling numerical simulation method, equipment and medium | |
| CN104778293A (en) | Volume integral Nystrom analysis method of inhomogeneous medium target electromagnetic scattering |
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 | ||
| CB03 | Change of inventor or designer information |
Inventor after: Dai Shikun Inventor after: Chen Qingrui Inventor before: Chen Qingrui Inventor before: Dai Shikun |
|
| CB03 | Change of inventor or designer information | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |