Numerical Calculation and Uncertain Optimization of Energy Conversion in Interior Ballistics Stage

: Gun ﬁring is a process that converts propellant chemical energy to projectile kinetic energy and other kinds of energies. In order to explore the energy conversion process, ﬁrstly, the interior ballistics mathematical model and the barrel-projectile ﬁnite element model are built and solved. Then, the related variable values and energy values are obtained and discussed. Finally, for improving energy e ﬃ ciency, the interval uncertainty optimization problem is modeled, and then solved using the two-layer nested optimization strategy and back-propagation (BP) neural network surrogate model. Calculation results show that, after optimization, the heat e ﬃ ciency raises from 31.13% to 33.05% and the max riﬂing stress decreases from 893.68 to 859.76 Mpa, which would improve the ﬁring performance and prolong the lifetime of the gun barrel.


Introduction
A gun is a kind of mechanical device which is used to fire projectiles through the gun firing system. The firing process of a gun can be briefly described as follows: firstly, the propellant combusts with generation of propellant gas, and then the propellant gas propels the projectile until it moves out of the gun barrel. This integral process is also known as the interior ballistics process. Besides movement of the projectile and combustion of the propellant, another three physicochemical phenomena also happen in the interior ballistics stage, namely friction and plastic deformation of the rotating band, recoil of the gun barrel and heat transfer of the gun barrel.
Many researchers have studied these physicochemical phenomena by numerical and experimental approaches. For the research on propellant combustion and the interior ballistics model, Liao et al. [1] established an innovative interior ballistics model which adopts actual combustion characteristics between propellants and pressure impulse. Monreal-González et al. [2,3] developed a program named UXGun to solve the one-dimensional interior ballistics model. They also compared numerical results with test date and the results showed a good agreement. An interior ballistics model based on conserved variables was built by  and calculation results were discussed by comparing with the model of Gough. A two-phase flow model of interior ballistics which contains igniters was given and four experiments were carried out by Jaramaz et al. [5] for investigating the interior ballistics process. To investigate detonation of granulated solid propellants, López et al. [6] used Rusanov and MacCormack-TVD (short for total variation diminishingschemes) numerical schemes to solve combustion equations and compared with calculation results by these two schemes. For the research on the movement of a projectile inside the gun barrel, Wu et al. [7][8][9][10][11] performed a series of engraving tests by employing a gas gun and various kinds of rotating bands. He discussed the effect of 2 of 21 shape, diameter, material and processing technology of rotating bands on engraving characteristics. Alexander et al. [12] described a finite element model of the firing system of a 155 mm gun. Solving this model by Abaqus software, he analyzed velocity, acceleration and displacement of the projectile. Sudarsan et al. [13] built a dynamics model of a 155 mm gun which contains a projectile gun barrel and rotating band, and then he conducted sensitivity analysis of the dependent parameters of shot start pressure. The influence of gun-projectile interactions as well as projectile-and gun-dependent factors on the dispersion were investigated by Dursun et al. [14]. Carlucci et al. [15] employed a simplified model to study the dynamic response of the projectile. His results indicated that accelerometer data should be used with caution. Ding et al. [16] proposed a mesh generation method for building a finite element model of the gun barrel. Moreover, he explored the relationship between a worn gun barrel and projectile velocity. For the research on heat transfer of the gun barrel, Zhang et al. [17] employed a method which combines the finite difference method with the sequential function specification method. His method can predict the heat flux of the gun barrel and has high accuracy by making comparisons of numerical tests. Considering material property to vary with temperature, Mathematica software was used by Hill et al. [18] for solving the one-dimensional transient heat transfer numerical model. In order to study the self-ignition phenomenon, Işık et al. [19] and Degirmenci et al. [20] used the finite element model of a 7.62 mm rifle to get temperature distribution of the gun barrel and discussed the effect of time interval, grain sizes and initial temperatures parameters. Talaee et al. [21] discussed the effect of different shooting rates in consecutive shooting on temperature distribution of the gun barrel. For the research on recoil of the gun barrel, a new mechanism is proposed by Sefiane et al. [22] to explain the occurrence of boiling crisis and it could improve recoil instability. Kathe and other researchers [23,24] suggested the concept of the soft recoil mechanism accelerating first in the forward direction in the armament technology. Mitchell et al. [25] studied the dynamics of the recoil system when firing projectiles with Mach 4.4 muzzle velocity. Hassaan et al. [26] investigated the dynamics of the barrel assembly-recoil mechanism of military cannons when using a hydraulic damper and a constant stiffness helical spring in their recoil mechanisms and evaluated the performance through the minimum and maximum displacements and the settling time of its response upon firing.
Heat machine is defined as a mechanical device which converts chemical energy to kinetic energy. From the energy point of view, guns also belong to the category of heat machine. Heat efficiency refers to the ratio of output kinetic energy to chemical energy of the heat machine. In engineering, heat efficiency is widely studied, such as heat efficiency of the gas turbine system [27], gasoline and ethanol mixture engine [28] and diesel engine [29]. Moreover, many methods have been employed to improve heat efficiency, such as optimizing structure and improving reliability [30,31]. For the gun, heat efficiency is also an important performance indicator, because higher heat efficiency means longer range and more penetration capability of the projectile. However, many researchers have built a simplified model and studied one or two physicochemical phenomena for getting several variable values, such as stress and velocity. Little attention has been paid to energy conversion in the interior ballistics stage.
In this paper, an energy calculation model of a fixed gun firing system is established and solved. Subsequently, energy values are obtained, and heat efficiency is calculated. Finally, the uncertain optimization is carried out for higher heat efficiency. This paper would provide an approach to investigate energy conversion of a fixed gun and would also be beneficial for investigating other kinds of thermal propulsion systems.

Interior Ballistics Mathematical Model
Before the propellant starts to combust, it is placed inside the gun barrel from breech to projectile bottom, as shown in Figure 1. When the propellant combusts, its products are gas-phase propellant and solid-phase propellant particles.
Energies 2020, 13 x FOR PEER REVIEW 3 of 22

