Модель напряженно-деформированного состояния окрестности трещины

Известия Тульского государственного университета Естественные науки. 2009. Вып. 2. С. 74-89

^ Механика ^

УДК 539.375

Модель напряженно-деформированного состояния окрестности трещины

М.В. Гаврилкина, В.В. Глаголев

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

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

Механика разрушения к настоящему времени имеет множество подходов к описанию напряженно-деформированного состояния (НДС) тел, ослабленных трещиной. Каждый из подходов включает в себя как модель трещины, так и соответствующий критерий ее продвижения. Наиболее распространенной моделью трещины в сплошной среде является представление в виде математического разреза [1, 2]. Используя критерий Гриффитса [3] в энергетическом или силовом варианте удастся прогнозировать трещиностойкость материалов, поведение которых описывается моделью линейно упругого тела. Если же материал проявляет пластические свойства, то модель математического разреза позволяет исследовать различные критерии начала образования новых материальных поверхностей [4] в зависимости от вида нагружения и степени пластической деформации. В этом случае ни одно из представлений не даст ответа на вопрос о начале пластического деформирования в концевой

Работа выполнена при поддержке РФФИ (проект №07-01-96402) и АВЦП «Развитие научного потенциала высшей школы» (проект №2.1.1/941).

области трещины Гриффитса. В частности, наиболее известная и используемая модель Леонова-Панасюка-Дагдейла [5, 6] предполагает наличие пластической зоны при сколь угодно малой внешней нагрузке. Однако определение начала пластического деформирования в вершине трещины является важным вопросом [7], в частности, при циклическом нагружении поврежденного материала. Таким образом, разработка математических моделей, адекватно описывающих зарождение и развитие пластических областей в вершине трещины является достаточно актуальным.

Далее трещину будем рассматривать как физический разрез. Толщину разреза выбираем как минимально допустимую с точки зрения выполнения гипотез механики сплошной среды. В совокупности с данным представлением в модель трещины включаем и материал среды, лежащий на мысленном продолжении физического разреза. Деформация слоя формирует граничные условия для сопряженной с ним среды. В работах [8-10] соответствующий материал определен как слой взаимодействия. Основными задачами при данном подходе являются определение численного значения толщины слоя, включая постановку возможных натурных экспериментов по его нахождению, а также исследование конкретных краевых задач механики разрушения.

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

Рассмотрим нагружение берегов трещины, моделируемой физическим разрезом в линейно упругой плоскости согласно схеме, показанной на рис. 1, соответствующей разрушению типа нормального отрыва. Предполагается возможным существование пластической области (088;0;) в пределах слоя взаимодействия.

Рис. 1. Схема нагружения плоскости

