Numerical algorithm for efficiently simulating elastic metasurface in vertical layered elastic model

Abstract

The invention provides a numerical algorithm for efficiently simulating an elastic metasurface in a vertical layered elastic model, which is particularly suitable for rapid and accurate simulation of the elastic metasurface. Figure 1 demonstrates the core of the problem in which metasurface thickness is negligible compared to wavelength, in horizontal orientation and in contact with the interface of Bloch Periodic Boundary Conditions (BPBC). According to the method, a reflection type elastic metasurface is introduced to control the phase leading edge of a reflected wave, a semi-analytical solver based on a numerical mode matching (NMM) method is developed, and a three-dimensional elastic wave problem is simplified into a two-dimensional characteristic value problem in the horizontal direction and an analytical expression in the vertical direction. Compared with a traditional three-dimensional numerical method, such as a finite element method (FEM), the NMM method has the advantages that the calculation cost is obviously reduced, and meanwhile high precision is kept.

Description

Technical Field

[0001] The present invention belongs to the technical field of elastic wave simulation, and specifically relates to a numerical algorithm for efficiently simulating elastic supersurfaces in a vertically layered elastic model, and is particularly suitable for accurately simulating elastic supersurfaces. Background Art

[0002] In recent years, metasurfaces have become an important research direction in the field of elastic wave engineering. Due to their unique physical properties, metasurfaces can flexibly and effectively control the phase, polarization state, and propagation mode of elastic waves. Therefore, fast and accurate design and simulation of metasurfaces are key to the successful development of related devices. To achieve this goal, there is an urgent need for a fast and accurate simulation tool to handle metasurface simulation problems.

[0003] The thickness of some metasurfaces is much smaller than their operating wavelength, so they can be considered as extremely thin surfaces attached to the interface of certain layers. In this case, the properties of these metasurfaces can be characterized by simulating the propagation and action of elastic waves in multilayer structures.

[0004] With the development of engineering technology, the propagation behavior of elastic waves in multilayer structures has become a research hotspot in the field of elastic wave engineering. Traditional numerical simulation methods, such as the finite element method (FEM) and the finite difference time domain method (FDTD), although they have certain universality, face huge challenges in dealing with complex multilayer structures. These methods usually require the discretization of the entire three-dimensional space, resulting in an exponential increase in the amount of calculation with the increase in structural complexity, especially when simulating thin-layer structures or high-frequency elastic wave propagation, the consumption of computing resources is more significant. Therefore, developing a simulation method that can maintain high accuracy and significantly reduce computing costs is of great significance for research in the field of elastic wave engineering. Summary of the invention

[0005] The present invention aims to solve the problems of low simulation efficiency and high consumption of computing resources in the vertical layered elastic model in the prior art, and proposes an efficient numerical algorithm based on the numerical pattern matching (NMM) method, which is specially used for simulating vertical layered elastic models and is also applicable to elastic metasurfaces. The algorithm decomposes the three-dimensional elastic wave problem into a two-dimensional eigenvalue problem in the horizontal direction and an analytical expression in the vertical direction, thereby greatly reducing the complexity of the calculation and improving the calculation efficiency.

[0006] Specifically, the core technical solution of the present invention includes the following aspects:

[0007] First, accurate processing of elastic metasurfaces: The present invention introduces elastic thin layer boundary conditions to simulate elastic metasurfaces whose thickness is much smaller than the wavelength. By analyzing the influence of changes in thin layer parameters on the phase of the reflected wave, the present invention can achieve accurate control of the phase front of the reflected wave, thereby meeting specific simulation requirements. This design based on changes in thin layer parameters allows the propagation direction of the wave to undergo an expected deflection, achieving an ideal reflection or transmission effect. The specific process is as follows:

[0008] For two half-spaces connected by an elastic thin layer boundary condition, when the P wave is incident vertically, the reflection coefficient can be expressed analytically as:

[0009]

