Атомистический подход к определению зависимости от температуры потенциальной энергии и равновесных межатомных расстояний для монокристаллических образцов с ГПУ-решеткой Текст научной статьи по специальности «Физика»

УДК538.911, 539.32

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

© 2012 И.Ю. Зубко, В.И . Кочуров, М.В. Симонов

Пермский национальный исследовательский политехнический университет

Поступила в редакцию 11.12.2012

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

Ключевые слова: ГПУ-монокристалл, метод атомарной статики, равновесные параметры решетки, температура

Учет температуры в дискретном подходе.

Одним из способов теоретического исследования термомеханических свойств конденсированных сред является дискретный подход, основанный на прямом моделировании движения частиц, вызванного их взаимодействием и приложенными внешними воздействиями. Наиболее распространенным является метод молекулярной динамики (МД), при использовании которого для каждой частицы (атома или молекулы) численно решаются уравнения движения (уравнения баланса количества движения и момента количества движения), содержащие силы и моменты взаимодействия выбранной частицы со всеми остальными частицами тела или с частицами из некоторой ее окрестности [4, 5]. Силы взаимодействия частиц вводятся с помощью потенциалов межатомного взаимодействия, которые являются ядром дискретных подходов и представляют собой приближенный способ описания взаимодействия частиц материала, качественно отражающий основные свойства атомов отталкиваться на малых и притягиваться на больших расстояниях. Получаемые количественные результаты зависят от числовых значений параметров потенциалов. Метод МД хорошо работает для тел малых размеров на малых интервалах физического времени. Ограничения на размеры моделируемых объектов накладываются возможностями вычислительной техники и применяемых алгоритмов для численного решения систем дифференциальных уравнений движения частиц. Интенсивно развиваемые высокопроизводительные методы

Зубко Иван Юрьевич, кандидат физико-математических наук, доцент кафедры «Механика композиционных материалов и конструкций». E-mail: [email protected] Кочуров Владимир Игоревич, старший преподаватель кафедры «Математическое моделирование систем и процессов». E-mail: [email protected] Симонов Максим Владимирович, студент

«распараллеливания» расчетов на многопроцессорных вычислительных устройствах расширяют границы применимости метода МД при исследовании больших систем.

В работе рассматривается предложенная ранее модификация [1, 2] метода частиц, позволяющая провести идентификацию параметров потенциала межатомного взаимодействия для монокристаллов с кубической решеткой (ГЦК- и ОЦК-типа) по их макроскопическим параметрам - равновесному межатомному расстоянию и упругому модулю сдвига. В основе этой методики лежит рассмотрение статики взаимодействующих частиц при явном задании структуры кристаллической решетки. Возможности подхода были продемонстрированы на примере использования потенциала Леннарда-Джонса. Это обосновано тем, что при расчете упругих модулей нет необходимости рассматривать динамические процессы, идущие при сверхнизких температурах или при скоростном деформировании, когда важную роль во взаимодействии атомов могут играть квантовые эффекты, учитываемые в потенциалах специального вида. Использование статической постановки позволяет получить аналитическое решение задачи определения равновесного межатомного расстояния для различных кристаллических решеток и образцов различного размера. При этом точные значения межатомного расстояния получаются для объемов материала с небольшим числом атомов п на ребре образца (от 5 до 50). Для получения макроскопических значений параметров по этим точным решениям можно сделать предельный переход, устремив число атомов п к бесконечности. В частности, это позволяет идентифицировать параметры потенциала по известным макроскопическим физико-механическим свойствам материалов. Отметим, что упругие модули, получаемые при использовании

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

