Development of a Method of Extended Cells for the Formulation of Boundary Conditions in Numerical Integration of Gas Dynamics Equations in the Domains of a Curvilinear Shape

A method of extended cells for the formulation of boundary conditions in the numerical integration of the Euler equations in domains with curvilinear boundaries for the cases of one-dimensional and two-dimensional compressible gas flows has been proposed in this study. The proposed method is based on the use of the explicit Godunov type finite volume method on a regular rectangular Cartesian grid. The essence of the extended cell method is that when integrating the basic equations of gas dynamics in a fractional cell, numerical fluxes are calculated through the sides of the new extended cell. This new cell is constructed tangentially to the curved boundary and has a size no less than the cell size in the regular domain. Parameters inside the new cell are calculated as mean integral values over the area of the neighboring regular cells included in it. In this case, when choosing a time step in accordance with the Courant condition, the stability of the method in the main computational domain is preserved both when integrating fractional cells and when integrating regular cells. Thanks to this approach, the proposed method has low requirements for computing resources, the ability to generalize for three-dimensional space and increase the order of accuracy without major modifications of the algorithm.<br><br>To test the proposed method, solutions were obtained for the generally accepted test problems of gas dynamics: normal and double Mach reflection of a shock wave from a plane wall. The choice of the time step was made in accordance with the Courant condition in regular finite volumes. The results obtained have made it possible to assess the convergence of the proposed method and their comparison with the results of calculations using other methods have shown a good quantitative and qualitative agreement.


Introduction
The process of designing modern technical equipment includes a large number of various stages. One of the key ones is obtaining the fields of physical parameters resulting from the operation of the designed product. This is made by solving problems of computational gas dynamics in the geometry of the domain of interest. This allows one to make more exact the product design and make an assumption about its characteristics. The solution to such problems is associated with numerical integration of basic equations of gas dynamics in domains of a complex shape with curvilinear boundaries. The main ways to solve such problems include the application of the finite volume method on structured and unstructured grids and solution of basic equations of gas dynamics written in curvilinear coordinates [1,2]. Each of them has its own advantages and disadvantages when modeling physical processes in domains with curvilinear boundaries but high complexity and low efficiency are their common disadvantages. This results in overestimated requirements for computing power and more longer calculation time which is unacceptable in the context of the rapid development of engineering and technology as well as high competition. Therefore, it is an urgent task to develop a simple and economical method for solving gas dynamics problems in the domains with curvilinear boundaries.

Literature review and problem statement
The study results presented in [3,4] are a reasonable alternative to the above approaches. They show that the use of the finite volume method in problems with curvilinear boundaries on regular Cartesian rectangular grids can give correct results within the computational domain. However, the complexity of this procedure lies in the fact that it is necessary to apply special relations in finite volumes outside the computational domain. Two approaches can be distinguished here: immersed particles [3] and fractional cells [4]. However, the absence of conservative property and implicit consideration of the body boundaries allow us to exclude the first of them from further consideration in this work [5].
It was shown in [4] that until recently, the issue of ensuring the stability of computations due to the presence of small cells in the computational domain still does not have an unambiguous solution. One of the common ways to ensure stability implies the application of the cell merging method characterized by a decrease in accuracy of the solution near the boundary and difficulty of obtaining a universal grid generation algorithm [6]. The "h-box" method can be an option of overcoming these difficulties [7,8]. What is proposed in this method is to integrate fractional finite volume with fluxes calculation performed using several additional "h-cells" constructed on each of the fractional finite volume faces. The solution stability is maintained here due to the property of reducing the difference in fluxes calculated with the use of "h-cells". However, the extension of the "h-box" method to the multidimensional case will still be ineffective because of the complexity of calculating fluxes through one face of the fractional finite volume. To this end, it is necessary to introduce a new coordinate system and build four "h-cells" in it and then find the flow parameters in each of them. After that, one should solve two Riemann problems along the corresponding directions and calculate flux through one face on this basis. Implementation of such an algorithm at each time step is resource-intensive even when using schemes of the first order of accuracy.
Thus, the analysis of existing studies shows that the existing methods do not ensure a universal economical algorithm for solving the computational gas dynamics problems on regular rectangular Cartesian grids in domains with curvilinear boundaries. All this makes it possible to assert that it is expedient to conduct a study devoted to the development of an effective method for formulation boundary conditions for solving problems of the computational gas dynamics.

