An efficient three-dimensional spectral element metamaterial electromagnetic simulation technique

Abstract

A highly efficient three-dimensional spectral element metasurface electromagnetic simulation technique belongs to the field of electromagnetic simulation. 1) Based on the desired function, find suitable metasurface scattering elements, solve for the scattering parameters, and synthesize the polarizability tensor required for the generalized thin-film transition conditions; 2) Establish a three-dimensional geometric model, selecting the model's material parameters, boundary conditions, and the size and form of the incident beam; 3) Mesh the geometric model using a hexahedral mesh, and mesh the zero-thickness equivalent thin film of the metasurface using a quadrilateral mesh; 4) Read the mesh information, perform preprocessing, set boundary conditions and material parameters, generate the system matrix using the spectral element method, and obtain the spatially discretized matrix equations; 5) Solve for the electric field value and observe the electric field distribution; 6) Determine if the electric field distribution meets expectations. If so, plot the electric field distribution; otherwise, correct the result by adjusting the scattering parameters and resolving the polarizability, repeating steps 4) to 5) until convergence.

Description

Technical Field

[0001] This invention belongs to the field of electromagnetic simulation, and in particular relates to an efficient three-dimensional spectral metasurface electromagnetic simulation technology that combines metasurface simulation methods, spectral element methods, and generalized thin-film transition conditions. Background Technology

[0002] Metasurfaces are subwavelength artificial two-dimensional structures composed of tiny metals or nanostructures. Through special structural design and arrangement, metasurfaces can induce anomalous reflection and transmission of electromagnetic waves within a specific frequency range, providing the possibility of manipulating electromagnetic waves. (CL Holloway, E.F. Kuester, and J.A. Gordon, “An overview of the theory and applications of metasurfaces: The two-dimensional equivalents of metamaterials,” IEEE Antennas and Propagation magazine, vol. 54, no. 2, pp. 10-35, 2012).

[0003] With the deepening research on metasurface structures, thin and compact metasurface structures are constantly being proposed. However, their complex structures and difficulties in selecting unit cells increase fabrication costs. (K. Achouri, C. Caloz, “Design, concepts, and applications of electromagnetic metasurfaces,” Nanophotonics, vol. 7, no. 6, pp. 11095-1116, 2018) Therefore, researching a method to quickly analyze the electromagnetic scattering characteristics of metasurfaces and perform high-fidelity simulations is one of the key technologies for current metasurface design.

[0004] In current commercial software numerical simulations of metasurfaces, for complex-shaped subwavelength units and ultrathin electrically large structures of metasurfaces, very fine meshes are required to divide the model and generate huge degrees of freedom, resulting in a huge amount of computer runtime and resource consumption.

[0005] A novel algorithm combining the three-dimensional spectral element method and generalized thin-plate transition conditions can effectively study the electromagnetic wave propagation characteristics of subwavelength-thickness metasurfaces. It also serves as an important method for integrating the structural and parameter design of metasurfaces based on pre-defined functional designs. This algorithm proposes a novel metasurface meshing method, treating subwavelength-thickness metasurfaces as equivalent to zero-thickness thin plates, significantly reducing mesh density and the number of degrees of freedom, thus improving computational time and space efficiency.

[0006] This numerical simulation technique enables rapid, multiple simulation analyses of complex metasurface structures. While electromagnetic numerical simulation experiments cannot completely replace real experiments, they can provide theoretical support for experimental parameters, thereby significantly reducing the total number of experiments and saving costs. In conclusion, research on computational electromagnetic methods for metasurfaces can aid in the design and optimization process of metasurfaces, possessing high academic research value and promising practical applications. Summary of the Invention

[0007] The purpose of this invention is to address the aforementioned problems in existing technologies by providing a highly efficient three-dimensional spectral element method for electromagnetic simulation of metasurfaces, thereby improving computational efficiency in electromagnetic simulation. By employing the three-dimensional spectral element method and generalized thin-plate transition conditions, a densely meshed metasurface is equated to a sparsely meshed, zero-thickness thin plate, reducing mesh degrees of freedom and increasing numerical solution speed. This improves simulation efficiency in terms of both simulation time and computational resource consumption, providing a theoretical foundation for metasurface design.