Interior Ballistics Mathematical Model
Before the propellant starts to combust, it is placed inside the gun barrel from breech to projectile bottom, as shown in Figure 1. When the propellant combusts, its products are gas-phase propellant and solid-phase propellant particles. breech propellant gun barrel projectile The interior ballistics mathematical model describes the law of propellant combustion and product's state from the moment the propellant is ignited to the moment the projectile leaves the muzzle. Based on multiphase fluid mechanics theory and combustion theory, a modified onedimensional two-phase flow of the interior ballistics model is established. In addition, for establishing this model, some hypotheses are stated below: (1) The gas-phase propellant and solid-phase propellant have independent element volume, and they are considered as homogeneous continuum. (2) The appropriately defined average of flow properties is used to describe the heterogeneous compound flow, which consists of two interacting continua. (3) The solid-phase propellant is incompressible and its density value is fixed.
(4) The Nobel-Abel state equation is used to describe the gas-phase state. (5) The geometrical combustion theory and exponential combustion theory are used to describe the law of propellant combustion.
This model takes the characteristic that the cross-section area of the barrel bore varies with its distance from the breech into consideration. Firstly, five main equations are given.
Mass conservation equation of gas-phase: where φ is the porosity, g ρ is the gas-phase propellant density, g u is the gas-phase propellant velocity, A is the cross-section area of gun bore and  c m is the generation mass of propellant gas per unit time.
Momentum conservation equation of gas-phase: where p is the propellant gas pressure, p u is the solid-phase propellant velocity and s f is the interphase resistance. The interior ballistics mathematical model describes the law of propellant combustion and product's state from the moment the propellant is ignited to the moment the projectile leaves the muzzle. Based on multiphase fluid mechanics theory and combustion theory, a modified one-dimensional two-phase flow of the interior ballistics model is established. In addition, for establishing this model, some hypotheses are stated below: (1) The gas-phase propellant and solid-phase propellant have independent element volume, and they are considered as homogeneous continuum. (2) The appropriately defined average of flow properties is used to describe the heterogeneous compound flow, which consists of two interacting continua. (3) The solid-phase propellant is incompressible and its density value is fixed.
(4) The Nobel-Abel state equation is used to describe the gas-phase state. (5) The geometrical combustion theory and exponential combustion theory are used to describe the law of propellant combustion.
This model takes the characteristic that the cross-section area of the barrel bore varies with its distance from the breech into consideration. Firstly, five main equations are given.
Mass conservation equation of gas-phase: where ϕ is the porosity, ρ g is the gas-phase propellant density, u g is the gas-phase propellant velocity, A is the cross-section area of gun bore and . m c is the generation mass of propellant gas per unit time. Momentum conservation equation of gas-phase: where p is the propellant gas pressure, u p is the solid-phase propellant velocity and f s is the interphase resistance. Energy conservation equation of gas-phase: where Q p is the interphase heat transfer, e g is the internal energy of gas-phase propellant, e p is the potential energy of solid-phase propellant and ρ p is the solid-phase propellant density. Mass conservation equation of solid-phase: Momentum conservation equation of solid-phase: where R p is the particle stress.
To solve these equations, it is necessary to add another seven auxiliary equations, which are given in the following.
Gas-phase state equation: where R is the gas constant, T is the propellant gas temperature and β is the co-volume. Particle stress equation: where c 1 is the sound velocity of the solid-phase particle when ϕ = ϕ 0 , k is the stress attenuation factor, ϕ 0 is the initial porosity constant by natural packing of the propellant and its value is 0.45. Additionally, Interphase resistance equation: where d p is the equivalent diameter of the propellant particle.

of 21
Interphase heat transfer equation: where S p is the superficial area of the propellant particle, M p is the mass of the propellant particle and q is the heat flux coefficient between two phases. Temperature equation of particle surface: where T ps is the temperature of the particle surface, a p is the temperature conductivity coefficient and λ is the thermal conductivity coefficient. Here, the temperature principle is adopted as combustion physical condition. Specifically, the initial temperature of the particle is 293 K. If the temperature of the particle surface is higher than 615 K, which is the combustion temperature, then the propellant begins to combust. Combustion rate equation: where r p is the equivalent radius of the propellant particle, . d is the linear combustion rate in porous bed, b is the combustion rate coefficient and n is the combustion rate exponent.
Projectile motion equation: where ν is the projectile velocity, A d is the projectile bottom area, p d is propellant gas pressure at projectile bottom, f is projectile motion resistance and m is projectile mass. In order to solve the equations above, the MacCormack difference scheme which has second order accuracy is adopted, as shown in Equations (13a)-(13c). Fortran language is used to develop a solving program called TFIBS (short for two-phase flow interior ballistics solver).

Barrel-Projectile Finite Element Model
The barrel-projectile finite element model contains four parts: gun barrel, rifling, projectile and rotating band. This model, which is established by Hypermesh software, is shown in Figure 2, with partially hidden gun barrel and riflings for better visualization. Related material parameters are listed in Table 2. Totally, there are 1,239,092 elements and 1,421,076 notes in the model. In addition, element sizes of the gun barrel, rifling and projectile are approximately 5 mm. As for the rotating band, considering the accuracy of the engraving process, its meshing size is approximately 0.5 mm. The contact algorithm between rotating band and bore surface is set as general contact in Abaqus software. The dynamics and partial energy values are calculated by Abaqus solver.
Energies 2020, 13, 5824 6 of 21 partially hidden gun barrel and riflings for better visualization. Related material parameters are listed in Table 2. Totally, there are 1,239,092 elements and 1,421,076 notes in the model. In addition, element sizes of the gun barrel, rifling and projectile are approximately 5 mm. As for the rotating band, considering the accuracy of the engraving process, its meshing size is approximately 0.5 mm. The contact algorithm between rotating band and bore surface is set as general contact in Abaqus software. The dynamics and partial energy values are calculated by Abaqus solver. gun barrel rotating band projectile rifling   Figure 3 depicts the calculation process of energies. Firstly, we solve the interior ballistic model by TFIBS program and get the related variable values of propellant gas. Then, we divide the calculation process into two routes. One is to integrate the values of pressure and porosity to calculate energy Q1 and E4. Another one is adopting Abaqus subroutine Vuamp and Film to set the values of pressure and temperature as initial load and boundary condition in the finite element model, and then to use Abaqus to calculate energy E1, E2, E3 and E5.   Figure 3 depicts the calculation process of energies. Firstly, we solve the interior ballistic model by TFIBS program and get the related variable values of propellant gas. Then, we divide the calculation process into two routes. One is to integrate the values of pressure and porosity to calculate energy Q 1 and E 4 . Another one is adopting Abaqus subroutine Vuamp and Film to set the values of pressure and temperature as initial load and boundary condition in the finite element model, and then to use Abaqus to calculate energy E 1 , E 2 , E 3 and E 5 .

