3D CAD integrated plasma simulation tool

Abstract

The present disclosure describes a 3D CAD integrated simulation tool for modeling plasma discharges and electromagnetic phenomena. The system integrates multiple computational physics solvers including finite element method (FEM) electromagnetic solvers, particle-in-cell (PIC) kinetic solvers, fluid dynamics solvers, and reaction solvers with coordinated multi-time-step integration schemes. The software system provides an integrated 3D CAD environment for user interaction and specialized capabilities for low-temperature plasma applications, plasma-enhanced chemical vapor deposition (PECVD), plasma etching, sputter deposition, spacecraft charging, vacuum discharges, and RF plasma processing.

Description

Docket No.3528169.000302 3D CAD INTEGRATED SIMULATION TOOL CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims priority to U.S. Provisional Patent Application No. 63 / 708,272, titled "3D CAD Integrated Simulation Tool for RF Plasma Discharges with Multi-Solver Interface", filed October 17, 2024, which is hereby incorporated by reference in its entirety. FIELD

[0002] The present disclosure relates to computational physics simulations, and more particularly to a 3D CAD integrated simulation tool for plasma discharges. BACKGROUND

[0003] Current commercial tools for plasma simulation are generally optimized for high-density or high-temperature plasma environments, such as those found in tokamaks and fusion reactors. These tools often require overly complex meshes, leading to inefficient simulations and extensive computational resources. Moreover, they lack the integration of advanced solvers and a comprehensive 3D CAD environment, making it difficult to create, clean, simulate, and analyze plasma systems effectively. SUMMARY

[0004] This summary is provided to introduce a selection of concepts in a simplified form that are further described below in the detailed description. This summary is not intended to identify key features or essential features of the claimed subject matter, nor is it intended to be used as an aid in determining the scope of the claimed subject matter.

[0005] According to an aspect of the present disclosure, a non-transitory computer- readable storage medium is provided. The non-transitory computer-readable storage medium stores instructions executable by a processor to simulate plasma discharge. The instructions comprise a plasma solver for predicting the behavior of plasma discharge. The instructions comprise a 3D CAD system having an interface configured to receive user input relating to a physical environment. The instructions comprise a simulation 1 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 engine configured to simulate plasma discharge within the physical environment using the plasma solver.

[0006] According to another aspect of the present disclosure, a method of simulating plasma in a simulation software is provided. The method comprises defining a physical environment using a 3D CAD module. The method comprises enumerating a plurality of physical parameters and initial conditions corresponding to the physical environment. The method comprises predicting, using a multi-physics solver framework comprising a finite element method electromagnetic solver and a particle- in-cell solver, the behavior of plasma within the physical environment.

[0007] According to another aspect of the present disclosure, a plasma simulation system is provided. The plasma simulation system comprises a processor. The plasma simulation system comprises a memory comprising a non-transitory computer-readable storage medium, the memory being coupled to the processor. The plasma simulation system comprises a multi-physics solver framework stored in the memory and executable by the processor. The framework comprises a finite element method electromagnetic solver. The framework comprises a particle-in-cell kinetic solver. The framework comprises a fluid dynamics solver. The framework comprises a reaction solver. The plasma simulation system comprises an integrated 3D CAD environment stored in the memory and executable by the processor, the environment configured to enable a user to interact with the multi-physics solver framework.

[0008] The foregoing general description of the illustrative embodiments and the following detailed description thereof are merely exemplary aspects of the teachings of this disclosure and are not restrictive. BRIEF DESCRIPTION OF DRAWINGS

[0009] Non-limiting and non-exhaustive examples are described with reference to the following figures.

[0010] FIGS. 1-2 show flow charts corresponding to example computational workflows for a software system according to an embodiment of the present disclosure.

[0011] FIG.3 is an example CAD design created in the software system of the present disclosure according to an embodiment.

[0012] FIG. 4 shows an RF signal for voltage boundary conditions displayed in the software system of the present disclosure according to an embodiment. 2 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0013] FIG.5 shows post processing performed by the software system of the present disclosure according to an embodiment.

[0014] FIG.6 shows simulation results acquired via the software system of the present disclosure according to an embodiment. DETAILED DESCRIPTION

[0015] The present disclosure describes a 3D CAD integrated simulation tool for modeling plasma discharges and related electromagnetic phenomena. The system (e.g., the software, the software system, the simulation system, the simulator) integrates multiple computational physics solvers into a single software environment, including finite element method (FEM) electromagnetic solvers, particle-in-cell (PIC) kinetic solvers, fluid dynamics solvers, and reaction solvers. These operate with coordinated multi-time-step integration schemes. The software system may provide an integrated 3D CAD environment for enabling user interaction with the software system. The disclosed system provides specialized capabilities for low-temperature plasma applications, plasma-enhanced chemical vapor deposition (PECVD), spacecraft charging, vacuum discharges, and RF plasma processing.

[0016] The present disclosure progresses through aspects of the software, beginning with the Finite Element Method (FEM) for electromagnetic field calculations using Maxwell's equations, natural boundary conditions, and gauge constraints. The disclosure then details the Boundary Element Method (BEM) for spacecraft charging applications using analytic plasma environment models. Next, the Particle-In-Cell (PIC) solver is described, followed by the fluid solver. The disclosure then describes a reaction solver for modeling chemical interactions between plasma species. Lastly, example implementations of the software are provided. FINITE ELEMENT METHOD (FEM)

[0001] FEM use cases for simulating electromagnetic systems will now be described. Background information such as Maxwell's equations, natural boundary conditions, and gauge constraints is covered first, then how space is discretized in the FEM, and how this discretization is used to derive the electromagnetic governing equations which can be implemented in code. This section of the present disclosure concludes by describing 3 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 how simulation setup informs the FEM solver, and how the governing equations are solved.

[0002] The four well-known Maxwell-Heaviside equations, commonly called Maxwell's equations, are provided below: ^^ ⋅ ^^ ൌ^^ ^^^1a^ ^^ ⋅ ^^ ൌ 0 ^1b^^^^^ ^^^^

[0003] (1a) is Faraday's law; and (1d) iscalled Ampere's law. Maxwell's equations describe how the electric and magnetic fields are related to one another, as well as how they respond to field sources. Field sources refer to physical quantities which generate the electromagnetic fields. Charged particles are sources for the electric field, and currents - i.e. moving charges - are sources for the magnetic field.

[0004] While Maxwell's equations are useful, they do not provide how the sources for the electromagnetic fields change over time. This description is provided by the continuity equation for electromagnetism: ^^^^ ^^ ⋅ ^^ ൌ െ^^^^

[0005] The charge continuity be used to track how the electromagneticsources evolve over time. These sources can then be plugged into Gauss's law, (1a), and Ampere's law, (1d), to solve for the electromagnetic fields that develop. These equations can be used to provide a description of the electromagnetic fields and their sources in the bulk. However, understanding of how the fields behave at boundaries is needed in order to achieve a complete description of electromagnetic systems in general.

[0006] Understanding of how the electromagnetic fields behave at boundaries comes from Maxwell's equations. When Maxwell's equations in integral form are applied to solve for the fields at the boundary between two different materials, the following is found: ^^^^^^^ െ ^^ ^ଶ^^ଶ ൌ ^^^ ^∀^^ ∈ ^^Ω^

[0007] The subscriptis the free surface charge density, ^^^is the free surface current density, and ^^Ω denotes an arbitrary 4 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 boundary. These boundary conditions show that the electromagnetic fields can be discontinuous over a boundary if one of two conditions are present in the system: 1. The two materials have different electromagnetic properties ( ^^ or ^^ ) 2. A surface charge / current density exists on the boundary ( ^^^or ^^^, respectively)

[0008] Only one of the conditions needs to be satisfied in order to get discontinuous electromagnetic fields.

[0009] What makes these boundary conditions "natural" is the fact that they are naturally satisfied when solving Maxwell's equations. This means that they only need to be explicitly imposed on the domain boundaries. Any internal boundary would satisfy the natural boundary conditions, as the simulation system of the present disclosure solves Maxwell's equations on these internal boundaries. How the natural boundary conditions are imposed on the domain boundaries is described in later sections.

[0010] A description of the electromagnetic fields and their sources for an arbitrary system has now been provided. If the fields themselves were to be solved for, then the introductory discussion could stop here. However, the software system of the present disclosure solves for the electromagnetic potentials as opposed to the fields. This is discussed in greater detail later in the disclosure. In order to get the solution for the potentials, additional constraints on the system are required. These are the gauge constraints.

[0011] The uniqueness of a solution is important in solving electromagnetic problems. A solution is considered to be unique if it is the only possible answer for a given set of inputs. In electromagnetism, the electromagnetic fields are guaranteed to be unique. This means that for a given set of sources (charge and / or current density), the same solution for the electromagnetic fields is always obtained. Furthermore, if the electromagnetic fields were to be uniformly offset - as in, their values are increased / decreased by some constant value everywhere in the solution domain - then it would be found that these new offset fields correspond to a unique set of sources - as in they differ from the original set of sources.

[0012] While the solution for the electromagnetic fields is guaranteed to be unique, the same cannot be said for the solution for potentials. It is known that in quasistatics the electric field only depends on the gradient of the scalar potential: ^^ ൌ െ^^^^. This isultimately what leads to the solution of the electromagnetic potentials not being 5 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 guaranteed to be unique; the ability to get the same electromagnetic fields with different sets of potentials. While this is physically acceptable, it makes converging on a solution when using the FEM very difficult. A way of reducing the number of possible solutions is needed. According the methods of the present disclosure, this is accomplished through gauge fixing.

[0013] Gauge fixing is performed by imposing a gauge condition on the electromagnetic potentials. The gauge condition leaves the solution for the fields unchanged while constraining the possible solutions for the associated potentials. The simplest gauge is the Lorenz gauge: 1 ^^^^ ^^ ⋅ ^^ ^^^ଶ^^^^ൌ 0

[0014] The Lorenz gauge possible solutions, but does notguarantee uniqueness of the there exists a set of additional constraints that can be applied to the potentials. These additional constraints still obey the Lorenz gauge, but further reduce the number of possible solutions. The most common of these additional constraints is provided by the gauge condition for the Coulomb gauge: ^^ ⋅ ^^ ൌ 0

[0015] The Coulomb gauge guarantees uniqueness, and will be the default gauge for the remainder of the present disclosure.

[0016] The spatial domain is discretized in the FEM through interpretation of the mesh. The present disclosure is generally directed toward tetrahedral mesh elements interpreted in the FEM, though the principles apply to hexahedral elements as well.

[0017] The Finite Element Method is built around basis functions. These basis functions describe the spatial dependence of the nodes in the mesh. For example, consider a reference element having nodes P1=(0,0,0), P2=(1,0,0), P3 =(0,1,0), and P4=(0,0,1), which is defined in the coordinate system ( ^^, Γ, ^^ ). This is analogous tothe Cartesian coordinate system, ( ^^, ^^, ^^ ). This reference element is convenient as ithas well-defined extents, simplifying mathematical calculations. The finite elements in the mesh are related to this reference element.

[0018] The basis function for node ^^ in the reference element is given by the linear equation: ^^^^ ൌ ^^^ ^ ^^^^^ ^ ^^^Γ ^ ^^^^^6 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0019] Where ^^^ , ^^^ , ^^^, and ^^^ are coefficients which depend on the node number ^^.The values of the coefficients are determined by making use of the following property of the basis functions: ^^^^ ൫^^^൯ ൌ ^^^^

[0020] Where ^^^points from ^^ in the element, and ^^^^is the Kronecker delta which is equal to0 otherwise. By evaluating the fourbasis functions - one for each node in the element - at the four nodal positions, a system of 16 equations is obtained. This system of equations can then be solved to find that the coefficients are: ^^^^^^^^^^^^ 1 െ1 െ1 െ1^^ଶ^^ଶ^^ଶ^^ଶ 1 0 0൪

[0021] The nodal with respect to theelement. A nodal element; not the firstnode in the mesh - which would be called the global index. Any given node in the mesh could be a part of multiple elements, and the local index of this node is not the same for all of these elements. For example, a particular node may be the first node in some particular element, but it may be the third node in some other element. Each node will have as many basis functions defined for it as there are elements which contain said node. If a node is included in six elements, then there are six basis functions which describe the node - one for each element. This is why there is the superscript ^^ in the equation above; because the basis function - or more specifically, the basis function coefficients - will depend on the element being considered.

[0022] Properties of the basis functions will be discussed. The first property is one that has already been shown: ^^^^ ൫^^^൯ ൌ ^^^^

[0023] The basis functioncan be used to expand physical variables: ே ^^^^^^ ^

[0024] The sum on the rightover the nodes which belong to the element that contains the point, ^^, where the physical variable is to be evaluated. In this case the physical variable is the scalar potential. For tetrahedral elements, ^^ ൌ7 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 4. To better understand this equation both sides can be evaluated at the position of node ^^: ே ^^൫^^^൯ ൌ ^   ^^^^^^^ ൌ ^^^

[0025] where the RHS be zero for all nodes ^^ except when ^^ ൌ ^^, which onlythat the coefficients, ^^^, arethe values of the scalar potential at the nodes. This extends to all other physical variables as well. For instance the coefficients for the vector potential, ^^^, are also the values of the vector potential at the nodes.

[0026] The basis functions are used to expand physical variables because it allows the concept of the weak derivative to be taken advantage of. A particular concern with electromagnetic problems is discontinuities in the physical variables. If a boundary has a nonzero surface charge density or current density, then there will be a discontinuous jump in the electric and magnetic field values across the boundary. While this discontinuity is physical, it can prove to be a problem should differentiation of the physical variable over the boundary be needed. The discontinuity makes the local derivative undefined. This is where expansion of the physical variables becomes relevant.

[0027] By expanding the physical variables the spatial dependence of the variable is given to the basis functions. The resulting coefficients, for instance ^^^, no longer contain any spatial dependence. This allows the spatial derivatives of the physical variables to be expressed as: ே ^^^^^^^^ ^^^^^^^^^^^^^^

[0028] The coefficients - inany discontinuities that might be present in the physical variable. However, since the spatial derivatives are applied to the basis function instead of the coefficients, this discontinuity no longer creates problems. This is because the basis function is always differentiable, as can be verified using the linear equation shown above.

[0029] By using the basis functions to expand physical variables, all of the spatial dependence is given to the basis function; which is always differentiable. This allows differentiation of physical variables even when they are discontinuous - which would 8 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 otherwise make taking the derivative of the physical variables impossible. This fact is what makes the derivative a so-called weak derivative.

[0030] With an understanding of what basis functions are - along with their properties - gradients and integrals of the basis functions can be derived which are important when implementing Maxwell's equations into the Finite Element Method. The linear basis function defined above was formulated in the reference coordinate system, ( ^^, Γ, ^^ ).However, the necessary gradients and integrals need to be derived in the Cartesian coordinate system, ( ^^,^^, ^^ ), that is actually used by the solver.

[0031] Two approaches can be taken to obtain the necessary gradients and integrals in the Cartesian coordinate system. One approach is to define the basis functions directly in the Cartesian coordinate system. However, the preferable path is to use the Jacobian to transform the gradients and integrals from the reference to the Cartesian coordinate system. The Jacobian is defined as: ^^^^ ^^^^ ^^^^ é ù ^^^^ ^^Γ ^^

[0032] A positions of thenodes in the element. These positions are given by ^^^ , ^^^, and ^^^ in the equation above,where ^^ is the local node indices of element ^^.

[0033] The Jacobian's determinant is related to the volume of the element. Specifically, for a tetrahedral element: det^^^^^ ൌ 6^^^. Where ^^^ is the volume of the element, andthe factor of 6 comes from the fact that the Jacobian treats the nodes as if they form a cube with the same side lengths as the tetrahedron. It is well-known that the volume of a cube is 6 times larger than that of a tetrahedron with the same side lengths. This property is leveraged when deriving integrals of the basis function.