The aim and objectives of the study
The study objective is to develop an effective method for the formulation of boundary conditions when solving problems of computational gas dynamics using the explicit finite volume method on regular rectangular Cartesian grids. The proposed method should be universally extended to multidimensional cases and allow an increase in the order of accuracy.
To achieve the objective, the following tasks were set: -to develop a method for integrating equations of gas dynamics in boundary cells; -to generalize this method for the case of two spatial variables; -to verify the proposed approach by solving a series of test problems of the dynamics of a inviscid compressible gas.

The method of numerical solution of the system of equations for the dynamics of a inviscid compressible gas
The system of two-dimensional unsteady Euler equations in an integral form closed by the equation of state [1,2,9] was used as the basic model of an ideal inviscid compressible gas. The vector form of the used system of equations is as follows: where U is the vector of conservative variables; F, G are the flow components of quantity U along axes of coordinates OX, OY, respectively; ρ is density; p is pressure; u, v are velocities of the flow along the axes of OX and OY coordinates, respectively; e is specific internal energy; E is specific total energy; k is the gas heat capacity ratio.
The slip wall boundary conditions were set on the solid wall [2,9]: are the velocities along tangential and normal directions to the wall.
The explicit Godunov type finite volume method was used for numerical integration (1) [1,2,9] which consists of the following sequential stages: 1) reconstruction of the flow parameters was performed by means of piecewise constant functions on boundaries of a finite volume by their mean values inside the volume. For example, reconstruction of the first order of accuracy for a finite volume with indices i, j will take the form: where U i,j is the mean value of the vector of conservative variables inside the finite volume; , , 2) solution of the Riemann problem using the Lax-Friedrichs relations [10] on the boundaries of a finite volume and calculating the fluxes through them. For the right face of the finite volume with indices i, j, the ratio takes the form: is the solution of the Riemann problem on the right face of a finite volume based on which flux ( ) through this face is calculated which is necessary for further search for parameters in a finite volume; λ max is the maximum wave propagation velocity in the problem under con-sideration. In this and the following notations, half-integer indices refer to the solution of the Riemann problem on the corresponding faces.
3) time integration using the Euler's explicit method was performed. In this case, equation (1), taking into account (3), (4) in a semi-discrete form, will be written as follows: where A i,j is the area of the final volume, Δy i,j , Δx i,j are the height and width of the final volume, respectively. Numerical modeling of slip wall boundary conditions (2) was carried out using "dummy" cells [2,9].

1. The one-dimensional case
Let us consider a one-dimensional problem of evolution of gas flow in a semi-bounded region Ω={0≤x≤L, 0≤t≤T} with an impenetrable wall at the point x=L.
Divide the domain Ω into N equal finite volumes according to the following relation in accordance with Fig. 1: where x i is coordinate of the right side of the i-th finite volume; h is the discretization step in space in the regular domain.
The spatial discretization step is chosen so that the following condition was fulfilled: As seen from Fig. 1, with this step choice, the last cell goes beyond the boundary of the calculation domain for the length h-h' where h' is the spatial discretization step near the wall and For the one-dimensional case, equation (5) for the i-th finite volume is transformed: As is known, a necessary condition for stability of explicit finite-volume methods is the Courant-Friedrichs-Levy (CFL) condition [10,11]: where C is the Courant number.
In accordance with the CFL condition, the time step Δt for the inner cells is determined as follows: where K≤1 is a certain stability factor determined for each problem being solved. It is easy to show that if the condition K/α>1 is satisfied, the CFL condition in the last cell will be violated leading to the algorithm instability. The idea of the proposed approach is to construct an extended cell N' of width h including the cell N. In this case, the right border of the cell N' coincides with the wall and the left border is inside the cell (N-1).

