CN118551630B

CN118551630B - Method for accelerating finite volume of grid sequence of frequency domain electromagnetic field of structural grid

Info

Publication number: CN118551630B
Application number: CN202411026342.0A
Authority: CN
Inventors: 许勇; 陈坚强; 牟斌; 江涛; 丁明松; 高铁锁; 董维中; 梅杰; 蒋军; 张凡
Original assignee: Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Current assignee: Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority date: 2024-07-30
Filing date: 2024-07-30
Publication date: 2024-10-18
Anticipated expiration: 2044-07-30
Also published as: CN118551630A

Abstract

The invention discloses a method for accelerating finite volume of a grid sequence of a frequency domain electromagnetic field of a structural grid, which belongs to the field of calculating the numerical value of an electromagnetic frequency domain and comprises the following steps: simulation modeling; dividing a two-dimensional quadrilateral or three-dimensional hexahedral structure grid of a numerical calculation region, encrypting at a wall surface and a geometric singular position, and gradually separating the grid from a scattering wall surface to gradually sparsity; outputting a grid data file and a boundary condition file; inputting target calculation electromagnetic parameters and numerical calculation control parameters; inputting grid data and boundary condition information files, and initializing and calculating a space electromagnetic field; based on virtual time step iteration propulsion, each iteration carries out electromagnetic flux residual implicit calculation in each level grid according to a specific loop mode of a grid sequence, adjacent grid layer flux residual, conservation variable and electromagnetic parameter interpolation transmission are carried out, and acceleration iteration solution is carried out on a frequency domain electromagnetic field; outputting the data. The invention can solve the problem of large-scale electromagnetic scattering of any complicated shape and high-frequency electric large-size targets.

Description

Method for accelerating finite volume of grid sequence of frequency domain electromagnetic field of structural grid

Technical Field

The invention relates to the field of calculation of electromagnetic frequency domain values, in particular to a method for accelerating finite volume of a grid sequence of a structural grid frequency domain electromagnetic field.

Background

Electromagnetic scattering of a target with a complex appearance and electromagnetic interference of a complex electromagnetic environment all need to calculate electromagnetic field spatial distribution, the electromagnetic field meets Maxwell's equations, and the equations can be solved directly along with development of computer technology. The hyperbolic mathematical features, which are the same as the euler equation, facilitate the application of computational fluid dynamics (Computational Fluid Dynamics, CFD) techniques in electromagnetic field computation, with time-Domain finite difference methods (FINITE DIFFERENCE TIME Domain, FDTD) and time-Domain finite volume methods (Finite Volume Time Domain, FVTD) being the most notable. K.S.Yee published a precursor time domain finite difference algorithm in the 60 th century of 20 th, directly and differentially calculated a time-varying Max Wei Weifen equation set, successfully simulated the time domain response of the electromagnetic pulse and the action of an ideal conductor, and opened a new electromagnetic field time domain calculation method. In the Yee algorithm, first, a cartesian orthogonal grid is generated in a region of interest (a target and a certain space around the target), and electric field and magnetic field components are placed in a crossing manner at a value point of a grid space, so that the periphery of each electric field component is surrounded by the magnetic field component on each coordinate plane, and the periphery of each magnetic field component is surrounded by the electric field component, so that the electromagnetic field configuration meets the requirements of faraday induction law and ampere loop law, and the grid is commonly called a Yee grid. The FDTD method is based on Ye grids, has high calculation efficiency and parallelism expandability, is very simple, convenient and effective on the structured grids of the regular regions, and has the defects that the conservation characteristic of a discrete equation is difficult to ensure, the adaptability to irregular regions is poor, the processing capacity to complex models is weak, and the calculation accuracy is poor. The time domain finite volume method (FVTD) directly applies the integral form of maxwell's equations conservation law to discrete grid cells, by integrating the divergence equation over the control volume, converting the volume integral into a control surface normal flux surface integral, i.e., using the conservation law integral form in physical space, allowing the discontinuous function to be calculated. The finite volume method keeps the diversity of the finite difference method in the format construction, and can conveniently utilize the design thought and theoretical results of almost all finite difference methods; and the limited volume method has no limitation on mesh subdivision and unit shape, is easy to process complex appearance, and is a common method for most commercial and engineering CFD software at present.

Time domain calculations are useful for simulating wideband pulsed electromagnetic wave signal radiation, scattering, but if the incident wave is a single frequency simple harmonic signal, the electromagnetic field can be calculated in the frequency domain. The traditional frequency domain method mainly comprises an analysis method, a high-frequency approximation method and a full-wave numerical method. The analysis method can only solve electromagnetic scattering of the simple target with a special geometric shape, and cannot be used for the actual complex geometric target. The high-frequency approximation method comprises geometrical Optics (Geometrical Optics, GO), geometrical diffraction theory (Geometrical Theory of Diffraction, GTD), physical Optics (PO), physical diffraction theory (Physical Theory of Diffraction, PTD), consistency geometrical diffraction theory (Uniform Theory of Diffraction, UTD), consistency progressive theory (Uniform Asymptotic Theory, UAT), equivalent fringe current method (Method of Equivalent Current, MEC) and the like, only the scattering field generated by a part or a tiny unit under an incident wave is considered based on the high-frequency field locality principle, the mutual coupling between the parts or the units is not considered, and the high-frequency method has poor precision in analyzing electromagnetic scattering of a complex structural target. The full-wave numerical method directly solves the Maxwell partial differential equation or the electromagnetic flow integral equation, does not perform any approximation, has higher calculation precision, and can solve the electromagnetic problem of any frequency under the condition of computer resource permission. The high-precision full-wave electromagnetic numerical methods are mainly divided into two types: one is to solve the integral equation with current as a variable, including a Moment Method (MOM) and a subsequently developed multipole Method (Fast Multipole Method, FMM), a multilayer fast multipole Method (Multi-LEVEL FAST Multipole Algorithm, MLFMA); another class of FDTD methods and finite element methods (FINITE ELEMENT methods, FEM) solve Maxwell differential or Helmholtz wave equations with electromagnetic fields as variables.

These full-wave simulation methods generally involve an iterative solution process, so that in addition to the calculation accuracy, the calculation efficiency of numerical simulation is also a key ring related to the applicability of the method. The iteration maximum time step is still limited by the local stability condition, and the slow convergence speed can restrict the application of numerical simulation in practical engineering to a great extent. As engineering problems become more and more complex, the mesh size of calculation becomes larger and larger, how to improve the electromagnetic field calculation efficiency and accelerate the iterative convergence speed of the electromagnetic field becomes important.