[0034] The Jacobian can be used to express the gradient of the basis function in the Cartesian coordinate system: ^^^^^ ் ି^^ ൌ ^^^^^ ൫^^^^̀^ ^ ^^^^^^ ^ ^^^^̀^൯

[0035] The quantityRHS is the gradient of the basis function in the reference coordinate system. While this gradient is given in terms of the reference coordinate axes, by multiplying by ^^^்^^ି^- which is the inverse of the 9 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 transposed Jacobian matrix - the axes can be transformed to those in the Cartesian coordinate system. In other words, the unit vectors ^^^,̀ ^^^ , ^̀^^ are replaced by ^^̀^, ^̀^, ^̀^^.

[0036] With the gradient of the basis function described, the integrals to be derived can be considered. Before doing so, how the volume of integration, ^^^^, can be transformed between the two coordinate systems will be examined. To do so, the integral is multiplied by det^^^^^.

[0037] The validity of this statement can be shown straightforwardly. To do so, an arbitrarily sized cube is considered. If ^^^^ were to be integrate^^^^d over the bounds of the cube, then ^^^would be obtained which is the volume of the cube in the Cartesian coordinate system. An equivalent reference cube is considered, which has a side length of 1 in all three components of the reference coordinate system. According to the statement made in the previous paragraph, the integral of ^^^^ would become: ^det^^^ ௧^^^^^^ . Where ^^^^௧ denotes that the integral is being performed in the referencecoordinate system. The determinant of the Jacobian does not depend on the coordinates of the reference coordinate system, and can thus be pulled out of the integral. Furthermore, the determinant is exactly equal to the volume of the cube, ^^^, since a tetrahedron is not being considered. Applying this to the integral leaves: ^^^^ ^^^^௧. Theintegral of ^^^^௧is equal to the volume of the cube in the reference coordinate system, which is 1. Therefore, just ^^^is obtained. This is clearly the same answer obtained when ^^^^ is directly integrated ^^^^over the volume in the Cartesian coordinate system.

[0038] Three integrals are of interest for derivation, and the reason for this will be seen in later sections. By using what was just covered, the three integrals of interest can be expressed as: ^^det^^^ ^ ^^^ ൌ ^ ൫^^^^^ ⋅ ^^^^^ ൯^^^^ ൌ ^൫^^^^^ ⋅ ^^^^^൯ ൌ ^^ ^^ ^ ^൫^^^^^ ⋅ ^^^^^^൯

[0039] the present disclosure as governing equations. As mentioned previously, the software system solves for the electromagnetic potentials, ^^ and ^^, rather than the fields, ^^ and ^^. The governing equation for the scalar potential is first considered, followed by the 10 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 vector potential. In both cases, a governing equation with the following matrix formalism is sought: ^^^^^^^ ൌ ^^^

[0040] Where ^^^^ is the mass matrix for node pair ^^^, ^^^,^^^ is the force matrix for node^^, and ^^^is the variable matrix - either ^^ or ^^ - for node ^^.

[0041] potentials are solved for instead of the fields because the solution for the potential converges faster than that for the fields. In light of the disclosure regarding the uniqueness of the solution, this might seem counterintuitive. In the case of having a nonunique solution it would be true that the electromagnetic potential solution might struggle to converge. However, this is not true if the uniqueness of the solution is guaranteed. Hence the importance of the gauge constraints as discussed above. Under these conditions, the solution for the potentials will converge faster than the solution of the fields due to ^^ being a scalar quantity as opposed to ^^, which is a vector quantity. It is easier to converge on the solution for a scalar value than for a vector quantity.

[0042] To obtain the governing equation for the scalar potential, Gauss's law is used as a starting point. This equation is then differentiated with respect to time, leading to a ^^^^ / ^^^^ term on the right-hand side. This term can then be re-expressed using the charge continuity equation. The equation is then re-expressed to be in terms of the electromagnetic potentials, after which the Coulomb gauge is applied by setting ^^ ⋅ ^^ ൌ0. This results in: ^^^^^ଶ^^^ ^^ ^ଶ^^^^ ^^ ^^^^^ ^ ^^ ൌ െ ^^^^^

[0043] Where ^^^^^ / ^^^^ ischarge density. This value comes from the current source, should one be used in the model. The current source can be thought of as adding charges into the system which were not originally there. This is why this term must be added into the governing equation; to account for these charges which do not originate from within the system.

[0044] To obtain the governing equation in a form that is similar to the matrix formalism, Galerkin's method is employed. This involves multiplying both sides of the governing equation by ^^^. Integration over the volume is then performed. This results in terms which involve ^ ^^ଶ^^^^^^. These terms can be re-expressed using integration11 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 by parts. The final step in applying Galerkin's method is to use the basis function expansion of the physical variables. This leads to: ^^^^ ^^^^ ^^^ ^ ^ ^^^^ ^^^^^ ⋅ ^^^^^ ^^^^ െ ∮ ^^^^^ ^ ^^^^^^ ൫ ^ ^൯ ^൬^^ ^^^^^ ^^^^^^^ ⋅ ^^^^ ൌ ^ ^,^^^ ^^^^^^^^^^^^ ^^ theformalism is now obtained. The only remaining step is to re-express the time derivative of ^^ using backwards differencing:^^ ^ ^ ^ ^ ^ ^^ ^^^ ^ ^^Δ^^^ ^^^^ ⋅ ^^^^ ^^^^ െ ∮ ^^ ^^^ ^ ^^Δ^^^^^^^ ⋅ ^^^^൫ ൯^ ^ ^ ^^ ^ ^ ^^^^^^

[0046] for,^ି^^^ is the solution from the previous timestep. A governing equation which can be^put into the form of the matrix formalism is now obtained. The governing equation for the vector potential is now considered.

[0047] The governing equation for the vector potential is obtained by starting with Ampere's law. The conduction current is then used to re-express the current density: ^^ ൌ ^^^^. The equation is then put in terms of the electromagnetic potentials. This leadsto a term on the LHS that appears as: ^^ ൈ ^^^ ൈ ^^^. This term is re-expressed using avector calculus identity. Finally, the Coulomb gauge is enforced by setting ^^ ⋅ ^^ ൌ 0.This results in: ଶ 1 ^^^^ ^^^^ ^^^^^^ଶ ^^ ^^ െ ^^ െ ^^ ൌ ^^^^^^ ^^^ ^^^^ ^^^^ ^^^^

[0048] Galerkin's was applied to the scalarpotential governing equation, which also leads to some surface integrals: ଶ ^^^^ ^^^^1 1^ ^^ ^ ^ ^^^ ⋅ ^^^^ ^ ^^ ^ ^^ ^^^^^^^^^^^^^^^^ ⋅ ^^^^ ^^^^ ^ ^^ ^ ^ ^ ^

[0049] timederivatives. Unlike with the scalar potential, the governing equation for the vector potential includes a second order time derivative. This can still be rewritten using backwards differencing: 12 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 Δ^^ଶΔଶ^^ ^^^^^ ^^^^^ ൫^^^^ ^^^ ⋅ ^^^^^ ൯ ^ ^^^ ^ ^^Δ^^^^^ ^^^^^ ^ ^^^^ െ ∮^^ ^^^ ^^^^^^^ ⋅ ^^^^^ ^^^^

[0050] seemthe governing equation for ^^ does not depend on the vector potential, both potentials can be solved independently of one another. This means that on a given timestep, ^^ is solved first, which can then be plugged into the RHS of the governing equation for ^^. Furthermore, backwards differencing is not used to express the time derivative of the scalar potential. This is because the time derivative can be calculated explicitly and plugged directly into the governing equation. Although the time derivative is still computed using backwards differencing, it need not be used to express the time derivative in the governing equation.

[0051] A governing equation for the vector potential, which can be put into the same form as equation the force matrix equation above, is now available. Before doing so, the surface integrals which appear in the governing equations will be examined more closely; specifically, how those surface integrals account for the natural boundary conditions discussed above.

[0052] The natural boundary conditions which arise from Maxwell's equations are given by the boundary condition equations. These natural boundary conditions are applied to the domain boundaries of the simulation. The domain boundaries are treated as being infinitely far away, which leads to ^^^ ൌ 0 and ^^^ ൌ ^^. Furthermore, both sidesof the domain boundary are treated as having the same material. Under these conditions, the natural boundary conditions simplify down to: ^^^^^ ൌ 0 ^∀^^ ∈ ^^Ω^1 ^^^^‖ ൌ ^^ ^∀^^ ∈ ^^Ω^

[0053] Before these boundary conditions can be applied to the governing equations, they need to be expressed in terms of the electromagnetic potentials. The quasistatic case is considered first. For such simulations there is no magnetic field, and so only the natural boundary condition for the electric field needs to be considered. To express this in terms of the scalar potential the constitutive relation is used: ^^ ൌ െ^^^^. Plugging thisinto the natural boundary condition on the electric field leads to: 13 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 ^^^^^ ൌ െ^^^^^^ ⋅ ^̀^ ൌ 0 ⇒ ^^^^ ⋅ ^̀^ ൌ 0 ^∀^^ ∈ ^^Ω^

[0054] The natural boundary conditions for a full-wave simulation will now be considered. In this case, both electromagnetic fields are present. Furthermore, the constitutive relation for the electric field changes to: ^^ ൌ െ^^^^ െ ^^^^ / ^^^^. Plugging thisinto the natural boundary condition on the electric field yields: ^^^^^^^^^ ⋅ ^̀^^^^^^^ ൌ െ^^ ൬^^^^ ^^^^^^ ⋅ ^̀^ ൌ 0 ⇒ ^^^^ ⋅ ^̀^ ൌ െ^^^^^∀^^ ∈ ^^Ω^

[0055] terms of theused: ^^ ൌ^^ ൈ ^^. Plugging this in gives:1 1 ^^‖ൌ ^^ ^^^^^ ൈ ^^^ ൈ ^̀^ ൌ ^^ ⇒ ^^^ ൈ ^^^ ൈ ^̀^ ൌ ^^ ^∀^^ ∈ ^^Ω^

[0056] One alongwith the natural is enforced in order to not excite the divergence of A, and is given by: ^^ ⋅ ^̀^ ൌ 0 ^∀^^ ∈ ^^Ω^

[0057] This additional constraint allows the natural boundary conditions to be rewritten as: ^^^^ ⋅ ^̀^ ൌ 0 ^∀^^ ∈ ^^Ω^^^^^ ⋅ ^̀^ ൌ ^^ ^∀^^ ∈ ^^Ω^

[0058] Having expressions for the natural boundary conditions in both the quasistatic and full-wave cases, the manner in which they are imposed onto the governing equations can be determined. In both the quasistatic and full-wave cases, the natural boundary condition for the scalar potential reads: ^^^^ ⋅ ^̀^ ൌ 0 ^∀^^ ∈ ^^Ω^

[0059] When examining the governing equation for the scalar potential, it can be observed that the surface integrals are proportional to ∮ ^^^^ ⋅ ^^^^. This expression canbe rewritten as ∮ ^^^^^ ⋅ ^̀^^^^^^^Ω^, where ^̀^ is the unit vector normal to the surface, ^^Ω,over which the surface integral is performed.

[0060] By expressing the surface integral in this manner, it becomes clear that the natural boundary condition implies that the surface integral is 0. The natural boundary condition is enforced by setting the surface integral to 0. This causes the governing equation for the scalar potential to read: 14 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 ^^^^ ^^^^^^^^ ^ ^^Δ^^^൫^^^^^ ^ ^,^^^^ ⋅ ^^^^^൯^^^^ ൌ ^ ^Δ^^^^ ^^^^^^^^ ^ ^^^^^^ି^൫^^^^^^ ⋅ ^^^^^^൯൨ ^^^^

[0061] ^^^^ ⋅ ^̀^ ൌ ^^ ^∀^^ ∈ ^^Ω^

[0062] The governing equation for the vector potential involves surface integrals which are propertional to ∮ ^^^^ ⋅ ^^^^. Just like with the scalar potential, this may be written as∮ ^^^^^ ⋅ ^̀^^^^^^^Ω^. To enforce the natural boundary condition for the vector potential,these surface integral terms are set to 0. This yields the following for the vector potential governing equation: ଶ ^^^^Δ^^ ^ ^ ^^൫^^^^^ ⋅ ^^^^^^ ^൯ ^ ^^^ ^ ^^Δ^^^^^^^^^^ ^^ ^^^^ ൌ ^ ^^^^൫2^^^ି^ െ ^^^ିଶ^ ^ ൯ ^ ^^Δ^^^^^^ି^൧^^^^^^ ^^gauge: ^^^^ ^^^^^^^ ^ ^^Δ^^^ ^^^^^ ⋅ ^^^^^ ^^^^ ൌ ^ Δ^^ ^,^^^^^^^ ^ ^^^^^ି^^^^^^ ⋅ ^^^^^ ^^^^൫ ൯ ^൫൯൨^ ^ ^ ^ ^ ^ ^ ^^^^^ ^ formof the matrix equation (reproduced below) will now be described: ^^ ^^ ൌ ^^^^ ^ ^

[0065] While the mass and force matrices are computed for the entire simulation domain, computation begins with computing them on an element-basis. The element- based matrices can then be summed to get the global matrices. To see how this is done, ^ the mass matrix is considered first. The element mass matrix, ^^ , is given by the LHS^^of the governing equation (after pulling ^^ out):^^ ^^^^ൌ ^ ^^^ ^ ^^Δ^^^ ^^^^ ⋅ ^^^^൫ ൯^^^^^^ ^ ^

[0066] By summing the element-based contributions to the mass matrix, the global mass matrix can be obtained: ^^^  ^^^^^ ^^15 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0067] Attention is now be shifted to the force matrix. The element force matrix, ^^^^, is given by summing the RHS of the governing equation over ^^ : ே ^^^^ ^^^ ^,^ ൌ ^   ^ ^Δ^^ ^^^^^^^ ^ ^^^^^ି^൫^^^^^^⋅ ^^^^^^^^^ ^ ^ ^ ^ ൯൨ ^^^^

[0068] The force matrix, just like with^^^ ൌ ^  ^^^^^

[0069] In practice, FEM solvers over the element-based matrices and instead directly add contributionsmatrices. This is done as it decreases overhead. In certain embodiments, the software system of the present disclosure skips over the element-based matrices and instead directly adds contributions to the global matrices.

[0070] Just like with the scalar potential, the mass and force matrices for the vector potential can first be computed on an element-basis. The element mass matrix for the vector potential is given by the LHS of the governing equation: ଶ ^ Δ^^ ^^^ ൌ ^ ^൫^^^^^ ⋅ ^^^^^൯ ^ ^^^ ^ ^^Δ^^^^^^^^^^^^^ ^ ^ ^ ^ ^^^^

[0071] while the of the governingequation over ^^ : ே ^^^^ ^^^^ ^   ^ ^^^ିଶ ^ Δ^^ଶ^ ^ ^ ^ ^^^^^^ ^ ^^ ^ ^^^^^^^^^^ ^^^^^^^^ ൌ ^   ^^ ^ ^^^

[0073] It can be noticed that thevector potential is a vector quantity. This comes from the fact that the vector potential itself is a vector quantity, and the components of the force matrix can be interpreted as the source for that component of the vector potential. Meaning, to get a solution for ^^௫, the x-component of the force matrix is used. This also allows for the components of ^^ to be solved independently of one another, which helps with solution convergence. Just like with the scalar potential, 16 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 the element-based matrices for the vector potential can also be skipped over and the global matrices can instead be computed directly.

