Метод решения краевой задачи релаксации остаточных напряжений в упрочненном слое цилиндрического образца при виброползучести Текст научной статьи по специальности «Физика»

УДК 539.376+621.787

МЕТОД РЕШЕНИЯ КРАЕВОЙ ЗАДАЧИ РЕЛАКСАЦИИ ОСТАТОЧНЫХ НАПРЯЖЕНИЙ В УПРОЧНЕННОМ СЛОЕ ЦИЛИНДРИЧЕСКОГО ОБРАЗЦА ПРИ ВИБРОПОЛЗУЧЕСТИ

М. Н. Саушкин, Е. В. Дубовова

Самарский государственный технический университет,

443100, Самара, ул. Молодогвардейская, 244.

E-mail: msaushkin@gmail.com

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

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

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

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

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

Михаил Николаевич Саушкин (к.ф.-м.н., доц.), докторант, каф. прикладной математики и информатики. Елена Валерьяновна Дубовова, аспирант, каф. прикладной математики и информатики.

- восстановление полей остаточных напряжений и деформаций по схеме сложного напряжённого состояния после применения процедуры упрочнения для цилиндрического образца;

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

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

1. Решение задачи восстановления полей остаточных напряжений и деформаций в цилиндрическом образце. Введём стандартную цилиндрическую систему координат г, в, z и обозначим через o^es, ^es, ^Zes окружную, радиальную и осевую компоненты тензора остаточных напряжений (соответственно) в цилиндрическом образце. Через qe, qr, qz обозначим пластические деформации, возникающие в цилиндре после процедуры ППД и инициирующие остаточные напряжения.

Касательные компоненты тензоров остаточных напряжений и деформаций считаются малыми, что справедливо для ряда технологических режимов упрочнения (гидродробеструйная обработка, бомбардировка микрошариками и другие), поэтому ими можно пренебречь Для нормальных компонент пластических деформаций, наведённых в результате упрочнения, выполняется соотношение

qe (r) = qz (r), (1)

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

Будем считать, что экспериментально определена окружная компонента 00es(r) и поставленная задача сводится к вычислению напряжений a^es(r)

и ^Zes(r), а также остаточных пластических деформаций с помощью инфор-

мации об измеренной компоненте 00es(r) [1-3].

Из уравнения равновесия

dnres

Л* r I res res (сл\

г — + «7,. = a, (2)

и условия ^res(a) = 0 (а — радиус цилиндрического образца), которое означает, что после упрочнения образец находится в естественном ненагруженном состоянии, нетрудно выразить через ^ees:

1 Г

оГ(г) = - aT{z)dz. (3)

r J о

Компонента ^Zes(r) может быть определена лишь с учётом остаточных пластических деформаций.

Пусть главные компоненты тензора полной деформации цилиндрического образца е° (i = r, в, z), приобретённые в результате упрочняющей обработки, представлены в виде

£°(r) = e°(r)+ qi(r), (4)

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

Из условия несжимаемости при пластическом деформировании qz + q$ + + qr = 0, соотношений (1), (4) легко установить, что

qr

qe = qz = -у, (5)

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

'£+*>=*-4-* (б>

Входящие в (6) упругие деформации нетрудно выразить через остаточные напряжения из закона Гука:

„О _ 1 [„.res „.res | _resM „О __ W res ../'„.res I „.res'll

er ~ Elar — ^\ав + az ee — Eiae — ^\ar + az )J (‘)

(p — коэффициент Пуассона, E — модуль Юнга).

В соотношениях (7) наряду с уже известными компонентами o^es и входит неизвестная компонента ^Zes. Для определения ^Zes требуется ввести гипотезу плоских сечений для цилиндрического образца, согласно которой поперечные сечения цилиндрического образца до упрочнения остаются плоскими и после упрочнения. При этом очевидно e°(r) + qz(r) = £* (e*z = const, r € [0; aj). Гипотеза может нарушаться лишь вблизи свободных торцов цилиндрического образца.

Выражая упругую осевую деформацию через закон Гука и подставляя её в формулу полных осевых деформаций, найдём

О) = E К - qz (r)) + p(<7“ (r) + fff(r)). (8)

С помощью полученного соотношения в законе Гука (7) можно исключить

CTres:

z

е° = К1 ~ ~ ^аТ] ~ Qz),

ев = Ч^ К1 “ ^)аТ ~ - К4 ~ Qz) ■

