Electromagnetic particle adaptive simulation method and system for non-uniform and conformal grid

Abstract

The application discloses a non-uniform and conformal grid-oriented electromagnetic particle adaptive simulation method and system, and relates to the technical field of electromagnetic particle simulation. The method comprises the following steps: performing non-uniform right-angle grid division based on the average free path of charged particles, and identifying conformal grid units; adopting macro particles for initialization; performing electromagnetic field propulsion, field interpolation, particle propulsion, charge and current backfilling and collision simulation in each time step; and realizing accelerated calculation through a CPU and GPU hybrid parallel architecture. Through the application, the numerical simulation precision can be improved, the adaptability of complex geometry can be enhanced, the stability of electromagnetic field and particle coupling can be ensured, the information interaction between the field and the particles can be improved, and the large-scale calculation efficiency can be improved.

Description

Technical Field

[0001] This invention relates to the field of electromagnetic particle simulation technology, and in particular to an adaptive simulation method and system for electromagnetic particles oriented towards non-uniform and conformal meshes. Background Technology

[0002] Numerical simulation of the interaction between charged particles and electromagnetic fields is an important research tool in fields such as plasma physics, vacuum electronics, micro / nano fabrication, and electric propulsion. Computer simulations of particle motion and electromagnetic field evolution can provide support for device design, process optimization, and fundamental mechanism research.

[0003] Among them, the Particle-In-Cell (PIC) method is a widely used numerical simulation method. This method utilizes macroparticles to represent a large number of real charged particles and solves the electromagnetic field equations on a discrete grid, thus achieving the coupling of particle dynamics and electromagnetic field evolution. Since its inception, the PIC method has achieved numerous applications under conditions of uniform grid partitioning, regular geometric boundaries, and medium-scale computation.

[0004] As research objects increasingly involve larger spatial scales, more complex geometric structures, and longer evolutionary processes, higher demands are placed on numerical simulation methods. In these applications, how to ensure the reliability of simulation results while maintaining computational efficiency has become a common concern in related fields. Summary of the Invention

[0005] To address the challenge of achieving reliable simulation of the interaction between electromagnetic fields and charged particles under complex geometries and multi-scale conditions, balancing simulation accuracy and computational efficiency, this invention proposes an adaptive simulation method and system for electromagnetic particles in non-uniform and conformal meshes.

[0006] The present invention achieves the above objectives through the following technical solutions:

[0007] An adaptive simulation method for electromagnetic particles in non-uniform and conformal meshes, the method comprising:

[0008] Based on the mean free path distribution of charged particles, the computational domain is partitioned into non-uniform rectangular grids. Grid cells that intersect with curved or inclined boundaries are identified as conformal grid cells, while the rest are regular grid cells.

[0009] Multiple charged particles with the same properties are represented by macroparticles. The spatial distribution and velocity distribution of the macroparticles are initialized based on a non-uniform rectangular grid, and the electromagnetic field is initialized.

[0010] Within each time step, the following iterative calculation process is executed sequentially:

[0011] The electromagnetic field is advanced by solving Maxwell's equations using the finite-difference time-domain method. Specifically, the regular grid cells are processed using a standard difference scheme, the conformal grid cells are processed using conformal techniques, and time filtering is applied to the electromagnetic field to suppress numerical noise.

[0012] The advanced electromagnetic field is mapped to the position of macroparticles; for macroparticles located in regular grid cells, Clouds-in-Cells or linear weighting methods are used for mapping, while for macroparticles located in conformal grid cells, an improved weighting function based on least squares or radial basis functions is used for mapping.

[0013] Based on the electric and magnetic fields obtained from the mapping, the velocity and position of the macroparticles are updated according to the Lorentz force law, and the updated charge and current density of the macroparticles are backfilled into the non-uniform right-angled grid.

[0014] The Monte Carlo collision method and pseudo-collision model are used to simulate the collision process between macroparticles and background neutral particles, and the motion state of the macroparticles that collide is updated.

[0015] Throughout the iterative computation process, a hybrid parallel architecture of CPU and GPU is used to accelerate the computation.

[0016] As a preferred embodiment of the present invention, the step of partitioning the computational domain into a non-uniform rectangular grid based on the mean free path distribution of charged particles includes:

[0017] Input the three-dimensional geometry of the device to be simulated and the gas pressure and gas temperature required for the calculation;

[0018] Calculate the mean free path distribution of different types of charged particles based on gas pressure and gas temperature;

[0019] Set the grid scale such that the size of each grid cell is less than or equal to the product of the empirical coefficient and the corresponding mean free path;

[0020] The mesh generation algorithm is called to perform non-uniform right-angle mesh partitioning on the computational domain. Dense mesh is set in the region where the mean free path of charged particles is less than or equal to a preset first threshold and the particle density gradient is greater than a preset second threshold, and sparse mesh is set in the region where the mean free path of charged particles is greater than the preset first threshold and the particle density gradient is less than or equal to the preset second threshold.

[0021] The specific grid generation algorithm is either the octree partitioning method or the forward push method.

[0022] As a preferred embodiment of the present invention, multiple charged particles with the same characteristics are equivalently represented by macroparticles, and the spatial distribution and velocity distribution of the macroparticles are initialized based on a non-uniform rectangular grid, including:

[0023] Real charged particles are equivalent to macroparticles, and the equivalent charge and equivalent mass of macroparticles are set according to the number density of charged particles in different grid cells.

[0024] On a non-uniform rectangular grid, the spatial distribution of macroparticles is initialized based on the charged particle number density using a quiet multi-beam initialization method; specifically, the target density function of the preset charged particle number density distribution is initialized. Normalization is performed to obtain the normalized distribution function. Calculate the cumulative distribution function of the normalized distribution function. And divide the number of macroparticles within the grid cell into The bundle is obtained by solving the cumulative distribution function. To determine the initial spatial position of each macroparticle beam. ,satisfy ,in The sequence number of the bundle. ;

