Numerical Schemes of Higher Approximation Orders for Dynamic Problems of Elastoviscoplastic Media Текст научной статьи по специальности «Математика»

Journal of Siberian Federal University. Mathematics & Physics 2024, 17(1), 8—17

EDN: BXNNJQ УДК 519.63

Numerical Schemes of Higher Approximation Orders for Dynamic Problems of Elastoviscoplastic Media

Vasily I. Golubev*

Moscow Institute of Physics and Technology Dolgoprudny, Moscow Region, Russian Federation Institute of Computer Aided Design of the RAS Moscow, Russian Federation

Ilia S. Nikitint

Institute of Computer Aided Design of the RAS Moscow, Russian Federation

Xin Mi*

Moscow Institute of Physics and Technology Dolgoprudny, Moscow Region, Russian Federation

Received 31.10.2023, received in revised form 17.11.2023, accepted 24.12.2023 Abstract. For a stable numerical solution of the constitutive system of an elastoviscoplastic model of a continuous medium with the von Mises yield condition and hardening, an explicit-implicit second-order scheme was proposed. It includes explicit approximation of the equations of motion and implicit approximation of the constitutive relations containing a small relaxation time parameter in the denominator of the non-linear free term. To match the approximation orders of the explicit elastic and implicit corrective steps, an implicit second-order approximation was constructed for isotropic elastoviscoplastic medium with hardening model. The obtained solutions with the second-order implicit approximation of the stress deviators of the elastoviscoplastic system of equations allow limiting case when relaxation time tends to zero. Correction formulas were obtained in this case, and they can be interpreted as regularizers of numerical solutions for elastoplastic systems with hardening.

Keywords: numerical simulation, elastoviscoplastic media, semi-linear hyperbolic systems, explicit-implicit schemes of higher orders.

Citation: V.I. Golubev, I.S. Nikitin, Xin Mi, Numerical Schemes of Higher Approximation Orders for Dynamic Problems of Elastoviscoplastic Media, J. Sib. Fed. Univ. Math. Phys., 2024, 17(1), 8-17. EDN: BXNNJQ.

An elastoviscoplastic (EVP) medium is a medium that exhibits viscous properties beyond the elastic limit in non-stationary loading processes, i.e., it has a dynamic yield strength that depends on the strain rate. The governing equations of the model are derived from the additive representation of the deviator part of the strain rate tensor as the sum of elastic and viscoplastic components [1,2]. The volumetric viscoplastic component is assumed to be zero. The differential part of the non-stationary system of equations for the elastoviscoplastic media coincides with the system of equations for the dynamic elasticity theory [2,3,4]. However, equations for the

*w.golubev@mail.ru, golubev.vi@mipt.ru https://orcid.org/0000-0003-3113-7299 ti_nikitin@list.ru https://orcid.org/0000-0003-3499-6910 ^ misin@phystech.edu © Siberian Federal University. All rights reserved

deviatoric stress components beyond the yield point include strongly non-linear free terms with a characteristic relaxation time t in the denominator.

In non-stationary processes when characteristic time is much bigger than t elastoviscoplastic media behave like elastoplastic (EP) media [4,5]. When t ^ 0 EVP systems of equations transform into systems of the Prandtl-Reuss type [5,6]. This fact was proved with the use of model formulations based on variational inequalities [7, 8].

For essentially non-stationary processes when characteristic time is less than t effects of highspeed hardening and sharp increase in the dynamic yield strength reveal themselves [4, 6]. To describe these effects it is required to use the full EVP system of equations.

Differential part of the EVP system of equations coincides with the equations of the dynamic elasticity. Then this part is a system of hyperbolic equations, and it can be reduced to a divergence form. The form allows one to build conservative numerical methods for solving initial-boundary value problems.

At the same time, the EP systems of equations, to which the theories of plastic flow are reduced, can not be represented in a divergence form. To describe strong discontinuities one has to construct generalized solutions [8].

On the one hand, formulation of essentially non-stationary problems of the inelastic deformation in the form of a EVP system of equations reflects the physics of dynamic processes. On the other hand, it ensures the regularization of the EP system of equations that does not have divergence form.

