On New Method and 3D Codes for Shock Wave Simulation in Fluids and

Solids in Euler Variables based on a Modified Godunov Scheme

M. H. ABUZIAROV, E. G. GLAZOVA, A. V. KOCHETKOV, S.V. KRYLOV

Research Institute for Mechanics,

Lobachevsly State University of Nizhni Novgorod,

RUSSIA

Abstract: A three-dimensional technique for modeling shock-wave processes both in fluids and solids and for

modeling fluid-structure interaction problems is proposed. The technique is based on a modified Godunov's scheme

of increased accuracy, which is the same for both fluids and solids, and uses Eulerian-Lagrangian multimesh

algorithms. Improving the accuracy of the scheme is achieved only by changing the "predictor" step of the original

Godunov scheme. A three-dimensional and time-dependent solution of Riemann's problem is used, which provides

a second-order approximation in time and space in the domain of smooth solutions. Monotonicity in the domain of

discontinuous solutions is ensured by the transition to the "predictor" step of the first-order scheme. A similar

solution of the Riemann problem is used at the contact "fluids - solids”. For each body, three types of computational

grids are used with an explicit Lagrangian choice of movable free and contact surfaces. The first type of mesh used

is a Lagrangian surface mesh in the form of a continuous set of triangles (STL file), which is used both to set the

initial geometry of an object and to accompany it in the calculation process, and two types of volumetric three-

dimensional meshes. These are the basic Cartesian fixed grid for each object, and auxiliary movable local Euler-

Lagrangian grids associated with each triangle of the surface Lagrangian grid. The results of numerical simulation

of the processes of the impact of ice fragments on a titanium plate, acceleration by detonation products of

deformable elastoplastic bodies of various shapes, and steel strikers piercing an aluminum plate are presented.

Key-Words: - Godunov scheme, increased accuracy, Riemann problem, multimesh approach, 3D problems,

detonation, elastic-plastic, hypoelastic, impact, penetration.

Received: December 8, 2022. Revised: October 15, 2023. Accepted: November 11, 2023. Published: December 14, 2023.

1 Introduction

Godunov scheme, [1], and its most famous

modifications, [2], [3], [4], [5], [6], are widely used

for solving nonlinear dynamic problems of fluid

dynamics (CFD) in Eulerian variables due to the

ability to distinguish and describe discontinuous

solutions without artificial viscosity. At present,

various modifications of this scheme are also used

to solve problems of the dynamics of a deformable

rigid body (CSD) in Eulerian and Eulerian-

Lagrangian variables. The main problem of

Godunov's scheme is the first order of

approximation and, as a consequence, significant

scheme viscosity, which leads to a fast decay of the

solution. Numerous attempts to eliminate this

drawback for CFDs, which are close in meaning to,

[3], [4], [6], increase the difference stencil of the

scheme and do not provide second-order accuracy

in time in the domain of smooth solutions in the

spatial case, and also create additional difficulties in

the implementation of boundary conditions. In CSD,

when modeling wave processes, the influence of the

schematic viscosity is even more significant, and in

many problems, it is necessary to use variants of the

scheme with an approximation order of at least two.

At present, for CSD, there are a large number of

modifications of the Godunov scheme of increased

accuracy, in particular, works, [7], [8], [9], [10],

[11], [12], [13], [14], [15], [16], [17], [18], [19],

[20]. In these works, modifications of the Godunov

scheme of increased accuracy are developed, based

on various versions of hyperelastic models of rigid

body dynamics. These models are hyperbolic,

invariant to the rotation of a rigid body,

thermodynamically compatible, and can be written

as a system of first-order differential equations in

the form of conservation laws. The resulting

modifications are laborious, are of more academic

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

173

Volume 18, 2023

interest, and have received limited distribution in

computational practice and commercial packages.

For hypoelastic CSD models, the solution to the

problem of increasing the accuracy of numerical

models was presented, [21], and further developed

by him for various models of nonlinear material

behavior in, [22]. It is shown that for hypoelastic

models of media, including those describing

irreversible deformations, and for “predictor-

corrector” schemes with splitting into physical

processes, it is sufficient at the “predictor” stage to

obtain a solution of the linearized equations in the

elastic approximation with the second order of

accuracy. The nonlinear behavior of the material is

taken into account at the “corrector” stage. In this

case, the second order of approximation of the

system of equations as a whole is preserved. In,

[23], [24], a modification of Godunov's scheme for

CSD of the second order of accuracy, monotonic on

discontinuities was proposed. In this case, the exact

solution of the Riemann problem in the elastic

formulation is used for the linearized equations of

the theory of plastic flow by the approach, [21], in a

two-dimensional formulation on a compact stencil.

This modification solved the problem of the

increased scheme viscosity and the problem of

boundary conditions. The increase in accuracy is

achieved due to the convergence of the areas of

influence of the differential and difference problems

of the Riemann problem, the monotonicity of

solutions in the area of discontinuities is ensured by

the transition to the “predictor” of the scheme of the

first order of accuracy. At the "fluid-elastic body"

contact, the exact solution of the Riemann problem

is also used. In, [25], [26], [27], [28], this

modification was generalized to a three-dimensional

case and three-dimensional problems of shock-wave

loading of elastoplastic bodies were solved.

Modeling three-dimensional dynamic processes of

fluid-solid interaction (FSI) in Eulerian variables

also requires an adequate description of complex

processes at moving contact boundaries. Therefore,

it is desirable to highlight and accompany the

moving boundaries in the process of calculations.

Currently, there are two approaches to describe the

spatial motion of free and contact boundaries in

Euler variables. The first approach is (the Sharp

Interface Method - SIM), and the second is (the

Diffusive Interface Method - DIM). The SIM

approach, [29], [30], [31], [32], [33], [34], [35],

[36], [37], involves precise selection and tracking of

the motion of the boundary surface. The best option

is the coincidence of the computational grid with the

boundaries of the body, which is not always

possible in the case of large displacements and

deformations, and, in practice, is possible only in

one-dimensional and two-dimensional cases.

Variants associated with the use of various

algorithms for tracking the location of the contact

surface within moving or stationary Euler grids,

often using the subgrid mesh technique to improve

the accuracy in the most interesting parts of the

computational domain (Adaptive Mesh Refinement

technology -AMR), are also complex and are

successfully applied only to solve two-dimensional

problems. In the 3D case, this approach causes

significant difficulties associated with tracking and

restoring the surface of the bodies itself, dynamic

non-Lagrangian rearrangement of meshes, and the

implementation of boundary conditions. In, [35], a

SIM approach for the 3D case using volume

fractions and solving the Riemann problem of a

"fluid-solid" discontinuity to restore and move the

contact boundary inside cells with a mixture was

proposed. The approach is conservative and

includes AMR, but due to its complexity, it did not

receive further development, even though it

indicated the way to solve the problem. In, [36],

also a SIM variant for the 3D case was proposed.

The authors solve the Riemann problem inside cells

with a mixture by interpolating and extrapolating

parameters from the surrounding cells without a

mixture to formulate and solve the Riemann

problem. Then this solution is used to move the

contact boundary inside the cells with the mixture

and to calculate the fluxes to the surrounding cells,

cutting them by the volume fractions in the cells

with the mixture. Several procedures are iterative.

Also, due to the complexity, the method did not find

further development and application. The work,

[38], is indicative in this respect, evolution from 2D

SIM, [29], [30], [31], [32], [33], to 3D DIM, [38],

[39]. It concludes the practical inapplicability of

SIM for 3D problems. The second DIM approach,

[40], [41], [42], [43], [44], [45], [46], [47], [48],

[49], [50], [51], [52], [53], which is used on

Eulerian grids, does not imply an exact selection of

the contact surface and allows the use of cells

containing mixtures of substances. With this

approach, one has to construct an artificial non-

physical equation of state for the mixture.

Accordingly, it is necessary to construct a solution

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

174

Volume 18, 2023

for the Riemann problem for the Godunov-type

schemes or special algorithms for determining

