Application of grid-characteristic method for numerical solution of deformable solid mechanics dynamical problems Текст научной статьи по специальности «Медицинские технологии»

UDC 519.6

Application of grid-characteristic method for numerical solution of deformable solid mechanics dynamical problems * I. B. Petrov*

Moscow Institute of Physics and Technology (State University), Moscow, Russian Federation

The grid-characteristic method is a promising numerical method for solving hyperbolic systems of equations, e.g., equations describing elastic and acoustic waves. This method has high precision and allows you to physically correctly simulate wave processes in heterogeneous media. The grid-characteristic method makes it possible to correctly take into account boundary conditions and conditions on surfaces with different physical characteristics. Most fully the advantages of the method are for one-dimensional equations, especially in combination with a fixed difference grid, as in conventional grid-based methods. However, in the multidimensional case using the algorithms of splitting with respect to spatial variables, the author has managed to preserve its positive qualities. The use of the method of Runge - Kutta type, or integro-interpolation method for hyperbolic equations makes it possible to effectively carry out the generalization of methods developed for linear equations, in nonlinear case, in particular, to enforce the difference analogs of the conservation laws, which is important for shock-capturing, for example, discontinuous solutions. Based on the author's variant of the grid-characteristic method were numerically solved several important problems of seismic prospecting, seismic resistance, global seismic studies on Earth and Mars, medical applications, nondestructive testing of railway lines, the simulation of the creation and characteristics of composite materials for the aerospace industry and other areas of practical application. A significant advantage of the constructed method is the preservation of its stability and precision at the strains of the environment. This article presents the results of numerical solution based on the grid-characteristic method to the problem of modeling elastic-plastic deformation in traumatic brain injury.

Key words: hyperbolic type equations, elastic and plastic deformation, grid-characteristic method, finite difference schemes, discontinuous solutions.

Introduction. The numerical solution to nonstationary dynamic problems of elastic and elastic-plastic deformations of isotropic solids in case of two or three spatial coordinates is the object of many works, a fairly detailed overview is given in [1], [2]. In solid mechanics various mathematical rheological models have been developed, among which an important place is occupied by those which are described by systems of equations of hyperbolic type, and, above all, the models of linear and non-linear theory of elasticity and plasticity.

The most important concept for such equations is their characteristic properties. The use of the characteristic records of the original equations when constructing a particular numerical

* The research is done within the frame of the independent R&D. **E-mail: petrov@mipt.ru

methods (i.e., approximation of not the original equations, but the equivalent to them conditions of compatibility along some characteristic directions) allows in the most natural way building computational algorithm on the borders of the region of integration, to a certain extent taking into account solutions dependence domain, etc. Most fully the advantages of the characteristic approaches occur for one-dimensional equations, especially in combination with a fixed difference grid, as in conventional grid-based methods (grid-characteristic methods [4]); however, in the multidimensional case using the algorithms of spatial splitting the variables are able to preserve their positive qualities. The use of the equations of = Runge - Kutta type or integro - interpolation method [6] makes it possible to effectively carry out the generalization of methods developed for linear equations, in nonlinear case, including by ensuring the implementation of the difference analogs of the conservation laws, which is important for shock capturing of discontinuous solutions.

Out of other approaches widely used for numerical solution to such problems, we will note the method of disintegration of discontinuities and various difference schemes of splitting. The choice of the coordinate system used is crucial at considerable strains of bodies. If for small displacements (and strains) the traditional for solid mechanics problems, Lagrangian finite-difference grid does not undergo significant changes then with large deformations in the area of integration heavy distortion on it until disturbances are possible (zero area of differential cells, their "eversion", etc.). Used in such cases methods of regularization of difference grids (reinterpolation on a new grid, etc.) do not always give the expected effect.

On the other hand, the use of a fixed in space Eulerian coordinates leads to serious difficulties in the approximation of boundary conditions. Therefore, it is important to compare calculations performed in various coordinate systems, using the same numerical scheme that is one of the objectives of the present work.

2. Using conventional kinematic equations for a symmetric tensor of velocities of deformations [7] in spatially fixed orthogonal curvilinear coordinate system x and x2, x3

choosing the defining relations in the form

do,. do- u, do- u do,, u3 do,. * . .

= —J + —--- + —--- + —--- = y q„m„emn, ,,J = 1,2,3,

at at Hj dxj H2 dx2 H, dx, „fj^,

we write a closed system of two-dimensional nonstationary equations in the form of

du ~ ôu ~ du

— + Aj — + A„-=/. (1)

dt dXj dx2

Strictly speaking, in the case of finite deformations, instead of substantial derivative

