Научная статья на тему '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'

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 Текст научной статьи по специальности «Медицинские технологии»

CC BY
43
8
Читать
i Надоели баннеры? Вы всегда можете отключить рекламу.
i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
Предварительный просмотр
i Надоели баннеры? Вы всегда можете отключить рекламу.

Текст научной работы на тему «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.

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

Keywords: computational gas dynamics, finite volume

method, Cartesian grids, extended cells -□ □-

Received date 11.09.2020 Accepted date 02.10.2020 Published date 23.10.2020

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

UDC 519.6+533

[DOI: 10.15587/1729-4061.2020.213795|

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

I. Dubrovskiy

Modeling Engineer FLIGHT CONTROL LLC Gagarina ave., 115, Dnipro, Ukraine, 49050 E-mail: iskolyar@gmail.com V. Bucharskyi PhD, Associate Professor Department of Engine Construction Oles Honchar Dnipro National University Gagarina ave., 72, Dnipro, Ukraine, 49010 E-mail: bucharsky1@gmail.com

Copyright © 2020,1. Dubrovskiy, V. Bucharskyi This is an open access article under the CC BY license (http://creativecommons.Org/licenses/by/4.0)

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.

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

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

-d J Udi + <J( F (U) + G (U ))d5 = 0,

'p ^ pu N

u = pu , F = pu2 + p

pv puv

,pe , (pE + p) u/

G =

(1)

pv

puv

pv2 + p v(pE+p) v7

E = e +1 (u2 + v2), P = (k -1) ep,

where U is the vector of conservative variables; F, G are the flow components of quantity U along axes of coordinates OX, OY, respectively; p 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]:

dV.

dT

= 0, V = 0,

(2)

where VT, Vn 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:

ul. = u* = ur. = ub. = u ■,

t,1 t,1 t,1 t,1 t, ■

(3)

where Ui,j is the mean value of the vector of conservative variables inside the finite volume;

U

u * ■,

tj'

UJ-, U *

values of the vector of conservative variables on the left, right, upper and lower faces of the finite volume with indices i,j, respectively;

2) solution of the Riemann problem using the Lax-Fried-richs 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 j, the ratio takes the form:

4. 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:

U^ = \(F (UR) +F (Ut, (UL1, j - UR ) (4)

where u*+ y is the solution of the Riemann problem on the

right face of a finite volume based on which flux F (u*+ ^)

through this face is calculated which is necessary for further search for parameters in a finite volume; Xmax 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:

3U.

dt

A

Ay'j ( F K >2 )-F (U >2 )

j ( G (—+ >2 )-G (u- >2 )

xl = i • h, i = 1...N,

L , L — < h <-.

N N -1

h =ah <h, a<1.

x<L (domain of the gas flow)

C =

At -I

-< 1,

(10)

where C is the Courant number.

In accordance with the CFL condition, the time step At for the inner cells is determined as follows:

At = K

h

(11)

(5)

where Ai,j is the area of the final volume, Ayi,j, Axi,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].

5. The method of integrating the equations of gas dynamics in boundary cells

5. 1. The one-dimensional case

Let us consider a one-dimensional problem of evolution of gas flow in a semi-bounded region fl={0<x<L, 0<t<T} with an impenetrable wall at the point x=L.

Divide the domain fl into N equal finite volumes according to the following relation in accordance with Fig. 1:

where K<1 is a certain stability factor determined for each problem being solved.

It is easy to show that if the condition K/a>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).

N-3

N-2

N-1 N

. h

(6)

where xi 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:

N N'+1

Fig. 2. Structure of an extended cell

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

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:

(7)

——N-=-1 If IUI

dt

h

- F|U

N '- K

(12)

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

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 Un-1 and Un:

(8)

Un,= 1 J U (x)dx =

hL-h

= h (h'UN +(1 - h')UN-1)

= aUn + (1 -a)UN-1.

! h h' h-h'