fluxes and contact parameters for other schemes. In

the most complicated versions, [39], [43], [44], a

multi-component mixture with dynamic equilibrium

is assumed in the cell, with possible sliding of the

mixture components inside the cell — a multi-

velocity continuum. This mixture must exhibit the

properties of a deformable solid, passing

continuously into a fluid or gas. The contact surface

is not explicitly defined in this approach. This

approach is convenient for solving 3D problems, but

it has a significant numerical viscosity and does not

have the required accuracy when describing

complex contact phenomena such as friction,

separation, cavitation, etc. In, [50], [51], similar

approaches that have made it possible to solve

several complex problems of the dynamic

interaction of structural elements with gases in a 3D

setting were developed. The multimesh approach

proposed in, [25], [26], [27], [28], which is close in

meaning to the gas dynamic numerical method

"Chimera", [54], has no drawbacks associated with

the difficulties in identifying and tracking contact

surfaces and describing complex equations of state

for dissimilar materials. This approach uses three

types of computational meshes for each body and

will be detailed below. The approach does not

require complex three-dimensional mesh generators,

it is enough to define the surfaces of bodies with

STL files, which significantly speeds up the process

of preparing data for calculation.

2 Governing Equations

The closed system of equations describing the

deformation of a continuous medium in the

approximation of a hypoelastic model of a

compressible elastoplastic body in the Cartesian

coordinate system has the following form, [21],

[28]:

 

,,

0

i

ti

x

u





(1)

 

,,

+ 0

j

i i j ij

tx

u u u

  



(2)

 

,,

0

j

t j i ij x

e eu u



  

(3)

ij

DS

+ S 2

Dt t ij

e





(4)

= (p, )

  

, (5)

Notation: t - time,

i

x

- spatial coordinates,

i

u

-

components of the velocity vector along the axes

i

x

respectively,



- density,

)2/( iiuue 



- total

energy per unit volume,



- internal energy per unit

mass, given by the equation of state (5),

ij



-

which is represented as a ball and deviator parts,

deviator of strain rate tensor

iiijijij pSp



3

1

, 

,

ij

e

- deviator of strain

rate tensor

)(2/1 ,

3

1,, ijjiijijkkijij uuwheree 



. The

symbol

Dt

D

denotes the Yaumann derivative,

which takes into account the rotation of the stress

tensor in Euler variables.

ikjkjkikktij SSuS









k

ij

,

ij

x

S

Dt

DS

, where

)(

2

1

,, ijjiij uu 



,



- is the material shear

modulus, and the index after the comma denotes

differentiation concerning the corresponding

variable (hypoelastic model). As a criterion for the

transition from the elastic to the plastic state, the

von Mises yield condition

2

3

2

Sijij SS





under

uniaxial tension is used, where

S



- the yield stress

under uniaxial tension. The parameter



must

remain positive during plastic deformation under the

condition of fluidity

2

3

S

ijij SS







. The plastic flow is

described by keeping the deviator on the surface of

fluidity, [55]. The system of equations (1) – (5) is

closed by equations of state with the corresponding

parameters. In the absence of shear stresses, the

system (1) – (5) goes over to the Euler equations for

the motion of a compressible fluid or gas.

3 Second order Godunov Modification

for CFD in Euler Variables on

Compact Stencil

Assume all hydrodynamics parameters (primitive

variables) to be linearly distributed in space

between the cell centers for considered time steps.

For example, for cell number i time step

tn

in X

direction, vector function

U

with scalar

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

175

Volume 18, 2023

components

321 ,,,, uuup



has the following

distribution, Figure 1. In this section, index “I”

defines the number of cells.

x

if

x

xx

x

U

x

U

x

UxU iii

ii

i



 



11

1

1 )(

)()(

x

if

x

xx

x

U

x

U

x

UxU iii

ii

i1

1

1 )(

)()(

)()( 









Interpolating these primitive variables (for

example, for cell number “ i “ with boundaries

xi2/1

and

xi2/1

) into points



ii xx ,

1

- values

)(),(1



ii xUxU

and into points





1

,ii xx

- values

)(),(1





ii xUxU

respectively.

Riemann problem is solved like in the original

Godunov method (exact Riemann solution ) but for

the cell's boundary

xi2/1

between values

)( and )( 1



ii xUxU

and for the boundary

2/1i

x

between values

)( and )( 1





ii xUxU

respectively.

For the original Godunov method, these Riemann

problems are solved for the boundary

xi2/1

between

)( and )( 1ii xUxU 

and for the boundary

xi2/1

between values

)( and )( 1ii xUxU

respectively. This is the main difference between

this scheme and the original Godunov method. The

coordinates of the points



ii xx ,

1





1

,ii xx

are not

known yet.

Fig. 1: Linear distribution of primitive parameters

between cell centers

Flux calculation and the corrector step

(integration) are the same as in the original

Godunov method.

So now the integrated values at the time step

1n

t

are functions of primitive variables interpolated

to these points



ii xx ,

1





1

,ii xx

. For linearized Euler

equations, it is possible to solve the Riemann

problem explicitly ( see APPENDIX 1 ) and obtain

hydrodynamics parameters at the new time step as

the explicit functions of coordinates of these points

for interpolation of the primitive variables.

Let's try to use the appropriate selection of these

points to optimize our numerical scheme (for

example, to exclude or introduce the required

viscosity terms). Let's analyze the one-dimensional

linearized Euler equations (one-dimensional are

considered only for simplicity, the analysis of three-

dimensional cases is in, [56]):

0*

2

*

0*

1

*

0



ux

c

p

u

p

u

xt

x

xt



,

(6)

where

00

0,, u

p



are points of linearization,

),( 0000



pcc 

sonic velocity, in our case as this

point of linearization our integrating cell

ii

iu

p,,



can be taken. Using formulas for the solution of the

modified Riemann problem from APPENDIX 1

with expressions for interpolated values and for

second order derivative of velocity in time from

APPENDIX 2 for Godunov scheme formulas for

these linearized Euler equations in APPENDIX 1, it

follows from eq.6 the following modified equation:

),,(*

2

*

0

*

2

****

),,(*

*2

*

2

**

1

**

2

0

2

1

2

0

1

0

2

0

2

1

0

1

0

h

tt

h

o

B

u

A

p

B

uc

A

p

u

p

h

tt

h

o

B

c

p

A

u

B

p

A

uuu

ii

xxxx

x

xt

ii

xx

x

xt

c











(7)

where

nn ttt  1

;

2/12/1   ii

ixx

h

;

2/32/1

1

 ii

ixx

h

;

2/12/3

1

 ii

ixx

h

 

































11

00

111

1/

2

1

iiii

ixxxx

uc

xxxx

h

A

(8)

 

































1100

111

1/

2

1

iiii

ixxxxcu

xxxx

h

B

(9)

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

176

Volume 18, 2023

2/

)(*)(

2

2/

)(

*

)(

*

)(

2

)(*

11

0

1

0

2

0

2

0

2























































iiii

i

ii

i

xxhh

h

c

xx

hh

xx

hh

h

u

uu

t

A

(10)

2/

)

(*)

(

)

(*

)(

2

2/

)(

*)