da{j I dt Qne ghouls use the so-called Jaumann derivative for the components of the stress tensor deviator of type

do S do S i

— -——Ift^+W

k'—l

dt dt

1

œa = — J 2

dVj

Ox,-

dv

j

dx,

A

1 3

sj = °J - sJi kk,

3 k=i

which provides a zero rate of change of the stress state of the particle during its rotation as a rigid whole [2]. However, as the calculations show, in the examined task it is almost insignificant.

The symbol

u = {ui ,u2 ,on a12 a22 ,o33 ,T}

means the vector of desired variables, including

the components u and u 2 of the velocity vector V (on the axes x1, x2 respectively) non-zero components atj of the symmetric stress tensor and the temperature T;

{/i> fi->fw> fm fiffr}— vector of right parts with the following components in the curvilinear orthogonal coordinate system x1, x2, x3

f = F +1

(a

11-a22 H1H 2

) dH2 (a

f2 = F2-

(a

22 -a33 H 2H3

_ qiJ11V2 OH

dx,

dx

J11 -a33 H1H3

)

) dd (a

'22 -a11 -> H 2H1

/■■ =

iJ H ¡H 2 dx2

2_ (qiJ 12 - qiJ 21 )

1 =1 1 pc

a11v2 dH1 H1H2 dx2

2H1H2

a 22 H1H2

dH

a 12 H„

a12 H,

dH

dH

H dx H dx

H1H2

dH

dH

1 Ox- 2 dx,

dH

OH,

H dx

H dx

H2H1

2 dH

dH,

2 cx, 1 dx

2

dH

V J--r V -

1 Ox o 2 Ox 7

qij22v2 OH

H H dx

2 \ iJ33

H

3

Vj dH

H j Ox j

3^ V2

2

dH31 4yq

'3

H2 Ox2

pc

OH,

OH 2'

1 Ox0 2 Ox,

a22V1 OH2 H H Ox

a 33 H,

v, dH

3

H Ox

2 sh£ H 2 Ox2

-Q

Here F1, F2 are the components of the vector of mass forces, Ç - internal energy for thermoelastic medium, с — the specific heat of the material, T — temperature, Q —the density of heat sources; Hi, i=1, 2, 3 — Lame coefficients that characterize the selectable orthogonal curvilinear coordinate system, P is the density determined by the equation of state of a solid body, for example

lnp = -(3K)-1 £ a, - 3LeT Po

i=1

P — density of material in the deformable condition, К - the coefficient of volumetric compression. In accordance with the intended here two-dimensionality of strain state, displacement of material points in the direction х3 and is absent, and the matrices

ct13 = a23 = v3 = 0, O / Ox3 = 0 thus have the form:

Ä -± 1 ~ H

Л2--П~ H2

u1 0 -1/p 0 0 0 0

0 u1 0 -1/p 0 0 0

~1llll ~(llll2 +41121)' '2 и 0 0 0 0

1l211 ~(ll212 + q 1221)' 2 01 U1 0 0 0

q2211 (q2211 + q2221) 20 0 U1 0 0

q3311 (q3312 + q3321) 20 0 0 U1 0

~G11f pc - °11/pc 0 0 0 0 u1

U2 0 0 -1/ p 0 0 0

0 U2 0 0 1 /p 0 0

~(llll2 + Il 121 ) ' '2 -q1122 U2 0 0 0 0

-(q1212 + q1221)' '2 -q1222 0 U2 0 0 0

~(q2212 + q2221 ) '2 -q2222 0 0 U2 0 0

q3311 (q3312 + q3321) 20 0 0 U2 0

- °12/pc ^22/pc 0 0 0 0 uv

(2)

For the adopted in the present work the model of Prandtl —Reuss components of a tensor of the fourth rank qiikl elastic-plastic material have the form: yo.A

qijkl ^ij^kl pc + & jl + & jk & il ) к where А, / — the parameters of the Lame, k is the yield stress in shear, Smn are the symbols of

Kronecker, the stress tensor deviator

1 з

Smn = °mn ~ & mn 1 SS

3 S=1

where у = (3А + 2/i)at (at is the coefficient of linear expansion of the material during heating), I is determined from the Mises plasticity J = |0 при S = S1 + й + й + Si2 < , |1 при S > 2k2

Defining relations for I=0 are linear Hooke's law for elastic isotropic materials. When у =0, Q=0 first 6 equations do not included temperature and can be solved independently from the energy equation. It should be noted that thermal effects are usually significant only when external heat is supplied, or for intensive computation. The above mentioned system of equations is a formal combination of known dynamic equations for small thermo-elastic (at I=0) and elastic-plastic (у =0, Q=0), without the energy equation deformations of isotropic materials. These two well-known models were used for the numerical solution is given in the work task. For small deformations, the differences between the Lagrangian and Eulerian description of the motion is insignificant and the system can be have further simplifications; in particular, it is possible to neglect the convective members (leaving out и, Ь on the main diagonals of the

matrices A, A, which was assumed in problems on small elastic and elastic-plastic deformations), i.e. in the matrix A and in this case were replaced by A = A = (u/HjJE, A = A = (v2/H2)EAk = {ak},k = 1,2,i;j = 1,.,7. To account for finite deformations and displacements, we introduce the moving coordinate

system connected with the fixed Eulerian х1, х2, х3 ratios

^Sij Skl

t=t, x1(t,Tjl,Tj2) = Tj1+jc1(t,Tj1,Tj2)dt

0

t

x2(t,1i,l2) = 1i + j x3 = 1-

(3)

in which 8xk / 8t\n=const = ck (t,1,12), ^=1, 2, — the components of the local velocity of the

points in a moving coordinate system with respect to a fixed in space coordinate system x1, x2, x3. The Jacobian matrix of the transformation has the form

't 1 0 0

B = xit xil! xil2 = Ci dxi / dii dxi /di2

x2t x2i2 x2l2 Ci dx2 / 5ii dx2> di2

(4)

is nonexceptional if A = DetB = x} x2 — x1 x2 . Then:

B = I A

A 0 0 h tx x2

(Cixil — Cix2i2) x2l2 —xil2 = in 1ixi

( Cix2ii —C2 xiii) —x2li xii 'h, h2 x 2

A-1 = DeW-1 = rj1xi^2x1 - Vixc^x =1 / A

and thereby it is established the connection between derivatives 81 / 8^ . = 1 and

8xk / 81 = xkm .

With this transformation, the system of equations takes the form:

Ut + A1Ui + A2U%= f Ak = 1E + A1 + 1kx2 A2, k = ^ (5)

moreover, the diagonal member of the matrices Xk = {«*} and Xk = {«*} source matrices are connected by ratios

1 , hxji , 1x2V2

H

H,

2 , l2xi Vi , 12X2 V2

H

H

V Hi

7 „

X2li

A

V H 2

V H 2 f „

xi1i

\

11 "2

and non-diagonal members take the form of:

~1 ~2

1 ~1 , ~2 X2n aij + Xln aij

V H1

— c

112

2li

A—1 = v

A—1 = v

A

~2 ~1

<*,] = + Vi,",, = -—

A

For the used in the work the Lagrangian coordinate system:

c = u / H¡,

c2 = u2/H 2

C$=Dk=0,

the coordinates of the nodal points x^, xn+1 of a uniform coordinate system n difference grid

Vim = mK m = 0,1,...,M, r = /h2, / = 0,1,.., hj = 1/Z

(6)

were determined (and was memorized as vectors u"ml, u"2 ) and ratios approximation, for example the simplest:

\ml 2 TC1ml ,

yn+1 _ vn . _,,«+1 X1ml = X1ml + lc1,

Yn+1 _ vn . n X2ml = X2ml + ic2ml,

and their derivatives dx / d^ were found by numerical differentiation using central differences in the internal nodal points 1 < m < N -1, 1 < l < L -1, for example

(x2n )ml = '

_ (X2 )m+1,l (X2 ) m—1,l x2r)ml = 2 ' (X2%2 2h

^ \ ( X2) m,l+1 — ( X2) m,l—1

, (X2r )ml =

and one-sided differences for the boundary node points. In some cases, and in the internal nodal points were used one-sided differences of form

(X2r )ml =

(X2r2 )ml , (X2n )ml < (X2n )ml

(X2r )ml , (X2n )ml > (X2r )ml

where

