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.
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]where is the ratio of the -fold ionized ions of the elemental species to the total heavy particles [12]. Apparently, determining the composition of plasma by solving the Saha equation actually means determining by solving (1). The function of is given by [12, 15, 18]where is the total number density of the heavy particles, is the ratio of the electron to the total heavy particles, denotes the partition functions which are dependent on the temperature , is the mass of the electron, is Boltzmann’s constant, is Planck’s constant, is the ionization potential. For the value of and necessary spectroscopic data used to calculate , one can find them from the database online [19]. is the lowering of the ionization potential which is introduced since the nonideal effect of plasma. is given by [4, 12]where is the permittivity of free space, is the electronic charge, is the de Broglie wave length, and is the standard Debye length which could be expressed aswhere is defined aswhere stands for the maximum allowed ionization stage of the elemental species . 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 bywhere (6) is the condition of electroneutrality of plasma. is the ratio of the elemental species to the total heavy particles. 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 or the mass density , and 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 , (1) could be expressed asNamely, could be expressed asFrom (7), can be expressed asBy substituting (9) into (10) with some rearrangements, one getsFirstly, for an easier case that , which is the ideal plasma, only in is unknown as one can find from (2), so that (11) could be expressed asEquation (12) indicates that could be expressed explicitly in terms of . By substituting (12) into (9), one getswhere . Similarly, as (11), only in the right side of (13) is unknown which indicates that can also be expressed explicitly in terms of . Then, together with (12), all of the element species could be expressed explicitly in terms of . Namely,Since (14) is actually effective for any element species j considered in a simulation case, (14) also means that of each element species considered in the simulation could be expressed explicitly in terms of . Then, by substituting (14) into (6), one could get an iteration equation about . Namely,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 to be the iterative initial value.(ii)Calculate (11) to obtain .(iii)Calculate (9) to obtain .(iv)Repeat steps (ii)-(iii) for all element species considered in the simulation case.(v)Calculate (6) to obtain a new denoted as .(vi)Set to be and repeat steps (ii)–(v) till the termination criteria are satisfied.
As one can find from procedures (i)–(vi), could be determined simultaneously with . Namely, for the case of , a relatively difficult multidimensional problem is reduced into a one-dimensional problem. For the present case that , as can be seen from (2)–(5), cannot be expressed explicitly in terms of unless is 1, so that both and are unknown in . In this case, (12) and (14) becomeThen, (15) becomesApparently, (17) cannot be solved with procedures (i)–(vi) directly. However, if is artificially set to be an empirical value, and might still be obtained by solving (17) with procedures (i)–(vi) as solving (15). Then, a new could be obtained by calculating (5). This is actually an iteration process about which could be expressed asThe bar over refers to the numerical nature of the function [12]. The procedure of solving (18) is as follows:(a)Set to be the iterative initial value.(b)Solve (17) with procedures (i)–(vi) to obtain .(c)Calculate (5) to obtain a new denoted as .(d)Set to be 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 which is easy to find its solution to any desired accuracy [12].
4. Solving the Saha Equation in the Simulation of Capillary Discharge
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, stands for time and stands for axial coordinate, is the mass density, is the velocity of the plasma flow, is the ablation rate which causes the increasing of , is the pressure, is the specific internal energy (SIE), is the conductivity, and 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 , so that one needs to solve the Saha equation in each computational cell to determine firstly, and, then, (19)–(23) could be solved for the plasma flow properties of next time step. However, (19)–(21) cannot give directly whereas is required as one of the inputting data pieces for solving the Saha equation. Similarly, as constructing (17) for the new emerged 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 while maintaining the one dimensionality of the method. The iteration equation about could be easily obtained by modifying (22) intowhere and represent the SIE and density determined by (19)–(21). These two could be viewed as inputting data for solving (24). stands for the result by calculating (22). The procedure of solving (24) is as follows:(A)Set to be the iterative initial value.(B)Solve (18) with procedures (a)–(d) to obtain .(C)Calculate (22) to obtain .(D)Calculate (24) to obtain a new denoted as .(E)Set to be 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 about . 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]where represents the right side of (24). From inequation (25), one can obtain the range of , which isIn the capillary discharge for the ET or ETC launch technology, the magnitude order of is usually (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 , where represents the true root in a cell, it might need least iterations to find the root for this cell. Namely,An ideal case is that if 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 in a cell, it means that, with getting closer to , would be more and more closer to 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 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 over all the cells at every time point. Since one cannot obtain 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 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 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, since 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 , , and , are dynamic spatial-temporal data, so that one may also need to set the iterative initial value of , , and dynamically. Otherwise, with the changing of inputting current, , , or 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 Figure 1, for solving the Saha equation or (24) in “present cell,” one could set the iterative initial value of , , or 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 , , and is set to be their average over all the cells at every time point from a previous simulated result.
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 [20]. For , the terminal criteria could be expressed as
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.
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].
(a)
(b)
Based on the above simulation, the efficiency of present methods is studied. In Figure 5, 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 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,” and , where could be , , , or other properties, actually mean the distance of the iterative initial value to the true root, so that these two methods connect and with the computational efficiency. It is probably that is generally larger than . From the simulation results, it is found that the magnitude order of is and the magnitude order of is , so that the magnitude order of is approximately 1 K. For the magnitude order of , it is varied from 0 at the closed end of the capillary to at the exit. Since is , the could reach approximately 180 K at most.
In Figure 6, 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 always increases along the axial direction in the capillary discharge, 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 in which 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, 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, 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.
(a)
(b)
A similarity in Figures 6(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 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 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 is plotted. It is found that could reach its peak of approximately at the exit of the capillary in the start of the simulation and reaches its lowest level of approximately 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 over all the cells at every time point is which is closer to the side of , and, then, the side of is the furthest point to the average of . For the reason that 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.
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 for “Initial-T” whereas is for “Initial-S.” The former one is not as close as that of “Initial-S” to the average of which is 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, is also in its relatively high level as shown in Figure 7, so that could be adjusted a little to the lower side of , which is , 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 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, 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 near the exit of the capillary might even happen to be , 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 of the corresponding spatial-temporal region. This might help to avoid the situation that the unified is far away from the true 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.”
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.