Calculation Results and Analysis of Energy Conversion
As mentioned above, several physicochemical phenomena happen in the interior ballistics stage, and appearance of these phenomena means occurrence of energy conversion. From the perspective of energy, the firing process can be described as: firstly, chemical energy of propellant converts into propellant gas energy through combustion; then, the propellant gas energy, as intermediate energy, converts into projectile kinetic energy, thermal energy of gun barrel and recoil energy, etc. According to the law of energy conservation, these energies can be given as: where Q is the chemical energy of total propellant, Q1 is the chemical energy of the uncombusted propellant, E1 is the translational kinetic energy of the projectile, E2 is the rotational kinetic energy of the projectile, E3 is the sum of frictional energy and plastic deformation energy, E4 is the recoil energy, E5 is the thermal energy of the gun barrel and E6 is the propellant gas energy. Among them, Q is a constant while Q1, E1, E2, E3, E4, E5 and E6 are time-varying. The energy conversion process is displayed in Figure 4 and will be explained in this section.

Calculation Results and Analysis of Energy Conversion
As mentioned above, several physicochemical phenomena happen in the interior ballistics stage, and appearance of these phenomena means occurrence of energy conversion. From the perspective of energy, the firing process can be described as: firstly, chemical energy of propellant converts into propellant gas energy through combustion; then, the propellant gas energy, as intermediate energy, converts into projectile kinetic energy, thermal energy of gun barrel and recoil energy, etc. According to the law of energy conservation, these energies can be given as: Energies 2020, 13, 5824 where Q is the chemical energy of total propellant, Q 1 is the chemical energy of the uncombusted propellant, E 1 is the translational kinetic energy of the projectile, E 2 is the rotational kinetic energy of the projectile, E 3 is the sum of frictional energy and plastic deformation energy, E 4 is the recoil energy, E 5 is the thermal energy of the gun barrel and E 6 is the propellant gas energy. Among them, Q is a constant while Q 1 , E 1 , E 2 , E 3 , E 4 , E 5 and E 6 are time-varying. The energy conversion process is displayed in Figure 4 and will be explained in this section. and appearance of these phenomena means occurrence of energy conversion. From the perspective of energy, the firing process can be described as: firstly, chemical energy of propellant converts into propellant gas energy through combustion; then, the propellant gas energy, as intermediate energy, converts into projectile kinetic energy, thermal energy of gun barrel and recoil energy, etc. According to the law of energy conservation, these energies can be given as: where Q is the chemical energy of total propellant, Q1 is the chemical energy of the uncombusted propellant, E1 is the translational kinetic energy of the projectile, E2 is the rotational kinetic energy of the projectile, E3 is the sum of frictional energy and plastic deformation energy, E4 is the recoil energy, E5 is the thermal energy of the gun barrel and E6 is the propellant gas energy. Among them, Q is a constant while Q1, E1, E2, E3, E4, E5 and E6 are time-varying. The energy conversion process is displayed in Figure 4 and will be explained in this section.
propellant gas gun barrel rotating band projectile

Numerical Calculation of E1 and E2
As is exhibited in Figure 4, combustion of the propellant produces propellant gas, which is filled in the space from breech to projectile bottom and pushes the projectile to do translational and rotational movement inside the gun barrel. At the same time, the projectile gets translational kinetic energy and rotational kinetic energy. E1 and E2 can be expressed by:

Numerical Calculation of E 1 and E 2
As is exhibited in Figure 4, combustion of the propellant produces propellant gas, which is filled in the space from breech to projectile bottom and pushes the projectile to do translational and rotational movement inside the gun barrel. At the same time, the projectile gets translational kinetic energy and rotational kinetic energy. E 1 and E 2 can be expressed by: where m is the projectile mass, v is the translational velocity of the projectile, j is the rotational inertia of the projectile and w is the rotational velocity of the projectile. m and j are constants and their values are listed in Table 1. Solving the interior ballistics model by the TFIBS program, propellant gas pressure is obtained and its data at some time points are shown in Figure 5a. In addition, Figure 5b shows propellant gas pressures at the projectile bottom.
Energies 2020, 13 x FOR PEER REVIEW 8 of 22 where m is the projectile mass, v is the translational velocity of the projectile, j is the rotational inertia of the projectile and w is the rotational velocity of the projectile. m and j are constants and their values are listed in Table 1. Solving the interior ballistics model by the TFIBS program, propellant gas pressure is obtained and its data at some time points are shown in Figure 5a. In addition, Figure 5b shows propellant gas pressures at the projectile bottom. By loading the pressure shown in Figure 5b to projectile bottom in the finite element model with adoption of Abaqus subroutine Vuamp, translational velocity and rotational velocity of the projectile can be obtained by Abaqus, as shown in Figure 6a, b. Obviously, curves of the projectile translational velocity and rotational velocity are highly similar. This is because the entanglement angle of the rifling is fixed, so that the translational velocity is proportional to the rotational velocity. E1 and E2 are calculated by Equations (15) and (16). Their values are shown in Figure 6a, b, respectively. From Figure 6a, the muzzle translational velocity is 909 m•s −1 . Table 3 gives muzzle translational velocity data of six firing tests, which are obtained by high-speed photography. Comparing E1 obtained by By loading the pressure shown in Figure 5b to projectile bottom in the finite element model with adoption of Abaqus subroutine Vuamp, translational velocity and rotational velocity of the projectile can be obtained by Abaqus, as shown in Figure 6a,b. Obviously, curves of the projectile translational velocity and rotational velocity are highly similar. This is because the entanglement angle of the rifling is fixed, so that the translational velocity is proportional to the rotational velocity. E 1 and E 2 are calculated by Equations (15) and (16). Their values are shown in Figure 6a,b, respectively. From Figure 6a, the muzzle translational velocity is 909 m·s −1 . Table 3 gives muzzle translational velocity data of six firing tests, which are obtained by high-speed photography. Comparing E 1 obtained by numerical simulation with E 1 measured by firing tests, the error is −2.58%. By loading the pressure shown in Figure 5b to projectile bottom in the finite element model with adoption of Abaqus subroutine Vuamp, translational velocity and rotational velocity of the projectile can be obtained by Abaqus, as shown in Figure 6a, b. Obviously, curves of the projectile translational velocity and rotational velocity are highly similar. This is because the entanglement angle of the rifling is fixed, so that the translational velocity is proportional to the rotational velocity. E1 and E2 are calculated by Equations (15) and (16). Their values are shown in Figure 6a, b, respectively. From Figure 6a, the muzzle translational velocity is 909 m•s −1 . Table 3 gives muzzle translational velocity data of six firing tests, which are obtained by high-speed photography. Comparing E1 obtained by numerical simulation with E1 measured by firing tests, the error is −2.58%.    As shown in Figure 6a, duration from the moment the propellant starts combusting to the moment the projectile leaves the muzzle is 18.9 ms, and this process is the interior ballistics stage. According to projectile bottom pressure, the generation period of E 1 can be divided into three periods. The time between 0 and 3.6 ms belongs to the first period. In this period, the propellant combusts apart, and the propellant at the projectile bottom has not yet combusted. Hence, the projectile bottom pressure is 0 MPa and the projectile is static. Taking the pressure data at 3.2 ms in Figure 5a as an example, the total propellant length is 1.08 m. From 0 to 0.7 m, the propellant gas pressure is greater than 0 Mpa, and the propellant gas pressure from 0.7 to 1.08 m is equal to 0 Mpa. It means that the propellant only combusts to 0.7 m and the propellant from 0.7 to 1.08 m has not yet started to combust. The time of 3.6-4.7 ms belongs to the second period. During this period, the projectile bottom pressure gradually increases from 0 to 25.6 Mpa. Meanwhile, the static motion resistance of the projectile also increases and is equal to it. So, the projectile is still stationary. The third period is from 4.7 to 18.9 ms. During this period, as the projectile bottom pressure continues to increase, the projectile begins to move, and the translational velocity, rotational velocity, translational kinetic energy and rotational energy of the projectile increase with time.

Numerical Calculation of E 3
The rotating band, which is ring-shaped and surrounds the projectile, moves along with the projectile inside the gun barrel. By partially hiding the gun barrel and riflings for better visualization, Figure 7 presents the size relationship between rotating band diameter (D 1 ) and gun bore diameter Energies 2020, 13, 5824 9 of 21 (D 2 ) before the projectile moves. Apparently, contact will occur between the rotating band, riflings and inner surface of the gun barrel after the projectile moves. The purpose of this structural design is to make the projectile rotate and prevent leakage of propellant gas. In addition, the type of projectile in this paper has two rotating bands.
projectile begins to move, and the translational velocity, rotational velocity, translational kinetic energy and rotational energy of the projectile increase with time.

