Comparison of Computational Methods and Approaches Applied in Formulation of Boundary Conditions in Lagrange’s Ballistic Problem

This paper compares algorithms for constructing a differential solution to the gas dynamics equations in the surrounding of a moving boundary. In order to compare the algorithms, the Lagrange problem, also known as the piston problem, was chosen as a test problem. A comparison was made between the author’s algorithm based on the method of characteristics and the classical approach with a fictitious cell using the solution of the Riemann problem. Calculations were carried out for the case of subsonic and supersonic gas motion. Comparisons were made on a fixed grid with a dynamically expanding cell (Euler grid) and on a uniformly expanding grid (Arbitrary Lagrangian-Eulerian – ALE grid). The problem was solved using numerical schemes of second order in time and space Lax - Wendroff type: the Richtmyer scheme and the MacCormack scheme. The results of calculations using the characteristics method were used as a benchmark solution. The comparative analysis carried out allows conclusions to be drawn regarding the accuracy of the individual approaches, as well as the scope of their applicability.


Introduction
Historically, internal ballistics calculations were performed using zero-dimensional thermodynamic models.With the development of computational fluid dynamics (CFD), distributed-parameter internal ballistics models were developed to include a number of additional phenomena.Nowadays, while many ways of solving the gas dynamic flow equations are known, there is little information in the literature concerned about how to perform calculations near the boundary of the computational domain (especially in case of moving boundaries).
In the paper [1], numerical methods in computational fluid dynamics are very extensively presented.Nevertheless, the method of determining the boundary conditions is only presented in a cursory manner.
Reports on gas dynamics models of internal ballistics [2][3][4][5][6][7][8][9][10][11] neglect the issue of boundary conditions or only mention them, rarely providing a numerical formulation.In paper [4], the authors proposed an original approach with uniformly expanding cells combined with dynamically generated new nodes.In the report [2] and in the paper [3], the solution of the internal ballistics problem for a multiphase flow model uses the method of characteristics to determine the solution near the boundaries.
In this paper, the classical Lagrange problem is solved in order to compare solution construction algorithms near a moving boundary.Calculations were performed on a uniformly expanding grid and a fixed grid with a dynamically additive cell.In addition, a model using the method of characteristics to determine the solution near a moving boundary are also be presented.

The ballistic Lagrange problem
The simplest example of propelling a projectile with a gaseous charge is the Lagrange problem, shown in Figure 1.In a cylindrical one-sided open tube with rigid, perfectly smooth and non-heat-conducting walls a projectile is positioned with the mass mp .The space between the projectile and the wall enclosing the tube is filled with a gas with density 01 and pressure p01 .The open side of the tube is filled with gas at a pressure p02 less than p01.The projectile is maintained in a fixed position by means of a diaphragm, which bursts when a forcing pressure equal to p01 is reached in the chamber.The projectile starts to move under the influence of the force generated by the difference in gas pressure values on both sides of the projectile.The initiation of projectile motion generates a mechanical disturbances in the gas.These disturbances propagate as a rarefaction (dilution) wave in the gas behind the projectile.Initially, the waves travel in the region of undisturbed gas.Such waves are called 'straight' waves.The rarefaction wave reaches the wall and reflects off it as a dilution wave.It interferes with the wave propagating from the bottom of the projectile.The spatial distributions of the gas velocity, pressure, density and temperature at any given time are determined by the wave interference processes.The classical Lagrange problem neglects diffusive transport processes in the gas (viscosity, thermal conductivity) and treats the flow as one-dimensional.It is assumed that during the whole process of propulsion, the pressure of the gas in front of the projectile is negligibly small compared to the pressure acting on the bottom of the projectile.Under these assumptions, the mathematical formulation of the Lagrange problem includes the differential equations of gas dynamics: • mass balance equation: !
• momentum balance equation: • energy balance equation: where is the density of the gas, ug is the velocity of the gas, p is the pressure of the gas, u is the specific energy of the gas, x is the spatial coordinate and t is time.Caloric equation of state of a perfect gas with constant specific heat: where is the adiabatic exponent.Initial conditions for the problem are following: where vp is the velocity of the projectile and xp is the position of the projectile.Boundary condition on the fixed wall can be formulated as: boundary condition at the bottom of the projectile is formulated as: which is fulfilled by the equation of motion of the projectile: where S is the area of the base of projectile, Projectile trajectory equation is described by the following equation:

Reference solution -the method of characteristics
For the reference solution of the classical Lagrange problem, a simple method based on the differential approximation of the relationships on the characteristics is suitable.The approximate solution of the flow equations is constructed on the grid of characteristics shown in Figure 2.
where ugA , ugB , ugC are the gas velocity values at points A, B and C respectively, pA , pB , pC are the gas pressure values at points A, B and C respectively and cA , cB , are the values of speed of sound at points A and B respectively.From the above equations, the velocity and pressure at point C can be determined.The values of the other parameters are determined from the formulae: In the corrector step, values of sound velocity should be substituted into equation (10) and equation (11) as average values on the given characteristic sections and the above parameters recalculated.
The position of point C can be determined from the approximation of the equations of characteristics: where xA , xB , xC are the spatial coordinates of points A, B and C and tA , tB , tC are the temporal coordinates of points A, B and C.
When determining the values of the parameters at the boundaries of the domain, the above procedure should be applied taking into account the boundary conditions.To determine the time coordinate of point C on the moving boundary, equation ( 14) and the differential approximation of the projectile trajectory equation should be used: where vBC is the average velocity of the projectile over the distance between points B and C.

Numerical formulations -grid and numerical boundary conditions
Numerical calculations of the different approaches to determine the parameters on the moving boundary were performed on the uniformly expanding ALE grid (Arbitraty Lagrangian-Eulerian grid) shown in Figure 3, and on the fixed grid (Eulerian grid) with dynamically expanding cell shown in Figure 4.The cell splits when its length exceeds 1.5 standard size values.In both approaches, the determination of the parameters at the moving boundary was based on the solution of the Riemann problem using a fictitious cell with parameters: The author's approach to determining the parameters at the moving edge was based on an Euler grid. Figure 5 shows a section of the grid near the projectile trajectory for two time layers j and j+1.To calculate the values of the parameters at n and n+1 nodes on layer j+1, relationships on the characteristics are used.For a one-dimensional thermally and calorically ideal gas flow, these have the following form: where c is the speed of sound.
The equations of characteristics are of the form: The numerical solution uses a differential approximation of the above equations on a section of the characteristics: The quantities in sharp-angled brackets make sense of the average values over a given section of the characteristic.The differential equations of characteristics have the form: Figure 6 shows a section of the characteristic reaching node (n+1,j+1) -point C.This characteristic intersects the time layer j at point A. This point lies between nodes (n,j) and (n+1,j).We can find the parameter values at this point by the linear interpolation.Figure 7 shows the characteristic sections arriving at node (n,j+1) -point C.  Parameter values at node C were calculated using the predictor-corrector procedure, the differential form of the relations on the characteristics (20) and the differential form of the characteristics trajectory equations ( 21), supplemented by the boundary condition equation ( 7), the projectile motion equation ( 8) and the projectile trajectory equation ( 9).In the predictor step, the parameters are assumed at individual nodes as values from the nearest nodes with known parameters.By calculating the coordinates of the nodes, the parameter values are found by linear interpolation.The gas velocity and pressure are thus determined.The density of the gases is then determined from equation ( 12) and the speed of sound from equation (13).In the corrector step, the above procedure is repeated, substituting the previously calculated values at the nodes.Finally, the internal energy is determined using the following expression: For any formulation of the boundary condition with a moving boundary, an analogous procedure was used with a stationary boundary.

Numerical diagrams
Lax-Wendroff-type schemes were used to find a solution within the computational domain: the Richmyer scheme and the MacCormack scheme.The Richtmyer scheme consists of determining the parameter values between nodes in a half time step in the first step and then using them in the next step: where U is a vector of conservative parameters and F is a vector of fluxes.The MacCormak scheme, on the other hand, uses a first-order forward differential scheme in the predictor step and a backward differential scheme in the predictor step: (27)

Test data
In order to test the algorithms and schemes presented above, the classical Lagrange problem was solved for the two initial data sets shown in Table 1.

Figure 1 :
Figure 1: Illustration for the formulation of the classical Lagrange problem.

Figure 2 :
Figure 2: Characteristics grid for determining the approximate solution of the classical Lagrange problem.

Figure 3 :
Figure 3: A uniformly expanding ALE grid.Figure 4: Euler grid with dynamically added cell.

Figure 4 :
Figure 3: A uniformly expanding ALE grid.Figure 4: Euler grid with dynamically added cell.

Figure 5 :
Figure 5: Illustration of a section of the grid near the moving boundary.

Figure 6 :
Figure 6: Illustration of the method for calculating parameter values at a node on the boundary.

Figure 7 :
Figure 7: Illustration of the method for calculating parameter values at a node near the boundary.(a), (b) -subsonic flow, (c) -supersonic flow.

Figure 8 :
Figure 8: Illustration of the method for calculating parameter values at a new node near the boundary.(a) -subsonic flow, (b), (c) -supersonic flow.
p01 (MPa) 01 ( mp ( l (m) S ( Reference solutions to the Lagrange problem for the subsonic case are shown in Figure9and Figure10, and for the supersonic case in Figure11and Figure12.Dashed lines denote sound speed distributions.

Figure 9 :
Figure 9: Gas pressure distribution for the subsonic reference case solution.

Figure 10 :
Figure 10: Gas velocity distribution for the subsonic reference case solution.

Figure 11 :
Figure 11: Gas pressure distribution for the supersonic reference case solution.

Figure 12 :
Figure 12: Gas velocity distribution for the supersonic reference case solution.