# A Study on the Efficient Method for the Solution of the Saha Equation in the Numerical Simulation of Capillary Discharge

Authors: Hengzhu Liu; Cheng Wu; Fengjiang An; Shasha Liao
Journal: Advances in Mathematical Physics
 (2017)

---

## Abstract

In the numerical simulation of capillary discharge for the Electrothermal (ET) or ET-chemical (ETC) launch system, one needs to solve the Saha equation to determine the composition of plasma for the transport coefficients. This would cause a relatively longer simulation runtime compared with conventional computational fluid dynamics problems. In this paper, the relatively difficult multidimensional problem of solving the Saha equation is transformed into a one-dimensional problem in the simulation of capillary discharge by constructing an iteration equation about temperature. A coefficient is introduced to ensure the convergence of this iteration equation. In order to improve the computational efficiency, this coefficient is further optimized and the methods of setting the iteration initial value dynamically are introduced. Several simulation tests are conducted to study the performance of these two methods. The results show that the simulation runtime could be significantly reduced with the methods presented in this paper.


---

## Body

## 1. Introduction

Electrothermal (ET) or ET-chemical (ETC) launch technology can achieve a muzzle velocity of more than 2 km/s [1]. It is the hypervelocity that conventional gunpowder driven launch technology can be barely achieved. As key physical processes in the ET or ETC launch system, the research in this field is always focused on the capillary discharge which generates high pressure and high temperature plasma by ablating the wall of the tube to serve as working medium or ignite the propellant [2]. During the last several decades, various mathematical models with different degree of simplification were established based on the Euler equations to study this physical process [3–9]. As one can find from these models, a characteristic is that those driver terms or transport coefficients, such as conductivity and ablation rate, are all related to the plasma composition, so that one needs to determine the composition of plasma in each computational cell firstly before solving the mathematical model for the plasma properties of next time step.In the ET or ETC launch technology, the typical range of temperature and pressure of the plasma generated by capillary discharge are 10000–40000 K and 100 MPa, respectively [10]. At such extreme conditions, plasmas are assumed to be in local thermodynamic equilibrium (LTE) and ablated material is fully dissociated into their elemental constituents and partially or fully ionized [11]. With these assumptions, the composition of plasma could be determined by solving the Saha equation. The Saha equation defines the relation of ions in two adjacent ionization stages. It is usually a multidimensional problem if one solves it directly. As mentioned by Zaghloul [12], this would be rather difficult to find the solution. In the work of Trayner and Głowacki [13] the Saha equation for the situation of ideal plasma with single element species was firstly reformulated into a one-dimensional equation and then solved with some numerical algorithm such as the bisection method. Later, Zaghloul extended this method for the situation of nonideal plasma with multielement species while maintaining the one dimensionality of the method [12]. The methods of Trayner and Głowacki and Zaghloul are prominent, since their methods enable the work of solving the Saha equation individually to become very easy.For the simulation of capillary discharge used in the ET or ETC launch system, a problem of employing the method of Trayner and Głowacki or Zaghloul is that those conversation equations in the mathematical model of capillary discharge cannot give temperature directly whereas it is required as one of the inputting data pieces for solving the Saha equation. In the series researches of Zoler et al. [5, 14], the energy conversation equation is transformed through the variable substitution which makes it able to give temperature directly. However, couples of partial differential terms would also be introduced into the model which would cause extra computation amount for calculating these terms. In other works [4, 9, 15], the Saha equation was directly solved with the specific energy equation together in the simulation. In this paper, this direct ideal is retained. Based on the methods of Trayner and Głowacki and Zaghloul, the multidimensional problem of solving the Saha equation in the simulation of capillary discharge is transformed into a one-dimensional problem by modifying the specific energy equation into an iteration equation about temperature.This method enables the work of solving the Saha equation to become rather simple in the simulation of capillary discharge; however, one still needs to solve an iteration equation with some numerical algorithm in each cell. Due to the computational amount, the simulation of capillary discharge usually takes longer time compared with conventional computational fluid dynamics problem. Particularly, in the cases that a long tube might be joint with the exit of the capillary to serve as the barrel or cathode or other purposes [4, 16, 17], the larger computational region would increase the simulation runtime further. In this paper, two methods are introduced for a better computational efficiency including optimizing the coefficient introduced in the iteration equation and setting the iteration initial value dynamically. The mechanism of accelerating the computation with these methods is also numerically investigated.


