CN114065586A - Three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method - Google Patents

Three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method Download PDF

Info

Publication number
CN114065586A
CN114065586A CN202111386046.8A CN202111386046A CN114065586A CN 114065586 A CN114065586 A CN 114065586A CN 202111386046 A CN202111386046 A CN 202111386046A CN 114065586 A CN114065586 A CN 114065586A
Authority
CN
China
Prior art keywords
field
dimensional
space
bit
domain
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
Application number
CN202111386046.8A
Other languages
Chinese (zh)
Other versions
CN114065586B (en
Inventor
戴世坤
陈轻蕊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202111386046.8A priority Critical patent/CN114065586B/en
Publication of CN114065586A publication Critical patent/CN114065586A/en
Application granted granted Critical
Publication of CN114065586B publication Critical patent/CN114065586B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • 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)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Measuring Magnetic Variables (AREA)

Abstract

The invention discloses a numerical simulation method of finite elements of a three-dimensional magnetotelluric space-wavenumber domain, which converts a Maxwell equation set into a control equation set which is satisfied by vector bit scalar positions through coulomb regulation and a secondary field method, and further obtains control equations of secondary vector positions and scalar positions. And performing two-dimensional Fourier transform on the control equations of the secondary vector bit and the scalar bit in the horizontal direction, and converting the three-dimensional control equation into a one-dimensional control equation under a plurality of wave numbers. The one-dimensional control equation is fitted by adopting a quadratic interpolation function, and is solved by adopting a one-dimensional finite element method. The three-dimensional finite element control equation is converted into a plurality of one-dimensional problems through two-dimensional Fourier transform, the one-dimensional problems are solved by adopting a one-dimensional finite element method, the calculation speed is high, the storage requirement is low, the efficiency is high, the high-efficiency and high-precision numerical simulation of the three-dimensional magnetotelluric is realized, and an important technical support is provided for the precise inversion of large-scale magnetotelluric measured data. The method is applied to the field of numerical simulation.

Description

