CN114065585B - 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
- CN114065585B CN114065585B CN202111386036.4A CN202111386036A CN114065585B CN 114065585 B CN114065585 B CN 114065585B CN 202111386036 A CN202111386036 A CN 202111386036A CN 114065585 B CN114065585 B CN 114065585B
- Authority
- CN
- China
- Prior art keywords
- space
- field
- domain
- wave number
- electric field
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 64
- 238000004088 simulation Methods 0.000 title claims abstract description 23
- 230000005684 electric field Effects 0.000 claims abstract description 85
- 238000004364 calculation method Methods 0.000 claims abstract description 38
- 230000008569 process Effects 0.000 claims abstract description 7
- 230000002159 abnormal effect Effects 0.000 claims description 11
- 125000003275 alpha amino acid group Chemical group 0.000 claims description 6
- 150000001875 compounds Chemical class 0.000 claims description 6
- 230000035699 permeability Effects 0.000 claims description 6
- 238000011160 research Methods 0.000 claims description 5
- 230000005672 electromagnetic field Effects 0.000 abstract description 12
- 230000009466 transformation Effects 0.000 abstract description 5
- 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 15
- 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 4
- 238000004458 analytical method Methods 0.000 description 3
- 230000007547 defect Effects 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
- 238000004891 communication Methods 0.000 description 1
- 230000008094 contradictory effect Effects 0.000 description 1
- 238000013461 design Methods 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
- 230000010354 integration Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000005259 measurement Methods 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
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 three-dimensional electrical source numerical simulation method based on coulomb specification, which is characterized in that when background field calculation is carried out, lamellar coefficients in bessel integral are calculated in advance, fourier transformation 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 tight operator is adopted to correct an electric field in the solving process. The method is based on a secondary field method, does not need to carry out divergence correction, utilizes the advantage, utilizes the coulomb specification to convert the Maxwell equation set into a control equation of a coulomb scalar bit vector position, adopts a quick method when calculating a background field, and utilizes a Fourier transformation technology when calculating the secondary field, thereby greatly improving the calculation efficiency of the numerical simulation of the three-dimensional electric source electromagnetic field. The invention 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 specification.
Background
The artificial source electromagnetic method can overcome the defect of weak signal of the magnetotelluric field source, has strong anti-interference capability, and is widely applied to the fields of mine resources, natural gas, geothermal and hydrologic environmental engineering exploration and the like. The analysis and solution of the electric source electromagnetic field in the complex medium are difficult to obtain, so the frequency domain three-dimensional forward modeling mainly takes numerical simulation as a main part, and the current commonly used three-dimensional forward modeling method mainly takes a finite difference method (Streich R.2009.3D finite-difference frequency-domain modeling ofcontrolled-source electromagnetic data:direct solution and optimization for high accuracy.Geophysics,74(5):
F95-f105.), finite element method (Liu Ying, li Yuguo, han Bo.2017. Controllable source electromagnetic field two-dimensional adaptive vector finite element forward modeling. Geophysical journal, 60 (12): 4874-4886), finite volume method (Peng Ronghua, hu Xiangyun, han Bo, etc. 2016. Frequency domain controlled source three-dimensional forward computation based on the mimicry finite volume method. Geophysical journal, 59 (10): 3927-3939.) and integral equation methods (Shang Jingtian, zhou Feng, ren Zhengyong, et al 2018. Controlled source electromagnetic integral equation forward for complex subsurface anomalies. Geophysical journal, 61 (4): 1549-1562). In recent years, the efficiency of three-dimensional forward modeling is still a research hot spot, and the development of high-efficiency and high-precision three-dimensional electromagnetic method forward modeling technology of an electrical source is still quite large.
Disclosure of Invention
Aiming at the defects in the prior art, the invention provides a three-dimensional electric source numerical simulation method based on coulomb specification, which solves the problems of large calculation amount, high storage requirement, low calculation efficiency and the like in the analysis process of an electric source electromagnetic field in a complex medium in the prior art.
In order to achieve the above purpose, the present invention provides a coulomb specification-based three-dimensional electrical source numerical simulation method, which comprises the following steps:
Step 1, carrying out hexahedral subdivision on a research area along the x, y and z directions, wherein the directions of the x and y are uniformly subdivided, the directions of the z are uniformly or non-uniformly subdivided, and the number of subdivision nodes in the directions of the x, y and z is N x、Ny、Nz respectively;
Step 2, obtaining parameters of a ground model, including parameters of a background layer, calculation frequency and abnormal conductivity, and obtaining a background electric field and a background magnetic field based on the calculation frequency, the position of a source and the given parameters of the background layer;
step 3, converting the Maxwell equation set into equations related to vector bits and scalar bits by adopting 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 scattering current of a right end item in the control equation;
Step 4, carrying out horizontal two-dimensional Fourier transform on the control equation to obtain a normal differential equation of a vector bit scalar position of a space-wave number domain, adopting Gaussian offset wave numbers, weighting the normal differential equation under each wave number by Galerkin's method, and solving by adopting finite elements to obtain a vector bit scalar position of a space-wave number domain secondary field under each wave number;
Step 5, obtaining a space-wave number domain secondary electric field based on the relation between the space-wave number domain secondary field vector bit scalar bit and the electric field in the space-wave number domain, and then obtaining the space domain secondary electric field by utilizing horizontal two-dimensional inverse Fourier transform;
Step 6, adding the background electric field and the space domain secondary electric field to obtain a space domain total electric field, and judging whether the space domain total electric field meets a given iteration convergence condition:
If so, obtaining a space-wave number domain secondary magnetic field based on the space-wave number domain secondary field vector bit scalar bit, obtaining a space domain secondary magnetic field after horizontal two-dimensional inverse Fourier transform, 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;
If not, modifying the space domain total field by adopting a compact operator, multiplying the space domain total field by the abnormal conductivity to obtain a new scattering current, and repeating the steps 4-6 after updating the scattering current of the right end term in the control equation.
As a further improvement of the above scheme, in 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, which specifically includes the following steps:
The bessel integral is in the form of:
Setting a background field to have N layers of media, wherein the background field is a background electric field or a background magnetic field; c l and D l in formula (1) are coefficients in each layered medium of the background field; m l is an exponential coefficient, k l is an amount related to frequency and conductivity, permittivity and permeability of each layer, l=1, 2,3 … N; j i is the i-order bessel function; /(I) R is the horizontal distance, (x, y) is the calculated point horizontal coordinate, (x s,ys) is the source point horizontal coordinate; m is an integral variable; z is the z-direction coordinate of the calculated point;
When the background field is layered into N layers, the coefficient C l、Dl of each layer is calculated, and when bessel integral is calculated, the coefficient corresponding to each node is multiplied into the integral expression.
As a further improvement of the above solution, in step 3, the control equation satisfied by the secondary field vector bit vector bits is:
Of the formula (I) The vector bit of the secondary field in the x, y and z directions, and phi s is the scalar bit of the secondary field; k b is the background field wave number, i is the imaginary unit, ω is the angular frequency, ω=2pi f, f is the calculated frequency, μ 0 is the vacuum permeability; /(I) For background admittance, σ b for background conductivity, ε for dielectric constant; The scattering currents in the x, y and z directions are σ a, the abnormal conductivity is shown, and Ex, ey and Ez are the total electric fields in the x, y and z directions.
As a further improvement of the above scheme, in step 3, the process of obtaining the space-wave number domain quadric field vector bit scalar bit is:
performing horizontal two-dimensional Fourier transform on the formula (2) to obtain:
In the formula (3), the symbol on the variable represents a space-wave number domain variable, the symbol without the symbol represents a space domain variable, and kx and ky are wave numbers corresponding to x and y directions in Fourier transform;
The Gaussian offset wave numbers are adopted, the ordinary differential equation under each wave number is weighted by the Galerkin method, and a finite element method equation set is formed by:
In the formula (4), the symbol on the variable represents a space-wave number domain variable, the symbol without the symbol represents a space domain variable, N e is a subdivision unit in the z direction, N i is a Galao Jin Quan function, the relation N z＝2×Ne+1,kp is wave number, and N z is a normal vector; and (3) loading an uplink boundary condition on the upper boundary and loading a downlink boundary condition on the lower boundary, and solving an equation set (4) to obtain the vector bit scalar of the space-wave number domain secondary field.
As a further improvement of the above scheme, step 5, the relationship between the space-wave number domain secondary field vector bit scalar bit and the electric field in the space-wave number domain is:
In the formula (5), the amino acid sequence of the compound, Is a secondary electric field in the x direction of the space-wave number domain,/>Is a secondary electric field in the y direction of the space-wave number domain,/>Is a space-wave number domain z-direction secondary electric field.
As a further improvement of the above scheme, step 6, the relationship between the vector bit scalar of the space-wave number domain secondary field and the space-wave number domain secondary magnetic field is:
in the formula (6), the amino acid sequence of the compound, Is a secondary magnetic field in the x direction of the space-wave number domain,/>Is a secondary magnetic field in the y direction of the space-wave number domain,Is a space-wave number domain z-direction secondary magnetic field.
As a further improvement of the above scheme, in step 6, the format of the compact operator corrected total electric field is:
In equation (7), E (n)(xi,yj,zk on the left of the equation) is the corrected spatial domain total electric field, E (n-1)(xi,yj,zk on the right of the equation is the last corrected spatial domain total electric field and E (n)(xi,yj,zk which is the new spatial domain total electric field calculated by the current forward modeling.
As a further improvement of the above solution, in step 6, the iterative convergence condition of the total electric field in the spatial domain is:
In the formula, if the formula is satisfied, the iteration convergence condition is satisfied, otherwise, the iteration convergence condition is not satisfied; err is error, E (n)(xi,yj,zk) is the new spatial domain total electric field obtained by solving the equation set for the nth time of the node with the coordinates of (x i,yj,zk) and adding the background electric field, E (n-1)(xi,yj,zk) is the spatial domain total electric field obtained by last calculation of the node with the coordinates of (x i,yj,zk), and when n=1, E (0)(xi,yj,zk) is the background electric field of the node with the coordinates of (x i,yj,zk).
According to the three-dimensional electric source numerical simulation method based on the coulomb specification, the divergence correction is not needed based on the secondary field method, the advantage is utilized, the Maxwell equation set is converted into the control equation of the coulomb scalar bit vector position by using the coulomb specification, a rapid filter coefficient method is adopted when a background field is calculated, and a Fourier transformation technology is utilized when the secondary field is calculated, so that the calculation efficiency of the 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 that are required in the embodiments or the description of the prior art will be briefly described, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and other drawings may be obtained according to the structures shown in these drawings without inventive effort for a person skilled in the art.
FIG. 1 is a flow chart of a method for simulating three-dimensional electrical source values based on coulomb specifications according to an embodiment of the invention;
FIG. 2 is a schematic view 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 a comparison of an electric field and a magnetic field of the method and the intel 3D software calculation line y=0m according to the embodiment of the present invention, where (a) is an electric field E x component, (b) is an electric field E y component, (c) is an electric field E z component, (D) is a magnetic field H x component, (E) is a magnetic field H y component, and (f) is a magnetic field H z component.
The achievement of the objects, functional features and advantages of the present invention will be further described with reference to the accompanying drawings, in conjunction with the embodiments.
Detailed Description
The following description of the embodiments of the present invention will be made clearly and fully with reference to the accompanying drawings, in which it is evident that the embodiments described are only some, but not all embodiments of the invention. All other embodiments, which can be made by those skilled in the art based on the embodiments of the invention without making any inventive effort, are intended to be within the scope of the invention.
It should be noted that all directional indicators (such as up, down, left, right, front, and rear … …) in the embodiments of the present invention are merely used to explain the relative positional relationship, movement, etc. between the components in a particular posture (as shown in the drawings), and if the particular posture is changed, the directional indicator is changed accordingly.
Furthermore, descriptions such as those referred to as "first," "second," and the like, are provided for descriptive purposes only and are not to be construed as indicating or implying a relative importance or implying an order of magnitude of the indicated technical features in the present disclosure. Thus, a feature defining "a first" or "a second" may explicitly or implicitly include at least one such feature. In the description of the present invention, the meaning of "plurality" means at least two, for example, two, three, etc., unless specifically defined otherwise.
In the present invention, unless specifically stated and limited otherwise, the terms "connected," "affixed," and the like are to be construed broadly, and for example, "affixed" may be a fixed connection, a removable connection, or an integral body; the device can be mechanically connected, electrically connected, physically connected or wirelessly connected; either directly or indirectly, through intermediaries, or both, may be in communication with each other or in interaction with each other, unless expressly defined otherwise. The specific meaning of the above terms in the present invention can be understood by those of ordinary skill in the art according to the specific circumstances.
In addition, the technical solutions of the embodiments of the present invention may be combined with each other, but it is necessary to be based on the fact that those skilled in the art can implement the technical solutions, and when the technical solutions are contradictory or cannot be implemented, the combination of the technical solutions should be considered as not existing, and not falling within the scope of protection claimed by the present invention.
The embodiment discloses a three-dimensional electrical source numerical simulation method based on coulomb specification, which calculates a laminar coefficient in bessel integral in advance when calculating a background field, so as to reduce the solution of the coefficient and improve the calculation efficiency. And Fourier transformation is adopted when solving the secondary electromagnetic field, 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 calculated quantity and the storage requirement are greatly reduced, and the method has higher calculation efficiency. And finally, correcting an electric field in the solving process by adopting a tight operator, and stably converging an algorithm, so that the method is applicable to any model with conductivity other than 0. After the background field calculation and the secondary field calculation are processed, the rapid and accurate numerical simulation of the three-dimensional electric source electromagnetic field is realized, a new technology is provided for the high-efficiency and high-precision numerical simulation forward modeling of the large-scale electric source, and the method has important invention value for three-dimensional high-efficiency fine inversion. Referring to fig. 1, the method includes the following steps 1-6.
And step 1, carrying out hexahedral subdivision on the research area along the x, y and z directions to obtain grid node coordinates of the research area. The split nodes in the x direction and the y direction are uniformly split and non-uniformly split in the z direction, and the number of split nodes in the x direction, the y direction and the z direction is N x、Ny、Nz respectively.
Step 2, obtaining parameters of the ground model, including background layer parameters σ b(zk), calculation frequency f, abnormal conductivity parameters σ a(xi,yj,zk),(xi,yj,zk), and representing coordinates of split nodes with numbers (i, j, k), i=1, 2, … N x,j＝1,2,…,Ny,k＝1,2,…,Nz. And deriving background electric and magnetic fields based on the calculated frequency, the source location, and given background layer parameters.
The analysis solution of the background electromagnetic field has bessel integral terms, the bessel integral terms in the background electric field and the background magnetic field are obtained based on a rapid filter coefficient method, and the bessel integral form is as follows:
Setting a background field to have N layers of media, wherein the background field is a background electric field or a background magnetic field; c l and D l in formula (1) are coefficients in each layered medium of the background field; m l is an exponential coefficient, k l is an amount related to frequency and conductivity, permittivity and permeability of each layer, l=1, 2,3 … N; j i is the i-order bessel function; /(I) R is the horizontal distance, (x, y) is the calculated point horizontal coordinate, (x s,ys) is the source point horizontal coordinate; 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, see in detail literature (Cai Cheng.2014. Fast Hanker's transform and its use in forward calculations. Geophysical progress, 29 (03): 1384-1390.).
Since z and r are different for each node, the conventional filter coefficient method calculates the integral once for each node, C l and D l. When the node numbers in the x, y and z directions are N x、Ny、Nz, the number of calculation times of the layer coefficient is N x×Ny×Nz. In this embodiment, considering the layering property of the background field, when the background field is layered into N layers, the coefficient C l、Dl of each layer is calculated first, and when bessel integration is calculated, the coefficient corresponding to each node is multiplied into the integral expression, and at this time, the coefficient of each layer is calculated only N x×Ny ×n times. In three-dimensional electromagnetic numerical simulation, the number of layers generally considered is rarely more than 10. In the three-dimensional forward modeling, the number of nodes calculated in the vertical direction is generally not less than 10 nodes. At this time, the rapid filter coefficient algorithm can reduce the number of calculation of the coefficients, thereby improving the calculation efficiency.
Finally calculating the background field to obtain the background electric field of each split nodeBackground magnetic field/>Where (x i,yj,zk) denotes the split node coordinates numbered (i, j, k), i=1, 2, … N x,j＝1,2,…,Ny,k＝1,2,…,Nz.
And 3, converting the Maxwell equation set into equations related to vector bits and scalar bits by adopting coulomb specifications, 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 scattering current of a right end item in the control equation. The control equation for the scalar bit of the secondary field vector bit is:
Of the formula (I) The vector bit of the secondary field in the x, y and z directions, and phi s is the scalar bit of the secondary field; k b is the background field wave number, i is the imaginary unit, ω is the angular frequency, ω=2pi f, f is the calculated frequency, μ 0 is the vacuum permeability; /(I) For background admittance, σ b for background conductivity, ε for dielectric constant; The scattering currents in the x, y and z directions are σ a, the abnormal conductivity is shown, and Ex, ey and Ez are the total electric fields in the x, y and z directions.
And 4, performing horizontal two-dimensional Fourier transformation on the control equation to obtain a normal differential equation of the vector bit scalar position of the space-wave number domain, weighting the normal differential equation under each wave number by using a Galerkin method by using Gaussian offset wave numbers, and solving by using finite elements to obtain the vector bit scalar position of the space-wave number domain secondary field under each wave number. The space-wave number domain secondary field vector bit scalar bit acquisition process comprises the following steps:
performing horizontal two-dimensional Fourier transform on the formula (2) to obtain:
in the expression (3), the symbol-on-variable represents a space-wave number domain variable, the symbol-not represents a space domain variable, and kx and ky are wave numbers corresponding to x and y directions in fourier transform.
In this embodiment, the gaussian offset wave number is used, and the gaussian points are 4 gaussian points. The selection of the wave number is related to the number of sampling points and sampling intervals, and the selection rule of the standard wave number is described as follows: when the sampling number is M and the sampling interval is Δx, the fundamental wave Δk=2pi/(m×Δx), and the number of wave numbers k p are also M.
When M is even, the selection rule is:
when M is an odd number, the wavenumber k p is selected as:
The Gaussian offset wavenumber is selected as follows:
kp＝(p+tg)·Δk,g＝1,2,3,4
The 4 gaussian points t g and the corresponding coefficients a g are respectively:
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 at this time was 4M.
The ordinary differential equation at each wave number is weighted by Galerkin's method, and the finite element method equation set is formed as follows:
In the formula (4), the symbol on the variable represents a space-wave number domain variable, the symbol on the variable is not a space domain variable, N e is a split unit in the z direction, N i is a galao Jin Quan function, and the embodiment adopts a quadratic interpolation shape function. The existence relation N z＝2×Ne+1,kp is the wave number, and N z is the normal vector; and (3) loading an uplink boundary condition on the upper boundary and loading a downlink boundary condition on the lower boundary, and solving an equation set (4) to obtain the vector bit scalar of the space-wave number domain secondary field.
And 5, obtaining a space-wave number domain secondary electric field based on a relation between the space-wave number domain secondary field vector bit scalar bit and the electric field in the space-wave number domain, and obtaining the space domain secondary electric field by utilizing horizontal two-dimensional inverse Fourier transform. Wherein, the relation between the vector bit scalar of the secondary field in the space-wave number domain and the electric field in the space-wave number domain is as follows:
In the formula (5), the amino acid sequence of the compound, Is a secondary electric field in the x direction of the space-wave number domain,/>Is a secondary electric field in the y direction of the space-wave number domain,/>The space-wave number domain secondary electric field is obtained by using the formula (5) as the space-wave number domain z-direction secondary electric field.
And 6, adding the background electric field and the space domain secondary electric field to obtain a space domain total electric field, and judging whether the space domain total electric field meets a given iteration convergence condition. The iterative convergence condition of the total electric field in the space domain is as follows:
In the formula, if the formula is satisfied, the iteration convergence condition is satisfied, otherwise, the iteration convergence condition is not satisfied; err is error, E (n)(xi,yj,zk) is the new spatial domain total electric field obtained by solving the equation set for the nth time of the node with the coordinates of (x i,yj,zk) and adding the background electric field, E (n-1)(xi,yj,zk) is the spatial domain total electric field obtained by last calculation of the node with the coordinates of (x i,yj,zk), and when n=1, E (0)(xi,yj,zk) is the background electric field of the node with the coordinates of (x i,yj,zk).
If so, obtaining a space-wave number domain secondary magnetic field based on a space-wave number domain secondary field vector bit scalar bit, obtaining a space domain secondary magnetic field after horizontal two-dimensional inverse Fourier transform, 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-wave number domain secondary field vector bit scalar bit and the space wave number domain secondary magnetic field is as follows:
in the formula (6), the amino acid sequence of the compound, Is a secondary magnetic field in the x direction of the space-wave number domain,/>Is a secondary magnetic field in the y direction of the space-wave number domain,Is a space-wave number domain z-direction secondary magnetic field.
If not, modifying the space domain total field by a compact operator, multiplying the space domain total field by the abnormal conductivity to obtain a new scattering current, and repeating the step 4-6 after updating the scattering current of the right end item in the control equation, wherein the format of the compact operator for modifying the total electric field is as follows:
In the formula (7), E (n)(xi,yj,zk) on the left side of the equation is the corrected total electric field of the spatial domain, E (n-1)(xi,yj,zk) on the right side of the equation is the corrected total electric field of the spatial domain of the last time and E (n)(xi,yj,zk) of the new spatial domain obtained by the forward calculation And obtaining a new scattering current, and then carrying out a new round of solving.
The accuracy and efficiency of the three-dimensional electrical source numerical simulation based on the coulomb specification provided by the embodiment are verified by combining with a design model example, and the tested computer is an Intel (R) Core (TM) i7-6700HQ CPU main frequency of 2.60GHz, and a memory of 16GB and 64-bit win10 system.
(1) Verification of correctness
As shown in FIG. 2, the model XOY plane projection is characterized in that the background is a uniform half-space medium, the upper half space is air, the air conductivity sigma 0＝10-12 S/m, the background conductivity of the lower half space is sigma b =0.01S/m, the frequency is 1Hz, the calculation range x-1000 m, the y-1000 m and the z-1000 m are respectively set, the number of split grid nodes is 101 x 101, the three directions are uniformly split, the split interval deltax in the x direction and the split interval deltay in the y direction are 20m, and the split interval deltaz in the z direction is 10m. The range of the abnormal body is-100 m in the x direction, -200 m in the y direction and 400-800 m in the z direction, and the conductivity sigma=0.1S/m of the abnormal body. The electrical source is 200m long along the x direction, the center coordinates of the source are (0, -2000,0.01), and the current is 10A. The accuracy of the method is verified by using the calculation result of three-dimensional forward software INTEM3D based on an integral equation method developed by university of Utah in the United states as a reference.
Fig. 3 is a comparison of the electric field and the magnetic field of the measurement line y=0m calculated by the method of the present embodiment and the intel 3D software, 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 of software INTEM3D, and as can be seen from fig. 3, the two numerical solutions are well matched, and zero value points and E z have large deviation due to small values, but the precision meets the requirement. In conclusion, the electromagnetic field has small calculation error, the correctness of the method of the embodiment is verified, and the method of the embodiment has high calculation accuracy.
(2) Primary and secondary field calculation efficiency
When the electromagnetic fields of 101 x 101 sections are calculated by using a traditional filter coefficient algorithm and a rapid filter coefficient algorithm respectively, the calculation time of the traditional algorithm is 70.3s, and the occupied memory is 102.5MB; the calculation time of the fast algorithm is 39.8s, and the occupied memory is 120.1MB; under the condition of almost the memory, the calculation time is 2/5 faster, so that the fast algorithm has high calculation speed and high calculation efficiency.
For the secondary field calculation, when the simulation calculation of the electric source electromagnetic field value of the embodiment reaches the expected convergence relative residual error of 10 -4, the iteration is needed for 12 times, the forward time of 101 x 101 nodes is 100.35s, and the memory 2034MB is occupied. The calculation time of the background field is 39.8s, and the total time spent by forward modeling is 140s, so that the algorithm is high in calculation speed, small in occupied memory and high in calculation efficiency.
The foregoing description is only of the preferred embodiments of the present invention and is not intended to limit the scope of the invention, and all equivalent structural changes made by the description of the present invention and the accompanying drawings or direct/indirect application in other related technical fields are included in the scope of the invention.
Claims (8)
1. The three-dimensional electrical source numerical simulation method based on coulomb specification is characterized by comprising the following steps:
Step 1, carrying out hexahedral subdivision on a research area along the x, y and z directions, wherein the directions of the x and y are uniformly subdivided, the directions of the z are uniformly or non-uniformly subdivided, and the number of subdivision nodes in the directions of the x, y and z is N x、Ny、Nz respectively;
Step 2, obtaining parameters of a ground model, including parameters of a background layer, calculation frequency and abnormal conductivity, and obtaining a background electric field and a background magnetic field based on the calculation frequency, the position of a source and the given parameters of the background layer;
step 3, converting the Maxwell equation set into equations related to vector bits and scalar bits by adopting 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 scattering current of a right end item in the control equation;
Step 4, carrying out horizontal two-dimensional Fourier transform on the control equation to obtain a normal differential equation of a vector bit scalar position of a space-wave number domain, adopting Gaussian offset wave numbers, weighting the normal differential equation under each wave number by Galerkin's method, and solving by adopting finite elements to obtain a vector bit scalar position of a space-wave number domain secondary field under each wave number;
Step 5, obtaining a space-wave number domain secondary electric field based on the relation between the space-wave number domain secondary field vector bit scalar bit and the electric field in the space-wave number domain, and then obtaining the space domain secondary electric field by utilizing horizontal two-dimensional inverse Fourier transform;
Step 6, adding the background electric field and the space domain secondary electric field to obtain a space domain total electric field, and judging whether the space domain total electric field meets a given iteration convergence condition:
If so, obtaining a space-wave number domain secondary magnetic field based on the space-wave number domain secondary field vector bit scalar bit, obtaining a space domain secondary magnetic field after horizontal two-dimensional inverse Fourier transform, 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;
If not, modifying the space domain total field by adopting a compact operator, multiplying the space domain total field by the abnormal conductivity to obtain a new scattering current, and repeating the steps 4-6 after updating the scattering current of the right end term in the control equation.
2. The method for simulating three-dimensional electrical source values based on coulomb specifications according to claim 1, wherein in 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 in the form of:
Setting a background field to have N layers of media, wherein the background field is a background electric field or a background magnetic field; c l and D l in formula (1) are coefficients in each layered medium of the background field; m l is an exponential coefficient, k l is an amount related to frequency and conductivity, permittivity and permeability of each layer, l=1, 2,3 … N; j i is the i-order bessel function; /(I) R is the horizontal direction distance, (x, y) is the calculated point horizontal direction coordinate, (x s,ys) is the source point horizontal direction coordinate; m is an integral variable; z is the vertical direction coordinate of the calculation point;
When the background field is layered into N layers, the coefficient C l、Dl of each layer is calculated, and when bessel integral is calculated, the coefficient corresponding to each node is multiplied into the integral expression.
3. The method for simulating three-dimensional electrical source numerical simulation based on coulomb specification according to claim 1, wherein in step 3, the control equation for the scalar bit of the secondary field vector bit is:
Of the formula (I) The vector bit of the secondary field in the x, y and z directions, and phi s is the scalar bit of the secondary field; /(I)K b is the background field wave number, i is the imaginary unit, ω is the angular frequency, ω=2pi f, f is the calculated frequency, μ 0 is the vacuum permeability; For background admittance, σ b for background conductivity, ε for dielectric constant; /(I) The scattering currents in the x, y and z directions are σ a, the abnormal conductivity is shown, and Ex, ey and Ez are the total electric fields in the x, y and z directions.
4. The method for simulating three-dimensional electrical source numerical simulation based on coulomb specification according to claim 3, wherein in step 3, the process of obtaining the vector bit scalar of the space-wave number domain secondary field is as follows:
performing horizontal two-dimensional Fourier transform on the formula (2) to obtain:
In the formula (3), the symbol on the variable represents a space-wave number domain variable, the symbol without the symbol represents a space domain variable, and kx and ky are wave numbers corresponding to x and y directions in Fourier transform;
The Gaussian offset wave numbers are adopted, the ordinary differential equation under each wave number is weighted by the Galerkin method, and a finite element method equation set is formed by:
In the formula (4), the symbol on the variable represents a space-wave number domain variable, the symbol without the symbol represents a space domain variable, N e is a subdivision unit in the z direction, N i is a Galao Jin Quan function, the relation N z＝2×Ne+1,kp is wave number, and N z is a normal vector; and (3) loading an uplink boundary condition on the upper boundary and loading a downlink boundary condition on the lower boundary, and solving an equation set (4) to obtain the vector bit scalar of the space-wave number domain secondary field.
5. The method of three-dimensional electrical source numerical simulation based on coulomb specification as recited in claim 4, wherein in step 5, the relationship between the vector bit position of the secondary field in the space-wave number domain and the electric field in the space-wave number domain is:
In the formula (5), the amino acid sequence of the compound, Is a secondary electric field in the x direction of the space-wave number domain,/>Is a secondary electric field in the y direction of the space-wave number domain,/>Is a space-wave number domain z-direction secondary electric field.
6. The method for simulating three-dimensional electrical source values based on coulomb specifications according to claim 4, wherein in step 6, the relationship between the vector bit vector bits of the secondary field in the space-wavenumber domain and the secondary magnetic field in the space-wavenumber domain is:
in the formula (6), the amino acid sequence of the compound, Is a secondary magnetic field in the x direction of the space-wave number domain,/>Is a secondary magnetic field in the y direction of the space-wave number domain,/>Is a space-wave number domain z-direction secondary magnetic field.
7. The method for modeling three-dimensional electrical source values based on coulomb specifications as defined in claim 4, wherein in step 6, the format of the compact operator modified total electric field is:
In equation (7), E (n)(xi,yj,zk on the left of the equation) is the corrected spatial domain total electric field, E (n-1)(xi,yj,zk on the right of the equation is the last corrected spatial domain total electric field and E (n)(xi,yj,zk which is the new spatial domain total electric field calculated by the current forward modeling.
8. The method for three-dimensional electrical source numerical simulation based on coulomb specification according to any one of claims 1 to 7, wherein in step 6, the iterative convergence condition of the total electric field in the spatial domain is:
In the formula, if the formula is satisfied, the iteration convergence condition is satisfied, otherwise, the iteration convergence condition is not satisfied; err is error, E (n)(xi,yj,zk) is the new spatial domain total electric field obtained by solving the equation set for the nth time of the node with the coordinates of (x i,yj,zk) and adding the background electric field, E (n-1)(xi,yj,zk) is the spatial domain total electric field obtained by last calculation of the node with the coordinates of (x i,yj,zk), and when n=1, E (0)(xi,yj,zk) is the background electric field of the node with the coordinates of (x i,yj,zk).
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 CN114065585A (en) | 2022-02-18 |
| CN114065585B true 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) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114297905B (en) * | 2022-03-10 | 2022-06-03 | 中南大学 | Quick numerical simulation method of two-dimensional earth electromagnetic field |
| CN114781307B (en) * | 2022-06-17 | 2022-08-23 | 北京智芯仿真科技有限公司 | Non-uniform sampling method and device for integrated circuit Hankel transform filter |
| CN115017782B (en) * | 2022-08-08 | 2022-12-16 | 中南大学 | Three-dimensional natural source electromagnetic field calculation method considering medium anisotropy |
| CN116244877B (en) * | 2022-09-05 | 2023-11-14 | 中南大学 | Three-dimensional magnetic field numerical simulation method and system based on 3D Fourier transform |
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 |
|---|
| CSAMT 2.5维有限元数值模拟;姚大为;王大勇;王刚;朱威;李永博;张振宇;;物探化探计算技术;20160515(第03期);第297-306页 * |
| Finite element modelling of marine controlled-source electromagnetic responses in two-dimensional dipping anisotropic conductivity structures;Yuguo Li等;IEEE;20110531;第622-636页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114065585A (en) | 2022-02-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114065585B (en) | Three-dimensional electrical source numerical simulation method based on coulomb specification | |
| CN106199742B (en) | A kind of dimension of Frequency-domain AEM 2.5 band landform inversion method | |
| Ansari et al. | 3D finite-element forward modeling of electromagnetic data using vector and scalar potentials and unstructured grids | |
| Key et al. | Adaptive finite-element modeling using unstructured grids: The 2D magnetotelluric example | |
| CN106980736B (en) | A kind of ocean controllable source electromagnetic method finite element forward modeling method of anisotropic medium | |
| CN110058315B (en) | Three-dimensional anisotropic radio frequency magnetotelluric adaptive finite element forward modeling method | |
| CN112949134B (en) | Earth-well transient electromagnetic inversion method based on non-structural finite element method | |
| CN114065586B (en) | Three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method | |
| CN107657137A (en) | A kind of unusual diffusion three-dimensional simulation method of the fractional order electromagnetism of Approximation by rational function | |
| CN113792445B (en) | Three-dimensional magnetotelluric numerical simulation method based on integral equation method | |
| CN105717547A (en) | Anisotropy medium magnetotelluric meshless value simulation method | |
| 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 | |
| CN114036745A (en) | Anisotropic medium three-dimensional magnetotelluric forward modeling method and system | |
| CN111638556A (en) | Magnetotelluric forward modeling method and device based on geospatial solution strategy and storage medium | |
| CN113158527A (en) | Method for calculating frequency domain electromagnetic field based on implicit FVFD | |
| Zhou et al. | Three-dimensional finite-element analysis of magnetotelluric data using Coulomb-gauged potentials in general anisotropic media | |
| Huang et al. | 3D anisotropic modeling and identification for airborne EM systems based on the spectral-element method | |
| CN115017782B (en) | Three-dimensional natural source electromagnetic field calculation method considering medium anisotropy | |
| Li et al. | Extension of the regularization technique to controlled-source electromagnetic modeling in general anisotropic conductivity media | |
| CN113627027B (en) | Method and system for simulating electromagnetic field value of non-trivial anisotropic medium | |
| CN112748471B (en) | Gravity-magnetic data continuation and conversion method of unstructured equivalent source | |
| CN113238285B (en) | Resistivity calculation method, system and terminal for geophysical charging method exploration | |
| CN113297526B (en) | Horizontal layered soil structure joint inversion method based on Wenner quadrupole and magnetotelluric data | |
| CN114970289A (en) | Three-dimensional magnetotelluric anisotropy forward modeling numerical simulation method, equipment and medium |
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 |