Numerical Calculation of E3
The rotating band, which is ring-shaped and surrounds the projectile, moves along with the projectile inside the gun barrel. By partially hiding the gun barrel and riflings for better visualization, Figure 7 presents the size relationship between rotating band diameter (D1) and gun bore diameter (D2) before the projectile moves. Apparently, contact will occur between the rotating band, riflings and inner surface of the gun barrel after the projectile moves. The purpose of this structural design is to make the projectile rotate and prevent leakage of propellant gas. In addition, the type of projectile in this paper has two rotating bands.  Again, the finite element model and the projectile bottom pressure in Section 3.1 are adopted for calculation, and Abaqus is used to output the Mises stress and morphology of the rotating band and riflings at 6.2 and 9.5 ms, as shown in Figures 8 and 9. According to the characteristics of energy, the generation period of E3 can be divided into two distinct parts. The engraving period is from 4.7 to 8.2 ms. During this period, the rotating band is squeezing into the riflings and is subjected to the cutting and squeeze of the riflings. When the Mises stress experienced by the rotating band is greater than the yield limit of its material, the rotating band will undergo plastic deformation and form rotating band grooves, which are shown in Figure 8. Meanwhile, the rotating band also produces friction with riflings and the inner surface of the gun barrel. The translation-rotation period is from 8.2 to 18.9 ms. During this period, rotating band grooves have been completely formed, and the rotating band no longer undergoes plastic deformation. The shape of the rotating band is shown in Figure 9. However, the rotating band still produces friction with riflings and the inner surface of the gun barrel. In conclusion, both frictional energy and plastic deformation energy generate in the engraving period, while only frictional energy generates in the translation-rotation period. Again, the finite element model and the projectile bottom pressure in Section 3.1 are adopted for calculation, and Abaqus is used to output the Mises stress and morphology of the rotating band and riflings at 6.2 and 9.5 ms, as shown in Figures 8 and 9. According to the characteristics of energy, the generation period of E 3 can be divided into two distinct parts. The engraving period is from 4.7 to 8.2 ms. During this period, the rotating band is squeezing into the riflings and is subjected to the cutting and squeeze of the riflings. When the Mises stress experienced by the rotating band is greater than the yield limit of its material, the rotating band will undergo plastic deformation and form rotating band grooves, which are shown in Figure 8. Meanwhile, the rotating band also produces friction with riflings and the inner surface of the gun barrel. The translation-rotation period is from 8.2 to 18.9 ms. During this period, rotating band grooves have been completely formed, and the rotating band no longer undergoes plastic deformation. The shape of the rotating band is shown in Figure 9. However, the rotating band still produces friction with riflings and the inner surface of the gun barrel. In conclusion, both frictional energy and plastic deformation energy generate in the engraving period, while only frictional energy generates in the translation-rotation period. The frictional energy and the plastic deformation energy are calculated by Abaqus energy output functions ALLFD (short for all frictional dissipated energy) and ALLPD (short for all plastic dissipated energy), respectively. The value of E3 is shown in Figure 10. From 8.2 to about 15 ms, the increase rate of E3 is less than its increase rate during the engraving period, which is due to the lack of plastic deformation energy during this period. From about 15 to 18.9 ms, the average increase rate of E3 is greater than its average increase rate during the engraving period, because after 15 ms, the The frictional energy and the plastic deformation energy are calculated by Abaqus energy output functions ALLFD (short for all frictional dissipated energy) and ALLPD (short for all plastic dissipated energy), respectively. The value of E3 is shown in Figure 10. From 8.2 to about 15 ms, the increase rate of E3 is less than its increase rate during the engraving period, which is due to the lack of plastic deformation energy during this period. From about 15 to 18.9 ms, the average increase rate of E3 is greater than its average increase rate during the engraving period, because after 15 ms, the The frictional energy and the plastic deformation energy are calculated by Abaqus energy output functions ALLFD (short for all frictional dissipated energy) and ALLPD (short for all plastic dissipated energy), respectively. The value of E 3 is shown in Figure 10. From 8.2 to about 15 ms, the increase rate of E 3 is less than its increase rate during the engraving period, which is due to the lack of plastic deformation energy during this period. From about 15 to 18.9 ms, the average increase rate of E 3 is greater than its average increase rate during the engraving period, because after 15 ms, the value of projectile velocity is increasing and relatively high, which makes the projectile displacement higher in unit time.
rotating band groove Figure 9. Mises stress and morphology at 9.5 ms in the translation-rotation period.
The frictional energy and the plastic deformation energy are calculated by Abaqus energy output functions ALLFD (short for all frictional dissipated energy) and ALLPD (short for all plastic dissipated energy), respectively. The value of E3 is shown in Figure 10. From 8.2 to about 15 ms, the increase rate of E3 is less than its increase rate during the engraving period, which is due to the lack of plastic deformation energy during this period. From about 15 to 18.9 ms, the average increase rate of E3 is greater than its average increase rate during the engraving period, because after 15 ms, the value of projectile velocity is increasing and relatively high, which makes the projectile displacement higher in unit time.
where k p is the coefficient of heat value, m p is the propellant mass and ρ p is the propellant density, and the values of these three constants are listed in Table 1. The calculated Q is 6.05 × 10 7 J. ϕ is the porosity and is obtained by the TFIBS program. The data of porosity at some time points are exhibited in Figure 11a. m 1 is the uncombusted propellant mass. b is the percentage of uncombusted propellant energy, and its value is shown in Figure 11b. As can be seen from Figure 11a, the porosity tends to increase roughly before 16.4 ms. This indicates that the propellant is combusting and producing propellant gas. At 16.4 ms, the porosity value of all spaces inside the gun barrel is 1, which means that the propellant has completely combusted and there is no solid-phase propellant. In other words, from 16.4 to 18.9 ms, no new propellant gas is produced, and the movement of the projectile depends on the propelling effect of the existing propellant gas. This propelling force is greater than the frictional resistance experienced by the projectile, so the projectile kinetic energy is still increasing. The absolute value of the slope of the curve in Figure 11b represents the combustion rate of the propellant, and the combustion rate generally shows a trend of increasing first and then decreasing. Referring to Equation (11b), this is because the combustion rate of the propellant is a monotonous increasing function of the propellant gas pressure, and the propellant gas pressure first increases and then decreases in the interior ballistics stage.
where p k is the coefficient of heat value, p m is the propellant mass and p ρ is the propellant density, and the values of these three constants are listed in Table 1. The calculated Q is 6.05 × 10 7 J. φ is the porosity and is obtained by the TFIBS program. The data of porosity at some time points are exhibited in Figure 11a. As can be seen from Figure 11a, the porosity tends to increase roughly before 16.4 ms. This indicates that the propellant is combusting and producing propellant gas. At 16.4 ms, the porosity value of all spaces inside the gun barrel is 1, which means that the propellant has completely combusted and there is no solid-phase propellant. In other words, from 16.4 to 18.9 ms, no new propellant gas is produced, and the movement of the projectile depends on the propelling effect of the existing propellant gas. This propelling force is greater than the frictional resistance experienced by the projectile, so the projectile kinetic energy is still increasing. The absolute value of the slope of the curve in Figure 11b represents the combustion rate of the propellant, and the combustion rate generally shows a trend of increasing first and then decreasing. Referring to Equation (11b), this is because the combustion rate of the propellant is a monotonous increasing function of the propellant gas pressure, and the propellant gas pressure first increases and then decreases in the interior ballistics stage.

Numerical Calculation of E4 and E5
As is exhibited in Figure 4, the propellant gas also exerts force on the breech, which propels the gun barrel to move with a direction opposite to the projectile movement. b p is the propellant gas pressure at the breech and is obtained by the TFIBS program, as shown in Figure 12a. Recoil force, z F , generates from the gun recoil mechanism and can be expressed by Equation (21). By loading b p and z F on the breech of the finite element model in Section 3.1, L is the gun barrel displacement and is obtained by Abaqus, as shown in Figure 12a. Finally, 4 E can be calculated by Equation (22) and its value is shown in Figure 12b.