## 2. The Saha Equation

In the simulation research of capillary discharge for the ET or ETC launch technology, a widely used form of the Saha equation is given by [15](1)ψji+1ψji=Kji,where ψji is the ratio of the i-fold ionized ions of the elemental species j to the total heavy particles [12]. Apparently, determining the composition of plasma by solving the Saha equation actually means determining ψji by solving (1). The function of Kji is given by [12, 15, 18](2)Kji=1nψe2Uji+1Uji2πmekTh23/2exp⁡-Iji-ΔIjikT,where n is the total number density of the heavy particles, ψe is the ratio of the electron to the total heavy particles, Uji denotes the partition functions which are dependent on the temperature T, me is the mass of the electron, k is Boltzmann’s constant, h is Planck’s constant, Iji is the ionization potential. For the value of Iji and necessary spectroscopic data used to calculate Uji, one can find them from the database online [19]. ΔIji is the lowering of the ionization potential which is introduced since the nonideal effect of plasma. ΔIji is given by [4, 12](3)ΔIji=i+1e24πε0λD+ΛB/8,where ε0 is the permittivity of free space, e is the electronic charge, ΛB is the de Broglie wave length, and λD is the standard Debye length which could be expressed as(4)λD=ε0kTe2nψe+ψa,where ψa is defined as(5)ψa=∑j∑i=1zji2ψji,where zj stands for the maximum allowed ionization stage of the elemental species j. zj could be estimated by the possible upper limit of the temperature in a specific simulation case. Two supplementary equations for solving the Saha equation are given by(6)ψe=∑j∑i=1zjiψji,(7)ψj=∑i=0zjψji,where (6) is the condition of electroneutrality of plasma. ψj is the ratio of the elemental species j to the total heavy particles. ψj could be obtained from the density of plasma with the chemical formula of the material ablated from the wall of the capillary. For solving the Saha equation expressed as (1), the density of plasma, whether the number density n or the mass density ρ, and T are required as the inputting data. The above equations are used to illustrate the method of solving the Saha equation in this paper. For more details about the Saha equations, one can find them in [4, 12, 15].


## 3. Method of Solving the Saha Equation Individually