[0074] In order to obtain a solution for the first timestep, the initial conditions for the model must be known. This is because the solution for the potentials depends on the values from the previous timestep. When solving for the solution of the very first timestep, values must be provided for the previous timestep. These values are the initial conditions. The initial conditions are defined as follows: ^^^∀^^ ∈ Ω, ^^^^ ൌ 0^^^^ ^^^^^∀^^ ∈ Ω, ^^^^ ൌ 0^^

[0075] To enforce in the force matrices areset to 0 at the very the arrays which store the previous timestep solutions to 0.

[0076] Boundary conditions other than natural boundary conditions may be considered. The types of boundary conditions include two basic types: Dirichlet-like and Neumann- like. Dirichlet-like boundary conditions involve setting a quantity to some value on a surface: ^^ ൌ ^^ ^∀^^ ∈ ^^Ω^

[0077] where ^^ is the physical variable to which the boundary condition is applied and could be ^^ or ^^, and ^^ is the value to which the physical variable is set.

[0078] Neumann-like boundary conditions set the normal derivative of a quantity to some value on a surface: ^^^^ ^^^^ൌ ^^ ^∀^^ ∈ ^^Ω^

[0079] Where ^^^^ / ^^^^derivative, whose exact expression depends on whether ^^ is a scalar or vector quantity.

[0080] Dirichlet-like boundary conditions are defined by the user whereas the Neumann-like boundary conditions are used to enforce the natural boundary conditions. Dirichlet-like boundary conditions may be utilized to set the scalar potential to some value on a surface. The same may be done with the vector potential. When applied to the vector potential, one, two, or all three components of ^^ may be set to some value on a surface.

[0081] Sources may be incorporated into the solver. The simplest source to consider is a voltage source. A voltage source is defined in a model by utilizing the scalar potential 17 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 boundary condition. The scalar potential boundary condition allows for both time- constant and time-varying definitions. The time-varying boundary condition is commonly referred to as the voltage source. Since the voltage source is defined by using boundary conditions, it is also enforced through the boundary conditions.

[0082] Another source that may be used is the current source. Unlike the voltage source, the current source definition does not use a boundary condition. The current source is enforced when computing contributions to the force matrix. This is accomplished by looping through the nodes with a current source defined on them, obtaining the corresponding global node index, and using all of that information to obtain ^^^^^ / ^^^^. The current source waveform is provided in units of Amps. The software system converts this to ^^^^^ / ^^^^ by using the volumes of the elements which contain the nodes which have a current source defined on them.

[0083] Additional sources to consider are Particle-In-Cell (PIC) and fluid sources. These sources use the charged particle densities which are simulated in PIC and / or fluid to act as a source for the electromagnetic potentials. The software system uses these charge densities to compute the ^^^^^ / ^^^^ value that arises from the PIC and / or fluid environments. This value is then added to the force matrix in the same manner as accomplished for the current source.

[0084] The matrix equation is solved by the software system. A particular concern with solving the matrix equation is the size of the matrices. Specifically, the size of the mass matrices may prove to be computationally expensive. The mass matrices have size (N ൈ N) where N is the number of nodes in the simulation domain. Advantageously,the mass matrices are sparse. This is a result of the fact that contributions are only computed for nodes that neighbor each other, meaning the nodes are a part of the same element. The fact that the mass matrices are sparse may be utilized to greatly reduce the computational requirements of solving the matrix equation. To exploit this fact, the BiConjugate Gradient Stabilized method, or BiCGStab as it is commonly called, is employed in the software for improved performance. BOUNDARY ELEMENT METHOD (BEM)

[0085] Now that the mathematical framework for embodiments of the software system has been established, the boundary element method (BEM) will be described in greater detail. 18 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0086] The software system has been recognized as particularly affective in simulating spacecraft charging (e.g., satellite, rocket, etc.). This refers to the charging the spacecraft experiences in the presence of the space-based plasma which is trapped by Earth's magnetic field lines. Such simulations involve representing the spacecraft using a surface mesh which includes all outer-most surfaces of the spacecraft. In general, spacecraft charging simulations focus on the potential and normal electric field which develops on the spacecraft. The end goal is to determine what the risk of a vacuum discharge occurring on the craft is.

[0087] The software system is configured to simulate (e.g., calculate) the potential and normal electric field on the spacecraft in two ways. One of these methods involves meshing the volume which surrounds the spacecraft, and using the FEM - as described above - to simulate the potential which arises from the surrounding plasma. The other way foregoes the volume mesh and instead directly solves for the potential and normal electric field on the surface of the spacecraft using analytic plasma environment definitions. This second method is referred to as the Boundary Element Method, or BEM. The general idea is to relate both the potential and normal electric field to the surface charge density on the surface of the spacecraft. These relations can then be used to directly relate the potential to the normal electric field. Finally, a charging equation can be constructed which relates the change in potential to a change in the surface charge density.

[0088] Unlike with the FEM, the BEM uses element-based quantities. To denote this difference, capitalized subscripts are used to represent element-based quantities.

[0089] The equation which relates the potential to the surface charge density stems from the well-known equation for the potential from a point charge: ^^ ൌ1 ^^ 1 ൌ^^^ ^^^^ 4^^^^^^^ 4^^^^^^^

[0090] The potential at an element, ^^ூ, is acquired based off of the surface charge density on all other elements, ^^^, and the distance between the elements, ^^ூ^. The distance is determined based off of the centers of the elements. These variables can be substituted into the equation above and summed over ^^ to get ^^ூthat results from the surface density on all elements. This can then be written as aequation: 1 ^^ ^^^   ^ ^ ^^^^^19 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0091] This equation does not include any sort of basis function. This is because the basis functions are only needed for nodal-based quantities. More explicitly, the basis functions are needed to describe the nodes because nodes will be shared between multiple elements, and a way of accounting for that is needed. There is only one way to represent the element so basis functions are thus not needed.

[0092] Now that an equation which relates the potential to the surface charge density has been established, the normal electric field can be addressed. The equation which relates the normal electric field to the surface charge density is derived from the equation for the electric field from a point charge: 1 ^^ 1 ^^ ^^ ൌ െ ^̀^ ൌ െ ^^̀^^^^^ 4^^^^^^^ଶ4^^^^^^^ଶ

[0093] To obtain the surface, the dot product ofboth sides with ^̀^ is for the potential equation is then applied to yield: ^^ ^^ െ ^^ ⋅ ^̀^^^^ ⋅ ^̀^^ ൌ െ1  ^ ^൫ ூ ^൯ூ ^ ^^^^ ൌ ^^ ^^^^ ^^ଷ ^ ூ^ ^

[0094] Where on the right-handside. Having obtained these two equations, the potential can be related to the normal electric field directly. Relating the equations for ^^ and ^^^ ⋅ ^̀^^ begins with expressingூ ூ^^ in terms of the potential:^ି^^^ ൌ ^^^^^ ூூ^

[0095] This can then be substituted for the normal electric field. Doingso results in a change of subscript variables: ି^^^^ ⋅ ^̀^^ ൌ ^^ ^^^^ூ ூ^ ^^^

[0096] Having related the normal electric field to the potential, an equation is needed which describes how the potential responds to a change in the surface charge density. Such an equation is commonly referred to as a charging equation. To obtain the charging equation, an expression for the total surface charge density on an element must first be defined. The total surface charge density includes the free charges which are sources for the normal electric field, as well as the bound charges which arise from capacitive effects. This expression is given as: 1்^^^^^ ⋅ ^   ^20 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0097] Where ^^ூ^is a capacitance matrix. Substituting the equation for the normal electric field allows the right-hand side to be expressed purely in terms of the potential: 1 ^^்ூ ൌ ^^ூ^^^^ି^^^^^^ ^^^   ^^ூ^൫^^ூ െ ^^^൯ ൌ ^^ூ^^^^ூ^

[0098] Where to time and using backwardshand side results in: ^^^^்^ூ ^^^ ^ ൌ^ ൌ ^ ^ െ ^^^ ^ ^ ^ ^ூ ூ^ ^ ^ → ^^ூ Δ^^ ൌ ^^ூ^൫^^^ െ ^^^ି^൯ ^^^^ Δ^^^

[0099] with the surface^^^^on the left- hand side, yielding the final charging equation: ^^^ ^ூ^^^^ ൌ ^^ூ Δ^^ ^ ^^ூ^^^^^ି^

[0100] This equation can potential on the surface of thespacecraft over time. are are currents that arise from the plasma environment, transport, and yield quantities.

[0101] Three broad sources of the surface current density are accounted for. These current sources can be accounted for by expressing the total surface current density as: ^^ூ ൌ ^^ூ,^^௩൭1 െ ^    ^^^^^^ ^ ^  ^^ூ,^

[0102] where ^^ூ,^^௩isplasma;^^^^^is related to the yield quantity, and the sum over ^^ is performed to account for all yield quantities; and ^^ூ,^is the current density that arises from transport, and the sum over ^^ is performed to account for all transport quantities.

[0103] Expressions for each individual current source term are now sought and will be described. The current that arises due to the plasma environment will now be discussed. At a high-level, this current can be expressed as: ^^ூ,^^௩ ൌ ^^^^ூ|^^ூ|

[0104] Where ^^ூis the plasma density on element ^^, and|^^ூ|is the magnitude of the plasma velocity on element ^^. The problem with this equation is the flux quantity, ^^ூ|^^ூ|, as it implies that all of the particles which comprise ^^ூhave the same velocity, |^^ூ|. This is not the case. In reality, there is a large spread in the particle velocities. To account for this, the Maxwellian distribution can be used: 21 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 ^^^^^,^^^ ൌ4^^ 1 ^^ ^ ^^^^ ^^ ^^ ^^exp൬െ^^^ ^^^^^

[0105] If the then the probability of finding a particle withwith an average energy - aka temperature - of ^^^^^ can be obtained. While this equation is very useful, the flux still needs to be known; not what the probability is of finding a particle at some specific energy. Fortunately the Maxwellian distribution may be modified in order to get this information. All that needs to be done is multiply it by the flux. However, since the Maxwellian distribution is given as a function of energy the expression for the plasma flux must be converted to be a function of energy as well. After doing so, the following is provided:: ^^^^^,^^^ ൌ ^8 ^^ ^^ ^^exp ൬െ ^

[0106] If this plasma flux may beobtained while

[0107] An expression for ^^ூ,^^௩is nearly ready to be defined. The only thing that is left to do is consider how the plasma particles' energies would interact with the potential on the surface of the spacecraft. As a plasma particle approaches the spacecraft, it will experience an acceleration due to the potential on the surface of the craft. This modifies the energy of the plasma particle, which must be accounted for when ^^^^^,^^^ is integrated over energy. To do so, ^^ூ,^^௩can be expressed as: ^ ^^ூ,^^௩^^^ூ ,^^^ ൌ ^^^  ^^ ^^ ^ ^^^^ ^^^^^ ^ ^^^^ூ ,^^^^^^^

[0108] wherethe acceleration the spacecraft potentials, and the fraction inside the integral is used to renormalize ^^^^^,^^^ to the updated energy, ^^ ^ ^^^^ூ.

[0109] The lower bound of integration, ^^, is chosen based off of whether the plasma particle will be repelled or attracted by the potential. If the particle will be repelled by the potential then ^^ ൌ |^^^^|. This is because if a particle has an energy less than thisvalue, then it will lose all of its energy due to the deceleration from the spacecraft potential, and is thus unable to reach the spacecraft. Hence, the particle will not contribute to the flux on the surface of the spacecraft. However, if the particle will be attracted by the potential then ^^ ൌ 0, because the particle will always impinge on the22 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 surface of the spacecraft in this case. While the expression inside of the integral is complex, it can actually be integrated analytically. Doing so gives an expression for ^^ூ,^^௩which can be easily implemented. The yield quantities will now be considered.

[0110] When a plasma particle impacts the surface of the craft, several interactions may occur. The particle can induce ionization, producing an additional electron, or the particle can be elastically scattered off the surface back into the plasma. All of these interactions can be described by a current, referred to as the yield quantity. The software system of the present disclosure is configured to track three yield quantities: 1. Secondary Electron Yield (SEY) 2. Proton Electron Yield (PEY) 3. Backscattered Electron Yield (BEY)

[0111] SEY is characterized by an electron impacting the spacecraft surface and inducing ionization. PEY describes the ionization that results from an ion impinging on the surface of the spacecraft. If an incident electron does not induce ionization, but is instead scattered back into the plasma, then that electron is a product of ^^^^^^. The probability of any particular yield event occurring is based on the material of a surface element, as well as the energy of the plasma particle that is incident on the element. These dependencies are accounted for by the variable ^^^^^^, which is a user-defined material dependent property. Multiplying this with the expression for the plasma flux allows the current density due to yield quantities to be expressed as: ^ ^^ ,^^,^^^ ^^^  ^^ ^^^^^ ^ ,^^^^^^^^

[0112] as before. To obtain this expression in the context of a previously described surface current equation, division by ^^ூ,^^௩is performed: ^ ^   ^^^^^^^^ ^^ ^^ ^ ^^^^ ^^^^^ ^ ^^^^ ,^^^^^^^^

[0113] Theby the spacecraft potential has been discussed above, as well as how, if the particle's energy is not high enough, it will not contribute to the plasma environment current. The same consideration can be applied to the secondary electrons - those produced by SEY and / or 23 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 PEY - but in the opposite direction. If the potential on a surface element is positive, then the secondary electrons will be attracted to the element. This means that for such an element, there is a chance that the secondary electrons get stuck within the surface as they are attracted by the potential. These secondary electrons would not contribute to the yield quantities. To account for these secondary electrons that do not manage to escape the surface of the spacecraft when the surface potential is positive, ^^^can be multipl^^^ied by a term that takes the form: exp ^െ^^^^ூ / ^^^^^^. The secondary electrons are assumed to have a temperature of ^^^^^ ൌ 2eV. Plugging this into the correctionfactor gives: ^^ exp ൬െூ^ 2

[0114] By multiplying ^^^by this the secondary electrons that areunable to escape surfaces with are accounted for.

[0115] The final sources of surface current density are transport-based. There are three transport-based contributions to the surface current density that are accounted for: 1. Photoemission 2. Conduction 3. High energy particles

[0116] When the spacecraft is being illuminated by the Sun, the solar light can induce currents on conductors via the photoelectric effect. The electrons that are produced by this process and manage to escape the surface of the spacecraft are referred to as photoemission, or photocurrent. Different materials experience different photocurrents for the same intensity of solar illumination. This may be characterized by a user-defined material property which defines the photocurrent under normal, full-intensity illumination. This value is then linearly reduced according to the user-defined illumination intensity and angle.

[0117] The illumination intensity is based on how far the spacecraft is from the sun. For any Earthbound orbits, an illumination intensity of 1.0 can be used. This means the spacecraft will roughly experience the same illumination as on Earth. The argument here being that the distance between Earth and the orbit of the Earth-bound spacecraft is negligible compared to the distance between the Earth and the Sun, which is true. However, the intensity can be reduced if the element is in the shadow of some other 24 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 feature on the surface of the spacecraft. The software system of the present disclosure is configured to account for such shadowing.

[0118] The illumination angle is the angle between the unit vector normal to a surface element and the unit vector which describes the illumination direction. The smaller the illumination angle is, the smaller the photocurrent will be. The reason for this is subtle, but comes from the fact that the intensity of the illumination is spread out over a larger area when it is incident on the spacecraft at an angle.

[0119] Something to consider with photoemission is the ability, or lack thereof, of a produced electron to escape the surface of the spacecraft. In the case of the surface having a positive potential, the electron will be attracted to the surface. This reduces the photoemission current. To account for this, the same correction factor that was used for the secondary electrons is used.

[0120] Currents which arise due to the plasma, interactions between the plasma and spacecraft, and photoemission have now been discussed. The current that would flow between elements also needs to be accounted for. This current is commonly the conduction current. There are two sources of the conduction current which are of concern: current that flows from an insulating material on the surface of the spacecraft to the underlying metallic chassis; and the current that flows along the surface of the spacecraft between elements. The current that flows from insulating materials to the underlying chassis is described by: Δ^^ ^^^^ൌ ^^^^