C учётом (Б) и (9) уравнение (б) можно записать в виде

(9)

где

r^ + IT^'Mr) = 9(r)' (1С1)

_<8(г)-(ТГ(г) I-L ^(laTir) ^<кт™(г) і ^

EE

dr dr

Решение дифференциального уравнения (10) с учётом граничного условия ііш qQ(г) = 0 имеет следующий вид:

r^r 0

qe{r) =--1 2/Х^ Г z& [аГ(^ + 2oT(z)] dz-

Е{ 1 + ц)г-'О

-^[а-дугм-^гм]- (12)

Теперь можно полностью восстановить поля остаточных пластических деформаций: согласно (12) вычисляется qQ, а затем в соответствии с (5) — дг и qz.

Для определения компоненты ст^^г), согласно (8), достаточно найти е*. Это можно сделать исходя из условия нулевого суммарного осевого усилия, действующего на образец:

/ гст^г) = 0.

Jo

Из последнего условия и (8) находится выражение для е*:

£*=^£г (,'<г) ■ |[<7“(г)+<7“(г)]) *■ (1з>

Вычислив согласно (13) величину е*, можно в соответствии с (8) однозначно определить распределение компоненты 0^68(г).

2. Методика расчета процесса релаксации напряжений в упрочнённом слое цилиндрического изделия при ползучести. В настоящей работе предлагается метод расчёта процесса релаксации остаточных напряжений в поверхностно упрочненном слое цилиндрического образца при продольной растягивающей нагрузке ст°(£), на которую наложена циклическая составляющая при постоянной температуре в условиях виброползучести, то есть когда амплитудная величина ста составляет 1-3% от статической компоненты ст°(£).

Для построения модели релаксации остаточных напряжений введём в рассмотрение следующие гипотезы и допущения [1].

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

В работах [1, 2] показано, что кинетику напряжённо-деформированного состояния в поверхностном слое можно рассматривать исходя из идей декомпозиции и агрегирования. При этом расчёт проводится исходя из предположений, что упрочнённый слой не оказывает существенного влияния на жёсткость и деформируемость самого цилиндра и что он (поверхностный слой) деформируется в режиме «жёсткого» нагружения при заданных значениях компонент полных деформаций ег (£), ег (£) и е$ (£) на поверхности цилиндра, определяемых одноосной виброползучестью цилиндрического образца при действии растягивающей нагрузки без учёта поверхностного слоя.

На основании вышеизложенного компоненты полной деформации в поверхностно-упрочненном слое можно представить в виде

е*(*) + е° = %(г)+ еГМ) + РГМ), г = г,0,*, (14)

где еД£) —полная деформация цилиндрического образца, рассчитываемая с помощью определяющих реологических уравнений для одноосного напряженного состояния; е° — величина полных остаточных осевых деформаций

после процедуры поверхностного пластического деформирования, определяемая соотношением (4); qi(r) —компоненты остаточных пластических деформаций, определяемые согласно методике п. 1; ej;es(r, t) —компоненты упругих деформаций и pjjes(r, t) —компоненты деформации ползучести в упрочнённом слое, рассчитываемые согласно реологической модели для сложного напряженного состояния по заданным компонентам ez (t), er (t) и ez (t).

Для вычисления осевой компоненты деформации ez (t) можно использовать любую одноосную модель неупругого деформирования в условиях виброползучести при растягивающей нагрузке. При этом поперечные деформации er (t) и eg (t) цилиндрического образца можно будет рассчитать исходя из пуассоновского сужения материала при растяжении.

Для вычисления кинетики остаточных напряжений при релаксации из (14) выразим

eres (r, t) = ei (t) + e0 (r) - qi(r) - pres(r, t), i = r, 0,z и из закона Гука получим

<s (r, t) - ^es(r, t) - ^ges(r, t)

-^es(r, t) + ^Zes(r, t) - ^ges(r, t)

-^es(r, t) - ^Zes(r, t) + CTges(r, t)

откуда

*) = ifeef* (r, *) + (1+/jff-2/i) (еГ(^ t) + еГ(г, t) + ef%r, t)). (15)

Соотношения (15) позволяют следить за процессом релаксации остаточных напряжений в упрочненном слое при неупругом реологическом растяжении цилиндрического образца. Для расчёта неупругих деформаций может быть применена любая модель реологического деформирования материала при виброползучести.