The method of solving the Saha equation individually listed below is from the work of Trayner and Głowacki and Zaghloul [12, 13]: for the element species j, (1) could be expressed as(8)ψj1=ψj0·Kj0,ψj2=ψj1·Kj1=ψj0·Kj0·Kj1,ψj3=ψj2·Kj2=ψj0·Kj0·Kj1·Kj2,⋮.Namely, ψji  i=1,2,…,zj could be expressed as(9)ψji=ψj0∏k=0i-1Kjk.From (7), ψj0 can be expressed as(10)ψj0=ψj-∑i=1zjψji.By substituting (9) into (10) with some rearrangements, one gets(11)ψj0=ψj1+∑i=1zj∏k=0i-1Kjk.Firstly, for an easier case that ΔIji=0, which is the ideal plasma, only ψe in Kjk is unknown as one can find from (2), so that (11) could be expressed as(12)ψj0=fj0ψe.Equation (12) indicates that ψj0 could be expressed explicitly in terms of ψe. By substituting (12) into (9), one gets(13)ψji=fj0ψe∏k=0i-1Kjk,where i=1,2,…,zj. Similarly, as (11), only ψe in the right side of (13) is unknown which indicates that ψji  i=1,2,…,zj can also be expressed explicitly in terms of ψe. Then, together with (12), all ψji  i=0,2,…,zj of the element species j could be expressed explicitly in terms of ψe. Namely,(14)ψji=fjiψe.Since (14) is actually effective for any element speciesj considered in a simulation case, (14) also means that ψji  i=0,2,…,zj of each element species j considered in the simulation could be expressed explicitly in terms of ψe. Then, by substituting (14) into (6), one could get an iteration equation about ψe. Namely,(15)ψe=∑j∑i=1zjifjiψe.In this paper, the simple fixed-point iteration method [20] is employed to solve the iteration equation as (15). The procedure is as follows:(i)
Setψe to be the iterative initial value.(ii)
Calculate (11) to obtain ψj0.(iii)
Calculate (9) to obtain ψji.(iv)
Repeat steps (ii)-(iii) for all element species considered in the simulation case.(v)
Calculate (6) to obtain a new ψe denoted as ψe_re.(vi)
Setψe to be ψe_re and repeat steps (ii)–(v) till the termination criteria are satisfied.As one can find from procedures (i)–(vi),ψji could be determined simultaneously with ψe. Namely, for the case of ΔIji=0, a relatively difficult multidimensional problem is reduced into a one-dimensional problem. For the present case that ΔIji≠0, as can be seen from (2)–(5), ψa cannot be expressed explicitly in terms of ψe unless zj is 1, so that both ψe and ψa are unknown in Kjk. In this case, (12) and (14) become(16)ψj0=gj0ψe,ψa,ψji=gjiψe,ψa.Then, (15) becomes(17)ψe=∑j∑i=1zjigjiψe,ψa.Apparently, (17) cannot be solved with procedures (i)–(vi) directly. However, if ψa is artificially set to be an empirical value, ψe and ψji might still be obtained by solving (17) with procedures (i)–(vi) as solving (15). Then, a new ψa could be obtained by calculating (5). This is actually an iteration process about ψa which could be expressed as(18)ψa=g¯ψa.The bar over g refers to the numerical nature of the function [12]. The procedure of solving (18) is as follows:(a)
Setψa to be the iterative initial value.(b)
Solve (17) with procedures (i)–(vi) to obtain ψji.(c)
Calculate (5) to obtain a new ψa denoted as ψa_re.(d)
Setψa to be ψa_re and repeat steps (b)-(c) till the termination criteria are satisfied.As mentioned by Zaghloul, the above method transforms the relatively difficult multidimensional problem into a one-dimensional problem aboutψa which is easy to find its solution to any desired accuracy [12].


## 4. Solving the Saha Equation in the Simulation of Capillary Discharge