[0121] Where Δ^^ is the potentialthe insulating surface element and the chassis; ^^ is the resistivity of the insulating material; and ^^ is the thickness of the insulating material. ^^ and ^^ are user-defined quantities for the material.

[0122] To obtain the conduction current that flows between surface elements, the conductance on the edge that is shared between a pair of surface elements is first determined. Any given pair of neighboring surface elements shares only one edge. The resistance of the edge is based on the surface resistivity of the two surface elements along with geometric quantities. The geometric variables are used to convert the surface resistivity to a resistance. The contributions to the edge's resistance from the two surface elements are: ^^ ^ ^^ூ^ ^^^^^^ௌ,^25 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0123] Where ^^ூand ^^^are distances from surface element centers to the center of the edge; ^^ூ^is the length of the edge; and ^^ௌ,ூand ^^ௌ,^are the surface resistivities of the two surface elements.

[0124] The conductance of the edge is given by the inverse of the sum of the two resistances given above: 1 ^^ூ^ൌ ^^ூ ^ ^^^

[0125] ^^ூ^can be multiplied with the potential difference between surface elements ^^ and ^^ to obtain the current that flows between the two elements. To obtain the current density on one of the elements, the current can be divided by the area of the surface element. By summing over ^^, the current that flows between all neighboring elements and element ^^ can be accounted for: ^^ Δ^^ ^^^ଶ,ூ ൌ ^   ூ^ ூ^^ ^^ூ

[0126] The final transport contribution to the current density comes from high energy particle sources. There are two kinds of high energy particle sources in the software system of the present disclosure: Geant4 and plumes. For both kinds, the contribution to the current density is computed in substantially the same way as the plasma environment contribution, ^^ூ,^^௩is computed. The general approach involves using user-defined quantities to obtain the particle flux on the surface of the spacecraft. This flux is then multiplied by the charge of the plasma species to obtain the contribution to the current density.

[0127] The first matrix to be considered is the charging matrix, ^^ூ^. The construction of this matrix proceeds by first looping through all elements in the mesh to obtain element ^^. The solver then performs another loop to obtain element ^^. What elements are included in this second loop depends on what the user specifies for the solver to include.

[0128] By default this second loop includes all elements which are located within two elements of element ^^. The default behavior is to include the elements which neighbor element ^^, as well as the elements which neighbor those neighboring elements. This is called the ^^^^^^ order nearest neighbor approximation

[0129] A user could alternatively choose to only include the elements which directly neighbor element ^^ - the nearest neighbor approximation. There is also the option to 26 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 choose the self-capacitance only approximation. In this configuration, the solver uses the nearest neighbor approximation, but only the diagonal terms of the capacitance matrix are computed. Only ^^ூூis computed when assembling ^^ூ^. The final method to consider is the full matrix method. In this configuration, the second loop is performed over all elements in the mesh.

[0130] While constructing the charging matrix, the solver also computes the contributions to the current density according to the equations provided above.

[0131] After performing the matrix assembly, the charging matrix equation can be constructed by using the charging matrix, ^^ூ^; the previous solution, ^^^^ି^; and the surface currents, ^^ூ. By default, the matrix equation is directly a usermay specify (e.g., via the user interface) that the software system use solver. The iterative solver is particularly helpful when the solver does not use the full matrix method when assembling ^^ூ^. In this case, the charging matrix will be sparse, which can be taken advantage of by the iterative solver to arrive at a solution faster than the direct solver can.

[0132] After solving for the potentials, the normal electric fields can then be obtained. The assembly of the matrices for this matrix equation is performed in the same way as was done to assemble the matrices in the charging equation. ANALYTIC PLASMA ENVIRONMENTS

[0133] The BEM solver utilizes analytic plasma environments. Such plasma environments are not directly simulated, but rather treated as a charge density boundary condition. The analytic plasma environments are derived from established computational methods, and will now be discussed in greater detail.

[0134] Three analytic plasma environments can be defined in a surface charging (BEM) simulation: 1. Laplace 2. Nonlinear 3. Barometric

[0135] Each analytic model estimates the space charge density of the plasma near the spacecraft. The manner in which each analytic model estimates the plasma space charge density will be described. 27 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0136] The Laplace environment assumes that the space charge density near the spacecraft is zero. This does not indicate an absence of plasma, but rather represents an approximation wherein the plasma space charge density is negligible compared to the charge density on the surface of the spacecraft. The space charge density for the environment is given as: ^^ ^^ൌ 0^

[0137] The Nonlinear environment assume that space charge density near the spacecraft is zero. Instead, the spacenear the craft is given by: é ù ^^ ^^^^êmax^1,^^^^^, ‖^^‖^^ú ú ú û

[0138] Where ^^ is the the plasma wake; ^^^^^, ‖^^‖^is the convergence plasma: ^^^^ ^^ ^ ^^ ^ ^^

[0139] When a spacecraft the plasma density is reducedbehind the spacecraft. This region of depleted plasma is known as the plasma wake, and can be considered analogous to the wake that a boat leaves behind as it travels through water. The factor by which the plasma density is reduced in the wake is given by ^^.

[0140] The convergence factor, ^^^^^, ‖^^‖^, accounts for how charged particles can beaccelerated in the presence of high potentials. Without the convergence factor, a small timestep size would be required in order to accurately resolve the motion of the plasma. By including the convergence factor, the timestep size need not be constrained as much, naturally resulting in faster simulations.

[0141] The final analytic environment to be considered is the Barometric environment. This environment also does not assume that the space charge density is zero near the spacecraft, like the Nonlinear environment. However, unlike the Nonlinear environment, the Barometric environment provides separate equations for the space charge density of the ions and electrons. These space charge densities are given as: 28 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 ^^^ ൌ ^^^^^^ì ï1 ^ ^^ െ ^^^^^ln ^^^^, if ^^ ^ ^^^^^ln ^^^^

[0142] Where space charge density. Ascharge density of theൌ ^^^^, reduced by the wake factor, ^^. It might appear that the ionspace charge density is constant in time, but that is only if the wake factor does not change over time, which is not strictly guaranteed.

[0143] As was described above, the Laplace environment assumes that the plasma space charge density near the spacecraft is negligible compared to the charge density on the surface of the spacecraft. As such, the Laplace environment works best when describing plasmas with low densities. It is most commonly used to describe the plasma environment found in a geosynchronous environment, as the plasma is sparse this far from the Earth.

[0144] Both the Barometric and Nonlinear models work well with plasma environments that have a high density, such as the ones that can be found in most low- Earth orbits. The primary difference between these models is what plasma temperature they work well with. If the plasma temperature is below the potentials on the surface of the spacecraft, then the Nonlinear model will work best. This is due to the convergence factor that is included in the Nonlinear model, but left out of the Barometric model. Conversely, when the plasma temperature is comparable to or above the potentials on the spacecraft surface, then the Barometric model works best. Under these conditions the convergence factor can be ignored, which would also allow for the simulation to run faster. The Barometric model is better under these conditions because it does not include the convergence factor.

[0145] Generally speaking, the Nonlinear model is optimized for use in spacecraft charging simulations in low-Earth orbits under eclipse conditions where the spacecraft is not being illuminated by the Sun. Under these conditions, the surface potentials on the spacecraft could reach large values. This necessitates the need to include the convergence factor which is excluded from the Barometric model. If the spacecraft is under illumination, however, then it is advisable to use the Barometric model. This is because the photoemission will counteract the other sources for the surface potential, leading to a potential which is comparable to the plasma temperature. Both the 29 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 Nonlinear and Barometric models can simulate plasma wakes, as can be inferred from the inclusion of the wake factor, ^^, in the equations for both models. PARTICLE-IN-CELL (PIC) SOLVER

[0146] While the analytic plasma environments may be used to simulate space-based plasmas, they are not equally affective for simulating all plasmas in general. In cases wherein the analytic plasma environments may not be affective- most notable of such cases are discharge and plasma enhanced chemical vapor deposition (PE-CVD) models - another way to model plasmas is needed.

[0147] To address the limitations of the analytic plasma environments, the software system includes a particle-in-cell (PIC) solver. The PIC solver directly models the plasma using what is known as the kinetic theory. Kinetic theory attempts to model plasmas by explicitly simulating the plasma particles. However, it is challenging to simulate each individual particle in a plasma. If the software system were to do so, then the time required to simulate just a single timestep would be prohibitively long.

[0148] To address this issue, kinetic theory utilizes so-called macroparticles. A macroparticle is defined as a group of particles which belong to the same species. The fact that the macroparticle only groups together like particles is important, and will be described in greater detail later in the present disclosure. By clumping particles together into macroparticles, all particles in the system can be effectively described without needing to model them individually. This reduces the number of particles which need to be explicitly simulated, as they are replaced with the macroparticles, which are far fewer in number. This enables the simulation of a plasma to run on a timescale that could be on the order of minutes, hours, or days.

[0149] The PIC solver is a mesh-less solver. That is, the PIC solver does not require a mesh to run by itself. Instead, the macroparticles can be thought of as being equivalent to mesh nodes. This presents some challenges, particularly with respect to obtaining the electromagnetic fields that a macroparticle experiences.

[0150] For each macroparticle, the software system is configured to track the position, velocity, and number count - the number of particles contained within the macroparticle. The equation used to govern how the velocity changes over time will be described first.

[0151] The first step to understanding how the velocity changes over time is to determine what the acceleration is that the macroparticle experiences. To obtain this 30 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 quantity, the forces that the macroparticle is subject to must be known. As charged particles, the macroparticles are subject to the Lorentz force. By dividing the Lorentz force law on both sides by the particle mass, the acceleration can be obtained: ^^ ^^^ൌ ^^^^^ ^ ^^^ ൈ ^^^^

[0152] There are three quantities on the RHS which depend on particle quantities: charge, ^^; mass, ^^; and velocity, ^^. Strictly speaking, these quantities belong to the particle and not the macroparticle. However, in practice there is no difference. Because the macroparticles only group together particles of the same type, the charge-mass ratio, ^^ / ^^, remains unchanged.

[0153] To see how the macroparticle velocity is the same as the particle velocity, it can be imagined that 100 electrons start at rest. If they were placed inside of a uniform field, then the acceleration they would experience would be the same for every electron. This would then lead to every electron possessing the same velocity.

[0154] In reality, the electromagnetic fields are very rarely uniform. In this case, it can be said that all of the particles contained within a macroparticle are located at approximately the same position. In which case they would all experience approximately the same acceleration. This means that the macroparticle velocity can be said to be approximately the average velocity of all particles contained within the macroparticle. And so the acceleration obtained from the Lorentz force law is the average acceleration experienced by all particles within the macroparticle.

[0155] Now that the macroparticle acceleration is known, it can be used to update the macroparticle velocity and position using well-known kinematic equations: 1 ^^^ ൌ ^^^ ^ ^^^Δ^^ ^^^ ൌ ^^^ ^ ^^^Δ^^ ^^^^Δ^^ଶ

[0156] The fact thattime difference, Δ^^, is what makes the PIC solver explicit in time.

[0157] The PIC solver is a mesh-less solver. This presents a challenge in that electromagnetic field data on the nodes of the FEM mesh must be used to determine the field strength that a macroparticle experiences. This interpolation is performed through a C 0 scheme. This involves determining what element a given macroparticle is in and computing weights - which range from 0 to 1 - for each node in the element. These weights correspond to how much influence a particular node has on the electromagnetic field strength at the macroparticle, and are determined based off of how 31 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 close the macroparticle is to the node. These weights are then utilized according to the C 0 interpolation scheme which defines a linear interpolation: ^^ெ^ ൌ ^^^^^^ ^ ^^^^^^ ^ ^^ଶ^^ଶ ^ ^^ଷ^^ଷ

[0158] Where ^^ெ^is the field strength at the macroparticle; ^^^is the weight for local node ^^; and ^^^is the field strength at local node ^^. The field that is plugged into ^^ could be either the electric or magnetic field.

[0159] When a user defines a PIC environment, an initial density and temperature are provided. These quantities are used along with the PIC statistic to determine the initial macroparticle quantities.

[0160] During initialization, the number of macroparticles to spawn for each environment in each element is first determined. This is done by taking the user-defined PIC statistic and dividing by the number of PIC environments in the model. This resulting quantity is then divided again by the number of mesh elements in the model. With the number of initial macroparticles known, the macroparticles can be spawned. There are three quantities which must be known to the software when doing so: the position, velocity, and number count. These macroparticle-based quantities are initialized so that when interpolated onto the nodes of the FEM mesh the following conditions are obtained: ^^^∀^^ ∈ Ω, ^^^^ ൌ ^^^ ^⃐^^ ൌ ^^^ ^^^∀^^ ∈ Ω, ^^^^ ൌ ^^

[0161] Where ^^^is the user-defined initial density, ^^^is the user-defined initial temperature, and ^^⃐^ denotes the temperature averaged over the entire model domain. It is also worth noting that the third condition does not mean that the temperature will be zero.

[0162] To initialize the macroparticles, the mesh elements in the model are looped through. The volume of the mesh element is then obtained and multiplied with the user- defined initial density. This determines how many individual particles there will be in the mesh element. This quantity is then divided by the number of macroparticles per element for the environment of interest to determine how many individual particles are contained within each macroparticle. This quantity is used to initialize the number count of each macroparticle that will be spawned into this element, which ensures that the first initial condition will be satisfied. 32 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0163] To obtain the initial positions of the macroparticles, the space inside of the mesh element is randomly sampled. This is done for each macroparticle within the element, so that each macroparticle is given a unique position.

[0164] When the macroparticle velocities are initialized, a Maxwellian distribution is first constructed using the user-defined initial temperature. This distribution is then randomly sampled to obtain an initial velocity. This velocity is then applied to a pair of macroparticles. One macroparticle in this pair is given the velocity directly; while the other macroparticle is given the initial velocity multiplied by -1. This is done so that the total velocity of the macroparticle pair is ^^. By ensuring that the number of macroparticle per element for each PIC environment is even, this method of initializing the macroparticle velocities ensures that the third initial condition is satisfied.

[0165] Because of how the macroparticle velocities are initialized, an initial temperature of ^^^on the nodes of the FEM mesh is not guaranteed. However, the average temperature across the whole model would be roughly equal to the user-defined initial temperature.

[0166] The software system supports the use of Dirichlet-like boundary conditions in the PIC solver. Specifically, for example, this includes: 1. Reflect / replace 2. Density 3. Velocity

[0167] The reflect / replace boundary condition is the most basic boundary condition. When defining a PIC environment, a user is enabled to specify how the solver should treat the macroparticles which impact the boundaries of the body the environment was defined for. The user is given at least two settings for this: the boundary type, and the loss fraction. The boundary type specifies whether the macroparticle(s) will be reflected off the boundary, or replaced somewhere randomly in the model. The loss fraction specifies how much the solver should reduce the macroparticle number count by when performing the reflect / replace operation.

[0168] When the boundary condition type is set to reflect, macroparticles will elastically bounce off of the boundaries so that their momentum is conserved. When the type is set to replace, the macroparticles will be removed from the boundary elements and placed somewhere randomly in the model. Even though this replacement 33 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 is random, it is biased so that the PIC environment's temperature and density is conserved over the whole model.

[0169] When performing these boundary operations, the solver multiplies the number count of the macroparticle by one minus the loss fraction. Thus, the loss fraction is the fraction of particles which are lost - for whatever reason - to the boundary.

[0170] The software system includes two types of density boundary conditions: inlets and outlets. The density inlet boundary condition is used to enforce a particular density on a boundary. This is accomplished by first determining what the number count should be in the nearby elements - so-called boundary elements - in order to get the user- specified density on the boundary. This number count is then enforced on the macroparticles which occupy the boundary elements.