, , _(X2)m+1,l —( X2)m—1,^ , .(^l —(X2)m—1,l

(X2rr )ml = 7 , (X2r2 )ml =

h

For moving Eulerian coordinate system

x2 - Xx -(7,Xj)

h

t = t, = x2, ?;2 = ■

X2-(t,xl)-X2 ■ (t,Xj) ^

where X2 ■ (Lxi) and X2 ■ (I.xi) are the equations of the upper and lower limits of integration we

have B 1 =

1 0

—P

Vi X

0 1

-r X

0 0

r X,

(here P = [X2i+t h{x;-t-X2i)\ q = [x^ +rj2(x:_Xi-X^j], ,hxi = (X2+-X^,).

This allows us to directly use the ratios, or to restrict ourselves instead of calculation (and

memorizing) only of values X2+ (t, xx) n X2 (7, Ay), with the help of e.g., approximation

X-2(t"+1,nlm) = X-2(f,r,lm) + rX-t(f\rjlm), X+2 (t"+1,n1m) = X2(t",nm) + tX(t",nm),

2t{ hJ 2t{ lJ H^o) H^o) ' 2t{ hJ 2t{ lJ H^l) H^J)

Then the necessary derivatives drfc //8xf, (либо cxI / ct,cxI let]-) are determined

through X2t±=dX2±(t,v)/8t,X2x±=SX2±(t,v)/d^1, using numerical differentiation of the

proportions.

Finally believing

Ci=0, C2= и 2/Н2,

we get a mixed Eulerian-Lagrangian coordinate system, in which you can restrict yourself to calculation (and memorizing) only value x^h = xn2ml + w"2ml / H^ml, different from zero and derivative units dx2/ , /'=1,2.

All three of the considered system of coordinates leaves the original rectangular region of integration

0<x(0,^1,^2) = V <Лъ 0<x2(0,v,v2) = щ < 1

and introduced in it uniform differential grid (6) stationary in variables tw , however, in variables t, xi, x2 the location of grid nodes for //,„,.//,, each of these coordinate systems is uneven and significantly affect the calculation of problems with large deformations.

When using for approximation of equations system (1) grid-characteristic method, the

design equations in the internal (Vm V), m = 1,2..., M-1,/ = 1,2,...,/-1,

2+1 _ 2 1 —-f2 . yn I yn

Um/ = Um/ + Jm/1 + b 1m/ + b 2m/'

n

\2 U„

nodal points are

the

ratios

h" - T h1ml - T

h" - T h2ml - 1

( "—^n—1 )" ( «^A )

У 1 1 1 'ml ^

" Umt1,l ~ Uml

( a 71a ;a—1 )

" U , -u ,

l—1,m ml

ml

h,

(a21A2a1 )

" Um+1,l~ Uml

ml

A,

(8)

in which the diagonal matrices

At

Ak +\Ak\

> Ak

Ak~ \Ak\

, Ak| - (Kk }

Kk = A }, i = 1,...,7;

— eigenvalues of the matrices Ak, k=1.2, determined from the characteristic equations Det(AkA E)=0; Q = \cok} — nonexceptional matrix, rows of which are linearly independent left

eigenvectors rak of matrices Ak, defined with accuracy to the length of the plurality of linear homogeneous systems of equations (AkT— AkE)ak =0, i=1,..., 7; Q-1 —the inverse k of the matrix; A[ — transposed matrix Ak).

For matrices A|-= {a* } we have A- = vk + y-,i = 1,...,7,

J_

yk — -y7 —

ak +

/ T \1/2

( a2- - )

1/2

y2k=-yk =

ak-

1/2

( ak - 4dk )

1/2

y3k = -y4k = y5 = 0.

ak — ai3a31 + ai4a41 + a24a42 + a25a52, dk — ai3a24 ((a3ia42 a32a41 )) + ai4a25 X ((a4ia32 a42a31 ))

+aoa225 ((ahak32 ~a\2a\i)), k = 1,2,

+

1 1 T12 T13 1 T14 Ti5

1 Ti1 1 ®23

0 0 „1 T32 „1 T34 1

0 0 1 T43 T44 0

0 0 1 T53 „1 T54 0

1 Ti1 1

1 1 T12 -Ti3 -t34 -Ti5

0 0 0 0 00

1 0

00 00

^2

2 T11 1 2 T13 „2 T14 „2 T15 0 0

1 _2 T22 _2 T23 21 2 T25 0 0

0 0 1 _2 T34 _2 T35 0 0

0 0 0 „2 T44 „2 T45 1 0

0 0 0 2 T54 _2 T55 0 1

1 2 T22 -T223 -T224 -T225 0 0

„2 T11 1 -Ti23 -T124 -T125 0 0

Here

i _( yi) -a1

®12 =

®11 =■

ai

a;

i _ a2

®21 — TTTK2 ~T,

( yi)2 -ai1

1 _( y22) -ai2

( yi2)2-ai

2

®12 =

a

ai ai3a31 + ai4a41,

v - aWi

y»

k _ a25©i2

5 =

3 =

i — 1,2,

y,

a\2ßii- a\iß\

(9)

V - aik4^ü + a2k4W»k2

k —1,2,

a32ßii aiißi2 2 _ a52ß2 - a51ß,^2

4 —

4 —

y,

a32ßii a31ßi2 a31ßi2 - ai2ßii k _ a42ßl -a41ß,^

a51ßi2 a52ß41

ßi — ai ;Ö35 + ai ;Ö36 + ai j©^ ßj — j2 + j2 +

®15 —

a2 ß2 a2 ß2 a52ßii a51ß42

j —1,2, i — 3,4,5.

, 1 — | pk | in this case for inverse matrices we have

2

1

3

4

1

5

6

7

Q-1

pf1 pk2 0 0 0 p1k2 pÎ1

p2k1 P22 0 0 0 p2k2 p2k1

pk31 P32 p3k3 0 0 -p3k2 -p3k1

P41 p42 p4k3 0 0 ^4k1

pks1 pk2 p5k3 0 0 -p5k2 -p5k1

A P62 p6k3 1 0 -p62 -p6k1

pn pk2 p7k3 0 1 -p7k2 -p7k1

where

i i (a2)

Pi = P22 =

2

1 _ p31 =

Pi2 =■

ii

»

2Ai

P21 = -

ii

»

21

2Ai

2 A1

1 _ »»L

p33 =

11

2 A1

1 »14

P12 = —

1 »23 P41 = —

2 A1

11

2A1

„1 _ »15»33 ~»113 p42 =

2A1

„1 _ »215»113 - ».5».3 p43 =

A1

„1 _ »214»313 p51 =

2A1

1 _»113»^4 -®114»^3

p52 =

2 A1

„1 _ »114»^3 -®113»^4 p53 =

A11

A2 = 1 -®112»^1-

p,1 =

pÎ2 =

p13 =

»1-2,3 (» - »»L ) - »»-2,4 (»^3 - ®^5®]3 )] (2 A1 ) ,

»/-2,4 (»B - »115»313 ) - »/-2,3 (»1^4 - »115»314 )] (2A1 ) 1, »)-2,3 (»114»^5 - »115»^4 ) - »)-2,4 (»1» - »»3 )] (A1 ) , A1 = »115 (»:>4»313 - »^»L ) + »25 (»113»114 - »114»313 ) + (»U»^ - »^»L )