3. Математическая модель виброползучести. Для решения краевой задачи о деформировании цилиндрического образца в условиях виброползучести (определение компонент er (t), ее (t) и ez (t)) и расчёта процесса релаксации остаточных напряжений в настоящей работе используется модель виброползучести энергетического типа [4,5], на основании которой полная деформация ez (t) в одноосном случае представляется в виде

ez (t) = ez (t) + pz (t), (16)

где ez(t) — полная деформация, ez(t) — упругая деформация, pz(t) — деформация ползучести. Причём деформация ползучести при ст0 = const и = const описывается в следующем виде:

P(t) = Uk (t) + Vk (t) + w(t);

k k Uk(t) = Ak [ak(a(t)/a*)n2 - Uk(t)];

Vk(t) = { Ak^ если ТЦ(t) > 0, Л(t) = bk(CT(t)/CT*)n2 - vk(t);

= eres (r, t), = ezes (r, t), = eges (r, t),

w(t) = c(a(t)/CT* )mi;

u(t) = a(<T°)<r(t)p(t) + gi(a0, a°, f)^rN; CT(t) = (1 + w(t))a0; ста (t) = (1 + w(t))CT0.

Здесь и(Ь), ^(Ь), эд(Ь) —вязкоупругая, вязкопластическая и вязкая составляющие деформации ползучести; ст(Ь) и ста(Ь) —соответственно статическая и амплитудная компоненты истинного напряжения; ст° и ст° — статическая и амплитудная компоненты номинального напряжения; Л&, а^, 6^, с, п, т, ст* —реологические константы материала, при помощи которых описываются первая и вторая стадии ползучести и обратимая часть деформации ползучести; ш(£) —параметр повреждённости материала; а(ст°) д1(ст°,ст°,/) —параметры материала, контролирующие процессы разупрочнения (описывающие третью стадию ползучести); N — число циклов в реализации; / — частота изменения ста.

Критерий разрушения материала имеет вид

где А**, А** — критические величины работ напряжений в условиях стационарной ползучести и виброползучести соответственно, при которых в момент времени Ь* происходит разрушение материала.

Модель (16), (17) может быть распространена и на случай сложного напряжённого состояния [4,5], в условиях которого рассчитывается релаксация остаточных напряжений в упрочнённом слое.

Поперечные деформации цилиндрического образца ег (Ь), ее (Ь), возникающие за счёт пуассоновского сужения материала, рассчитываются на основании формул

где р' и р" — коэффициенты Пуассона для обратимой u и необратимой v компонент деформации ползучести соответственно.

4. Методика идентификации параметров модели виброползучести. Для

идентификации параметров предложенной модели используется серия кривых стационарной ползучести от начала нагружения материала до момента разрушения при ст0 = const и серия кривых циклической ползучести от начала нагружения до момента разрушения при ст0 = const, ст^ = const.

Отметим, что приведённые в п. 3 определяющие соотношения для деформации ползучести при ста = 0 полностью совпадают с соотношениями для стационарной ползучести [1,6]. Поэтому, считая константы модели, при помощи которых описываются первая и вторая стадии ползучести, известными, необходимо определить параметры модели, описывающие процессы разупрочнения в условиях стационарной и циклической ползучести.

Сначала определяются параметры а и A* для стационарной ползучести. Из экспериментальных данных по серии кривых ползучести материала для заданного напряжения ст0 = const определяются время разрушения t* материала и величина неупругой деформации р* в момент разрушения. Для

(17)

ei(t) = -рст0/Е - p'uz(t) - p''vz(t) - wz(t)/2, i = r, 0, (18)

каждой экспериментальной кривой (пары значений £*, р*) целенаправленно варьируется параметр а и осуществляется численный расчёт по модели (16), (17) до достижения минимума функционала

где рР, £р — расчётные значения по модели (16), (17), полученные при заданном параметре а.

Используя значение параметра а, при котором выполняется (19), вычисляем критическую величину работы напряжений на деформациях ползучести:

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

В случае, когда а и А* для материала являются константами, задача их определения сводится к минимизации функционала отклонения расчётных значений рР и £Р от соответствующих экспериментальных значений р*, £*:

где M — количество экспериментальных кривых стационарной ползучести.

Далее, считая, что все параметры модели для стационарной ползучести известны, определяем параметры gi и A** для виброползучести. Как и в предыдущем случае, из кривых виброползучести для ст0 = const, СТ0 = const определяется пара величин t* и р*. Варьированием параметра gi для каждого ст0 проводится минимизация функционала аналогичного (19) и определяется значение этого параметра. Отметим, что при варьировании gi в модели (16), (17) в качестве а берётся значение, которое найдено именно при минимизации функционала (19) при заданном ст0 и ст0 = 0. В итоге находится величина критической работы напряжений на деформации циклической ползучести:

Считая параметры $1 и А* для материала константами, их можно взять как средние значения полученных величин при заданных ст^ с использованной выше методики.

Алгоритм минимизации функционала (19) и определения значений параметров а и $1 для каждой экспериментальной кривой приведён в [7].

(19)

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

5. Расчёт релаксации остаточных напряжений в цилиндрическом образце при виброползучести. Численная реализация предложенного в п. 3 метода проводилась на модельном примере исходного распределения полей остаточных напряжений в поверхностном слое цилиндрического образца радиуса а = 3,76 мм из сплава ЭП 742 при ползучести с температурой Т = 650С при следующих режимах вибронагружения: ст0 = 690 МПа, ст° € {0, 50, 75} МПа.

Из методики, изложенной в п. 2, следует, что для моделирования полей остаточных напряжений в упрочнённом слое образца достаточно знать экспериментальные данные для компоненты ст#68 в слое. В работах [1-3] показано, что для аппроксимации этих данных хорошо подходит функция

где ст0, ст1, Ь — параметры аппроксимации.

Моделирование полей остаточных напряжений проводилось при следующих заданных параметрах аппроксимации (20): сто = 815,54 МПа, ст1 = = 15,54 МПа, Ь = 0,08 мм. Для расчёта релаксации остаточных напряжений при вибронагружении при Т = 650 С по серии стационарных кривых и кривых виброползучести по методике п. 4 определены параметры разупрочнения материала для модели (16), (17): а = 0,0195, А* = 98,013; 71 = 3,23 ■ 10-4, А* = 3,656 ■ 103. Остальные параметры модели для этого материала при Т = 650С взяты из работ [1,6]: ст* = 490,5 МПа, к = 2, А1 = 0,022, а1 = = 7,32 ■ 10-4, Ь1 = 5,37 ■ 10-3, А2 = 0,2, а2 = Ь2 = 0, с = 7,22 ■ 10-7, Ш1 = 14,3.

Согласно общей методике решалось две задачи:

1) при заданном нагружении рассчитывалась кинетика компонент полных