[0008] This invention includes the following steps:

[0009] 1) Determine the function to be achieved by the metasurface and find the metasurface scattering unit that meets the requirements; use electromagnetic simulation software to obtain its scattering parameters and input the scattering parameters into the polarizability tensor that needs to be determined for the generalized thin-film transition conditions to synthesize the polarizability tensor.

[0010] 2) Establish a three-dimensional geometric model, select the computational domain size, material parameters, boundary conditions, and incident light beam in the model;

[0011] 3) Perform hexahedral meshing on the geometric model, and convert the thick metasurface into a thin sheet with zero thickness. Use quadrilateral meshes to mesh the thin sheet.

[0012] 4) Read the mesh information, perform preprocessing, set boundary conditions and material parameters, and use the spectral element method to generate the system matrix, thus obtaining the spatially discretized matrix equations.

[0013] 5) Solve the equations to obtain the electric field value;

[0014] 6) Determine if the electric field result achieves the expected function of the metasurface. If so, plot the electric field distribution; otherwise, re-solve the polarizability tensor by adjusting the scattering parameters and return to repeat steps 4) to 5) until the simulation result meets expectations.

[0015] In step 1), a suitable metasurface scattering unit is found. By adjusting the material and size of the scattering unit, a full-band simulation is performed under the periodic boundary to solve the scattering parameters of scattering units with different structures. A suitable band and scattering unit structure are found, the scattering parameters of the scattering unit are obtained, and the electromagnetic polarizability tensor of the homogenized metasurface model is synthesized.

[0016] In step 2), the three-dimensional geometric region is established, and its shape, size, boundary conditions, material parameters, etc., are confirmed; the amplitude, frequency, and propagation direction of the incident beam are determined. The size and placement of the metasurface equivalent sheet are determined.

[0017] In step 4), the system matrix is ​​generated using the spectral element method. The spectral element method is combined with the generalized thin-plate transition condition to generate a three-dimensional system matrix equation. (K. Achouri, AM Salem, and C. Caloz, “General metasurface synthesis based on susceptibility tensors,” IEEE Transactions on Antennas and Propagation, vol. 63, no. 7, pp. 2977-2991, 2015) When setting the boundary conditions of the equivalent thin plate of the metasurface, the generalized thin-plate transition condition is used.

[0018] Compared with the prior art, the advantages of the present invention are as follows:

[0019] 1. This invention employs a generalized thin-film transition condition for the subwavelength thickness metasurface, combining the structure of the scattering unit with the electromagnetic polarizability tensor after the metasurface is homogenized. By calculating electromagnetic methods, the process of electromagnetic waves entering and propagating on the metasurface can be simulated, thereby verifying whether the scattering unit can achieve the expected function.

[0020] 2. This invention employs the spectral element method to discretize and solve the Helmholtz equation. The GLL polynomials used in the spectral element method have spectral accuracy, and the interpolation error decreases exponentially with increasing polynomial order. Therefore, the spectral element method can improve computational accuracy by appropriately increasing the order of the interpolation basis functions.

[0021] 3. This invention employs a combination of the three-dimensional spectral element method and generalized thin-film transition conditions to equate subwavelength-thickness metasurfaces to zero-thickness thin films. This novel metasurface mesh generation method improves the running efficiency and memory usage of simulation software without sacrificing accuracy, thereby enhancing simulation efficiency. Attached Figure Description

[0022] Figure 1 This is a flowchart illustrating an embodiment of the present invention.

[0023] Figure 2 This is a schematic diagram of the scattering unit structure according to an embodiment of the present invention.

[0024] Figure 3 This is a phase diagram showing how the scattering parameters change with the radius of the nanopillar in an embodiment of the present invention.

[0025] Figure 4 This is an amplitude diagram showing the variation of scattering parameters with the radius of the nanopillar in an embodiment of the present invention.

