A two-dimensional finite element mesh non-conformal h-type adaptive refinement method

Abstract

The application belongs to the field of finite element grid encryption, and specifically provides a two-dimensional finite element grid non-conformal h-type adaptive encryption method to overcome the shortcoming that the traditional finite element method cannot perform non-conformal encryption. The application adopts a normalized edge-based function to process a non-conformal area in a triangular grid, so that the finite element method can solve unknown electric field intensity by using a non-conformal grid, and then calculate a numerical solution of a radar scattering cross section by using the obtained electric field intensity. In combination with an error estimation method, the application can freely encrypt a grid with a large error contribution, and can efficiently reduce the error between the numerical solution and the true solution. The application effectively avoids the increase of distorted grids in the encryption process without introducing additional unknown quantities, and ensures the grid quality. Meanwhile, the application realizes grid adaptive encryption, and can reduce unknown quantities and improve the solving precision compared with uniform encryption.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of finite element mesh refinement, and particularly provides a two-dimensional finite element mesh non-conformal h-type adaptive refinement method based on normalized edge basis functions. BACKGROUND

[0002] In the field of computational electromagnetics, the finite element method is often used to solve electromagnetic fields under specific boundary value problems. In general, in order to improve the accuracy of the calculation results, more refined meshes are needed to discretize the electromagnetic field. However, the use of too refined meshes will significantly increase the calculation degrees of freedom, thereby increasing the calculation time and memory consumption. In order to solve this problem, researchers have proposed mesh adaptive refinement technology, which selectively refines meshes with larger error contributions. However, the traditional finite element method requires the mesh to be conformal, and the nodes at both ends of the common edge of adjacent elements are required to be the same in the refinement process. This will cause the model to have more and more distorted triangular elements as the mesh is continuously refined, which reduces the error elimination efficiency and even increases the error. The use of domain decomposition method to realize non-conformal mesh requires to impose continuity conditions on the interface, which introduces additional unknown quantities and increases the consumption of computing resources. SUMMARY

[0003] The purpose of the present application is to provide a two-dimensional finite element mesh non-conformal h-type adaptive refinement method based on normalized edge basis functions, which overcomes the shortcomings of the traditional finite element method that cannot perform non-conformal refinement. Moreover, the present application effectively avoids the increase of distorted meshes in the refinement process without introducing additional unknown quantities, thereby ensuring the quality of the mesh. In addition, the present application can realize mesh adaptive refinement, reduce unknown quantities and improve the solving accuracy in combination with the error estimation method. Furthermore, the present application is easy to implement, effectively saves computing resources and improves the efficiency of the adaptive refinement method.

[0004] To achieve the above object, the technical scheme adopted by the present application is as follows:

[0005] A two-dimensional finite element mesh non-conformal h-type adaptive refinement method, characterized in that it comprises the following steps:

[0006] Step 1. Discretize the model using a triangular mesh to obtain the node information and element information of all triangular elements. The node information includes the node number and node coordinates, and the element information includes the element number and the node numbers of the three vertices of the element.

[0007] Step 2. Establish a matrix equation for the electric field intensity to be solved based on the mesh information, and solve to obtain the electric field intensity coefficient to calculate the numerical solution of the radar cross section.

[0008] Step 3. The error estimation method is used to calculate the error between the numerical solution and the true solution of the radar cross section, and when the error is less than the preset threshold, the iteration is stopped, and the numerical solution of the radar cross section calculated in step 2 is output; otherwise, step 4 is executed;

[0009] Step 4. Select N e triangular elements with the largest absolute value of error contribution, and perform mesh refinement on each triangular element;

[0010] Step 5. Take the refined mesh information as input and return to step 2 for iteration.

[0011] Further, in step 2, the matrix equation for the electric field intensity to be solved is specifically:

[0012] [K]·{E Ω}={b},

[0013] Where K is a sparse symmetric matrix, E Ω is the electric field intensity coefficient vector to be solved, b is the right end excitation vector, [·] represents a square matrix, and {·} represents a vector.

[0014] The element in the nth row and mth column of K nm is K

[0015]

[0016] Where represents the curl of the function, N m and N n represent the normalized edge basis functions of edges m and n, μ r and ε r represent the relative permeability and relative permittivity of the triangular element, k0 represents the wave number in vacuum, Ω0 represents the solution domain, j represents the imaginary unit, is the outward normal vector on the boundary, represents the absorbing boundary of the solution domain.