[0171] The density outlet boundary condition, on the other hand, is used to overwrite the reflect / replace boundary condition for particular boundaries. A user specifies whether the density outlet should use reflect or replace operations, as well as what the loss fraction is. This density outlet boundary condition is enforced in the same manner as the reflect / replace boundary condition.

[0172] The velocity boundary condition may be used in conjunction with the density inlet boundary condition. There are a few different ways in which a user can define the velocity boundary condition, but in the end, the software system obtains or computes a vector which describes what the velocity components should be on the boundary. These values are then enforced on the macroparticles in the boundary element when the solver enforces the density inlet boundary condition.

[0173] Sources to the PIC solver may now be discussed. These sources are defined as features which act as a source of density, velocity, and / or energy for the PIC solver. Aside from the boundary conditions - which can act as sources for all three quantities - and the electromagnetic solver - which acts as a source of velocity and energy - there are two other sources to be considered: 1. Emission 2. Reactions

[0174] With respect to emission, under certain conditions, conductors can emit electrons from their surface. This process is similar to the photoelectric effect, but proceeds through different mechanisms. The energy which binds an electron to a conductor is described by the work function, ^^, which has units of ^^^^. If an electron 34 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 inside the conductor is provided energy equal to or greater than the work function, then it will be liberated - or emitted - from the conductor. With the photoelectric effect, this energy is provided by a photon. For the emission which is considered, this energy is provided by the electric field strength on the surface of the conductor and / or the temperature of the conductor. If the conductor is not at an elevated temperature, then emission can be described by Fowler-Nordheim-type field emission: ^^^^ଶ^^^^^^^^^^ଷ / ଶ^^^ൌ exp ^െ ^ ^^ ^^

[0175] Where ^^^is the ^^ and ^^ are constants;^^ is the electric field field strength and work function. If the conductor is at an elevated temperature, then Schottky emission (also called field-enhanced thermionic emission) can be used in conjunction with the Fowler- Nordheim-type field emission: ^^ െ Δ^^^^ ൌ ^^ ^^ଶ் ீ exp ൬െ^ ^^^^^

[0176] Where ^^்is the ^^ீis a constant; ^^ isthe temperature of the conductor; and Δ^^ is the energy supplied by the surface electric field.

[0177] In certain embodiments, the software system further comprises a module configured to compute electron injection velocity for particles emitted from a conductive surface under space-charge-limited conditions. The system may determine, at each emission surface element, the local normal electric field ^^^^^^and the nearest interior node distance ^^^. From these quantities, the system calculates the potential change experienced by an emitted electron over the first finite-element cell, represented as ∆^^ ൌ ^^^^^^^^^ .

[0178] The system uses this local potential change to adjust the velocity of emitted electrons according to relativistic energy and momentum invariants. In particular, a Lorentz factor ^^^ൌ ^ is determined for the pre-emission velocity ^^, and a post- ^^ି^^మ ^మemission Lorentz ᇱ^^థ^^ ൌ ^^^ െ^^^మ^ is computed to represent the energy change resulting from the localThe tangential momentum of the electron is conserved such that ^^௧ ൌ ^^^^^^^௧ ൌ ^^ᇱ^^^^ᇱ௧ , where ^^௧ and ^^ᇱ௧ᇲdenote the tangential 35 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 velocity components before and after the update. The new normal velocity component is determined from conservation of energy and momentum as: ଶ ^^ᇱ^ ൌ ^^^^^^^^^ ^^^ ^^^ᇱଶ െ 1 െ^^^^^ ^௧^^^ᇱ^ ^^

[0179] The total ^^^ᇱ^^^ where ^^^ and ^ ^^ are thelocal normal andଶ

[0180] If the computed condition ^^ᇱଶ െ 1 ^ ^ఊబ^^^^ ^ is satisfied, the system determines that the local potential rise is sufficient to emitted electron. In this case, the emission is suppressed. This conditionthat the electrons lacking sufficient kinetic energy to overcome the local electrostatic barrier are not injected into the simulation domain.

[0181] This approach ensures exact energy conservation across the emission boundary and maintains physically correct tangential momentum. It further prevents unphysical potential increases that may arise when the simulation mesh does not fully resolve the near-surface space-charge region. The formulation remains local, depending only on ^^^^^^and ^^^, and can be implemented directly within finite-element or particle-in-cell (PIC) solvers without performing a line integral of the potential.

[0182] In certain embodiments, the emission model described above is implemented within the simulation framework as a distinct emission source referred to as the Space- Charge-Limiting (SCL) emission source. The SCL emission source represents emission that occurs under space-charge-limited conditions and incorporates a non-zero initial kinetic energy corresponding to the injected electron energy at the point of emission. The SCL source utilizes the relativistic velocity adjustment and reflection criteria to determine whether emitted particles are transmitted or reflected and to compute their injection velocities accordingly. In certain embodiments, the software system is configured to utilize the SCL emission source in conjunction with, or independently from the other two emission sources. That way emission is fully characterized at all times.

[0183] With respect to reactions, certain models may include multiple PIC environments and / or fluid environments. Such models - which are typically used for PE-CVD or discharge simulations - could have interactions between the various environments which lead to a change in density, velocity, and / or energy. How exactly 36 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 these reactions can be used as a source of PIC quantities will be discussed will be discussed later in the present disclosure.

[0184] The manner in which the electromagnetic solver can be influenced by the PIC solver through the governing equations has been described. The influence of the PIC solver on the electromagnetic solver can also be modeled, depending on how the user configures the modeling of such influence.

[0185] When PIC source is enabled in a domain settings portion of the software system, the charge density of the PIC environments is used as a source for the electromagnetic fields. A more subtle way of having the PIC solver influence the electromagnetic solver is through the conductivity of the PIC environments, which can be enabled using a PIC conductivity source domain setting. When this is enabled, quantities from the reaction solver are used to compute the conductivity of the PIC environments. These conductivities are then summed and provided to the electromagnetic solver. Since the conductivity calculation depends on reaction solver quantities, the manner in which the conductivity is computed will be described in connection with the reaction solver.

[0186] The manner in which PIC quantities are updated over a given timestep may now be described. This update is carried out through the following steps: ^ Reaction updates are computed if reactions are enabled ^ Reaction updates are used to update macroparticle quantities such as density, velocity, and energy, if reactions are enabled ^ Density and velocity boundary conditions are applied ^ Reflection / replacement is checked for and enforced ^ Macroparticles are propagated (macroparticle positions are updated) ^ Element-based PIC quantities are updated for visualization ^ Field interpolation weights are computed ^ Electromagnetic fields are updated ^ Macroparticles are accelerated (macroparticle velocities are updated) FLUID SOLVER

[0187] While the PIC solver is very useful for simulating plasmas, more efficient ways to simulate neutral species may be employed. This fact, along with the recognition that plasmas can also be treated as a fluid, motivated the inclusion of a fluid solver inside 37 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 of the software system. Like with conventional CFD software, the fluid solver included in the software system uses the FEM. As such, a matrix equation of the following form is sought: ^^^^^^^ ^ ൫^^^ ⋅ ^^^^൯^^^ ൌ ^^^

[0188] In this case, two hand side as the fluid equations solved by the softwareterms. The stiffness matrix, ^^^^, is used to represent these non-linear terms.

[0189] In this section of the present disclosure, the construction of continuity equations for the fluid mass density, temperature, and velocity is described; how these equations are put into a form that can be solved with the Finite Element Method is described; and how the matrices in the equation above are constructed and solved is described.

[0190] All of the fluid continuity equations come from a single, generic continuity equation: ^^^^ ^^^^^ ^^ ⋅ ^^^ ^ ^^^^^ െ ^^ ൌ 0

[0191] Where ^^ is the transported flux; ^^ is the fluxin the absence of fluid transport; and ^^ is a term that represents the source / sink of fluid quantity ^^.

[0192] This generic continuity equation can be used to derive continuity equations for the fluid mass density, temperature, and velocity. The software system is configured to use both a compressible and incompressible fluid description, but for the sake of the present disclosure, the focus is on the incompressible fluid equations.

[0193] A continuity equation for the mass density of an incompressible fluid is first derived. In this case, ^^ ൌ ^^, where ^^ is the mass density. Furthermore, mass can onlybe transported by a fluid flux, and thus ^^ ൌ ^^ for the mass density continuity equation.The transported flux term becomes ^^ ⋅ ^^^^^^ ൌ ^^ ⋅ ^^^^ ^ ^^^^^ ⋅ ^^^. An incompressiblefluid is defined as one in which the fluid velocity is divergence free - i.e. ^^ ⋅ ^^ ൌ 0 -and the gradient of the density is approximately zero. This leaves only the time derivative and the source / sink term, which is written as ^^ and moved over to the right- hand side: ^^^^ ^^38 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0194] Next, the continuity equation for the fluid velocity is considered. While this equation could be derived by hand, the Navier-Stokes equation for incompressible fluid flow is used: ^^^^ଶ^^^^^ ^^^ ⋅ ^^^^^ െ ^^^^ ^^ ൌ ^^

[0195] Where ^^ is the and ^^ is the kinematic viscosity. It is noted thathand side of this equation is nonlinear in ^^.

[0196] To get the continuity equation for the temperature, ^^ ൌ ^^^^^^^ is plugged intothe generic continuity equation. Where ^^^is the heat capacity of the fluid, and the quantity ^^^^^^^ is the amount of heat per unit volume. Unlike with density, heat can flow in the absence of fluid transport. To account for this, ^^ ൌ െ^^^^^^ is used, where ^^ is thethermal conductivity of the fluid. After plugging everything in and rearranging the equation to get it purely in terms of the temperature, ^^, the following is obtained: ^^^^ଶ^^^^^ ^^^ ⋅ ^^^^^ െ ^^^^ ^^ ൌ ℎ

[0197] Where ^^:ൌ ℎ is the source / sink term forthe fluid temperature.

[0198] The continuity equations are transformed into a form that can be implemented into the Finite Element Method. This transformation is accomplished by utilizing Galerkin's method and backwards differencing, similar to the approach employed for the electromagnetic governing equations described previously.

[0199] The governing equation for the fluid mass density of an incompressible fluid is given by: ^^^ ^^^^^^^^^^ ^^^^ ൌ ^ ൫^^^^ Δ^^ ^ ^^^^ି^൯^^^^^^^^^^^^

[0200] For the velocity of an incompressible fluid: ^^^^^^^^^^^^^^ ^ ^^Δ^^൫^^^^^^ ⋅ ^^^^^^ ൯൧^^^^ െ ∮ ^^^^^^Δ^^^^^^^^^^^^ ⋅ ^^^^ ^ ^ Δ^^^^^^^^^^^ ⋅ ൫^^^^^^^^^^^^^^൯൧^^^^^^^^ ^^^^^ ^^^^ ^ ^ Δ^^^^൫^^^^^ ⋅ ^^^^^^ ^ ൯൧^^^^ െ ∮ ^^^Δ^^^^^^^^^^^^ ⋅ ^^^^ ^ ^ ^^^Δ^^^^^^ ^ ^ ^ ^ ⋅ ൫^^^^^^^^^^^^^^൯^^^^to all domain boundaries. Furthermore, it can be difficult for the solver to determine which boundaries should be given a natural boundary condition a priori. As such, the 39 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 fluid solver may rely on the user specifying boundary conditions for all boundaries in order to enforce the natural boundary conditions. How exactly these boundary conditions are defined and enforced is described below.

[0203] The final forms of the governing equations for the fluid mass density, velocity, and temperature have been established. How these governing equations can be used to construct the below matrix equation will be described: ^^^^^^^ ^ ൫^^^ ⋅ ^^^^൯^^^ ൌ ^^^

[0204] Where the stiffness for nonlinear terms. Ingeneral, the matrix follows roughly the same steps the electromagnetic solver took. Most notably, element-based versions of the matrices above could be defined, constructed first, and then summed over elements to get the global matrices. While this is the formulation that is used in the present section of the disclosure, the fluid solver does directly compute the global matrices; just like the electromagnetic solver does. The fluid solver also uses the same finite element mesh as the electromagnetic solver. Therefore the fluid solver can use the same basis function integrals during matrix assembly.

[0205] The governing equation for the incompressible fluid mass density does not have a nonlinear term. Therefore only the mass and force matrices will have nonzero entries. The element-based contributions to the matrices can be directly constructed from the governing equation. These contributions can then be summed over all elements to get the global matrices: ^^^ ^^^ ൌ ^ ^^^^^^^ ^^^^ → ^^^^ ൌ ^   ^^^ ^^^ ^

[0206] Theincludes a nonlinear term. The element-based contributions to the matrices can be directly constructed from the governing equation, and these contributions can then be summed over all elements to obtain the global matrices: 40 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 ^^^^^ ൌ ^ ^^^^^^^^^ ^ ^^Δ^^൫^^^^^ ^^ ⋅ ^^^^^ ൯൧^^^^ െ ∮ ^^Δ^^^^^^^^^^^^ ⋅ ^^^^ → ^^^^ ൌ ^   ^^^ ^^^ ^^^^ thethe surface integral in the mass matrix is computed only on the boundaries present in the model.

[0208] The governing equation for the incompressible fluid temperature also has a nonlinear term. The construction of global matrices is performed using the same process as described previously: ^^^ ^ ^ ^^^ ൌ ^ ^^^^^^^^ ^ ^^Δ^^൫^^^^^ ⋅ ^^^^^ ൯൧^^^^ െ ∮ ^^Δ^^^^^^^^^^^^ ⋅ ^^^^ → ^^^^ ൌ ^  ^^^^^ ^ matrixis computed only on the boundaries included in the model.

[0210] The calculation of the mass, stiffness, and force matrices in the fluid solver proceeds through the same steps that are used for this calculation in the electromagnetic solver. Two primary differences exist: 1. A stiffness matrix is computed in addition to the mass and force matrices 2. The stiffness matrix uses a new basis function integral which must be computed

[0211] The initial conditions for the fluid quantities are as follows: ^^^^ ^^^∀^^ ∈ Ω, ^^^^ ൌ ^^^^^^^^∀^^ ∈ Ω, ^^^^ ൌ 0^^^∀^^ ∈ Ω, ^^ ^ ൌ ^^^^^^ ^ ^^^^^∀^^ ∈ Ω, ^^^^ ൌ ^^^^^∀^^ ∈ Ω, ^^^^ ൌ ^^^ ∈ Ω, ^^^^ ൌ 0

[0212] Where ^^^and ^^^are the user-defineddensity and temperature, respectively. The velocity is set to zero initially so that the solution for the first timestep is static, thereby ensuring stability. 41 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0213] The initial conditions are enforced by initializing the arrays which store the values of the fluid quantities from the previous timestep so that the initial conditions are satisfied everywhere.

[0214] Dirichlet-like boundary conditions are supported by the fluid solver for use on the fluid quantities: ^^^∀^^ ∈ ^^Ω^ ൌ ^^^^^∀^^ ∈ ^^Ω^ ൌ ^^^^

[0215] Where ^^,^^, and ^^ are conditions on surface ^^Ω for the fluid mass density, velocity,A density outlet boundary condition is also supported by the fluid solver, which allows specification of how many fluid particles will leave the system through the boundary. This is specified by using a loss fraction. The density outlet boundary condition can be considered as a time-varying Dirichlet-like density boundary condition. The boundary condition for a given timestep is given by multiplying the current density on the boundary by the loss fraction: ^^^^∀^^ ∈ ^^Ω ^ି^outlet ^ ൌ ^^ ∗ ^^ ^^^^