Three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method
Technical Field
The invention relates to the technical field of numerical simulation, in particular to a numerical simulation method of a finite element in a three-dimensional magnetotelluric space-wavenumber domain.
Background
The frequency domain electromagnetic method is widely applied to the fields of earth deep detection, mineral oil and gas resource exploration, engineering exploration and the like. The three-dimensional electromagnetic field inversion technology is a final means of quantitative interpretation, and the improvement of the inversion calculation speed and accuracy has an important promoting effect on the detection potential of the method. In the inversion process, the forward algorithm occupies most of the calculation time, so the final foothold for improving the application effect of the method lies in improving the calculation speed of the forward algorithm.
The finite element method has strong adaptability to complex models, and in recent years, researches on a finite element method of nodes (GrayverAV, Burg M,2014, Robust and scalable 3-D geo-electromagnetic modeling using the finite element method, geographic Journal International,198(1), 110-125), a finite element method of vectors (consider the WeChat, Wushi, Stone inkstone; 2020. magnetotelluric rapid three-dimensional forward research, geophysical exploration and chemical exploration, 44(6), 1387-1398) and the like improve the efficiency of the three-dimensional electromagnetic method to a certain extent, but the efficiency still has a great improvement space.
Disclosure of Invention
Aiming at the efficiency problem of the finite element method in magnetotelluric numerical simulation in the prior art, the invention provides the three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method, the three-dimensional finite element control equation is converted into a plurality of one-dimensional problems through two-dimensional Fourier transform, and the one-dimensional problems are solved by adopting the one-dimensional finite element method, so that the calculation speed is high, the storage requirement is low, and the efficiency is high.
In order to achieve the above object, the present invention provides a finite element numerical simulation method for three-dimensional magnetotelluric space-wavenumber domain, comprising the following steps:
step 1, establishing a cuboid target area in a rectangular coordinate system, and dispersing the target area in the directions of x, y and z, wherein a three-dimensional conductivity abnormal body is in the target area;
step 2, setting the polarization direction of a geodetic electromagnetic field source, obtaining a background electric field and a background magnetic field of a sampling node in a target area based on a background field equation set, background conductivity and frequency, and taking the background electric field as an initial electric field;
step 3, multiplying the abnormal conductivity of the sampling node of the target area by the initial electric field to obtain the scattering current of the sampling node, and performing two-dimensional Fourier transform on the scattering current in the horizontal direction to obtain two-dimensional wavenumber domain scattering current;
step 4, substituting the two-dimensional wavenumber domain scattering current into a space-wavenumber domain secondary field control equation set, obtaining a space-wavenumber domain secondary field vector bit and a scalar bit by adopting a one-dimensional finite element method based on secondary interpolation, and obtaining a space-wavenumber domain secondary electric field by utilizing the relation between the space-wavenumber domain secondary field vector bit and the scalar bit and the electric field;
step 5, performing two-dimensional Fourier inverse transformation on the space-wavenumber domain secondary electric field to obtain a space domain secondary electric field, and adding the space domain secondary electric field and a background electric field to obtain a space domain total electric field;
step 6, judging whether the total electric field of the spatial domain meets a given iteration convergence condition:
if the magnetic field is satisfied, obtaining a space-wavenumber domain secondary magnetic field based on the relation between the space-wavenumber domain secondary field and the vector bit scalar bit and the magnetic field, obtaining the space domain secondary magnetic field after the two-dimensional inverse Fourier transform in the horizontal direction, adding a background magnetic field to obtain a space domain total magnetic field, and then outputting the space domain total electric field and the space domain total magnetic field;
if not, taking the total electric field of the spatial domain as an initial electric field, and repeating the step 3-6;
and 7, obtaining an electromagnetic field obtained by forward modeling of the earth electromagnetic field source polarized in the x direction and the y direction based on the steps 2-6, and further obtaining apparent resistivity and phase at the measuring point.
As a further improvement of the above technical solution, in step 1, in the process of discretizing the target area in the x, y, and z directions, the x and y directions are uniform subdivision, the z direction is uniform subdivision or non-uniform subdivision, and the numbers of subdivision nodes in the x, y, and z directions are respectively set as Nx、Ny、Nz。
As a further improvement of the above technical solution, the background field equation set in step 2 and the space-wavenumber domain quadratic field control equation set in step 4 have an acquisition process:
introducing vector bits and scalar bits under the coulomb specification, and converting a Maxwell equation system related to the electromagnetic field into a control equation system related to the vector bits and the scalar bits based on the coulomb specification;
based on a quadratic field method, dividing a control equation set into a background field equation set and a control equation set with a quadratic field vector bit and a scalar bit;
and performing two-dimensional Fourier transform on the control equation set of the secondary field vector bit and the scalar bit in the horizontal direction to obtain a one-dimensional space domain secondary field control equation set under a two-dimensional wavenumber domain, namely a space-wavenumber domain secondary field control equation set.
As a further improvement of the above technical scheme, the time constant is set as eiωtWherein i is an imaginary unit, ω is an angular frequency, ω is 2 pi f, and f is a calculation frequency, and a control equation system satisfied by the vector bit a and the scalar bit Φ based on the coulomb specification can be obtained as follows:
Figure BDA0003367143340000021
in the formula, mu0In order to achieve a magnetic permeability in a vacuum,
Figure BDA0003367143340000022
k is the wave number of the electromagnetic field,
Figure BDA0003367143340000023
for admittance, ε is dielectric constant, σ is conductivity, JsIs the source item.
As a further improvement of the above technical solution, based on the quadratic field method, equation (1) is split into a background field equation set and a control equation set of a quadratic field vector bit and a scalar bit, where the control equations of the quadratic field vector bit and the scalar bit are:
Figure BDA0003367143340000031
in the formula, AsIs a secondary field vector bit, phisIn the form of a quadratic field scalar bit,
Figure BDA0003367143340000032
kbnumber of background electromagnetic field, JaIn order to scatter the current, the current is,
Figure BDA0003367143340000033
in order to have an abnormal admittance ratio,
Figure BDA0003367143340000034
in order to be the background admittance ratio,
Figure BDA00033671433400000311
as Hamiltonian, σbFor background conductivity, E is the total electric field;
performing horizontal two-dimensional Fourier transform on the formula (2), and converting the vector bit AsThe space-wavenumber domain secondary field control equation set obtained by expanding the vector bit in the x, y and z directions is as follows:
Figure BDA0003367143340000035
in the formula (I), the compound is shown in the specification,
Figure BDA0003367143340000036
is a space-wavenumber domain quadratic field vector bit,
Figure BDA0003367143340000037
is a space-wavenumber domain quadratic field scalar bit,
Figure BDA0003367143340000038
Figure BDA0003367143340000039
is the scattering current in the x, y and z directions, kx、kyThe wave numbers corresponding to the x and y directions in Fourier transform.
As a further improvement of the above technical solution, the two-dimensional fourier transform in step 3 and the two-dimensional inverse fourier transform in step 5 are independent of each other between the layered media, and are performed in parallel by Open MP.
As a further improvement of the above technical solution, in step 6, the iteration convergence condition is:
Figure BDA00033671433400000310
wherein E is an error, E(n)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) The nth iteration of the node of (a) to obtain the 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 in the spatial domain obtained in the (n-1) th iteration(0)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) The background electric field of the node of (a);
when e is<10-4If so, the iterative convergence condition is satisfied, otherwise, the iterative convergence condition is not satisfied.
The invention provides a three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method, which has the following beneficial technical effects:
1. by utilizing horizontal two-dimensional Fourier transform, a control equation met by a secondary vector bit scalar bit of a three-dimensional space domain is split into one-dimensional ordinary differential control equations under a plurality of wave numbers, the calculation amount and the storage requirement are greatly reduced, and the calculation efficiency of numerical simulation can be improved;
2. the one-dimensional control equation is fitted by adopting a quadratic interpolation function, and is solved by adopting a one-dimensional finite element method, so that the calculation speed is high;
3. the forward and reverse Fourier transformation process is parallel by adopting Open MP, so that the utilization of computer multithreading is improved, and the algorithm efficiency is further improved;
4. the three-dimensional magnetotelluric high-efficiency and high-precision numerical simulation is realized, and important technical support is provided for the precise inversion of large-scale magnetotelluric measured data.
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 finite element numerical simulation method of a three-dimensional magnetotelluric space-wavenumber domain according to an embodiment of the present 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 calculated magnetotelluric resistivity in an example of an embodiment of the present invention, where (a) is an apparent resistivity in xx mode, (b) is an apparent resistivity in xy mode, (c) is an apparent resistivity in yx mode, and (d) is an apparent resistivity in yy mode;
fig. 4 shows exemplary calculated values of magnetotelluric phases in an embodiment of the present invention, where (a) is a phase in xx mode, (b) is a phase in xy mode, (c) is a phase in yx mode, and (d) is a phase in yy mode.
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 numerical simulation method of finite elements of a three-dimensional magnetotelluric space-wavenumber domain. The control equations of the secondary vector bit and the scalar bit are subjected to two-dimensional Fourier transform in the horizontal direction, and the three-dimensional control equation is converted into one-dimensional control equations under a plurality of wave numbers, so that the calculation amount and the storage requirement are reduced, and the calculation efficiency is improved. The one-dimensional control equation is fitted by adopting a quadratic interpolation function, is solved by adopting a one-dimensional finite element method to obtain a 23-diagonal equation, is solved by adopting a catch-up method, and is high in calculation speed. The forward and reverse Fourier transformation process is parallel by adopting Open MP, so that the utilization of computer multithreading is improved, and the algorithm efficiency is further improved. The three-dimensional magnetotelluric high-efficiency and high-precision numerical simulation is realized, and important technical support is provided for the precise inversion of large-scale magnetotelluric measured data.
Specifically, the model of the background electromagnetic field is a uniform layered medium, and the conversion process of the Maxwell equation set is as follows:
introducing the vector bit and the scalar bit under the coulomb specification converts the Maxwell equation system about the electromagnetic field into a control equation system about the vector bit and the scalar bit based on the coulomb specification. The implementation process is as follows:
let the time constant be eiωtWherein i is an imaginary unit, ω is an angular frequency, ω is 2 pi f, and f is a calculation frequency, and a control equation system satisfied by the vector bit a and the scalar bit Φ based on the coulomb specification can be obtained as follows:
Figure BDA0003367143340000051
in the formula, mu0In order to achieve a magnetic permeability in a vacuum,
Figure BDA0003367143340000061
k is the wave number of the electromagnetic field,
Figure BDA0003367143340000062
for admittance, ε is dielectric constant, σ is conductivity, JsIs the source item.
Based on the quadratic field method, the control equation set in formula (1) is divided into a background field equation set and a control equation set with a quadratic field vector bit and a scalar bit, wherein the background medium is a uniform layered medium, the background electromagnetic field can be obtained by an analytic solution, and the implementation process is a conventional means, so that details are not repeated in this embodiment. The governing equations for the quadratic field vector bit and scalar bit are:
Figure BDA0003367143340000063
in the formula, AsIs a secondary field vector bit, phisIn the form of a quadratic field scalar bit,
Figure BDA0003367143340000064
kbnumber of background electromagnetic field, JaIn order to scatter the current, the current is,
Figure BDA0003367143340000065
in order to have an abnormal admittance ratio,
Figure BDA0003367143340000066
in order to be the background admittance ratio,
Figure BDA00033671433400000613
as Hamiltonian, σbFor background conductivity, E is the total electric field;
then, the control equation system of the secondary field vector position and the scalar position in the formula (2) is subjected to horizontal two-dimensional Fourier transform, and the vector position AsThe vector bit in the three directions of x, y and z is expanded to obtain a one-dimensional space domain secondary field control equation set under a two-dimensional wavenumber domain, namely a space-wavenumber domain secondary field control equation set, which is as follows:
Figure BDA0003367143340000067
in the formula (I), the compound is shown in the specification,
Figure BDA0003367143340000068
is a space-wavenumber domain quadratic field vector bit,
Figure BDA0003367143340000069
is a space-wavenumber domain quadratic field scalar bit,
Figure BDA00033671433400000610
Figure BDA00033671433400000611
is the scattering current in the x, y and z directions, kx、kyIs the wave number corresponding to the x and y directions in Fourier transform, whereinUsing standard Fourier transform sampling law, and sampling point number N in horizontal direction of space domainxAnd NyAnd when the sampling interval is delta x and delta y, calculating the wave number k corresponding to the x direction and the y direction according to the literature (Chenlongwei, Daishikun, Wumeiping, discrete frequency calculation when applying FFT algorithm of any sampling point number, geophysical progress 2016,31(01):164-169.)xAnd ky。
The relationship between the secondary electric field, the secondary magnetic field, and the space-wavenumber domain secondary field vector bit and scalar bit is:
Figure BDA00033671433400000612
in the formula, EsIs a secondary electric field, HsIs a secondary magnetic field.
Referring to fig. 1, based on the above equation transformation, the method for numerical simulation of finite elements in three-dimensional magnetotelluric space-wavenumber domain in this embodiment specifically includes the following steps:
step 1, establishing a cuboid target area in a rectangular coordinate system, and dispersing the target area in the x, y and z directions, wherein a three-dimensional conductivity abnormal body is in the target area. In the process of dispersing the target areas in the x direction, the y direction and the z direction, the x direction and the y direction are uniformly divided, the z direction is uniformly or non-uniformly divided, and the number of dividing nodes in the x direction, the y direction and the z direction is respectively set as Nx、Ny、Nz。
And 2, giving the conductivity of each sampling node as sigma, wherein the abnormal body conductivity is different from the surrounding conductivity. The electric conductivity without an abnormal body node is the background electric conductivity sigmabIt does not change in the horizontal direction but only in the vertical z direction, and can be found
Figure BDA0003367143340000075
And setting the polarization direction of the earth electromagnetic field source, obtaining a background electric field and a background magnetic field of a sampling node in the target region based on the background field equation set, the background conductivity and the frequency, and taking the background electric field as an initial electric field. Wherein, according to the frequencyf and background conductivity σbCalculating background electromagnetic fields of x-polarization and y-polarization directions respectively: eb(xi,yj,zk) Three components E ofbx(xi,yj,zk)、Eby(xi,yj,zk)、Ebz(xi,yj,zk)，Hb(xi,yj,zk) Three components of (H)bx(xi,yj,zk)、Hby(xi,yj,zk)、Hbz(xi,yj,zk) 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。
Step 3, multiplying the abnormal conductivity of the sampling node of the target area by the initial electric field to obtain the scattering current of the sampling node, and performing two-dimensional Fourier transform on the scattering current in the horizontal direction to obtain the two-dimensional wave number domain scattering current in the x, y and z directions
Figure BDA0003367143340000071
And 4, substituting the two-dimensional wavenumber domain scattering current into a space-wavenumber domain secondary field control equation set, obtaining a space-wavenumber domain secondary field vector bit and a scalar bit by adopting a one-dimensional finite element method based on secondary interpolation, and obtaining a space-wavenumber domain secondary electric field by utilizing the relation between the space-wavenumber domain secondary field vector bit and the scalar bit and the electric field. Wherein, the diagonal equation system formed by the one-dimensional finite element method based on the quadratic interpolation function is CAX ═ B, where CAIs a coefficient matrix which is a 23 diagonal matrix; x is an unknown quantity for solving space-wavenumber domain, comprising
Figure BDA0003367143340000072
And
Figure BDA0003367143340000073
b is the right end term. The number of such 23 diagonal equations sets is Nx×NyAnd (4) respectively. This embodiment converts a three-dimensional large-scale system of equations into Nx×NyThe one-dimensional finite element equation set has small calculated amount and less storage requirement, and greatly improves the calculation speed.
Step 5, performing two-dimensional Fourier inverse transformation on the space-wavenumber domain secondary electric field to obtain a space domain secondary electric field, and adding the space domain secondary electric field and a background electric field to obtain a space domain total electric field;
step 6, judging whether the total electric field of the spatial domain meets a given iterative convergence condition, specifically:
Figure BDA0003367143340000074
wherein E is an error, E(n)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) The nth iteration of the node of (a) to obtain the 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 in the spatial domain obtained in the (n-1) th iteration(0)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) The background electric field of the node of (a);
if e<10-4If the space-wavenumber domain secondary magnetic field and the vector bit scalar bit are in a linear relationship with each other, obtaining a space-wavenumber domain secondary magnetic field, 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;
if e<10-4If the iterative convergence condition is not satisfied, the total electric field of the spatial domain is taken as the initial electric field, and the step 3-6 is repeated.
And 7, obtaining an electromagnetic field obtained by forward modeling of the earth electromagnetic field source polarized in the x direction and the y direction based on the steps 2-6, and further obtaining apparent resistivity and phase at the measuring point.
It should be noted that the two-dimensional fourier transform in step 3 and the two-dimensional inverse fourier transform in step 5 are independent from each other between the layered media, and are performed in parallel by the Open MP, so as to improve the utilization of computer multithreading and further improve the computational efficiency.
In the following, a design model example is combined to perform magnetotelluric three-dimensional numerical simulation on the space-wavenumber domain finite element method provided by this embodiment, so as to verify the correctness and efficiency of the method. The computer tested was an Intel (R) core (TM) i7-6700HQ CPU with a master frequency of 2.60GHz and a 16GB 64-bit win10 system as memory.
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 σbThree-dimensional numerical simulation of the earth electromagnetic field was performed at a frequency f of 10Hz and 0.01S/m. The calculation range is-1000 m in the x direction, 1000-1000 m in the y direction, the calculation range is 0-1000 m in the z direction, the number of nodes of the subdivision grid is 101 multiplied by 101, the three directions are uniformly subdivided, both the delta x, the delta y and the delta z are 10m, the abnormal body range is-100 m in the x direction, 200-200 m in the y direction, 200-400 m in the z direction, and the conductivity sigma of the abnormal body is 0.1S/m. 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 graph comparing the magnetotelluric resistivity calculated by the method of this embodiment with the software of INTEM3D in the z-0 plane, and fig. 4 is a graph comparing the magnetotelluric resistivity calculated by the method of this embodiment with the software of INTEM3D in the z-0 plane. Rhoxx、ρxy、ρyx、ρyy、
Figure BDA0003367143340000081
The average relative errors in the plane of the eight components are 1.02%, 0.45%, 0.43%, 0.99%, 1.21%, 0.012%, 0.008%, 1.05%, ρxy、ρyx、
Figure BDA0003367143340000082
The relative errors are all less than 1 percent,the relative error of other components is about 1%, and the correctness of the method of the embodiment is verified. In this embodiment, the number of the computing nodes is 101 × 101 × 101-1030301, when the space-wavenumber domain finite element numerical simulation iteration is terminated, the X polarization direction is iterated 13 times, the Y polarization direction is iterated 19 times, Open MP parallel is adopted for forward and reverse fourier transform, when the parallel threads are 8 threads, the total computing time is 60.12s, and the occupied memory is 2032.5 MB.
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 (7)