В [1-2] не обсуждался вопрос об учете температуры при расчете упругих модулей монокристаллических твердых тел. В данной работе ставится вопрос о способе задания температуры в статическом подходе. При моделировании поведения механических свойств различных тел методом молекулярной динамики температура вводится через кинетическую энергию движения атомов относительно среднего движения тела. При стабилизации системы энергия крупномасштабного движения атомов переходит в мелкомасштабные колебания с меньшей амплитудой. Эти мелкомасштабные колебания относительно неменяющихся со временем центров отождествляются с тепловым движением и с помощью выражения Е=3кТ, где Е - амплитуда пульсаций кинетической энергии, к - постоянная Больцмана, Т - температура тела. При построении равновесных кривых зависимости физико--механических свойств исследуемых материалов (например, коэффициента теплового расширения, внутренней энергии, упругих модулей) от температуры при динамическом моделировании после подвода к системе частиц энергии, вызывающей повышение температуры, необходимо проводить расчеты по переходу системы в равновесное состояние, и затем определять значение требуемого свойства в равновесном состоянии. При этом желательно для каждой точки получаемой зависимости проводить осреднение по достаточной для статистической обработки выборке реализаций рассматриваемой системы частиц. В статическом подходе нет возможности вводить температуру как характеристику теплового движения атомов, поскольку само движение при этом исключается. Можно было бы считать, что статический подход дает возможность определять равновесные положения центров, относительно которых атомы совершают колебательное движение, независимое от перемещений этих центров, связанных с механическим движением как жесткого целого или деформированием тела. Тогда температуру в статическом подходе можно учесть с помощью введения некоторой искусственной дополнительной силы, приводящей к увеличению при нагреве межатомного расстояния в равновесной конфигурации, свободной от внешних механических воздействий.

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

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

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

б)

Рис. 1. Виды распределений 5000 точек по единичной сфере, полученные: а) при равномерном разбиении интервала для угла 0 е [-% /2; % /2], б) при равномерном разбиении интервала [-1;1] значений функции sm0

Алгоритм получения равномерного распределения точек по сфере. В сферической системе координат равномерное разбиение интервалов для угла фе[0;2л), откладываемого от оси 0x1 в

плоскости х10х2 и угла 9 е [-п /2; п /2], откладываемого от плоскости х]0х2, приводит к концентрации точек на полюсах сферы (рис. 1а). Связь декартовых и сферических координат имеет вид: XI = г cos9 cosф, х2=г cos9 sinф, х3 = г sin9. Для того чтобы не происходило концентрации точек вблизи полюсов сферы 9 = ±л/2, разбивать необходимо не интервал для угла 9е[-п/2; п/2], а интервал [-1;1] значений функции sin9. Функция cos9 (положительная для рассматриваемого интервала значений угла) выражается через полученное случайное значение как -\Д - ят2 9 . В результате получается равномерное покрытие сферы точками (рис. 1б).

Для проверки равномерности распределения точек по сфере проводится разбиение сферы на области равной площади, совокупность которых полностью покрывает сферу. С помощью вершин вписанного в единичную сферу икосаэдра, соединенных геодезическими линиями, она покрывается 20 одинаковыми «равносторонними» криволинейными треугольниками. Для исследования статистических свойств получаемых случайных распределений икосаэдр случайным образом поворачивается относительно сферы и при каждом новом его положении для одного и того же распределения определяется число точек, попавших в каждый криволинейный треугольник. Это позволяет получить представительную выборку попаданий точек на сфере в криволинейные треугольники, достаточную для анализа статистических характеристик распределения точек на сфере. Анализ первого алгоритма распределения показывает наличие треугольников с большим числом точек, которые соответствуют областям вблизи полюсов сферы (рис. 1а). При использовании второго алгоритма получается, что все накладываемые на сферу треугольники примерно одинаково заполняются точками распределения. Соответствующие виды функций распределения вероятностей показаны на рис. 2, который демонстрирует узкую локализацию при использовании второго алгоритма. Функция, соответствующая первому алгоритму имеет длинный «хвост».

Таким образом, второй алгоритм, используемый для размещения точек по сфере, действительно позволяет получить равномерное распределение направления смещений атомов в пространстве.

Для оценки преимущественного направления в распределении, аналогичном рис. 1а, используется тензор Р = , Р, - единичный вектор,

I

соответствующий выбранному из данного распределения направлению, причем ^р, ~ 0. Получен-

I

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

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