(

)(

*

)(

2

1

0

1

0

00

2























































ii

i

ii

i

xx

hh

xx

hh

h

c

xx

hh

xx

hh

h

u

t

cu

B

(11)

So (7) are modified equations of (6) for our scheme

and (7) have parameters to vary for obtaining a

better scheme. The difference between (6) and (7) is

the approximation error of our numerical scheme in

simulating equations (6).

In formulas (8-11) coordinates



ii xx ,

1

and





1

,ii xx

are taken relatively to

i

x

; this means

0

i

x

.

Further for simplicity



ii xx ,

1

are used relative to

2/1i

x

(relative to the border of the left cell), and





1

,ii xx

are used relative to

2/1i

x

(relative to the

border of the right cell’s boundary) respectively. For

system (6) to have an approximation error of the

first order of accuracy, it is necessary that:

1;1 11  BA

(12)

From eq. (6-9) and (12) it follows that:



ii xx 1

,





1ii xx

(13)

So choosing the interpolated points according to

(13) is sufficient for first-order accuracy. Now the

difference between (6) and (7) is only in the right

terms (7). For system (6) to have the second order

of accuracy, it is necessary in addition to (12) to

have:

0;0 22  BA

(14)

Then it follows from (10 – 14) that:















































1

11

1100

1

;

4

2

1/

/

82

)(

ii

i

iii

ii

i

xx

h

hhh

hhcu

t

x

(15)











































ii

i

iii

ii

i

xx

h

hhh

hhcu

t

x

1

11

1100

;

4

2

1/

/

82

)(

(16)

Formulas (15-16) determine the coordinates of such

points, interpolation to which provides an

approximation error of the second order of accuracy

of system (6) for an irregular mesh for a given

numerical method. For an irregular grid, it follows

from (15-16) that the resulting scheme is not

conservative, since the coordinates of these points

are functions of the length of the integrating cell and

its left and right neighbors. But for the commonly

used uniform mesh or arithmetic progression mesh,

where

  iiii hhhh 11

, it is conservative

and formulas (15-16) take the following form:

























1

00

1 ;

42

)(

iii xx

cu

t

x

; (17)























iii xx

c

ut

x1

0

0 ;

42

)(

(18)

It is obvious from (17-18) that for an irregular mesh

there is a displacement to a larger cell.

Thus, the second-order approximation error

algorithm for the 1D linearized Euler equation is the

following:

1. definition of coordinates of points



ii xx ,

1

and





1

,ii xx

using formulas (15-16) or 17-18;

2. definition of

)( and )( 1



ii xUxU

and

)( and )( 1





ii xUxU

respectively (Figure 1);

3. solving the Riemann problems between

)( and )( 1



ii xUxU

for boundary

2/1i

x

and

between

)( and )( 1





ii xUxU

for boundary

2/1i

x

respectively;

4. flux calculation and corrector step (integration) as

in the original Godunov method.

Formulas (17) and (18) for obtaining the

coordinates of interpolation points are very physical

and have an obvious geometric interpretation. In

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

177

Volume 18, 2023

Figure 2 the coordinates of points



ii xx and

1

for

the cell’s boundary

2/1i

x

are presented. The zones

of influence on the solution of the Riemann problem

for this boundary of this cell at

2/t

are restricted

by these interpolation points. The total length of this

zone is

tc 

0

. All possible variants are presented in

Figure 2. It is quite obvious that for this scheme

there is no difference in the Riemann problem for

the case of Lagrangian or Eulerian variables - there

is always a Lagrangian choice of the solution, there

is no analysis of the plane of characteristics. There

is also no difference between subsonic and

supersonic cases. The algorithm is obviously

extended for the case of moving meshes – the

interpolation coordinates are calculated like the

zone of influence of the moving boundary or the

coordinates of interpolation are calculated for the

position of the boundary

2/1i

x

at time step level

2/t

. In other words, the coordinates of the zone of

influence with the length

tc 

0

are calculated for

the particle arriving at the boundary of the cell

2/1i

x

at the time step level

2/t

.

The “predictor” step of this scheme is

practically the characteristics method for the time

step level

2/t

for the interpolated values. For

acoustic equations for Courant number 1 for regular

meshes for 1D case, it is the exact characteristics

method for the predictor step and the scheme

elaborated in this case is equivalent to the original

Godunov method and has the same second-order

accuracy. In general case choosing the interpolation

points in the cell’s centers changes this modification

into the exact original Godunov method. Moving

points of interpolation from values calculated by

(17) or (18) one can introduce viscosity into this

scheme, the value of this numerical viscosity can be

regulated, by choosing the points of interpolation.

The main advantage of this approach is that 3 cell

stencil is enough for the second order accuracy,

another very important thing – switching to the

Riemann solution without interpolation (to the

original Godunov method) provides the monotony

of this modification on the same 3 cell stencil.

The algorithm for the exact second-order

accuracy for linearized Euler equations in 2D, [57],

and 3D, [56], cases is the following:

1. for each boundary of the cell it is assumed that

the normal direction to the center of the surface of

the boundary is X direction;

2. the coordinates of these interpolation points are

the same as for the 1D case;

3. primitive variables at the cell's centers are to be

corrected with the influence of the tangential

gradients at

2/t

;

4. interpolation of these corrected primitive

variables to appropriate points;

5. Riemann problem between these interpolated

primitive variables;

6. calculation of the velocities tangential to the

boundary surface (Riemann solver for tangential

velocities) with the higher order accuracy at

2/t

like in the Lax-Wendroff scheme;

7. flux calculation and integration as in the original

Godunov's method.

Fig. 2: Zone of influence on the Riemann solution at

the cell boundary at

2/t

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

178

Volume 18, 2023

The scheme obtained has second-order accuracy

and it is not monotone. It is possible to introduce the

numerical viscosity into this scheme moving the

points of interpolation from values calculated by

(17) and (18) (zero viscosity points) to the

appropriate cell’s centers – to the viscosity of the

original Godunov scheme, which is enough for

monotony. Analysis of the numerical results shows

that the origin of nonmonotonicity is the big rate of

the gradients. So the nonmonotonicity problems are

coupled with the regions where the behavior of the

solution mainly depends on the second-order

derivatives. Physically these regions are at the top

or the bottom of acoustic or shock waves, density

discontinuities, and rarefaction waves. From the

analysis of numerical results, obtained for this

scheme, for the monotony of the numerical solution

in such a problem-generating region it is enough to

use the original Godunov predictor step for the one

or two cell's boundaries. The main problem is how

to detect such boundaries. It was suggested to use

two parabolic splines, constructed from the

primitive variables (pressure and density) in the

vicinity of the analyzed cell’s boundary, using three

points. The two points of the centers of the cells, are

used for the Riemann problem, and the nearest point

to the left for the left spline and the nearest point to

the right for the right spline. Such parabolic

functions represent more precisely the behavior of

numerical results. For example, these spline

functions can be not monotone in these appropriate

intervals when the numerical ones are still

monotone. So these spline functions can be used

even to predict the nonmonotone behavior of the

numerical one. If the coordinate of the maximum or

minimum of this spline function is not in

appropriate interval it is possible to use Riemann

solver for zero viscosity. The advantage of such

spline function analysis is that it does not depend on

the absolute values, only the analysis of the

coordinate of maximum or minimum is enough for

the choice of interpolation points. In the regions of

constant solutions, the behavior of these spline

functions does not have meaning for switching,

because in these regions it does not matter to use a

first or second-order Riemann problem solver. In

the regions of smooth solutions or even in the

regions of big gradients these spline functions are

monotone, but when there are regions of big rates of

gradients, these spline functions are not monotone.

4 Second-order Godunov Modification

for CSD in Euler variables on

Compact Stencil

In, [21], it was showed that for numerical modeling

of dynamic elastoplastic equations with a second

order of approximation for schemes of the

"predictor-corrector" type at the "predictor" stage, it

is sufficient to solve elastic equations

0

t



with a

second order of approximation using linearized

equations (1) - (5 ). In this case, the plastic behavior

of the medium is taken into account at the

“corrector” stage after the integration of the

linearized equations and is reduced to the “landing”

of deviators on the yield surface, [21]. Let us denote







)(

1

; )(

2

p

f

p

cs









, as a result, after

linearizing the energy equation in the system of

equations (1) - (5) and excluding the terms

containing derivatives with respect to

2

x

and

3

x

(splitting in spatial variables is performed), we

obtain a one-dimensional system of equations:

0 1,11,1,  uu

t



(19)

0*/1*/1 1,

111,1,11,1  Spuuu t



(20)

0*/1+ 1,

121,21,2  Suuu t



(21)

0*/1+ 1,

131,31,3  Suuu t



(22)

0

)(

1,11,3131,212

1,111

2

,





puufSufS

ufScpt



(23)

0

11,31,2

4/3- 1,1113121,1,11  SuuSuSuS t



(24)

0 -3/2 1,2211,2121,1,22  SuuSuS t



(25)

0 -3/2 1,3311,3131,1,33  SuuSuS t



(26)

00.5 )0.5(2 1,1211,3231,22211,12  SuuSuSSS t



(27)

0)0.5(2 0.5 1,1311,333111,223,13  SuuSSuSS t



(28)

00.5 -0.5 - 1,2311,3121,213,23  SuuSuSS t

(29)

In matrix form, equations (19-29) are:

0







x

U

A

t

U

where

)],,,,,,,,,,[( 231312332211321 SSSSSSpuuuU





A

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

179

Volume 18, 2023

11213

1

1133

23

123

1122

113

112

11312

11

2

1

0000005.05.000

000000

2

5.000

0000005.0

2

00

0000000

3

2

0

0000000

3

2

0

000000

3

4

0

0000000

0

1

00000000

00

1

0000000

00000

11

000

000000000

uSS

u

SS

S

uS

SS

uS

uSS

ufSfS

fS

с

u

















In the case of constant coefficients of matrix A

(linearized case), the system is hyperbolic and can

be written in the form of 11 transport equations (in

invariant form), [23], [24].

11,,1 ,0 







l

x

R

c

t

Rl

l

,

Where

l

R

are Riemann invariants constant for the

corresponding characteristic velocity

l

c

. In the

formulas below and Figure 3, the index "1" for the

velocity

1

u

and coordinate

1

х

is omitted.

; ; ;

; ; ; ;

111098765

4321

ucccccucuc

ucucaucauc

zz

yy







Where





11

22 3/4 fS

сa



,

2

23

2

3322

11

y

)(*25.0

5.0

4/3







SSS

S







2

23

2

3322

11 )(*25.0

5.0

4/3







SSS

S

z







In the plane,

tх,

1

the trajectories of

discontinuities (characteristics) are depicted by rays

emanating from the point x = x0, and divide the half-

plane t> 0 into 8 zones.

Fig. 3: Riemann solution zones for an elastic case, 1

order.

The ratios on these characteristics

11,1, lRl

have

the following form:

 

13

2

23

2

3

2

22

2312

2

13

12

2

23

2

3

2

22

2313

2

3

2

12

3

2

23

2

3

2

2312

2

13

2

23

2

3

2

2313

2

3

2

12

1111

25.0))((

5.0*)(

)1(

25.0))((

5.0*)(

)1(

/25.0))((

/5.0*)(

)1(

/25.0))((

/5.0*)(

)1(

S

Sbaba

SSbS

f

S

Sbaba

SSbaS

f

u

Sbaba

SSbaS

fa

u

Sbaba

SSbaS

fa

SpuaR























































































 

13

2

23

2

3

2

22

2312

2

13

12

2

23

2

3

2

22

2313

2

3

2

12

3

2

23

2

3

2

22

2312

2

13

2

23

2

3

2

22

2313

2

3

2

12

1112

25.0))((

5.0*)(

)1(

25.0))((

5.0*)(

)1(

25.0))((

5.0*)(

)1(

25.0))((

5.0*)(

)1(

S

Sbaba

SSbaS

f

S

Sbaba

SSbaS

f

u

Sbaba

SSbS

fa

u

Sbaba

SSbaS

fa

SpuaR























































































 

1312323 ][][ SCSuCuR yy 



 

1312324 ][][ SCSuCuR yy 



where









)(5.0

;

)(5.0

3311

2

3

2211

2

SS

b

SS

b







     

1312325 SSCuuCR zz 



WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

180

Volume 18, 2023

where

10 ;

)(

4

11

2

1

2

3322

2

23







 С

SS

S

С

   

1312326 SСSuuСRzz 



 

13

2

23

2

3

2

122313

2

12

2

23

2

3

2

1323

2

312

11

2

7

25.0

)5.0)(1(

25.0

)5.0)(1(

S

Sbb

SSSbf

S

Sbb

SSbSf

SpaR















































13

2

23

2

3

2

122313

2

12

2

23

2

3

2

1323

2

312

118

25.0

5.0

25.0

5.0

3/4

S

Sbb

SSSb

S

Sbb

SSbS

SR



























































13

2

23

2

3

2

1223

12

2

23

2

3

2

312

229

25.0

5.0

25.0

3/2

S

Sbb

SS

S

Sbb

bS

SR





















































13

2

23

2

3

2

13

2

12

2

23

2

3

2

1323

3310

25.0

5.0

3/2

S

Sbb

Sb

S

Sbb

SS

SR

























































2313

2

23

2

3

2

132312

2

12

2

23

2

3

2

1223

2

313

11

25.0

25.05.0

25.0

25.05.0

SS

Sbb

SSSb

S

Sbb

SSbS

R















































Here, the values in square brackets are the

linearization coefficients determined from the

averaged parameters in the cells. For a scheme of

the first order of accuracy, for the zone where the

solution is sought (where the corresponding face of

the cell is located), the corresponding invariants are

determined and the primitive parameters necessary

for calculating the fluxes are determined from them.

For a second-order accurate scheme, in contrast to

the equations of gas dynamics, the Riemann

invariants, [23], [24] obtained from primitive

parameters are interpolated from the centers of the

cells and are to be corrected with the influence of

the tangential gradients at

2/t

. The coordinates

of the interpolation points are determined as the

boundaries of the areas of influence of the

corresponding invariants on the position of the

boundary face at the time

2/t

, where

w

edge

velocity (indicated by a dotted line in Figure 4) is as

follows:

11,..,1 ,

2

)(

2

)( 1







n

t

wc

xx

xn

ii

n

,

where

; ; ;

; ; ; ;

111098765

4321

ucccccucuc

ucucaucauc

zz

yy







Let us denote the interpolated invariants by the

index “ '', respectively

11..,,1 , )( 1

1

1



 



nxx

xx

RR

RR in

ii

i

n

i

n

i

n

m

n

.

Fig. 4: The coordinates of interpolation points of

Riemann invariants for a second-order accurate

scheme

The obtained values of the Riemann invariants

are used to determine the primitive parameters and

"stream" values. The stage of numerical integration

of the equations ("corrector" step) remains

unchanged (as for the first-order scheme by, [21]).

The solution according to the second-order scheme

will experience dispersive oscillations at the

discontinuities. To ensure monotonicity, it is

necessary to construct the solution in a hybrid

manner - in the regions of smoothness according to

the relations providing the second order of

approximation, and on discontinuities - according to

the relations of the first order scheme. For

elastoplastic flows, in contrast to the problems of

gas dynamics, where it is necessary to analyze the

pressure and density fields, to obtain monotonic

solutions, it is sufficient to analyze the normal stress

m

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

181

Volume 18, 2023

field and use the obtained viscosity as a relative one

to determine the interpolation coordinates of the

remaining invariants, [23], [24].

5 Riemann Solver for FSI

The solution of the contact problem between fluid

and solid is realized at the "predictor" stage of the

Godunov scheme (at the stage of solving the

Riemann problem). From the side of the deformable

body, 8 invariants arrive at the boundary and 3

boundary conditions are used. On the fluid side,

three relationships are used for nonlinear

compression and rarefaction waves. To improve the

accuracy in the domain of smooth solutions,

extrapolation of parameters from the boundary and

boundary cells is used. An iterative algorithm for

obtaining a joint solution, in the general case, looks

like this:

1. for the fluid parameters, boundary conditions of

the movable "rigid wall" type are realized, the

normal velocity of which is equal to the velocity in

the adjoining cell of the deformable body; as a

result, we obtain the pressure at the "fluid -

deformable body" interface.

2. for obtaining the parameters of the boundary cell

of a deformable body, the following equations are

used: Riemann invariants number 1,3,5,7,8,9,10,11

are taken from the center of the boundary cell,

possibly with interpolation from inside the body to

improve accuracy, instead of the Riemann invariant

number 2 the normal stress equal to pressure with

the opposite sign, obtained in item 1, is used, and

instead of Riemann invariants number 4 and 6 - zero

tangential stresses or stresses taking into account

friction are used.

3. The new normal boundary speed defined in item

2 is used cyclically in item 1 as the normal velocity

of the rigid wall. The process continues until the

convergence of this normal velocity. As a rule, 3-4

iterations are sufficient until convergence with a

relative accuracy of 0.01.

6 Multimesh SIM for 3D Nonmoving

Euler Mesh

This SIM approach is multimesh and uses three

types of computational meshes, [25]. The first type

of meshes is Lagrangian meshes in the form of STL

files that define and accompany the deforming

surfaces of bodies. Fixed regular Euler meshes with

cubic cells are used inside homogeneous regions.

The third type of grid is auxiliary local movable

Euler-Lagrangian meshes associated with the

surfaces of bodies. In general, the algorithm for

calculating the contact interaction of fluids and

solids in nonmoving Euler consists of a sequence of

the following steps:

1) Fluids and solids are specified in the form of

surfaces from sets of triangles with the required

precision - in the form of STL files containing

external normals and coordinates of triangle

vertices. In Figure 5 such objects are marked in red

and blue.

Fig. 5: STL surfaces of touching objects and local

meshes

2) Each computational domain (fluid, solid)