And, similarly

A! =

»

2A

pu =

2

2A 2

2 2 2 2 2 _ »24»35 »225»34 p31 =

2 A

2 22

„2 _ »25 - »23»35 p 31

2A

p222 =

»

2A

2

A2 = 1 -»,,».

22

2 2 2 2 2 _ »15»34 »»1

p32 =

2 A

2 2 2 2 _ »13 »35 - »15 p 31

2A

2 2 2 2 2 _ »4»25 -»5»24

p33 =

A2

2

2

2

2 2 2 2 2 2 2 2 _ P15P23 ~ P13P25 1 P23P34 ~ P24

Аз A2 ' 2A2 '

2 2 2 2 2 2 2 2 Щ14 - Щ1ЗЩЗА , _ ®1_з®^-®23®1_1

Р52 2А2 ' p53 А2 '

2 I 2 / 2 2 2 \ 2 / 2 2 2 \ I/* 2\ 1

p,1 = L®i- 2,3 (Р24 _ ®25®34 ) " Р2,4 (Р25 " Р23Р35 ) j (2 А ) ,

2 I 2 / 2 22\ 2 / 2 22 \ I i 1 \ 1 Pi2 = L®i-2,4 (®15 - P13P35 ) " "-2,5 ("14 - Р13Р34 )j j2 А 2 ) ,

Р* = \_tf-2,4 (Р1зР25 - ""З ) - Pf-2,5 (""З - Р13Р24 )] ( А1 ) ,

i = 6,7,

2 2 { 2 2 2 2 \ 2 / 2 2 2 2 \ / 2 2 2 2 \

А j = p3 (ю24 ®35 - ю25 о34) + p (®15®34 - о14 ®35) + (о14 ю25 - ®15®24 ).

When constructing calculation formulas on the borders of the rectangular (in coordinates ) region of integration will only consider the upper (щ = 0) and lower (щ = 1) borders, meaning that the remaining boundaries (щ = 0, щ) are often a plane (or axis) of symmetry or periodicity of the decision or chosen so that during the time t < tj the disturbance did not reach these boundaries. The generalization for more complicated conditions on the boundaries щ = 0,щ present no major difficulties and are similarly described below. Scalar multiplying by the eigenvectors (pf)"m}, we obtain ratios

К )>rn+1=B2=K L (u"ml+T/n):

T

V

approximating with the first order of accuracy the conditions of compatibility along the lines of intersection of the characteristic surfaces of the system and the plane coordinates щ =щт (with

equations dq2 = l2dt)

Q2ut + Л2ю;Ч?3 = ш2 (f - ^V^ i = 1,...,7.

As it is known, the number of boundary conditions for hyperbolic system of equations of the type is determined by the number of negative (positive) eigenvalues of the matrix А1 on the upper (respectively, lower) boundary of the region of integration. In the examined tasks at the upper boundary щ = 1 occurs X < Л < 0 at the lower boundary щ = 0, Л2 > Д2 > 0 respectively,

and; consequently, each of these boundaries requires the setting of two boundary conditions, which we write in the form (f)i(t,n1,u1,...,u7} = 0, i = 1,2 при n = 0, (pi (t,n;,u;,...,u7) = 0, i = 1,6,7 при n = 1,

and it is necessary that Det^ * 0, DetQ * 0, where respectively

Q_ = IIoo^^CJo^Oc»2...W^H , Q+= Цсо^.с^о^Ц .

+

ml "ml V 'ml1....... .......

№ km+«mi), I=10

Here юг- = {дф /ды},...,дф /du7}, i = 1,2,6,7, а i = 1,...,7, - the eigenvectors of the matrices А2. For the considered tasks and the boundary conditions were chosen semilinear and after approximation had the form of

ф = Шm ) Ugt ( m ) = 0,

i = 1,2 при 7 = 0, i = 6,7 при щ= 1.

Attracting the differential ratio at /=3,.., 7 for the border 7 =0 and /=1,.., 5 for border 7 =1, we get all the necessary formulas for points belonging to these boundaries. For the calculation of points belonging to the boundaries 7 =0, 7 = 7, we used calculation formulas with additional "rays" 7 = , (m=-1) and 7 = 7 +h (m=M+1) for which the components of the sought vector u are determined according to the data inside the region of integration with appropriate symmetry or periodicity of the decision or extrapolation (depending on the type of boundary).

Below there are some results of calculations illustrating the capabilities of the considered difference scheme. The data were nondimensionalized in the following way: the linear dimensions were related to some characteristic value х0, the density - to the initial density of the material p0, the components of the stress tensor - to some value a0, temperature T - to the

characteristic value T0; the other dimensionless quantities were defined using the dimension (index р) thus

= vP

V

Po

-1/2

k = 1,2, t = tp

_,p (oo /Po )1/2

xn

2 Xp

-, X =—, p =

On

pp On

J kP Kp k = —, K =-, c =

Oo O O

PoTo - ~y P To, a = aPT0, Q = ■

^ Y = ■

Qp x0

o

o

O

(Oo /Po )

1/2

H-0' HO, (11)

The examples of numerical modeling of various complex dynamic problems are presented in [8-19].

The results of the numerical solution of some urgent problems. Traumatic brain injuries (the task was set by the specialists of the Institute named after Sklifosovsky and Ambulance Central clinical military hospital named after Burdenko)

Fig. 1

b

а

c

Working together with the neurosurgeons we formulated three mechanical-mathematical models of human head. The simplest of these is a two-component model (Fig. 1a), in which bone tissue and brain are described by homogeneous isotropic materials with average mechanical properties; more complex models take into account the presence of ventricle (Fig. 1b) and two membranes (Fig. 1c). The external load is set as the collision of skull-brain system with a completely fixed rigid barrier with a given initial velocity (1-3 m/s). It is known that a significant contribution to the formation of a wave pattern of disturbance propagation in an elastic medium makes heterogeneity. Therefore, a two-component model was further developed in the form of a three-component model, consisting of bones and brain and ventricle. The dural membranes have an inhibiting effect on the movement of brain inside skull. Therefore, subsequently, it was decided to introduce a model of the falx cerebri, with a vertical membrane separating the hemispheres in the parietal region. Fig. 2 a, b, c shows the corresponding computational grids (3 a - quadrangular, 3b, 3c triangular) for these models, which were based on numerical calculations. Rheological properties of biomaterials were also varied. Thus, the rheology of the medulla varied from lianopoulos to viscoelastic. The behavior of bone material was simulated by lianopoulos continuous medium with average properties of lamellar and trabecular bone. Modelling the interaction of skull - brain is challenging due to the fact that brain has several different mechanical properties of the membranes with cavities, including fluid-filled. This work addresses the selection of the contact gap with conditions that ranged from full adhesion to slip with the possibility of detachment.

The distribution of the shear stress in the formulation of various conditions is shown in Fig. 3, 5.

ab c

Fig. 2

Model of traumatic brain sections: a - two-component model, b - model of ventricles, c - model with ventricles and membrane (falx cerebri).

a b c

Fig. 3. The distribution of shear stresses a - complete adhesion; b - slide; c - two-component model.

a b

Fig. 4. a - CT of a patient with brain contusion of severe degree, left impact; b - the corresponding distribution of maximum shear stress.

a = -90c

a = -30c

Fig. 5. The dependence of the distribution of the maximum shear stresses (0 ^ 0.04 ATM) on the impact direction (dark regions correspond to larger shears) a - two-component model; b - model with ventricles;

c - model with membrane

Fig. 3 shows the overall characteristics of the mechanical effects on the brain at a side impact, obtained using different models with free-slip on the border of the skull - brain. The most dangerous are concentrations of maximum tensile (positive) and shear stress. The use of conditions for the complete adhesion on the border of the skull - brain leads to concentration of the shear stress along the contact boundary on the lateral surfaces, while the sliding contact completely removes them. Accounting for the presence of ventricles almost has no effect on the distribution of the regions of maximum compression and extension, but significantly affects the distribution of shear loads. The presence of membrane is more essential to localize areas of compression-tension in side impacts. Fig. 4a shows an example of a CT image of a patient with brain contusion of severe degree with left impact from an accident (data provided by the Main military hospital named after N. N. Burdenko). The arrows on the CT scan indicate the lesion of the brain substance. Fig. 4b shows the corresponding distribution of maximum shear stress. When reducing the area of impact (similar to impact with a pointed object) there is a region of concentration of the shear stress, which coincides with the center of hematoma. Fig. 5 shows the distribution fields of shear stress in human brain depending on the angle of application of the

shock load ^a = -90° ^ 90° j for the three considered models of human head, which gives the

neurosurgeon the information about the localization and the volume of region of brain damage in traumatic brain injury.

Direct problems of seismic prospecting

a b

c

d

Fig. 6

Fig. 6 shows the isosurfaces of the modulus of the velocity for ten-layer geological environment at four points in time (see Fig. 6a, 6b, 6c, 6d for 500 ms, 600 ms, 700 ms, 800 ms, respectively) effects of seismic loading at the upper boundary. It is seen that the used in the paper numerical scheme gives well the complex wave pattern resulting from the interaction of seismic waves with geological layered rock. On the calculated seismogram the reflected from all contact boundaries waves are clearly visible. The obtained calculations show the possibility of using a full numerical simulation of complex wave processes in layered geological media to obtain data and study of the structure of rocks.

Fig. 7

Let's consider the problem of reflection of plane Ricker pulse, the direction of propagation of which is 14 degrees with vertical (this corresponds to the displacement of the point of the explosion along the surface to 0.5 km) from the reservoir 5 vertical filled cracks. The pattern of reflected waves at various distances between the cracks is examined. The q value will change from 0.5 to 4.0 (q = 0.5, 1.0, 1.5, 2.0, 3.0, 4.0). Fig. 7 shows the corresponding velocity

fields. A comparison of the results of calculations of the geological problems on two-dimensional triangular grid on a three-dimensional cubic grid is done. This comparison shows that the discrepancy of the results obtained is not more than 30%. Thus, the two-dimensional calculation gives not only a qualitatively correct result, but also a good quantitative assessment. Therefore, in the study of three dimensional dynamic processes in geological environments it is possible to use two-dimensional calculations. This can significantly speed up computation time and reduce storage costs and data processing.

Impact task

a b c d

Fig. 8 (a, b, c, d) presents the results of calculations (velocity fields) for the impacts on five-layer barriers with gaps (spaced obstacles) and without gaps between layers.

Similar designs are used as protective ones (e.g., body armor). The lower (in the figures on the right) calculated layer corresponds to the protected from hitting human body. As it can be seen from the calculations, after hitting the upper boundary of the five-layer construction, the effect of elastic-plastic compression waves leads to a collision of layers from top to bottom, after which the compression waves go to a protected environment. By simulating attacks on the multilayer barrier we obtained wave pattern corresponding to the distribution of the secondary waves generated by the collisions of obstacles and directed against the impact. In addition, the layers in such structures face then go away from each other, depending on the sign of the normal stress (negative or positive) on the contact boundaries. These processes (collision and detachment) initiate the emergence of secondary waves of compression-stretching in the design and the protected body.

Fig. 10.

References.

1. Kukudzhanov, V.N. Chislennoe reshenie neodnomernykh zadach rasprostraneniya voln napryazheniy v tverdykh telakh. [Numerical solution of multidimensional problems of propagation of stress waves in solid bodies.] Soobshch. po prikl. matem., iss. 6, Moscow, VTS AN USSR, 1976 (in Russian).

2. Novatskiy, V. K.Volnovye zadachi teoriiplastichnosti. [Wave problems of the theory of plasticity.] Moscow, Mir, 1978 (in Russian).

3. Magomedov, K.M., Kholodov, A.S. O postroenii raznostnykh skhem dlya uravneniy giperbolicheskogo tipa na osnove kharakteristicheskikh sootnosheniy. [On construction of difference schemes for equations of hyperbolic type on the basis of the characteristic ratios.] Zh. vychisl. matem. i matem. fiz., 1969, vol. 9, no. 2, pp. 373-386 (in Russian).

4. Kholodov, A.S. Opostroenii raznostnykh skhem spolozhitel'noy approksimatsiey dlya uravneniy giperbolicheskogo tipa. [On construction of difference schemes with positive approximation for equations of hyperbolic type.] Zh. vychisl. matem. i matem. fiz., 1978, vol. 18, no. 6, pp. 1476-1492 (in Russian).

5. Kholodov, A.S. O postroenii raznostnykh skhem povyshennogo poryadka tochnosti dlya uravneniy giperbolicheskogo tipa. [On construction of difference schemes of increased order of accuracy for equations of hyperbolic type.] Zh. vychisl. matem. i matem. fiz., 1980, vol. 20, no. 6, pp. 1601-1620 (in Russian).

6. Tikhonov, A.I., Samarskiy, A.A. Ob odnorodnykh raznostnykh skhemakh. [On homogeneous difference schemes.] — Zh. vychisl. matem. i matem. fiz., 1961, vol. 1, no. 1, pp. 5-63 (in Russian).

7. Sedov, L.I. Mekhanika sploshnykh sred. [Continuum mechanics.] Moscow, Nauka, vol. I, II, 1973 (in Russian).

8. Rozhdestvensky, V.L., Yanenko, N. N. Sistemy kvazilineynykh uravneniy. [Systems of quasilinear equations.] Moscow, Nauka, 1978 (in Russian).

9. Chislennoe reshenie mnogomernykh zadach gazovoy dinamiki. [Numerical solution to multidimensional problems of gas dynamics.] Ed. Godunov, S.K., Moscow, Nauka, 1976 (in Rusian).

10. Redi, J. Deystvie moshchnogo lazernogo izlucheniya. [The action of powerful laser radiation.] Moscow, Mir, 1974 (in Russian).

11. Bullough, R., Gilman, I. 1. Elastic explosions in solid caused by radiation. J. Appl. Phys., 1966, vol. 37, no. 6, pp. 2283-2287.

12. Maynchen, D., Sak, S. Metod rascheta "Tenzor". [Calculation method "Tensor".] Vychisl. metody v gidrodinamike, Moscow, Mir, 1967, pp. 185-211 (in Russian).

13. Favorskaya, A.V., Petrov, I.B. O volnovykh otklikakh ot neftesoderzhashchikh rezervuarov v shel'fovykh zone Arktiki. [On the wave responses from petroliferous resevoirs in the Arctic shelf zone.] DAN, 2016, vol. 466, no. 6, pp. 722-725 (in Russian).

14. Kvasov, I.E., Levyant, V.B., Petrov, I.B. Chislennoe issledovanie volnovykh protsessov v poristoy srede s ispol'zovaniem setochno-kharakteristicheskogo metoda. [Numerical study of wave processes in a porous medium using grid-characteristic method.] ZHVM i MF, 2016, vol. 56, no. 9, pp. 1645-1656 (in Russian).

15. Astanin, A.V., Dashkevich, A.D., Petrov, I.B., Petrov, M.N., Utyuzhnikov, S.V., Khokhlov N.I. Modelirovanie vliyaniya golovnoy udarnoy volny chelyabinskogo meteorite na poverkhnost' Zemli. [Modeling of the influence of bow shock waves of the Chelyabinsk

meteorite on the Earth's surface.] Matematicheskoe modelirovanie, 2016, vol. 28, no. 8, pp. 3345 (in Russian).

16. Sannikov, A.V., Miryaha, V.A., Petrov, I.V. Chislennoe modelirovanie eksperimentov po issledovaniyu prochnostnykh kharakteristik l'da. [Numerical modeling experiments on the study of strength characteristics of ice.] Cb. Polyarnaya mekhanika. Materialy 3-ey mezhdunarodnoy nauchnoy konferentsii 27-30 sentyabrya 2016g. [Col. Polar mechanics. The materials of the 3rd international scientific conference 27-30 September 2016.] Vladivostok, pp. 141-150 (in Russian).

17. Petrov, I., Vasyukov, A., Ermakov, A., Favorskaya, A. Numerical modeling of nondestructive testing of composities. Science Processing Computer Science 96 (2016), pp. 930-938.

18. Vassilevski, Y.V., Beklemysheva, K.A., Grigoriev, G.K., Salametova, A.O., Vasyukov, A.V. Transcranial ultrasound of concept. Russian journal of Numerical Analysis and Mathematical Modeling, vol. 31, no. 5, pp. 239-328, 2016.

19. Levyant, V.B., Kvasov, I.E., Petrov, I.B. Otsenka vozmozhnosti obnaruzheniya kartirovaniya zon razvitiya mezotreshchin v plastakh pri ispol 'zovanii obmennykh rasseyannykh voln. [Evaluation of the detection capability of the mapping areas for the development of mezocracks in layers using converted scattered waves.] Tekhnologii seysmorazvedki, 2016, no. 1, pp. 14-30 (in Russian).

Autor:

Petrov, Igor B., Correspondent Member of RAS, Dr.Sci. (Phys. -Math.), Professor, Moscow Institute of Physics and Technology (State University), Moscow, Russian Federation

УДК 519.6

Применение сеточно-характеристического метода для численного решения динамических задач механики деформируемых твердых сред*

И. Б. Петров **

Московский физико-технический институт

Сеточно-характеристический метод - перспективный численный метод решения гиперболических систем уравнений, например, уравнений, описывающих упругие и акустические волны. Этот метод обладает высокой точностью и позволяет физически правильно моделировать волновые процессы в гетерогенных средах. Сеточно-характеристический метод позволяет корректно учитывать граничные условия, а также условия на поверхностях раздела сред с различными физическими характеристиками. Наиболее полно преимущества метода проявляются для одномерных уравнений, особенно в сочетании с фиксированной разностной сеткой, как в обычных сеточных методах. Однако, и в многомерном случае при использовании алгоритмов расщепления по пространственным переменным, удалось сохранить его положительные качества. Использование для гиперболических уравнений методов типа Рунге - Кутта, либо интегро-интерполяционного метода позволяет эффективно проводить обобщение методов, развитых для линейных уравнений, на нелинейный случай, в том числе, для обеспечения выполнения разностных аналогов законов сохранения, что важно при сквозном счете, например, разрывных решений. На основе разработанного автором варианта сеточно-характеристического метода был численно решен ряд важных задач сейсморазведки, сейсмостойкости, глобальной сейсмики на Земле и Марсе, медицинских приложений, неразрушающего контроля железнодорожных путей, моделирования процессов создания и характеристик композитных материалов для аэрокосмической отрасли и в других областях практического применения. Существенным преимуществом построенного метода является сохранение его устойчивости и точности при значительных деформациях среды. В данной статье представлены результаты численного решения на основе сеточно-характеристического метода задачи о моделировании упруго-пластической деформации при черепно-мозговых травмах.

Ключевые слова: уравнения гиперболического типа, упругие и пластические деформации, сеточно-характеристический метод, разностные схемы, разрывные решения.

Автор:

Петров Игорь Борисович, член-корреспондент РАН, доктор физико-математических наук, профессор, заведующий кафедрой «Вычислительная математика и информатика» Московского физико-технического института (государственного университета), г. Москва, РФ.

* Работа выполнена в рамках инициативной НИР

** E-mail: pctrov@mipt.ru