[0017] The nth element in b n is b

[0018]

[0019] Where U inc represents the excitation term.

[0020] The mth element in E Ω is E m , and E m represents the electric field intensity coefficient of the mth edge edge basis function.

[0021] Furthermore, in step 2, the numerical solution of the radar cross section is expressed as:

[0022]

[0023] Where φ represents the receiving direction in polar coordinates, σ(φ) represents the radar cross section in the receiving direction φ, and J x 、J y They represent the equivalent current J in the Cartesian coordinate system respectively Γ The x and y components of M z Denotes the equivalent magnetic current M in the Cartesian coordinate system Γ The z component of ; η represents the wave impedance in vacuum, and x' and y' represent the coordinates of the source point;

[0024] Equivalent current J Γ and the equivalent magnetic current M Γ Expressed as:

[0025]

[0026]

[0027]

[0028]

[0029] in, represents the outward normal vector of the cylindrical surface, ω is the angular frequency, and u is the magnetic permeability.

[0030] Furthermore, in step 3, an error estimation method based on an adjoint operator is adopted, where the QOI is constructed as:

[0031]

[0032] The error estimate between the numerical solution and the true solution of QOI is used as the error between the numerical solution and the true solution of radar cross section, which is expressed as:

[0033]

[0034] Among them, R E represents the error estimate between the QOI numerical solution and the true solution, R e Indicates the error contribution of the triangle element with sequence number e, U e is the total number of triangle elements.

[0035] Furthermore, in step 4, N e The triangle elements satisfy the conditions:

[0036]

[0037] wherein C1 is an empirical coefficient between 0 and 1, R e denotes the error contribution of the triangular element with index e, U e is the total number of triangular elements.

[0038] Based on the above technical solution, the application has the following beneficial effects:

[0039] The application provides a two-dimensional finite element mesh non-conformal h-type self-adaptive refinement method based on normalized edge basis functions, which adopts normalized edge basis functions to process non-conformal regions in triangular meshes, so that the finite element method can solve unknown electric field intensity by using non-conformal meshes. The method provided by the application allows the introduction of non-conformal regions in the refinement process, and can divide any triangular element into four smaller triangular elements similar to the original triangular element without additional processing of adjacent triangular elements. Moreover, the application can well ensure mesh quality in the refinement process and avoid the appearance of abnormal triangular elements; in combination with an error estimation method, the application can freely refine meshes with large error contributions, so that the application can efficiently reduce the error between a numerical solution and a true solution. Meanwhile, the normalized edge basis function is simple to operate in processing non-conformal regions, does not need to impose continuity, and avoids the introduction of additional unknown quantities. Finally, the application has the advantages of simple operation, guaranteed mesh quality in the refinement process, avoidance of the introduction of additional unknown quantities, combination with an error estimation method, realization of mesh self-adaptive non-conformal refinement, improved solving precision, saved computing resources, and improved self-adaptive refinement efficiency. BRIEF DESCRIPTION OF DRAWINGS

[0040] Figure 1 FIG. 1 is a schematic diagram of four-refinement of a triangular element in an embodiment of the application.

[0041] Figure 2 FIG. 2 is a schematic diagram of a non-conformal structure in an embodiment of the application.

[0042] Figure 3 FIG. 3 is an edge basis function vector diagram on a boundary edge of a non-conformal region in an embodiment of the application.

[0043] Figure 4 FIG. 4 is a mesh diagram of a non-conformal refinement of a uniform medium cylindrical mesh in an embodiment of the application.

[0044] Figure 5 FIG. 5 is a relative error reduction diagram of a monostatic radar scattering cross section of an infinite uniform medium cylinder in an embodiment of the application. DETAILED DESCRIPTION

[0045] To make the purpose, technical solutions and beneficial effects of the application clearer, the application is further described in detail below with reference to the drawings and embodiments.

[0046] The embodiment provides a two-dimensional finite element mesh non-conformal h-type self-adaptive encryption method based on normalized edge-based functions, and takes an infinite long uniform medium cylinder XOY section as an example to adaptively encrypt the infinite long uniform medium cylinder, adopts an absorbing boundary condition, and calculates single station RCS (Radar Cross Section) of the infinite long uniform medium cylinder.

[0047] Specifically, the method comprises the following steps:

[0048] Step 1, triangular mesh discrete models are adopted to obtain node information and element information of all triangular elements, wherein the node information comprises a node serial number and node coordinates, and the element information comprises an element serial number and node serial numbers of three vertexes of the element;