with curved boundaries is enclosed in a bordering

rectangular parallelepiped and covered with a

Cartesian mesh. Figure 6 shows a cross-section of

such a computational domain (black solid curve)

with a bordering parallelepiped. We get four types

of cells for the computational domain, the first type

- cells cut by triangles of the surface or boundary

cells, are colored green, the second type - is cells

outside the surface, the third and fourth types - are

cells inside the surface, for which there are enough

(brown) or not enough to integrate (marked with

black dots) of a difference stencil of whole cells

located inside the surface.

Fig. 6: Cross-section of the 3D computational

domain

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

182

Volume 18, 2023

3) On each triangle of the surface, inside the

volume along the normal for this surface, an

auxiliary local Cartesian grid 3 * 3 * 3 is built , in

Figure 5 these meshes are yellow. The cell sizes of

this local three-dimensional mesh are taken close to

the cell sizes of the inner mesh. In case of contact

with another subdomain or boundary conditions that

require additional parameters, the local mesh is

symmetrically completed in this subarea from the

plane of the triangle, in Figure 5 this completed

mesh is marked in green. Figure 7 shows this

auxiliary mesh for one triangle.

Fig. 7: Local mesh for one triangle of STL mesh

Figure 8 shows one triangle and the adjacent

central cells of the auxiliary mesh, the centers of the