Classical iteration methods such as Gauss-Seidel iteration, jaccobi iteration, SOR iteration and the like start to iterate the error attenuation between the approximate solution and the true solution for several times very fast, but the later iteration error attenuation is very slow, and the higher the condition number of the coefficient matrix is, the more iteration times are needed to reach a certain precision. The iteration is subjected to Fourier analysis to find that the high-frequency error component decays rapidly in the iteration process, the low-frequency error component decays very slowly, and the high-frequency oscillation error is recognized to be a local behavior and is derived from mutual coupling among a plurality of nearby grid points, and is irrelevant to the boundary or grid point information with a longer distance; whereas the low frequency smoothing error is global behavior, mainly derived from boundary information. The traditional relaxation method is a method with stronger locality, so that the local high-frequency oscillation error can be quickly smoothed, but the global low-frequency smoothing error is slowly attenuated, so that the error is smooth after the initial several iterations. The first proposed grid sequence or multiple grid method (Multigrid Method) in Fdeorenko, "The Speed of Convergence of One Iterative Process," USSR Computational Mathematics and Mathematical Physics, vol. 4, no. 3,1964. in 1964 is a very efficient accelerated convergence method, and in 1979 Brnadt , "Multi-Level Adaptive Solutions to Boundary-Value Problems," Mathematics of Computation, vol. 31, no. 138, 1977. the Fedorenko result was used to solve the elliptic equation, the high frequency error can be quickly eliminated by the iterative format, and the remaining low frequency error is the limiting convergence speed. Brnadt adopts a series of grids which are coarsened step by step, firstly iterates on the denser grids to eliminate high-frequency errors, then transmits the solution to the coarse grids through interpolation, and the low-frequency errors remained by the fine grids form high-frequency errors relative to the coarse grids, so that a part of errors can be eliminated further through iteration. Thus, the error of different frequencies is eliminated by iteration on a series of grids with different thicknesses, and then the value with the low frequency error eliminated is interpolated on the fine grids, so that the convergence speed is greatly improved. Jameson in 1991 applied the multiple mesh technique to the structural mesh finite volume method in CFD at AIAA 91-1956,"Time Dependent Calculations Using Multigrid, with Applications to Unsteady Flows Past Airfoils and Wings",. The geometric multiple grid has the advantages of clear structure level, very clear and regular structural relation among points, elements and elements, very clear topological relation among grids of each layer, very easy control of element number ratio between thick and thin grids and easy realization of data transmission relation between two adjacent grids, so that the multiple grid method on the structural grid is widely applied to CFD.

However, in the field of computational electromagnetics, the grid sequence method or the geometric multiple grid method is not much researched, and is currently applied to the finite element method, for example, ,IEEE TRANSACTIONS ON MAGNETICS, VOL. 42, NO. 4 "A Multigrid Solver for Time Harmonic Three-Dimensional Electromagnetic Wave Problems", electric wave science report 2013 1 st year in 2006, the known method is applied to the electromagnetic scattering problem in the composite grid method, the region is divided into two regions of coarse grids and fine grids, and the finite element is adopted to calculate the electromagnetic scattering of the fine grid region containing the microstructure. The moment method multiple grid method adopting the basis function is applied less, a small amount of documents such as the MGFFT method for analyzing electromagnetic scattering problem of the university of electronic technology journal 1999 3, long Yi are adopted, and the moment method is combined with the fast moment vector product of Fourier transformation to calculate the solution of the metal strip integral equation. In the paper of An Yuyuan, the radar scattering section of a rectangular flat plate and an H-plane corner reflector model is calculated by adopting a multi-grid pretreatment and RWG basis function moment method through two subdivision of a triangular coarse surface element. In 2008 'multi-layer grid successive approximation calculation electromagnetic field based on FDTD', the method for calculating the electromagnetic field based on FDTD by successive approximation of coarse to fine grids is applied to electromagnetic field calculation of a built-in metal resonant cavity with cracks. Chatterjee in 2015: "A Multilevel Numerical Approach with Application in Time-Domain Electromagnetics", proposes an algebraic multiple grid composed of different precision spatial formats combined with explicit RK time-marching FVTD and applied to two-dimensional cylindrical and airfoil electromagnetic scattering.

Electromagnetic computing FVFD in combination with geometric grid sequence acceleration has not been reported so far. The frequency domain micro-classification limited volume method in the literature comprises the following steps: in Huh AIAA-92-0453："a compact high-order finite-volume time-domain/frequency-domain method for electromagnetic scattering" of 1992, a Huh method is proposed, the Huh method adopts a compact differential combined filtering artificial viscous construction flux, a point hidden 4-step Runge-Kutta method is adopted for time iteration, and the process is complex. Secondly, the Bonnet method was proposed in "Frequency-Domain Finite Volume Method for Electromagnetic Scattering" in 1998, which adopted the BICGSTAB (1) method of solving a linear algebraic equation set. Thirdly, in the prior art, some calculation states in practical application are found, and the phenomenon of slow convergence due to long-time oscillation exists.

In summary, in order to improve the electromagnetic field calculation efficiency of large grid quantity and complex electromagnetic engineering problems, accelerate the iterative convergence speed of the electromagnetic field and achieve the purpose that FVFD numerical tools solve more practical electromagnetic engineering requirements, the invention provides a corresponding solution.

Disclosure of Invention

The invention aims to overcome the defects of the prior art, provides a method for accelerating a finite volume of a structural grid frequency domain electromagnetic field grid sequence, realizes the solving of the large-scale electromagnetic scattering problem of a complex-appearance, high-frequency and electric large-size conductor target by a high-efficiency frequency domain electromagnetic numerical method, and improves the calculation efficiency on the premise of ensuring high calculation precision; the invention provides a new electromagnetic numerical method for eliminating numerical errors of different frequencies in an iteration process through grid sequences with different geometric dimensions and thicknesses, and achieving an accelerating convergence effect through the grid sequence circulation.

The invention aims at realizing the following scheme:

a method for accelerating finite volume by using a structural grid frequency domain electromagnetic field grid sequence comprises the following steps:

step 1: according to the physical background of the electromagnetic problem simulated by the target, combining boundary condition information to perform simulation modeling;

step 2: meshing the simulation model by adopting a two-dimensional quadrilateral or three-dimensional hexahedral structure, encrypting the grid at the wall surface and geometric singular positions, gradually keeping the grid away from the scattering wall surface and gradually sparsely; calculating grids of the corresponding areas by numerical values, outputting grid data files, and setting and outputting boundary condition files;

step 3: inputting target calculation electromagnetic parameters and numerical calculation control parameters;

step 4: initializing a calculation space electromagnetic field;

Step 5: performing iterative solution on the frequency domain electromagnetic field of the Maxwell equation set based on virtual time propulsion and space flux residual error division calculation, and performing steady virtual time step circulation on the simulation model until convergence is finished; in each virtual time iteration process, according to a grid sequence loop mode, carrying out space flux and residual calculation on each grid block grid and each grid unit of each structure grid of the hierarchical grid in sequence, carrying out implicit iterative calculation, carrying out interpolation transfer on adjacent grid layer flux residual, conservation variable, forcing function and electromagnetic parameter, and updating the conservation electromagnetic field value of the next-stage virtual time sub-iteration step number;

Step 6: and outputting the real part and the imaginary part of the electromagnetic field, and outputting the surface induced current and the radar scattering cross section spatial distribution data.

Further, in the step 2, the density of the grid ensures 13-20 grid points per wavelength, the density of the wall surface is >300 points/wavelength, and the geometric singular point is encrypted to 50-100 grid points/wavelength.

Further, in the step 2, the following processing is performed for the two-dimensional grids in the grids: and (3) pushing a layer of the two-dimensional grid on a plane perpendicular to the two-dimensional grid according to a right-hand rule, and taking the layer of the two-dimensional grid as a three-dimensional problem special case for unified calculation.

Further, in the step 2, the grid data file includes the number of structural grid blocks and three lower dimensions of the curve coordinate system of each block.

Further, in the step3, if there is an external flow field condition of the plasma, the step3 further includes the following sub-steps: inputting flow field parameters corresponding to the external flow field of the plasma.

Further, in the step 5, a steady virtual time step is circulated for the simulation model until convergence is finished; in each virtual time iteration process, according to a grid sequence loop mode, space flux and residual calculation are sequentially carried out on each grid block grid and each grid unit of each structure grid of the hierarchical grid, implicit iteration calculation is carried out, adjacent grid layer flux residual, conservation variable, forcing function and electromagnetic parameter interpolation transmission are carried out, and the next-stage virtual time sub-iteration step number conservation electromagnetic field value is updated, and the method specifically comprises the following sub-steps:

