阅读说明：本技术 一种基于超收敛插值逼近的三维电磁慢扩散数值模拟方法 (Three-dimensional electromagnetic slow diffusion numerical simulation method based on hyperconvergence interpolation approximation ) 是由 嵇艳鞠 王世鹏 赵雪娇 吴琼 黎东升 关珊珊 于 2019-11-12 设计创作，主要内容包括：本发明涉及一种基于超收敛插值逼近的三维电磁慢扩散数值模拟方法,通过将复电导率模型引入频域Maxwell方程组后,电磁场扩散方程中含有复频变量的负分数次幂项,先进行频-时变换得到含有Caputo分数阶微分项的时间域控制方程；再采用Alikhanov超收敛插值逼近方法,对电场控制方程中Caputo分数阶导数进行超收敛逼近,获得分数阶微分项的非均匀步长离散近似表达式,从而完成时间域分数阶微分项的稳定、高精度直接求解；最后基于有限差分算法对控制方程进行离散,推导出电场和磁场迭代方程,最终实现了三维时域电磁慢扩散的高精度数值模拟。本发明目的在于克服分数阶微分求解的弱奇异不稳定性及误差较大问题,实现三维时域电磁慢扩散的高精度数值模拟。(The invention relates to a three-dimensional electromagnetic slow diffusion numerical simulation method based on hyperconvergence interpolation approximation.A complex conductivity model is introduced into a frequency domain Maxwell equation set, an electromagnetic field diffusion equation contains a negative fraction power term of a complex frequency variable, and frequency-time transformation is carried out to obtain a time domain control equation containing a Caputo fractional order differential term; then, performing hyperconvergence approximation on the Caputo fractional derivative in the electric field control equation by adopting an Alikhanov hyperconvergence interpolation approximation method to obtain a non-uniform step discrete approximate expression of a fractional derivative term, thereby completing the stable and high-precision direct solution of the time domain fractional derivative term; and finally, dispersing the control equation based on a finite difference algorithm, deducing an electric field and magnetic field iterative equation, and finally realizing high-precision numerical simulation of three-dimensional time domain electromagnetic slow diffusion. The invention aims to solve the problems of weak singular instability and large error of fractional order differential solution and realize high-precision numerical simulation of three-dimensional time domain electromagnetic slow diffusion.)

1. A three-dimensional electromagnetic slow diffusion numerical simulation method based on hyperconvergence interpolation approximation is characterized by comprising the following steps:

1) introducing a complex conductivity expression into a Maxwell equation set as an electromagnetic field iterative control equation, wherein the early stage of the electromagnetic diffusion process after introducing the complex conductivity is consistent with the classical diffusion, the late stage is slowly attenuated, and the electromagnetic trailing phenomenon in actual observation is met;

2) converting the negative fractional power of the complex frequency variable in the control equation into a positive fractional power, and performing frequency-time conversion to obtain a time domain electromagnetic field control equation containing a Caputo fractional order differential term;

3) performing quadratic interpolation polynomial approximation on electric field components in an electric field control equation, then deriving the two ends of the equation, substituting the derivation into a Caputo fractional derivative hyperconvergence interpolation approximation formula to obtain a fractional differential term non-uniform step discrete approximation formula, and completing direct solution of a fractional differential equation;

4) dispersing each partial derivative item of the control equation by adopting a finite difference method, and recursively obtaining the control equation of each component of the electromagnetic field;

5) subdividing the calculation area by adopting a non-uniform three-dimensional Yee's grid, setting the conductivity, the permeability and the artificial dielectric constant of the calculation area, and calculating an initial field;

6) loading a C-PML boundary condition, accelerating by using a GPU, and carrying out iterative operation on each component of the electric field and the magnetic field in observation time;

7) after the iterative computation is finished, extracting each component response of the electromagnetic field, mapping, and analyzing and processing the result;

in the step 1), the complex conductivity is mainly used for simulating the electromagnetic field slow diffusion phenomenon, the diffusion process of a complex medium usually spans multiple spatial scales, and the heterogeneity and anisotropy of the medium determine that the electromagnetic diffusion needs to be subjected to multi-scale composite simulation; introducing a weight coefficient which is more than 0 and less than 1, describing the change of the rock along with the frequency conductivity by adjusting the weight coefficient, and degenerating the electromagnetic diffusion process into a classical electromagnetic diffusion process when the weight coefficient is equal to zero;