Наряду с напряжением <хц (жг) в слое учитываем напряжение (Т22 (жг), обусловленное наличием касательных напряжений а21 (ж 2) вдоль границы с полуплоскостью.

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

где а — девиатор тензора истинных напряжений; ё — девиаторная составляющая тензора деформаций; 9 = £ •• Е, К — модуль объемного сжатия; р = |<7 •• Е — гидростатическая составляющая тензора напряжений; э — параметр упрочнения; б'(э) — сдвиговой модуль; С{э) = Су при а   а ^ т^п: С(э) = при а  а > т£г; тпг — предел текучести.

В силу симметрии задачи рассмотрим только верхнюю полуплоскость (х\ ^ <5о/2) (рис. 2), а действие слоя заменим нагрузкой <?(ж) = — (ацвх + <?21ё2) (здесь и далее ж = Х2/5о — безразмерная координата; с= с• Р г, .7 = 1,2 — безразмерные напряжения; (3 = ^ — параметр

материала, Е — модуль упругости; р — коэффициент Пуассона).

Основными неизвестными компонентами слоя считаем средние напряжения, определяемые следующим образом:

В силу антисимметрии касательных компонент тензора напряжений относительно ОСИ ОХ2 имеем <712 () = <?21 () = 0. СрСДНСС напряжение <722 (ж) слоя будет связано с касательными напряжениями, действующими вдоль границы слоя соотношением

Да = 2(7(э)Дё, Ар = К АО,

г II

Рис. 2. Эпюра нагрузки

Соотношения Фламана [11] связывают внешние нагрузки <тц и а\2 с перемещениями границы безразмерными выражениями

щ(х) = -Р 1п + £ап(0 1п ^

^ оЫО 1п ^

здесь щ = и^/5о, г = 1,2 — безразмерные перемещения; Р = Р(3/8о — безразмерная сила на единицу толщины; Ь — удаленная точка с нулевым перемещением; Ь — расстояние от начала координат до Ь.

Перемещения границ слоя определяются в следующем виде:

«1

РХ

«2

Напряжения до достижения предела текучести связаны с деформациями законом Гука

£ц = Асгц — В <7 221

£22 = ^4 О" 22 ~~ В<Т 11,

О733 = V

где А = §, В = 2(1-и) — безразмерные постоянные.

Поведение материала слоя на стадии пластического деформирования согласно (1) определяется следующими выражениями:

£п = С + Арап - Вра22,

£ 22 = С\ + АрТт 22 — ВрСТц,

*33 = ^33 + иР

где

0 ’ /?

ЗА' - -Юр ЭА'Ср

• Еп =

р 6А' + ’Юр ’ р ЗА+С/

С = (Л - -(В- Вр)ак22; ^1 = (Л - Лр)а|2 - (В - Яр)<7*1;

*11 > *22! *!з — критические напряжения, соответствующие переходу материала слоя в пластическое состояние.

Подставим в формулу (4) выражения (6), (8) и (11) в результате получим уравнения относительно ?ц и сг 22

В левую часть (15) подставим (9) и (12) и получим следующие уравнения относительно СТц, (722 И <721:

Таким образом, имеем систему интегральных уравнений (14), (16), дополняемую связью (3).

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

где <7 = <7 • (3 — безразмерный девиатор тензора напряжений; тт — безразмерный предел текучести.

Если <7Т = <7Т • 0 — безразмерный предел текучести при одноосном растяжении ((7ц = <7Т, <722 = <733 = 0), ТО ПрСДСЛ ТСКуЧССТИ Тт ОПрСДСЛЯСТСЯ следующим соотношением:

Тт = \^2аТ.

Запишем условие перехода в состояние пластичности через главные компоненты тензора напряжений <тц и (722 для случая плоского деформирования

где С = у (и — 1).

Используя условие (19), соотношения (14), (16), (3), запишем соответствующие системы уравнений для пластической области, области перехода из упругого в пластическое состояние и упругой области.

ИЛИ

Система уравнений, описывающая поведение пластической области слоя взаимодействия

1

х + а

- {С + Арсгц - Врсг2г) = -Р 1п

А ь + а

+ [ *п(£)1п J о

(^,

С\ + АрСГ 22

Вра ц =

ло

*2і

(х - О

да 22

~д7

= —2521.

Области слоя, переходящей из упругого в пластическое состояние

Г 1

х + а

(Аац - Ва22) = ~Р\п

А Ь ~т &

+ [ *іі(£)1п ./о

Аа

22

Ва п = 52і

ло

(х - О

ЙСГ22

= —2(721,

„ ^*11 + *22 “ *11*22 + С ((7ц + (Т22)2 = *т-

Система уравнений для упругого поведения слоя

1 , т— й— \ й, х + а -(Асти — Ва 22) = — Р ш

Л

г , •[ *п(С)1п^7—

Ь + а Уо Ь - £

^ Л(Т22 - Бац = а2і(0

Jo

(Уа22

{х - О

дх

= —2521.

Результатом решения систем (20)—(22) является нахождение полей средних напряжений ац и (722! касательных напряжений по границе со слоем 521, СИЛЫ Р С учетом граничного условия 522105=0 = 0.

Поведение материала слоя взаимодействия будем описывать в рамках дискретной модели, представляя его набором взаимодействующих квадратных в плане ^о-элементов. Основным постулатом данной модели является положение об однородности НДС в каждом из элементов.

Для решения задачи разобьем границу полуплоскости ОЬ на N граничных элементов [12]. Считаем, что каждый граничный элемент характеризуется постоянным (средним по элементу) значением напряжения а^\ где к = 1, N, ; = 1,2.

Построим дискретные выражения интегральных операторов

Проинтегрируем соотношение (3) по 2 элементу

Подставив (23)-(25) в (20)—(22), получим дискретную модель упругопластического поведения материала слоя взаимодействия.

Для пластической области

= -Р1п

* С\ + АрСТ 22 — ВрСГ ^

2

Для элемента, переходящего из упругого в пластическое состояние

2

Система уравнений для упругой области слоя

N с

Ы» - _ Р . Ь, ХІ + а а. Гк+1 1п -3

к=1

я

ж і

~0) _ 22 22

- 2 '

Система уравнений содержит Зп + 1 уравнение от Зи + 1 неизвестного. Неизвестными являются Зп обобщенных напряжений и сила Р, обеспечивающая соответствующее напряженное состояние.

Нахождение основных характеристик сводится к решению соответствующих нелинейных систем уравнений. На первом этапе решения одной из важных задач является определение корректного начального приближения. Для начального приближения используется решение соответствующих линейных систем уравнений, в которых в качестве критерия перехода в пластическое состояние вместо условия (19) использовался критерий Трсска-Ссн-Венана.

Из решения соответствующей упругой задачи имеем: сгц > <733 ^ <722; следовательно, критерий Трсска-Ссн-Венана имеет вид

*11 - *22 = 2т*,

где п = Т* • 0, П — предел текучести.

Или, если 2п ассоциировать с пределом текучести при одноосном растяжении, то

*11 - *22 = *т-

Образование пластической зоны рассматривается как последовательный переход элементов слоя в пластическое состояние. На каждом этапе решения упругопластичсской задачи учитывается перераспределение напряжений, вызываемое переходом нового элемента в состояние пластического деформирования. При решении задачи учитывался вариант, при котором *11 > *22 ^ *33- В этом случае в качестве приближения к критерию (19) использовалось условие

*11 - *33 = *Т-

Было проведено исследование влияния числа граничных элементов на результаты расчетов. Установлено, что при N = 200 элементов отношение

\ап (М) — <7^ (А7' — 100)|

= --------ттт----------------

<7^

приблизительно равно 5.842 • 10-5, при N > 1000 элементов — 1.264 • 10-6.

На рис. 3 построены эпюры распределения напряжений в слое взаимодействия на первых 6 элементах при следующих расчетных характеристиках: N = 1000; а = 10; = О.Шу; Е = 2.1 • 105 МПа; стт = 600 МПа; и = 0.25 для

случая пластического деформирования первого элемента и перехода второго в состояние текучести. Непрерывные линии соответствуют решению нелинейных систем (26)—(28), а штриховые — начальному приближению. Кривые 1 и 4 определяют напряжение стц, 2 и 5 — <722, 3 и 6 — <733.

1

ч\ ч\ ч\ ч\ Л/ 5 4

/ чуу/ , ^ /3 2 г

/

0 -----------------------------------------------------------------------------------------------------------

1 2 3 4 5 Л

Рис. 3. Распределение напряжений для случая пластического деформирования первого элемента и перехода второго в состояние текучести

Проведено исследование влияния параметра упрочнения на распределение напряжений в зоне предразрушения. На рис. 4 построены эпюры распределения напряжений в слое взаимодействия на первых 6 элементах для случая перехода третьего элемента в состояние пластичности при пластическом течении па первых двух. Непрерывные липни соответствуют Ор = О.Юу, а штриховые — Ор = 0.001Сгу. Кривые 1 и 4 определяют напряжение <тц, 2 и 5 — ^22, 3 и 6 — ст33.

Из рисунка видно, что данный параметр не существенно влияет на распределение напряжений в зоне предразрушения.

Исследовался характер распространения пластической зоны в слое взаимодействия при идеально упругопластичсском поведении материала [13] для двух критериев перехода в состояние пластичности: Трсска-Ссн-Вснана и Губсра-Мизсса. В этом случае предполагалось, что справедливо следующее разложение приращения полной деформации через упругую и пластическую составляющие:

Де = Дер + Дее,

— 3

&ІІ

Тт

2.5

1.5

0.5

1

\ У £ 4 /

лЯЧ V / У

1 2 3 4 5 ^

Рис. 4. Распределение напряжений при переходе третьего элемента слоя в

состояние пластичности

где є — тензор полных деформаций, єр — тензор пластических деформаций, єе — тензор упругих деформаций.

Перепишем (31) в виде

єіі = єіі + єііі

где I = 1, 2.

В случае плоского деформирования будем предполагать, что £р3 = е|3 = 0. Таким образом, из условия пластической несжимаемости имеем

+ Є99 — 0-

Подставив в формулу (32) выражения (4), (6), (8), (9), (15) и (33), получим уравнения относительно <7ц, <722 И <721:

х + а

(А - В) (<7ц + <722) = -2РІГ

її "Ь сі

2 / <711 (О 1п , Jo

I

*2іШт------77

о ( - О

Л<7

22

В а ц =

Ґ

Jo

*2іШ

1

Отмстим, что первое уравнение системы (34) описывает упругое изменение относительного объема элемента слоя и остается справедливым как на стадии упругопластичсского, так и упругого деформирования. В качестве критерия перехода из упругого в пластическое состояние, используем соотношение (19), а для начального приближения решения нелинейной системы

критерий (29) или (30) в зависимости от значения НДС на пластически деформируемых элементах.

Таким образом, имеем соответствующие системы интегро-дифференциальных уравнений.

Система уравнений для поведения материала слоя в пластической области

ж + а

{Л-В) (стц + стгг) = -2Р1

да 22

~д7

П

Ь/ -{- (X

рь I т _ СI

-2/ агц(С) 1п ~7------------

J о

= —2(721,

[ *21(0 J О

1

(х - О

к V*11 + *22 - *11*22 + с ((7ц + (722)2 = 5Т.

Для области, переходящей из упругого в пластическое состояние

х + а

А — В) ((7ц + (722) = —2Р 1п

да22

~д7

= — 2^21,

Л(722 - Ва\\ = ^21(0/ ,ч

Уо ( - о

\ZcTii + *22 - *11*22 + с {ап + а22)2 = а?

1

1

1

*21(0, V

о ( - О

Система уравнений для упругой области слоя

х + а Ь + а

2/ стц(С) 1п -^7—Т^+ [ *21(0 .) о ^ € Уо

А — В) ((7ц + (722) = —2Р 1п

#(722 ~дх~

(ж - О

= —2521,

Л(7

22

5*11 = *21(0

Jo

(ж - О

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

(23), (24) и соотношения (25), получим дискретное представление для пластической области

и) , _ о 3 1„ хз + а Л.

А - В) { ] = -2 Р 1п

"Ь а

2

+2 С1п ('^тУ+£ Г ^£'

,аи) - аи^1}

~0) _ 22 22

- 2 ’

, 1^з)2 , ^О')2 -=(з)-=(з) , Р АД?) « — я

к у а11 + 22 ~~ *11 *22 + С I *11 + *22 ) — *т-

Для элемента, переходящего в пластическое состояние {Л-В)(а^ + аЩ = -2Р\п ^ +

+ 25^^п / 1п(^ £ + / 7"^7^’

^22} - [ 7^7^’ ^

^ ^ ^ - *

,аи) - а{^1}

~0) _ *22 *22

- 2 ’

, 1^з)2 , ^О')2 7=г(^)-=гО') , Р ^0') , 7^0'Л

у сг11 + сг22 - сг11 сг22 + С ^сг11 + сг22 ^

Для упругой области слоя

' {А-В) (а^ + аЩ = -2 Р 1п '

+2 Іґ* С'1п (^§?й’ Г ^

^22} - [ 7^7^’

•/£*. хз - ?

к=1

.ст

~0) _ *22 *22

- 2 '

Из распределения напряжений на элементах слоя взаимодействия, можно определить НДС в области, прилегающей к слою со стороны плоскости. Соответствующую область, показанную на рис. 5, представим рядом квадратных элементов с единичной стороной. Считаем, что в пределах каждого

элемента НДС однородно и соответствует осреднснным напряжениям [14] определяемым следующим образом:

1 гх

к+1

10

0%

где Д — площадь квадратного элемента с единичной толщиной.

Рис. -5. Дискретное представление области, прилегающей к слою со стороны

плоскости

В этом случае, зная нагрузку на полуплоскость в виде распределенной (ж) и с/2 (ж) со стороны слоя и сосредоточенной Р, можно определить напряженное состояние в любой точке полуплоскости, в том числе и в точках с координатами ж™ и ж2, следуя принципу суперпозиции

2

7Г

Р (ж?)3 гЬ

\2 п

[ Я1Ш Jo

Ж

2ч2

I

Ь (грП _ С\^ _га

Я2 Ш , (Ж2 Ж2

2

7Г

((?-02 + (?)2)2 42(Ж?,Ж?) =

Р (ж2 + а)2 ж” [ь (х%-£)*х

((ж2 + а)2 + (ж”)2)2 ! О

[ 41 и) Jo

2 „га 1

2\ 2

(!£,

Iі

J о

42 (О

О (ї

2\2

2

7Г

Р (ж? + а) (ж?) ((ж£ + а)2 + (ж?)2)2

21\‘П’Л2. 2 />Ь

[ чі(0 J о

(2 - О (ж?У

2\ 2

/'

Jo

Я2 {О

Подставив (42)—(44) в (41), получим осрсднсннос распределение напряжений в квадратном элементе с единичной стороной.

На рис. 6 построены эпюры интенсивности дсвиатора напряжений и максимального касательного напряжения в области, прилегающей к слою со стороны плоскости, на первых 6 элементах для случая достижения предела текучести на первом элементе. Непрерывные линии соответствуют интенсивности дсвиатора напряжений (критерий Губера—Мизеса), а штриховые — максимальному касательному напряжению (критерий Трсска-Ссн-Вснана). Кривые 1 и 3 определяют характеристики напряжений в слое взаимодействия,

2 и 4 — в области, прилегающей к слою со стороны плоскости. К — критерий перехода в пластическое состояние (по интенсивности или максимальному касательному напряжению), тк — соответствующий данному критерию предел текучести.

<Тп

Л

Тт

Рис. 6. Распределение интенсивности девиатора напряжений и максимального касательного напряжения в слое и смежной ему области при достижении предела текучести на первом элементе

Из рисунка видно, что интенсивность напряжений в слое превышает интенсивность напряжений в приграничной области. Следовательно, развитие пластической области (по крайней мерс, се начальная стадия) будет проходить в пределах слоя взаимодействия. В этом случае напряженное состояние слоя в зоне пластического деформирования может ассоциироваться с силами сцепления [15]. Однако в данной модели, их распределение получается из

решения конкретной краевой задачи механики сплошной среды, а не вводится априорно. Как показывают расчеты, величина напряжения отрыва для случая плоской деформации в окрестности конца трещины может существенно превышать предел текучести (см. рис. 4), что является следствием существенной гидростатической составляющей тензора напряжений. Данный результат является известным экспериментальным фактом [16, 17], но его описание невозможно в рамках классического представления трещины Гриффитса, где в концевой точке <тц = и22 = ос и подавление соответствующей сингулярности силами сцепления возможно при конкретном задании закона их распределения, а не получить его из решения конкретной задачи.

Список литературы

1. Черепанов Г.П. Механика хрупкого разрушения. М.: Наука, 1974. 640 с.

2. Панасюк В. В. Механика квазихрункого разрушения материалов. Киев: Наук, думка, 1991. 416 с.

3. Гольдштейн Р.В., Перелъмутер М.Н. Рост трещин но границе соединения материалов // В кн.: Проблемы механики. Сб. статей. М.: ФИЗМАТЛИТ. 2003. С. 221-239.

4. Клевцов Г.В., Ботвина Л.Р. Микро- и макрозона пластической деформации как критерии предельного состояния материала при разрушении // Проблемы прочности. 1984. № 4. С. 24-28.

5. Леонов М.Я., Панасюк В.В. Развитие мельчайших трещин в твердом теле // Прикладная механика. 1959. Т. 5, № 4. С. 391-401.

6. Dugdale D.S. Yielding of steel sheets containing slits // J. Mech. and Phys. Solids. 1960. V. 8, № 2. P. 100-108.

7. Новожилов В.В., Рыбакина О.В. О перспективах построения критерия прочности при сложном нагружении // Прочность при малом числе циклов нагружения. М.: Наука. 1969. С. 71-80.

8. Глаголев В.В., Маркин А.А. Определение термомеханических характеристик процесса разделения // Изв. РАН. Механика твердого тела. 2007. № 6. С. 101-112.

9. Глаголев В.В., Маркин А.А. Оценка толщины слоя взаимодействия как универсального параметра материала // Изв. РАН. Механика твердого тела. 2006. Ш 5. С. 194-203.

10. Гаврилкина М.В., Глаголев В.В., Маркин А.А. К решению одной задачи механики разрушения // ПМТФ. 2007. № 4. С. 121-127.

11. Лурье А.И. Теория упругости. М.: Наука, 1970. 939 с.

12. Крауч С., Старфилд А. Методы граничных элементов в механике твердого тела. М.: Мир, 1987. 328 с.

13. Ишнинский А.Ю., Ивлев Д.Д. Математическая теория пластичности. М: ФИЗМАТЛИТ, 2001. 701 с.

14. Новожилов В. В. О необходимом и достаточном критерии хрупкой прочности // ПММ. 1969. № 2. С. 212-222.

15. Баренблатт Г.И., Христианович С.А. О модуле сцепления в теории трещин // Инженерный журнал. Механика твердого тела. 1968. № 2. С. 69-75.

16. Нотт Дж.Ф. Основы механики разрушения. М: Металлургия, 1978. 256 с.

17. Irvin G.R. Linear fracture mechanics, fracture transition, and fracture control // Engn. Fracture Mechanics. 1968. V. 1. P. 241-257.

Поступило 06.03.2009

Гаврилкина Мария Владимировна (gavrilkina-mv@ramblcr.ru), аспирант, кафедра математического моделирования, Тульский государственный университет.

Глаголев Вадим Вадимович (vadim@tsu.tula.ru), д.ф.-м.н., профессор, кафедра математического моделирования, Тульский государственный университет.

Model of the strain state of the crack’s vicinity

M.V. Gavrilkina, V.V. Glagolev

Abstract. An approach to the resolution of fracture mechanics problems is offered. The fraction model in this case includes a crack in the form of a physical cut and a material laying on an imaginary continuation of a physical cut in the continuum. The given material is defined as an interaction layer and its state of strain forms boundary conditions for adjacent environment. Possible existence of plastic area within the limits of the given layer is supposed. In the given work influence of possible material’s hardening on development of thin plastic area in an interaction layer is investigated. The Gubcr-Von Miscs plastic-yicld criterion is accepted as the condition of the material transition from the clastic state to the plastic one. The case when the layer’s material is described in the perfectly clastic-plastic model is considered.

Keywords: characteristic dimension, interaction layer, thin plastic zone, perfectly clastic-plastic model, material hardening, the Gubcr-Von Miscs plastic-yicld criterion.

Gavrilkina Maria (gavrilkina-mv@ramblcr.ru), postgraduate student, department of mathematical modeling, Tula State University.

Glagolev Vadim (vadim@tsu.tula.ru), doctor of physical and mathematical sciences, professor, department of mathematical modeling, Tula State University.