[0010] Where Z1=ρ1v p,1 and Z2=ρ2v p,2 They represent the P wave impedance of regions 1 and 2 respectively, and ω is the angular frequency. At this time, the reflection coefficient R PP The amplitude and phase can be expressed as:

[0011] R PP =|R PP |[cos(θ1-θ2)+jsin(θ1-θ2)]

[0012] in, and Therefore, the incident wave is reflected from the surface with a phase shift of Therefore, the design of the elastic metasurface is based on the thin layer parameter Z = diag ([Z t ,Z t ,Z n ]) can deflect the wave propagation direction to achieve the reflection or refraction effect of the target. The reflection phenomenon caused by the phase change includes abnormal reflection, and its reflection direction follows the generalized Snell's law:

[0013]

[0014]

[0015] in, and are the reflection phase gradients of the surface in the x and y directions, respectively. In order to achieve the desired reflection angle (θ r ,φ r ), the analytical expression of ψ(x,y) can be calculated according to the above formula. The thin layer parameters include elastic modulus and thin layer thickness, and these parameters vary in spatial coordinates.

[0016]

[0017] where ψ(x,y) is the desired phase change.

[0018] The second aspect is an efficient simulation method: the present invention introduces the high-order spectral element method (SEM) and the numerical pattern matching (NMM) method. The NMM method simplifies the three-dimensional problem into a two-dimensional eigenvalue problem and an analytical expression by decomposing the solution process. This method not only reduces the number of degrees of freedom in the solution process, but also effectively improves the computational efficiency of the simulation by calculating the eigenmode using SEM. Compared with traditional three-dimensional numerical methods, the method of the present invention significantly reduces the requirements for memory and computing time, while ensuring high precision, greatly improving the overall efficiency of the simulation. The algorithm specifically includes:

[0019] Each layer of the region needs to obtain the characteristic mode pair through the spectral element method where γ z is the propagation constant, is the corresponding characteristic function:

[0020]

[0021] Among them, Φ(x,y) is the row vector composed of node basis functions, u j is a column vector containing the values ​​of the characteristic function at discrete points. The stress field is discretized as:

[0022]

[0023] Using the basis function test on the stress expression, we can get the column vector expression of the stress at discrete points:

[0024] T ij =M0- 1 M ijk u k

[0025] Where, M0=<Φ,Φ>, Where <,> represents the inner product.

[0026] The incident displacement field and the stress normal component of the incident field in region 1 can be expressed as:

[0027]

[0028] The reflected displacement field and the stress normal component of the reflected field in region 1 can be expressed as:

[0029]

[0030] The transmitted displacement field and the stress normal component of the transmitted field in region 2 can be expressed as:

[0031]

[0032] The total field in region 1 can be expressed as the incident field plus the reflected field. Therefore, the expressions for the normal component of the displacement field and stress in region 1 are as follows:

[0033]

[0034] The total field in region 2 can be expressed as the transmitted field. Therefore, the expressions for the displacement field and the normal component of the stress in region 2 are as follows:

[0035]

[0036] in, and is a matrix consisting of a set of characteristic functions of region N, Γ N,z It is a diagonal matrix of propagation constants corresponding to the characteristic function. The superscripts “+” and “-” indicate the propagation direction. Define the unknown quantity to be solved: reflection matrix R 1,2 , the transmission matrix T 1,2 , excitation coefficient a1.

[0037] According to the given incident field expression u0(x,y,z), the excitation coefficient a1 corresponding to the incident field from region 1 is calculated. At the incident surface z=z0, we get:

[0038]

[0039] After testing this expression with each characteristic function in turn, the excitation coefficients are obtained by solving the following system of linear equations:

[0040]