The stable integration of the constitutive relations between stresses and strains in the EVP system of equations by an explicit scheme for small relaxation times requires a stronger constraint on the time step than the usual Courant condition. This additional limitation can be eliminated with the help of implicit schemes for calculating "rigid" constitutive relations. Note that explicit-implicit schemes for solving constitutive equations for EVP media described below do not require numerical solution of algebraic systems of equation since all implicit relations can be solved analytically. Each time step is done explicitly with the time step defined by the Courant condition. If it is possible to make limit transition t ^ 0 in the difference approximation for the EVP system then the obtained numerical solution can be considered as the solution of the corresponding elastoplastic problem. Further this approach is systematically used.

1. Isotropic model of elastoviscoplastic continuum media with hardening

In Cartesian rectangular coordinate system, the system of equations for an elastoviscoplastic (EVP) medium with hardening has the form [2]

dt

(1)

(2)

(3)

e = (V ® v + v ®V) /2,

e' = e — (e : I) I/3, 3 = a — al, a = (a : I) /3,

M f(

yjs : s

as(k)

- 1

/t,

(5)

(6) (7)

(8)

where v is the velocity vector, a is the stress tensor, s is the deviator part of the stress tensor, a is the mean stress, e is the strain rate tensor, e' is the deviator part of the strain rate tensor, s : s is the second invariant of the deviator part of the stress tensor, F(x) is the non-linear viscosity function, describing the rate of hardening, F > 0, F(0) = 0, (F) = FH(x), H(x) is the Heaviside function, t is the characteristic relaxation time of components of the deviator part of the stress tensor on the yield surface, p is the medium density, A and n are the Lame constants. Here additional notations are introduced: e'vp is the deviator part of the viscoplastic strain rate tensor, as is the variable yield strength depending on the hardening parameter k, a0 is the initial value of the yield strength (before the hardening process), f (k) is the hardening function such that f (0) = 1.

Taking into account the expression for e'vp, the evolutionary equation for hardening takes the form

dk dt

H

as(k)

1

/t.

(9)

2. Implicit 2nd order approximation scheme

for the elastoviscoplastic medium with hardening

2.1. General solution for an arbitrary hardening function

Consider linear equations for the velocity components and for the mean stress from system (1)-(9). Any explicit scheme of the required 2nd order can be constructed for these equations without difficulties [9]. It is assumed further that such an approximation is carried out. Then taking into account the initial and boundary conditions, values of velocities and the mean stress on the next time layer are determined.

Consider now the equation for the deviator part of the stress tensor from EVP system with hardening (1)-(9). Let us construct an implicit 2nd order approximation of this equation. The right part of the equation is approximated as a half-sum of terms on the next and current time layers:

sn+l _ sn

At M I

t

= M (

sn+1

e'n+1 + e'n) —

\/sn+1

sn+1

«

^/sn+1

sn+1

aof(kn+1)

1

+

V

sn : sn

Mfi—0»,

kn+1 _ kr'

1

t

V/s"+T7s"+T

1

+< A •s"—i

/2.

(10)

(11)

At r\r \ ac f (kn+1) -J/ ' \ ^ \aof (kn) Indices n +1 and n mark the values on the next and current time layers, At is the time step.

e

n

s

This non-linear system of equations with respect to s"+1 can be written as

^n+l = sn+1-

( sn+1 (F(V^+^T^ - \\ _s^ I fVtñJZ - (12)

^ f {kn+i) J J (kn) j/j'

:<f ( - ■)>+(F (

- -kn=¿ Hf^-i)HF[fM--i ■ -

Here new notations are introduced s"+1 = sn+1/a0, sn = sn/a0, srn+1 = sn+1/v0, where

srn+1 = = sn + ¡ (e'n+1 + e'n) At is the deviator part of the stress tensor after the "elastic"

time step, S = and ¡32 = -r— are dimensionless parameters of the system of equations,

At ¡ At

S/^i = &0/ (2¡). In this construction, except for specially stated cases, restrictions on the values of S and ¡i2 are not set.

Taking into account that velocity values on the next time layer are known, ~srn+1 and values s on the n-th time layer are also known.

The non-linear system for the unknown values of sn+1 and hardening parameter kn+1 on the next time layer can be written as follows

+ (f/f tkn+-) - !))■ +

F fórTF/f (kn) - i))

_sn : sn

(14)

Ssn+\

kn+1 = kn + mfñ+Sf - i)) + (f№^V (15)

Performing convolutions of equations (14), (15) sequentially with s"+1, sn and s^1, the system of equations for convolutions X = v/Sn+T7Sn+T, Y = , Z = is

obtained, where values of T = V^Ts^ S = y/SrTSn*1, E = ^s^1 : s^1 are known. Then the following system is obtained

¿X + (F (X/f (hn+1) - 1)) = SS, (16)

^kn+1 - (F (X/f (kn+1) - 1)) = fokn, (17)

where £n+1 = sn+1 - — (F(T/f (kn - ^, kn = kn + — (F (T/f (kn) - 1)), S = VSn+T7sn+T, __T S P2

T = V sn : sn.

It is easy to see that intermediate value of £n+1 and its convolution S are calculated on the basis of results of the "elastic" time step.

Then components of the deviator part of the stress tensor and hardening parameter on the next time layer can be found from Eqs. (14)-(17) as follows

SSn+1 sn+1 ^n+l. = _^_ = £_X (18)

s S + F(X/f(X"+1)-1)) S X (18)

+

kn+1 = S(S - X)/ß2 + kn. (19)

Thus, to completely determine components of sn+1 and hardening parameter kn+1 Eqs. (16), (17) should be solved, and the result should be substituted into Eqs. (18), (19).

Now, it is also possible to consider a limiting case S ^ 0 in the formula for the intermediate value sn+1 which contains a small parameter in the denominator:

Sn+1 = sn+i - T {F (T/f k) - T)). (20)

To do this Eq. (16) that is valid on the n + 1-time layer is rewritten for the n-th time layer:

ST + {F (T/f (kn) - 1)) = S\/sn : sn. (21)

From here

(F (T/f (kn) - 1)) = S^Vsn : sn /sn : sn^ . (22)

Taking into account that S/$2 = a0/ (2^), it is possible to remove the singularity when S ^ 0 and obtain

-sn I \

sn+1 = sn+1 - , (V sn : sn -VSrTSr), (23)

e Vsn : sn \ /

'sn : sn

kn+1 = kn + g - X^, kn = kn + 2°(^Sn : Sn - Vs?1 : sn^ . (24)

This technique was proposed in [10,11], and it was widely used in computational practice (see, for example, [12-15]). However, it was interpreted as a first-order correction based on physical splitting of the elastoplastic process into elastic step and bringing the deviator part of the stress tensor to the yield circle.

2.2. Linear viscosity function and linear hardening

In order to solve non-linear algebraic system (16), (17) and obtain explicit expressions for the unknown quantities s?+1 and kn+1, it is necessary to specify the viscosity function F (x) and the hardening function f (k). To obtain an analytical solution in a system with hardening linear viscosity function F(x) = x is accepted. In this case, values of the dimensionless parameters ¡32 ^ 0, S ^ 0 are not bounded. The hardening function is often taken in the form f (k) = 1 + ak, where parameter a > 0 regulates the mode and rate of change of the yield strength. For a > 0 we have an increase in the yield strength (hardening), for a = 0 we have the ideal viscoplasticity.

Under these assumptions Eqs. (16)-(19) have the form

SX + - t) =SS (25)

1 + akn+1

ß*kn+1 -( Y^U - Y) = ß2kn (26)

X, kn+1 = s(S - X^/ß2 + kn. (27)

sn+1

sn+l. = f_ X kn+1

S

Excluding kn+1 from this system, we obtain the equation for X:

X = (l + S (S - X)) (l + aS(S - X + akn^j

(28)

Let us introduce a new variable Z = S — X and notation y = a/ft2. Then equation for Z takes the form

ys2z2 + (1 + S (1 + Y + akn^ Z — (s — 1 — akn^ = 0. (29)

The exact solution of this equation that satisfies condition Z > 0 is

\/(l + S (l+ y + afc™))2 +4yS2 (S — 1 — a kn) — (l + S (l + y + afc"))

Z = ^------—----. (30)

Hence,

\/Î1 + S (1 + y + a kA) 2 +4yS2 (S — 1 — a — (1 + S (1 + y + a

X = S — ^---—---,-------—. (31)

2ys2

Relations for the deviator part of the stress tensor and hardening parameter on the next time layer take the form

Sn+1 = X

S

/ 1 2 \

v(1 + S (1 + Y + a kn)) +4yS2 (s — 1 — a — (1 + S (1 + y + a k"))

S 2YS2

(32)

/

kn+i = A k = ß2

+ kn.

^ \J(l + S (l + Y + a kn)) 2 +47A2 (S - 1 - aJk^j - (1 + S (1 + 7 + a kn)) ^

27s2

V )

+

(33)

These relations can be simplified if we assume that hardening coefficient a is small, or more precisely, parameter y ^ 1-

Then relations for the invariant, deviator part of the stress tensor and hardening parameter on the next time layer take a compact form

1 + a kn + SS (1 + y + a kn)

X «-^-, (34)

1 + S i 1 + y + akn)

Sn+1 1 + akn + SS (1 + y + a kn) s^1 « ^--V-_J_, (35)

S 1 + S (l + 7 + a kn j

kn+i ^ O0 s- 1 - ^ ) +kn. (36)

1 + S (1 + Y + akn j

These relations allow one to consider a limiting case S ^ 0 (pure elastoplastic material):

sn+1

S^1 - 1 + a kn

sk

(1 + a , (37)

kn+1 ~ 2J (S - 1 - akn) + kn ■ (38)

Condition of the EVP mode with hardening is S - 1 - alkn ^ 0.

3. Numerical simulation results

Using the described explicit-implicit scheme, the wave propagation in an isotropic elas-toviscoplastic medium initiated by instantly applied load was simulated for various values of the relaxation time (parameter S). The medium has the following parameters: density p = 2500 kg/m3 P-wave velocity Cp = -y/(A + 2^)/p = 4500 m/s, S-wave velocity Cs is varied, yield strength as = 112500 Pa. The computational domain was a parallelepiped with dimensions of 50 x 50 x 10000 m. Computational grid was introduced with space step of 5 m. To correctly simulate the process of loading a half-space the condition of zero spatial derivative of the solution was used on the side faces of the computational domain. For the stability of the explicit grid-characteristic scheme [16,17,18,19] used for solution of the elastic problem, time step At = 0.001 s was used, and it satisfies the Courant condition. Totally, 2 s of the physical experiment was simulated.

The problem of instantaneous application of compressive load an = -3as = -337500 Pa along the axis x3 to an isotropic elastoviscoplastic half-space was considered in three-dimensional case. The velocity of transverse wave in the medium is 2250 m/s. For the case of an elastoplastic medium, analytical solution includes longitudinal plane wave and plastic deformation wave that propagates with lower speed Cf = ^(A + 2^/3)/p. In this case, a zone of constant stress is established between two discontinuities of the solution. The value of this stress acan be also calculated analytically. Using the von Mises yield condition for certain ratio between stress components a11 = a22 = A/(A + 2p)a33, one can obtain

(an - a/3)2 + (a22 - a/3)2 + (a33 - a/3)2 = a2s,a = an + a22 + a33 = 3A+2^a33. (39)

A + 2^

Finally, we have a33 = as. The computational algorithm developed in this work can

be also successfully used if relaxation time t is not small. The calculations were carried out for the problem of deformation of an elastoviscoplastic medium with a linear viscosity condition for non-zero values of dimensionless parameter S. The results are presented in Fig. 1. It can be seen that propagation velocities of both waves and the value of stress a|3 are correctly restored. An increase in viscosity of the medium leads to broadening of the plastic wave front, and it does not affect the front of the elastic precursor.

Further, hardening effect of the EVP medium was studied. Hardening parameter is varied in a wide rage, however, the condition y << 1 was always fulfilled. The simulation results are presented in Fig. 2. For large values of a the effect of hardening significantly influences the propagation velocity of the plastic deformation wave.

Conclusions

The dynamic behaviour of isotropic elastoviscoplastic media with hardening effect under the action of an external load was studied in this paper. The resulting system of equations is semi-linear hyperbolic system. Its structure doesn't allow us to directly apply standard explicit numerical schemes. An approach based on the implicit approximation of the constitutive relations

Fig. 1. The distribution of stress a33 along the axis x3 after 2 s from the beginning of the calculation, an = -3as, medium without hardening

---Viscosity <5 = 10 ..... Viscosity <5 = 10 and hardening a = : 105 -

---Viscosity <5 = 10-2 - Viscosity <5 = 10~2 and hardening a = 102

-

H-1-1-1-r

1.5 1.6 1.7 1.8 1.9 2.0 2.1

u=vTt

Fig. 2. The distribution of stress a33 along the axis x3 after 2 s from the beginning of the calculation, an = -3as, y = 0.044

containing a small relaxation time parameter in the denominator of the non-linear free term was used. In this work, the grid-characteristic method on structured computational grids was used. It allows one to increase approximation order and to avoid strong oscillations of discontinuous solutions. The explicit computational algorithm that has second-order approximation in time and space was developed.

The obtained solutions with the second-order implicit approximation for the deviator part of the stress tensor of the elastoviscoplastic system of equations allow one to consider the limiting case when relaxation time tends to zero. Correction formulas obtained in this case can be interpreted as regularizers of numerical solutions of elastoplastic systems.

The developed algorithms were applied to the three-dimensional direct dynamic problems for elastoviscoplastic media. The obtained results may be useful for the inverse seismic problems [20].

The construction of numerical algorithms for elastoviscoplastic media was carried out by Golubev V. I. and Nikitin I. S. with the financial support of the Russian Science Foundation (project

no. 19-71-10060, https://rscf.ru/project/19-71-10060/). The simulation results were obtained by

Mi Xin with the financial support by the Fellowship of China Scholarship Council.

References

[1] P.Pezhina, Basic issues of viscoplasticity (Osnovnie voprosi vyazkoplastichnosti), Moscow, Mir, 1968 (in Russian).

[2] V.N.Kukudzhanov, Computational continuum mechanics (Vichislitelnaya mechanica splosh-nikh sred), Moscow, Fizmatlit, 2008 (in Russian).

[3] V.K.Novatskiy, Computational continuum mechanics (Vichislitelnaya mechanica sploshnikh sred), Moscow, Mir, 1978 (in Russian).

[4] A.Freidental, Kh.Geiringer, Mathematical theories of inelastic continuum (Matematicheskie teorii neuprugoy sploshnoy sredi), Moscow, Fizmatgiz, 1962 (in Russian).

[5] V.N.Kukudzhanov, Rasprostranenie voln dv uprugovyazkoplasticheskih materialah s dia-grammoi obchego vida, Mehanika Tverdogo Tela, 5(2001), 96-111 (in Russian).

[6] D.Kolarov, A.Baltov, N.Boncheva, Mechanics of plastic media (Mehanika plasticheskih sred), Moscow, Mir, 1979 (in Russian).

[7] G.Duvo, N.Lions,Inequalities in mechanics and physics (Neravenstva v mechanike i fizike), Moscow, Nauka, 1980 (in Russian).

[8] V.M.Sadovskii, , Discontinuous solutions in problems of dynamics of elastoplastic media (Razrivnie recheniya v zadachah dinamiki uprugoplasticheskih sred), Moscow, Nauka, 1997 (in Russian)

[9] A.G.Kulikovskii, N.V.Pogorelov, A.Yu.Semenov, Semenov, Mathematical issues of numerical solution of hyperbolic systems of equations (Matematicheskie voprosi chislennogo resheniya giperbolicheskih system uravnenii), Moscow, Fizmatlit, 2001 (in Russian).

[10] M.L.Wilkins, Calculation of elastoplastic flows. Computational methods in hydrodynamics (Raschet uprugoplasticheskih techenii. Vichislitelnie metodi v gidrodinamike), Moscow, Mir, 1967, 163-212 (in Russian).

[11] M.L.Wilkins, Computer simulation of dynamic phenomena, Berlin-Heidelberg-New-York: Springer, 1999.

[12] V.N.Kukudzhanov, Method for splitting elastoplastic equations (Metod rascheppleniya uprugoplasticheskih uravneniy), Mechanika tverdogo tela, 1(2004), 98-108 (in Russian).

[13] M.Kh.Abuzyarov, V.G.Bazhenov, V.L.Kotov, A.V Kochetkov, S.V.Krylov, V.R.Feldgun, The method of the decay of discontinuities in the dynamics of elastoplastic media, Comput. Math. Math. Phys., 40(2000), no. 6, 900-913.

[14] B.D.Annin, V.M.Sadovskii, On the numerical realization of a variational inequality in problems of the dynamics of elastoplastic bodies, Comput. Math. Math. Phys., 36(1996), no. 9, 1313-1324.

[15] N.G.Burago, Modeling the destruction of elastoplastic bodies (Modelirovanie razrusheniya uprugoplasticheskih tel), Computational continuum mechanics, 1(2008), no. 4, 5-20 (in Russian).

[16] V.I.Golubev, A.V.Shevchenko, I.B.Petrov, Raising convergence order of grid-characteristic schemes for 2D linear elasticity problems using operator splitting, Computer Research and Modeling, 14(2022), no. 4, 899-910. DOI: 10.20537/2076-7633-2022-14-4-899-910

[17] I.B.Petrov, V.I.Golubev, A.V.Shevchenko, Higher-Order Grid-Characteristic Schemes for the Acoustic System, Proceedings - 2021 Ivannikov Memorial Workshop, IVMEM 2021, 2021, 61-65.

[18] V.I.Golubev, A.V.Shevchenko, N.I Khokhlov, I.B.Petrov, M.S.Malovichko, Compact Grid-Characteristic Scheme for the Acoustic System with the Piece-Wise Constant Coefficients, International Journal of Applied Mechanics, (2022), paper 2250002.

DOI: 10.1142/S1758825122500028

[19] A.S.Kholodov, Construction of difference schemes with positive approximation for hyperbolic equations, USSR Computational Mathematics and Mathematical Physics, 18(1978), no. 6, 116-132. DOI: 10.1016/0041-5553(78)90141-6

[20] V.I.Golubev, I.S Nikitin, A.V.Vasyukov, A.D.Nikitin, Fractured inclusion localization and characterization based on deep convolutional neural networks, Procedia Structural Integrity, 43(2023), 29-34. DOI: 10.1016/j.prostr.2022.12.230

Численные схемы повышенной аппроксимации для задач динамики упруговязкопластических сред

Василий И.Голубев

Московский физико-технический институт Долгопрудный, Московская область, Российская Федерация Институт автоматизированного проектирования РАН Москва, Российская Федерация

Илья С.Никитин

Институт автоматизированного проектирования РАН Москва, Российская Федерация

Синь Ми

Московский физико-технический институт Долгопрудный, Московская область, Российская Федерация

Аннотация. Для устойчивого численного решения определяющей системы упруговязкопластиче-ской модели сплошной среды с условием текучести Мизеса и с учетом упрочнения предложена явно-неявная схема 2-го порядка с явной аппроксимацией уравнений движения и неявной аппроксимацией определяющих соотношений, содержащих малый параметр времени релаксации в знаменателе нелинейных свободных членов. Для согласования порядков аппроксимации явного упругого и неявного корректировочного шагов построена неявная аппроксимация второго порядка для изотропной упрочняющейся упруговязкопластической модели сплошной среды. Полученные решения неявной аппроксимации 2-го порядка для девиаторов напряжений упруговязкопластической системы уравнений допускают предельный переход при стремлении времени релаксации к нулю. Корректировочные формулы, полученные таким предельным переходом, можно трактовать как регуляризаторы численных решений упругопластических систем с упрочнением.

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