(19)
∂
ρ
∂
t
+
w
∂
ρ
∂
z
+
ρ
∂
w
∂
z
=
ρ
α
′
,
(20)
ρ
∂
w
∂
t
+
ρ
w
∂
w
∂
z
+
∂
p
∂
z
=
-
w
ρ
α
′
,
(21)
ρ
∂
E
∂
t
+
ρ
w
∂
E
∂
z
+
p
∂
w
∂
z
=
J
2
σ
-
ρ
α
′
E
-
w
2
2
,
(22)
E
=
E
t
r
+
E
e
l
+
E
c
h
+
E
v
+
E
d
,
(23)
p
=
n
+
n
e
k
T
.Equations (19)–(23) are the time-dependent 1D mathematical model used to simulate capillary discharge for the ET or ETC launch technology [4, 11, 21]. In this set of models, t stands for time and z stands for axial coordinate, ρ is the mass density, w is the velocity of the plasma flow, ρα′ is the ablation rate which causes the increasing of ρ, p is the pressure, E is the specific internal energy (SIE), σ is the conductivity, and J2/σ is the inputting energy in the form of joule heating. σ is determined based on the electron collisions with both neutral atoms and ions in which the collision frequency with ions is evaluated with the model of Zollweg and Liebermann model [22]. For detailed theories or methods of calculating σ as well as transport coefficient such as ρα′ and those terms in the SIE equation (22), one can find them from [4, 15]. Again, as mentioned before, these transport coefficients are all related on ψji, so that one needs to solve the Saha equation in each computational cell to determine ψji firstly, and, then, (19)–(23) could be solved for the plasma flow properties of next time step. However, (19)–(21) cannot give T directly whereas T is required as one of the inputting data pieces for solving the Saha equation. Similarly, as constructing (17) for the new emerged ψa in the case of nonideal plasma mentioned before, the method of Trayner and Głowacki and Zaghloul could be further extended by constructing an iteration equation about T while maintaining the one dimensionality of the method. The iteration equation about T could be easily obtained by modifying (22) into(24)T=Ein-ET,ρin·δ+T,where Ein and ρin represent the SIE and density determined by (19)–(21). These two could be viewed as inputting data for solving (24). ET,ρin stands for the result by calculating (22). The procedure of solving (24) is as follows:(A)
SetT to be the iterative initial value.(B)
Solve (18) with procedures (a)–(d) to obtain ψji.(C)
Calculate (22) to obtain ET,ρin.(D)
Calculate (24) to obtain a new T denoted as Tre.(E)
SetT to be Tre and repeat steps (B)-(C) till the termination criteria are satisfied.By the above method, the relatively difficult multidimensional problem of solving the Saha equation in the simulation of capillary discharge is transformed into a one-dimensional problem aboutT. In (24), a coefficient δ is introduced. It is used to enable iteration equation (24) to satisfy the condition of convergence. For (24), this convergence condition could be expressed as [20](25)h′T<1,where h(T) represents the right side of (24). From inequation (25), one can obtain the range of δ, which is(26)0<δ<2∂E/∂T.In the capillary discharge for the ET or ETC launch technology, the magnitude order of ∂E/∂T is usually 1×104 (in SI units). It is obvious that (24) could not satisfy the convergence condition if δ is not introduced or set to be 1. For the specific value of δ in a single cell, one can set it for maximum computational efficiency. From inequation (25), one can find that if δ enables h′Tt=0, where Tt represents the true root in a cell, it might need least iterations to find the root for this cell. Namely,(27)δ=∂E∂TTt-1.An ideal case is that if h(T) is a linear function and δ is set by (27), only one time of iteration would be needed to find the solution and the convergence rate would be infinitely great. For the nonlinear function, if δ enables h′Tt=0 in a cell, it means that, with T getting closer to Tt, h′T would be more and more closer to h′Tt which is 0 and the convergence rate would also increase with this process. It is more like an acceleration effect of the convergence rate. Unlike the Newton-Raphson method which is quadratically convergent, the simple fixed-point iteration method itself is linearly convergent [20], yet one might still obtain some acceleration effect by optimizing δ. Namely, if δ is set to be closer to ∂E/∂T-1 for a cell, less iterations would be needed to find the root.In the actual simulation of capillary discharge, it would be tedious and unnecessary to setδ for each cell at different time point. A more practical way is setting δ to be a unified value for the whole simulation as long as the condition of convergence could be satisfied. Based on the above analysis, one can find that, for the minimum runtime, this unified δ should be chosen to be very close to the average value of ∂E/∂T-1 over all the cells at every time point. Since one cannot obtain ∂E/∂T-1 nor its average beforehand, δ could be set to be a relatively small value to ensure converging as one sees from inequation (26) for the first time of simulation; then, one can obtain the average value of ∂E/∂T-1 from the simulation result. Since δ is set as a unified value, the exact value of δ for the minimum runtime might also be affected by other factors such as the method of setting the iterative initial value. One may need to adjust δ based on this average value of ∂E/∂T-1 for extra rounds of simulation tests till the minimum simulation runtime is achieved.


## 5. Setting Iterative Initial Value