edges of these cells coincide with the center of the

triangle.

Fig. 8. Triangle and two adjacent cells of the local

mesh.

Parameter values in the generated local mesh

are determined by interpolation of parameters from

the previous local mesh and parameters from the

main mesh. This stencil is sufficient to integrate the

central cells (in Figure 7 it is highlighted in dark and

in Figure 9 these cells are shown separately and

marked with crosses) adjacent to the Lagrangian

contact surface, with increased accuracy according

to the modified Godunov scheme. For the central

cells from Figure 8, the Riemann problem at the

contact boundaries of the media is solved. Its

solution results in velocities and forces at a half-

time level in the center of the corresponding

triangle. We move this contact boundary with

normal speed and get a local mesh on a new time

layer. Figure 9 shows these center cells before (left

side) and after moving the contact boundary (right

side). We carry out the standard integration of the

parameters of these central cells in movable meshes

(ALE integration). Since the motion of the local

grid is one-dimensional, the volume and surface

integrals over the cells are calculated exactly in this

case.

Fig. 9: Central cells before and after integration

with the motion of the contact boundary in the

normal direction (ALE integration).

4) Using the velocities in the center of each

triangle obtained from solving the Riemann problem

in step 3), we calculate the velocities at the vertices

of the triangles in the STL file with weights

proportional to the areas of the triangles. With these

speeds, we move the vertices and get the position of

the surface on a new time layer (the new position of

the STL file).

5) We rebuild the bordering parallelepiped with

the possible addition or reduction of layers of cells

by the new position of the STL file. The cross-

section of this surface with a bordering

parallelepiped is shown in Figure 10. The red curve

corresponds to the position of the computational

domain on the new time layer. To the cells of the

fourth type, which remained in the computational

domain on a new temporary layer (marked with

crosses on a white background), add non-integrated

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

183

Volume 18, 2023

cells captured during surface movement (marked

with crosses on a green background).

Fig. 10: New STL surface position with additional

non-integrated cells.

6) In the cells of the fourth type (marked with

crosses on a white and green background), we

interpolate the parameters from the integrated cells

of the third type and the integrated cells of the local

meshes, thus completing the calculations.

7 Modeling the Interaction of Ice

Impactors with a Clamped

Titanium Plate

The problem statement is shown in Figure 11. An

ice ball on the left and a cylinder on the right hit the

titanium plate. It is required to describe in a coupled

formulation the processes of interaction of ice

projectiles with an elastoplastic deformable plate.

Plate dimensions from the origin (embedment

plane): length along the X-axis 150 mm, width

along the Y-axis 50 mm, height along the Z-axis 5

mm, the plate is rigidly fixed in the plane X = 0,

material - VT6 titanium, [58], with mechanical

characteristics: density 4.51 g/cm3, the bulk

compression modulus is 112 GPa, the shear

modulus is 41GPa, yield strength 1.45 GPa. At the

moment of impact, the coordinates of the point of

contact between the ball and the axis of symmetry

of the cylinder are X = 134 mm, Y = 25 mm, and Z

= 5 mm. An ice ball has a diameter of 28 mm, a

cylinder - diameter of 17.1 mm, height of 50 mm;

weights are the same as 10.35 g, mechanical

characteristics of ice: density 0.9g / cm3, the bulk

compression modulus is 9.196 GPa, the shear

modulus is 3.5263 GPa, yield point 5.2 MPa. The

initial vertical velocity of the impacters in both

cases is 350 m/s. At the outer boundaries of the

metal barrier and ice impacters, the conditions are

fulfilled as at the “free boundary” with a given

pressure p = 0.1 MPa. The dimensions of the cells

along the ice and the obstacle were taken 0.5 mm,

which required about six hundred thousand cells of

the main grid. In the process of impact, the hailston

and the plate undergo significant displacements and

deformations, and the destruction of the hailston

occurs. Figure 12 and Figure 13 also show the

positions of the bodies and the computational mesh

in the longitudinal plane of symmetry along the ball

and cylinder at an instant of 800 μs. Figure 14

shows the velocities of the plate versus time on the

axis of symmetry of the sphere, marked in red, and

of the cylinder, marked in green, from the beginning

of the interaction to 120 μs. Figure 15 shows the

same parameters from the beginning of the

interaction to 800 μs. The difference in the initial

moment is due to the larger initial interaction area

for a cylindrical object. In the process of interaction,

the strikers spread out over the plate, completely

losing their shape, and starting from the moment of

120 μs, the influence of the initial geometry ceases

to affect the movement of the plate.

Fig. 11: Initial position of impactor and plate

Fig. 12: Shapes of objects 800 μs after impact

Fig. 13: Computational mesh for spherical and

cylindrical impactors, a cross-section along the X

axis, 800 μs after impact

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

184

Volume 18, 2023

Fig. 14: The vertical velocity of the plate at the

point of impact versus time, red for the ball, green

for the cylinder (0-120 μs)

Fig. 15: The vertical velocity of the plate at the

point of impact versus time, red for the ball, green

for the cylinder (0-800 μs)

8 Interaction of Elastic-Plastic Bodies

with the Product of Detonation

In Figure 16 the statement of the problem is given.

Charge of TNT-RDX TR 8530 of the spherical

shape, radius 5 cm, the source of initial detonation

was set in a region of radius less than 0.2 cm. The