the complex conductivity expression is:

σ(ω)＝σ₀+ησ₀(iω)^-β(1)

σ in formula (1)₀The direct current conductivity is adopted, η is a weight coefficient, i is an imaginary number unit, omega is an angular frequency, β is roughness, and the values of β and η are both larger than 0 and smaller than 1;

in step 2), the negative fractional power of the complex frequency variable in the control equation is converted into a positive fractional power, and then frequency-time conversion is carried out to obtain:

where ε is the dielectric constant in vacuum, μ is the permeability, E is the electric field, H_mIs a magnetic field and is used for generating a magnetic field,

wherein the process of performing super-convergence interpolation approximation solving on the electric field component in the electric field control equation in the step 3) is as follows:

firstly, in the interval [ t_k-1,t_k]Upper, using three points (t)_k-1,E(t_k-1)),(t_k,E(t_k)),(t_k+1,E(t_k+1) Approximation of a quadratic interpolation polynomial to e (t) can be obtained:

and then carrying out derivation on the formula (4):

wherein E (t)_k-1)、E(t_k)、E(t_k+1) The value of the electric field, Δ t, at different moments of the electric field_k＝t_k-t_k-1，Δt_k+1＝t_k+1-t_k；

The Caputo fractional order hyperconvergence format is:

in the formula (6)

Finally, the formula (5) is substituted into the formula (6) to obtain an ultra-convergence discrete approximation format of fractional order differential under any step length:

in the step 4), a finite difference method is adopted to disperse the partial derivative term in the formula (3), an electromagnetic field control equation is expanded, difference discrete approximation is carried out, and iteration equations of all components of the electromagnetic field are obtained, wherein the iteration equation of the electromagnetic field in the x direction is as follows:

Technical Field

The invention relates to a three-dimensional electromagnetic slow diffusion numerical simulation method based on hyperconvergence interpolation approximation, which is suitable for time domain electromagnetic slow diffusion numerical simulation, in particular to numerical simulation of magnetic source electromagnetic slow diffusion.

Background

The magnetic source transient electromagnetic detection method selects a coil as an emission source, loads step current to excite an underground good conductor to generate eddy current, and acquires abundant underground electrical information by a receiver through acquiring a secondary induction field. The gradual maturity and the practical application of Superconducting Quantum interferometers (SQUIDs, Superconducting Quantum Interference devices) and atomic magnetometers obviously improve the exploration precision, obtain more accurate electromagnetic attenuation signals in actual geology, gradually observe the transient electromagnetic slow diffusion phenomenon which does not conform to the theoretical-5/2 power law, and along with the continuous research on the characteristics of geological structures, the traditional single-scale and single-parameter geological conductivity model can not meet the requirements of high-precision geological detection. Establishing a physical model more conforming to the geological structure is one of the key technologies for improving the electromagnetic detection resolution.

The Three-dimensional Finite-difference Time-domain method is one of the important methods for electromagnetic field numerical calculation, and the basic idea is that the first-order partial derivatives of field quantity to Time and space are approximated by central difference, and the propagation process of wave is simulated by recursion in Time domain, so as to obtain field distribution.

Chester J.Weiss et al observed a time domain electromagnetic slow diffusion response with TEM47 in Brassuzshire, USA, and excluded the possibility that the slow diffusion is due to a layered structure or lateral anisotropy, giving the slow diffusion response a fractional order diffusion caused by a rough medium. Mark e.everett et al propose that the medium has a non-uniformity, i.e. roughness, and introduce a roughness parameter into the expression for the frequency domain conductivity. Detwiler et al studied the problem of solute migration in the fracture and found that solute migration in the fracture medium is a two-scale coupled diffusion of Taylor diffusion inside the fracture and normal diffusion at the fracture surface; in geological problems, the propagation of electromagnetic waves in large-scale rough medium fracture media also has the same double-coupling property, namely double-scale coupling diffusion of common diffusion and secondary diffusion. The structure of the random change of the conductivity along with the space is mainly researched in China, and the slow diffusion of the change of the conductivity along with the time is rarely researched.