[0049] Step 2, mesh information is read, an electric field intensity to be solved is taken as an unknown quantity, a vector-based function is adopted to expand the unknown electric field intensity, a matrix equation is established by testing a wave equation and a boundary condition equation according to a Galerkin method, the matrix equation is solved, and the unknown electric field intensity is obtained, and then a numerical solution of the radar scattering cross section is calculated through the obtained electric field intensity;

[0050] The specific process is as follows:

[0051] Firstly, a boundary value problem is given, the wave equation is formula (1), and the boundary condition takes an absorbing boundary condition as an example, as shown in formula (2);

[0052]

[0053]

[0054] Wherein, is a Hamiltonian operator, for any vector function G, represents a curl of the function G, represents a divergence of the function G, for any scalar function g, represents a divergence of the function g; μ r and ε r respectively represent relative permeability and relative permittivity of a discrete element, k0 represents a wave number in vacuum, E represents the electric field intensity to be solved, and Ω0 represents a solution domain;

[0055] j represents an imaginary unit, is a normal vector outward of the boundary, represents an absorbing boundary of the solution domain; U inc represents an excitation term, E inc represents an incident electric field;

[0056] In the embodiment, the incident electric field is the incident magnetic field is H inc =ez e -jkx For example, in order to facilitate the derivation of the following formula (16), wherein, e y is the unit vector along the positive direction of the y-axis in the Cartesian rectangular coordinate system, e z is the unit vector along the positive direction of the z-axis in the Cartesian rectangular coordinate system, η0 is the numerical value of the wave impedance in the vacuum, and x represents the x-axis coordinate in the Cartesian rectangular coordinate system.

[0057] Next, the electric field intensity is expanded by using the basis function. For any triangular unit, the normalized vector edge basis function form used on the edge m is as formula (3). The edge basis function of this form ensures that the tangent vector on the unit boundary is 1, which is one of the key points of the non-conformal encryption method of the application.

[0058]

[0059] Wherein, N m represents the normalized vector edge basis function of the edge m, and m represents the serial number of the edge; p and q represent two nodes of the edge m, represents the gradient of the scalar function L q , represents the gradient of the scalar function L p , l m is the length of the edge m; L p represents the shape function corresponding to the node p in the triangular unit; L p is as formula (4):

[0060]

[0061] Wherein, (x p , y p ), (x q , y q ), (x r , y r ) are the coordinates of the nodes p, q, r in the triangular unit, (x, y) represents the coordinates of any point in the triangular unit, and S represents the area of the triangular unit.

[0062] Similarly, L q represents the shape function corresponding to the node q in the triangular unit:

[0063] Next, the Galerkin method is used to select the same weighted function as the basis function. The inner product of the weighted function on both sides of formula (1) and formula (2) is obtained, that is, the weak form of the original edge value problem, as formula (5) and formula (6).

[0064]

[0065]

[0066] Among them, N n represents the normalized vector edge basis function of edge n, which is also the nth weighting function;

[0067] Derivation can be obtained as formula (7);

[0068]

[0069] When solving the above equation using the finite element method, the electric field E can be approximately expanded into a linear combination of basis functions in the finite element space: Among them, E m is the electric field intensity coefficient of the edge basis function of the mth edge, and M is the total number of edge basis functions;

[0070] We can get formula (8):

[0071]

[0072] Among them, K nm The specific expression of is shown in formula (9), b n The specific expression of is shown in formula (10)

[0073]

[0074] Equation (8) can be written more compactly as a matrix equation, as shown in Equation (11):

[0075] [K]·{E Ω}={b} (11)

[0076] Among them, K is a sparse symmetric matrix, E Ω is the electric field strength coefficient vector to be determined, b is the right-hand excitation vector; [·] represents a square matrix, {·} represents a vector; the element in the nth row and mth column of K is K nm , the specific expression is shown in formula (9), E Ω The mth element in is E m , the nth element in b is b n , the specific expression is shown in formula (10);

[0077] For the edge basis functions at the boundary of the non-conformal region, their electric field strength coefficients need to be the same to ensure the continuity of the tangential electric field at the boundary of the unit. In this embodiment, this can be achieved in the following two ways:

[0078] 1) Let edge m be the boundary edge of the non-conformal region, as Figure 2 Taking the boundary edge of the non-conformal region shown as an example, the edge basis functions on the boundary edge of the non-conformal region are The corresponding electric field strength coefficient is In the process of solving, let the electric field strength coefficient vector to be solved be