( • • • • • • •

x>L (wall)

(13)

N-6 N-5 N-4 N-3 N-2 N-1

Fig. 1. Discretization of the calculation domain

For the one-dimensional case, equation (5) for the i-th finite volume is transformed:

f = -h (F K >2 )-F K

(9)

ow, 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:

As is known, a necessary condition for stability of explicit finite-volume methods is the Courant-Friedrichs-Levy (CFL) condition [10, 11]:

u' = r[u] = r

pu

pE

P

-pu

PE.

(14)

1

X

X

L

x

N

Then the value of the vector of variables in the dummy cell is found from the expression Un+1=^[Un'].

Using (13) and (14), find solutions of the Riemann problem on faces of the cell N^ from formulas (3) and (4)

ur -u • UL = U • Ul = U • UR -U

^n '-,j n' • n,-fj n' • n '+1 -,J n '+1; n-1 = ^ N-1,

u;+>2 -2(F(UR)+F(ut'+1 ))-^^(uL'+1 -uR,), (15) u; _yt - 2 (f (un-1)+f (un. ))-^^(ul' - ur-1).

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

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:

A - K hy; (17)

2) fractional cells: the cells in which one part of the area lies in the calculation domain (Ac) and the other is inside the solid wall (Aw). The following relations are valid for them:

Af - Ae + Aw,

Ac < Ar; (18)

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]:

5. 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<LX, 0<y<Ly, 0<t<T}, a part of which under the straight-line 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:

xi- ihx, y,- jK,

At <-

i -1...N, j -1... Ny,

L Ly

hx - NX' hy - N'

(16)

where hx is the step of discretization along the OX axis; hy is the sampling step along the OY axis; hx, hy are the numbers of finite volumes along the OX and OY axes, respectively.

"^J

X \

1 1'

02 \

\ >V \\ \

\ \\

\ \ \\

0 L„

Fig. 3. Two-dimensional calculation domain: regular cells fractional cells 2; empty cells 3

(19)

where , are the maximum velocities of wave propagation along the OX and OY axes, respectively.

In this case, the standard method becomes unstable only when integrating fractional cells in which Ai,j=Ac<Ar. 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 A' > A .

'j r

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'x 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 cn 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:

^[cn, hx, hy ],

C

>/(C)2 + (cr)2 ^ max[hx,hy],cfx = 0.5max[h,,h], Cn =-0.5max [h,, hy ], (20)

^(C )2 +(cr )2 > max [hx, hy ] ,C

c. = c. , c. = c. ;

T T ' T T '

6) build an extended cell in a new coordinate system as a rectangle with the lower-left coordinate (c™n, 0) and the upper right coordinate (cm™, cn). In this case, the area of the constructed cell will be calculated as follows:

U* u denotes solution of the Riemann problem on the same face. The value U/j is determined using the piecewise-con-

stant reconstruction u/j = U'-, and U^ by calculating the mean integral value along the length of the face from the piecewise constant functions that fell on it:

U-=(CFbr) !U <I)dI=

1

\ (li+1,jUi+1,j + li,j+1Ui

i, j+1/

(26)

4,,j' = (cr+C11 )cn

It is obvious that: A. . > A ;

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

i, j C'

(21)

(22)

7) determine value of the vector of variables in the new cell by calculating the mean integral value over area of the new cell from the piecewise constant functions Uij, Ui+ij, Ui.1j, Uij+1, Ui j.1 found inside:

Ui,j = AT -U (S )d S =

Ai j S

A.

Si+1,jUi+1, j + Si, j'+1Ui, j+1 + +Si-1, jUi-1, j + Si, j-Ui ,j-1 + Si, jUi,

(23)

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.

Fig. 4. The construction of a new turned cell for a finite volume with indices i, j

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:

1

i+1, j

-%-(u L+1 j - U Rj )

au

i j dt

cn ( F K >2 )-F K >2 )) +

where symbols UR

(C" + c"

U L' i+1, j

G(U *+K )-G (U * -X

(24)

, (25)

where li+1yj (li 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:

R [U ] = R

P

PUT PUn

PE

P

PUT -Pun

. PE

(27)

denote values of the vectors of

variables on the right face of the new cell and the symbol

where ux, Un are the velocities along the O't, 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 R [Ui'+1]. Next, solve the Riemann problem on the face which lies on the wall using values of the vectors of variables Ufj = U'j and R [Ui'+1 - ]. Then integrate (25) and save the result obtained as a parameter in a fractional cell with indices i, j. Since integration uses the area A'- of the extended cell and the fluxes corresponding to it, the resulting solution remains stable under condition (19).

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.

6. The results obtained in test calculations

To verify the proposed method of formulation of boundary conditions, two test problems of dynamics of inviscid

1

2

1

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.

6. 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<i<0.3| shown in Fig. 5. The part of the domain below the straight line

f ( X ) = f tg|] X - 0.09624,

was occupied by the wall. Initial conditions in the entire domain were set equal to parameters of the gas behind the shock wave with the Mach number of the shock wave Mw=10:

U (0, x, y ) =

' Psh ' f 7.9 ^

PSHUSH 32.34

PSHVSH -56.21

KpSHESH , v 551.7 ,