[0041] The third aspect is the processing of boundary conditions: at the interface between different regions, sliding boundary conditions are used to simulate the characteristics of the elastic thin layer. The introduction of this boundary condition allows the reflection and transmission behavior of the elastic wave to be accurately described when the thickness of the thin layer is much smaller than the wavelength. In addition, the accuracy and reliability of the simulation results can be further improved by numerically solving the reflection matrix and the transmission matrix. The specific process is as follows: By applying boundary conditions at the interface between the two layers, the reflection matrix and the transmission matrix are determined. Since the thickness of the elastic thin layer is much smaller than the wavelength, we introduce sliding boundary conditions to simulate the elastic thin layer as follows:

[0042]

[0043] Among them, K is the thin layer material parameter, h is the thin layer thickness, and the elastic thin layer parameter Z = K / h. At the interface z = z1, we get the relationship:

[0044]

[0045] Using the characteristic function of region 2 as the test function for each relation we obtain the linear system:

[0046]

[0047] The expression of the coefficient matrix is ​​as follows:

[0048]

[0049] By solving the linear equations, the reflection matrix and transmission matrix can be obtained.

[0050] Fourthly, a wide range of applications: the method of the present invention is not only applicable to elastic wave simulation in multi-layer structures, but can also be extended to other similar wave phenomena, such as the propagation of sound waves, etc. The versatility of this method makes it an effective tool for solving complex multi-layer structure simulation problems.

[0051] Compared with the prior art, the beneficial results of the present invention are:

[0052] Compared with existing numerical simulation methods, the present invention shows significant advantages in many aspects. First, by introducing the NMM method, the present invention successfully reduces the computational complexity of three-dimensional problems to two dimensions, which significantly reduces the computing time and resource requirements. Secondly, when dealing with the simulation of elastic metasurfaces, the present invention can more accurately control the phase change and avoid the numerical errors in traditional methods. In addition, the method of the present invention shows excellent adaptability and versatility when simulating complex vertical layered structures, and is suitable for a variety of engineering application scenarios. DETAILED DESCRIPTION

[0053] The present invention will explain in detail the application of numerical pattern matching (NMM) method in simulating vertical layered elastic model through multiple embodiments. Each embodiment covers the whole process from model construction, meshing to final simulation result analysis, ensuring that users can clearly understand and correctly apply the technical solution of the present invention.

[0054] like Figure 1 As shown, the present invention proposes a numerical algorithm for efficiently simulating elastic metasurfaces in vertical layered elastic models. The whole method comprises the following steps:

[0055] Based on the theoretical analysis of the phase change caused by the change of thin layer parameters, an elastic metasurface that meets specific simulation requirements is designed. By adjusting the thin layer parameters so that they change with the spatial coordinates, a specific refraction or reflection effect can be achieved.

[0056] In CUBIT software, the cross section of each layer is accurately modeled to generate a two-dimensional grid suitable for calculation. Second-order quadrilateral meshing is used to ensure calculation accuracy and efficiency. Finally, the mesh data is exported for subsequent calculations.

[0057] Based on the elastic wave equation, the model is solved by a numerical pattern matching method. This method significantly improves the simulation efficiency by decomposing the three-dimensional elastic wave problem into a two-dimensional eigenvalue problem in the horizontal direction and an analytical expression in the vertical direction.

[0058] First, the spectral element method (SEM) is used to analyze the waveguide problem of each layer model to obtain the eigenmodes of each layer. By expressing the incident wavefield as a linear combination of these eigenmodes, the propagation characteristics of the wavefield can be described more accurately.

[0059] Secondly, the excitation coefficients are calculated based on the given incident field representation. This step is crucial and determines the accuracy of the simulation results.

[0060] Boundary conditions are applied at the interface of the elastic metasurface, and the behavior of the wave field on the interface is described by solving the reflection matrix and the transmission matrix.

[0061] The total field expression of each layer is calculated using the excitation coefficients and matrices obtained in the previous steps. In this process, the superposition and propagation of the wave fields of each layer are strictly controlled to ensure the simulation accuracy.