деформаций (£), (£), е# (£) по соотношениям (16), (18);

2) по известным из решения первой задачи значениям компонент полных деформаций согласно (15) рассчитывалась кинетика остаточных напря-

Результаты расчётов при указанных нагрузках приведены на рис. 1 и 2. На рис. 1 приведена кинетика компонент ст#68, стГГ68 на поверхности упрочнённого слоя. На рис. 2 приведены эпюры компоненты ст^68 по упрочнённому слою при нагрузке ст0 = 690 МПа, ст0 = 50 МПа в различные моменты времени. Отметим, что на рис. 2 имеется характерный «упругий скачок» компоненты ст^68 в момент приложения нагрузки на величину ст0. Аналогичная кинетика эпюр по толщине слоя наблюдается и у компоненты ст#68, но без «упругого скачка».

Анализ расчётных данных показывает, что амплитудная составляющая ст0 оказывает существенное влияние на процесс релаксации остаточных напряжений — происходит «расслоение» кривых релаксации, что хорошо видно на рис. 1.

Выводы. По результатам работы можно сделать следующие выводы:

1) на основе идей декомпозиции и агрегирования разработан метод расчёта релаксации остаточных напряжений в цилиндрическом образце на фоне

(20)

жений.

Рис. 1. Кривые релаксации остаточных напряжений на поверхности цилиндрического образца из сплава ЭП 742 при T = 450°Спри виброползучести (сплошная линия — ales(a,t); штриховая линия — Uges(a,t)): 1 — a° = 690 МПа, a° = 0; 2 —а° = 690 МПа, а° = 50 МПа; 3 —

а° = 690 МПа, 7° = 75 МПа

Рис. 2. Кинетика эпюр а^ев в упрочнённом слое цилиндрического образца из сплава ЭП 742 при Т = 450°Спри виброползучести (сплошная а0 = 690 МПа, а° = 50 МПа): 1 — Ь = 0 — 0 (эпюра после упрочнения); 2 — Ь = 0 + 0 (эпюра в момент приложения нагрузки); 3 — Ь = 50 ч; 4 —