[0025] The velocity distribution of macroparticles is initialized on a non-uniform rectangular grid by: calculating the thermal velocity based on the gas temperature and the mass of the charged particles, obtaining a standard normal distribution variable using a random number generation method, scaling the standard normal distribution variable according to the thermal velocity to obtain the velocity component, so that the macroparticle velocity of regular grid cells follows a Maxwell-Boltzmann distribution, and the macroparticle velocity of conformal grid cells follows a drift Maxwell distribution;

[0026] Initialize the electromagnetic field on a non-uniform rectangular grid: initialize the electric and magnetic fields to zero, or add an external static field or a preset field distribution, and set boundary conditions.

[0027] As a preferred embodiment of the present invention, the regular mesh cells are processed using a standard difference format, and the conformal mesh cells are processed using conformal techniques, including:

[0028] On the regular grid cells, the electric and magnetic fields are updated using a central difference scheme based on the Yee grid structure. The update coefficients of the electric and magnetic fields are determined based on the dielectric constant, conductivity, permeability and magnetic loss parameters.

[0029] Within a conformal mesh cell, the effective edge length and effective surface area intersecting the boundary are introduced to correct and update the magnetic field. When the conformal mesh cell contains multiple media, the permittivity and permeability are weighted and averaged according to the area occupied by each medium. In conformal mesh cells where the geometric mismatch exceeds a preset third threshold, a backward weighting scheme is used to correct the magnetic field update. The weighting factor is... The value range is 0.5–0.9.

[0030] As a preferred embodiment of the present invention, the step of applying time filtering to the electromagnetic field to suppress numerical noise specifically includes:

[0031] After updating the electric and magnetic fields at each time step, a Friedman-Rambo time-law filter is applied to the electric and / or magnetic field components. This Friedman-Rambo time-law filter introduces a weighted combination of historical electric field quantities into the electromagnetic field propagation formula, where the filter parameters... The value range is 0.1–0.3.

[0032] As a preferred embodiment of the present invention, when mapping the propelled electromagnetic field to the macroparticle position, the macroparticles located in the conformal grid cells are mapped using an improved weighting function constructed based on the least squares method or radial basis functions, including:

[0033] Centered on the macroparticle's location, a local support region is selected, encompassing the conformal mesh cell containing the macroparticle and its adjacent mesh cells. N mesh nodes are then selected from this local support region as interpolation nodes, denoted as the coordinate set. And obtain the electric field component value and magnetic field component value of each interpolation node;

[0034] When using the radial basis function method:

[0035] Let the radial basis functions be ,in Indicates the position to be interpolated With the Nodes The Euclidean distance;

[0036] Constructing the kernel matrix , element is ;

[0037] Let the nodal field value vector be ;

[0038] By solving the system of linear equations Obtain the coefficient vector Then calculate the field value at the position of the macroparticle. ,in Indicates the electric field at the th Nodes Interpolated value of position;

[0039] When using the least squares method:

[0040] Let the basis function vector be The interpolation function is: ,in Let be the vector of coefficients to be determined;

[0041] Constructing a matrix , No. Behavior ;

[0042] Let the nodal field value vector be The weight matrix is ,in , It is a positive constant;

[0043] By solving equations Obtain the coefficient vector Then calculate the field value at the position of the macroparticle. .

[0044] As a preferred embodiment of the present invention, the step of updating the velocity and position of the macroparticle based on the mapped electric and magnetic fields according to the Lorentz force law includes:

[0045] With macroparticles at time speed Using the initial value, first perform half-step electric field acceleration to obtain the intermediate velocity. ;

[0046] Calculate the rotation vector and rotation factor Based on the rotation vector and rotation factor For intermediate speed Rotate the object to obtain its speed. ;

[0047] Perform half-step electric field acceleration again to obtain the updated velocity. ;

[0048] Utilize the updated speed Proceeding the macroparticle's position to time ,Right now ;

[0049] In the formula, , They represent macroparticles respectively At time step Time step Spatial location vector; For time step.

[0050] As a preferred embodiment of the present invention, when the updated charge and current density of the macroparticles are backfilled into the non-uniform right-angled grid, the charge carried by the macroparticles and the current component corresponding to the displacement are allocated to the corresponding nodes according to the weight function between the macroparticles and the grid nodes. For macroparticles located in regular grid cells, the weight function is still calculated using the Clouds-in-Cells method or the linear weight allocation method.

[0051] For macroparticles located within conformal mesh cells, the weighting function is still calculated using an improved weighting function constructed based on the least squares method or radial basis functions.

[0052] As a preferred embodiment of the present invention, the Monte Carlo collision method and a pseudo-collision model are used to simulate the collision process between macroparticles and background neutral particles, and the motion state of the macroparticles that collide is updated, including:

[0053] The total collision probability is calculated based on the macroparticle energy and the preset energy-dependent collision cross section.

[0054] A pseudo-collision model is used to determine the candidate macro-particle set, and the collision probability of the candidate macro-particles is compared with a random number to determine whether it is a real collision or a pseudo-collision.

[0055] For macroparticles that undergo real collisions, the collision type is determined using random numbers, and the motion state of the macroparticles is updated according to the collision type, wherein:

[0056] In an elastic collision, the magnitude of the velocity remains constant, only the direction of the velocity changes;

[0057] When a collision is triggered, the excitation energy is subtracted from the macroparticle's energy, and the velocity magnitude is reduced.

[0058] During ionization collisions, ionization energy is subtracted from the macroparticle's energy, and new electronic macroparticles and ion macroparticles are generated.

[0059] An adaptive simulation system for electromagnetic particles in non-uniform and conformal meshes, the system comprising:

[0060] The mesh generation module is used to perform non-uniform rectangular meshing of the computational domain based on the mean free path distribution of charged particles, and to identify mesh cells that intersect with curved or inclined boundaries as conformal mesh cells.

[0061] The particle initialization module is used to represent multiple charged particles with the same properties using macroparticles as equivalents, and to initialize the spatial distribution and velocity distribution of macroparticles based on the non-uniform rectangular grid, while also initializing the electromagnetic field.

[0062] The electromagnetic field propulsion module is used to solve Maxwell's equations using the finite-difference time-domain method in each time step. It uses the standard difference scheme for regular mesh cells, the conformal technique for conformal mesh cells, and applies time filtering to the electromagnetic field.