[0216] Where ^^ is the loss fraction, boundary is an outlet, and^^ ∗ ^^^ି^^^^^ is the time-varying value of the boundary condition which can beconsidered as ^^^^^^. Enforcement of these Dirichlet-like boundary conditions on the fluid solver is performed in the same manner as in the electromagnetic solver. After construction of all matrices needed to solve for a given fluid quantity, nodes on boundaries with specified boundary conditions are processed. For each node, designated as node ^^, the force matrix is reset so that the entry for the node is equal to the value of the boundary condition: ^^^^^^^ ൌ ^^. Where ^^ is the condition for theboundary that node ^^ lies on. Additionally, the row for node ^^ in the mass matrix is reset so that the only nonzero value is: ^^^^^^^, ^^^ ൌ 1.0. The row for node ^^ in the stiffnessmatrix is reset so that all entries in the row are zero. By this process, the matrix equation for node ^^ becomes: ^^^^^^^, ^^^^^^ ൌ 1.0 ∗ ^^^ ൌ ^^

[0217] The solution for the^^^, is thereby determined to be equal to the value of the boundary condition, ^^.

[0218] Depending on the setup of the model, the fluid quantities can receive contributions from two sources: 1. Electromagnetic fields 42 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 2. Reactions

[0219] If a fluid species is charged, then the fluid species will be accelerated by the electromagnetic fields. In such cases, the electromagnetic fields act as a source of fluid temperature and velocity. The electric force is treated as a source with contributions calculated from Coulomb's electric force law and kinematic equations. The contribution to the velocity from the electric force is: ^^ ^^ாெ,^ൌ ^^ ^^^

[0220] Where ^^ாெ,^is the force, ^^, from the electric force onnode ^^ and ^^^is the electric field contribution to the temperature from the electric force is: ^^ଶℎாெ,^ൌ ^^ଶ^^^^^^

[0221] Where ℎாெ,^is the contribution to the heat source, ℎ, from the electric force on node ^^. If reactions are defined in the model, then the reactions can act as sources to the fluid density, velocity, and temperature. In general, the contributions to the fluid quantities are accumulated when the software system constructs the force matrices, ^^^.

[0222] The steps taken by the fluid solver during matrix assembly are substantially the same as those taken by the electromagnetic solver. The fluid quantities - ^^, ^^, and ^^ - can be solved independent of each other. The order in which the fluid quantities are solved is as follows: 1. Temperature 2. Density 3. Velocity

[0223] Since the temperature is solved for first, the velocity that is used in the nonlinear term in the temperature governing equation comes from the previous timestep.

[0224] The mass and stiffness matrices for the fluid quantities are sparse, similar to the mass matrices in the electromagnetic solver. This allows BiCGStab (an iterative method / algorithm for finding the numerical solution of nonsymmetric linear systems) to be used to solve the matrix equations which are constructed for the fluid quantities.

[0225] Similar to the PIC solver, the fluid solver can influence the electromagnetic solver. There are two ways in which this influence can be handled: 1. Charge density 43 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 2. Conductivity

[0226] If fluid source is enabled in the domain settings, then the charge density of the fluid environments - should the fluid environments be charged - will be used as a source for the electromagnetic fields. If fluid conductivity source is enabled in the domain settings, then the reaction solver will compute the conductivity of the fluid environments. After summing up the contributions from all fluid environments, the conductivity is then provided to the electromagnetic solver. REACTION SOLVER

[0227] When simulations are run with multiple PIC and / or fluid environments, it is possible for the environments to interact with each other, depending on what the environments are and what is being simulated. These interactions can lead to a change in the density of the environments, their energies, or both. Depending on what is to be observed in the simulation, accounting for these interactions and their effects can be very important. This is accomplished through the reaction solver.

[0228] Reactions may be included for a variety of reasons. In the case of plasma enhanced chemical vapor deposition (PE-CVD) - a semiconductor manufacturing technique - the reactions govern the energy of the plasma to be used for the etching process. The reactions also determine the rate of etching as the plasma interacts with the substrate. When the breakdown of a gas / liquid is simulated, the reactions govern the formation of the discharge. The reactions then determine what the arc will do after breakdown; whether it will sustain or extinguish. The reactions even determine what the conductivity of the plasma is.

[0229] It is common to categorize reactions as elastic or inelastic. In elastic reactions, the kinetic energy of the interacting particles is conserved. Inelastic reactions do not conserve the kinetic energy; instead some kinetic energy is "lost" to a variety of mechanisms. In some cases the lost energy goes towards breaking chemical, molecular, or nuclear bonds. In other cases the energy can be lost in the form of heat, or by causing one of the interacting molecules to vibrate / rotate around its molecular bond.

[0230] Reactions can also be discussed in terms of how many participating particles, or species, there are. The most common type is a two-body reaction, in which two species interact with each other. However, it is also possible to have a three-body 44 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 reaction, which has three participating species. There are a few ways to write out reactions, but the format used by the software system is: ^^ ^ ^^ → ^^ ^ ^^ ^ ^^loss

[0231] Where ^^ and ^^ are the participating species, ^^ and ^^ are the product species, and ^^lossis the energy that is lost by inelastic reactions in units of electron-Volts, ^^^^. When this format is used to define an elastic reaction ^^lossis zero, since no energy is lost. The product species can be the same as the participating species - in which case the reaction equation would read: ^^ ^ ^^ → ^^ ^ ^^ ^ ^^loss - but they do not have to be.It is also possible to have only one of the product species be the same as one of the participating species. Furthermore, it is possible to define a reaction with more product species than participating species.

[0232] It can also be helpful to have names to describe the participating species. One of the participating species is referred to as the incident species, whereas the other species is dubbed the target species. While either participating species could technically be the incident species, for the purposes of the present disclosure, it is assumed that the incident species is either a charged particle, or the participating species with the smaller mass - in the case that both of the participating species are charged particles, or neither of them are.

[0233] As will be seen, most reactions are considered to be inelastic. This can make discussion of inelastic reactions somewhat complex, as simply calling them "inelastic" does not provide any clarity on what reaction is specifically being discussed. As such, the inelastic reactions may be broken apart into their "types." There are two types of elastic reactions, which are called momentum transfer and Coulomb scattering. Both of which are written as: ^^ ^ ^^ → ^^ ^ ^^ ^ 0eV

[0234] The primary difference between these two elastic reaction types is what the participating species are. If one or both participating species is neutral, then that would be called a momentum transfer reaction. Meanwhile, if both participating species are charged particles, then a Coulomb scattering reaction is present.

[0235] For most ionization reactions, the incident species, ^^, would be the electrons and the target species, ^^, could be either a neutral species or an ion. Although it is much more common for the target species to be a neutral than an ion. Furthermore, it is 45 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 possible for ions to be the incident species. Assuming the incident species is an electron, ionization reactions may be written as: ^^ ^ ^^ → 2^^ ^ ^^ା ^ ^^loss

[0236] Where an additional ^^ appears on the right hand side - denoted by the preceding 2 - because an electron is produced by the ionization event, and the original electron is not consumed in the process. ^^losswill be nonzero for ionization reactions, and this lost energy is spent freeing the produced electron from the target species. In the context of ionization reactions, ^^lossis commonly called the ionization energy.

[0237] There also exists a reaction which is the inverse of ionization, commonly referred to as recombination reactions. For this reaction, an electron (incident species) will be "captured" by an ion (target species) to produce a neutral (product species): ^^ ^ ^^ା → ^^ ^ ^^loss

[0238] Both the electrons and ions are "consumed" in this reaction, which is why they do not appear on the RHS of the equation. ^^losscould be nonzero, but it is actually typically zero. The reason why this reaction is still considered to be inelastic - even when ^^loss ൌ 0 - is due to a subtlety with how energy conservation is defined. Strictlyspeaking, energy conservation is defined in terms of the energy of the participating species. Since both of the participating species are consumed in the recombination reaction, their energy is inherently not conserved. Of course, their energy is given to the product species ^^, so energy is conserved over the whole system. Which is what leads to ^^lossbeing zero. In the event that ^^lossis nonzero, it is likely because the target molecule actually prefers being an ion due to orbital filling and the Pauli exclusion principle. However, this is unlikely and typically only occurs for polyatomic ions with a lot of atoms.

[0239] Closely related to the ionization reaction type is the dissociative ionization reaction type. This reaction can occur when the incident species is a charged particle and the target species is a diatomic / polyatomic neutral molecule. This means the molecular structure of the target species includes more than one atom, such as for Nitrogen gas ^^^ଶ^. Note that this reaction type can only occur if the target species is diatomic / polyatomic. This reaction type ends up having more product species than participating species: ^^ ^ ^^ → 2^^ ^ ^^ ^ ^^ା ^ ^^loss46 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0240] In this reaction, the target species will dissociate - or break apart - into the product species ^^ and ^^. One of which ends up ionizing. ^^losswill be nonzero, and this energy is spent doing two things: breaking apart the target molecule, and ionizing one of the resulting products. In general, dissociative ionization reactions will lose more energy than plain ionization reactions. As additional energy is needed to break the molecular bond of the target species.

[0241] It is also possible to have dissociation without any ionization. Just like with the dissociative ionization reaction type, the target species for dissociation reactions has to be a diatomic / polyatomic molecule. The reaction equation for dissociation reactions can be written as: ^^ ^ ^^ → ^^ ^ ^^ ^ ^^ ^ ^^loss

[0242] In this case, no electron is produced as no ionization occurs. ^^losswill still be nonzero, and is spent breaking the molecular bond in the target molecule. In some cases, it is also possible to produce negatively-charged ions instead of positively-charged ions. Such reactions are referred to as attachment reactions as an electron will "attach" to some neutral molecule: ^^ି ^ ^^ → ^^ି ^ ^^loss

[0243] Like with recombination reactions, attachment reactions consume the electron that participates in them. For most target species, ^^losswill be nonzero. However, the attachment reaction is unlikely to occur when ^^lossis nonzero. As a result, it is typically only important to include attachment reactions when ^^loss ൌ 0. As such, attachmentreactions usually have ^^loss ൌ 0 although this is not strictly true. Even when ^^loss ൌ 0,this reaction is still considered to be inelastic owing to the fact that the participating species do not appear on the right-hand side of the equation.

[0244] Electrons can also be "detached" from negatively-charged ions. Such reactions are called detachment reactions: ^^ ^ ^^ି → ^^ ^ ^^ ^ ^^ି ^ ^^^^^^

[0245] In most cases the incident species will also be an electron, which would allow the reaction to be rewritten as: ^^ି ^ ^^ି → 2^^ି ^ ^^ ^ ^^loss . Detachment reactions canhave a zero ^^loss, but it is more commonly nonzero as energy has to be given to the produced electron in order to liberate it from the participating negative ion.

[0246] Negatively-charged ions can participate in recombination reactions, but such reactions are typically referred to as mutual neutralization reactions: 47 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 ^^ି ^ ^^ା → ^^ ^ ^^ ^ ^^loss

[0247] It is fairly common for the negative and positive ions to come from the same neutral species. In which case the equation can be rewritten as: ^^ି ^ ^^ା → 2^^ ^ ^^loss .These reactions, just like with recombination, will also typically have ^^loss ൌ 0. Thisis because the electric field produced by both ions usually provides enough energy to free the extra electron from the negative ion and to "give" it to the positive ion.

[0248] The final reaction type to be considered are excitation reactions. In these reactions a charged particle - typically an electron - will impact the target particle - typically neutral - putting it into an excited state: ^^ ^ ^^ → ^^ ^ ^^∗ ^ ^^^^^^

[0249] Where the superscript * denotes that the species is in an excited state, and ^^lossis spent putting the target species into said excited state. For these excitation reactions, it is common to call ^^lossthe excitation energy.

[0250] The excitation reactions can be further categorized by defining sub-types that relate to how exactly the target species is being excited, although the format used to write the reaction equation would be the same for all sub-types. In doing so three sub- types of excitation reactions are defined: rotational, vibrational, and electronic. Rotational and vibrational excitation can only occur when the target species is a diatomic / polyatomic molecule. Rotational excitation reactions will result in the target molecule rotating about its molecular bond(s). The axis of rotation as well as the frequency will determine the excitation energy of these rotational excitation reactions. Vibrational excitation reactions are similar, but result in the target molecule vibrating about its molecular bond(s) instead of rotating. The vibrational frequency will determine the excitation energy for these reactions. The last excitation sub-type, electronic excitation, can occur for any target species.

[0251] In electronic excitation reactions, the energy of the incident charged particle is given to an electron in one of the orbitals of the target species. This results in a "promotion" of the electron to a higher orbital. The excitation energy for these reactions is based off of the energy difference between the initial and final orbitals.

[0252] The manner in which the reactions are modeled will now be described. The end goal is to be able to describe the change in density, energy, and / or momentum of all participating and product species for a given reaction. Before that point is reached, how the reaction rate of a particular reaction can be obtained will be described. The reaction 48 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 rate is defined as the number of reaction events that are occurring per unit volume, per unit time.

[0253] Two formalisms for computing the reaction rate are provided by the software system: the analytic formalism, and the cross-section formalism. Data for either formalism can generally be found in literature.

[0254] In the analytic formalism, the reaction rate is computed using a parameterized equation. Three equations that can be used in this formalism are supported by the software system. These equations, along with which values are supplied by the user, are given in the table below. Formula EquationParameters

[0255] ere^s e reac on ra e or reac on n un s o ;^s the number density of the incident species in units of ^^ିଷ;^^௧is the number density of the target species in units of ^^ିଷ; and ^^^is the temperature of the incident species in units of Kelvin. For the parameters, the units of ^^^are ^^ଷ^^ି^;^^^is unitless; and ^^^has units of Kelvin.

[0256] When reactions are defined using the analytic formalism, the user need not specify which equation is being used to the software system. Since each equation has a unique number of parameters, the equation to use can be determined by the number of parameters that are provided when defining the reaction.

[0257] While the analytic formalism is convenient for defining reactions in an easy manner, the parameterization may lead to a loss of fidelity. For reactions whose rate does not depend heavily on the energy of the incident species, this loss is minimal. However, not all reactions are so insensitive to the incident species' energy.

[0258] When the reaction rate depends strongly on the energy of the incident species, the cross-section formalism can be used to obtain the reaction rate with a high 49 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 fidelity. Cross-section data is provided as a function of incident species energy with units of area, ^^ିଶ. Conceptually, the cross-section represents the probability of a reaction occurring; a larger cross-section represents a larger target with a higher probability of hitting said target.

[0259] To obtain the reaction rate from the cross-section data, the following equation is employed by the software system: ^^^ ൌ ^^^^^௧^^^^^^^^^^^^^^^

[0260] Where ^^^^^^^^is the^^^^^is the magnitude of the incidentspecies' velocity; and value of the cross-section datamultiplied with the incident species. The expected value is used here asnot every particle in possess the same velocity. Computing the expected value is performed as follows: ^ ^^^^^^^^^^^^^^^ ൌ ^  ^^^^^^^^^^^^^^^^^^^ ^

[0261] Where ^^^^^^^ is provides the probability offinding a particle with energy ^^^out of a group of particles with temperature - i.e. average energy - ^^^. The use of the Maxwellian distribution is supported by the software system.

[0262] In order to use either formalism, information on the participating species is needed: namely, the densities and temperatures of these species. When one of the participating species is modeled as a PIC environment a question arises: how should PIC quantities be used to compute the reaction rate?

[0263] This question comes from the fact that the PIC quantities are stored in two forms. The PIC quantities are stored on a macroparticle basis, which is then used to construct element-based PIC quantities. So the question can be phrased as: should the macroparticle-based or element-based quantities be used to compute reaction rates?