y >| tg^ Ix - 0.09624

U (0, x, y ) =

Psh " f

PSHUSH 0

PSHVSH 0

PSHESH , v 0;

(28)

y <| tg-6 I x - 0.09624,

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;l , L2 [13]. Tables 1, 2 show dependences of the relative calculation error A of density, the internal energy of the gas in the norms C, L\, L2, 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 1

Relative error of the proposed approach calculated from the flow density

N 1D

AC nAC An nAL1 Al2 nAL2

40 0.032 - 0.0188 - 0.0207 -

80 0.022 0.56 0.0075 1.32 0.0106 0.97

160 0.015 0.53 0.0038 0.98 0.0063 0.75

320 0.010 0.52 0.0019 1.04 0.0037 0.78

- 2D

40 0.038 - 0.0316 - 0.0319 -

80 0.028 0.46 0.0204 0.63 0.0210 0.61

160 0.020 0.45 0.0131 0.63 0.0139 0.60

320 0.016 0.37 0.0092 0.51 0.0098 0.50

Table 2

Relative error of the proposed approach calculated from the flow energy

N 1D

AC nAC AL1 nAL1 Al2 nAL2

40 0.0316 - 0.0181 - 0.0203 -

80 0.0217 0.54 0.0075 1.28 0.0106 0.94

160 0.0150 0.53 0.0038 0.98 0.0063 0.75

320 0.0105 0.52 0.0018 1.04 0.0037 0.78

- 2D

40 0.0446 - 0.0174 - 0.0182 -

80 0.0141 1.67 0.0083 1.07 0.0093 0.96

160 0.0090 0.64 0.0037 1.18 0.0047 0.99

320 0.0055 0.70 0.0023 0.69 0.0026 0.84

Fig. 5. Formulation of the problem of normal reflection of a shock wave

As can be seen from the values in Tables 1, 2, values of the errors in different norms, and the values of the order of approximation were close in magnitude. The order of accuracy

did not deteriorate when the proposed method was used for the formulation of the boundary conditions.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

6. 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 420x180 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. 8. Initial conditions in the entire domain were set according to the Mach number of the shock wave Mw=10:

Fig. 7. Density domain when solving the problem of double Mach reflection in the standard formulation

u (0, x, y ) =

Po

Pouo

Po V0

vPo E0 J

0 o

v 2.5j

y >| tg- Ix - 0.09624,x > 0.1667

u (0, x, y ) =

x < 0.1667

U (0, x, y ) =

Psh

pshush pshvsh vpshesh

Psh

pshush pshvsh vpshesh /

y <| tg- |x-0.09624

7.9 64.9 0

551.7

1.4 0 0

v 2.5j

Fig. 8. Formulation of the problem of double Mach reflection for verification of the proposed method

(29)

The density domain obtained in the calculation is shown in Fig. 9.

4.56 7.60 10.64 _____ _____ _____

Fig. 9. Density domain in solving the problem of double Mach reflection using extended cells

Fig. 6. Formulation of the problem of double Mach reflection [14]

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.

I p2~pl I

0.0 0.5 1.0 1.5 2.0 2.5 3.0

0.18 0.90 1.62 2.34 3.06 3.78 4.50

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

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.

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.

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

The differences between the proposed method and existing methods [4-8] are as follows.

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.

8. Conclusions

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.

References

1. Moukalled, F., Mangani, L., Darwish, M. (2015). The finite volume method in computational fluid dynamics. An Advanced Introduction with OpenFOAM® and Matlab. Springer. doi: https://doi.org/10.1007/978-3-319-16874-6

2. Blazek, J. (2015). Computational Fluid Dynamics: Principles and Applications. Elsevier, 466. doi: https://doi.org/10.1016/c2013-0-19038-1

3. Chi, C., Lee, B. J., Im, H. G. (2016). An improved ghost-cell immersed boundary method for compressible flow simulations. International Journal for Numerical Methods in Fluids, 83 (2), 132-148. doi: https://doi.org/10.1002/fld.4262

4. Berger, M. (2017). Cut Cells. Handbook of Numerical Methods for Hyperbolic Problems - Applied and Modern Issues, 1-22. doi: https://doi.org/10.1016/bs.hna.2016.10.008

5. Sidorenko, D. A., Utkin, P. S. (2018). A Cartesian grid method for the numerical modeling of shock wave propagation in domains of complex shape. Num. Meth. Prog., 17 (4), 353-364.