Ь =100 ч; 5 — Ь = 200 ч

виброползучести самого изделия;

2) построена модель виброползучести материала и методика идентификации её параметров;

3) проведён численный расчёт модельной задачи по релаксации остаточных напряжений в цилиндрическом образце при виброползучести.

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

Работа выполнена при поддержке Федерального агентства по образованию (проект РНП 2.1.1/3397).

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Радченко В. П., Саушкин М.Н. Ползучесть и релаксация остаточных напряжений в упрочнённых конструкциях. — М.: Машиностроение-1, 2005. — 226 с.

2. Радченко В. П., Саушкин М. Н. Математические модели восстановления и релаксации остаточных напряжений в поверхностно упрочнённом слое цилиндрических элементов конструкций при ползучести// Изв. вузов. Машиностроение, 2004. — №11. — С. 3-17.

3. Радченко В. П., Саушкин М. Н. Расчёт релаксации остаточных напряжений в поверхностно упрочнённом слое цилиндрического изделия в условиях ползучести // Вестн.

Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2001. — №12. — C. 61-73.

4. Радченко В. П., Кичаев П. Е., Симонов А. В. Энергетический вариант модели реологического деформирования и разрушения металлов при совместном действии статических и циклических нагрузок // ПмТФ, 2000. — №3. — C. 169-176; англ. пер. Radchenko V. P., Kichaev E. K., Simonov A. V. Energetic version of the model of rheological deformation and destruction of metal under a joint action of static and cyclic loads // J. Appl. Mech. Tech. Phys., 2000. — Vol. 41, No. 3. — P. 531-537.

5. Радченко В. П., Кичаев П. Е., Гагаринский В. С. Уравнения состояния и критерий разрушения материалов в условиях циклической ползучести // Изв. вузов. Машиностроение, 2006. — №5. — C. 9-18.

6. Радченко В. П., Ерёмин Ю. А. Реологическое деформирование и разрушение материалов и элементов конструкций. — М.: Машиностроение-1, 2004. — 256 с.

7. Саушкин М. Н., Дубовова Е. В. Идентификация параметров, контролирующих процессы разупрочнения материала в условиях ползучести и виброползучести / В сб.: Труды Третьей Всероссийской научной конференции (29-31 мая 2008 г.). Часть 1: Математические модели механики, прочности и надёжности элементов конструкций/ Матем. моделирование и краев. задачи. — Самара: СамГТУ, 2008. — C. 266-272.

8. Радченко В. П., Симонов А. В. Разработка автоматизированной системы построения моделей неупругого деформирования металлов на основе методов непараметрического выравнивания экспериментальных данных// Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 1999. — № 7. — C. 51-62.

Поступила в редакцию 01 /II/2010; в окончательном варианте — 10/III/2010.

MSC: 74A10, 74C05

A METHOD OF SOLVING BOUNDARY VALUE PROBLEM OF RESIDUAL STRESSES RELAXATION IN THE HARDENED LAYER OF CYLINDRICAL SPECIMEN DURING VIBROCREEP

M. N. Saushkin, E. V. Dubovova

Samara State Technical University,

244, Molodogvardeyskaya str., Samara, 443100.

E-mail: msaushkin@gmail.com

Presents the method of calculating the relaxation of residual stresses under vibrocreep.

Based on the ideas of decomposition and aggregation of the calculation of the kinetics of relaxation of residual stresses in the surface layer is reduced to gluing solutions of two boundary value problems. In the process of solving the first boundary value problem defined by the stress-strain state of a cylindrical specimen during creep without taking into account the surface hardened layer. The second boundary value problem is investigated relaxation of residual stresses in surface hardening layer, deformable in the “hard” loading at given values of the strain tensor components, which are determined by solving the first boundary value problem. It is shown that the imposition of vibroloads landlines load significantly accelerates the process of relaxation of residual stresses. Results of model calculations are presented.

Key words: vibropcreep, parameter identification, residual stresses, residual stresses relaxation.

Original article submitted 01/II/2010; revision submitted 10/III/2010.

Mikhail N. Saushkin (Ph. D. (Phys. & Math.)), Doctoral Candidate, Dept. of Applied Mathematics & Computer Science. Elena V. Dubovova, Postgraduate Student, Dept. of Applied

Mathematics & Computer Science.