Рис. 2. Распределение вероятности для двух способов размещения 5000 точек по треугольникам: 1) гистограмма с «длинным хвостом» соответствует рис. 1а; 2) локализованная гистограмма соответствует рис. 1 б

б)

Рис. 3. Геометрия образца с ГПУ-решеткой, атомы слоя А изображены серым цветом, атомы слоя В -коричневым цветом: а) вид спереди, б) вид сверху; п=6

Влияние температуры на свойства ГПУ-монокристалла. При исследовании зависимости равновесных межатомных расстояний ГПУ-решетки от амплитуды возмущений рассматривается монокристалл в форме прямой призмы с правильным треугольником в основании (рис. 3), симметрия которого соответствует симметрии ГПУ-решетки. Такой выбор формы образца обоснован в работе [3]. Для каждого заданного значения амплитуды A при произвольных параметрах решетки a, Ь и произвольных параметрах а, р потенциала Лен-нарда-Джонса:

Ф(г) = р((а/г)12 -2(а/г)6)

(1)

вычисляется полная потенциальная энергия атомов образца.

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

5U = 5A + 5Q §A ^ ^-8ц ds + ^-8и dv

S = 81 и ^ 8Ч = 0

(2)

значений a и Ь (то есть объема), 5S=0, и приток тепла также не зависит от варьируемых параметров решетки в равновесном состоянии: 5Q=0.

Будем считать, что в состоянии равновесия энергия равномерно распределена по всем атомам кристалла. Таким образом, для определения равновесных параметров a и Ь кристаллической решетки при заданной температуре (амплитуде и частоте тепловых колебаний атомов) система должна удовлетворять условию 5^0, и = шКй, где и -

массовая плотность внутренней энергии образца, осредненная по значительному промежутку времени, N - количество атомов, m - масса атома. Массовая плотность внутренней энергии u складывается из плотности кинетической энергии теплового движения атомов и плотности потенциальной энергии их взаимодействия, то есть в каждый момент времени она равна выражению

и = ]^у,2/2К + £ £ ф(г. - г. )/шК

¿=1 j=i+1

(3)

где 5U - вариация внутренней энергии образца, 5A - вариация работы внешних поверхностных ^ и объемных рГ сил, 5Q - изменение тепла в рассматриваемом состоянии равновесия, 5и - кинематически допустимые вариации перемещений, Sl - часть поверхности тела, на которой задаются статические граничные условия, S2 - часть поверхности, на которой задаются кинематические граничные условия, S - полная (внешняя) поверхность тела. Приток тепла к телу определяется интегралом

8Q = J р Т 8Sdv, где T - температура, S - энтропия

V

системы. Объемные силы в рассматриваемом случае на образец не действуют. Поверхностные силы в зависимости от постановки задачи могут быть заданы любыми, например, соответствовать гидростатическому давлению при исследовании его влияния на кристалл. В данном случае поверхностные силы равны нулю, поэтому справедливо 5A=0.

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

где Vi - скорость 7-го атома, ri - его радиус-вектор, ф - парный потенциал межатомного взаимодействия. Осредняя записанную сумму по большому периоду времени т и принимая, что каждый атом в состоянии термодинамического равновесия совершает близкие к гармоническим колебания с заданной амплитудой A и частотой ю, получим

т N-1 N

и = 3А2ю2 / 2 + т-1£ ф(г (1) - г (1)) / шNdt

(4)

0 ¿=1 .=¿+1

Здесь учтено, что кинетическая энергия тепловых колебаний одинаково распределена по трем взаимно ортогональным направлениям в пространстве. Осреднение по времени в последнем слагаемом в дискретном подходе заменяется осреднением по некоторому большому числу реализаций теплового возмущенного состояния образца. Заметим, что среднее значение потенциальной энергии взаимодействия атомов не равно средней кинетической энергии, как для случая классических осцилляторов. Минимальный уровень потенциальной энергии выражается через параметр р потенциала Леннард-Джонса, который является независимой постоянной для каждого материала.