[0264] The answer to this question depends on the conditions present in the simulation. For PECVD simulations where the plasma is fairly sparse - as in the density is less than 1E ^ 016^^ିଷ - then using the macroparticle-based quantities provides the mostaccurate results. However, for discharge simulations where the plasma is moderately dense the macroparticle-based quantities no longer provide an accurate description of the plasma. Instead, using the element-based quantities provides the best results. 50 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0265] The software system allows the user to specify whether the macroparticle-based or element-based PIC quantities are to be used. The scheme which uses macroparticle- based quantities is referred to as particle-based, whereas the scheme that uses element- based quantities is referred to as continuum-based. To use the continuum-based reaction rate calculations, the user can enable a "PIC Continuum Scattering" setting.

[0266] Once the reaction rate has been calculated for a given reaction, it can be used to update the density, energy, and / or momentum of the participating and product species.

[0267] The density update is first considered. The first step in computing the density update is to determine the change in particle count from the reaction equation. This is most easily understood by considering examples. For the density updates for the ionization of ^^ଶ, the equation used to define this reaction would be: ^^ି ^ ^^ଶ → 2^^ି ^ ^^ାଶ ^ ^^loss

[0268] To determine the change in particle counts, the LHS can be subtracted from the RHS. In doing so, it is seen that the particle count increases by 1 for the electrons and ^^ଶା- since those particles are produced - and decreases by 1 for ^^ଶ- since it is consumed. Δ^,^is used to express the change in the particle count of species ^^ for reaction ^^. Once all of the changes in particle count have been computed, obtaining the rate of change of the density is straightforward: ^^^^ ൬^^ ^^^^ൌ Δ^,^^^^^

[0269] It is noted that allwill have a change in density; assuming Δ^,^is nonzero for all species.

[0270] To obtain the rate of change of a species energy, the change in energy per reaction, Δ^^col, is first determined. For inelastic reactions the change in energy is given by ^^loss, and for elastic reactions it can be computed from the conservation of energy and momentum: 2^^ Δ^ ൌ^^^ ^௧^^^ ^^^ ^^ ^ ^^ ^ଶ ^^^ െ ^^௧^௧

[0271] Where ^^^is the mass of the incident^^௧is the mass of the target species; ^^^is the energy of the incident species; and ^^௧is the energy of the target species. The rate of change of a participating particle's energy can then be expressed as: ^^^^^Δ^^ ^ ^^^^^^51 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0272] The sign of the right-hand side is determined by which species the energy rate is applied to. Since Δ^^colis defined to be the energy lost / gained by the target species, the sign of the right-hand side is positive when applying the energy rate to the target species, and negative when applying it to the incident species. The energy rate is only applied to the target species when reaction ^^ is elastic. Otherwise the energy lost by the incident species goes towards doing work and is irreversibly lost.

[0273] The quantity േΔ^^col^^^represents the change in total energy density of the species to which the energy rate is applied. However, the average change in a single particle's energy is obtained by dividing the quantity by the number density of the species.

[0274] The final reaction update to be considered is that for the species' momentum. This quantity is computed for elastic reactions, in which case the change in a particle's velocity can be computed using conservation of momentum and energy: 2^^ Δ^^ ൌ௧ ^^ െ ^^^^^ ^^ ^ ^^ ^^ ^ ൬ ^ ௧௧^ ^^ ^ ^^ െ 1^ ^^^^

[0275] Where ^^ is ^^௧is the velocity oftarget species. The rate of change of a species' momentum can then be expressed as: ^^^^ ^^ Δ^^ ^^ ൬^^ ^ ^ൌ േ ^ ^^^ ^^^ ^ ^^^

[0276] Like with the energy hand side is determined by thespecies to which the momentum rate is applied. The difference comes from the fact that Δ^^colis defined as the velocity lost / gained by the incident species. So the sign convention is opposite to what it was for the energy rate; the sign is positive for the incident species, and negative for the target species.

[0277] The quantity that appears in the numerator on the right-hand side is the change in the total momentum density of species ^^. However, the average change in a particle's momentum is obtained by dividing the right-hand side by the number density of species ^^.

[0278] The momentum rate is computed for elastic reactions. The momentum rate calculation can be further limited to only when the elastic reaction has an ion as the incident species and neutral as the target species. The momentum rate need not be calculated when one of the participating species is an electron. This is due to the randomness of the electrons' motion as compared with the ions' motion. Because the electron mass is small, electrons are heavily scattered during elastic reactions. This 52 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 leads to their motion being much more randomized than it would be for the ions. Due to this randomness in the electron motion, if Δ^^^^^is averaged over multiple collisions, it is found to be approximately zero.

[0279] Once the rate of change of the density, energy, and / or momentum has been determined, those rates can be applied to the PIC and / or fluid environments. A distinction is made between PIC and fluid environments since each environment type tracks different quantities.

[0280] PIC quantities are tracked on a macroparticle-basis. These macroparticle-based quantities are then summed / averaged in order to get element-based quantities. When applying reaction updates to a PIC environment, the updates are directly applied to the macroparticles.

[0281] The macroparticle quantities that are tracked are the number count - the number of particles contained within a macroparticle - and the velocity. Therefore the reaction updates are framed in terms of these quantities. The manner in which the density update is applied to the macroparticle number count is described first.

[0282] The change in the density is obtained by multiplying the rate of change of the density with the timestep size. This can then be divided by the density of the element to get how much the density changes as a fraction of the original density: ∑ ^^^^ ^   ൬ ^^^^^ ^Δ^^ ^

[0283] Where ^^dens ,^is theof PIC species ^^, and the sum over ^^ on the RHS is performed to get the total rate of change in the density of PIC species ^^.

[0284] Subsequent processing depends on whether ^^dens ,^is positive or not. If it is negative, then all macroparticles inside of the element for which ^^dens ,^was computed are processed. The number count is then set to be ൫1 ^ ^^dens ,^൯ times the original value.When ^^dens ,^is positive, then the number of particles that are created is determined, which is given by the numerator of ^^dens ,^. New macroparticles are then created to hold these new particles.

[0285] The manner in which the energy and momentum rates are applied is described. The fractional change in the particle's average energy is computed, similar to what was performed with the density: 53 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 ∑ ^^^^ ^   ൬ ^^^^^ ^Δ^^ ^

[0286] The macroparticles element for which ^^^^,^was computed are then processed.each macroparticle is computed using: ì ï^1 ^ ^abs൫^^^^,^൯^^^^^ௗ,^, if ^^^^,^ ^ 00

[0287] This for the scatteringangle. After contributions from the momentum rates are added. Before doing so, the momentum rates are converted to a change in velocity, which can then be added to the updated macroparticle velocity. The expression for the updated velocity which accounts for both the energy and momentum rates is: Δ^^ ^^^^ ^^^^௪,^ ൌ ^1 േ ^abs൫^^ ^^^,^൯^^^^^ௗ,^ ^ ^  ൬ ^ ^^^^ ^^^^ ^

[0288] ^^^^,^. The randomreorientation of the velocity is only performed on the first term in the equation above; the change in velocity that comes from the momentum rates is assumed to be highly directed. Hence the reorientation is not applied to that term.

[0289] When using the particle-based rate calculation, these reaction rates are used to update the macroparticle quantities as they are computed. In other words, the reaction rates are applied on a reaction-by-reaction basis when using the particle-based calculation. However, when using the continuum-based rate calculation, the reaction rates are only applied after computing the contributions from all reactions.

[0290] Fluid quantities are tracked on a nodal-basis. The first step in applying the reaction updates to the fluid involves converting element-based rates to being nodal- based. This conversion is performed by averaging the reaction updates over the elements which contain a particular node. Once this averaging is performed, the reaction updates can be applied to the fluid species.

[0291] The fluid solver tracks the mass density, temperature, and velocity of the fluid species. To obtain the change in mass density, the rate of change in the density is first 54 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 multiplied by the mass of the fluid species. This can then be incorporated into ^^ on the right-hand side of the fluid mass continuity equation: ^^^^ ^^ ൌ ^^^ ^  ൬^^ ^ ^^^^ ^

[0292] To obtain the change is treated as an ideal gas. This treatment allows the averageto the temperature of a group of particles. The ideal gas law can be used to convert the energy rate to a temperature rate. This can then be incorporated into ℎ on the right-hand side of the fluid temperature continuity equation: 2^^ ^^^^ ℎൌ ^^  ൬^^ 3 ^^^^ ^ ^

[0293] Where the prefactor, relates the average particleenergy of an ideal gas to the

[0294] Obtaining the change in fluid velocity is straightforward. The momentum rate is divided by the mass of the fluid species. This quantity can then be incorporated into ^^ on the right-hand side of the fluid velocity continuity equation: ^^ ൌ1 ^^^^^^^^  ൬ ^ ^ ^^^^ ^

[0295] After computing the ^^ for all nodes, the fluid solveraccounts for the reaction updates when solving for the timestep's solution.

[0296] How the reaction updates are applied to the fluid environments does not depend on whether the rate calculation was particle-based or continuum-based. In both cases, the contributions to the reaction updates are summed for all reactions and used to construct ^^, ℎ, and ^^.

[0297] When PIC conductivity and / or fluid conductivity is enabled (e.g., by a user), the reaction solver performs additional calculations which are used to update the plasma conductivity. Under these conditions, the solver computes the collision frequency of each reaction. This can be computed by dividing the reaction rate, ^^^, by the density of the incident species: ^^ ^^^^,^ൌ ^^^

[0298] Where ^^^,^is the collision frequency of species ^^ from reaction ^^. These collisioncan then be summed over all reactions and incorporated into the 55 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 Drude model in order to obtain the conductivity of the species. By summing over all species, the total plasma conductivity is obtained: ଶ ^^ ൌ  ^^^^^ ^^^^^^^൫∑^    ^^^,^൯

[0299] This calculation is basis so that the plasma conductivity can be directlysolver. Since the fluid quantities are nodal-based, some data manipulation is required in order to obtain the element-based fluid conductivity. However, this is accomplished by averaging the nodal-based quantities over the nodes in an element to obtain the corresponding element-based quantity.

[0300] The system is configured to compute plasma conductivity based on how the plasma interacts with itself as well as the surrounding medium. This approach allows the plasma conductivity method to be generalized, enabling the system to simulate novel discharge problems which traditional plasma conductivity methods may struggle with due to lacking experimental data. The conductivity calculation accounts for the dynamic interactions between plasma species and their environment. EXAMPLES

[0301] Referring to FIG. 1, an example computational workflow for the software system is shown in the form of a flow chart, and is generally indicated at reference number 1000 (e.g., the method 1000).

[0302] The system is initialized at step 1002: simulation settings are established, finite element method matrix assembly is performed, MPI / Metis domain divisions are implemented, shape functions are defined, particle-in-cell parameters are initialized, scattering rates are set, and boundary and field initial conditions are established.

[0303] The physical environment may be defined by a user through the integrated 3D CAD module during the initialization process at step 1002. In some aspects, the user may create or import three-dimensional geometries that represent the simulation domain (e.g., vacuum chambers, plasma reactors, discharge vessels, spacecraft, etc.). The 3D CAD module may enable the user to specify geometric boundaries, material regions, and surface characteristics that define the physical space where plasma discharge phenomena will be simulated. During step 1002, the user may enumerate various physical parameters and initial conditions corresponding to the defined physical environment (e.g., electromagnetic properties of materials, thermal conductivities, 56 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 surface emission characteristics, plasma densities, temperatures, boundary voltages, positions, etc.). In some cases, the user may also define mesh parameters that control the finite element discretization of the CAD geometry (e.g., to ensure that the computational mesh accurately represents the physical environment while maintaining numerical efficiency for the subsequent plasma simulation calculations).

[0304] Following initialization, at step 1004, a central coordination point for finite element method operations is provided including boundary condition updates, application of particle-in-cell and fluid charge distributions, potential updates, and implicit field calculations. Step 1004 is a point of convergence for three distinct computational cycles that operate at different temporal scales and address different physical phenomena within the plasma simulation.

[0305] The cycle time-stepping loop is operated with an accelerated time-step and encompasses steps 1006, 1008, 1010, 1012, and 1014 before returning to step 1004. At step 1006, finite element method average values are calculated over a single cycle consisting of N standard time-steps. At step 1008, these average values are compared to previous cycle results to determine whether the accelerated time-step remains acceptable for maintaining numerical accuracy and stability. At step 1010, the cycle averages are utilized to calculate mass and heat rates through the reaction solver, incorporating chemical kinetics and energy transfer processes. At step 1012, ion species are updated based on the cycle averages and electrons are subsequently updated based on the ion changes, maintaining charge neutrality and species coupling throughout the plasma discharge region. At step 1014, an implicit solution of Navier-Stokes equations is found and fluid charge is calculated.

[0306] The fluid time-stepping cycle is operated with a fluid time-step and includes steps 1024 and 1026 before returning to step 1004. At step 1024, implicit solution of Navier-Stokes equations is performed and fluid charge distributions are calculated for species treated using continuum fluid modeling approaches. At step 1026, particle-in- cell macroparticle renormalization is implemented, particle positions are updated, and charge interpolation to the finite element method mesh is performed using shape functions.

[0307] The standard time-stepping cycle is operated with an enforced time-step determined by plasma frequency constraints and encompasses steps 1016, 1018, 1020, and 1022, followed by either step 1024 for fluid simulation or step 1026 for particle-in- cell simulation, before returning to step 1004. At step 1016, particle velocities and 57 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 positions are updated using the enforced time-step that resolves plasma frequency phenomena. At step 1018, charge conservation field corrections and energy conservation particle corrections are performed to maintain conservation laws throughout the simulation. At step 1020, emission sources, particle-in-cell loss mechanisms, and surface reactions are updated.

[0308] At step 1022, reaction solver operations are represented that calculate reaction rates, update particle-in-cell macroparticles, and apply force terms to fluid equations for coupled fluid-kinetic modeling.

[0309] The three computational cycles are operated with independent time-stepping schemes that accommodate the different temporal scales of electromagnetic, kinetic, and fluid phenomena. The cycle time-stepping is utilized with accelerated time-steps that span multiple plasma periods, enabling efficient simulation of long-term plasma evolution. The fluid time-stepping is operated at time scales appropriate for fluid transport phenomena (e.g., in the millisecond to seconds range). The standard time- stepping is used to resolve plasma frequency phenomena using time-steps that are fractions of the plasma period, ensuring accurate representation of high-frequency plasma dynamics.

[0310] In an embodiment, a simulation processor 1028 executes the software systems of FIG.1.

[0311] Referring to FIG. 2, a more general example execution of the software system is provided in the form of a flow chart, and is generally indicated at reference number 2000 (e.g., the method 2000). The method 2000 provides a comprehensive workflow for plasma simulation from initial setup through final analysis and encompasses four primary steps that guide users through the complete simulation process.

[0312] At step 2002, the foundational elements of the plasma simulation are established. During this step, CAD definition is performed to create the geometric representation of the simulation domain, which may include complex three- dimensional structures. Domain setup is implemented to define the computational boundaries and regions where different physics solvers will be applied. Environment setup is conducted to specify the plasma environments (e.g., particle-in-cell environments for kinetic modeling and fluid environments for continuum modeling). Material definitions are established to assign electromagnetic properties, thermal properties, and surface characteristics to different regions of the simulation domain. Boundary conditions are configured to specify voltage sources, current sources, surface 58 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 charging conditions, and natural boundary conditions as described in the electromagnetic solver sections. Emission sources are defined to account for field emission, thermionic emission, and photoemission processes that may occur on conductor surfaces. Meshing parameters are set to control the finite element mesh generation, including element size and shape.

[0313] At step 2004, the configured simulation parameters are prepared for execution by the computational solvers. During this step, any adaptive mesh elements programmatically modify the mesh based on the defined meshing parameters and geometric complexity requirements. The mesh is converted and exported. Simulation definitions are also exported (e.g., environment parameters, material properties, boundary conditions, solver settings).