[0062] By observing and analyzing the simulation results, it is verified whether the designed metasurface has the expected refractive or reflective properties. The simulation results show that the method of the present invention not only significantly improves the simulation efficiency, but also maintains high accuracy in complex layered structures. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] The accompanying drawings are included to provide a further understanding of the embodiments and are incorporated into and constitute a part of this specification. The accompanying drawings illustrate the embodiments and are used together with the description to explain the principles of the present invention. It will be easy to recognize other embodiments and many expected advantages of the embodiments because they become better understood by reference to the following detailed description. The elements of the drawings are not necessarily to scale with each other. The same reference numerals refer to corresponding similar parts.

[0064] Figure 1 A schematic flow chart of a numerical algorithm for efficiently simulating an elastic metasurface in a vertically layered elastic model according to an embodiment of the present invention;

[0065] Figure 2 A schematic diagram of a double-layer non-uniform model structure in a numerical algorithm for efficiently simulating an elastic hypersurface in a vertically layered elastic model according to an embodiment of the present invention;

[0066] Figure 3 Comparison of displacement field components along the x-axis at y=0 and z=0.5 using the NMM method, a hybrid spectral element method, and the commercial software COMSOL for a double-layer non-uniform half-space structure of an embodiment of the present invention;

[0067] Figure 4 This is a comparison diagram of phase shifts calculated using the NMM method and the analytical method when the impedance parameter of the elastic thin layer connecting the two layers changes in the range of [-300, 300] GPa in a double-layer uniform half-space of an embodiment of the present invention;

[0068] Figure 5 Schematic diagram of the phase shift and impedance of the gradient metasurface of an embodiment of the present invention as a function of coordinate x. The blue dashed line shows the schematic diagram of the designed phase shift as a function of coordinate x. The red line shows the schematic diagram of the impedance as a function of coordinate x calculated using the corresponding equation in the present invention;

[0069] Figure 6 Comparison of displacement field components along the z-axis at x=0 and y=0 using the NMM method, a hybrid spectral element method, and the commercial software COMSOL for a structure containing a gradient metasurface according to an embodiment of the present invention;

[0070] Figure 7 FIG. 1 is a schematic diagram of the divergence of the reflected displacement field and the transmitted displacement field in the XOZ plane view at y=0 for a structure containing a gradient metasurface according to an embodiment of the present invention.

[0072] Example 1: Two non-uniform half-spaces

[0073] In this implementation, we studied Figure 2A double-layer inhomogeneous non-spatial model is shown, which is a soli-pile system to verify the performance of the NMM method. The calculation area contains two solid inhomogeneous regions. Region 1 is a cladding block with a cylindrical core. The cladding block geometry is [-0.5,0.5]m by [-0.5,0.5]m by [-4,0]m, the P-wave velocity is 2378.7m / s, the S-wave velocity is 1373.3m / s, and the density is 2400kg / m^3. The radius of the cylindrical core is 0.2m, the center is located at (-0.1,0,-2), the P-wave velocity is 220.51m / s, the S-wave velocity is 90.02m / s, and the density is 1939kg / m^3. Region 2 is also a cladding block with a cylindrical core. The cladding block geometry is [-0.5,0.5]m by [-0.5,0.5]m by [0,4]m, the P wave velocity is 1802.2m / s, the S wave velocity is 1103.6m / s, and the density is 650kg / m^3. The radius of the cylindrical core is 0.2m, the center is located at (0.1,0,-2), the P wave velocity is 245.26m / s, the S wave velocity is 100.12m / s, and the density is 1888kg / m^3. Figure 3 The displacement field component results obtained by the three methods when the model is incident are shown. There is good consistency between the three methods, which verifies the accuracy of the NMM method in simulating the half-space structure with inhomogeneous double layers. At the same time, the following table shows the computational cost comparison of the three different methods. The results show that the NMM method is very efficient when simulating simple layered inhomogeneous regions, with the shortest computation time and the lowest memory usage.