Numerical Calculation of E 4 and E 5
As is exhibited in Figure 4, the propellant gas also exerts force on the breech, which propels the gun barrel to move with a direction opposite to the projectile movement. p b is the propellant gas pressure at the breech and is obtained by the TFIBS program, as shown in Figure 12a. Recoil force, F z , generates from the gun recoil mechanism and can be expressed by Equation (21). By loading p b and F z on the breech of the finite element model in Section 3.1, L is the gun barrel displacement and is obtained by Abaqus, as shown in Figure 12a. Finally, E 4 can be calculated by Equation (22) and its value is shown in Figure 12b.
where K 1 is the hydraulic resistance coefficient, Ω 1 is the minimum cross-sectional area of the hydraulic branch, α x is the area of fluid orifice, A f j is the cross-sectional area of piston, v l is the recoil velocity and A b is the breech area.
Energies 2020, 13 x FOR PEER REVIEW 12 of 22 where 1 K is the hydraulic resistance coefficient, 1 Ω is the minimum cross-sectional area of the hydraulic branch, x α is the area of fluid orifice, fj A is the cross-sectional area of piston, l v is the recoil velocity and b A is the breech area. As is exhibited in Figure 4, due to the temperature difference between the propellant gas and the gun barrel, heat transfer occurs at the contact surface between them. This process causes the gun barrel temperature to increase and the gun barrel to obtain thermal energy. The average temperature of propellant gas calculated by the TFIBS program is shown in Figure 13a. By using Abaqus subroutine Film to apply the average temperature value of propellant gas as the boundary condition As is exhibited in Figure 4, due to the temperature difference between the propellant gas and the gun barrel, heat transfer occurs at the contact surface between them. This process causes the gun barrel temperature to increase and the gun barrel to obtain thermal energy. The average temperature of propellant gas calculated by the TFIBS program is shown in Figure 13a. By using Abaqus subroutine Film to apply the average temperature value of propellant gas as the boundary condition to the inner surface of the gun barrel, the temperature distribution is calculated by Abaqus, as shown in Figure 14. Furthermore, E 5 is calculated by Abaqus energy output function ALLIE (short for all internal energy), and the result is shown in Figure 13b. As is exhibited in Figure 4, due to the temperature difference between the propellant gas and the gun barrel, heat transfer occurs at the contact surface between them. This process causes the gun barrel temperature to increase and the gun barrel to obtain thermal energy. The average temperature of propellant gas calculated by the TFIBS program is shown in Figure 13a. By using Abaqus subroutine Film to apply the average temperature value of propellant gas as the boundary condition to the inner surface of the gun barrel, the temperature distribution is calculated by Abaqus, as shown in Figure 14. Furthermore, E5 is calculated by Abaqus energy output function ALLIE (short for all internal energy), and the result is shown in Figure 13b. At 18.9 ms, the ratio between E4 and Q is 0.56%, and the ratio between E4 and Q is 1.18%. The values of both E4 and E5 are relatively low. This is the reason that E4 and E5 are ignored sometimes. In engineering, for convenience, E4 and E5 can be calculated roughly in the way that Q is multiplied by the ratio. This ratio is determined by the type of gun and it can be found in the interior ballistics manual.

Numerical Calculation of E6
The propellant gas energy, which contains thermal energy of propellant gas and kinetic energy of propellant gas, is intermediate energy and could convert to other kinds of energies. After the projectile leaves the muzzle, the propellant gas will emit out of the gun barrel as exhaust gas. In Equation (14), the other kinds of energies have been obtained, except E6. So, E6 can be given by Equation (23) and its value is shown in Figure 15. At 18.9 ms, the ratio between E 4 and Q is 0.56%, and the ratio between E 4 and Q is 1.18%. The values of both E 4 and E 5 are relatively low. This is the reason that E 4 and E 5 are ignored sometimes. In engineering, for convenience, E 4 and E 5 can be calculated roughly in the way that Q is multiplied by the ratio. This ratio is determined by the type of gun and it can be found in the interior ballistics manual.

Numerical Calculation of E 6
The propellant gas energy, which contains thermal energy of propellant gas and kinetic energy of propellant gas, is intermediate energy and could convert to other kinds of energies. After the projectile leaves the muzzle, the propellant gas will emit out of the gun barrel as exhaust gas. In Equation (14), the other kinds of energies have been obtained, except E 6 . So, E 6 can be given by Equation (23) and its value is shown in Figure 15.
In the early stage of interior ballistics, since the combustion rate of the propellant is increasing, the values of E 1 , E 2 , E 3 , E 4 and E 5 are relatively low. So, E 6 increases with time. In the later stage of interior ballistics, conversely, because the combustion rate of the propellant continuously decreases until zero, the values of E 1 , E 2 , E 3 , E 4 and E 5 are relatively high. So, E 6 decreases with time.
Energies 2020, 13, 5824 13 of 21 of propellant gas, is intermediate energy and could convert to other kinds of energies. After the projectile leaves the muzzle, the propellant gas will emit out of the gun barrel as exhaust gas. In Equation (14), the other kinds of energies have been obtained, except E6. So, E6 can be given by Equation (23) and its value is shown in Figure 15. In the early stage of interior ballistics, since the combustion rate of the propellant is increasing, the values of E1, E2, E3, E4 and E5 are relatively low. So, E6 increases with time. In the later stage of interior ballistics, conversely, because the combustion rate of the propellant continuously decreases until zero, the values of E1, E2, E3, E4 and E5 are relatively high. So, E6 decreases with time.

Energy Distribution Proportion
According to application characteristics, E1 and E2 belong to the effective energy, which designers aim to get, E3, E4 and E5 belong to the secondary energy, which appears inevitably, and E6 belongs to the potential energy, which is able to convert to other kinds of energies. Based on the calculation results above, Table 4 summarizes energy values at 18.9 ms. Additionally, Figure 16 illustrates the proportionate relationship between each kind of energy and total propellant chemical energy. For instance, the percentage of effective energy, η 1 , can be given as:

Energy Distribution Proportion
According to application characteristics, E 1 and E 2 belong to the effective energy, which designers aim to get, E 3 , E 4 and E 5 belong to the secondary energy, which appears inevitably, and E 6 belongs to the potential energy, which is able to convert to other kinds of energies. Based on the calculation results above, Table 4 summarizes energy values at 18.9 ms. Additionally, Figure 16 illustrates the proportionate relationship between each kind of energy and total propellant chemical energy. For instance, the percentage of effective energy, η 1 , can be given as:   For this type of fixed gun, its heat efficiency, namely percentage of effective energy, is 31.13%. That is to say, the projectile gets 31.13% of Q as kinetic energy, while the rest (68.87%) of Q converts to other kinds of energies. More specifically, 4.09% of Q dissipates directly and it is a relatively low percentage. E6 accounts for 64.78% of Q and its percentage is relatively high. So, it seems that there is a possibility to decrease the percentage of E6 further. Generally, two methods can be employed to pursue higher heat efficiency. One is recycling the propellant exhaust gas, and the other one is optimizing the parameters of the gun firing system. Considering the difficulty of recycling propellant exhaust gas, an optimization work will be carried out and discussed in the next section.