Fig. 2. Structure of an extended cell
Obviously, the CFL condition (10) will be fulfilled in cell N' at the selected time step (11).
Equation (9) for the cell N' has the form: As seen from Fig. 2, the state in the new cell N' can be determined by calculating the mean integral value along the length h of the new cell from the piecewise constant functions U N-1 and U N : Now, when knowing the value of the vector of variables in the extended cell, in order to calculate (12), it remains to find the value of the vector of variables in the dummy cell with index (N '+1) using boundary conditions (2). In accordance with [7], introduce the reflection operator ℜ which will have the following form in the one-dimensional case: Then the value of the vector of variables in the dummy cell is found from the expression Using (13) and (14), find solutions of the Riemann problem on faces of the cell N` from formulas (3) and (4) ; Following substitution of (15) into (12), calculate fluxes, integrate (12), and store the parameters obtained in the new cell in the last cell with index N.
Since the integration of the last finite volume uses the size h of the new cell and corresponding fluxes, the resulting solution remains stable if condition (10) is fulfilled. Note also that calculations by (13) and (14) are performed at the first stage of the finite volume method and the parameter α is determined only once when generating the grid. Due to all this, the proposed method can be modified to obtain a higher order of solution accuracy.

1. The two-dimensional case
Let us generalize the proposed approach for the two-dimensional case. Consider, for example, the two-dimensional flow of a inviscid compressible gas in a rectangular domain {0≤x≤L x , 0≤y≤L y , 0≤t≤T}, a part of which under the straightline f(x)=kx+b is occupied by a solid wall. Divide the calculation domain into finite volumes with a regular rectangular Cartesian grid as shown in Fig. 3: where h x is the step of discretization along the ОХ axis; h y is the sampling step along the ОY axis; h x , h y are the numbers of finite volumes along the ОХ and ОY axes, respectively.
As seen in Fig. 3, the grid cells can be divided into three types: 1) regular cells: the cells entirely lying in the calculation domain. Their area is calculated by the formula: 2) fractional cells: the cells in which one part of the area lies in the calculation domain (A c ) and the other is inside the solid wall (A w ). The following relations are valid for them: 3) cells that are completely behind the solid wall. These cells are not processed in the calculation and will not be mentioned further.
In the two-dimensional case, semi-discrete equation (5) is used for numerical integration over a regular finite volume with indices i, j. The time step is chosen so that it satisfies the CFL condition for the two-dimensional problems in the regular domain [10,11]: are the maximum velocities of wave propagation along the ОХ and ОY axes, respectively.
In this case, the standard method becomes unstable only when integrating fractional cells in which A i,j =A c <A r . Obviously, to fix this, one needs to increase the size of the fractional cells. This can be done by constructing a new extended finite volume with an area , .
i j r A A ≥ ′ Let us consider some arbitrary fractional cell in Fig. 4 with indices i, j. As in the one-dimensional case, construct an extended cell for the selected fractional cell in accordance with the following technique: 1) find coordinates of points of intersection of the wall line and the fractional cell (points D, F in Fig. 4); 2) find a midpoint between the obtained coordinates and select the found point as origin O' of the new coordinate system; 3) select the tangent to the wall line as the horizontal axis O'τ of the new coordinate system (in the case under consideration, this is the straight line f(x)=kx+b) and the direction perpendicular to it as the vertical axis O'n; 4) define in the new coordinate system the cell vertices for which the coordinate in the O'n direction will be strictly positive (points A, B, C in Fig. 4). Among them, find the maximum coordinate c n along the normal to the wall (point B') as well as the maximum (point C ') and minimum (point A') along the wall. All points considered in this section are shown in Fig. 4; 5) redefine the found coordinates according to the following relations: 6) build an extended cell in a new coordinate system as a rectangle with the lower-left coordinate ( min , c τ 0) and the upper right coordinate ( max , c τ c n ). In this case, the area of the constructed cell will be calculated as follows: where S is the area of the polygons resulting from the intersection of cells in the calculation domain with a new cell. Fig. 4 shows the constructed extended cell and the corresponding fictitious one inside the wall. Let us write the relation (4) for the right face of the extended cell and then the semi-discrete equation (5) for the entire extended cell: where symbols , , where l i+1,j (l i,j+1 ) is the length of the section of the face of the turned cell which is inside the finite volume with indices (i+1, j) (i, j+1).
Reconstruction of the vectors of variables on the remaining faces of the new extended cell, except for the face that lies on the wall, is performed in a similar way.
As in the one-dimensional case, to calculate the value of the vector of variables in a dummy cell, introduce a two-dimensional reflection operator ℜ which will take the form: where u τ , u n are the velocities along the O'τ, O'n axes of the new coordinate system, respectively. With its help, determine the value of the vector of variables in a dummy cell as 1, Let us note the advantages of the proposed approach: 1) in contrast to the fractional cell method, the stability condition is not violated when calculating parameters in the boundary cells; 2) the results of integration over the extended finite volume are saved in the corresponding fractional volume while calculating four fluxes instead of five, as in the "h-box" method; 3) regardless of the geometry of the fractional cell, one extended finite volume is built for it while the number varies from 9 to 17 in the "h-box" method; 4) at each time step, the transition to a new coordinate system occurs only once when calculating parameters of a dummy cell compared with the "h-box" method in which this is done twice for each face of a fractional cell for which the flux coming through it is calculated; 5) calculation of parameters in the extended cells is performed at the reconstruction stage and their number can be easily increased by making larger the template. This ensures the direct application of the methods with a high-order accuracy without significant modification of the main algorithm.