В отличие от объектов статистической физики рассматриваемая система не является достаточно большой для применения методов статистического анализа и подвергается непосредственной обработке. Реализации возмущенного состояния кристалла генерируются на компьютере, и в каждой из них с помощью потенциала Леннард-Джонса вычисляется потенциальная энергия взаимодействия всех атомов кристалла. В силу выпуклости потенциальной энергии кристаллического образца требуется найти ее минимум по двум переменным параметрам решетки a и Ь в каждой полученной реализации. В работе для решения задачи минимизации использовался метод Хука-

Дживса. Зависимость равновесных параметров кристаллической решетки и плотности потенциальной энергии от амплитуды тепловых колебаний A получается осреднением найденных параметров решетки и потенциальной энергии образца по всем реализациям при каждом заданном значении амплитуды тепловых колебаний.

Для каждого значения амплитуды A было получено по 1000 реализаций случайных конфигураций кристаллической решетки, для каждой из которых определены равновесные межатомные расстояния, и по этим 1000 испытаниям определены медианы распределений значений параметров решетки. Также в каждой реализации определялась

N-1 N

Ф = ^^ф(Т -Tj)/N потенциальная энергия всей

i=i j=i+i

системы атомов, отнесенная к числу атомов образца. Осреднением по всем реализациям определялось значение ф.

Для различного числа атомов на ребре образца, как и в работе [3], получены различные значения безразмерных равновесных межатомных расстояний а/а, b/а и их отношения b/a, а также найдены их зависимости от безразмерной амплитуды A/а наложенных тепловых возмущений (рис. 4а, б). Наименьший образец имел 4 атома на ребре основания. Кривые для параметров решетки монотонно возрастают с ростом амплитуды тепловых колебаний, то есть предложенный подход качественно верно описывает явление теплового расширения кристалла. Все кривые для зависимости параметра а/а от температуры пересекаются при значении А/а«0,2. При этом же значении пересекаются все кривые для зависимости параметра b/a от температуры. Для получения предельной кривой, соответствующей макроскопическому телу (при n^-да), для каждого фиксированного значения амплитуды возмущений строилась гиперболическая аппроксимация зависимости значений параметров I ПУ-решетки от числа атомов на ребре образца. Заметим, что предельная кривая для отношения параметров решетки b/a показывает, что оно не зависит от температуры. Для наночастиц с ГПУ-решеткой отношение b/a начинает уменьшаться после достижения относительной амплитудой теплового движения значения А/а«0,3. С увеличением размеров образца отношение b/a перестает зависеть от температуры.

а/cr, Ь/а

_1_i_i_._i_1_i_i_i_i_I_i_i_._i_1_._i_i_i_1_i_i_._i_l_ ill ffg

0.0 0.1 0.2 0.3 0.4 0.5 a)

b/a

1.4 -1.3 -1.2 -1.1 ;

.:---------_ А,а

0.0 0.1 0.2 0.3 0.4 0.5 б)

Рис. 4. Зависимости от амплитуды тепловых возмущений при различных п: а) параметров решетки а и Ь, б) их отношения Ь/а. Пунктирная линия -предельная кривая при п^-да

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

а/а, Ь/а

.........................А/а

0.0 0.1 0.2 0.3 0.4 0.5

а)

Ь/а

1.4 -1.3 -1.2 -1.1 ;

-i————— А/а 0.0 0.1 0.2 0.3 0.4 0.5

б)

Рис. 4. Зависимости от амплитуды тепловых возмущений при различных n: а) параметров решетки a и b, б) их отношения b/a. Пунктирная линия -предельная кривая при n^-да

Для установления связи между температурой и амплитудой теплового движения атомов в простейшем случае можно принять, что частота колебаний при тепловом движении атомов не зависит от амплитуды. Принимая, как и выше, что атомы кристалла совершают близкие к гармоническим колебания с амплитудой А и частотой ю, получим осредненное по периоду колебаний выражение 3шА2ю2/2 для кинетической энергии отдельного атома. Приравнивая эту величину энергии 3кТ/2 теплового движения атома идеального газа (или используя «газовый термометр»), получим равенство шА2ю2, то есть температура Т рассматриваемой системы при фиксированном значении

