Advances in Mathematical Physics

Volume 2017, Article ID 6473234, 9 pages

https://doi.org/10.1155/2017/6473234

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

^{1}School of Mechatronical Engineering, Beijing Institute of Technology, Beijing 100081, China^{2}State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, Beijing 100081, China

Correspondence should be addressed to Cheng Wu; nc.ude.tib@uwgnehc

Received 14 June 2017; Revised 17 August 2017; Accepted 22 October 2017; Published 20 November 2017

Academic Editor: Prabir Daripa

Copyright © 2017 Hengzhu Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.