[0074] Solver Degrees of Freedom CPU time(s) Memory (GB) NMM 13380 68 0.26 Hybridization method 250749 224 3.54 COMSOL 2627607 908 166

[0075] Example 2: Two uniform half-spaces containing thin layers

[0076] In this embodiment, the model consists of two uniform half-spaces with a constant impedance Z in between. n The center is at the origin, and the cross-sectional dimensions are 1.2 meters by 1.2 meters. The P wave velocity in region 1 is 1201m / s, the S wave velocity is 532.53m / s, and the density is 1481kg / m^3. The P wave velocity in region 2 is 3500m / s, the S wave velocity is 2000m / s, and the density is 2400kg / m^3. The operating frequency is 1201Hz. Figure 4 As shown in Figure 2, when the impedance of the thin layer changes in the range [-300, 300] GPa, the phase calculated by the NMM method is compared with the phase calculated by the analytical solution formula, showing an excellent match. Therefore, we can design an elastic gradient metasurface that can be obtained by Z n The change of the coordinates in the cross section can produce a phase shift change in the range of (-π,π).

[0077] Example 3: Two uniform half-spaces with gradient metasurfaces

[0078] Assuming that the P wave is incident vertically, a gradient metasurface is designed to obtain a reflection angle (θ r ,φ r )=(50.28°,0) The expected phase shift changes with the x coordinate as Figure 5 Substituting the expression of ψ(x) into the elastic parameter expression, we get Z n The detailed expression of (x) and Figure 5 The drawing is done in red.

[0079] Next, we verified the accuracy of the numerical method provided in simulating the gradient metasurface, and observed whether abnormal reflection phenomena occurred on the designed gradient metasurface through numerical solution. Figure 6 As shown in , we used three different methods to simulate the metasurface: numerical pattern matching method, a hybrid spectral element method and traditional finite element method for comparison. The displacement field components obtained by the above three methods have achieved good consistency, which verifies the accuracy of our proposed numerical pattern matching method in simulating the vertical layered model containing gradient metasurface. Finally, as Figure 7 As shown, the divergence diagram of the displacement field of the reflected wave shows the distribution of the reflected P wave, and its reflection angle is The calculation results achieved the expected goals and verified the correctness of our design method.

Claims

1. A numerical algorithm for efficiently simulating elastic metasurfaces in vertically layered elastic models, characterized in that the method comprises the following steps: 1.1 analyzing the phase change caused by the change of thin layer parameters, and designing the elastic metasurface based on this; 1.2 constructing a numerical mode matching (NMM) method to improve the simulation efficiency by simplifying the three-dimensional elastic wave problem into a two-dimensional eigenvalue problem in the horizontal direction and an analytical expression in the vertical direction; 1.3 using the spectral element method (SEM) to calculate the eigenmodes of each layer; 1.4 The incident wave field is represented as a linear combination of the eigenmode and the excitation coefficient and propagator in the vertical direction, so as to achieve an accurate description of the wave field; 1.5 The boundary conditions are applied at the interface to determine the reflection matrix and the transmission matrix.

2. The numerical algorithm for efficiently simulating elastic metasurfaces in vertically layered elastic models according to claim 1, characterized in that: For two half-spaces connected by a thin elastic layer, when the P wave is incident vertically, the reflection coefficient can be expressed analytically as: Where Z1=ρ1v p,1 and Z2=ρ2v p,2 They represent the P wave impedance of regions 1 and 2 respectively, and ω is the angular frequency. At this time, the reflection coefficient R PP The amplitude and phase can be expressed as: R PP =|R PP |[cos(θ1-θ2)+jsin(θ1-θ2)] in, and Therefore, the incident wave is reflected from the surface with a phase shift of Therefore, the design of the elastic metasurface is based on the thin layer parameter Z = diag ([Z t ,Z t ,Z n ]) can deflect the direction of wave propagation to achieve the reflection or refraction effect of the target.