[0079] 2) Let edge m be the boundary edge of the non-conformal region, as shown in Figure 2 Taking the boundary edge of the non-conformal region shown as an example, the edge basis functions on the boundary edge of the non-conformal region are The corresponding electric field strength coefficient is In the solution process, the elements in the sparse symmetric matrix K are merged, the rows m1, m2, and m3 are merged into the mth row, and the columns m1, m2, and m3 are merged into the mth column. Let

[0080] Solve the above matrix equation (the solution process is well known in the art and will not be repeated here), and get the electric field intensity coefficient vector E Ω , that is, the electric field intensity coefficient of the edge basis function of each edge is obtained, and then the numerical solution E of the unknown electric field intensity is calculated by the electric field intensity expansion formula (12) h ;

[0081]

[0082] The numerical solution of the electric field intensity can be obtained by calculating the numerical solution of the magnetic field intensity through formula (13);

[0083]

[0084] Where ω is the angular frequency and u is the magnetic permeability;

[0085] Γ represents the cross section of the cylindrical surface, J Γ Represents the equivalent current on the cylindrical surface, M Γ Represents the equivalent magnetic flux on the cylindrical surface; through equations (14) and (15), J can be calculated respectively Γ and M Γ :

[0086]

[0087]

[0088] in, Represents the outward normal vector of the cylindrical surface;

[0089] When the amplitude of the incident magnetic field intensity is 1 Tesla, according to Huygens principle, we can calculate the numerical solution of the radar cross section through equation (16):

[0090]

[0091] where, φ represents the receiving direction in polar coordinates, σ(φ) represents the radar scattering cross section in the receiving direction φ, J x , J y respectively represent the x component and y component of the equivalent current J Γ in the Cartesian rectangular coordinate system, M z represents the z component of the equivalent magnetic current M Γ in the Cartesian rectangular coordinate system; η represents the wave impedance in vacuum, and x' and y' represent the source point coordinates;

[0092] Step 3. An error estimation method is used, which can estimate the error between the numerical solution and the true solution;

[0093] The error of the QOI (quantity of interest) can be estimated based on the accompanying error estimation method (reference: Adjoint-Based Accelerated Adaptive Refinement in Frequency Domain 3-D Finite Element Method Scattering Problems). The present application takes the error estimation method based on the accompanying as an example, constructs the QOI related to the radar scattering cross section, estimates the calculation error of the QOI, and reduces the QOI error through adaptive encryption, so as to realize the effect of reducing the calculation error of the radar scattering cross section;

[0094] The form of the QOI constructed by the present application is shown in formula (17)

[0095]

[0096] It can be known from formula (16) and formula (17) that the smaller the error of the QOI is, the smaller the error of the radar scattering cross section is;

[0097] R E represents the estimated value of the error between the numerical solution and the true solution of the QOI, R e is the error contribution of the triangular element with the serial number e, and the relationship between the two is shown in formula (18);

[0098]

[0099] where, U e is the total number of triangular elements, and the error estimation of the QOI is obtained by calculating R e of each element and accumulating;

[0100] If the error estimation of the QOI is less than a threshold value, the iteration is stopped, and the numerical solution of the radar scattering cross section calculated in step 2 is output; otherwise, step 4 is executed, and the grid non-conformal encryption is performed;

[0101] Step 4. Ne The cumulative error contribution modulus of each unit satisfies formula (19):

[0102]

[0103] Wherein, C1 is an empirical coefficient between 0 and 1;

[0104] And the N e units are meshed, and the meshing method is shown in Figure 1 , four-partition meshing is adopted, the midpoint of the three edges of the unit i to be meshed is connected, the i-th unit is divided into four smaller units i1, i2, i3, i4, and adjacent units are not meshed; since adjacent units are not meshed, the adjacent units (assuming the number is j) and i1, i2 have overlapping edges, as shown in Figure 2 , j has one same node A and another different node with i1, j also has one same node C and another different node with i2, the edge AC partially overlaps with the edge AB and the edge BC, at this time the mesh is non-conformal; the electric field intensity coefficients of the edge basis functions corresponding to the three edges AB, AC and BC are equal, which can ensure the tangential continuity of the interpolation electric field on the edge AC, and the normalized edge basis functions of the edges AB, AC and BC are as shown in Figure 3 ; in this embodiment, the mesh of the uniform medium cylindrical mesh after non-conformal meshing once is as shown in Figure 4 ;