[0063] The field mapping module is used to interpolate electromagnetic fields to the macroparticle positions. For macroparticles located in regular grid cells, Clouds-in-Cells or linear weighting methods are used, while for macroparticles located in conformal grid cells, an improved weighting function based on least squares or radial basis functions is used.

[0064] The particle propulsion module is used to update the velocity and position of macroparticles according to the Lorentz force law based on the mapped electric and magnetic fields.

[0065] The charge and current backfill module is used to backfill the updated macroparticle charge and current density into the non-uniform right-angled grid.

[0066] The collision simulation module is used to simulate the collision process between macroparticles and background neutral particles using the Monte Carlo collision method and pseudo-collision model, and to update the motion state of the macroparticles that collide.

[0067] The parallel computing module is used to accelerate the calculation of electromagnetic field evolution and particle motion using a hybrid parallel architecture of CPU and GPU throughout the iterative calculation process.

[0068] The beneficial effects of this invention are as follows: Based on the non-uniform rectangular mesh partitioning of the mean free path distribution of charged particles, regions with large particle density gradients or small mean free paths achieve higher resolution, thereby enhancing the ability to characterize the physical processes in key regions without significantly increasing the overall computational load. By identifying cells intersecting curved or inclined boundaries as conformal mesh cells and combining conformal techniques to process the electromagnetic field, the fitting effect of electromagnetic field boundary conditions can be effectively improved, reducing errors caused by geometric discretization. Introducing a time-filtering operator during electromagnetic field propagation can suppress the accumulation of numerical noise during long-term iterations, improving the stability of the electromagnetic field-particle interaction process. By using Clouds-in-Cells or linear weighting methods in regular mesh cells and improved weighting functions based on least squares or radial basis functions in conformal mesh cells, accurate mapping of the electromagnetic field at particle positions is achieved, ensuring accurate backfilling of charge and current densities. The use of a hybrid parallel architecture of CPU and GPU enables parallel accelerated computation of electromagnetic field evolution and particle motion, which is beneficial for large-scale three-dimensional long-term simulations. Attached Figure Description

[0069] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. Wherein:

[0070] Figure 1 This is a flowchart of the method of the present invention;

[0071] Figure 2 This is a schematic diagram of the modular structure of the system in an embodiment of the present invention. Detailed Implementation

[0072] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the described embodiments of the present invention are within the scope of protection of the present invention.

[0073] like Figure 1 The illustration shows an embodiment of the present invention, which provides an adaptive simulation method for electromagnetic particles in non-uniform and conformal meshes, enabling efficient and accurate numerical simulations in semiconductor processes and vacuum electronic device design. For ease of understanding, the implementation process of the present invention can be divided into a preprocessing and initialization stage, a main loop stage, and a parallel computation and post-processing stage, as detailed below:

[0074] Phase 1: Preprocessing and Initialization

[0075] S1: Computational Domain Discretization and Non-Uniform Mesh Generation

[0076] S11: Perform geometric modeling of the device to be simulated. Input parameters such as the three-dimensional CAD geometry of the device to be simulated, initial plasma density, gas type, gas pressure, temperature, and electromagnetic field (electrostatic field and static magnetic field) as prerequisites for the simulation method.

[0077] "Device to be simulated" refers to semiconductor process devices or vacuum electronic devices that require simulation analysis. Examples include plasma etching chambers, reaction chambers in sputtering coating equipment, vacuum tubes, microwave power devices, electron beam devices, and other three-dimensional device structures involving the interaction of charged particles with electromagnetic fields. It is typically defined by inputting its three-dimensional CAD model and material and electrode boundary conditions. Gas pressure can be obtained from the process parameters of the vacuum chamber, in Pa; temperature can be determined by experimental settings or the external environment, in K.

[0078] S12: Non-uniform mesh partitioning based on particle mean free path

[0079] The mesh size should be appropriate for the mean free path of the charged particles. The mean free path of the charged particles is calculated using gas dynamics theory based on gas pressure and temperature. The formula is as follows:

[0080] ;

[0081] In the formula, Boltzmann's constant; The gas temperature; The effective diameter of the molecule (which can be obtained by looking up the scattering cross-section of different gas molecules in a database); This represents gas pressure.

[0082] In actual simulations, the mean free path distribution is calculated for different types of charged particles (such as electrons and ions), and a correspondence with spatial coordinates is established. This yields a mean free path field covering the entire computational domain.

[0083] After obtaining the mean free path distribution, set the constraints on the mesh size: ,in, The side length of the grid cell; This is the mean free path of the region; This is an empirical coefficient, usually taken as... This constraint ensures that the mesh can distinguish the motion characteristics of charged particles. During mesh generation, the software automatically checks the size of each mesh according to this formula and refines the mesh as necessary.

[0084] Under the condition of satisfying the above scale constraints, the octree partitioning method or the forward push method is used to mesh the three-dimensional computational domain:

[0085] Octree partitioning: recursively partitions a 3D cube into elements, when the element size does not meet the requirements. Continue subdividing until all elements satisfy the constraints;

[0086] Advance-front method: Start by generating the mesh from the boundary step by step, determine the mesh size based on the mean free path and density distribution, and advance inward step by step until the entire computational domain is filled.

[0087] The generated non-uniform rectangular mesh has the following characteristics:

[0088] When the mean free path is less than or equal to a preset first threshold And the particle density gradient is greater than the preset second threshold. In the area, a dense grid is arranged;

[0089] When the mean free path is greater than And the particle density gradient is less than or equal to In the area, a sparse grid is arranged.

[0090] After completing the non-uniform rectangular mesh generation, check the relationship between each mesh cell and the device boundary:

[0091] If a cell is located entirely inside or outside the device, it is directly identified as a regular grid cell.

[0092] If the boundary of a cell intersects with the curved surface or inclined boundary of the device, the cell is identified as a conformal mesh cell.

[0093] For the identified conformal mesh cells, the effective edge length and effective surface area need to be calculated further. These parameters will be used to correct the difference formula during the electromagnetic field iteration process in order to achieve conformal FDTD calculation.

[0094] S2: Particle Initialization

[0095] S21: Macroparticle Generation