1. A three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method is characterized by comprising the following steps of:
step 1, establishing a cuboid target area in a rectangular coordinate system, and dispersing the target area in the directions of x, y and z, wherein a three-dimensional conductivity abnormal body is in the target area;
step 2, setting the polarization direction of a geodetic electromagnetic field source, obtaining a background electric field and a background magnetic field of a sampling node in a target area based on a background field equation set, background conductivity and frequency, and taking the background electric field as an initial electric field;
step 3, multiplying the abnormal conductivity of the sampling node of the target area by the initial electric field to obtain the scattering current of the sampling node, and performing two-dimensional Fourier transform on the scattering current in the horizontal direction to obtain two-dimensional wavenumber domain scattering current;
step 4, substituting the two-dimensional wavenumber domain scattering current into a space-wavenumber domain secondary field control equation set, obtaining a space-wavenumber domain secondary field vector bit and a scalar bit by adopting a one-dimensional finite element method based on secondary interpolation, and obtaining a space-wavenumber domain secondary electric field by utilizing the relation between the space-wavenumber domain secondary field vector bit and the scalar bit and the electric field;
step 5, performing two-dimensional Fourier inverse transformation on the space-wavenumber domain secondary electric field to obtain a space domain secondary electric field, and adding the space domain secondary electric field and a background electric field to obtain a space domain total electric field;
step 6, judging whether the total electric field of the spatial domain meets a given iteration convergence condition:
if the magnetic field is satisfied, obtaining a space-wavenumber domain secondary magnetic field based on the relation between the space-wavenumber domain secondary field and the vector bit scalar bit and the magnetic field, obtaining the space domain secondary magnetic field after the two-dimensional inverse Fourier transform in the horizontal direction, adding a background magnetic field to obtain a space domain total magnetic field, and then outputting the space domain total electric field and the space domain total magnetic field;
if not, taking the total electric field of the spatial domain as an initial electric field, and repeating the step 3-6;
and 7, obtaining an electromagnetic field obtained by forward modeling of the earth electromagnetic field source polarized in the x direction and the y direction based on the steps 2-6, and further obtaining apparent resistivity and phase at the measuring point.
2. The method according to claim 1, wherein in step 1, in the process of discretizing the target region in x, y, and z directions, the x and y directions are uniformly subdivided, the z direction is uniformly or non-uniformly subdivided, and the numbers of nodes of subdivision in x, y, and z directions are respectively set to be Nx、Ny、Nz。
3. The method for numerical simulation of finite element in three-dimensional magnetotelluric space-wavenumber domain as defined in claim 1, wherein the background field equations in step 2 and the space-wavenumber domain secondary field control equations in step 4 are obtained by:
introducing vector bits and scalar bits under the coulomb specification, and converting a Maxwell equation system related to the electromagnetic field into a control equation system related to the vector bits and the scalar bits based on the coulomb specification;
based on a quadratic field method, dividing a control equation set into a background field equation set and a control equation set with a quadratic field vector bit and a scalar bit;
and performing two-dimensional Fourier transform on the control equation set of the secondary field vector bit and the scalar bit in the horizontal direction to obtain a one-dimensional space domain secondary field control equation set under a two-dimensional wavenumber domain, namely a space-wavenumber domain secondary field control equation set.
4. The method of claim 3, wherein the time constant is eiωtWherein i is an imaginary unit, ω is an angular frequency, ω is 2 pi f, and f is a calculation frequency, and a control equation system satisfied by the vector bit a and the scalar bit Φ based on the coulomb specification can be obtained as follows:
Figure FDA0003367143330000021
in the formula, mu0In order to achieve a magnetic permeability in a vacuum,
Figure FDA0003367143330000022
k is the wave number of the electromagnetic field,
Figure FDA0003367143330000023
for admittance, ε is dielectric constant, σ is conductivity, JsIs the source item.
5. The three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method of claim 4, wherein equation (1) is split into a background field equation set and a control equation set of a secondary field vector bit and a scalar bit based on a secondary field method, the control equations of the secondary field vector bit and the scalar bit being:
Figure FDA0003367143330000024
in the formula, AsIs a secondary field vector bit, phisIn the form of a quadratic field scalar bit,
Figure FDA0003367143330000025
kbnumber of background electromagnetic field, JaIn order to scatter the current, the current is,
Figure FDA0003367143330000026
Figure FDA0003367143330000027
in order to have an abnormal admittance ratio,
Figure FDA0003367143330000028
for background admittance ratio,. v is the Hamiltonian, σbFor background conductivity, E is the total electric field;
performing horizontal two-dimensional Fourier transform on the formula (2), and converting the vector bit AsExpanded into vector bits in the x, y and z directions, and marked as Axs、Ays、AzsThe obtained space-wavenumber domain secondary field control equation set is as follows:
Figure FDA0003367143330000029
in the formula (I), the compound is shown in the specification,
Figure FDA00033671433300000210
respectively, space-wavenumber domain secondary field vector bits,
Figure FDA00033671433300000211
is a space-wavenumber domain quadratic field scalar bit,
Figure FDA00033671433300000212
is the scattering current in the x, y and z directions, kx、kyThe wave numbers corresponding to the x and y directions in Fourier transform.
6. A three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method according to any one of claims 1 to 5, wherein the two-dimensional Fourier transform in step 3 and the two-dimensional inverse Fourier transform in step 5 are independent of each other between the respective layered media and are performed in parallel by Open MP.
7. The three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method according to any one of claims 1 to 5, wherein in step 6, the iterative convergence condition is:
Figure FDA0003367143330000031
Figure FDA0003367143330000032
wherein E is an error, E(n)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) The nth iteration of the node of (a) to obtain the 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 in the spatial domain obtained in the (n-1) th iteration(0)(xi,yj,zk) Is a coordinate of (x)i,yj,zk) The background electric field of the node of (a);
when e is<10-4If so, the iterative convergence condition is satisfied, otherwise, the iterative convergence condition is not satisfied.
CN202111386046.8A 2021-11-22 2021-11-22 Three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method Active CN114065586B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111386046.8A CN114065586B (en) 2021-11-22 2021-11-22 Three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111386046.8A CN114065586B (en) 2021-11-22 2021-11-22 Three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method