[0314] At step 2006, the actual plasma simulation calculations are performed using the multi-physics solver framework described throughout the present disclosure. The simulation is initialized and memory is allocated therefore. Time stepping begins using the multi-time-step integration scheme that coordinates the various solvers at their respective optimal time scales. Data is written out at a user-defined frequency in binary format to capture the temporal evolution of electromagnetic fields, particle distributions, fluid properties, and reaction rates. Upon completion of the specified simulation time, the simulation is completed and memory is freed.

[0315] At step 2008, the simulation results are processed and analyzed. Simulation data is imported and processed from the binary output files generated during execution. All major simulation data variables are made available to the user (e.g., electromagnetic potentials and fields, particle densities and velocities, fluid properties such as temperature and pressure, reaction rates and species concentrations, surface charging characteristics, etc.). Data may be presented to the user via three-dimensional visualizations, two-demensional visualizations, or presentation schemes (plots, graphs, tables, etc.).

[0316] The method 2000 accommodates the complex aspectsof plasma systems while maintaining user accessibility through the structured workflow. The method enables simulation of diverse plasma applications including spacecraft charging, plasma- enhanced chemical vapor deposition, discharge phenomena, and vacuum discharge systems through the integration several distinct modeling modalities.

[0317] Referring to FIG. 3, an example CAD design created in the software system is illustrated according to an embodiment. The figure shows a three-dimensional 59 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 cylindrical chamber with a flange attachment displayed within the software's integrated CAD interface (e.g., user interface). The chamber features a transparent outer wall that allows visualization of internal components. The CAD model demonstrates the software system's capability to handle complex three-dimensional geometries that may be encountered in plasma simulation applications such as vacuum chambers, plasma reactors, or discharge vessels. The CAD design may be created internally by the user via the software system of the present disclosure, but may also be imported.

[0318] Referring to FIG.4, an RF signal for voltage boundary conditions displayed by the software system is shown according to an embodiment.

[0319] Referring to FIG. 5, post-processing performed by the software system is illustrated according to an embodiment. The figure shows a contour plot displaying electric potential distributions, with a curved boundary region exhibiting varying potential values represented by different grayscale intensities.

[0320] Referring to FIG. 6, simulation results from an example Gaseous Electronics Conference (GEC) reference cell capacitively coupled plasma (CCP) simulation are shown. Six contour plots are arranged in a 2x3 grid, each showing different aspects of the plasma simulation within a cross-sectional view of the GEC cell geometry. The GEC cell configuration includes two electrodes separated by a distance (e.g., up to 6.35 cm), with peak-to-peak voltages ranging from tens to hundreds of volts applied to generate the plasma discharge. A neutral pressure gradient from high on the left to low on the right results in an asymmetric plasma.

[0321] The top left plot shows the density distribution of Argon neutrals modeled using the fluid solver, exhibiting a four-fold decrease across the chamber when the flow rate is set to 80 cc / s using pressure boundary conditions.

[0322] The top right plot displays the velocity field of the Argon neutrals, showing increased flow speeds in the narrower region between the electrodes.

[0323] The middle left plot presents the electron density distribution generated when 100 V peak-to-peak voltage is applied to the electrodes. The electrons are modeled using the particle-in-cell solver after the background neutral species have reached steady state, creating a spatial imbalance in the electron density that corresponds to the neutral pressure gradient.

[0324] The middle right plot shows the ion density distribution, which follows the electron generation patterns due to the chemical reactions solved through the reaction solver. 60 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0325] The bottom left plot illustrates the electric potential distribution modeled using the finite element method, accounting for space charge contributions from each plasma species and the applied voltage boundary conditions.

[0326] The bottom right plot demonstrates surface reaction modeling through Chemkin integration, showing sputtered Aluminum generated from Argon ion bombardment. The Aluminum species are tracked using the particle-in-cell solver and contribute to the overall plasma composition.

[0327] The simulation from chamber flow dynamics to sputtering phenomena is demonstrated through these six simulation result snapshots, utilizing multi-physics modeling, multi-scale time stepping, and simulation restart capabilities from initial conditions. Beyond reproducing the established results of the GEC Cell reference configuration, the software system enables users to modify chamber parameters including geometry, power settings, flow conditions, and material properties to analyze their effects on surface etching, sputtering, and deposition processes.

[0328] As used herein, the term "solver" is used to refer to a portion of the software system configured to perform calculations in order to enable the simulation of physical phenomena. A solver may be implemented as a computational module, algorithm, or mathematical framework that processes input data and generates output data. The software system may incorporate various types of solvers such as those described above including, but not limited to, a plasma solver, a finite element method (FEM) electromagnetic solver, a particle-in-cell (PIC) solver, a kinetic solver, a fluid solver, a reaction solver, a boundary element method (BEM) solver, a matrix solver (e.g., a BiConjugate Gradient Stabilized (BiCGStab) solver), a hybrid fluid-kinetic solver, a Navier-Stokes equation solver, a Maxwell's equations solver, etc.

[0329] As used herein, “plasma solver” refers generically to any of the aforementioned solvers which, alone or in combination, enable the simulation of plasma discharge and related electromagnetic phenomena. In some aspects, the plasma solver comprises multiple individual solvers working in coordination to predict the behavior of plasma within a defined physical environment, such as combining the finite element method electromagnetic solver with the particle-in-cell solver to model the coupled electromagnetic-kinetic behavior of charged particles in a plasma discharge. In cases wherein only a single solver is used, that solver may be referred to as a plasma solver.

[0330] A simulation engine may be implemented as a computational framework that coordinates the execution of multiple physics solvers within the software system as 61 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 described above. In some aspects, the simulation engine may manage the multi-time- step integration schemes that enable different solvers to operate at their respective optimal temporal scales, such as the accelerated time-steps for cycle time-stepping, fluid time-steps for continuum modeling, and enforced time-steps for plasma frequency resolution. The simulation engine may control the initialization and memory allocation processes, coordinate data exchange between the various solvers (e.g., the finite element method electromagnetic solver, particle-in-cell kinetic solver, fluid dynamics solver, reaction solver, a plasma solver in general, etc.) and manage the temporal evolution of various parameters (e.g., electromagnetic fields, particle distributions, fluid properties, reaction rates, etc.). In general, the simulation engine a computational portion of the software system which coordinates and executes functions thereof.

[0331] One or more functions of the software system (e.g., the plasma solver, the simulation engine) may be executed by a simulation processor (e.g., the simulation processor 1028 of FIG.1). In general, the simulation processor is configured to execute the operations required for plasma simulation (e.g., to execute the plasma solver, to execute the simulation engine). In some aspects, the simulation processor may comprise one or more central processing units (CPUs) or graphics processing units (GPUs). For example, the software system may be configured to be executed on at the CPU of a consumer device (e.g., a computer, a laptop, a desktop, a tablet, a mobile device, a cell phone). In some cases, the simulation processor may comprise distributed computing across multiple nodes (e.g., multiple CPUs, multiple GPUs, multiple CPUs and GPUs, etc.).

[0332] In summary, the software system of the present disclosure includes a plurality of advantageous capabilities / functions, including but not limited to:

[0333] Comprehensive Multiphysics Simulation: The software system provides integrated solvers including full-wave time-domain Electromagnetics (EM), Particle- in-Cell (PIC), Computational Fluid Dynamics (CFD), Bulk and Surface Particle Collision Modeling (Reaction), and Hybridized PIC / CFD solutions.

[0334] Advanced Numerical Features: The software system incorporates analytic and 1D-PIC sheath models for accurate boundary condition simulation. It employs semi- implicit, exactly energy-conserving time-stepping methods for PIC, along with RF- cycling time-stepping for plasma frequencies and long-timestepping for fluid and circuit solvers (ms to seconds). The system supports both structured and unstructured meshes and provides Power, Voltage, and Current excitation models. 62 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0335] Parallelization and Scalability: The software system features full parallelization across multiple nodes and is compatible with various operating system environments. It provides simulation flexibility with the ability to restart simulations from previous results and supports an unlimited number of species and reactions in simulations.

[0336] 3D CAD Integration and Workflow Efficiency: The software system is fully integrated with a three-dimensional (3D) computer aided engineering (CAE) user interface (UI), supporting all CAD formats. It includes built-in post-processing for 3D analysis of EM, PIC, and fluid results, including species temperature, density, momentum, reaction rates, potentials, fields, surface deposition, and more. Data can be exported to common unstructured grid formats (e.g., .vtd files) for further analysis in tools like ParaView.

[0337] Coupling with External Tools: The software system provides interoperability with the Ansys Suite, including coupling with Ansys Maxwell for magnet simulations, integration with Chemkin for plasma libraries, and access to Granta for materials libraries. It also features advanced field imports from Ansys HFSS and Ansys EMC Plus.

[0338] Advanced Boundary and Physical Models: The software system includes field emission, secondary emission, thermionic emission, and space charge limiting emission models as boundary conditions. It incorporates auxiliary circuit boundary conditions integrated with a full SPICE circuit modeler.

[0339] Sophisticated Reaction Modeling: The software system provides comprehensive reaction models with advanced gas-gas and gas-solid interaction models for elastic scattering, ionization, excited states, etching, and deposition. It can simulate secondary particle generation at surfaces and maintains a growing gas reaction library to support an ever-increasing number of reactions and species.

[0340] While the systems and methods above have been described and disclosed in certain terms and have disclosed certain embodiments or modifications, persons skilled in the art who have acquainted themselves with the disclosure, will appreciate that it is not necessarily limited by such terms, nor to the specific embodiments and modification disclosed herein. Thus, a wide variety of alternatives, suggested by the teachings herein, can be practiced without departing from the spirit of the disclosure, and rights to such alternatives are particularly reserved and considered within the scope of the disclosure. 63 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302

[0341] When introducing elements, the articles “a,” “an,” “the,” and “said” are intended to mean that there are one or more of the elements. The terms “comprising,” “including,” and “having” are intended to be inclusive and mean that there may be additional elements other than the listed elements.

[0342] Not all of the depicted components illustrated or described may be required. In addition, some implementations and embodiments may include additional components. Variations in the arrangement and type of the components may be made without departing from the spirit or scope of the claims as set forth herein. Additional, different or fewer components may be provided and components may be combined. Alternatively, or in addition, a component may be implemented by several components.

[0343] The above description illustrates embodiments by way of example and not by way of limitation. This description enables one skilled in the art to make and use aspects of the disclosure, and describes several embodiments, adaptations, variations, alternatives and uses of the aspects of the disclosure, including what is presently believed to be the best mode of carrying out the aspects of the disclosure. Additionally, it is to be understood that the aspects of the disclosure are not limited in its application to the details of construction and the arrangement of components set forth in the following description or illustrated in the drawings. The aspects of the disclosure are capable of other embodiments and of being practiced or carried out in various ways. Also, it will be understood that the phraseology and terminology used herein is for the purpose of description and should not be regarded as limiting.

[0344] It will be apparent that modifications and variations are possible without departing from the scope of the disclosure defined in the appended claims. As various changes could be made in the above constructions and methods without departing from the scope of the disclosure, it is intended that all matter contained in the above description and shown in the accompanying drawings shall be interpreted as illustrative and not in a limiting sense.

[0345] In view of the above, it will be seen that several advantages of the aspects of the disclosure are achieved and other advantageous results attained.

[0346] The Abstract and Summary are provided to help the reader quickly ascertain the nature of the technical disclosure. They are submitted with the understanding that they will not be used to interpret or limit the scope or meaning of the claims. The Summary is provided to introduce a selection of concepts in simplified form that are further 64 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 described in the Detailed Description. The Summary is not intended to identify key features or essential features of the claimed subject matter, nor is it intended to be used as an aid in determining the claimed subject matter. 65 CORE / 3528169.000303 / 232777312.1

Claims

Docket No.3528169.000302 WHAT IS CLAIMED IS:

1. A non-transitory computer-readable storage medium storing instructions executable by a processor to simulate plasma discharge, the instructions comprising; a plasma solver for predicting the behavior of a plasma discharge; a 3D CAD system having an interface configured to receive user input relating to a physical environment in which the plasma discharge is predicted to occur, the 3D CAD system configured to generate a model of the physical environment in response to the user input; and a simulation engine coupled to the plasma solver and the 3D CAD system, the simulation engine configured to simulate the plasma discharge within the physical environment based on the predicted behavior from the plasma solver.

2. The non-transitory computer-readable storage medium of claim 1, wherein the plasma solver comprises two or more of: a simultaneous charge conserving and energy conserving semi- implicit timestepping method; a full-wave electromagnetic equation solver; a dynamic integration of kinetic solver; or a hybrid fluid-kinetic solver.

3. The non-transitory computer-readable storage medium of claim 1, wherein the 3D CAD system is further configured to render the simulation.

4. The non-transitory computer-readable storage medium of claim 1, wherein the plasma solver comprises; a particle-in-cell solver; and a finite element method electromagnetic solver configured to interact with the particle-in-cell solver, the finite element method electromagnetic solver responsive to the user input for deriving a plurality of boundary conditions.

5. The non-transitory computer-readable storage medium of claim 1, wherein the plasma solver comprises a finite element method electromagnetic solver configured to solve Maxwell's equations using basis functions and Galerkin's method.

6. The non-transitory computer-readable storage medium of claim 1, wherein the plasma solver comprises a reaction solver configured to model interactions between multiple plasma species using at least one of analytic formalism or cross-section formalism. 66 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 7. The non-transitory computer-readable storage medium of claim 1, wherein the simulation engine is configured to operate one or more time-stepping schemes including standard time-stepping, fluid time-stepping, or cycle time-stepping with different temporal scales.

8. The non-transitory computer-readable storage medium of claim 1, wherein the plasma solver is configured to model low-temperature plasma applications.

9. The non-transitory computer-readable storage medium of claim 1, wherein the plasma solver is configured to model plasma-enhanced chemical vapor deposition.

10. The non-transitory computer-readable storage medium of claim 1, wherein the plasma solver comprises a fluid solver configured to solve Navier-Stokes equations for fluid density, velocity, and temperature using finite element method discretization.

11. A method of simulating plasma in a simulation software, comprising: defining a physical environment using a 3D CAD module; enumerating a plurality of physical parameters and initial conditions corresponding to the defined physical environment; executing a multi-physics solver framework comprising a finite element method electromagnetic solver and a particle-in-cell solver; and predicting a behavior of plasma within the defined physical environment based on the executed multi-physics solver framework.

12. The method of claim 11, wherein the particle-in-cell solver is exactly energy- conserving.

13. The method of claim 11, wherein the plasma is low-temperature plasma.

14. The method of claim 11, further comprising representing the behavior of the plasma within the defined physical environment as data in an unstructured grid format.

15. The method of claim 14, further comprising exporting said data.

16. The method of claim 11, further comprising solving electromagnetic fields using a finite element method solver that applies natural boundary conditions and gauge constraints.

17. The method of claim 11, further comprising modeling chemical reactions between plasma species using reaction rates.

18. The method of claim 11, further comprising coordinating multiple computational physics solvers using multi-time-step integration schemes with different temporal scales for electromagnetic, kinetic, and fluid phenomena.

19. The method of claim 11, further comprising visualizing simulation results through a three-dimensional rendering within the 3D CAD module interface. 67 CORE / 3528169.000303 / 232777312.1Docket No.3528169.000302 20. A plasma simulation system comprising: a simulation processor; a memory comprising a non-transitory computer-readable storage medium, the memory being coupled to the simulation processor; a multi-physics solver framework stored in the memory and executable by the simulation processor, the framework comprising: a finite element method electromagnetic solver; a particle-in-cell kinetic solver; a fluid dynamics solver; and a reaction solver; and an integrated 3D CAD environment stored in the memory and executable by the simulation processor, the environment configured to enable a user to interact with the multi-physics solver framework. 68 CORE / 3528169.000303 / 232777312.1

Images