[0026] Figure 5 This is a model diagram and its mesh subdivision diagram in COMSOL of an embodiment of the present invention.

[0027] Figure 6 This is the equivalent model diagram and its mesh subdivision diagram of an embodiment of the present invention.

[0028] Figure 7 This is an electric field distribution diagram according to an embodiment of the present invention. Detailed Implementation

[0029] The following embodiments will further illustrate the technical solution of the present invention with reference to the accompanying drawings.

[0030] like Figure 1 The embodiments of the present invention include the following steps:

[0031] S1: This embodiment of the invention aims to find a metasurface capable of achieving a 14° generalized deflection of the x-component of the electric field. First, it is necessary to find transmission-type scattering units capable of achieving phase changes within the range [0, 2π] with similar amplitudes. Based on the refraction formula of the generalized Snell's law, a suitable wavelength for the incident beam and the structure of the eight scattering units are determined. The eight scattering units are then arranged sequentially. The eight scattering parameters are substituted into the polarizability tensor required to determine the generalized thin-film transition conditions, thus synthesizing the polarizability tensor of the equivalent model for each scattering unit.

[0032] S2: Establish a three-dimensional geometric model, select the computational domain size, material parameters, and boundary conditions in the model, and select the incident beam. In this embodiment of the invention, the computational domain length in the x-direction is 2000 nm, and the computational domain length in the y-direction is 250 nm. Periodic boundary conditions are applied to the front, back, left, and right boundaries; absorbing boundary conditions are applied to the top and bottom; and a generalized thin-film transition condition is applied to the zero-thickness thin film equivalent to the metasurface. The incident electric field is a y-polarized plane wave.

[0033] S3: Perform hexahedral meshing on the geometric model. Then, treat the thick metasurface as a zero-thickness sheet and mesh the sheet using quadrilateral meshes.

[0034] S4: Read the mesh information, perform preprocessing, set boundary conditions and material parameters, generate the system matrix using the spectral element method, and obtain the spatially discretized matrix equations.

[0035] S5: Use the sparse matrix solver Pardiso to solve the equations and obtain the electric field distribution;

[0036] S6: Determine if the electric field result achieves the desired function of the metasurface. If yes, plot the electric field distribution; otherwise, resolve the scattering parameters, synthesize the polarizability tensor, and return to repeat steps 4) to 5) until the simulation result meets expectations.

[0037] S1 shows a schematic diagram of the scattering unit structure for implementing the present invention, as follows: Figure 2 The bottom hexahedron is a SiO2 substrate with a width of 250 nm in both the x and y directions; the top cylinder is a TiO2 nanopillar with a height of 600 nm and a radius of R. The relationship between the nanopillar radius and scattering parameters was measured using COMSOL simulation software.

[0038] When the height of the nanopillar is 600 nm and the incident wavelength is 500 nm, the scattering parameters are... The results of the phase and amplitude variations with the nanopillar radius are as follows: Figure 3 and 4 It can be seen that the phase of the scattering unit can vary from [0, 2π], and the amplitudes are similar. Based on the requirements of the metasurface in this embodiment of the invention, eight scattering units were selected, and their corresponding phase changes and nanopillar radii R are shown in Table 1.

[0039] Table 1

[0040] Phase change 0° 45° 90° 135° 180° 225° 270° 315° Nanopillar radius (nm) 87.6 81.6 75.6 69.3 98.5 98.1 97.3 93.7

[0041] The statement in S1 regarding the generalized Snell's law states that metasurfaces, by introducing abrupt phase changes within the wavelength range, exhibit anomalous reflection and refraction phenomena when this phase changes linearly along the interface. The refraction formula of the generalized Snell's law is:

[0042]

[0043] θ in the formula t and θ i n represents the transmission angle and the incident angle, respectively. t and n i Let represent the refractive index of the transmitted field and the refractive index of the incident field, respectively, and λ0 be the wavelength of the incident beam. It is the phase gradient of phase Φ as position x changes.