The basic principle of setting the iterative initial value is to set it as close as possible to the true root so that less iteration would be needed to find it. Besides, sinceh(T) is nonlinear function, the above method of setting δ also indicates that if the iterative initial value is too far away from the true root, the iteration equation cannot satisfy the condition of convergence. In the time-dependent simulation of capillary discharge, the plasma flow properties, including ψe, ψa, and T, are dynamic spatial-temporal data, so that one may also need to set the iterative initial value of ψe, ψa, and T dynamically. Otherwise, with the changing of inputting current, ψe, ψa, or T in a cell at present time point might change greatly from these properties in this cell at a previous time point. This might also cause large amount of iteration to find the solution or even make the iteration diverging.In the simulation of capillary discharge, in order to achieve a high quality solution or satisfy the Courant-Friedrichs-Lewy (CFL) condition to maintain the stability of the simulation, space step or the time step is usually set to be relatively small value, so that plasma properties from two spatial or temporal adjacent cells might be very close to each other. This actually provides two easy approaches to set the iterative initial value dynamically. As shown in Figure1, for solving the Saha equation or (24) in “present cell,” one could set the iterative initial value of ψe, ψa, or T to be the solution of (24) correspondingly from “last time step” or “last cell.” In this paper, for convenience, setting the iterative initial value from the solution of “last time step” is named as “Initial-T” and setting the iterative initial value from the solution of “last cell” is named as “Initial-S”; since there is no “last cell” for the cell in the closed end of the capillary, “Initial-T” is still used for this one cell in the method of “Initial-S.” Additionally, for comparison research, setting the iterative initial value fixedly is named as “Initial-F.” In this paper, the fixed iterative initial value of ψe, ψa, and T is set to be their average over all the cells at every time point from a previous simulated result.Figure 1
Schematic diagram of the method of setting iterative initial value.


## 6. Results and Discussion