амплитуды А может быть оценена как Т=шА2ю2/к. Принимая, что амплитуда возмущений Ас соответствует некоторой критической температуре Тс, получим оценку частоты колебаний атомов го = Ас -1^/кТс /т. Вопрос выбора значения Тс остается открытым.

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

= д(тй) _ а(3кТ/2 + Ф(А))=зк /2 аф(А) аА

р ат ~ +"

ат

аА ат

Учитывая, что Т=шА2ю2/к, получим

аф(А)

С =3к/2+-

р аА

ат

аф(А) к

— = 3к / 2 + - ,

аА I аА 2тАю2

Вводя безразмерную удельную теплоемкость С*у, окончательно получим:

р (А/ а) а

С у - (Ср/к - 3/2)

2т а2ю2

Р а(А/а) А

где ф* - ф / Р - безразмерная величина, соответствующая осредненной потенциальной энергии ф.

-1

Зависимость введенной теплоемкости от амплитуды тепловых колебаний имеет максимум вблизи точки А/а«0,3, а затем начинает уменьшаться (рис. 5б).

Ф/Р

б)

Рис. 5. Зависимость от амплитуды тепловых возмущений при различных п: а) удельной потенциальной энергии кристалла, б) его теплоемкости. Пунктирная линия - предельная кривая при п^-да

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

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

Работа выполнена при финансовой поддержке РФФИ (грант №11-01-96033-р-Урал-а).

СПИСОК ЛИТЕРАТУРЫ:

1. Зубко, И.Ю. Определение упругих постоянных ГЦК-монокристаллов с помощью потенциала межатомного взаимодействия / И.Ю. Зубко, П.В. Трусов // Вестник ПНИПУ. Механика. 2011. № 1. С. 147-169.

3.

Зубко, И.Ю. Вывод упругого закона монокристаллов металлов из потенциала межатомного взаимодействия / И.Ю. Зубко, О.В. Мелентьева, В.П. Морозова, В.И. Кочуров // Вестник Нижегородского университета им. Н.И. Лобачевского. 2011. №4. Ч.5. С. 2181-2183. Симонов, М.В. Определение равновесных параметров решетки различных ГПУ-монокристаллов с помощью потенциала межатомного взаимодействия Ми / М.В.

Симонов, И.Ю. Зубко // Вестник ПНИПУ. Механика. 2012. № 3. С. 205-218.

Кривцов, АМ. Деформирование и разрушение твердых тел с микроструктурой. - М.: Физматлит, 2007. 304 с.

Метод молекулярной динамики в физической химии. Под ред. Ю.К Товбина. - М.: Наука, 1996. 334 с.

ATOMISTIC APPROACH TO FINDING THE DEPENDENCE OF POTENTIAL ENERGY AND EQUILIBRIUM LATTICE PARAMETERS FOR HCP-MONOCRYSTALS ON TEMPERATURE

5

© 2012 I.Yu. Zubko, V.I. Kochurov, M.V. Simonov Perm National Research Polytechnical University

The method of specifying temperature for estimating the physical properties of monocrystals using the atomistic statics approach was considered. Thermal oscillations of atoms in the crystal lattice are emulating by imposed perturbations on every atom positions with prescribed amplitude A. The directions of the perturbations are homogeneously distributed in space. The frequency of atoms oscillations is considered as independent parameter of model that has to be identified as well as interatomic interaction potential parameters for each chosen material. The dependencies of potential energy per atom and two parameters of crystal lattice spacing on the perturbation amplitude A were found for HCP-monocrystals of several sizes.

Key words: HCP-monoctystal, atomistic statics method, equilibrium lattice parameters, temperature

Ivan Zubko, Candidate of Physics and Mathematics, Associate Professor at the Department "Mechanics of Composite Materials and Constructions". E-mail: [email protected]

Vladimir Kochurov, Senior Teacher at the Department "MathematicalModeling of Systems and Processes". E-mail: [email protected] Maksim Simonov, Student