step 5-1: and (3) stabilizing the virtual time step circulation until the calculation convergence is finished:

；(1)

；(2)

；(3)

Wherein, Is a conservation variable of the scattered electromagnetic field,Is a virtual time period of time and,Is the frequency of the incident simple harmonic electromagnetic wave,、、Is electromagnetic flux under rectangular coordinate system、、The component(s) of the composition,Is a complex type scattering magnetic induction intensity vector,Is the complex type scattered field electric displacement vector,Is the complex type scattered field electric field intensity vector,Is a complex type scattered field magnetic field intensity vector and contains subscripts、、Rectangular coordinate system with scalar values respectively corresponding to vectors、、A component;

grid space of curve coordinate system, and steady iteration adopts implicit algorithm:

；（4）

Wherein, For the structural grid curve coordinate system direction 1,For the structural grid curve coordinate system direction 2,The direction 3 of a structural grid curve coordinate system; are respectively corresponding to a curve coordinate system Directional electromagnetic flux; Is implicit control parameter, fetch Other parameters correspond to the explicit and implicit mixed formats; subscript ofIs the grid cell number and,Is the firstGrid cell numberThe electromagnetic conservation variable of the virtual time step,Is the firstGrid cell numberThe electromagnetic conservation variable of the virtual time step,Is the firstThe electromagnetic conservation variable of the virtual time step,Is the firstThe electromagnetic conservation variable of the virtual time step,Is the firstGrid cell numberThe spatial flux residual of the virtual time step,Is the firstGrid cell numberThe spatial flux residual of the virtual time step,The virtual time step is controlled by stability and calculated by CFL number, local grid cell geometric scale and characteristic value;

step 5-2: in each virtual time iteration process, in the structural grid sequence loop, each level grid carries out space flux calculation and implicit iteration solution calculation by grid blocks and grid units, and the value of the conservation electromagnetic field of the iteration step number of the next level virtual time is updated; in the structural grid sequence loop, the spatial flux in each level grid is calculated as follows:

Calculating the interface flux of the grid cells by adopting Steger-Warming splitting;

；

；(5)

In the subscript Respectively taking curve coordinate systemsOne of the directions, corresponding toNamely is corresponding toThe electromagnetic flux in the direction is directed to,Representing a curved coordinate systemIn the corresponding direction Steger-Warming splitting, the electromagnetic flux obtained after the positive eigenvalue is split; Representing a curved coordinate system In the corresponding direction Steger-Warming splitting, the electromagnetic flux obtained after the negative eigenvalue is split;， in the form of a similarity matrix, Respectively a diagonal matrix formed by positive and negative eigenvalues,，Representing left and right state variables at the interface respectively, which can adopt MUSCL format to reach the highest third-order precision;

；

Wherein, Is a limiter, subscriptIs the grid cell number and,Corresponding to the cell interface plane,Is a control parameter in a 3-order precision format,AndThe post-difference and pre-difference operators, respectively; Represented in grid cells The left state electromagnetic field conservation variable at the interface,Represented in grid cellsA right state electromagnetic field conservation variable at the interface; Is the first The individual grid cells scatter electromagnetic field conservation variables,Is the firstA grid cell scattered electromagnetic field conservation variable;

In the sequence loop of the structural grid, space flux implicit iteration and the iterative solution of the Jacobian coefficient matrix before and after splitting are adopted in each level grid, the Jacobian coefficient Steger-Warming splitting generated by flux bias conservation variable is adopted to obtain,

；(6)

；

Wherein, Is a matrix of coefficients after the splitting and,Is the spatial flux residual calculated in the last iteration time step,Is the implicit virtual time iterative electromagnetic field difference;

Will be Expressed as LDU approximate factorization

； (7)

Wherein the subscriptIs the grid cell number and,Is the maximum eigenvalue split parameter of the jacobian coefficient matrix,Is the biggest eigenvalue of the Jacobi coefficient matrix; is a matrix of unit diagonals, Is a diagonal matrix of the type,In the form of an upper triangular matrix, the upper triangular matrix,In the form of a lower triangular matrix,Is the difference of the electromagnetic conservation variables corresponding to the upper triangular matrix,Is the difference value of the electromagnetic conservation variable corresponding to the lower triangular matrix; Finger means Adjacent iteration time steps of the grid unit are electromagnetic conservation variable difference values; Finger means Adjacent iteration time steps of the grid unit are electromagnetic conservation variable difference values; Finger means Adjacent iteration time steps of the grid unit are electromagnetic conservation variable difference values; Finger means Adjacent iteration time steps of the grid unit are electromagnetic conservation variable difference values; Finger means Adjacent iteration time steps of the grid unit are electromagnetic conservation variable difference values; Finger means Adjacent iteration time steps of the grid unit are electromagnetic conservation variable difference values; Refers to adjacent A coefficient matrix after grid cell splitting; Refers to Coefficient matrix after splitting adjacent grid cells; Refers to adjacent A coefficient matrix after grid cell splitting; Refers to adjacent A coefficient matrix after grid cell splitting; Refers to Coefficient matrix after splitting adjacent grid cells; Refers to Coefficient matrix after splitting adjacent grid cells;

finally, the step of obtaining the product, obtaining the difference value of the forward and backward iterative calculation electromagnetic field ：

；

Wherein, Is a diagonal matrixIs used for the inverse matrix of (a),、Respectively according to、An upper triangular matrix and a lower triangular matrix are calculated;

forward cycle: ；

Backward circulation: ；

Wherein, Is the intermediate transition variable of the electromagnetic conservation variable difference value;

In the structural grid sequence circulation process, the grid sequence parameters which are set in advance are applied in the process, wherein the grid sequence parameters comprise the maximum number of multiple grids, the iteration number of each grid and the grid sequence circulation mode; wherein, the adjacent level grid calculating step:

Step one, solving a discrete equation in the iteration of the closest layer grid to obtain electromagnetic field quantity Increment of electromagnetic field quantityAnd residual errors, storing the amounts and transmitting the electromagnetic field amounts and residual errors to the coarse grid through a limiting operator;

；(8)

； (9)

； (10)

In the subscript Representing the geometric dimensions of the dense grid, subscriptsRepresenting the geometric scale of the thin grid; subscript ofTaking 4 in two dimensions and 8 in three dimensions for obtaining sum variables; Is that The triangular matrix on the level grid is used for generating a three-dimensional matrix,Is thatLevel grid diagonal matrixIs used for the inverse matrix of (a),Is thatA triangular matrix under the grid of the level,Is a virtual time step size of the time,Is thatStage grid space flux residual; Is that The initial value of the electromagnetic field quantity on the stage grid,Is thatStage grid toThe volume weighted interpolation operator of the stage mesh,Is thatThe volumes of grid cells on the grid of stages that participate in the summation,To pass the resulting residual from the most dense grid,Is thatThe level grid spatial flux residual is calculated,Is thatResidual errors participating in summation on the stage grid;

step (II), transferring to the iteration of the sub-dense grid layer, and obtaining a forcing function in the first step of iteration The meaning is that the difference between the primary field and the transmitted residual error on the obtained secondary dense grid is limited:

；(13)

In the method, in the process of the invention, To pass the resulting residual from the most dense grid,Iteratively calculating a residual error calculated by an explicit format for the first step of coarse mesh; in an iteration of the layer of the present system,Will remain unchanged; if a W loop is employed, then, in the downward calculation,The previous value will be taken; the actual residual of the iteration on the coarse grid is the residual calculated by the format plus a forcing term:

； (11)

Is that On the level gridThe actual residual of the next sub-iteration,Is thatResidual errors obtained by calculation of space format flux on the grade grid, and electromagnetic field variables after coarse grid iteration are recorded asRepeating the step (I);

Step (III), transferring to a thinner grid to continue calculation, and repeating the steps (one) - (three) until the most thin grid is reached;

Step four, after iterative calculation of the thinnest grid, sequentially inserting correction amounts back into the upper grid until the grid is the thickest; the back-inserting process carries out 'post iteration' on the grid of the layer, and 'forcing item' is added to the iteration at the moment and is used for avoiding the uploading of high-frequency error components;

；(12)

； (13)

In the method, in the process of the invention, Is thatThe amount of electromagnetic field on the mesh of the stage,To transfer interpolation operators for correction from the lower mesh to the upper mesh,Is thatElectromagnetic field volume increment on a stage grid，Is thatThe first iteration on the mesh of stages is followed by the electromagnetic field quantity,Is thatA primary value of the electromagnetic field quantity of the grade grid; the method adopts a tri-linear interpolation mode, and is characterized in that the correction quantity is firstly interpolated from the grid center to the nodes through volume weighting, and then the correction quantity is linearly interpolated to the dense grid center.

Further, subscriptsDilute grid geometry scale relative subscriptIs a dense grid geometry, each dimension of the grid is halved.

Further, the lower layer grid is more specifically a thin grid to a dense grid.

Further, in the step 5, the steady virtual time step loop is implicit, and the CFL number is not limited by the explicit stability requirement.

Further, in the step 6, the iterative convergence criterion is that the second moment of the magnetic field drops to a standard designated in advance when the absolute value of the spatial amplitude difference is calculated, 8.e-5 is selected in two dimensions, and 5.e-4 is selected in three dimensions.

The beneficial effects of the invention include:

1. The invention can solve the problem of large-scale electromagnetic scattering of any complicated shape and high-frequency electric large-size targets. The method adopts a multi-block structural grid and geometric grid sequence iterative convergence acceleration algorithm in the specific conception, and can solve the technical problem that the prior application scheme has long-time oscillation and slow convergence in some calculation states in practical application.

Specifically, the invention directly numerically simulates a frequency domain electromagnetic field Maxwell equation set, carries out flux residual error, conservation variable, forcing function and electromagnetic parameter interpolation transmission of adjacent grid layers at each grid level in a specific circulation mode in virtual time iteration, firstly iterates at the most dense grid to eliminate high-frequency errors, and then eliminates low-frequency errors at the coarse grid, so that the errors of different frequencies are iterated and eliminated at a series of grids with different thicknesses, and the purpose of accelerating convergence is achieved.

The calculation example of the embodiment of the invention shows that the convergence curve display adopts multiple grid sequence iteration to converge faster than single grid calculation and the W shape is optimal in the whole complex electromagnetic field steady calculation regardless of the V shape or the W shape circulation. The two-dimensional metal square column calculation example shows that: the V-shaped grid sequence is accelerated by 60.8% compared with a single grid, so that the efficiency is improved by 155%; the W-shaped grid sequence is 73.4% faster than single grid convergence, and the efficiency is improved by 227.8%; three-dimensional football examples show that: the V-shaped grid sequence is accelerated by 70.3% compared with a single grid, so that the efficiency is improved by 236.8%; the W-shaped grid sequence is 75.5% faster than single grid convergence, and the efficiency is improved by 308.5%.

Specifically, the calculation accuracy of the method FVFD is not affected by the iterative convergence mode, and typical examples show that FVFD has high calculation accuracy. Because the V-shaped grid sequence and the W-shaped grid sequence have quick dissipation effect on oscillation errors of various frequencies, the electromagnetic fields at the same space position are calculated and compared, and the oscillation amplitude of the numerical value is smaller and the oscillation amplitude is faster to converge than that of the single grid, wherein the W-shaped grid sequence is optimal.

Specifically, the structural grid and geometric grid sequence in the method of the invention are simple to realize, have clear structure, can select a proper number of grid levels, and can also arbitrarily select the circulation mode of the grid sequence, such as W-shaped circulation or V-shaped circulation, so as to achieve the optimal acceleration efficiency.

Specifically, in the method, an incident field is given by analysis, a scattered field form Maxwell equation set is solved, propagation of an incident electromagnetic wave is not required to be calculated in the whole calculation grid space, dissipation and dispersion of the incident electromagnetic field are avoided, and numerical accuracy is maintained.

2. According to the invention, the conservation electromagnetic field increment is calculated by adopting a forward and backward iteration implicit inversion algorithm, on one hand, the full hidden format relaxation stability limit can obtain a large iteration time step, and in addition, the two scans are simplified to replace sparse matrix inversion, so that the memory and the operand are saved, the defect of large calculation amount caused by the traditional explicit limited-volume global time step limit is overcome, the engineering is simple and easy, and the calculation performance is improved.

3. The method can solve the large-scale electromagnetic scattering problem of any complicated shape and high-frequency electric large-size targets or solve the electromagnetic problem that the multi-scale MOM containing a complicated electronic detail structure is difficult to process; and, support the structural grid of the curve coordinate system and multizone decompose and parallel algorithm; for the simple harmonic single-frequency incidence condition (the time signal is a continuous periodic signal at the moment), in order to improve the calculation efficiency and the precision and reduce the grid quantity, the time is directly derived and converted into a frequency domain in a Maxwell equation set, so that the variable dimension is reduced from 4 dimensions in space time to 3 dimensions in only residual space, the variable is reduced from an unsteady algorithm to a steady algorithm, and the variable is complex and graceful and more efficient, thereby reducing Fourier transformation links required by calculating the target electromagnetic characteristics in the time domain.

4. The implicit frequency domain Finite Volume (FVFD) method in the method conception is different from FDTD, the FDTD adopts Cartesian orthogonal grid to simulate that the wall surface has a ladder effect to influence the numerical accuracy, and artificial viscosity is added to a 2-order central differential format by using electromagnetic field component space-time cross placement, the method FVFD adopts a body-attached curve coordinate system grid to better fit an object plane and encrypt the grid at a geometric singular position, and electromagnetic field quantity is placed in the grid unit center in the grid space, and the artificial viscosity is maintained by adopting a windward format, so that the accuracy and algorithm design are more favorably maintained.

5. Different from the finite element method, FEM and FVFD can both adopt arbitrary shape grid units to simulate discrete calculation space, FEM adopts a basic function to simulate one of node or edge vector electric field and magnetic field vector, and adopts a variation method or residual value weighting to construct a matrix form equation set, so as to obtain a strip discrete full-space matrix and solve the linear algebraic equation set.

6. Unlike FVFD of AIAA-92-0453："a compact high-order finite-volume time-domain/frequency-domain method for electromagnetic scattering" Huh in 1992, the Huh method adopts a compact differential combined filtering artificial viscosity structure flux, the time iteration adopts a point hidden 4-step Runge-Kutta method, the process is complex and complicated, and the FVFD method adopts windward interpolation Steger-Warming split structure flux and implicit iteration to calculate the electromagnetic field quantity.

7. The invention is different from FVFD of Bonnet in 'Frequency-Domain Finite Volume Method for Electromagnetic Scattering' in 1998, wherein first, the FVFD flux of the invention adopts a Steger-Warming split non-simple phase geometric relationship, second, the invention FVFD adds a virtual time step, and adopts an implicit forward and backward split matrix iterative computation, and the Bonnet adopts a method BICGSTAB (1) of solving a linear algebraic equation set.