Several numerical tests are conducted to study the performance of methods presented in this paper. For the accuracy of the present methods, two tests are conducted. The first one is for procedures (a)–(d). In this test, the Saha equation or (18) for polyethylene plasma is solved individually by employing procedures (a)–(d) with different pressure and temperature. In this paper, the terminal criteria for solving the Saha equation are constructed with the approximate percent relative error εa [20]. For T, the terminal criteria could be expressed as(28)εa=Tcurrent-TpreviousTcurrent<5×10-6.By this terminal criteria, 5 significant digits could be achieved and this is found to be sufficiently accurate for the present research. Then, with the solution, the conductivity is calculated and compared with the results of Kim [11]. In Figure 2, calculated conductivity for the polyethylene plasma as a function of temperature and pressure from present results and [11] are plotted. As can be observed, present results are found to agree perfectly with the results of Kim [11] which could prove the accuracy of procedures (a)–(d) as well as the correctness of the coding work.Figure 2
Calculated conductivity for the polyethylene plasma as a function of temperature and pressure. Solid line represents the results from present calculations and symbols represent the results from [11].The second test is for procedures (A)–(D). In this test, the capillary discharge used in the ET or ETC launch system is numerically simulated with the experimental data from shot 8 of series experiments conducted by Powell and Zielinski [4] and procedures (A)–(D) are employed to solve the Saha equation in each cell in the simulation. Also, the simulated plasma flow properties are compared with the simulated results presented in the work of Powell and Zielinski [4]. In Figure 3, the inputting current and the geometry of the capillary used in the simulations are plotted. The inputting current data in Figure 3(a) is collected from [11] with image recognition software and then is fitted to a ninth-order polynomial in time. The simulation is conducted from 200 μs to 700 μs to avoid the initial period dominated by the electrical explosion used to initialize the arc [4, 11]. The mathematical model used for the simulation is (19)–(23). The Lax scheme is employed to solve this model. The grid points and CFL number are fixed to be 100 and 0.5, respectively. The method of setting the iterative initial value is “Initial-T.” Results from simulation with “Initial-S” or “Initial-F” can barely be distinguished from each other. In Figure 4, simulated plasma flow properties along the axial direction at 500 μs are plotted. As can be observed, simulated plasma flow properties from the present simulation are a little higher than those of Powell and Zielinski [4]. One of the possible reasons could be that the present inputting current collected with image recognition software is not exactly the same as that used in the simulation of Powell and Zielinski. Other possible reasons for these minor differences could be different setting of grid points or time steps, different scheme used for the simulation, and so forth. For more detailed information about the simulation and the discussion about the simulated plasma flow properties, one can found them in [4, 11, 21].Figure 3
Inputting current profile and the geometry of the capillary discharge [4].
(a)
(b)Figure 4
Simulated plasma flow properties along the axial direction at 500µs. L is the total length of computational region (L=60.9 mm). Solid lines represent the results from present calculations and symbols represent the results from [4].Based on the above simulation, the efficiency of present methods is studied. In Figure5, simulation runtime with δ varied from 0.4 to 0.8 and different methods of setting the iterative initial value are plotted in order to select δ for the minimum runtime. The simulation program is coded with C++ and run on Windows 7 with the CPU of Intel Core i7-6700HQ. Most of those file writing functions are turned off in the simulation program for efficiency. Each final recorded runtime is the average runtime over 50 times repeated simulation. As can be observed, there is a significant advantage of using “Initial-S” or “Initial-T” over “Initial-F.” The runtime with “Initial-T” and “Initial-S” is approximately arranged from 9 to 13 sec and 20 to 29 sec, respectively, whereas the runtime with “Initial-F” is approximately arranged from 45 to 69 sec. In Figure 5, it is also found that even if δ is set to be 0.7×10-4 for “Initial-T,” the total runtime is still approximately 50% less than the minimum runtime of the simulation with “Initial-S.” In the methods of “Initial-S” and “Initial-T,” ∂ζ/∂z×Δz and ∂ζ/∂t×Δt, where ζ could be ψe, ψa, T, or other properties, actually mean the distance of the iterative initial value to the true root, so that these two methods connect ∂ζ/∂z×Δz and ∂ζ/∂t×Δt with the computational efficiency. It is probably that ∂ζ/∂z×Δz is generally larger than ∂ζ/∂t×Δt. From the simulation results, it is found that the magnitude order of ∂T/∂t is 1×107 and the magnitude order of Δt is 1×10-8, so that the magnitude order of ∂T/∂t×Δt is approximately 1 K. For the magnitude order of ∂T/∂z, it is varied from 0 at the closed end of the capillary to 3×105 at the exit. Since Δz is 60.9×10-5, the ∂T/∂z×Δz could reach approximately 180 K at most.Figure 5
Simulation runtime with differentδ and different method of setting the iterative initial value.In Figure6, the spatial-temporal distribution of the loop count of procedures (i)–(vi) from the simulations with “Initial-S” and “Initial-T” are plotted, respectively. Since the structure of procedures (A)–(E) is actually a three layers’ loop and procedures (i)–(vi) are the innermost loop, the loop count of procedures (i)–(vi) could indicate the total iterations in one cell. The value of δ used in the simulations is set from Figure 5 for minimum runtime. In Figure 6(a), since ∂ζ/∂z always increases along the axial direction in the capillary discharge, ∂ζ/∂z×Δz would be larger in the region near the exit of the capillary. This causes that the iterative initial value obtained from “last cell” becomes further to the true root along the axial direction and the loop count would also increase and reach its highest level in the region near the exit of the capillary. In Figure 6(b), the distribution is generally consistent with ∂I/∂t in which I is the inputting current. In the early and late period of the simulation, current increases and decreases rapidly as can be observed from Figure 3(a); namely, ∂I/∂t is in a relatively high level in these two periods. This makes the plasma flow properties also changing rapidly, so that the iterative initial value obtained from “last time step” might be relatively far away from the true root in the early and late period of the simulation. In the middle period from 430 μs to 530 μs when current is near its peak, ∂I/∂t is about to be 0, so that the plasma flow properties would remain relatively constant in this period. This makes that the iterative initial value obtained for a present cell from “last time step” becomes very close to the true root, so that the loop count could reach its lowest level in the middle period.Figure 6
Spatial-temporal distribution of the loop count of procedures (i)–(vi) with (a) “Initial-S” andδ=0.7×10-4. (b) “Initial-T” and δ=0.54×10-4.
(a)
(b)A similarity in Figures6(a) and 6(b) is that all the loop count can reach its own peak near the exit of the capillary in the start of the simulation. As mentioned before, for a single cell, the solution of (24) could be found with least iterations by setting δ as close as possible to the true ∂E/∂T-1 of this cell. In this paper, δ is set to be a unified value for the whole simulation, it might happen that the unified δ is relatively far away from their true ∂E/∂T-1 in some spatial-temporal regions. This causes that the convergence rate might remain in a relatively low level in the iteration process which leads to a relatively high level of the loop count in these spatial-temporal regions. In Figure 7, the spatial-temporal distribution of ∂E/∂T is plotted. It is found that ∂E/∂T could reach its peak of approximately 2.14×104 at the exit of the capillary in the start of the simulation and reaches its lowest level of approximately 1.3×104 at the closed end of the capillary near 430 μs to 530 μs when current also reaches its peak. From the simulation results, it is found that the average of ∂E/∂T-1 over all the cells at every time point is 0.69921×10-4 which is closer to the side of 1.3×104-1, and, then, the side of 2.14×104-1 is the furthest point to the average of ∂E/∂T-1. For the reason that ∂E/∂T could be in its relatively high level near the exit in the start or late period of the capillary discharge, one can find it from the definition of SIE or by analyzing the dimension of SIE.Figure 7
Spatial-temporal distribution of∂E/∂T.As mentioned above, the exact value ofδ when minimum runtime is achieved might also be affected by the method of setting the iterative initial value. In Figure 5, when minimum runtime is achieved, δ is 0.54×10-4 for “Initial-T” whereas δ is 0.7×10-4 for “Initial-S.” The former one is not as close as that of “Initial-S” to the average of ∂E/∂T-1 which is 0.69921×10-4 as mentioned above. An approximate explanation for this is as follows: in Figure 6(b), the loop count is in its higher level and is in the early and late period of the simulation due to the profile of the inputting current. Meanwhile, in these two periods, ∂E/∂T is also in its relatively high level as shown in Figure 7, so that δ could be adjusted a little to the lower side of ∂E/∂T-1, which is 2.17×104-1, to decreasing the loop count in the early and late period of the simulation. This adjustment also causes that δ become further to the side of 1.3×104-1 in the middle period of the simulation. However, the loop count in this period could still be very low since current also reaches its peak and the plasma flow properties could remain constant approximately as mentioned before. For the simulation with “Initial-S,” the total runtime cannot be decreased anymore by adjusting δ to any side. The region where the loop count is in the relatively high level is near the exit of the capillary as shown in Figure 6(a). In this region, ∂E/∂T decreases and increases along the time direction and it is not consistently high or low as shown in Figure 7. It is possible that the average of ∂E/∂T near the exit of the capillary might even happen to be 0.7×10-4, so that there is not such a side that δ could be adjusted to for a less total runtime in the simulation with “Initial-S.” In order to reduce the runtime further, a possible approach might be dividing the whole spatial-temporal regions into 3-4 small regions and setting δ for these regions individually, so that δ could be closer to the average ∂E/∂T-1 of the corresponding spatial-temporal region. This might help to avoid the situation that the unified δ is far away from the true ∂E/∂T-1 in those regions such as the exit of the capillary in the start of the simulation.


## 7. Conclusion

A method based on the work of Trayner and Głowacki and Zaghloul discharge is introduced to solve the Saha equation in the simulation of capillary discharge used in the ET or ETC launch system. By modifying the specific energy equation, the multidimensional problem of solving the Saha equation is reduced into a one-dimensional iteration equation about the temperature. A coefficientδ is introduced into the iteration equation to satisfy the convergence condition of iteration method. In order to improve the computational efficiency, two methods including optimizing δ and setting the iteration initial value dynamically are introduced. Several numerical simulations of capillary discharge used in the ET or ETC launch system are conducted to test the detailed performance of the methods presented in this paper. The loop count of procedures (i)–(vi) and simulation run time is recorded to study the efficiency of the methods. The results show that the simulation runtime could be significantly reduced with an optimized δ and “Initial-T.” Firstly, “Initial-T” enables that the iteration initial value could be closer to be the true root in a cell; Secondly, “Initial-T” could make use of the pulse profile of the inputting current which enables that δ could be optimized further for the minimum run time compared with “Initial-S.”



---

*Source: 6473234-2017-11-20.xml*