Numerical simulation of thermoelastic waves arising in materials under the action of different physical factors Текст научной статьи по специальности «Физика»

UDC 519.633: 536.212

Numerical Simulation of Thermoelastic Waves Arising in Materials Under the Action of Different Physical Factors

I. V. Amirkhanov*, I. Sarhadov*, A. F. Ghalebt, N. H. Sweilamt

* Laboratory of Information Technologies Joint Institute for Nuclear Research Joliot-Curie, 66, 141980 Dubna, Moscow Region, Russia ^ Department of Mathematics, Faculty of Science Cairo University at Giza, Egypt, 126613

In this work, a system of equations of one-dimensional thermoelasticity is presented, which takes into account the nonlinear dependence between stresses and deformation. Different problems are considered, which correspond to different external actions on the sample. The numerical simulation of thermoelastic waves arising in a metallic sample under a variable external pressure on the sample's boundary is considered, when the maximum of the absolute value of pressure changes in the interval of observation or violation of the linear Hooke's Law. The numerical investigation of the dynamics of these waves in two-layered structures is also performed. It is shown that the nonlinear dependence of stresses on deformation strongly affects the form dynamics of the thermoelastic wave.

Key words and phrases: thermoelasticity, thermoelastic waves, numerical simulation, nonlinearity, Hooks law, deformation, stress.

1. Introduction

The appearance of optical quantum generators, accelerators of intensive current electron beams and powerful pulsed ion beams has provided a unique opportunity for the study of the irradiation effect of concentrated energy streams on materials (c.f. [1] and the references therein). This is one of the most promising trends of the beam technology from the viewpoint of modification of material properties [1-6] (microsolidity, wear resistance, corrosion resistance, heat resistance, etc.), as the modification provides a way for changing purposely the element composition and the structural-phase state of pre-surface layers by means of a pulsed thermal treatment. The interaction of the charged particle beams with the surface of metals and alloys can cause elastic stresses in the whole irradiated area. This phenomenon evidently plays an important role in the process of the radiation treatment of the surface. When the charged particles affect the surface of the solid, there are different causes leading to the formation of elastic stresses [3]. One of these causes is the local overheating occurring around charges particles trajectories, which leads to the generation of thermoelastic stress [7], and this can exceed the mechanical material strength. The material temperature changes not only because of heat supplied from external sources, but also because of the deformation processes. Elastic and thermal waves are produced by deformation. These phenomena are usually studied in the framework of the theory of thermoelasticity [8]. In previous papers we studied thermoelastic waves arising in materials under the action of pulsed ion beams in the frame of the system of equations of thermoelasticity when the linear Hookes law is valid. In this paper, we shall study the dynamics of the thermoelastic waves arising in materials under the action of different factors for the case when dependence between stresses and deformation is nonlinear.

It is necessary to emphasized, that in the majority of practical cases thermal processes happen under the influence of exterior factors during very short times and in microsize areas, for example for our cases: t ~ 10-7 s, I ~ 10-5 m. It is complicatedly or impossible to spend physical measurements at such sizes and times in the nature.

Received 18th December, 2012. This work is performed under the financial support of the RFBR, grants N 12-01- 00396-a, 11-01-00278-a, Project N203, ARE-JINR cooperation.

In such cases researches resorts to numerical modelling which allows to study physical processes in the extraordinary substance conditions.

2. Basic equations

The basic equations of thermoelasticity are defined by following algorithm:

1. Write down the free energy function of the deformed body F(T,UUij). Here, T is the temperature as measured from a reference temperature and Uij are the components of the strain tensor.

2. Through this function define the stress components 7ij and the entropy of the deformed body S by the formulas:

_ dF_ __ dF 7iJ _ dU~j; S _ - dT •

3. From the known entropy and stress components, set the equations of motion and the equation of heat conduction as:

d2ui dai-j mdS „a^^a ,