Drawings

In order to more clearly illustrate the embodiments of the invention or the technical solutions of the prior art, the drawings which are used in the description of the embodiments or the prior art will be briefly described, it being obvious that the drawings in the description below are only some embodiments of the invention, and that other drawings can be obtained according to these drawings without inventive faculty for a person skilled in the art.

FIG. 1 is a software flow diagram of a method of an embodiment of the invention;

FIGS. 2 a-2 b are two main grid sequence loops (V-shaped or complex W-shaped), where FIG. 2a is a V-shaped grid sequence loop; FIG. 2b is a complex W-shaped mesh sequence loop pattern;

Fig. 3 a-3 d are two-dimensional metal square column 4-heavy structure grid sequence distributions, wherein fig. 3a is an ign=1-nearest network, fig. 3b is an ign=2 network, fig. 3c is an ign=3 network, and fig. 3d is an ign=4-nearest network;

FIG. 4 is a comparison of the two-dimensional metal square column electromagnetic scattering calculation 4-fold grid sequence convergence characteristics; wherein, the ordinate represents the total Residual error (Residual) of the whole grid electromagnetic field iterative calculation process, and the abscissa represents the iterative calculation Step number (Step); lu-SGS SINGLE GRID represents the calculation situation adopting a single-grid and Lu-Sgs algorithm, 4 Level V-MultiGrid represents the calculation situation adopting a 4-grid and V-shaped grid sequence circulation algorithm, and 4 Level W-MultiGrid represents the calculation situation adopting a 4-grid and W-shaped grid sequence circulation algorithm;

FIGS. 5 a-5 b are real and imaginary contour cloud diagrams of a two-dimensional metal square column scattered electromagnetic field; wherein FIG. 5a is a real contour cloud of the two-dimensional metal square-column scattered electromagnetic field, FIG. 5b is a virtual contour cloud of the two-dimensional metal square-column scattered electromagnetic field, bsz representing the z-component of the scattered magnetic induction vector;

FIG. 6 is a comparison of radar cross-section distributions calculated in 3 iterative ways for two-dimensional metal square columns; wherein the abscissa is Indicating the azimuth angle of scattering, 0 in brackets on the abscissa representing in degrees; the ordinate sigma/lambda represents the two-dimensional radar scattering width normalized by wavelength, the dB inside the ordinate bracket represents in decibels, and MOM represents the moment method;

fig. 7 a-7 c are three-dimensional football (f=1.18 GHz) 3-weight structural grid sequence distributions; wherein fig. 7a is the ign=1 nearest network, fig. 7b is the ign=2 network, and fig. 7c is the ign=3 nearest network;

Fig. 8 is a 3-fold grid sequence convergence characteristic comparison of three-dimensional football (f=1.18 GHz) electromagnetic scattering calculations; wherein, the ordinate represents the Residual error (residual_L2_Norm) represented by L2 Norm, 3 Level V-MultiGrid represents the calculation situation of 3-fold grid and V-shaped grid sequence circulation algorithm, and 3 Level W-MultiGrid represents the calculation situation of 3-fold grid and W-shaped grid sequence circulation algorithm;

Fig. 9 is a plot of the convergence history of three-dimensional football (f=1.18 GHz) at the same location in 3 iterative ways; wherein, the ordinate represents the electromagnetic field Amplitude (Amplitude), dy (Lu-SGS SINGLE GRID) represents the situation of an electric displacement vector y component calculated by adopting a single-grid and Lu-Sgs algorithm, dy (3 Level V-MultiGrid) represents the situation of an electric displacement vector y component calculated by adopting a 3-grid and V-shaped grid sequence circulation algorithm, and Dy Level W-MultiGrid represents the situation of an electric displacement vector y component calculated by adopting a 3-grid and W-shaped grid sequence circulation algorithm;

fig. 10 is a three-dimensional rugby (f=1.18 GHz) electromagnetic scattering surface induced current cloud plot; wherein Jz represents the z-component of the target surface induced current density vector;

Fig. 11 is a comparison of radar cross-section distributions (horizontal polarization) calculated in 3 iterative ways for three-dimensional rugby (f=1.18 GHz); wherein, the ordinate represents Radar Cross Section (RCS), the abscissa represents angle, SINGLE GRID FVFD represents calculation by adopting single-grid and FVFD algorithm, 3 level V_MultiGrid FVFD represents calculation by adopting 3-grid and V-grid sequence circulation FVFD algorithm, and 3 level W_MultiGrid FVFD represents calculation by adopting 3-grid and W-grid sequence circulation FVFD algorithm.

Detailed Description

All of the features disclosed in all of the embodiments of this specification, or all of the steps in any method or process disclosed implicitly, except for the mutually exclusive features and/or steps, may be combined and/or expanded and substituted in any way. An embodiment of a method for accelerating finite volume of a grid sequence of a frequency domain electromagnetic field of a structural grid according to an embodiment of the present invention is further described below with reference to the accompanying drawings.

Referring to fig. 1, the whole structure grid frequency domain electromagnetic field finite volume convergence acceleration method software can be divided into: pretreatment, electromagnetic field calculation and post-treatment. The preprocessing mainly comprises three modules, namely grid data input, calculation parameter data input and control parameter input, and is mainly used for reading the grid data, the calculation parameter data input and the control parameter file, and preprocessing is carried out on the basis, so as to provide calculation support for electromagnetic field calculation; the electromagnetic field calculation includes: the space electromagnetic field MUSCL format interpolation, unit interface flux calculation and correction, time pushing and convergence judgment module; the post-treatment is mainly used for outputting the space real part and imaginary part distribution of an electromagnetic field, the induced current density of the target surface and the radar scattering cross section output.

The two rotation equations (time factors) of the frequency domain Maxwell equation set to be numerically simulated are combined as follows) Faraday (Faraday) law of electromagnetic induction: ; ampere (amp) theorem: An implicit frequency domain finite volume method numerical computation process is introduced. Wherein the method comprises the steps of Is the imaginary symbol corresponding to the complex variable,Is the frequency of the simple harmonic electromagnetic wave,Is a complex type scattering magnetic induction intensity vector,Is the complex type scattered field electric displacement vector,Is the complex type scattered field electric field intensity vector,Is a complex type scattered field magnetic field intensity vector,Is forced current applied.

The rectangular coordinate system conservation form of the two rotation equations of the frequency domain Maxwell equation set under the passive condition is as follows:

；

Wherein, Is a conservation variable of the scattered electromagnetic field,Is a virtual time period of time and,Is the sign of the imaginary part,Is the rate of the circumference of the circle,Is the frequency of the incident simple harmonic electromagnetic wave,、、Is electromagnetic flux under rectangular coordinate system、、The component(s) of the composition,Is a complex type scattering magnetic induction intensity vector,Is the complex type scattered field electric displacement vector,Is the complex type scattered field electric field intensity vector,Is a complex type scattered field magnetic field intensity vector and contains subscripts、、Rectangular coordinate system with scalar values respectively corresponding to vectors、、A component.

Clearly visible whenThe system of equations is equivalent to the original system of equations when converging.

For complex-appearance objects, a computational space body-attached multi-block structural grid is adopted, so that coordinate transformation exists:

Wherein, Representing a curved coordinate systemThree directions are respectively takenOne of them. Obtaining the conservation shape of Maxwell equation set under the curve coordinate system of the required numerical simulation:

；