Publications (2)

Publication Number Publication Date
CN114065586A true CN114065586A (en) 2022-02-18
CN114065586B CN114065586B (en) 2022-09-02

Family

ID=80278890

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111386046.8A Active CN114065586B (en) 2021-11-22 2021-11-22 Three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method

Country Status (1)

Country Link
CN (1) CN114065586B (en)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115017782A (en) * 2022-08-08 2022-09-06 中南大学 Three-dimensional natural source electromagnetic field calculation method considering medium anisotropy
CN115292973A (en) * 2022-10-09 2022-11-04 中南大学 Arbitrarily sampled space wave number domain three-dimensional magnetic field numerical simulation method and system
CN115795231A (en) * 2022-10-09 2023-03-14 中南大学 Space wave number mixed domain three-dimensional high-intensity magnetic field iteration method and system
CN116244877A (en) * 2022-09-05 2023-06-09 中南大学 Three-dimensional magnetic field numerical simulation method and system based on 3D AS-FT

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090006053A1 (en) * 2006-03-08 2009-01-01 Carazzone James J Efficient Computation Method for Electromagnetic Modeling
US20100298967A1 (en) * 2009-05-19 2010-11-25 Frisken Sarah F Method for Reconstructing a Distance Field of a Swept Volume at a Sample Point
CN108509693A (en) * 2018-03-13 2018-09-07 中南大学 Three-dimensional frequency domain controllable source method for numerical simulation
WO2021068527A1 (en) * 2019-10-12 2021-04-15 中南大学 Numerical simulation method for footprint-guided high-efficiency airborne electromagnetic survey
CN113656750A (en) * 2021-10-20 2021-11-16 中南大学 Magnetic induction intensity calculation method of strong magnetic medium based on space wave number domain

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090006053A1 (en) * 2006-03-08 2009-01-01 Carazzone James J Efficient Computation Method for Electromagnetic Modeling
US20100298967A1 (en) * 2009-05-19 2010-11-25 Frisken Sarah F Method for Reconstructing a Distance Field of a Swept Volume at a Sample Point
CN108509693A (en) * 2018-03-13 2018-09-07 中南大学 Three-dimensional frequency domain controllable source method for numerical simulation
WO2021068527A1 (en) * 2019-10-12 2021-04-15 中南大学 Numerical simulation method for footprint-guided high-efficiency airborne electromagnetic survey
CN113656750A (en) * 2021-10-20 2021-11-16 中南大学 Magnetic induction intensity calculation method of strong magnetic medium based on space wave number domain

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
SHIKUN DAI 等: "Three-dimensional numerical modeling of gravity and magnetic anomaly in a mixed space-wavenumber domain", 《GEOPHYSICS》 *
朱峰 等: "一种求解复波数域波动弥散方程的通用数值方法", 《应用力学学报》 *
苏晓波 等: "大地电磁三维矢量有限元正演研究", 《地球物理学进展》 *
陈龙 等: "应用任意采样点数FFT算法时离散频率计算", 《地球物理学进展》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
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
CN115292973A (en) * 2022-10-09 2022-11-04 中南大学 Arbitrarily sampled space wave number domain three-dimensional magnetic field numerical simulation method and system
CN115795231A (en) * 2022-10-09 2023-03-14 中南大学 Space wave number mixed domain three-dimensional high-intensity magnetic field iteration method and system
CN115795231B (en) * 2022-10-09 2023-08-04 中南大学 Space wave number mixed domain three-dimensional strong magnetic field iteration method and system