»H¥ _ Hxf; T~et _ V(XVT) + A(r' '>•

where X(T) is the coefficient of heat conduction of the solid body, » — its density and A is a source function.

The free energy function of the deformed body has the form [?]: F _ Fo(T) - Ka(T - To)Un + » (l7ik - 3SiU^ + y U2, Uik _ 2 (f| + f^) ,

E E s 1

K _ 3(1—2^) ' » _2(TT^)' -1 <*< 2 '

where Ô is Poisson's ratio, K, » are respectively the modulus of uniform compression and shear, E is Young's modulus, ui is the i-th component of the displacement vector u, Uij is deformation tensor, a _ 3at is the coefficient of thermal volume expansion ( at is the coefficient of thermal linear expansion), 5ik are the Kronecker delta symbols and Fo and To are respectively the free energy and the initial temperature of the body in the undeformed state.

Following the above algorithm, one obtains the stress tensor components and the entropy as:

dF ( 1 \

&ik _ dUj^ _ —Ka(T - To) ôik + 2» (Uik - 3 SikUuj + KUiAk,

S _ -dF _So(T) + KaUii, d T

where S0(T) is the body entropy in the undeformed state. The obtained 7ik and S leads to the next move and heat equations:

d2Ui _ d<7ij _ KadT i fK i »\ d2ui i d2i

K d1 (k »^

dxi V 3/ dxidx] »

dt2 dxj dxi V 3/ dxidxl dx2 '

rrdS rrfdSodTl duii\ dSo Cv

Tm _T{wm +Ka^T)_ V(AVT) + A(-'W _ T'

Here Cv is the heat capacity per unit volume of the body under constant volume — p, ^v iS the specific heat capacity).

Thus in the case of linear Hooke's law, the system of equations of thermoelasticity takes the form:

p^ — -KaVT + [K + |) Vdivu + ^V2u, (1)

dT d

Cv — + Ka-div u — V(XVT ) + A(f, t). (2)

at at

In [9,10] the system of equations (1)-(2) with proper initial and boundary conditions was used for the numerical modeling of thermoelastic waves arising in materials under the action of pulsed ion beams for the one dimensional case. The source function A(r, t) and initial and boundary conditions accompanying system (1)-(2) are defined according to the statement of the physical problem under consideration.

But this system of equations is obtained within the frame of the linearized theory. When this is not the case, the system of equations (1)-(2) cannot be used anymore and must be replaced accordingly, in order to take into account the existing nonlin-earities, for example the nonlinear dependence of stress on strain. In [11], a model of thermo-magnetoelasticity is investigated, which involves the dependence of the specific heat on strain. The resulting equations were used to study the propagation of nonlinear thermo-magnetoelastic waves in materials. A restriction of this model to thermoelasticity is investigated in [12] to study the propagation of nonlinear thermoelastic waves in a half-space. The phenomenon of shock wave formation was put in evidence.

In the present paper, following the above-mentioned algorithm, we shall consider a system of nonlinear equations of thermoelasticity, along the guidelines proposed in [11], with different types of nonlinearities in the one-dimensional case. The free energy function of the deformed body is taken in the form:

F — Fo(T) + ^ + ^Eul + fEui + E 5u| + ■■■-

- uxEaT(T - To)(1 + T\ux/2 + 72^/3) - Ea2T(T - Tq)2(7i/2 + 72^), (3)

where ai, a2, ■■■, n, r2, 71, 72 are dimensionless constants, ax is the coefficient of linear thermal expansion. The application of the above-mentioned algorithm to this free energy function leads to the following equations of motion and heat conduction:

OP+W — o- (4)

/ dv + dv \ da; Ou du Ou + Ou (5)

\dt dx) dx' a dx' di dt dx'

a — Eux(1+aiux+a2u2x+a3u3x +----)-EaT (T-To)(1+nux + T2u2x)-72Ea^ (T-To )2,

p cv + TEa? (71 + 272 ux)

^ + TaTE at

1 + TlUx + T2U2X + 272aT(T - To)

dux ~df

dx^\ dx)

+ A(x, t), X(T ) = A(T ) + [p 1 Ux + P2U2x + P3UxaT (T -To )]A(To). (6)

Here p1, p2, p3 are dimensionless constants, A(T), A(T), A(T0) are respectively the thermal conduction of the material in the presence and in the absence of the deformation, and at room temperature T0. Equation (6) shows that the elastic deformation contributes to the thermal conduction and heat capacity of the material. The deformation thus leads to a change in the material's temperature. In the sequel, it is

assumed that any function f(ux) depending only on the deformation ux in the above equations may be expanded as:

f(ux) = Eux(1 + a\ux + a.2u2x + a3usx +-----+).

For large enough deformations, the material density may change essentially and to take this into account, the continuity equation (4) is added to the system of equations.

If the sample is composed of different superposed material layers, then at the points of contact with coordinates xk, k = 1, 2, ■ ■ ■ ,n (n is the number of contact surfaces) the following coupling conditions must be satisfied:

[7]x=xk _ 0, [T]x=xk _ 0, [u]x=xk _ 0

A—

d x

0.

Here [ /]x=xfc means the usual jump in the value of f at the interface x _ xk passing from one layer to the next. In order to avoid complications in the numerical satisfaction of these conditions at the layer interfaces, we introduce a continuous function which changes smoothly on the few nets points from one value to the other. For example if Young's modulus for one layer has the value Ei and in other layer has the value E2 then for Young's modulus one can introduce the coordinate dependence:

E(x) _ Ei + E2 + E2 2 Ei tanh(»(x - xk)), » > 1,

here tanh(x) is hyperbolic tangent which changes from -1 to 1 when its argument x changes from -œ to to. The coefficient » has no concrete value. Its value is chosen in such a way that on the calculation grids the interval of passing from -1 to 1 of the hyperbolic tangent has enough number of grid points. For example if the space step is h _ 0.001 then may take » _ 200. In this case the mentioned interval contains more than 10 points.

3. Initial and boundary conditions

To the system (4)-(5) we will add the following sets of initial and boundary conditions:

1. »(x, 0)_»o, u(x, 0)_0, v(x, 0)_0, T(x, 0) _ To, A(x, t) _ 0,

Ani

7(0, t)_ -P(t), 7(1, t)_0, —

d x

_0, %

x=o d x

0.

x =

2. »(x, 0) _ »o, u(x, 0) _ 0, v(x, 0) _ 0, T(x, 0) _ To, A(x, t) _ 0, 7(0, t) _ 0

7<'• (>_°- »f

o d T

x=o _g()• 9( ) _ 1 + exp[a*(i - io)] ^ dx

x =

3. »(x, 0)_»o, u(x, 0)_0, v(x, 0)_0, T(x, 0) _ To, A(x, t) _ 0,

O/Ti

7(0, t) _ 0, 7(1, t) _ 0, T(0, t)_ fi (t),—

d x

0.

x =

4. »(x, 0)_»o, u(x, 0)_0, v(x, 0)_0, T(x, 0) _ To, A(x, t) _ 0,

d T

7(0, i)_0, 7(1, t)_0, dx

_0, f

x=0 d x

x =

0,

A(x, t) _ Aq(x)q(t), q(x) _

1

1 + exp[a*(x - xo)]

x = xfc

0

Here I is the thickness of the sample and go, a, to, xo are constant parameters. In all of these four problems the initial conditions are the same and mean that the sample at the initial time moment was in an strain-free state and its temperature was To. Excluding the first problem, in all the other three problems the boundary conditions for the stresses are the same and mean that the sample boundaries are stress free. In the first problem, a variable in time pressure P(t) acts on the boundary x = 0, while the boundary x = lis stress free. The boundary conditions for the sample temperature for the first and fourth problems mean that both boundaries are thermally insulated, while for the second and third problems only the boundary x = is thermally insulated. At the other boundary x = 0 there is a thermal flux of intensity g(t) for the second problem and the temperature varies with time according to a given law fi (t) for the third problem. In the first three problems there is no volume energy supply, while in the fourth problem the sample is heated by a volume source with known intensity. In the papers [9,10] such sources were used for the simulation of the thermoelastic waves arising in metals under an action of pulsed ion beams. The choice of other initial and boundary conditions to accompany the system equations (4)-(5) allows to treat a number of other various problems, but we will confine our attention here only to the above-mentioned four problems.

4. The original field equations in dimensionless variables

For the numerical simulation of thermoelastic waves it is convenient to transform the initial field equations to dimensionless unknown functions and dimensionless independent variables according to the formulas:

- t x _ u ^ T a ^ E _ p

i=—, x = 7, u = y. T = yr. a =—,E =—, p =—, to I I To ao ao po

v = - c (T) = ^^ A(T) = ^ A(T) = ^

Vo ' ^ (T ) = Cv (To)' A(T ) = A (To)' A(T ) = A(To) • In order to preserve the form of the equations, one introduces the parameters

pol2 I

ao = -2- , Vo = -T ■

o o

With this, the equations in dimensionless form reduce to:

I + 1 H =0. (7)

, , /dv dv\ da _ du du

p(£){m + vd~x) = dx;e = us = dx, * = dv (8)

a = /(e) - EeaT(T - 1)(1 + ne + t2s2) - jiEa^(T - 1)2, /(e) = Ee (1 + aie + a2?2 + ^e3 +----),

pcv + p*a?T (71 + 272 e)

dT - ,, _ + 1 p aT

at

1 + ne + T2e2 + 272«T(T - 1) 2

de dl

8 /-dT\ - - - - -= ko(A-8^) +A(x,t), A(T) = A(T)+ pie + P2£2+p3eäT(T - 1). (9)

The constants ko, aT, p* are defined by the formulas:

fcn =

A(To )t 0 Po Cv (To) I2

, aT = ®tT0, p* =

E

Po Cv (To) To

5. The initial and boundary conditions in dimensionless

variables

The original initial and boundary conditions in dimensionless variables take the form:

1. »(x, 0) _ 1, u(x, 0)_0, V(x, 0)_0,T(x, 0)_1, A(x, i)_0,

d T

7(0, i)_ -P(t), 7(1, t)_0,—

x=o dx

x=i

2. »(x, 0) _ 1, u(x, 0) _ 0, v(x, 0) _ 0, T(x, 0) _ 1, A(x, t) _ 0, 7(0, t) _ 0,

7(l,i)_0, -A(T)f

9o dT

x=o _ ^ g(t) _ 1 + exp[a*(t -to)] ^ ~dx

0.

x=i

3. »(x, 0) _ 1, u(x, 0)_0, -(x, 0)_0, T(x, 0) _ 1, A(x, -) _ 0,

- - - - dT 7(0,-)_0, 7(1,-) _0, T(0, -)_/i(-), dx

x=i

4. »(x, 0) _ 1, u(x, 0)_0, v(x, 0)_0, T(x, 0) _ 1, A(x, t)_0,

d T

7(0, t)_0, 7(1, t)_0, —

A(x, ) _ A (x) ( ), (x) _

_0. £

x=o dx

1

x=i

0,

1 + exp[a*(x - xo)]

Here a* dimensionless parameter, while the dimensionless source A(x, t) and energy flux go are obtained from their original functions A(x, t), go by the formulas:

^ A(x^ fto

A(x, t) _ ^ /rr ^ m , go

cv (To )poT0' X(T0 )T0-

In what follows the bar will be omitted and all variables are meant as dimensionless.

6. The method of numerical solution

For the numerical simulation we will use the dynamical adaptation method [13,14]. The accuracy of the numerical solution of the partial differential equations strongly depends on how good the distribution of the grid nodes is in agreement with the character features of the solution.

The construction of optimal grids can be achieved by means of various adaptation methods, to the development of which special attention is paid in the past few decades [15,16]. The method of dynamic adaptation relies on the procedure of transition to any arbitrary non-stationary system of coordinates. Such a procedure allows to formulate the problem of construction and adaptation of a computational grid on the differential level, so that in obtained mathematical model one part of the differential equations describes the physical processes, while the other part describes the behavior of the nodes in the physical grid [13,14]. The transition to an arbitrary system of coordinates is carried out by means of an automatic transformation of coordinates, with the help of the required solution.

Following [13,14], we will perform the transition from the physical space Qx,t to the computational space with arbitrary non-stationary system of coordinates Qq,T and variable coordinates ( q, r). The mentioned transition may be performed by using a

0

0

change of variables of general form x = ( , ), = , for which there exists a unique non-degenerate transformation q = ip(x, t), t = t.

The partial derivatives of the dependent variables for the transition from one system to the other are calculated from the following expressions:

d d Q d d 1 d dx ^ dx

— =--h — —, — =--, tf = —, Q =--,

dt dr tf dq' dx tf dq ' dq ' dr'

(10)

where ^ is the Jacobian of the inverse transformation. Using these transformations, the original system of equations in dimensionless form (7)-(9) may be rewritten in the variables and in the conservative form:

d_ dr

W

d_ d

(Q + v)p

0,

d_ dr

(pv tf)

d_ d

da

(Q + v)vp = —, a = a( £ ,T),

=

1 du tf dq '

d_ dr

(utf)

d_ d

(Q )u

= v,

(11)

(12) (13)

[p cv + p*a^T (71 + 272 e)j

Ldr (*T ) d

d

+ Tp*aT[1 + rie + T2£2 + 272ŒT(T - 1)] tf— + (Q + v) —

d

+ d

d£ 1

d

= ± I A dT )

dq V tf dq z1

(14)

0, 0 <q< 1, t> 0. (15)

To this system we add the initial and boundary conditions for the first problem:

dtf dQQ

d d

p(q, 0) = 1, u(q, 0) = 0, v(q, 0) = 0, T(q, 0) = 1,

a(0, r) = -P (t),

d T

a<l, r) = 0, -

q=„

0, f

d

0,

q=i

(16) (17)

The functions a(e,T), p(e) and A(T) are the same functions which were introduced before with a "bar". Here the bar has been omitted.

Thus, after passing to an arbitrary non-stationary system of coordinates, the original differential model transforms into an expanded differential model in which there appears an additional equation of of the type (15). Its type, properties and the form of the boundary conditions depend on the concrete form of function Q. The function Q at this stage is left unspecified. Once it has been specified, the equation (12) will be used for the construction of a grid that adapts to the solution. The ways according to which function Q can be chosen are explained in [13,14] for different problems. As seen, for the choice Q = -v the system (12)-(15) takes the form:

d d a d u

d? = d? u = = a(e'T),

rt^ rt^ ft / ft l \

[cv + p*o?TT(71 +272£)jir + Tp*Oct[1 + Tie + t2s2 +272aT(T- 1)]ir = k0— À— , (19)

=

d u d

(18)

de

d f-dT^

dr

dr dq\ dq J '

u(q, 0) = 0, v(q, 0) = 0, T(q, 0) = 1,

(20)

7(0, T)_ -P(t), 7(1, r) _ 0, dT

q=o

0, f

d

q=i

0,

(21)

Such a choice of transformation of coordinates corresponds to the Lagrangian mass variable. The coordinate q determines the position of the body in the undeformed state and u defines the displacements relative to this state.

As an application, the numerical simulation will be carried out for the iron sample used in [9,10]. The model parameters cv o, A (To), po, E, ar are those of iron at room temperature T0 = 300 K.

Cv (To) = 456 J / (kgK), A (To) = 78, 2 W/(mK), p0 = 7870 kg/m3,

E _ 2.02 ■ 10ii Pa, aT _ 3.6 ■ 10—5K-i,t0_3 ■ 10—7 s, I _ 10

-5 r^-i

1-7

1-5

m.

For these values, the dimensionless parameters ko, aT, p*, E take on the values:

k0 ~ 6, 54 ■ 10

2

ut _ 1.08 ■ 10-2, p* ~ 187.6, E ~ 2.31 ■ 104

The dimension scale for the stress is cto = 8.744 ■ 106 Pa. Thereby all calculable dimensionless constants are determined. There remains only to determine the values of the remaining constants 71, 72, n, r2, pi, p2, p3. The dimensionless functions /(e), cv(T), A(T), P(t) will be given using the experimental data for these relations. In order to avoid the false numerical oscillations appearing during the numerical simulation, we introduce artificial viscosity to the stress function a(e,T). This is zero for positive values of the velocity gradient and differ from zero otherwise:

7(e, T ) ^7(£ ,T ) + i7( % -

d

d

7. Numerical scheme

For the numerical simulation of the original problem we shall use the method of finite volumes [?,17]. In the computational space fiq,T we introduce a computational grid in which for convenience we use nodes with integer and half-integer indices:

= {(qi,^), ( %+i/2,tJ); qi+i = qi + Aq,

qi+i/2 = qi + A q/2, i = 0,1, ••• ,m; r^1 =r3 + Ar, fk =/(^, V), j = 0,1, ••• ,n}.

The values of the displacement and velocity will be taken relative to the integer cells, while those of temperature, deformation and stresses are taken relative to the halfinteger cells. Explicit or implicit numerical schemes can be used. Here we use an explicit scheme:

- 4 = %2 -^/2 - ; . = 0 (22)

Ar = 3A q ; ' = ° (22)

i+i j J _ J

= ^+1/2A "i-1/2; ¿ = 1,2,... m - 1, (23)

A A

- < _ —i/2 - 7rn-3/2 - 8^

m

; « _ m, 7m _ 0, (24)

Ar 3Aq

IV.

i_0,1, ■■■ ,m (25)

7 + i Î 7 + i i 7

K - u _ vy + vi ^ . _

A 2

Jij + l _ rpj

— J+l — J r0 ' 1/2 T1/2 + 1/2 1/2 _ , _

°i/2 Ar + ^/2 At 0 Aç

W -W1; ,_0,

(26)

rpj + 1 _ rpj rj Ti+1/2 Ti+1/2 + Fj °i+1/2 Ar + »+1/2

eÎ+1/2 ei+1/2

A

_ fco

^ _ w3

' '%+1 " t ,

Aq :

i_1,2, ••• , m - 2, (27)

TJ+1 — Ti PJ+1 — ^

T m-1/2 T m-1/2 + fcm-1/2 m-1/2

— 1 /2 Â r " m

"m-1/2 Ar + m-1/2 Ar

wm -wm 1

_ fco m A m-1 ; î_m - 1, (28)

T T

,,,] u i+1/2 ¿-1/2 • -i o -i

W, = —A-;—--—, i = 1, 2, ••• ,m — 1

u° = 0; = 0; Tf+1/2 = 1; i = 0,1, ••• ,m — 1. (29)

r° = —Pj; < = 0; V(0) = 0; W\m) = 0; (30)

^^ = 4 + P*^ (71 +272£ j), H = P*«TTj [1 + n ej + T2 ej2 + 272«t (Tj' — 1)], k = i + 1/2.

rP = Tj ) ^ = u^+1 — U

"¿+1/2 r(t¿+1 /2' T ¿+1 /2)' '¿+1/2 ^ '

To solve the system (22)-(28) with initial and boundary conditions (29)-(30), we first solve the equations (22)-(25) relative to v3i+1, u{+1 and then from the known values of £i+1/2 and £i+1/2 we solve equations (26)-(28).

8. Discussion of the numerical results

Investigation of the mechanical waves arising in the sample under the action of the variable pressure with time on the boundary x = 0 will be carried out for the uniform sample at first, then for the two-layered sample. The function expressing the nonlinear dependence of stress on the deformation we will taken in the form:

f(e)_a*( 1 - e-.

As we see for ^ 1 the given dependence is that of the linear Hooke's law f(e) = £E. Figure 1 shows the mentioned dependence for r* = 100.

For the first case we will investigate the stress waves without taking into account the change in temperature, i.e. we assume T = T°.

Now we consider the particular case when the time dependence of pressure on the boundary is defined by Gauss's formula and the maximum value of the pressure lies, either in the region of validity of the linear Hooke's law, or in the region of its nonlinearity:

P(t) = P° exp[—(t — r°)2/2r2].

Figure 2 shows the given dependence for ro = 0.01, Po = 10, aT = 2 ■ 10 4 which satisfies Hooke's linear law.

Figure 1. Nonlinear dependence of Figure 2. Time dependence of pressure stress on deformation f (e) for on the boundary x =0 for tq = 0.01,

a = 100 Pq = 10, aT = 2 • 10-4, within the domain

of linearity of Hooke's law

Figure 3 shows the profiles of the stress wave a(x, t) at the different time moments t = Ti, i = 1,2, ■ ■ ■ , 6 (t! = 0.011, Ti = Ti + (i - 1) ■ 0.002), which arise in the sample under the action of an external variable pressure on the boundary x = 0. As seen, the localized stress wave arising in the sample moves from the boundary x = 0 to the boundary x = 1 and getting this boundary the stress wave changes its polarity and returns to the boundary x = 0. At this stage, our numerical simulation confirms the results of the wave equation solution.

As mentioned earlier, the amplitude of the external pressure was taken not very large, in order not to violate the linear Hooke's law. Now let us remove this restriction and choose a large amplitude of the external pressure so as to work in the region of nonlinear dependence of stress on deformation (Po = 90), and watch the dynamics of the wave form at the same previous time moments. Fig. 4 shows the evolution of the wave form at those time moments. As seen, the nonlinearity strongly distorts the form of the wave: The amplitude of the wave decreases and its width increases. The change of the wave form happens during some interval of time, beyond which no further changes of the wave form are noticed.

t-,. „ t-, „ ,, , Figure 4. Profiles of the stress wave

Figure 3. Profiles of the stress wave N , , ,.

/ \ ^ ^ ^ a(x,T) at different time moments

alx.T) at different time moments v ' . , __ _ ,

v ' 7 t = Ti, i =1, 2, • • • , 6 (n = 0.011, Ti =

Ti + (i - 1) • 0.002, i = 2, 3, • • • , 6) when the linear Hooke's law is violated

T = Ti, i =1, 2, • • • , 6 (Ti = 0.011, Ti = Ti + (i - 1) • 0.002, i = 2, 3, • • • , 6)

Now we follow the evolution of the wave form arising in the sample under the action of the same external pressure, when our sample is formed of two layers of the same thickness, but with different Young's moduli Ei = 2.31 ■ 104, E2 = 4.62 ■ 104. Fig. 5 shows these profiles at the same previous times ri, i = 1,2, ■ ■ ■ , 6.

a(x,x)

.0 02 0.4 0.6 0.8 1.0

.0 42 0.4 0.6 0.8 1.0

<J(X,T)

0.2 0.4 0.6 0.8

12 ст(Хд) 8 4 0

-4-8

-12--

0.0

0.2 0.4 0.6 0.8

Figure 5. Profiles of the stress wave a(x,T) at different time moments in a two-layered sample for the case of linear Hooke's law

It is thus seen that after the wave gets to the common boundary of the layers, a part of the wave passes to the second layer and some part of the wave reflects back to the boundary x = 0.

9. Conclusions

In this paper, we carried out an investigation of the propagation of thermoelastic waves within the frame of a nonlinear model. Using numerical modeling, the following results were reached:

1. In the domain of validity of the linear Hooke's law the form of a thermoelastic wave almost does not vary with time.

2. Under violation of the linear Hooke's law, there is an essential modification of the waveform with time.

3. In two-layered samples, partial reflexion of the wave takes place at the contact interface.

References

1. Бойко В. И., Валяев А. Н., Погребняк А. Д. Модификация металлических материалов импульсными мощными пучками частиц // УФН. — 1999. — Т. 169, № 11. — С. 1243-1271. [Modification of metallic materials by intense pulsed particle beams / V.I. Boyko, A.N. Valyaev and A. D. Pogrebnyak // UFN. — 1999. — Vol. 169, No 11. — P. 1243-1271 ]

2. Stinnett R. W. et al. Thermal Surface Treatment using Intense, Pulsed Ion Beams // Materials Research Society Symp. Proc.: Materials Synthesis and Processing Using Ion Beams. — 1993. — Vol. 316. — Pp. 521-532.

3. Диденко А. Н., Лигачев А. Е., Куракин И. Б. Воздействие пучков заряженных частиц на поверхность металлов и сплавов. — Москва, Энергоатомиздат, 1987. — 184 с. [Didenko A.N., Ligachev A.E., Kurakin I. B. The Interaction of Charged Particle Beams with the Metals and Alloys Surfaces, 1987.— 184 p. ]

4. Коренев С. А. // Препринт ОИЯИ P13-89-615. — Дубна, 1989. — 16 с. [Korenev S. A. Preprint of JINR P13-89-615, Dubna, 1989, 16 p. ]

5. Действие излучения большой мощности на металлы / С. И. Анисимов, Я. А. Имас, Г. С. Романов, Ю. В. Ходько. — М.: Наука, 1970. — 272 с. [Anisimov S. I., Imas Ya. A., Romanov G. S., Khodiko Yu. V. Action of High-Power Radiation on Metals. — M.: Nauka.— 1970.— 272 p. ]

6. Быковский Ю. А., Неволин В. Н., Фоминский В. Ю. Ионная и лазерная имплантация металлических материалов. — Энергоатомиздат, 1991. — 237 с. [Bikovskiy Yu. A., Nevolin V. N. Fominskiy V. Yu. Ion and Laser Implantation of Metallic Materials. — Energoatomizdat, 1991.— 237 p. ]

7. Калиниченко А. И., Лазурик В. Т. Возбуждение акустических колебаний пучком заряженных частиц малой плотности // ЖЭТФ. — 1973. — Т. 65, № 6. — С. 2364-2368. [Excitation of Acoustic Oscillations by a Low Density Beam of Charged Particles / A.I. Kalinichenko and V.T. Lazurik //ZhTPh — 1974. — Vol.38, No 6, — P. 2364-2368 ]

8. Ландау Л. Д., Лифшиц Е. М. Теория упругости. — 1970. — 248 с. [Landau L. D., Lifshitz E. M. Theory of Elasticity. — 1987. — Vol. 7. — 248 p. ]

9. Numerical Simulation of the Thermoelastic Effects in Metals Irradiated by Pulsed Ion Beam / I. V. Amirkhanov, E. V. Zemlyanaya, I. V. Puzynin et al. // JCSME. — 2002. — Vol. 2, No 1s-2s. — Pp. 213-224.

10. Численное моделирование влияния вязкости на распространение термоупругих волн в образце металла, облучаемого импульсными пучками ионов / И. В. Амирханов, Е. В. Земляная, И. В. Пузынин и др. // Поверхность. — 2004. — № 10. — С. 92-99. [Numerical Simulation of Influence of Viscosity on the Propagation of Thermoelastic Waves in Metallic Samples Irradiated by Pulsed Ion Beams / I. V. Amirkhanov, E. V. Zemlyanaya, I. V. Puzynin at al. // Surface. X-ray, Synchrotron and Neutron Investigations. — 2004. — Vol.10. — P. 92-99 ]

11. Ghaleb A. F., Ayad M. M. Nonlinear Waves in Thermo-Magnetoelasticity // Int. J. Appl. Electromagn. Mat. Mech. — 1998. — Vol. 9, No 4. — Pp. 339-357.

12. Rawy E. K., Iskandar L., Ghaleb A. F. Numerical Solution of a Nonlinear, One-Dimensional Problem of Thermoelasticity // J. Comp. Appl. Math. (C.A.M.). — 1998. — Vol. 100. — Pp. 53-76.

13. Мажукин А. В., Мажукин В. И. Динамическая адаптация в параболических уравнениях // ЖВМ и МФ. — 2007. — Т. 47, № 11. — С. 1913-1936. [The Dynamical Adaptation in Parabolic Equations / A. V. Mazhukin, V. I. Mazhukin // ZhVM and MPh. — 2007. — Vol.47. -N11.,— P. 1913-1936 ]

14. Бреславский П. В., Мажукин В. И. Алгоритм численного решения гидродинамического варианта задачи Стефана при помощи динамически адаптирующихся сеток // Математическое моделирование. — 1991. — Т. 3, № 10. — С. 104-115. [Computational Algorithm of a Hydrodynamical Version of Stephan Problem by Dynamic Adapting Grid / P. V. Breslavskiy and V. I. Mazhukin // Mathematical Modelling. — 1991. — Vol. 3, No 10. — P. 104-115 ]

15. Thomson J. F., Warsi Z. U. A., Mastin C. W. Boundary-Fitted Coordinate Systems for Numerical Simulation of Partial Differental Equations-A Review // J. Comput. Phys. — 1982. — Vol. 47, No 1. — P. 108.

16. Thomson J. F., Warsi Z. U. A., Mastin C. W. Numerical grid generation. Foundation and Application., 1982. — New-York:Nort-Holland, 1985. — 483 p.

17. Ильин В. П. Методы конечных разностей и конечных объемов для элип-тических уравнений. — Новосибирск, Издательство Института математики, 2000. — 346 с. [Ilin V.P. Royak M.E., Njvjsibirsk. Finite Differences and Finite Volumes Methods for Elliptic Equations — 2000—346p. ]

УДК 519.633: 536.212

Численное моделирование термоупругих волн, возникающих в материалах под действием различных физических факторов

И. В. Амирханов*, И. Сархадов*, А. Ф. Галеб^ Н. Х. Суейлам^

* Лаборатория информационных технологий Объединённый институт ядерных исследований ул. Жолио-Кюри, д.6, Дубна, Московская область, 141980, Россия ^ Кафедра математики, Факультет науки Каирский университет, Гиза, Египет, 12613

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

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