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

ВЕСТНИК ТГГПУ. 2010. №3

ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ

УДК 534.2:532

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

© А.А.Аганин, Д.Ю.Топорков, Т.Ф.Халитова, Н.А.Хисматуллина

В статье рассматриваются малые (осесимметричные) искажения сферической формы кавитационного (парового) пузырька при его сверхсильном сжатии в условиях экспериментов по акустической кавитации дейтерированного ацетона (Taleyarkhan et al, Science. 2002. V.295). Изучается применимость относительно простой и очень экономичной в использовании модели (основанной на расщеплении двумерного движения жидкости и пара на сферическую составляющую и ее малое несферическое возмущение) для оценки максимально возможного роста амплитуды таких искажений. При описании эволюции искажения в этой модели используется ряд допущений, справедливость которых в финальной высокоскоростной стадии сжатия строго не обоснована. Показано, что и при таких допущениях получаемые оценки являются вполне удовлетворительными, что для низкочастотных искажений сферичности подтверждено сравнением с результатами применения двумерной гидродинамической модели.

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

1. Введение

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

В настоящей работе изучается возможность получения оценок максимально возможного нарастания амплитуды малых искажений сферичности пузырька при его сжатии в условиях экспериментов [1] с применением относительно простой и экономичной модели работы [4]. Исследование проводится путем сравнения результатов, рассчитанных с использованием упрощенной модели [4], с результатами, которые позволяет получить полная модель [2].

В модели [2] учитываются влияние неоднородности распределения параметров в паре и жидкости на изменение формы и объема пузырька, эффекты нестационарной теплопроводности

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

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

Этим и вызвана необходимость в проведении данного исследования.

2. Постановка задачи

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

Задача сверхсильного сжатия кавитационного пузырька рассматривается в следующей постановке. Давление в жидкости pX в процессе сжатия считается постоянным и равным 15 бар (в экспериментах [1] оно на фазе сжатия изменяется в пределах от 14 до 16 бар). Полагается, что в начальный момент времени (t = 0) радиус пузырька R равен Rm = 500 мкм, радиальная скорость R (здесь и далее точка сверху означает производную по времени) равна нулю, температура T пара и жидкости равна T0 = 273 K, давление p пара равно давлению насыщения pS0, плотность р жидкости равна рю = 858 кг/м3.

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

rs (0, t) = R(t) + X ai (t)Pi (cos 0).

i=2

Здесь ai - амплитуда (со знаком) отклонения формы пузырька от сферической в виде сферической гармоники, определяемой полиномом Лежандра Pi (cos 0) степени i. Величина искажения сферичности в виде Pi (cos 0) характеризуется амплитудой (со знаком) вг- = a JR. В рассматриваемом случае малой несферичности пузырька амплитуда искажения |вг| мала при всех i (|вг- | << 1). При малых |ег-| взаимодействием искажений формы с разными номерами i можно пренебречь. В силу этого эволюцию искажений сферичности пузырька можно изучать для разных номеров i отдельно, что и делается в настоящей работе. С учетом этого уравнение поверхности пузырька в начале сжатия имеет вид

rs(0) = Rm + ai,mPi (cos 0) (ai,m - начальная амплитуда), который практически остается таким же и в процессе всего сжатия:

rs(0,t) = R(t) + ai(t)Pi (cos 0). Принимается, что в начале сжатия вг- = 0.0003, в г- = 0.

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

3. Полная модель

В этой модели движение как пара в пузырьке, так и окружающей жидкости описывается следующими уравнениями

dt

-pV-u = 0,

p ~Г+Vp = 0, dt

dE

p— + V-(pu - kVT ) = 0.

Здесь и - вектор скорости частицы среды, E = e + и2/2 - удельная полная энергия, e -удельная внутренняя энергия, к - коэффициент теплопроводности.

Применяются реалистичные уравнения состояния жидкого и парообразного дейтерирован-ного ацетона в форме Ми-Грюнайзена [5].

Граничные условия на поверхности пузырька г = Ы0, 0 имеют вид:

р+^Б - и+)-п = р-(б - и-^-п = у,

P + = Р~,

pVT • n )+-pVT • n) = jl,

где D = дг8 /дt =Dn0 - скорость элемента поверхности, п°- внешняя единичная нормаль, у - интенсивность фазовых преобразований, l - теплота парообразования. Знак плюс относится к стороне жидкости, а минус - к стороне пара. Интенсивность фазовых преобразований у определяется по формуле Герца - Кнудсена - Ленгмюра [3,5]:

у = | (Р5(Т-)-х

'4-

2nRvT

где

= exp(-Q2) - qVtc

( 2 ° ^

1----;= I exp Р-Х2 ) dx

Q =

л/Л

ч

jyRT V2 p

0

Здесь аас - коэффициент аккомодации, -

газовая постоянная для пара. Посредством пара-

метра х учитывается подвижность поверхности пузырька. Зависимости k+(T), k(T), pS(T), l(T) представляют собой аппроксимации экспериментальных данных [3]. Коэффициент аккомодации в настоящей работе выбирается равным aac = 1.

На внешней поверхности жидкости r = rf (0, t) задаются давление и температура:

p = Px, T = T0 •

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

Метод численного решения системы уравнений (1) с граничными условиями (2)-(4) строится на основе принципа расщепления по физическим процессам [6]. На первом этапе игнорируется теплопроводность газа и жидкости. В уравнениях газовой динамики тепловые потоки не учитываются. Решение находится численно с помощью модификации схемы Годунова второго порядка точности по пространству и времени на основе UNO-схемы (UNO - uniformly nonoscillatory) [7]. Влияние теплопроводности учитывается на втором этапе. Для этого поле полной энергии корректируется с учетом тепловых потоков, определяемых из уравнения теплопроводности по неявной схеме переменных направлений. Эффективность применения данного метода численного решения системы (1) для изучения искажений сферичности газового пузырька при его сильном сжатии показана в [8].

4. Упрощенная модель

С учетом малости искажений сферичности пузырька движение жидкости и пара в упрощенной модели представляется как суперпозиция сферического (радиального) движения и его малого несферического возмущения. Радиальное движение описывается соответствующими системе (1)-(4) одномерными уравнениями.

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

(1 , )а , зК. (/-1)(1 -д{)К 0

(1+^) а +—а1------------ ------а = 0,

К К

(i + 1)р g

i pl0

Здесь ря = ря(0 - плотность пара в пузырьке. В рамках этой модели жидкость предполагается несжимаемой. В ходе всего сжатия при описании несферического возмущения плотность жидкости полагается равной р10.

Аналогичная упрощенная модель применяется в [4]. В отличие от модели настоящей работы в ней учитывается влияние вязкости жидкости, поверхностного натяжения, вихревого движения жидкости, но пар полагается гомобарическим и пренебрегается влиянием его плотности.

Обыкновенные дифференциальные уравнения (5) решаются методом Рунге-Кутта высокого порядка точности.

5. Радиальное движение На рис.1 показана зависимость К(0 при чисто сферическом сжатии пузырька (а) и изменение радиального распределения давления в его финальной стадии (Ь). Радиус пузырька при сжатии меняется от Кт = 500 мкм до Кса1 = 18.8 мкм в момент коллапса (под моментом коллапса понимается момент экстремального сжатия пара в пузырьке). Распределения давления, температуры и плотности в финальной стадии сжатия оказываются сильно неоднородными, и в конце сжатия в пузырьке образуется радиально сходящаяся ударная волна. В силу малости искажений сферичности межфазной поверхности изменение радиуса пузырька К(0 в двумерном и одномерном случаях практически не различается (рис.1а).

Рис.1. Изменение радиуса пузырька К (кружком отмечен момент коллапса) (а) и радиального распределения давления (в четыре последовательных момента времени (кривые 1-4) финальной стадии) при сжатии чисто сферического пузырька (точкой указана граница пузырька, слева от нее - пар, справа - жидкость) (Ь).

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

Рис.2. Изменение средней по объему плотности пара в пузырьке р (кривая 1) и ее граничного значения р-(кривая 2) в ходе чисто сферического сжатия пузырька.

На рис.2 представлено изменение средней по объему плотности пара в пузырьке р, и его граничного значения р-. Следует отметить, что различие парового и пустого пузырьков проявляется в первую очередь в том, что масса пара в значительной мере определяет глубину коллапса. Схлопывание пустого сферического пузырька было бы полным (до К = 0) и происходило бы со скоростью К ^ да. При наличии же сколь угодно малой начальной несферичности пустой пузырек при таком сжатии разрушается.

Представленные на рис.2 кривые определяют различия в эволюции искажений при использовании двух вариантов упрощенной модели (с ря = р и ря = р-). На низкоскоростной стадии сжатия при К/Кт > 0.2 плотность пара остается малой. В начале высокоскоростной стадии сжатия при К/Кт < 0.2 плотность пара заметно возрастает. При этом величина р- увеличивается значительно быстрее, чем р , поскольку на этом временном интервале плотность пара растет быстрее в приграничной области. Вблизи фокусировки ударной волны соотношение между величинами р и р- меняется. Среднее значение плотности по объему пузырька р начинает превосходить ее граничное значение р-, так как плотность в центральной зоне пузырька оказывается выше, чем у ее поверхности, в силу схождения сферической волны сжатия.

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

6. Расчет величины малых искажений

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

При использовании полной модели численная сходимость результатов расчетов эволюции малых искажений сферичности пузырька в процессе всего сжатия достигается за приемлемое время счета на компьютере лишь при 2 < / < 6.

Поэтому при 7 > 7 величина малых искажений сферичности пузырька в конце сжатия в настоящей работе определяется приближенно.

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

Рис.3. Эволюция амплитуды искажения в финале сжатия при 7 = 7, рассчитанная по гидродинамической модели на сетках: 1000x50 (1), 1000*100 (2), 1000x200 (3).

На рис.3 показано изменение величины е7 при 7 = 7 в финальной стадии сжатия. Представлены результаты вычислений на сетках с разбиением по г на 1000 ячеек (300 - в паре, 700 - в жидкости), по 0 - на 50, 100, 200 ячеек. Видно, что по мере измельчения сетки численное решение довольно быстро сходится (что типично для разностных схем второго порядка точности, к классу которых относится схема метода [2]). Разница величины максимального значения амплитуды искажения на соседних сетках с увеличением разбиения по углу 0 в 2 раза уменьшается примерно в 3 раза. Приблизительно такая скорость сходимости наблюдается и при 2 < 7 < 6. С учетом этого предполагается, что такая скорость сходимости имеет место как при дальнейшем измельчении сетки в случае 7 = 7, так и при ее измельчении при 7 > 8. Данное допущение применялось при определении приведенных ниже

приближенных значений для 7 < 7 < 10.

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

На рис.4 показаны зависимости величины

col / ~

S; / si>m, определяющей рост амплитуды иска-

жения пузырька при сжатии, от номера i в интервале 2 < i < 10, полученные как по обоим вариантам упрощенной модели

значение si в процессе сжатия. При 2 < i < 10 оно достигается либо в момент коллапса, либо незадолго до него. Как видно, зависимости sco1 / вг-, m от номера i имеют колебательный характер, что обусловлено тем, что величина Si в процессе сжатия изменяется, как известно, в режиме колебаний.

200 I----------1---------1----------

2 4 6 8 * 10

Рис.4. Зависимости величины / е,т от 7, полученные по упрощенной модели при рг = р (1) и рг = р- (2) и полной модели (символы).

По рис.4 можно заключить, что на отрезке

2 < 7 < 10 упрощенная модель дает весьма удовлетворительное приближение степени роста амплитуды искажения пузырька при сжатии. При этом в варианте рг = р- результаты получаются несколько лучше. При 7 > 10 применение полной гидродинамической модели не представляется возможным в силу чрезвычайно больших потребностей компьютерного времени. При 7 > 10 оказывается возможным применять только упрощенную модель, так как она позволяет практически одинаково быстро рассчитывать как низкочастотные, так и высокочастотные возмущения сферичности пузырька.

Методика применения упрощенной модели

[4] для оценки максимально возможного роста амплитуды малых искажений сферичности пузырька при сжатии основана на том, что зависимости амплитуды искажения в момент экс-

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

^ і col / і

максимальных значений зависимости \ Sy / є^І

от i. Такие огибающие, построенные по результатам применения упрощенной модели с pg = р и

pg = p-, приведены на рис.5. Видно, что в области низкочастотных возмущений (i = 2-5) огибающие для обоих вариантов упрощенной модели практически совпадают. При увеличении i в области i > 5 они начинают различаться, причем в варианте с pg = p- оценки максимально возможного роста получаются несколько более высокими. В интервале 2 < i < 10 результаты применения этого варианта упрощенной модели лучше согласуются с результатами применения полной модели (рис.4, 5). С учетом этого можно заключить, что для оценки максимально возможного роста амплитуды искажений сферичности пузырька при сжатии вариант упрощенной модели с = р является более предпочтительным.

Рис.5. Огибающие зависимостей | sco1 / е,т| от номера i, полученные по упрощенной модели при Pg = р-(кривая 1) и pg = р (кривая 2), и фрагмент такой зависимости, рассчитанной при pg = р- (кривая 3); симво-

I col / I ___

лами показаны значения | Si / е^т|, полученные по

полной модели.

Заключение

Исследована возможность получения оценок максимально возможного роста амплитуды малых искажений сферичности кавитационного пузырька при его сверхсжатии в условиях экспериментов по акустической кавитации дейтериро-ванного ацетона [1] с применением относительно простой и очень экономичной модели [4] (с учетом влияния плотности пара). Данная упрощенная модель основана на расщеплении двумерного движения жидкости и пара на сферическую составляющую и ее малое несферическое возмущение. Ряд допущений в ней является в конце

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

о правильности получаемых при таких допущениях результатов делается путем их сравнения с результатами применения полной двумерной гидродинамической модели [2], являющейся обобщением одномерной модели Р.И.Нигмату-лина [3].

Установлено, что для получения оценок максимально возможного роста амплитуды малых искажений сферичности кавитационного пузырька при его сверхсжатии в условиях экспериментов [1] упрощенная модель [4] является вполне пригодной. При этом вариант упрощенной модели с выбором в качестве средней плотности пара ее граничного значения более предпочтителен. При таком выборе оценки максимально возможного роста амплитуды малых низкочастотных (2 < i < 10) искажений получаются близкими к тому, что дает полная гидродинамическая модель.

Провести обоснование применимости упрощенной модели сравнением ее результатов с результатами полной гидродинамической модели при i > 10 не представляется возможным, так как потребности компьютерного времени при применении полной модели становятся неприемлемо большими. По этой же причине максимальную амплитуду искажений сферичности пузырька при 7 < i < 10 при использовании полной модели удается рассчитать только приближенно.

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

Работа выполнена в рамках программы РАН и при поддержке РФФИ (номер проекта 08-01-00215).

1. Taleyarkhan R.P., West C.D., Cho J.S., Lahey R.T., Jr., Nigmatulin R.I., BlockR.C. Evidence for nuclear emissions during acoustic cavitation // Science. -2002. - V.295. - P.1868-1873.

2. Аганин А.А., Халитова Т.Ф., Хисматуллина Н.А. Метод численного решения задач сильного сжатия несферического кавитационного пузырька // Вычислительные технологии, 2010, Т.15, №1, С.14-32.

3. Nigmatulin R.I, Akhatov I.Sh., Topolnikov A.S., Bolotnova R.Kh., Vakhitova N.K., Lahey R.T. (Jr.), Taleyarkhan R. The theory of supercompression of

vapor bubbles and nano-scale thermonuclear fusion // Phys. Fluids, 2005, V.17, №10, 107106.

4. Аганин А.А., Ильгамов М.А., Нигматулин Р.И., Топорков Д.Ю. Эволюция искажений сферичности кавитационного пузырька при акустическом сверхсжатии // МЖГ, 2010, №1, С.57-69.

5. Нигматулин Р.И. Динамика многофазных сред, т.

1 и 2. М.: Наука, 1987.

6. Kondic L., Gersten J.I., Yuan C. Theoretical studies of sonoluminescence radiation: Radiative transfer and parametric dependence // Phys. Rev., 1995, V.52, P.4976-4990.

7. Harten A., Engquist B., Osher S., Chakravarthy S.R. Uniformly high order accurate essentially non-oscillatory schemes III // J. Comp. Phys., 1987, V.71, P.231-303.

8. Аганин А.А., Ильгамов М.А., Халитова Т.Ф. Моделирование сильного сжатия газовой полости в жидкости // Математическое моделирование, 2008, T.20, №11, C. 89-103.

9. Lin H., Storey B.D., Szeri A.J. Inertially driven inhomogeneities in violently collapsing bubbles: the validity of the Rayleigh-Plesset equation //J. Fluid Mech., 2002, V.452, P.145-162.

SIMULATION OF SMALL SPHERICITY DISTORTION OF A CAVITATION BUBBLE DURING ITS SUPERCOMPRESSION

A.A.Aganin, D.Yu.Toporkov, T.F.Khalitiva, N.A.Khismatullina

Small (axially symmetric) distortions of the spherical shape of a cavitation (vapor) bubble during its superstrong compression under conditions of experiments on acoustic cavitation of deuterated acetone (Taleyarkhan et al, Science. 2002. V.295) are considered in this article. Applicability of relatively simple and very efficient-in-application model (based on splitting the two-dimensional flow of the liquid and the vapor into the spherical part and its small non-spherical perturbation) to estimate maximum growth of the amplitude of such distortions is studied. While describing evolution of the distortion, a number of assumptions are used in this model, the validity of which has not been proved rigorously in the high-speed final compression stage. It is shown that even under such assumptions the evaluated estimates are quite satisfactory, which in the case of low-frequency distortions is supported by the comparison with the results of application of two-dimensional hydrodynamic model.

Key words: bubble dynamics, acoustic cavitation, cavitation bubble, bubble collapse, stability of the spherical shape.

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

E-mail: aganin@kfti.knc.ru

Топорков Дмитрий Юрьевич - кандидат физико-математических наук, старший научный сотрудник Института механики и машиностроения Казанского научного центра РАН.

E-mail: dem17@rambler.ru

Халитова Талия Фаритовна - кандидат физико-математических наук, научный сотрудник Института механики и машиностроения Казанского научного центра РАН.

E-mail: taliny@mail.ru

Хисматуллина Наиля Абдулхаевна - кандидат физико-математических наук, доцент кафедры вычислительной математики и информатики Татарского государственного гуманитарно-педагогического университета.

E-mail: nailya_hism@mail.ru