含可变温度T的电子、中性及离子粒子连续性方程组求解方法咨询

Hey there! Let's walk through how to tackle this system of continuity equations for electrons ($e$) and neutral/ionic particles ($p$), especially since variable temperature $T$ adds a layer of complexity to both diffusion and reaction terms.

First, let's restate your equations clearly for reference:

$$
\frac{\partial n_{e,p}}{\partial t}+\nabla \cdot J_{e,p}=S_{e,p} \tag{1}
$$

$$
J_{e,p}=-\nabla n_{e,p} \cdot D(T)_{e,p} \tag{2}
$$

$$
S_{p,e}=\sum k(T){i}n{p}n_{e} \tag{3}
$$

Step 1: Reduce the System to Coupled Reaction-Diffusion Equations #

Start by substituting the flux term (equation 2) into the continuity equation (equation 1). This gives us a pair of coupled partial differential equations (PDEs):

$$
\frac{\partial n_e}{\partial t} = \nabla \cdot \left(D_e(T) \nabla n_e\right) + \sum k_i(T) n_p n_e
$$

$$
\frac{\partial n_p}{\partial t} = \nabla \cdot \left(D_p(T) \nabla n_p\right) + \sum k_i(T) n_p n_e
$$

Quick note: Double-check the sign and specific reaction rates for each species in $S_{p,e}$! For example, recombination reactions would decrease both $n_e$ and $n_p$, while ionization might increase one and decrease the other—getting these signs right is critical for physical accuracy.

Step 2: Handle the Variable Temperature $T$ #

The temperature dependence of $D(T)$ (diffusion coefficients) and $k(T)$ (reaction rate constants) is the key challenge here. You’ll fall into one of two scenarios:

• Case 1: $T(x,t)$ is known: If you have a precomputed temperature field (from experiments or a separate simulation), treat $D_e(T)$, $D_p(T)$, and $k_i(T)$ as spatiotemporally varying coefficients. Just plug in the corresponding $T$ values at each grid point and time step during solving.
• Case 2: $T(x,t)$ is unknown: You’ll need to couple these continuity equations to energy conservation equations for electrons and ions/neutral particles. Temperature affects diffusion and reactions, and reactions (like recombination) release/absorb energy that changes $T$. This turns the problem into a multi-physics coupling task—solve the continuity and energy equations iteratively until convergence.

Step 3: Choose a Solution Method #

Analytical solutions are only feasible for extremely simple cases (e.g., steady-state, 1D, constant $T$, trivial boundary conditions). For most real-world scenarios, go with a numerical approach:

• Finite Difference Method (FDM): Discretize space and time into grids. Use explicit schemes (forward Euler for time, central difference for space) for easy implementation, or implicit schemes (Crank-Nicolson) for better stability (though they require solving linear systems at each step).
• Finite Element Method (FEM): Ideal for complex geometries. Convert PDEs to their weak form, then use basis functions to approximate solutions. FEM handles irregular boundaries smoothly and integrates variable coefficients easily.
• Finite Volume Method (FVM): A conservative scheme popular in fluid dynamics and plasma physics. Discretize the domain into control volumes to enforce local mass conservation—perfect for problems where conservation properties are non-negotiable.

Step 4: Tackle Nonlinearity #

The source term $\sum k_i(T) n_p n_e$ is nonlinear (it’s a product of two unknowns). When solving numerically, you have two common approaches:

• Linearization: Use $n_e$ and $n_p$ values from the previous time step to approximate the product term, turning the equation linear for the current step. Simple to implement, but may require smaller time steps for stability.
• Newton-Raphson Method: For strong nonlinearities, use iterative Newton updates to solve the nonlinear system at each time step. Converges faster, but requires computing Jacobian matrices.

Step 5: Define Initial & Boundary Conditions #

PDEs can’t be solved without these:

• Initial conditions: Specify $n_e(x,0)$ and $n_p(x,0)$ for all points in your domain.
• Boundary conditions: Common options include:
  • Dirichlet: Fixed concentration at the boundary ($n_e = n_{e,0}$)
  • Neumann: Fixed flux (e.g., insulated boundary where $\nabla n_e \cdot \mathbf{n} = 0$, meaning no particle flow across the boundary)
  • Robin: A mix of concentration and flux (e.g., particle loss to a wall with a finite reaction rate)

Example Simplification: Steady-State, Constant $T$ #

If you’re in steady state ($\partial/\partial t = 0$) and $T$ is constant, the equations reduce to:

$$
0 = D_e \nabla^2 n_e + k n_p n_e
$$

$$
0 = D_p \nabla^2 n_p + k n_p n_e
$$

In 1D, this might be solvable with analytical methods (e.g., symmetry assumptions or separation of variables), but the nonlinear source term still makes it tricky—numerical methods are often safer here too.

内容的提问来源于stack exchange，提问作者Sirine