[0105] Step 5. Take the mesh information after meshing as input, and return to execute step 2 for iteration.

[0106] In this embodiment, the relative error reduction diagram of the single-station radar scattering cross section of the infinite uniform medium cylinder is as shown in Figure 5 , the relative error refers to the relative error of the numerical solution and the analytical solution. Among them, the two curves show the unknown quantity increased and the relative error of the radar scattering cross section reduced by each iteration of uniform meshing and non-conformal adaptive meshing, the first data point of the two curves represents the initial relative error of the radar scattering cross section, and each data point after that represents the relative error of the radar scattering cross section calculated after each meshing. As can be seen from the figure, compared with uniform meshing, non-conformal adaptive meshing can reduce more error while increasing fewer unknown quantities.

[0107] The above is only a specific embodiment of the present application, any feature disclosed in the specification can be replaced by other equivalent or similar purpose alternative features unless specifically described; all features disclosed, or steps in all methods or processes, except for mutually exclusive features and / or steps, can be combined in any way.

Claims

1. A two-dimensional finite element mesh non-conforming h-type adaptive refinement method, characterized in that, Comprising the following steps: Step 1. Discrete model is adopted by using triangular mesh, node information and element information of all triangular elements are obtained; wherein, the node information contains node serial number and node coordinates, and the element information contains element serial number and node serial number of three vertexes of the element; Step 2. Matrix equation for the electric field intensity to be solved is established according to the mesh information, the electric field intensity coefficient is solved, and the numerical solution of the radar cross section is calculated; The matrix equation for the electric field intensity to be solved is specifically: [K] · {E Ω} = {b}, wherein K is a sparse symmetric matrix, E Ω is the electric field intensity coefficient vector to be solved, b is the right end excitation vector, [·] represents a square matrix, and {·} represents a vector. The element in the nth row and the mth column of K is K nm is expressed as: where ∇× denotes the curl of a function, N m , N n denote the normalized vector edge basis functions of edges m, n, μ r and ε r denote the relative permeability and relative permittivity of the triangular element, k0denotes the wave number in vacuum, Ω0denotes the solution domain, j denotes the imaginary unit, is the outward normal vector on the boundary, denotes the absorbing boundary of the solution domain; the nth element of b n is represented as: wherein U inc represents an excitation term; E Ω the mth element is E m , E m denotes the edge-based function of the electric field intensity of the mth edge; The numerical solution of the radar cross section is expressed as: where φ denotes the receiving direction in polar coordinates, σ(φ) denotes the radar cross section at the receiving direction φ, J x , J y denote the x and y components of the equivalent current J Γ in the Cartesian rectangular coordinate system, respectively, M z denotes the z component of the equivalent magnetic current M Γ in the Cartesian rectangular coordinate system; η denotes the wave impedance in vacuum, and x' and y' denote the source point coordinates. Equivalent current J Γ Equivalent magnetic current M Γ is expressed as: wherein represents the outward normal vector of the cylindrical surface, ω is the angular frequency, and u is the magnetic permeability. Step 3. The error estimation value of the numerical solution of the radar cross section and the true solution is calculated by using the error estimation method, when the error estimation value is less than the preset threshold value, the iteration is stopped, and the numerical solution of the radar cross section calculated in step 2 is output; otherwise, step 4 is executed; Step 4. Select N triangle elements with the largest absolute values of error contributions and refine each triangle element. e Step 4. Select N triangle elements with the largest absolute values of error contributions and refine each triangle element. Step 5. The encrypted mesh information is taken as input, and step 2 is executed iteratively.

2. The two-dimensional finite element mesh non-conforming h-type adaptive refinement method of claim 1, wherein, In step 3, the error estimation method based on the adjoint operator is adopted, wherein the QOI is constructed as: The error estimation value of the numerical solution of the QOI and the true solution is taken as the error of the numerical solution of the radar cross section and the true solution, and is expressed as: where R E represents the error estimate of the QOI numerical solution from the true solution, R e represents the error contribution of the triangular element with index e, U e is the total number of triangular elements.

3. The two-dimensional finite element mesh non-conforming h-type adaptive refinement method of claim 1, wherein, In step 4, N e triangular cells satisfy the condition: where C1 is an empirical coefficient between 0 and 1, R e denotes the error contribution of the triangle cell with index e, U e is the total number of triangle cells.

Images