cross-section passes through the center of the charge

and the centers of the cubes. In the charge and the

air, the main mesh was used with the side of the

cubic cell 0.15 cm. For the explosive, the

parameters of the equation of state of the JWL type

are taken from, [13], as A=708.6056 GPa,

B=13.16457 GPa, C=1.058238 GPa, r1=4.94,

r2=1.35,



=0.28. The same state equation is used

for air, air is assumed like the product of detonation

for Euler simulation without tracking the interface

"detonation product" – "air". Three steel cubes with

sides 1 cm are marked with red color, the density is

7.8 g/cm3, the bulk compression modulus is 175

GPa, the shear modulus is 80.77 GPa, the hardening

modulus is 240 MPa, and the yield strength is 340

MPa, weight 7.8 grams. The step of the main grid in

cubes is 0.7 mm. In Figure 16 the grids for steel and

TNT are given. In Figure 17 the density distribution

is presented at the moment of 55 μs. The cubes are

strongly and irreversibly deformed, the streams of

detonation products move much faster, and gas jets

are formed around the cubes. In Figure 18 the shape

of the cubes, respectively, at the initial moment and

for 13 and 55 μs. By 13 μs the fragments practically

acquired a residual form by slightly changing their

position. Figure 19 shows the velocity versus time

on the surface of the cube, accelerated in the vertical

direction, the lower curve corresponds to the center

of the upper surface (far from the charge),

respectively, the upper curve - to the center of the

lower surface of the cube (near to the charge).

Figure 20 shows the velocity of the center of mass

versus time for the cube, accelerated in the vertical

direction, marked by "1" and for the cube,

accelerated in a 45-degree direction, marked by "2".

Dashed lines indicate these dependencies, calculated

in an elastic formulation. Taking into account the

plastic deformation of the bodies increases the

maximum velocity of the body by more than 20

percent. The maximum velocity of a body strongly

depends on its shape and position relative to the

surface of the charge.

Fig. 16: Statement of the problem. Mesh for steel,

TNT, and air

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

185

Volume 18, 2023

Fig. 17: Density distribution at 55 μs

Fig. 18: Shape of the cubes, respectively, at the

initial moment, at 13 and 55 μs

Fig. 19: Velocity versus time on the bottom and top

of the cube

Fig. 20: Velocity of the center of mass versus time,

initial positions in vertical and 45-degrees. Dashed

lines show the calculations of the elastic setting.

9 Perforation of Aluminum Plate with

Ogive Nose Steel Rod at Oblique

Impact

Figure 21 shows the formulation of the problem of

punching with a steel impactor, velocity of 400 m /

s, with an aluminum plate, [59], at an angle of 30

degrees. The surfaces (STL files) of the impactor

and target are shown in detail in Figure 22. Steel,

the density is 7.85 g/cm3, the bulk compression

modulus is 175 GPa, the shear modulus is 80.77

GPa, the yield strength is 3.4 GPa and the hardening

modulus is 2.4 GPa. Aluminum, the density is 2.71

g/cm3, the bulk compression modulus is 67.64 GPa,

the shear modulus is 26 GPa, the ideal plasticity,

and the yield strength is 0.262 GPa. The dimensions

of the main mesh in both bodies are 0.01 cm In the

experiment, [59], the plate is 55x55x2.63 cm; in the

calculations, the slab of smaller dimensions is

10x10x2.63 cm with free boundaries, rests along the

perimeter on a rigid frame (Figure 23) with a width

of 1 cm. In Figure 24 and Figure 25, there are

computational grids at 200 and 540 μs, respectively.

Figure 26 shows the punched plate in the direction

of the initial velocity vector of the impactor at the

moment of 540 μs. The Figure 27 shows the shapes

of the striker at the moment of 40 µs, with a

deflection with a tendency to ricochet, and at the

moment of departure 540 µs with the opposite

deflection, which was noted in experiments, [59].

The numerical values of the residual velocity of the

impactor vary in the range from 195 to 205 m / s

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

186

Volume 18, 2023

due to elastic vibrations, which is close to the

measured average values - 200 m / s, [59].

Fig. 21: Oblique impact, rod velocity 400 m/s. STL

surfaces

Fig. 22: Part of STL surfaces

Fig. 23: Rigid frame supporting the aluminum plate

Fig. 24: Computational mesh at 200 μs

Fig. 25: Computational mesh at 540 μs

Fig. 26: Perforated aluminum plate at 540 μs

Fig. 27: Steel rod at 40 and 540 μs (STL surfaces)

10 Conclusions

The multimesh numerical technique developed and

described in this article for solving three-

dimensional problems of propagation of wave

processes in fluids and solids and the interaction of

the fluid with elastoplastic deformable bodies at

large displacements and deformations allows

obtaining reliable results with sufficiently high

accuracy. The results of numerical studies are in

good agreement with the known experimental data.

The use of STL files for constructing all types of

meshes and describing the motion of Lagrangian

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

187

Volume 18, 2023

contact and boundary surfaces within the framework

of this multimesh approach, as well as the use of

regular fixed Euler mesh, makes it possible to

dramatically simplify the preparation of initial data

for solving complex problems and to increase the

efficiency and accuracy of calculations by

eliminating distortions and rebuilding of

computational meshes in traditional Euler-

Lagrangian techniques.

Acknowledgments:

This work has been supported by the grants of the

Russian Science Foundation, RSF 22-79-10076.

References:

[1] Godunov S.K.. A finite difference method for

the computation of discontinuous solutions of

the equations of fluid dynamics.

Math.Sbornik, Vol.47, 1959, pp.271–306.

[2] S.K. Godunov, A.V. Zabrodin, and G.P.

Prokopov. A difference scheme for two-

dimensional nonstationary fluid dynamics

problems and flow with a bow shock wave,

Zh. Vychisl. Mat. Mat. Fiz. Vol. 1, 1961,

pp.1020–1050.

[3] V.P. Kolgan. Application of the minimum

derivative principle for constructing finite

difference schemes for the calculation of

discontinuous solutions in fluid dynamics,

Uchen. Zapiski of TsAGI Vol.3, No.6, 1972,

pp.68–77.

[4] B. Van Leer, Towards. The ultimate

conservative difference scheme V: A second-

order sequel to Godunov’s method, J.

Comput. Phys. Vol.32, 1979, pp. 101–136.

[5] P.L. Roe. Approximate Riemann solvers,

parameter vectors and difference schemes,

Journal of Computational Physics, Vol. 43,

1981, pp.357-372.

[6] P. Colella and P.R. Woodward. The

piecewise-parabolic method (PPM) for gas-

dynamics simulations, J. Comput. Phys.Vol.

54, 1984, pp.174–201.

[7] Hill D.J., Pullin D., Ortiz M., Meiron D.. An

Eulerian hybrid WENO centered-difference

solver for elastic–plastic solids, Journal of

Computational Physics. Vol. 229, 2010, pp.

9053–9072.

[8] Trangenstein J.A., Colella P.. A higher-order

Godunov method for modeling finite

deformation in elastic–plastic solids,

Commun. Pure Appl. Math. Vol. 44, 1991, pp.

41–100.

[9] Miller G.H., Colella P.. A high-order

Eulerian Godunov method for elastic-plastic

flow in solids, J. Comput. Phys., Vol. 167,

2001, pp. 131–176.

[10] Gavrilyuk S.L., Favrie N., Saurel R..

Modelling wave dynamics of compressible

elastic materials, J. Comput. Phys. Vol. 227,

2008, pp. 2941–2969.

[11] Favrie N., Gavrilyuk S.L.. Dynamics of shock

waves in elastic–plastic solids, in: ESAIM

Proceedings. 2010, pp. 1–18.

[12] Dumbser M., Peshkov I., Romenski E.,

Zanotti O. High order ADER schemes for a

unified first order hyperbolic formulation of

continuum mechanics: viscous heat-

conducting fluids and elastic solids, J.

Comput. Phys. Vol. 314, 2016, pp. 824–862.