[0044] In this embodiment of the invention, eight scattering units were selected. Each scattering unit has a length of 250 nm in the x-direction, and the phase changes of each scattering unit are 0°, 45°, 90°, 135°, 180°, 215°, 270°, and 315°, respectively. Therefore, the total computational domain in the x-direction is 2000 nm, and the total phase change in the x-direction is 2π. The incident angle θ of the incident field... i =0°, the incident region is air, and the refractive index of the medium is n i =1; the transmission region is SiO2, and its medium refractive index n t =1.47; According to the refraction formula (1) of the generalized Snell's law, the theoretical transmission angle θ can be obtained. t ≈14°.

[0045] In S1, the scattering parameters need to be incorporated into the polarizability tensor required for the generalized thin-film transition condition. For a uniaxial, uniaxial anisotropic surface, the relationship between the polarizability components and the scattering parameters is as follows:

[0046]

[0047] in, and For the surface polarizability component of the metasurface, the scattering parameters The value of is defined as the ratio of the electric field when an x-polarized wave is incident from the incident port 1 and a y-polarized wave is transmitted from the transmission port 2. Other scattering parameters are defined similarly.

[0048] The specific steps for establishing the system matrix equation in S4 are as follows:

[0049] For a three-dimensional electromagnetic field, the governing equations are the Helmholtz equations, which are vector equations:

[0050]

[0051] Where Ω is the computational domain for solving, E is the total field within the computational domain, k is the wave number, and μ r Let ε be the relative permeability tensor. r is the relative permittivity tensor.

[0052] Using Galerkin's weighted residual method and the finite element method, the residual weighted integral of element e can be obtained as follows:

[0053]

[0054] in, The test function on the unit, E e The electric field vector represents the unit cell.

[0055] By introducing Green's identity, equation (4) can be written as:

[0056]

[0057] In the formula, Γ e Indicates the enclosure of V e face, Indicates Γ e The unit vector of the outward normal.

[0058] The spectral element method uses hexahedral elements for mesh generation and employs mixed-order Gauss-Lobatto-Legendre (GLL) polynomials to construct interpolation functions on the elements.

[0059] The expression for the Nth-order interpolation function in the one-dimensional case is as follows:

[0060]

[0061] In the formula, ξ∈[-1,1] is the position variable on the reference unit, L N (ξ) and L′ N (ξ) is an Nth-order Legendre polynomial and its first derivative. j The j-th GLL sampling point on the reference element is defined by equation (1-ξ). 2 )L′ N The root of (ξ) = 0.

[0062] By introducing GLL element discretization into the governing equations, the vector basis functions in the reference element (ξ,η,ζ)∈[-1,1]×[-1,1]×[-1,1] can be expressed as:

[0063]

[0064] in, respectively under the reference frame N-order vector basis functions in the direction.

[0065] In the transformation between physical units and reference units, a corresponding covariant mapping needs to be introduced:

[0066]

[0067] in, Let Φ represent the vector basis functions on the reference element, Φ represent the vector basis functions on the physical element, and J be the Jacobian matrix. Its expression is as follows:

[0068]

[0069] Then the electric field E on the physical unit e It can be written as:

[0070]

[0071] Substituting equation (10) into equation (5), we get:

[0072]

[0073] Where F = 3N(N+1) 2 The number of degrees of freedom of the edge of the unit cell at order N.

[0074] Writing (11) in matrix form, we have:

[0075]

[0076] in:

[0077]

[0078] Combining all units and setting the weighted integral of the residuals to zero, we obtain the system of equations:

[0079]

[0080] The above formula can also be written as:

[0081] [S-M+B]{e}={b} (15)

[0082] By solving the above equations, the degree of freedom for each problem can be obtained.

[0083] In S4, using the spectral element method, a generalized thin-film transition condition is applied to the metasurface equivalent thin film. The specific steps are as follows:

[0084] The generalized thin-plate transition condition is used to solve the problem of discontinuous electromagnetic wave distribution at the interface of a metasurface. Assuming the metasurface is placed in the xoy plane and neglecting the normal component of the electromagnetic polarization density, the generalized thin-plate transition condition is:

[0085]

[0086]

[0087] For unidirectional anisotropic materials, the electric polarization density P and magnetic polarization density M can be represented by the average of the polarizability and the electromagnetic field:

[0088]

[0089] in: and These four polarizability components have been determined in S2.

[0090] Substituting equation (18) into equation (17), and expressing the form of each component of the magnetic field in terms of the components of the electric field, we have:

[0091]

[0092] In this context, subscript 1 represents the field on the incident side of the metasurface, and subscript 2 represents the field on the transmission side of the metasurface. The expressions for A1 to A4 are shown below:

[0093]

[0094]

[0095]

[0096]

[0097] According to the process of establishing the system matrix equation in S5, in the spectral element method, for the boundary... And Substituting it into the following categories:

[0098]

[0099] Substituting equation (19) into equation (21) combines the spectral element method with the generalized thin section transition condition.

[0100] The above is the specific implementation process of the present invention.

[0101] The supercell in the method of this invention is modeled in COMSOL as follows: Figure 5 As shown in (a), the meshing of this model in COMSOL is as follows: Figure 5 As shown in (b), it can be seen that due to the small radius of the nanopillars and the close spacing between each nanopillar, a large number of fine meshes are required for subdivision. The equivalent model of the supercell of the method of this invention is as follows: Figure 6 As shown in (a), this method equates a thick metasurface to a thin sheet of zero thickness; therefore, the height of the computational domain of the equivalent model is reduced by the metasurface portion. The meshing of this equivalent model is as follows: Figure 6 As shown in (b), this model employs hexahedral meshing in the spectral element method, reducing the mesh density. Through... Figure 5 and 6 By comparison, it can be found that the method of the present invention treats the metasurface as a sheet of zero thickness, which greatly reduces the mesh density and the degree of freedom.

[0102] The generalized refraction angle deflection achieved by the method of this invention is as follows: Figure 7 As shown. By observing the electric field distribution, it can be found that the x-component of the electric field achieves the transmission angle θ. t =14° deflection.

[0103] Experiments show that this invention develops a highly efficient three-dimensional spectral element metasurface electromagnetic simulation technique. By combining the metasurface's scattering parameters with the polarizability of the generalized thin-plate transition condition, it provides practical assistance for metasurface structural design. This invention solves the computational domain using the spectral element method, reducing mesh density and computational errors. Furthermore, by equating a subwavelength-thickness metasurface to a zero-thickness thin plate, this invention significantly reduces the mesh's degrees of freedom, improves simulation efficiency, and has significant application value for metasurface design and simulation.

Claims

1. An efficient three-dimensional spectral element metasurface electromagnetic simulation technique, characterized by Includes the following steps: 1) Determine the function that the metasurface needs to achieve and find a metasurface scattering unit that meets the requirements; The scattering parameters were obtained by electromagnetic simulation software, and the scattering parameters were then substituted into the polarizability tensor that needs to be determined for the generalized thin sheet transition condition to synthesize the polarizability tensor. The process involves finding a suitable metasurface scattering unit, adjusting the material and size of the scattering unit, performing full-band simulation under periodic boundaries, and solving for the scattering parameters of scattering units with different structures. The process also involves finding the suitable band and scattering unit structure, obtaining the scattering parameters of the scattering unit, and synthesizing the electromagnetic polarizability tensor of the homogenized metasurface model. The process of incorporating the scattering parameters into the polarizability tensor required to determine the generalized thin-film transition conditions, for a uniaxial, unidirectional anisotropic surface, yields the following relationship between the polarizability components and the scattering parameters: in, , , and For the surface polarizability component of the metasurface, the scattering parameters The definition is that it is incident from port 1. Polarized wave is transmitted from transmission port 2 The electric field ratio for polarized waves is similar for other scattering parameters; 2) Establish a three-dimensional geometric model, select the computational domain size, material parameters, boundary conditions, and incident light beam in the model; 3) Perform hexahedral meshing on the geometric model, and convert the thick metasurface into a thin sheet with zero thickness. Use quadrilateral meshes to mesh the thin sheet. 4) Read the mesh information, perform preprocessing, set boundary conditions and material parameters, and use the spectral element method to generate the system matrix, thus obtaining the spatially discretized matrix equations. 5) Solve the equations to obtain the electric field value; 6) Determine if the electric field result achieves the expected function of the metasurface. If so, plot the electric field distribution; otherwise, re-solve the polarizability tensor by adjusting the scattering parameters and return to repeat steps 4) to 5) until the simulation result meets expectations.