The results obtained in test calculations
To verify the proposed method of formulation of boundary conditions, two test problems of dynamics of inviscid compressible gas were solved. These are the problem of normal reflection of a shock wave from a solid wall and the problem of double Mach reflection of a shock wave from an inclined wall. Note that when solving all problems, the time step was chosen according to relation (14) at C=0.9 and the calculation was stable in all cells of the computational domain.

1. Normal reflection of the shock wave
To assess the effect of using the formulation of boundary conditions with the help of extended cells on the convergence of the method, the problem of normal reflection of a shock wave from a plane wall was solved. The solution was sought in a two-dimensional region {0≤x≤4.5, 0≤y≤3.0, 0≤t≤0.3} shown in Fig. 5 where subscript SH refers to parameters of the gas behind the shock wave.
The calculation was carried out on regular grids with discretization steps in the OX, OY directions corresponding to 40, 80, 160, and 320 cells along each coordinate axis. This problem is one-dimensional in the Cartesian coordinate system associated with a solid wall. This has made it possible to find an exact solution to the selected problem using the one-dimensional theory of shock waves [12]. To reduce the results of a numerical solution to the one-dimensional case, a square normal to the wall with a side equal to 1.5 units was selected in the middle of the calculation domain with its base lying on the wall. Further, density and internal energy domains were averaged in this square, along the direction tangential to the wall in such a way as to obtain their dependences on the distance along the normal to the wall. The convergence rate was estimated by comparing the numerical and analytical values of the gas parameters behind the reflected wave in finite-dimensional analogs of the norms C, L 1 , L 2 [13]. Tables 1, 2 show dependences of the relative calculation error Δ of density, the internal energy of the gas in the norms C, L 1 , L 2, and the estimate of the order of accuracy n of the method of the numerical solution. There are also given values of errors of flow parameters when solving this problem in a one-dimensional formulation for the classical case of formulation of the boundary conditions (the solid wall coincides with the boundary of a finite volume).  Table 2 Relative error of the proposed approach calculated from the flow energy N 1D