In the method, in the process of the invention, Is a jacobian matrix of coordinate transformation, and the corresponding ≡superscript variable represents a value under a curve coordinate system and is obtained by the coordinate transformation.Is an electromagnetic conservation variable under a curve coordinate system; for the structural grid curve coordinate system direction 1, For the structural grid curve coordinate system direction 2,The direction 3 of a structural grid curve coordinate system; is the electromagnetic flux in a curved coordinate system. Respectively take outOne of the directions under the three curved coordinate systems.

In order to get rid of the large calculation amount defect caused by the limitation of the traditional explicit finite volume global on the iteration time step, the finite volume convergence acceleration method of the structural grid frequency domain electromagnetic field, which is conceived by the invention, utilizes the constant iteration promotion of the local time step and combines the implicit space flux residual calculation and correction to obtain a stable and efficient calculation flow, and comprises the following steps:

Step 1: and carrying out simulation modeling according to the physical background of the electromagnetic problem simulated by the target and combining boundary condition information.

Step 2: and meshing the simulation model by adopting a quadrilateral (2-dimensional) or hexahedral (3-dimensional) structure, encrypting the grid at the wall surface and geometric singular positions, and gradually thinning the grid along with gradually being far away from the scattering wall surface. And (5) calculating grids of the corresponding areas by numerical values, outputting grid data files, and setting and outputting boundary condition files. The grid density ensures 13-20 grid points per wavelength, the wall density is greater than 300 points per wavelength, the geometric singular point is encrypted to 50-100 grid points per wavelength, the two-dimensional grid is pushed to a layer on the plane perpendicular to the two-dimensional grid according to the right-hand method, and the three-dimensional grid is used as the unified calculation of the three-dimensional problem special cases. The grid data file includes the number of structured grid blocks and three curvilinear coordinate system downdimensions per block.

Step 3: and the preprocessing part inputs the target calculation electromagnetic parameters, the numerical calculation control parameters and the corresponding flow field parameter file under the condition of the external flow field with plasma. The virtual time iteration is implicit in that its CFL number is not constrained by explicit stability requirements.

Step 4: and inputting grid data and boundary condition information files, and initializing and calculating a space electromagnetic field.

Step 5: carrying out iterative solution on the frequency domain electromagnetic field of the Maxwell equation set based on virtual time propulsion and space flux residual error division calculation and correction; performing a steady virtual time step circulation on the simulation model until convergence is finished; in each virtual time iteration process, according to a grid sequence loop mode, space flux and residual calculation are sequentially carried out on each grid block grid and each grid unit of the structure grid of the hierarchy grid, implicit iteration calculation is carried out, adjacent grid layer flux residual, conservation variable, forcing function and electromagnetic parameter interpolation transmission are carried out, and the conservation electromagnetic field value of the next-stage virtual time sub-iteration step number is updated.

In step 5, the method specifically comprises the following sub-steps:

step 5-1: and (5) stabilizing the virtual time step circulation until the calculation convergence is finished.

；(1)

；(2)

；(3)

Wherein, Is a conservation variable of the scattered electromagnetic field,Is a virtual time period of time and,Is the frequency of the incident simple harmonic electromagnetic wave,、、Is electromagnetic flux under rectangular coordinate system、、The component(s) of the composition,Is a complex type scattering magnetic induction intensity vector,Is the complex type scattered field electric displacement vector,Is the complex type scattered field electric field intensity vector,Is a complex type scattered field magnetic field intensity vector and contains subscripts、、Rectangular coordinate system with scalar values respectively corresponding to vectors、、A component.

；(4)

Wherein, For the structural grid curve coordinate system direction 1,For the structural grid curve coordinate system direction 2,The direction 3 of a structural grid curve coordinate system; are respectively corresponding to a curve coordinate system Directional electromagnetic flux; Is implicit control parameter, fetch Other parameters correspond to the explicit and implicit mixed formats; subscript ofIs the grid cell number and,Is the firstGrid cell numberThe electromagnetic conservation variable of the virtual time step,Is the firstGrid cell numberThe electromagnetic conservation variable of the virtual time step,Is the firstThe electromagnetic conservation variable of the virtual time step,Is the firstThe electromagnetic conservation variable of the virtual time step,Is the firstGrid cell numberThe spatial flux residual of the virtual time step,Is the firstGrid cell numberThe spatial flux residual of the virtual time step,The method is a virtual time step controlled by stability, is calculated by CFL number, local grid cell geometric scale and characteristic value, and is obviously different from an explicit method, and different grid cells which are calculated constantly have different local virtual time sub-iteration step sizes so as to accelerate the convergence of the electromagnetic field of the cell.

Step 5-2: in each virtual time iteration process, in the structural grid sequence loop, each hierarchical grid performs space flux calculation and implicit iteration solution calculation on a grid-by-grid block basis and a grid-by-grid unit basis, and the value of the conservation electromagnetic field of the iteration step number of the next virtual time is updated.

In the structural grid sequence loop, the spatial flux within each hierarchical grid is calculated as follows:

grid cell interface flux was calculated using the Steger-Warming split.

；

；(5)

In the subscriptRespectively taking curve coordinate systemsOne of the directions, corresponding toNamely is corresponding toThe electromagnetic flux in the direction is directed to,Representing a curved coordinate systemIn the corresponding direction Steger-Warming splitting, the electromagnetic flux obtained after the positive eigenvalue is split; Representing a curved coordinate system In the corresponding direction Steger-Warming splitting, the electromagnetic flux obtained after the negative eigenvalue is split;， in the form of a similarity matrix, Respectively a diagonal matrix formed by positive and negative eigenvalues,，Representing the left and right state variables at the interface, respectively, which can be used in MUSCL format to achieve the highest third order accuracy.

；

Wherein the method comprises the steps ofIs a limiter, subscriptIs the grid cell number and,Corresponding to the cell interface plane,Is a control parameter in a 3-order precision format,AndThe post-difference and pre-difference operators, respectively; Represented in grid cells The left state electromagnetic field conservation variable at the interface,Represented in grid cellsA right state electromagnetic field conservation variable at the interface; Is the first The individual grid cells scatter electromagnetic field conservation variables,Is the firstThe individual grid cells scatter electromagnetic field conservation variables.

；(6)

；

Wherein the method comprises the steps ofIs a matrix of coefficients after the splitting and,Is the spatial flux residual calculated in the last iteration time step,Is an implicit virtual time iterative electromagnetic field difference.

Will beExpressed as LDU approximate factorization

；(7)

Wherein the subscriptIs the grid cell number and,Is the maximum eigenvalue split parameter of the jacobian coefficient matrix,Is the biggest eigenvalue of the jacobian coefficient matrix.Is a matrix of unit diagonals,Is a diagonal matrix of the type,In the form of an upper triangular matrix, the upper triangular matrix,In the form of a lower triangular matrix,Is the difference of the electromagnetic conservation variables corresponding to the upper triangular matrix,Is the difference value of the electromagnetic conservation variable corresponding to the lower triangular matrix; Finger means Adjacent iteration time steps of the grid unit are electromagnetic conservation variable difference values; Finger means Adjacent iteration time steps of the grid unit are electromagnetic conservation variable difference values; Finger means Adjacent iteration time steps of the grid unit are electromagnetic conservation variable difference values; Finger means Adjacent iteration time steps of the grid unit are electromagnetic conservation variable difference values; Finger means Adjacent iteration time steps of the grid unit are electromagnetic conservation variable difference values; Finger means Adjacent iteration time steps of the grid unit are electromagnetic conservation variable difference values; Refers to adjacent A coefficient matrix after grid cell splitting; Refers to Coefficient matrix after splitting adjacent grid cells; Refers to adjacent A coefficient matrix after grid cell splitting; Refers to adjacent A coefficient matrix after grid cell splitting; Refers to Coefficient matrix after splitting adjacent grid cells; Refers to Coefficient matrix after splitting adjacent grid cells;