2. The efficient three-dimensional spectral element metasurface electromagnetic simulation technique of claim 1, wherein In step 2), a three-dimensional geometric model is established to confirm the shape, size, boundary conditions, and material parameters of the region; the amplitude, frequency, and propagation direction of the incident beam are determined; and the size and placement of the metasurface equivalent sheet are determined.

3. The efficient three-dimensional spectral element metasurface electromagnetic simulation technique of claim 1, wherein In step 4), the system matrix is ​​generated using the spectral element method, which is combined with the generalized thin-plate transition condition to generate a three-dimensional system matrix equation.

4. The efficient three-dimensional spectral element metasurface electromagnetic simulation technique of claim 1, wherein In step 4), the establishment of the system matrix equation involves the following steps: For a three-dimensional electromagnetic field, the governing equations are the Helmholtz equations, which are vector equations: wherein, is the computational domain, is the total field within the computational domain, is the wave number, is the relative permeability tensor, is the relative permittivity tensor; By the Galerkin's weighted residual method and the finite element method, the weighted integral of the residual of the element is obtained as wherein representing a test function on the cell, representing an electric field vector on the cell; By introducing the Green's identity, the equation can be written as: wherein denotes the plane that encloses , denotes the unit vector of the outer normal The spectral element method uses hexahedral elements for mesh generation and constructs interpolation functions on the elements using mixed-order Gauss-Legend-Lobata polynomials. In one dimension The expression of the order interpolation function is as follows: wherein is a position variable on the reference unit, and are Legendre polynomials of order N and their first derivative; is called the k-th GLL sample point on the reference element and is the root of the equation .​ The GLL element is introduced into the control equation for discretization, and in the reference element The vector base function is expressed as: in, , , respectively under the reference frame Direction First-order vector basis functions; In the transformation between physical units and reference units, a corresponding covariant mapping needs to be introduced: wherein denotes a vector basis function on the reference element, denotes a vector basis function on the physical element, is the Jacobian matrix, which is expressed as follows: the electric field on the physical unit is written as: Substituting formula into formula then: wherein, is the number of degrees of freedom of the free edges of the prisms on the time unit. Write In matrix form, we have: in: Combining all units and setting the weighted integral of the residuals to zero, we obtain the system of equations: The above formula can be written as: By solving the above equations, each degree of freedom to be solved is obtained.

5. The efficient three-dimensional spectral element metasurface electromagnetic simulation technique of claim 4, wherein By combining the spectral element method with the generalized thin-plate transition condition, the generalized thin-plate transition condition is used when setting the boundary conditions of the equivalent thin plate of the metasurface. The specific steps are as follows: The generalized sheet transition condition is used to solve the problem of electromagnetic wave distribution discontinuity at the interface of a super surface; assuming that the super surface is placed in a plane, the generalized sheet transition condition is: For a single anisotropic material, the electric polarization density and the magnetic polarization density are expressed in terms of the polarizability and the average of the electromagnetic field: wherein: , , and These four polarizability components have been determined in the preceding steps; Substituting this into the expression for In terms of the components of the electric field, the components of the magnetic field are given by:​ wherein subscript 1 represents the field at the super surface incident side, and subscript 2 represents the field at the super surface transmission side; The expression of the above is shown as follows: In the spectral cell method, for the boundary has , where is introduced into which has: Substituting the expression for into the expression for i.e. the spectral element method is combined with the generalized thin strip transition condition.

Images