2. Double Mach reflection of the shock wave
To compare qualitative patterns of supersonic flow, the problem of double Mach reflection was solved in two formulations.
In the first case, classical formulation of the problem described in [14] was used: a shock wave runs at a certain angle onto a horizontal wall which coincides with the lower boundary of the calculation domain and starts at the point x=0.1667. Moreover, all finite volumes in the region are regular. This formulation of the problem is shown in Fig. 6.
The density domain obtained as a result of calculations on a 420×180 grid up to the time moment t=0.2 s is shown in Fig. 7.
Using extended cells, this problem was solved in the domain {0≤x≤3.2, 0≤y≤1.8, 0≤t≤0.3} the part of which located below the straight line was occupied by the wall, as shown in Fig The density domain obtained in the calculation is shown in Fig. 9. For a quantitative comparison of the obtained density domains, the results of the extended cell method were combined with the results obtained by the method [14]. Fig. 10 shows the isolines of the domain of modulus of the difference of these densities. Fig. 10. The modulus of difference between the density domains obtained in solving the problem of double Mach reflection in the classical formulation of the problem and when using extended cells As seen from Fig. 10, the results differed only on the discontinuity surfaces. The geometry of the disturbed region was the same in both cases and the small discrepancy of the parameters within it was within the error limits of both methods.

Discussion of the results obtained in using the method of extended cells in the explicit method of finite volumes
Persistence of stability of the explicit finite volume method when integrating boundary cells is explained by the fact that size of the extended cell was used instead of size of the fractional cell. In this case, the CFL condition calculated over the regular region was not violated.
1. New extended cells with a size greater than or equal to the size of cells in the regular domain were used to integrate fractional boundary cells.
2. Compliance with the algorithm of the explicit finite volume method. The use of relations of the proposed method at the stage of reconstruction makes it possible to apply the same methods to the boundary cells for increasing the order of accuracy, calculating fluxes, and integrating as for cells in a regular domain.
3. Reduction of calculations: instead of nine new cells [7,8], it is necessary to build and process only one extended cell.
It should be noted that only flat boundaries were considered in the current study. Therefore, the application of the proposed method for convex or concave geometry requires further development efforts.
Also, the disadvantage of the proposed method consists in that after the integration stage, parameters of the extended cell are written into a fractional cell. Consequently, calculation results may be less physically real in the case of using a more complex geometry of the problem. The direction of the study associated with the elimination of this drawback should be focused on finding values in fractional cells after the integration of the gas dynamics equation in extended cells.
In the next studies, it is planned to generalize the proposed approach to domains with curvilinear boundaries, increase the order of accuracy of the method, and generalize it to a three-dimensional space.

1.
A new method was proposed for the formulation of boundary conditions of numerical solution of gas dynamics equations in domains of a complex shape on Cartesian regular grids. The use of extended cells requires several times less amount of computations to determine flow parameters in a fractional cell compared to the existing approaches.
2. A technique has been developed for using the extended cell method in calculating two-dimensional flows of a inviscid compressible gas. The proposed method was generalized to the case of two spatial variables by using additional geometric relationships on the faces of extended cells in a local coordinate system.
3. Verification of the proposed method was carried out by solving the problems of double Mach reflections and normal reflections of shock waves from a flat wall. The qualitative pattern of the solution was in good agreement with the results obtained by applying earlier known methods to solve these problems. The magnitude of the error introduced by the proposed method did not exceed in orders of magnitude the error that occurs when using existing approaches. Since the volume of computations, in this case, is significantly less (as noted, up to nine times in the two-dimensional case), it can be concluded that the extended cell method has a good efficiency.

Introduction
It is known that the resistance of submarines at full immersion consists of friction resistance, shape, roughness, cuts, protruding parts, and fencing of retractable devices. Until recently, the only means of determining resistance was a model experiment in research basins. Recalculation of the results of such experiments should be carried out according to the