Also Published As

Publication number Publication date
CN114065586B (en) 2022-09-02

Similar Documents

Publication Publication Date Title
CN114065586B (en) Three-dimensional magnetotelluric space-wavenumber domain finite element numerical simulation method
Jin et al. A meshless method for some inverse problems associated with the Helmholtz equation
CN106980736B (en) A kind of ocean controllable source electromagnetic method finite element forward modeling method of anisotropic medium
Li et al. Hermite–Cloud: a novel true meshless method
Xiao et al. A local Heaviside weighted meshless method for two-dimensional solids using radial basis functions
CN114065585B (en) Three-dimensional electrical source numerical simulation method based on coulomb specification
Liu et al. A meshfree Hermite‐type radial point interpolation method for Kirchhoff plate problems
CN113051779B (en) Numerical simulation method of three-dimensional direct-current resistivity method
CN113792445B (en) Three-dimensional magnetotelluric numerical simulation method based on integral equation method
CN111079278B (en) Processing method for three-dimensional time domain hybridization discontinuous Galerkin method with additional electromagnetic source item
Sun et al. An Improved Interpolating Element‐Free Galerkin Method Based on Nonsingular Weight Functions
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
Ren et al. Static responses of magneto-electro-elastic structures in moisture field using stabilized node-based smoothed radial point interpolation method
Xiao Local Heaviside weighted MLPG meshless method for two-dimensional solids using compactly supported radial basis functions
CN115017782B (en) Three-dimensional natural source electromagnetic field calculation method considering medium anisotropy
Malhotra et al. Efficient Convergent Boundary Integral Methods for Slender Bodies
CN115659579B (en) Three-dimensional gravity field numerical simulation method and system based on 3D AS-FT
CN114970289B (en) Three-dimensional magnetotelluric anisotropy forward modeling numerical simulation method, equipment and medium
Yang et al. Two‐Dimensional Multiple‐Snapshot Grid‐Free Compressive Beamforming Using Alternating Direction Method of Multipliers
IT202100015602A1 (en) PROCEDURE FOR CALCULATING PHYSICAL QUANTITIES OF A CONDUCTIVE BODY, CORRESPONDING PROCESSING SYSTEM AND COMPUTER PRODUCT
Lv et al. A Kriging interpolation-based boundary face method for 3D potential problems
Zhou et al. An efficient hybrid direct-iterative solver for three-dimensional higher-order edge-based finite element simulation for magnetotelluric data in anisotropic media
Baxansky et al. Calculating geometric properties of three-dimensional objects from the spherical harmonic representation
Chen et al. A coupled RBF method for the solution of elastostatic problems

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
GR01 Patent grant
GR01 Patent grant