；

forward cycle: ；

Backward circulation: ；

Wherein, Is the intermediate transition variable of the electromagnetic conservation variable difference.

In the process of the structural grid sequence loop (V loop), the grid sequence parameters which are set in advance are applied in the process, wherein the grid sequence parameters comprise the maximum number of multiple grids (as shown in fig. 3 a-3 d, the maximum effect is generally good up to 4), the iteration number of each grid, and the grid sequence loop mode (as shown in fig. 2 a-2 b, V-shaped or complex W-shaped). And calculating adjacent level grids:

Firstly, solving a discrete equation in the iteration of the closest-layer grid to obtain electromagnetic field quantity Increment of electromagnetic field quantityAnd a residual, storing these quantities and passing the electromagnetic field quantities and residual to the coarse grid through a "limiting operator".

；(8)

； (9)

；(10)

In the subscriptRepresenting the geometric dimensions of the mesh (dense mesh), subscriptsRepresents the geometric scale of the mesh (thin mesh, relative to the preceding subscriptDense grid, halving each dimension of the grid), subscriptTaking 4 in two dimensions and 8 in three dimensions for obtaining sum variables; Is that The triangular matrix on the level grid is used for generating a three-dimensional matrix,Is thatLevel grid diagonal matrixIs used for the inverse matrix of (a),Is thatA triangular matrix under the grid of the level,Is a virtual time step size of the time,Is thatStage grid space flux residual; Is that The initial value of the electromagnetic field quantity on the stage grid,Is thatStage grid toThe volume weighted interpolation operator of the stage mesh,Is thatThe volumes of grid cells on the grid of stages that participate in the summation,To pass the resulting residual from the most dense grid,Is thatThe level grid spatial flux residual is calculated,Is thatResidual errors on the stage grid that participate in the summation.

(II) transferring to a next-secret grid layer iteration, and in the first step of iteration, obtaining a forcing functionThe meaning is that the difference between the primary field and the transmitted residual error on the obtained secondary dense grid is limited:

；(13)

In the method, in the process of the invention, To pass the resulting residual from the most dense grid,The residuals calculated by the explicit format are iteratively calculated for the first step of coarse mesh. In an iteration of the layer of the present system,Will remain unchanged; if a W loop is employed, then, in the downward calculation,The previous values will be used. The actual residual of iteration on the coarse grid is the residual calculated by the format plus a forcing term

； (11)

Is thatOn the level gridThe actual residual of the next sub-iteration,Is thatResidual errors obtained by calculation of space format flux on the grade grid, and electromagnetic field variables after coarse grid iteration are recorded asRepeating the step (one).

And (III) turning to thinner grid calculation, repeating the steps (I) - (III) until the most thin grid is reached.

And (IV) after iterative calculation of the most dilute grid, sequentially inserting correction amounts back to the upper grid until the most dense grid. The back-interpolation process may further perform "post-iteration" on the grid of the layer in order to avoid uploading the high-frequency error component, and the "forcing term" is still added to the iteration.

；(12)

； (13)

In the method, in the process of the invention,Is thatThe amount of electromagnetic field on the mesh of the stage,To transfer interpolation operators for modifier amounts from the lower grid to the upper grid (thin grid to dense grid),Is thatElectromagnetic field volume increment on a stage grid，Is thatThe first iteration on the mesh of stages is followed by the electromagnetic field quantity,Is thatPrimary electromagnetic field magnitude of the stage mesh.The method adopts a tri-linear interpolation mode, and is characterized in that the correction quantity is firstly interpolated from the grid center to the nodes through the volume weighting, and then the correction quantity is linearly interpolated to the grid center of the dense grid.

The above is an iterative process of the structure grid, grid sequence acceleration implicit FVFD to calculate the frequency domain electromagnetic field controlled by the Maxwell equation set.

Step 6: and (3) performing convergence judgment, performing post-processing, outputting spatial distribution of a real part and an imaginary part of an electromagnetic field, outputting spatial distribution data of surface induced current and radar scattering cross sections, and the like.

As shown in fig. 2 a-2 b, two main grid sequence calculation cycle modes, V-shaped cycle and W-shaped cycle, are respectively implemented by first rapidly filtering high-frequency errors from the denser grid, and filtering smooth high-frequency errors by calculating the geometric coarse-fine grid generated by semi-coarsening.

As shown in FIG. 3 a-FIG. 3d, the two-dimensional metal square column 4 is distributed in a grid sequence with a heavy structure, the length of the scatterer is 10 times of the wavelength, the width of the scatterer is 1 wavelength, and the incident electromagnetic wave forms an angle of 20 degrees with the X axis and is polarized by TE. The grid data dimension meets the half coarsening requirement of the grid sequence, namely the two space direction dimensions of the space grid block are，The two-dimensional 3 rd direction is a layer of units. In this example, the two-dimensional space has 2 grid blocks, the grid dimensions are all 201X353X2, the far-field boundary is outside 10 wavelengths, the normal direction is totally 201 grid points, the wall is encrypted to 300 grid points per wavelength, then gradually sparse to 4 grid points per wavelength along the radiation direction, and the other two directions are uniformly distributed, and the average of the wall is about 30 grid points per wavelength.

As shown in fig. 4, the two-dimensional metal square column electromagnetic scattering calculation is performed to obtain a 4-fold grid sequence convergence characteristic comparison, the grid sequence adopts grids shown in fig. 3 a-3 d, the calculation condition is that the maximum number of the grids=4, the cfl number=10, and each level of grids is calculated for 5 times in an iterative manner. The calculation shows that the 4-fold grid sequence iteration is faster than the single grid calculation in convergence no matter the V-shaped or W-shaped loop, the whole calculation is finished after the convergence standard (residual error < 8.E-5) is met, 357 iterations are needed for single grid calculation, 140 iterations are needed for the V-shaped grid sequence, 115 iterations are needed for the W-shaped grid sequence, and the V-shaped grid sequence in the calculation example is accelerated by 60.8% compared with the single grid convergence, and the efficiency is improved by 155%; the W-shaped grid sequence is 73.4% faster than single grid convergence, and the efficiency is improved by 227.8%.

As shown in fig. 5 a-5 b, the real part and the imaginary part of the electromagnetic scattering field are the real part and the imaginary part of the contour cloud pictures of the two-dimensional metal square column scattering electromagnetic field, and the real part and the imaginary part cloud pictures of the grid space complex electromagnetic scattering field can be distributed.

As shown in fig. 6, comparing the radar cross section distribution calculated by the two-dimensional metal square column 3 iteration modes, whether the grid sequence is adopted to accelerate the iteration convergence or not can be seen, comparing the double-station radar cross section distribution with the moment method, the peak and trough and the numerical value are well matched, the calculation accuracy is ensured without being influenced by the iteration convergence mode, and the high calculation accuracy is also shown by FVFD.

As shown in fig. 7 a-7 c, the three-dimensional rugby (f=1.18 GHz) 3-fold grid sequence distribution is an EMCC standard example, and the frequency of the incident electromagnetic wave is 1.18GHz, and the grid sequence distribution is horizontally polarized. The calculation example comprises 1 grid block, the dimension is 41x49x65, a full implicit format is adopted, wherein 20 grid points per wavelength are selected by object plane grids, the far field boundary is outside 3 wavelengths, and the radial grid wall surface is encrypted.