3. The numerical algorithm for efficiently simulating elastic metasurfaces in vertically layered elastic models according to claim 2, characterized in that: The reflection phenomenon caused by the phase change includes abnormal reflection, and its reflection direction follows the generalized Snell's law: in, and are the reflection phase gradients of the surface in the x and y directions, respectively. In order to achieve the desired reflection angle (θ r ,φ r ), the analytical expression of ψ(x,y) can be calculated according to the above formula. The thin layer parameters include elastic modulus and thin layer thickness, and these parameters vary in spatial coordinates. where ψ(x,y) is the desired phase change.

4. The numerical algorithm for efficiently simulating elastic metasurfaces in vertically layered elastic models according to claim 1, characterized in that: Each layer of the region needs to obtain the characteristic mode pair through the spectral element method where γ z is the propagation constant, is the corresponding characteristic function: Among them, Φ(x,y) is the row vector composed of node basis functions, u j is a column vector containing the values ​​of the characteristic function at discrete points. The stress field is discretized as: Using the basis function test on the stress expression, we can get the column vector expression of the stress at discrete points: T ij =M0 -1 M ijk u k Where, M0=<Φ,Φ>, Where <,> represents the inner product.

5. The numerical algorithm for efficiently simulating elastic metasurfaces in vertically layered elastic models according to claim 4, characterized in that: definition and is a matrix consisting of a set of characteristic functions of region N, Γ N,z It is a diagonal matrix of propagation constants corresponding to the characteristic function. The superscripts "+" and "-" indicate the propagation direction. Define the unknown quantity to be solved: reflection matrix R 1,2 , the transmission matrix T 1,2 , excitation coefficient a1.

6. The numerical algorithm for efficiently simulating elastic metasurfaces in vertically layered elastic models according to claim 1, characterized in that: The incident displacement field and the stress normal component of the incident field in region 1 can be expressed as: The reflected displacement field and the stress normal component of the reflected field in region 1 can be expressed as: The transmitted displacement field and the stress normal component of the transmitted field in region 2 can be expressed as: Therefore, the total field in region 1 can be expressed as the incident field plus the reflected field. Therefore, the expressions for the normal component of the displacement field and stress in region 1 are as follows: The total field in region 2 can be expressed as the transmitted field. Therefore, the expressions for the displacement field and the normal component of the stress in region 2 are as follows:

7. The numerical algorithm for efficiently simulating elastic metasurfaces in vertically layered elastic models according to claim 6, characterized in that: We need to calculate the excitation coefficient a1 corresponding to the incident field from region 1 based on the given incident field expression u0(x,y,z). At the incident surface z=z0, we get: After testing this expression with each characteristic function in turn, the excitation coefficients are obtained by solving the following system of linear equations:

8. The numerical algorithm for efficiently simulating elastic metasurfaces in vertically layered elastic models according to claim 6, characterized in that: We need to determine the reflection matrix and the transmission matrix by applying boundary conditions at the interface between the two layers. Since the thickness of the elastic thin layer is much smaller than the wavelength, we introduce a sliding boundary condition to simulate the elastic thin layer as follows: Wherein, K is the thin layer material parameter, h is the thin layer thickness, and the elastic thin layer parameter Z=K / h.

9. The numerical algorithm for efficiently simulating elastic metasurfaces in vertically layered elastic models according to claim 8, characterized in that: At the interface z=z1, we obtain the relationship: Using the characteristic function of region 2 as the test function for each relation we obtain the linear system: The expression of the coefficient matrix is ​​as follows: By solving the linear equations, the reflection matrix and transmission matrix can be obtained.

10. The numerical algorithm for efficiently simulating elastic metasurfaces in vertically layered elastic models according to claim 6, characterized in that: The numerical pattern matching method includes utilizing multiple eigenmodes to improve the simulation accuracy of propagation characteristics of elastic waves in complex multilayer structures.

Images