Uncertain Optimization of Energy Conversion
A gun is a firing device, and as a result, the projectile kinetic energy is an important performance indicator. Moreover, propellant gas pressure and rifling stress should also be considered in design since they have a significant effect on the lifetime of the gun barrel. Due to the uncertainty in the interior ballistics stage, an uncertain optimization is performed to obtain higher heat efficiency in this section. This optimization takes projectile kinetic energy as the objective function and takes  For this type of fixed gun, its heat efficiency, namely percentage of effective energy, is 31.13%. That is to say, the projectile gets 31.13% of Q as kinetic energy, while the rest (68.87%) of Q converts to other kinds of energies. More specifically, 4.09% of Q dissipates directly and it is a relatively low percentage. E 6 accounts for 64.78% of Q and its percentage is relatively high. So, it seems that there is a possibility to decrease the percentage of E 6 further. Generally, two methods can be employed to pursue higher heat efficiency. One is recycling the propellant exhaust gas, and the other one is optimizing the parameters of the gun firing system. Considering the difficulty of recycling propellant exhaust gas, an optimization work will be carried out and discussed in the next section.

Uncertain Optimization of Energy Conversion
A gun is a firing device, and as a result, the projectile kinetic energy is an important performance indicator. Moreover, propellant gas pressure and rifling stress should also be considered in design since they have a significant effect on the lifetime of the gun barrel. Due to the uncertainty in the interior ballistics stage, an uncertain optimization is performed to obtain higher heat efficiency in this section. This optimization takes projectile kinetic energy as the objective function and takes maximum propellant gas pressure and rifling stress as constraint functions.

Algorithm Model of Uncertain Optimization
Without loss of generality, the interval method is employed to model the uncertainty, hence the objective function and the constraints with respect to the interval uncertain vector can be written in a form of nonlinear continuous function. By this, an interval optimization problem can be formulated as follows: where U is a n-dimensional interval uncertainty vector, and the superscripts I, L,R, c and w are the interval, the lower bounds, the upper bounds, the center and the radius of an interval number, respectively. f and g are the interval objective function and the interval constraints, respectively. b I i is the permissible interval of the i-th constraint, and it will degenerate into a real number when b I i = b R i . Obviously, the traditional optimization algorithms cannot directly solve Equation (25), and some special treatments are required.
The interval order relations can judge that an interval number is better than another, and thus they are usually used to treat the uncertainty objective function. The interval order relation ≤ cw , proposed by Ishibuchi et al. [32], is herein adopted to rank the interval numbers. Assuming that there are two interval numbers, A and B, if and only if the center and the radius of B are both smaller than that of A, then interval B is better than interval A, i.e., A ≤ cw B (Jiang et al. [33]). Using ≤ cw , we can transform the interval objective function in Equation (25) into a multi-objective optimization problem: The lower and upper bounds of the interval objective function, i.e., f L (U) and f R (U), can be computed through two deterministic optimization processes, which can be formulated as follows: The above equations transform an uncertain objective function into two deterministic ones. Similar to the stochastic programming, f c (U) reflects the average performance of the objective function with respect to the uncertainty, and f w (U) is the fluctuation range of the objective function caused by the uncertainties. Therefore, we can consider the robustness in this optimization model by minimizing f w (U).
Besides, the interval probability degree model in Reference [34] is adopted to complete the deterministic transformation of uncertain constraints. By this, the inequality constraint, g i (U) ≤ b I i , can be rewritten as a deterministic one: where λ i is the presupposed threshold value of interval possibility degree. The larger λ i means the more rigorous constraint. The equality constraint g I i (U) = b I i can be transformed as follows: It can be further transformed through two deterministic inequality constraints: The interval of the constraint, g I i (U), can also be computed through two deterministic optimization processes: Through Equations (26)-(31), the interval uncertain optimization problem as in Equation (25) can be rewritten as a multi-objective optimization model: Through the above modeling, the obtained Equation (32) is generally solved by the two-layer nested optimization strategy. The inner optimizer is used to compute the interval bounds of the objective function and constraints with respect to the uncertainty, and the outer optimizer is adopted to search for the optimal solution of the design variables. The multi-objective optimization algorithm, namely the non-dominated sorting genetic algorithm-II (NSGA-II), is employed as the outer optimizer since Equation (32) is a typical multi-objective optimization problem. Moreover, the inner optimization is completed by the genetic algorithm (GA), which has fine global search ability so it avoids the local optimal solution. However, it takes more computational costs than the gradient-based algorithms. For more details of this solving strategy, the reader is referred to Reference [34]. Finally, the solving process of the uncertain optimization is depicted in Figure 17.

Optimization Model and Result Analysis
In design work, the maximum gas pressure that the gun barrel can bear depends on the material of the gun barrel. Once the maximum gas pressure exceeds the yield strength of the gun barrel material, the plastic expansion of the gun barrel would happen. For a large-caliber gun, the lifetime of the gun barrel usually allows several hundred uses due to the wear of the riflings. So, the rifling stress should be in a proper range. By maximizing the interval center of projectile kinetic energy and minimizing its interval radius, and meanwhile considering the constraints of the gas pressure and rifling stress, the optimization model is formulated as follows: , where A E is the projectile kinetic energy and is equal to the sum of E1 and E2, bmax P is the max propellant gas pressure at the breech and rmax S is the max rifling stress at the initial position.
Moreover, in actual engineering practice, smaller interval radius of the uncertain design variables usually represents higher manufacturing precision and higher production cost. Therefore, the best design scheme is to achieve the required performances under the maximum parameter errors, which