As shown in fig. 8, the three-dimensional rugby (f=1.18 GHz) electromagnetic scattering calculation is performed to calculate the convergence characteristic comparison of the 3-fold grid sequence, the grid sequence adopts the three-dimensional rugby (f=1.18 GHz) 3-fold structure grid sequence to distribute grids, the calculation condition is that the maximum number of the grids=3, the cfl number=100, and each level of grids is calculated for 8 times in an iterative manner. The calculation shows that the 3-fold grid sequence iteration is faster than the single grid calculation in convergence no matter the V-shaped or W-shaped loop, the whole calculation is finished after the convergence standard (residual error <5. E-4) is met, the single grid calculation needs 3072 iterations, the V-shaped grid sequence needs 912 iterations, the W-shaped grid sequence needs 752 iterations, and the V-shaped grid sequence in the calculation example is accelerated by 70.3% compared with the single grid convergence, and the efficiency is improved by 236.8%; the W-shaped grid sequence is 75.5% faster than single grid convergence, and the efficiency is improved by 308.5%.

As shown in fig. 9, the convergence histories of the induction intensities Dy of the three-dimensional rugby (f=1.18 GHz) in the same position in the 3 iterative manner are obviously shown, and due to the rapid dissipation effect of the V-shaped and W-shaped grid sequences on the oscillation errors of various frequencies, the oscillation amplitudes of the V-shaped and W-shaped grid sequences are smaller and the convergence is faster than that of the single grid case without the iteration of the grid sequences, wherein the W-shaped grid sequences are optimal.

Fig. 10 is a three-dimensional football (f=1.18 GHz) electromagnetic scattering surface induced current, one of the outputs of the present invention.

Fig. 11 is a comparison of radar cross-section distributions (horizontal polarization) calculated in 3 iterative ways for three-dimensional rugby (f=1.18 GHz). The three are almost completely matched, and compared with the moment method with highest precision at present, even if the radar scattering cross section is extremely small, most angles of errors are smaller than 1 dB, the calculation precision of the output of the whole radar scattering cross section can be ensured by the iteration of the V-shaped grid sequence and the W-shaped grid sequence, and the reliability is ensured.

The two-dimensional and three-dimensional numerical value examples show that the corresponding structure grid frequency domain electromagnetic field grid sequence acceleration FVFD algorithm is combined with the total hidden forward and backward iteration implicit inversion algorithm, so that the format numerical value precision can be ensured, the space distribution of the target complex electromagnetic scattering field can be accurately obtained while the virtual time step limit is relaxed, the virtual time propulsion iteration convergence efficiency is improved, and the calculation cost is saved.

In summary, the invention relates to improvement based on an implicit frequency domain finite volume method FVFD (Finite Volume Frequency Domain, FVFD), provides a novel numerical simulation method for calculating the electromagnetic field distribution of a grid space frequency domain, and contributes to the field how to improve the grid sequence iterative acceleration skill of the steady complex type scattered field iterative convergence efficiency. The method of the corresponding embodiment provides a high-efficiency grid sequence acceleration method for directly solving the geometric structure grid of the frequency domain Maxwell equation set, and forms a new method containing acceleration convergence skill FVFD. Unlike finite element and moment methods, the present invention FVFD has no large storage requirement and large dense matrix iterative computation, FVFD only requires iterative computation adjacent to several grid cell variables; unlike the FDTD multilayer grid method, the method adopts an implicit finite volume method with conformal grids and co-arranged electromagnetic field quantities; unlike the algebraic multiple grid of Chatterjee in 2015, the invention combines successful experience in CFD field, and accelerates convergence only on the geometric sequence of structural grid, thereby greatly improving the application efficiency of FVFD in various large-scale electromagnetic problems.

The units involved in the embodiments of the present invention may be implemented by software, or may be implemented by hardware, and the described units may also be provided in a processor. Wherein the names of the units do not constitute a limitation of the units themselves in some cases.

According to other aspects of embodiments of the present invention, there is provided a computer program product or computer program comprising computer instructions stored in a computer readable storage medium. The computer instructions are read from the computer-readable storage medium by a processor of a computer device, and executed by the processor, cause the computer device to perform the methods provided in the various alternative implementations described above.

In addition to the foregoing examples, those skilled in the art will recognize from the foregoing disclosure that other embodiments can be made and in which various features of the embodiments can be interchanged or substituted, and that such modifications and changes can be made without departing from the spirit and scope of the invention as defined in the appended claims.

Claims

1. The method for accelerating the finite volume of the grid sequence of the frequency domain electromagnetic field of the structural grid is characterized by comprising the following steps of:

step 4: initializing a calculation space electromagnetic field;

Step 6: outputting the real part and the imaginary part of the electromagnetic field, and outputting the surface induced current and the radar cross section spatial distribution data;

In the step 5, the simulation model is circulated in a steady virtual time step until convergence is finished; in each virtual time iteration process, according to a grid sequence loop mode, space flux and residual calculation are sequentially carried out on each grid block grid and each grid unit of each structure grid of the hierarchical grid, implicit iteration calculation is carried out, adjacent grid layer flux residual, conservation variable, forcing function and electromagnetic parameter interpolation transmission are carried out, and the next-stage virtual time sub-iteration step number conservation electromagnetic field value is updated, and the method specifically comprises the following sub-steps:

；(1)

；(2)

；(3)

；（4）

；

；(5)

；

；(6)

；

Will be Expressed as LDU approximate factorization

；(7)

；

forward cycle: ；

Backward circulation: ；

；(8)

； (9)

； (10)

；(13)

； (11)

step four, after iterative calculation of the thinnest grid, sequentially inserting correction amounts back into the upper grid until the grid is the thickest; the back-inserting process is carried out on the grid of the layer, and a forcing term is added to the iteration at the time and is used for avoiding uploading of high-frequency error components;

；(12)

； (13)

2. The method for accelerating finite volume according to claim 1, wherein in the step 2, the density of the grids ensures 13-20 grid points per wavelength, the density of the wall surface is >300 points/wavelength, and the geometric singular point is encrypted to 50-100 grid points/wavelength.

3. The method for accelerating finite volume according to claim 1, wherein in step 2, the following is performed for two-dimensional grids of the grids: and (3) pushing a layer of the two-dimensional grid on a plane perpendicular to the two-dimensional grid according to a right-hand rule, and taking the layer of the two-dimensional grid as a three-dimensional problem special case for unified calculation.

4. The method for accelerating finite volume according to claim 1, wherein in step 2, the grid data file includes the number of grid blocks of the structure and the lower dimension of three curve coordinate systems per block.

5. The method for accelerating finite volume according to claim 1, wherein in step 3, if there is an external flow field condition of the plasma, step 3 further comprises the following sub-steps: inputting flow field parameters corresponding to the external flow field of the plasma.

6. A method for accelerating finite volume of a structured grid frequency domain electromagnetic field grid sequence as set forth in claim 1 wherein the subscriptsDilute grid geometry scale relative subscriptIs a dense grid geometry, each dimension of the grid is halved.

7. The method for accelerating finite volume of a structured grid frequency domain electromagnetic field grid sequence according to claim 1, wherein the lower grid is a thin grid to a dense grid.

8. The method for accelerating finite volume according to claim 1, wherein in step 5, the steady virtual time step cycle is implicit, and the CFL number is not constrained by explicit stability requirements.

9. The method for accelerating finite volume of the electromagnetic field grid sequence in the structural grid frequency domain according to claim 1, wherein in the step 6, the iterative convergence criterion is that the absolute value second moment of the magnetic field in the calculated spatial amplitude difference falls to the pre-specified standard, two-dimensionally selecting 8.e and three-dimensionally selecting 5.e-4.