[13] Barton P.T., Drikakis D., Romenski E.,

Titarev V.A. Exact and approximate

solutions of Riemann problems in non-linear

elasticity, J. Comput. Phys. Vol. 228, 2009,

pp. 7046–7068.

[14] Barton P., Romenski E. On computational

modelling of strain-hardening material

dynamics, Commun. Comput. Phys. Vol.11,

2012, pp. 1525–1546.

[15] Menshov I., Mischenko A., Serezkin A.

Numerical modeling of elastoplastic flows by

the Godunov method on moving Eulerian

grids. J. Mathematical Models and Computer

Simulations, Vol. 6, No. 2, 2013, pp. 127-135.

[16] Barton P.T., Drikakis D. An Eulerian method

for multi-component problems in non-linear

elasticity with sliding interfaces, J. Cornput,

Phys. Vol.229, No. 15, 2010, pp. 5518-5540.

[17] Barton P.T., Drikakis D., Romenski E. An

Eulerian finite-volume scheme for large

elastoplastic deformations in solids, Int. J.

Numer. Methods Eng. Vol. 81, 2010, pp. 453-

461.

[18] A.López Ortega, Lombardini M., Pullin D.I.,

Meiron D.I. Numerical simulation of elastic-

plastic solid mechanics using an Eulerian

stretch tensor approach and HLLD Riemann

solver, J. Comput. Phys. Vol. 257, 2014, pp.

414–441.

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

188

Volume 18, 2023

[19] Schoch S., Nordin-Bates K., Nikiforakis N.

An Eulerian algorithm for coupled

simulations of elastoplastic-solids and

condensed-phase explosives, J.Com-put.

Phys. Vol. 252, 2013, pp. 163–194.

[20] Titarev V.A., Romenski E., Toro E.F.

MUSTA-type upwind fluxes for non-linear

elasticity, Int. J. Numer. Methods Eng. Vol.

73, 2008, pp. 897-911.

[21] Kukudzhanov V.N. Decomposition method

for elastoplastic equations Mech. Solids Vol.

39, No. 1, 2004, pp. 73–80.

[22] Kukudzhanov V. N., Coupled models of

elastoplasticity and damage and their

integration, Mechanics of Solids Vol. 41, No.

6, 2006, pp. 83–109.

[23] Abouziarov M., Aiso H., Takahashi T. An

application of conservative scheme to

structure problems, Series from Research

Institute of Mathematics of Kyoto University.

Mathematical Analysis in Fluid and Gas

Dynamics. No. 1353, 2004, pp. 192-201.

[24] Abouziarov М.Х., Aiso H. An application of

retroactive characteristic method to

conservative scheme for structure problems

(elastic-plastic flows), Hyperbolic Problems,

Theories, Numerics, Applications. Tenth

International Conference in Osaka.

Yokohama Publishers, Inc., pp. 223-230.

[25] Abuziarov K.M., Abuziarov M.H., Kochetkov

A.V. 3d fluid-structure interaction problem

solving method in Euler variables based on

the modified Godunov scheme, Materials

Physics and Mechanics. Vol.28, No.1-2,

2016, pp. 1-5.

[26] Abuziarov К.М. The method of

decomposition of gapes in the three-

dimensional dynamics of elastoplastic media.

Problems of Strength and Plasticity. Vol. 82.

No 3, 2020, pp. 5-17.

[27] K.М. Abuziarov, М. Abuziarov, Е.G.

Glazova, А.V. Kochetkov, S.V. Krylov, Е.Е.

Maslov, В.I. Romanov. Numerically

modeling 3d processes of explosive

acceleration of elastoplastic bodies. Problems

of Strength and Plasticity. Vol. 80. No 2,

2018, pp. 255-266.

[28] Abuziarov M, Glazova Е.G., Kochetkov А.V.,

Krylov S.V.. Numerical method for solving

three-dimensional problems of interaction of

high-speed gas jets with elastoplastic barriers.

VANT Series Mathematical modeling of

physical processes. No. 4, 2021, pp. 24-40.

[29] Barton P., Deiterding R., Meiron D., Pullin D.

Eulerian adaptive finite-difference method for

high-velocity impact and penetration

problems, J. Comput. Phys.,Vol. 240, 2013,

pp. 76–99.

[30] Barton P.T., Drikakis D. An Eulerian method

for multi-component problems in non-linear

elasticity with sliding interfaces, J. Cornput,

Phys. Vol. 229, 2010, pp. 5518-5540.

[31] Barton P.T., Drikakis D., Romenski E. An

Eulerian finite-volume scheme for large

elastoplastic deformations in solids, Int. J.

Numer. Methods Eng. Vol. 81, 2010, pp. 453-

463.

[32] Schoch S., Nordin-Bates K., Nikiforakis N.

An Eulerian algorithm for coupled

simulations of elastoplastic solids and

condensed-phase explosives, J.Com-put.

Phys. Vol. 252, 2013, pp. 163–194.

[33] X. Zeng, C. Farhat, A systematic approach for

constructing higher-order immersed boundary

and ghost fluid methods for fluid and fluid-

structure interaction problems, Journal of

Computational Physics, Vol. 231, 2012, pp.

2892–2923.

[34] C. Farhat, J.F.D.R. Gerbeau, and A.Rallu.

FIVER: A finite volume method based on

exact two-phase Riemann problems and

sparse grids for multi-material flows with

large density jumps. Journal of

Computational Physics, Vol. 231, No. 19,

2012, pp. 6360–6379.

[35] Miller G.H., Colella P. A conservative three-

dimensional Eulerian method for coupled

fluid-solid shock capturing, J. Comput. Phys.

Vol. 183, 2002, pp. 26-82.

[36] Barton P.T., Obadia B., Drikakis D. A

conservative level-set based method for

compressible solid/fluid problems on fixed

grids, J.Comput.Phys. Vol. 230, 2011, pp.

7867–7890.

[37] Udaykumar H.S., Tran L., Belk D.M., Vanden

K.J. An Eulerian method for computation of

multimatereal impact with ENO shock-

capturing and sharp interfaces, J.

Comput.Phys. Vol. 186, 2003, pp. 136–177.

[38] Barton P.T. An interface-capturing Godunov

method for the simulation of compressible

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

189

Volume 18, 2023

solid-fluid problems, J. Comput.Phys. Vol.

390, 2019, pp. 25-50.

[39] Wallis Tim, T. Barton Philip, Nikiforakis

Nikolaos A flux-enriched Godunov method

for multi-material problems with interface

slide and void opening, Journal of

Computational Physics. Vol. 442, 2021, pp.

110499-110511.

[40] Michael L., Nikiforakis N. A multi-physics

methodology for the simulation of reactive

flow and elastoplastic structural response,

J.Comput.Phys. Vol.367, 2018, pp.1–27.

[41] Favrie N., Gavrilyuk S.L., Saurel R. Solid-

fluid diffuse interface model in cases of

extreme deformations, J.Comput.Phys. Vol.

228, 2009, pp. 6037-6077.

[42] Favrie N., Gavrilyuk S.L. Diffuse interface

model for compressible fluid - compressible

elastic-plastic solid interaction.

J.Comput.Phys. Vol. 231, 2012, pp. 2695-

2723.

[43] Ndanou S., Favrie N., Gavrilyuk S. Multi-

solid and multi-fluid diffuse interface model:

applications to dynamic fracture and

fragmentation, J. Comput. Phys. Vol. 295,

2015, pp. 523–555.

[44] Michael L., Nikiforakis N. A hybrid

formulation for the numerical simulation of

condensed phase explosives, Journal of

Computational Physics. Vol. 316, 2016, pp.

193–217.

[45] Jackson Haran, Nikiforakis Nikos. A unified

Eulerian framework for multimaterial

continuum mechanics, Journal of

Computational Physics. Vol. 401, 2020, pp.

109022.

[46] Adler M.C., Jain S.S., West J.R., Mani A. and

Lele S.K. Diffuse-interface capturing