[0096] To reduce computational complexity, "macroparticles" are used to represent a large number of real charged particles with identical properties. Each macroparticle can be considered as a group of particles with the same coordinates and velocities, and a sufficient number of macroparticles are needed in each grid to ensure computational accuracy. For each grid cell, the total number of charged particles in the region is first calculated (charged particle number density × cell volume). Then, according to the set number of macroparticles, the number of real particles represented by each macroparticle is determined (i.e., the weight, obtained by dividing the total number of charged particles by the number of macroparticles). This determines the equivalent charge and equivalent mass of the macroparticle.

[0097] For example, for each particle species Set the target number of macroparticles (For example (The ion concentration is also set to 64).

[0098] For each grid cell With particle species The actual number of particles in the calculation unit ,in For particle species In grid cells number density on For grid cells Volume;

[0099] Setting macro particle weights ;

[0100] Setting the equivalent charge of macroparticles With equivalent quality : ;in, The valence state of the particle species (electrons are -1, monovalent ions are +1). For elementary charge, For particle species mass;

[0101] For each grid cell With particle species generate Instance of macro particles and will The value is assigned to the equivalent quantity given in the above formula.

[0102] S22: Initial Coordinate Assignment

[0103] On a non-uniform rectangular grid, the spatial distribution of macroparticles is initialized based on the charged particle number density using a quiet-beam initialization method (QuietStart's Multiple-Beam method); specifically:

[0104] Let the total number of macroparticles be Divided into Bundles, each bundle contains One macroparticle;

[0105] Define the target density function of the charged particle number density distribution. ,For example (A typical "parabolic distribution") For the center of maximum number density, The length of the region;

[0106] Normalization yields And calculate the cumulative distribution function. ;

[0107] Divide the interval into equal parts Part, No. Position of macroparticles From the equation , The position of the macroparticle is determined by inverse calculation to ensure that it accurately follows the target density function. ;

[0108] In the formula, It is a grid cell Spatial range; and In Both represent spatial coordinates, which in the one-dimensional case are along a certain direction (e.g.) The actual position variable (direction). It is a virtual integral variable in integration, and its physical meaning is... Both are spatial coordinates, but they are used in the integral sign. This indicates an attempt to avoid exceeding the points cap. Confusing.

[0109] Its function is to map "probability value" to "spatial location", thus performing a coordinate reverse calculation function.

[0110] Calculate the weighted equivalent total number of macroparticles across the entire computational domain, for all grid cells and all particle species. and the total number of actual charged particles If the two are equal, it means that the initialization correctly reproduced the distribution of real charged particles.

[0111] S23: Initial Velocity Distribution

[0112] First, based on the gas temperature and the mass of charged particles Calculate the thermal velocity of this particle species: ,in Boltzmann's constant;

[0113] Subsequently, standard normally distributed variables are generated using random number generation methods. For example, when using the Box-Muller transform, two independent uniformly distributed random numbers are first generated. Two independent standard normal distribution samples were obtained through calculation: ;

[0114] For macroparticles in regular mesh cells, their three-dimensional velocity components are obtained by thermal velocity scaling:

[0115] ;

[0116] in, It can be generated by the next Box-Muller transformation. The velocity components obtained from this are... They both obey the Maxwell-Boltzmann distribution, ensuring that the velocity field of macroparticles is consistent with the thermal equilibrium state of the gas.

[0117] For macroparticles in conformal mesh elements, their velocity distribution adopts a drift Maxwell distribution. Drift velocity is introduced based on the above thermal velocities: ,in The drift coefficient, The unit vector representing the drift direction can be either the direction of the applied electric field or the boundary normal, thus making the overall velocity distribution of macroparticles conform to a drift Maxwell distribution, more accurately reflecting the dynamic characteristics of particles near the boundary.

[0118] Through the above steps, the macroparticle velocities in regular grid cells are initialized to the Maxwell-Boltzmann distribution, and the macroparticle velocities in conformal grid cells are initialized to the drifting Maxwell distribution, thereby ensuring the consistency and physical rationality of the velocity distribution in the entire non-uniform rectangular grid system.

[0119] S3: Electromagnetic Field and Boundary Condition Initialization

[0120] S31: Field Initialization

[0121] At each grid point of the completed non-uniform right-angled grid, the electric field is... and magnetic field Assign a value.

[0122] Common initialization methods include:

[0123] Zero-field initialization: Electric field of all grid points and magnetic field Set to zero;

[0124] External field initialization: Depending on the problem setting, import the static (or sometimes dynamic) electromagnetic field distribution generated by an external source. For example, initialize a constant external magnetic field. , or an electrostatic field generated by electrode potential.

[0125] S32: Boundary Condition Settings

[0126] At the boundary of the computational domain, set boundary conditions for the electromagnetic field, for example:

[0127] Ideal electric conductor (PEC): The tangential component of the electric field is zero;

[0128] Ideal magnetic conductor (PMC): The tangential component of the magnetic field is zero;

[0129] Absorbing Boundary (ABC): Simulates an open boundary to prevent non-physical reflections;

[0130] Periodic boundary (PBC): Simulates an infinitely large periodic structure.

[0131] Phase Two: Main Loop Phase (each time step) (Execution steps)

[0132] S4: Electromagnetic Field Propulsion

[0133] In this embodiment, the finite-difference time-domain method is used to numerically advance Maxwell's equations, combining three steps: regular grid cells, conformal grid cells, and time filtering. The specific steps are as follows:

[0134] S41: On regular grid cells, a Yee grid structure is used for spatial discretization. Electric and magnetic field components are stored alternately in space and progressively updated in time. The update formula for the electric field components is:

[0135] ;

[0136] Among them, the update coefficient , Determined by dielectric constant and conductivity:

[0137] ;

[0138] The update formula for the magnetic field components is:

[0139] ;

[0140] Among them, the update coefficient , Determined by permeability and magnetic loss parameters:

[0141] ;

[0142] In the formula, Indicates time step Time and space The electric field components; These are the magnetic field components at adjacent locations; The grid step size; For unit conductivity, For unit The dielectric constant; For time step; Indicates time step Time and space Magnetic field components; For unit magnetic permeability, For unit The magnetic loss coefficient (equivalent magnetic conductivity).

[0143] The above formulas guarantee numerical stability and consistency of energy dissipation in the presence of conductivity and magnetic loss.

[0144] S42: For conformal mesh elements intersecting curved or inclined boundaries, effective geometric quantities must be calculated:

[0145] Effective edge length , represents the actual edge length located inside the computational domain in the cell;

[0146] Effective surface area , which represents the actual area within the computational domain of the cell.

[0147] The above effective geometric quantity correction formula is introduced into the magnetic field update, for example:

[0148] ;

[0149] In the formula, Indicates and Effective cross-sectional area corresponding to the direction The effective edge length of the cell edge within the computational domain; This represents the electric field along the edge tangentially; It still represents magnetic permeability;

[0150] When a conformal mesh element contains multiple media, the equivalent parameters are obtained by weighting the dielectric constant, permeability, and other parameters according to the area occupied by each medium on the element surface:

[0151] ;

[0152] In the formula, Indicates medium The area occupied by the unit surface , as medium The dielectric constant and permeability; For the quantity of media;

[0153] For conformal mesh elements whose geometric mismatch exceeds a preset third threshold, this embodiment employs backward weighting correction, and the updated magnetic field value is given by the following formula:

[0154] ;

[0155] in, The magnetic field value is calculated based on conformal correction. This is a weighting factor, with a value range of 0.5–0.9; this range was determined in numerical experiments by comprehensively considering stability and physical accuracy. If the correction is too weak, it will be insufficient to suppress numerical instability; if If the correction is too strong, it will cause the energy to decay too quickly.

[0156] S43: After updating the electric and magnetic fields at each time step, a time-domain filter is applied to the field quantities. This embodiment uses the Friedman-Rambo time-domain filter operator, with magnetic induction intensity as the primary factor. For example:

[0157] ;

[0158] in The filter parameter has a value range of 0.1–0.3; the selection criteria are as follows:

[0159] when At that time, the filtering effect was insufficient and could not effectively suppress high-frequency numerical noise;

[0160] when At times, excessive filtering may interfere with low-frequency physics.

[0161] Therefore By limiting it to the range of 0.1–0.3, high-frequency noise can be effectively suppressed while ensuring low-frequency fidelity.

[0162] The filtering target can be the electric field component, the magnetic field component, or both simultaneously.

[0163] S5: Electromagnetic field interpolation (from grid to particle)

[0164] The advanced electromagnetic field is mapped to the macroparticle position; for macroparticles located in regular grid cells, the Clouds-in-Cells or linear weighting method is used for mapping, while for macroparticles located in conformal grid cells, an improved weighting function is used for mapping.

[0165] In this embodiment, for macroparticles located within conformal mesh cells, the improved weighting function includes two implementation methods: the radial basis function method and the least squares method, as detailed below:

[0166] S51: Based on the position of macroparticles Centered on the macroparticle, a conformal mesh cell containing the macroparticle and its adjacent mesh cells is selected as the local support region. N mesh nodes are selected from the local support region as interpolation nodes, denoted as the coordinate set. The corresponding electric or magnetic field component values ​​are denoted as ;

[0167] S52: When using the radial basis function method:

[0168] Choose radial basis functions as ,in Indicates the position to be interpolated With the Nodes The Euclidean distance; the radial basis functions are preferably Gaussian functions. ,in For shape parameters, the range of values ​​is related to the local average mesh size. Proportional, that is proportionality coefficient Between 0.5 and 2;

[0169] Constructing the kernel matrix Its elements are defined as , ;

[0170] Let the nodal field value vector be ;

[0171] By solving the system of linear equations Obtain the coefficient vector ;

[0172] The field value for calculating the position of macroparticles ,in Indicates the electric field at the th Nodes Interpolated value of position.

[0173] S53: When using the least squares method:

[0174] Let the basis function vector be For example, basis functions for first- or second-order polynomials;

[0175] The interpolation function is: ,in Let be the vector of coefficients to be determined;

[0176] Constructing a matrix , its first Behavior ;

[0177] Let the nodal field value vector be The weight matrix is ,in , It is a positive constant used to avoid division by zero;

[0178] By solving equations Obtain the coefficient vector ;

[0179] Calculate the field value at the position of macroparticles .

[0180] S6: Particle Propulsion (Migration)

[0181] Based on the electric and magnetic fields obtained from the mapping, the velocity and position of the macroparticles are updated according to the Lorentz force law.

[0182] In this embodiment, after the electromagnetic field is interpolated and mapped to the macroparticle's position, the macroparticle's motion state is updated using the Boris algorithm. This algorithm is a commonly used numerical integration method that ensures energy conservation during the rotational process under the influence of a magnetic field. The specific steps are as follows:

[0183] S61: Input macroparticle state, including: macroparticle state at time step Spatial position vector Macroparticles at time speed The charge of macroparticles and quality Electric field at the position of macroparticle With magnetic field and time step These quantities serve as inputs for the calculation.

[0184] S62: Half-step electric field acceleration (first acceleration)

[0185] First, consider the effect of the electric field on the particle. To implement the leap-frog scheme, the electric field acceleration is divided into two half-steps:

[0186] ;

[0187] This is equivalent to the electric field acting on the macroparticle within half a time step, thereby updating and obtaining an intermediate velocity. At this point, the influence of the magnetic field has not yet been taken into account.

[0188] S63: Magnetic field rotation treatment: The effect of a magnetic field on a particle is to cause its velocity to rotate without changing its magnitude.

[0189] Define the rotation vector: ;

[0190] Define the rotation factor ;

[0191] Calculate the intermediate speed: ;

[0192] Calculate the velocity after rotation: ;

[0193] This step is equivalent to changing the velocity vector. Rotate by an angle around the direction of the magnetic field, the angle of rotation being... The decision is made. The magnetic field only changes the direction of the velocity, not its magnitude; therefore, the energy remains constant.

[0194] S64: Second-stage electric field acceleration (second acceleration)

[0195] After the magnetic field rotates, an electric field is applied again for half-step acceleration: ;

[0196] Thus, the complete effect of the electric field is divided into the first half-step and the second half-step, ensuring time symmetry and numerical stability.

[0197] S65: Update macroparticle positions

[0198] Advance position based on the updated half-step speed: This step allows the macroparticle's position to move into the next integer time step. .

[0199] S7: Charge-current density distribution (from particle to grid)

[0200] The updated charge and current densities of macroparticles are then backfilled into the non-uniform rectangular grid.

[0201] In this embodiment, after the macroparticles update their velocity and position under the influence of electric and magnetic fields, their charge and current density need to be backfilled into the non-uniform rectangular grid to ensure self-consistent coupling between the electromagnetic field and the particles. The specific steps are as follows:

[0202] S71: For the first There are macroparticles, whose equivalent charge is... ,in, This represents the actual number of charged particles that the macroparticle represents. This represents the charge of a single real charged particle, which remains constant throughout the calculation.

[0203] S72: Charge Backfilling

[0204] Based on the position of macroparticles Determine the non-uniform rectangular grid cell in which it is located and its adjacent grid nodes;

[0205] For regular mesh cells, the Clouds-in-Cells (CIC) or linear weighting method is still used, and the weights are calculated based on the geometric distance between the macroparticle and its neighboring nodes. For conformal mesh cells, the improved weighting function based on the least squares method or radial basis function is still used, and its calculation method is the same as in S5.

[0206] Macroparticle charge The charge density increments of each node are obtained by assigning them to the nodes of the corresponding grid cell according to a weighting function. For example, for the th... Nodes: ,in, For a moment node charge density, macroparticles For nodes The weight, For nodes The corresponding grid control volume.

[0207] S73: Current density backfill

[0208] time step Inside, the trajectory of macroparticles is approximately linear: ;

[0209] Current density of macroparticles Given by the kinetic charge flow rate: ;

[0210] Similarly, using the CIC / linear weighting function for regular mesh elements and the improved weighting function for conformal mesh elements, Assigned to grid nodes: ;

[0211] After the charge and current density backfilling is completed, check that the discrete continuity equation is satisfied: This is to ensure the self-consistency of the numerical simulation.

[0212] S8: Particle Collision (MCC Method)

[0213] In this embodiment, the Monte Carlo collision method is used to simulate the interaction between macroparticles and background gas molecules, and a pseudo-collision model is combined to improve computational efficiency.

[0214] S81: Calculate the total collision probability based on the macroparticle energy and the preset energy-dependent collision cross-section, specifically:

[0215] For each macroparticle First, calculate the kinetic energy: ,in, For macroparticle mass, For macroparticle velocity vectors;

[0216] According to kinetic energy Read different types of collision cross sections from the database The total cross-section is obtained as follows: ,in, Indicates the type of collision (including elastic, excitation, ionization, etc.);

[0217] Macroparticles at time step The total collision probability within is: ,in For gas number density, This represents the speed magnitude.

[0218] S82: A pseudo-collision model is used to determine the candidate macroparticle set, and the collision probability of the candidate macroparticles is compared with a random number to determine whether it is a real collision or a pseudo-collision. Specifically:

[0219] Predetermine the global maximum collision frequency: ;

[0220] The corresponding maximum collision probability is: ;

[0221] In total macroparticle number In, random sampling One macroparticle is selected as a candidate macroparticle;

[0222] For each candidate macroparticle, generate a first random number. ,like If the collision is real, the particle is considered to have experienced a true collision; otherwise, it is considered a false collision and its state remains unchanged.

[0223] S83: Generate a second random number for macroparticles that undergo real collisions. The specific collision type is determined based on the proportion of various cross sections:

[0224] like Then it is the first Collision type.

[0225] S84: Update macroparticle state

[0226] Elastic collision: The velocity direction is recalculated, while the velocity magnitude remains approximately unchanged;

[0227] Inducing Collisions: From Kinetic Energy Subtracting the excitation energy, the velocity magnitude decreases;

[0228] Ionization collision: from kinetic energy The ionization energy is subtracted, and new electron and ion macroparticles are generated and given corresponding initial velocities.

[0229] For all macroparticles in each time step, execute S81–S84 sequentially until the entire simulation process is completed.

[0230] Phase 3: Parallel Computing and Post-processing

[0231] S9: Parallel Computing

[0232] Throughout the iterative computation process, a hybrid parallel architecture combining CPU and GPU is employed to accelerate computation, including:

[0233] Domain decomposition: The entire computing grid is divided into multiple subdomains, with each computing process (MPI) responsible for one subdomain.

[0234] Task assignment:

[0235] CPU side (MPI+OpenMP): Responsible for main logic control, FDTD field calculation, and inter-process communication. Within each process, OpenMP can be used to parallelize the particle loop through multiple threads.

[0236] GPU-based CUDA: The computationally intensive and logically simple parts (such as particle-related operations: field interpolation, force calculation, particle propulsion, charge and current distribution, and collision detection) are ported to the GPU, using thousands of threads to process each particle in parallel.

[0237] Data synchronization: After each time step, data exchange is performed between processes (transferring particle information across subdomain boundaries and synchronizing field and source term information at the subdomain boundaries).

[0238] S10: Output and Post-processing

[0239] Timed output: Writes data such as particle distribution, electromagnetic field distribution, and plasma parameters (density, temperature) to a file at set time intervals.

[0240] Visualization: Post-processing software is used to visualize the output data, showing results such as particle motion, electromagnetic field evolution, and plasma formation processes.

[0241] like Figure 2 As shown, another embodiment of the present invention provides an adaptive simulation system for electromagnetic particles in non-uniform and conformal meshes, comprising:

[0242] The mesh generation module is used to perform non-uniform rectangular meshing of the computational domain based on the mean free path distribution of charged particles, and to identify mesh cells that intersect with curved or inclined boundaries as conformal mesh cells.

[0243] The particle initialization module is used to represent multiple charged particles with the same properties using macroparticles as equivalents, and to initialize the spatial and velocity distribution of macroparticles based on a non-uniform rectangular grid, while also initializing the electromagnetic field.

[0244] The electromagnetic field propulsion module is used to solve Maxwell's equations using the finite-difference time-domain method in each time step. It uses the standard difference scheme for regular mesh cells, the conformal technique for conformal mesh cells, and applies time filtering to the electromagnetic field.

[0245] The field mapping module is used to interpolate electromagnetic fields to the macroparticle positions. For macroparticles located in regular grid cells, Clouds-in-Cells or linear weighting methods are used, while for macroparticles located in conformal grid cells, an improved weighting function based on least squares or radial basis functions is used.

[0246] The particle propulsion module is used to update the velocity and position of macroparticles according to the Lorentz force law based on the mapped electric and magnetic fields.

[0247] The charge and current backfill module is used to backfill the updated macroparticle charge and current density into a non-uniform right-angled grid.

[0248] The collision simulation module is used to simulate the collision process between macroparticles and background neutral particles using the Monte Carlo collision method and pseudo-collision model, and to update the motion state of the macroparticles that collide.

[0249] The parallel computing module is used to accelerate the calculation of electromagnetic field evolution and particle motion using a hybrid parallel architecture of CPU and GPU throughout the iterative calculation process.

[0250] The parallel computing module connects to all the above modules and is used to accelerate the execution of electromagnetic field updates, field mapping, particle propulsion, backfilling, and collision calculations in a CPU and GPU hybrid architecture.

[0251] In summary, this invention enables high-precision, adaptive, and efficient numerical simulation of the interaction process between electromagnetic fields and charged particles under complex geometric and multi-scale conditions.

[0252] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any person skilled in the art can easily conceive of various variations or substitutions within the technical scope disclosed in this application, and these should all be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. An adaptive simulation method for electromagnetic particles in non-uniform and conformal meshes, characterized in that, The method includes: Based on the mean free path distribution of charged particles, the computational domain is partitioned into non-uniform rectangular grids. Grid cells that intersect with curved or inclined boundaries are identified as conformal grid cells, while the rest are regular grid cells. Multiple charged particles with the same properties are represented by macroparticles. The spatial distribution and velocity distribution of the macroparticles are initialized based on a non-uniform rectangular grid, and the electromagnetic field is initialized. Within each time step, the following iterative calculation process is executed sequentially: The electromagnetic field is advanced by solving Maxwell's equations using the finite-difference time-domain method. Specifically, the regular grid cells are processed using a standard difference scheme, the conformal grid cells are processed using conformal techniques, and time filtering is applied to the electromagnetic field to suppress numerical noise. The advanced electromagnetic field is mapped to the position of macroparticles; for macroparticles located in regular grid cells, Clouds-in-Cells or linear weighting methods are used for mapping, while for macroparticles located in conformal grid cells, an improved weighting function based on least squares or radial basis functions is used for mapping. Based on the electric and magnetic fields obtained from the mapping, the velocity and position of the macroparticles are updated according to the Lorentz force law, and the updated charge and current density of the macroparticles are backfilled into the non-uniform right-angled grid. The Monte Carlo collision method and pseudo-collision model are used to simulate the collision process between macroparticles and background neutral particles, and the motion state of the macroparticles that collide is updated. Throughout the iterative computation process, a hybrid parallel architecture of CPU and GPU is used to accelerate the computation.

2. The adaptive simulation method for electromagnetic particles in non-uniform and conformal meshes according to claim 1, characterized in that, The method of partitioning the computational domain into a non-uniform rectangular grid based on the mean free path distribution of charged particles includes: Input the three-dimensional geometry of the device to be simulated and the gas pressure and gas temperature required for the calculation; Calculate the mean free path distribution of different types of charged particles based on gas pressure and gas temperature; Set the grid scale such that the size of each grid cell is less than or equal to the product of the empirical coefficient and the corresponding mean free path; The mesh generation algorithm is called to perform non-uniform right-angle mesh partitioning on the computational domain. Dense mesh is set in the region where the mean free path of charged particles is less than or equal to a preset first threshold and the particle density gradient is greater than a preset second threshold, and sparse mesh is set in the region where the mean free path of charged particles is greater than the preset first threshold and the particle density gradient is less than or equal to the preset second threshold. The specific grid generation algorithm is either the octree partitioning method or the forward push method.

3. The adaptive simulation method for electromagnetic particles in non-uniform and conformal meshes according to claim 2, characterized in that, Multiple charged particles with identical properties are represented by macroparticles. The spatial and velocity distributions of the macroparticles are initialized based on a non-uniform rectangular grid, including: Real charged particles are equivalent to macroparticles, and the equivalent charge and equivalent mass of macroparticles are set according to the number density of charged particles in different grid cells. On a non-uniform rectangular grid, the spatial distribution of macroparticles is initialized based on the charged particle number density using a quiet multi-beam initialization method; specifically, the target density function of the preset charged particle number density distribution is initialized. Normalization is performed to obtain the normalized distribution function. Calculate the cumulative distribution function of the normalized distribution function. And divide the number of macroparticles within the grid cell into The bundle is obtained by solving the cumulative distribution function. To determine the initial spatial position of each macroparticle beam. ,satisfy ,in The sequence number of the bundle. ; The velocity distribution of macroparticles is initialized on a non-uniform rectangular grid by: calculating the thermal velocity based on the gas temperature and the mass of the charged particles, obtaining a standard normal distribution variable using a random number generation method, scaling the standard normal distribution variable according to the thermal velocity to obtain the velocity component, so that the macroparticle velocity of regular grid cells follows a Maxwell-Boltzmann distribution, and the macroparticle velocity of conformal grid cells follows a drift Maxwell distribution; Initialize the electromagnetic field on a non-uniform rectangular grid: initialize the electric and magnetic fields to zero, or add an external static field or a preset field distribution, and set boundary conditions.

4. The adaptive simulation method for electromagnetic particles in non-uniform and conformal meshes according to claim 1, characterized in that, The regular mesh cells are processed using a standard difference scheme, and the conformal mesh cells are processed using conformal techniques, including: On the regular grid cells, the electric and magnetic fields are updated using a central difference scheme based on the Yee grid structure. The update coefficients of the electric and magnetic fields are determined based on the dielectric constant, conductivity, permeability and magnetic loss parameters. Within a conformal mesh cell, the effective edge length and effective surface area intersecting the boundary are introduced to correct and update the magnetic field. When the conformal mesh cell contains multiple media, the permittivity and permeability are weighted and averaged according to the area occupied by each medium. In conformal mesh cells where the geometric mismatch exceeds a preset third threshold, a backward weighting scheme is used to correct the magnetic field update. The weighting factor is... The value range is 0.5–0.

9.

5. The adaptive simulation method for electromagnetic particles in non-uniform and conformal meshes according to claim 4, characterized in that, The application of time filtering to the electromagnetic field to suppress numerical noise specifically includes: After updating the electric and magnetic fields at each time step, a Friedman-Rambo time-law filter is applied to the electric and / or magnetic field components. This Friedman-Rambo time-law filter introduces a weighted combination of historical electric field quantities into the electromagnetic field propagation formula, where the filter parameters... The value range is 0.1–0.

3.

6. The adaptive simulation method for electromagnetic particles in non-uniform and conformal meshes according to claim 1, characterized in that, When mapping the advanced electromagnetic field to the macroparticle positions, the macroparticles located in conformal mesh elements are mapped using an improved weighting function based on least squares or radial basis functions, including: Centered on the macroparticle's location, a local support region is selected, encompassing the conformal mesh cell containing the macroparticle and its adjacent mesh cells. N mesh nodes are then selected from this local support region as interpolation nodes, denoted as the coordinate set. And obtain the electric field component value and magnetic field component value of each interpolation node; When using the radial basis function method: Let the radial basis functions be ,in Indicates the position to be interpolated With the Nodes The Euclidean distance; Constructing the kernel matrix , element is ; Let the nodal field value vector be ; By solving the system of linear equations Obtain the coefficient vector Then calculate the field value at the position of the macroparticle. ,in Indicates the electric field at the th Nodes Position interpolation; When using the least squares method: Let the basis function vector be The interpolation function is: ,in Let be the vector of coefficients to be determined; Constructing a matrix , No. Behavior ; Let the nodal field value vector be The weight matrix is ,in , It is a positive constant; By solving equations Obtain the coefficient vector Then calculate the field value at the position of the macroparticle. .

7. The adaptive simulation method for electromagnetic particles in non-uniform and conformal meshes according to claim 1, characterized in that, The updating of the velocity and position of macroparticles based on the mapped electric and magnetic fields, according to the Lorentz force law, includes: With macroparticles at time speed Using the initial value, first perform half-step electric field acceleration to obtain the intermediate velocity. ; Calculate the rotation vector and rotation factor Based on the rotation vector and rotation factor For intermediate speed Rotate the object to obtain its speed. ; Perform half-step electric field acceleration again to obtain the updated velocity. ; Utilize the updated speed Proceeding the macroparticle's position to time ,Right now ; In the formula, , They represent macroparticles respectively At time step Time step Spatial location vector; For time step.

8. The adaptive simulation method for electromagnetic particles in non-uniform and conformal meshes according to claim 6, characterized in that, When the updated charge and current density of the macroparticles are backfilled into the non-uniform right-angled grid, the charge carried by the macroparticles and the current components corresponding to the displacement are allocated to the corresponding nodes according to the weight function between the macroparticles and the grid nodes. For macroparticles located in regular grid cells, the weight function is still calculated using the Clouds-in-Cells method or the linear weight allocation method. For macroparticles located within conformal mesh cells, the weighting function is still calculated using an improved weighting function constructed based on the least squares method or radial basis functions.

9. The adaptive simulation method for electromagnetic particles in non-uniform and conformal meshes according to claim 6, characterized in that, The Monte Carlo collision method and pseudo-collision model are used to simulate the collision process between macroparticles and background neutral particles, and the motion state of the macroparticles that collide is updated, including: The total collision probability is calculated based on the macroparticle energy and the preset energy-dependent collision cross section. A pseudo-collision model is used to determine the candidate macro-particle set, and the collision probability of the candidate macro-particles is compared with a random number to determine whether it is a real collision or a pseudo-collision. For macroparticles that undergo real collisions, the collision type is determined using random numbers, and the motion state of the macroparticles is updated according to the collision type, wherein: In an elastic collision, the magnitude of the velocity remains constant, only the direction of the velocity changes; When a collision is triggered, the excitation energy is subtracted from the macroparticle's energy, and the velocity magnitude is reduced. During ionization collisions, ionization energy is subtracted from the macroparticle's energy, and new electronic macroparticles and ion macroparticles are generated.

10. An adaptive simulation system for electromagnetic particles in non-uniform and conformal meshes, applied to the adaptive simulation method for electromagnetic particles in non-uniform and conformal meshes as described in claim 1, characterized in that, The system includes: The mesh generation module is used to perform non-uniform rectangular meshing of the computational domain based on the mean free path distribution of charged particles, and to identify mesh cells that intersect with curved or inclined boundaries as conformal mesh cells. The particle initialization module is used to represent multiple charged particles with the same properties using macroparticles as equivalents, and to initialize the spatial distribution and velocity distribution of macroparticles based on the non-uniform rectangular grid, while also initializing the electromagnetic field. The electromagnetic field propulsion module is used to solve Maxwell's equations using the finite-difference time-domain method in each time step. It uses the standard difference scheme for regular mesh cells, the conformal technique for conformal mesh cells, and applies time filtering to the electromagnetic field. The field mapping module is used to interpolate electromagnetic fields to the macroparticle positions. For macroparticles located in regular grid cells, Clouds-in-Cells or linear weighting methods are used, while for macroparticles located in conformal grid cells, an improved weighting function based on least squares or radial basis functions is used. The particle propulsion module is used to update the velocity and position of macroparticles according to the Lorentz force law based on the mapped electric and magnetic fields. The charge and current backfill module is used to backfill the updated macroparticle charge and current density into the non-uniform right-angled grid. The collision simulation module is used to simulate the collision process between macroparticles and background neutral particles using the Monte Carlo collision method and pseudo-collision model, and to update the motion state of the macroparticles that collide. The parallel computing module is used to accelerate the calculation of electromagnetic field evolution and particle motion using a hybrid parallel architecture of CPU and GPU throughout the iterative calculation process.

Images