Optimization Model and Result Analysis
In design work, the maximum gas pressure that the gun barrel can bear depends on the material of the gun barrel. Once the maximum gas pressure exceeds the yield strength of the gun barrel material, the plastic expansion of the gun barrel would happen. For a large-caliber gun, the lifetime of the gun barrel usually allows several hundred uses due to the wear of the riflings. So, the rifling stress should be in a proper range. By maximizing the interval center of projectile kinetic energy and minimizing its interval radius, and meanwhile considering the constraints of the gas pressure and rifling stress, the optimization model is formulated as follows: where E A is the projectile kinetic energy and is equal to the sum of E 1 and E 2 , P bmax is the max propellant gas pressure at the breech and S rmax is the max rifling stress at the initial position. Moreover, in actual engineering practice, smaller interval radius of the uncertain design variables usually represents higher manufacturing precision and higher production cost. Therefore, the best design scheme is to achieve the required performances under the maximum parameter errors, which means the low manufacturing costs. Thus, ζ in Equation (33) is a non-negative and dimensionless evaluation coefficient which is designed herein to evaluate the overall uncertainty of all design variables: where A X w i is the preset maximum possible error of the design variable, and n is the number of design variables. The bigger ζ has the higher overall error level and the lower cost. The initial values of design variables are listed in Table 5. In Table 5, x 1 is the external diameter of the rotating band, x 2 is the width of the rotating band, x 3 is the rifling width, x 4 is the rifling height, x 5 is the projectile mass, x 6 is the forcing cone taper, x 7 is the propellant thickness and x 8 is the combustion rate coefficient. Figure 18 depicts optimization results, i.e., pareto front. Here, E c A and E w A denote the center and the radius of the projectile kinetic energy interval, respectively. All solutions in pareto front are feasible and satisfy the constraints. The choice of a final solution depends on the preference of the decision maker, and the maximum projectile kinetic energy is herein taken as the decision preference. Subsequently, Table 6 lists the optimal combination and intervals of design variables. Unlike deterministic optimization, the specific tolerance range of uncertain parameters are obtained, which provides a theoretical reference for the error formulation of design variables. By substituting the data in Table 6 to the original model, Figure 19 shows the specific, optimal ranges of objective function and constraint functions. Detailed optimization results are listed in Table 7.
Energies 2020, 13 x FOR PEER REVIEW 18 of 22 means the low manufacturing costs. Thus, ζ in Equation (33) is a non-negative and dimensionless evaluation coefficient which is designed herein to evaluate the overall uncertainty of all design variables: where w i X A is the preset maximum possible error of the design variable, and n is the number of design variables. The bigger ζ has the higher overall error level and the lower cost. The initial values of design variables are listed in Table 5. In Table 5, 1 x is the external diameter of the rotating band, 2 x is the width of the rotating band, 3 x is the rifling width, 4 x is the rifling height, 5 x is the projectile mass, 6 x is the forcing cone taper, 7 x is the propellant thickness and 8 x is the combustion rate coefficient.  Table 6 lists the optimal combination and intervals of design variables. Unlike deterministic optimization, the specific tolerance range of uncertain parameters are obtained, which provides a theoretical reference for the error formulation of design variables. By substituting the data in Table 6 to the original model, Figure 19 shows the specific, optimal ranges of objective function and constraint functions. Detailed optimization results are listed in Table 7.     Figure 19 shows that the performance indicator from the uncertain optimization method is also an interval value. This fully reflects the fluctuation of these dynamic responses which is caused by uncertainty. Hence, the upper and lower bounds of these performance indicators can provide a good reference for designers.  Table 7, the projectile kinetic energy increases significantly, and the increment is 6.78% with constraint conditions satisfied, which would improve the heat efficiency of the gun firing system.   Figure 19 shows that the performance indicator from the uncertain optimization method is also an interval value. This fully reflects the fluctuation of these dynamic responses which is caused by uncertainty. Hence, the upper and lower bounds of these performance indicators can provide a good reference for designers. From Table 7, the projectile kinetic energy increases significantly, and the increment is 6.78% with constraint conditions satisfied, which would improve the heat efficiency of the gun firing system. Furthermore, the max rifling stress decreases from 893 to 859.76 Mpa, which conduces to prolong the lifetime of the gun barrel. Table 8 shows a comparison of energies before and after optimization. As can be seen, after optimization, the percentage of E 6 decreased from 64.78% to 61.91%, while the other kinds of energy increased. It means that the optimization improves utilization of E 6 and the reduced part of E 6 distributes to other kinds of energies. The ideal optimization result is that of increasing E 1 and E 2 , as well as decreasing the other energies. But, it is hard to realize this result. For example, the higher propellant gas caused the increase in not only E 1 and E 2 but also E 4 .

Conclusions
In this paper, an interior ballistics mathematical model and barrel-projectile finite element model were established and solved. Then, the values of related variables such as pressure and porosity were obtained, the generation process of energies was explained, and especially, the energy versus time data were obtained and discussed. Furthermore, the uncertainty optimization was carried out for higher heat efficiency. From the calculation results, the following conclusions can be drawn: (1) For this type of fixed gun, interior ballistics duration was 18.9 ms when using initial parameters. At 18.9 ms, percentage of effective energy was 31.13%, percentage of secondary energy was 4.09% and percentage of potential energy was 64.78%. It indicates that a while small proportion of propellant chemical energy dissipates directly in the interior ballistics stage, most of it converts to propellant gas energy. (2) After the uncertain optimization, the percentage of the projectile kinetic energy rose from 31.13% to 33.05%, and the percentage of the improvement was 6.78%. It means that heat efficiency and firing performance improved. In addition, the max rifling stress decreased by 3.79%, which is a benefit for prolonging the lifetime of the gun barrel.
The gun firing model and optimization model built in this paper could calculate and optimize not only energies but also other dynamical variables such as stress and strain. They would be a benefit for the improvement of the comprehensive performance of the gun firing system.

Conflicts of Interest:
The authors declare no conflict of interest.