6. Bennett, W. P., Nikiforakis, N., Klein, R. (2018). A moving boundary flux stabilization method for Cartesian cut-cell grids using directional operator splitting. Journal of Computational Physics, 368, 333-358. doi: https://doi.org/10.1016/jjcp.2018.04.048

7. Berger, M. J., LeVeque, R. J. (1990). Stable boundary conditions for Cartesian grid calculations. Computing Systems in Engineering, 1 (2-4), 305-311. doi: https://doi.org/10.1016/0956-0521(90)90016-e

8. Berger, M., Helzel, C. (2012). A Simplified h-box Method for Embedded Boundary Grids. SIAM Journal on Scientific Computing, 34 (2), A861-A888. doi: https://doi.org/10.1137/110829398

9. Ferziger, J. H., Milovan, P., Street, R. L. (2020). Computational Methods for Fluid Dynamics. Springer. doi: https://doi.org/ 10.1007/978-3-319-99693-6

10. Chung, T. J. (2010). Computational Fluid Dynamics. Cambridge University Press. doi: https://doi.org/10.1017/cbo9780511780066

11. Courant, R., Friedrichs, K., Lewy, H. (1967). On the Partial Difference Equations of Mathematical Physics. IBM Journal of Research and Development, 11 (2), 215-234. doi: https://doi.org/10.1147/rd.112.0215

12. Loytsyanskiy, L. G. (2003). Mehanika zhidkosti i gaza. Moscow: Drofa, 840.

13. Bucharskiy, V. L. (2007). Metod sovmestnoy approksimatsii postroeniya raznostnyh shem dlya resheniya uravneniy v chastnyh proizvodnyh. Tehnicheskaya mehanika, 1, 129-138.

14. Woodward, P., Colella, P. (1984). The numerical simulation of two-dimensional fluid flow with strong shocks. Journal of Computational Physics, 54 (1), 115-173. doi: https://doi.org/10.1016/0021-9991(84)90142-6

-□ □-

A new method of reducing the resistance of submarines is presented, which consists in installing special circular recesses on its surface in the stern. It is found that during the movement, in the recesses there is a macro-vortex flow, in which pressure decreases significantly. This phenomenon affects the characteristics of the boundary layer and in general the pressure distribution on the surface of the hull, i. e. the resistance of the submarine. Using the methods of computational fluid dynamics, the influence of the number and size of the recesses at their fixed location on the resistance of two types of "Lira" and "Gepard" submarines is investigated. The results show that the decrease in resistance increases with increasing Reynolds number and reaches 6 % for "Lira" with 4 recesses with a diameter of d=0.01 D at Re=1.55.108 and 2 % at Re=1.35.108 for "Gepard" with 7 recesses with a diameter of d=0.01 D. The effect of the number of cells of the computational grid on the results of calculations in the Flow Simulation (USA, France, Canada) and Flow Vision (Russian Federation) software packages was also studied. The effect of resistance reduction obtained in both software packages is approximately the same, but the absolute values differ due to the small number of cells in Flow Vision, which is due to the limited capabilities of the used 2nd version of this complex. There was also a slight effect of resistance reduction on the model of the "Persia-110" (Iran) submarine with recesses during towing tests in the research basin at significantly lower Reynolds numbers. Unlike most resistance reduction means, the use of this method does not require significant changes in the design of the housing. This makes it possible to use it both on new facilities and on facilities that have already been commissioned

Keywords: submarine, macro-vortex means of reducing motion resistance, computational fluid dynamics -□ □-

Received date 29.07.2020 Accepted date 08.09.2020 Published date 23.10.2020

UDC 629.5.01

| DOI: 10.15587/1729-4061.2020.212Ö05|

A STUDY OF THE EFFECT OF RECESSES ON THE MOTION RESISTANCE OF SUBMARINES BY METHODS OF COMPUTATIONAL FLUID

DYNAMICS

J. Bod n arc h u k

Postgraduate Student* E-mail: julka.bodnar4uk@gmail.com Yu. Korol PhD, Associate Professor* E-mail: yuriy.korol@nuos.edu.ua M. Moonesun PhD in Engineering Science, Senior Lecturer Center for Underwater Research Malek Ashtar University of Technology Shahinshahr in Isfahan Iran Isfahan, Isfahan 17469-37181, Iran E-mail: m.moonesun@gmail.com *Educational and Scientific Center Hydromechanics Admiral Makarov National University of Shipbuilding Heroiv Ukrainy ave., 9, Mykolaiv, Ukraine, 54025

Copyright © 2020, J. Bodnarchuk, Yu. Korol, M. Moonesun This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)

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

i Надоели баннеры? Вы всегда можете отключить рекламу.