methods for compressible multiphase fluid

flows and elastic-plastic deformation in

solids: Part I. Methods, Center for Turbulence

Research, Annual Research Briefs. 2020,

http://dx.doi.org/10.13140/RG.2.2.15129.651

24.

[47] Wallis Tim, T. Barton Philip, Nikiforakis

Nikolaos. A flux-enriched Godunov method

for multi-material problems with interface

slide and void opening, Journal of

Computational Physics. Vol. 442, 2021, pp.

110499.

[48] Yashraj Bhosale, Tejaswin Parthasarathy,

Mattia Gazzol A remeshed vortex method for

mixed rigid/soft body fluid–structure

interaction, Journal of Computational

Physics. Vol. 444, 2021, pp. 110577.

[49] Krayukhin A.A., Stadnik L.N., Yanilkin

Yu.V. Numerical simulation of the motion of

rigid projectiles in elastoplastic media on a

fixed computational grid using the “Egak”

method, VANT Series Mathematical modeling

of physical processes. No.1, 2019, pp. 19-32.

[50] Yanilkin Yu.V. Closure models for the

equations of Lagrangian gas dynamics and

elastoplastics in multicomponent cells Part 1.

Isotropic models. VANT Series Mathematical

modeling of physical processes. No.3, 2017,

pp. 4-21.

[51] Yanilkin Yu.V., Toporova O.O., Kolobanin

V.Yu. Closure models for the equations of

Lagrangian gas dynamics and elastoplastics in

multicomponent cell Part 2. Anisotropic

models. VANT Series Mathematical modeling

of physical processes. No. 3, 2017, pp. 22-38.

[52] Hank S., Favrie N., Massoni Jacques

Modeling hyperelasticity in non-equilibrium

multiphase flows, J. Comput. Phys. Vol. 330,

2017, pp. 65–91.

[53] Saurel R., Pantano C. Diffuse-interface

capturing methods for compressible two-

phase flows, Annu. Rev. Fluid Mech. Vol.50,

2018, pp. 105–130.

[54] Luis Ramírez, Xesús Nogueira, Pablo Ouro,

Fermín Navarrina, Sofiane Khelladi, Ignasi

Colominas. A Higher-Order Chimera Method

for Finite Volume Schemes. Archives of

Computational Methods in Engineering,

Springer Verlag. Vol. 25, No.3, 2017, pp.

691-706.

[55] M.L. Wilkins, Calculation of elastic-plastic

flow, in Methods in Computational Physics,

edited by B.Alder, S.Fernbach, and M.

Rotenbeg, Academic, New York, Vol. 3,

pp.211-263.

[56] Abuziarov M. On the increase of the accuracy

of Godunov's method for solving the

problems of dynamics of gases and liquids,

XIII conference of young scientists of Moscow

Physical-Technical Institute, Vol.2. 1988,

pp.30-37.

[57] Abuziarov M., Bazhenov V.G., Kochetkov

A.V. On the monotonization of the Godunov

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

190

Volume 18, 2023

scheme of the second order of accuracy by

introducing the scheme viscosity, Applied

problems of strength and plasticity. Research

and optimization of structures: 1987. pp. 85-

90.

[58] Merzhievsky L.A., Paletsky. AV.

Calculations for diagrams of dynamic

deformation of metals and alloys. Physical

Mesomechanics, Vol.4, No 3, pp.85-96.

[59] Piekutowski AJ, Forrestal MJ, Poormon KL,

Warren TL. Perforation of aluminum plates

with ogive nose steel rods at normal and

oblique impacts. Int J Impact Eng., Vol.18,

1996, pp. 877-887.

Appendix 1

Godunov's scheme

0

12/12/1

0

2/12/1

0





























x

pp

x

uu

u

t

uu iiiii

i



0

2/12/1

2

00

2/12/1

0





























x

uu

c

x

pp

u

t

pp iiiii

i



Godunov's Riemann's problem

 

;

2

1

21

00

1

2/1 ii

ii

ipp

c

uu

u









 

;

22 1

001

2/1 ii

ii

iuu

cpp

p









 

;

2

1

21

00

1

2/1 







ii

ii

ipp

c

uu

u



 

;

22 1

001

2/1 







ii

ii

iuu

cpp

p



Riemann's problem for this scheme

 

; )()(

2

1

2

)()(

1

00

1

2/1













ii

ii

ixpxp

c

xuxu

u



 

; )()(

22

)()(

1

001

2/1













ii

ii

ixuxu

cxpxp

p



 

; )()(

2

1

2

)()(

4

00

4

2/1 xpxp

c

xuxu

ui

i

i









 

; )()(

22

)()(

1

001

2/1













ii

ii

ixuxu

cxpxp

p



Appendix 2

 

2

-

tt t otu

2

1





t

i

iu

t

uu

 

2

-

tt t otp

2

1





t

i

ip

t

pp

6 sys. of eq. second theof ;0*

2

*

6 sys. of eq.first theof ;0*

1

*

6 sys. of eq. second theof ;0*

2

*

6 sys. of eq.first theof

t

;0*

1

*

0

x

u

c

p

u

p

x

p

uuu

t

u

c

p

u

p

uuu

xx

xxtx

xx

xxtx

xx

xttt

xt

xttt



















 

xxxxtt ucup

u

u2

0

2

0

0 2 



   

t o 2

2

-

2

0

2

0

0





















xxxxt

i

iucup

u

t

u

t

uu



 

xxxxtt pcuucup 2

0

2

0

2

000 2 



    

t o 2

2

-

2

0

2

0

2

000 









xxxxt

i

ipcuucu

t

p

t

pp



 







 





 2/

2/

)()(

)()( 1

1

1ii

ii

ii xh

hh

xpxp

 

2/

)o(h

4*2

)(

2

)()()(

)( 1

1

-

2

i

2

11



























 ii

ii

ixx

ii

ixii

ixh

hh

xp

hh

xpxpxp

xp

 

)(2//2)(

2

1

2/)()()( 2

1111 iiiiiixxiiixii hoxhhhxpxhxpxpxp 













 

)(2//2)(

2

1

2/)()()( 2

1111 iiiiiixxiiixii hoxhhhxuxhxuxuxu 













 

)(2//2)(

2

1

2/)()()( 2

1iiiiiixxiiixii hoxhhhxpxhxpxpxp 



 

 

)(2//2)(

2

1

2/)()()( 2

1iiiiiixxiiixii hoxhhhxuxhxuxuxu 



 

 

)(2//2)(

2

1

2/)()()( 2

1iiiiiixxiiixii hoxhhhxpxhxpxpxp 



 

 

)(2//2)(

2

1

2/)()()( 2

1iiiiiixxiiixii hoxhhhxuxhxuxuxu 



 

 

)(2//2)(

2

1

2/)()()( 2

1111 iiiiiixxiiixii hoxhhhxpxhxpxpxp 













 

)(2//2)(

2

1

2/)()()( 2

1111 iiiiiixxiiixii hoxhhhxuxhxuxuxu 













WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

191

Volume 18, 2023

The authors equally contributed to the present

research, at all stages from the formulation of the

problem to the final findings and solution.

Sources of Funding for Research Presented in a

Scientific Article or Scientific Article Itself

This work has been supported by the grants of the

Russian Science Foundation, RSF 22-79-10076.

Conflict of Interest

The authors have no conflicts of interest to declare.

Creative Commons Attribution License 4.0

(Attribution 4.0 International, CC BY 4.0)

This article is published under the terms of the

Creative Commons Attribution License 4.0

https://creativecommons.org/licenses/by/4.0/deed.en

_US

WSEAS TRANSACTIONS on FLUID MECHANICS

DOI: 10.37394/232013.2023.18.17

M. H. Abuziarov, E. G. Glazova,

A. V. Kochetkov, S. V. Krylov

E-ISSN: 2224-347X

192

Volume 18, 2023

Contribution of Individual Authors to the

Creation of a Scientific Article (Ghostwriting

Policy)