Chinese patent CN107766666A discloses a three-dimensional time domain electromagnetic anomalous diffusion simulation method based on a fractional order differential method, which obtains a time domain fractional order differential-integral expression of an electromagnetic field through frequency-time conversion; performing time domain dispersion on integral and differential terms of a diffusion equation by adopting a Riemann-Liouville fractional order integral and finite difference method to construct a time domain iterative formula of an electric field and a magnetic field; and loading initial conditions and boundary conditions to realize numerical simulation of three-dimensional time domain electromagnetic anomalous diffusion. The numerical simulation of the three-dimensional time domain electromagnetic anomalous diffusion is realized.

Chinese patent CN 106776478A discloses a discrete fractional order differential method based on step-by-step calculation in anomalous diffusion, which takes the ratio of two gamma functions in the discrete format of a time-space anomalous diffusion equation as a parameter, utilizes the recurrence relation of the gamma functions, combines the thought of distribution calculation in numerical simulation and the thought of integral consideration, reduces the calculation limit of the gamma functions, expands the number of points of simulation, and improves the simulation efficiency.

Chinese patent CN108897052A discloses a three-dimensional time domain electromagnetic slow diffusion simulation method based on fractional order linear approximation, which introduces generalized conductivity into Maxwell equation set, carries out linear approximation on fractional order differential terms, and carries out frequency-time conversion to obtain integer order differential of approximate time domain, thus completing conversion from fractional order to integer order and achieving the purpose of simulating three-dimensional time domain electromagnetic slow diffusion phenomenon.

The method discloses a research method for abnormal diffusion at home and abroad. However, for the slow diffusion phenomenon in electromagnetic anomalous diffusion, at present, few researches on fractional order three-dimensional finite difference operation in a time domain are available, and for the problems that a fractional order operator has weak singularity in direct operation in the time domain, results are unstable, errors are large and the like, how to accurately and stably perform fractional order difference operation in the electromagnetic slow diffusion field is a technical problem which is urgently solved by technical personnel in the field.

Disclosure of Invention

The invention aims to solve the technical problem of providing a three-dimensional electromagnetic slow diffusion numerical simulation method based on hyperconvergence interpolation approximation. Performing hyperconvergence approximation on the Caputo fractional derivative in the electric field control equation by adopting an Alikhanov hyperconvergence interpolation approximation method to obtain a non-uniform step discrete approximation formula of a fractional derivative term, and completing stable and high-precision direct solution of the time domain fractional derivative term; and finally, discretizing the space and time partial derivatives of the control equation based on a finite difference algorithm, and deducing an iterative equation of each component of the electric field and the magnetic field. Finally, stable and high-precision numerical simulation of three-dimensional time domain electromagnetic slow diffusion is realized. The invention aims to solve the problems of weak singular instability and large error of fractional order differential solution and realize high-precision numerical simulation of three-dimensional time domain electromagnetic slow diffusion.

The invention is realized in such a way that a three-dimensional electromagnetic slow diffusion numerical simulation method based on hyperconvergence interpolation approximation:

1) introducing a complex conductivity expression into a Maxwell equation set as an electromagnetic field iterative control equation, wherein the early stage of the electromagnetic diffusion process after introducing the complex conductivity is consistent with the classical diffusion, the late stage is slowly attenuated, and the electromagnetic trailing phenomenon in actual observation is met;

2) converting the negative fractional power of the complex frequency variable in the control equation into a positive fractional power, and performing frequency-time conversion to obtain a time domain electromagnetic field control equation containing a Caputo fractional order differential term;

3) performing quadratic interpolation polynomial approximation on electric field components in an electric field control equation, then deriving the two ends of the equation, substituting the derivation into a Caputo fractional derivative hyperconvergence interpolation approximation formula to obtain a non-uniform step length discrete approximation formula of a fractional differential term, and completing the direct solution of the fractional differential equation;

4) dispersing each partial derivative item of the control equation by adopting a finite difference method, and recursively obtaining the control equation of each component of the electromagnetic field;

5) subdividing the calculation area by adopting a non-uniform three-dimensional Yee's grid, setting the conductivity, the permeability and the artificial dielectric constant of the calculation area, and calculating an initial field;

6) loading a C-PML boundary condition, accelerating by using a GPU, and carrying out iterative operation on each component of the electric field and the magnetic field in observation time;

7) after the iterative computation is finished, extracting each component response of the electromagnetic field, mapping, and analyzing and processing the result;

in the step 1), the complex conductivity expression is as follows:

σ(ω)＝σ₀+ησ₀(iω)^-β(1)

σ in formula (1)₀The conductivity is defined to simulate the electromagnetic field slow diffusion phenomenon, and when the weight coefficient is equal to zero, the electromagnetic diffusion process degenerates to the classical electromagnetic diffusion process;

the control equation of the electromagnetic field after introducing the complex conductivity expression can be expressed as:

wherein epsilon is the dielectric constant in vacuum, mu is the magnetic conductivity, E is the electric field, H is the magnetic field,

is Hamiltonian;

further, the method for converting the negative fractional power of the complex frequency variable in the control equation into the positive fractional power in step 2) is to perform differential processing on both ends of the formula (2) at the same time, that is, multiplying by i ω at the same time, so as to obtain:

wherein, the time domain control equation obtained by performing frequency-time conversion on the frequency domain control equation in the step 2) is as follows:

wherein H_mIs a magnetic field and is used for generating a magnetic field,

for the Caputo fractional derivative, α is a non-negative real number, 0 is the lower integration limit, t is the upper integration limit, and the Caputo fractional derivative is expressed as:

equation (7) is called the left α order Caputo fractional derivative, a is the lower integration limit, x is the upper integration limit, f (x) is a function, n is an integer, n-1 ≦ α < n, and Γ is a Gamma function expressed as:

wherein, the process of carrying out quadratic interpolation polynomial approximation and derivation on the electric field component in the electric field control equation in the step 3) comprises the following steps:

in the interval [ t_k-1,t_k]Upper, using three points (t)_k-1,E(t_k-1)),(t_k,E(t_k)),(t_k+1,E(t_k+1) Approximation of a quadratic interpolation polynomial to e (t) can be obtained:

the derivative of equation (9) is:

wherein E (t)_k-1)、E(t_k)、E(t_k+1) The value of the electric field, Δ t, at different moments of the electric field_k＝t_k-t_k-1，Δt_k+1＝t_k+1-t_k；

Wherein, the process of carrying out super-convergence approximate direct solution on the fractional order differential term in the control equation in the step 3) comprises the following steps:

the Caputo fractional order hyperconvergence format is:

in the formula (11)

Substituting the formula (10) into the formula (11) to obtain an ultra-convergence discrete approximation format of fractional order differential under any step length:

wherein, in the step 4), the finite difference method is adopted to disperse the partial derivative term in the formula (6) to obtain:

where Δ t_n+1＝t_n+1-t_n，Δt_n＝t_n-t_n-1；

Expanding an electromagnetic field control equation, and performing difference discrete approximation to obtain an iterative equation of each component of the electromagnetic field, wherein the iterative equation of the electromagnetic field in the x direction is as follows:

compared with the prior art, the invention has the beneficial effects that: a complex conductivity model is introduced into a time domain electromagnetic three-dimensional numerical simulation control equation, so that the electromagnetic slow diffusion phenomenon can be effectively simulated; aiming at the problems that the fractional order differential item has weak singular instability and large approximate solution integral error when being directly solved in a time domain, the fractional order differential item is directly solved by adopting a super-convergence interpolation approximation method, so that the problem of oscillation in the solving process is solved, and the purpose of high-precision numerical simulation of the three-dimensional time domain electromagnetic slow diffusion phenomenon is achieved.

Drawings

FIG. 1 is a schematic diagram of a three-dimensional electromagnetic slow diffusion numerical simulation method based on hyperconvergence interpolation approximation;

FIG. 2 is a comparison graph and an error curve of induced electromotive force attenuation curves of a receiving coil and numerical integration solutions obtained by three-dimensional time domain electromagnetic slow diffusion simulation based on fractional order differential term hyperconvergence direct solution;

FIG. 3 is a graph comparing an induced electromotive force attenuation curve of a receiving coil obtained by a three-dimensional time domain electromagnetic slow diffusion simulation based on a complex conductivity model with induced electromotive force attenuation curves of a classical model and a generalized conductivity model;

FIG. 4 is a ground two-dimensional electromagnetic response plan obtained by three-dimensional time domain electromagnetic slow diffusion numerical simulation based on a complex conductivity model;

FIG. 5 is a aerial two-dimensional electromagnetic response plan view obtained by three-dimensional time domain electromagnetic slow diffusion numerical simulation based on a complex conductivity model;

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is further described in detail with reference to the following embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.