ОПТИМИЗАЦИЯ ПРОГРАММНОГО БЛОКА В ЗАДАЧАХ МЕХАНИКИ И ЭЛЕКТРОДИНАМИКИ ПРИСТЕНОЧНОЙ ПЛАЗМЫ Текст научной статьи по специальности «Физика»

Электронный журнал «Труды МАИ». Выпуск № 50

www.mai.ru/science/trudy/

УДК 533

Оптимизация программного блока в задачах механики и электродинамики пристеночной плазмы М.В. Котельников, Нгуен Суан Тхау

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

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

Особое внимание уделяется возможности применения классического условия устойчивости Куранта-Фридрихса-Леви к задачам данного класса.

Ключевые слова: пристеночная плазма; уравнения неразрывности, уравнение Пуассона, метод крупных частиц, условие устойчивости Куранта-Фридрихса-Леви, расчетная сетка, шаг по времени, метод Рунге-Кутта.

Введение

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

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

Математическая модель задачи

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

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

В качестве примера рассмотрим бесконечно длинный цилиндр радиуса ^ и потенциала фр,, помещенный в поперечный поток слабоионизованной плотной плазмы. Скорость потока Ua,. Магнитное поле предполагается отсутствующим.

Математическая модель задачи включает:

- систему уравнений Навье-Стокса (или Эйлера) для нейтральной компоненты плазмы;

- уравнения неразрывности для заряженных компонент;

- уравнения движения для заряженных компонент;

- уравнения Пуассона для самосогласованного электрического поля.

Уравнение энергии для ионов не использовалось, поскольку в условиях слабой степени ионизации температуры ионов и нейтральных атомов совпадают (^ = Уравнение энергии для электронов заменяется модельным соотношением 8=^/^=^^^ при этом величина в варьируется.

В представленной ниже математической модели учтены также следующие допущения:

- химические реакции заморожены;

- вязкостью пренебрегается;

- собственное магнитное поле мало. сП-

-н- + div(n; v;) = 0 ot

Эп

—e + div(ne ve) = 0 ot

mi dVL = - Vni + Ze(E + Vi x B) - ^ Via (Vi - Va ) dt n;

dv IcT

me dVe =--e Vne - e(E + Ve x B) - ^a Vea (Ve - Va )

dt ne

Аф = e(ne-Zn;)/£o, E = -Уф 1)

^ + div(Pa Va ) = 0

at

0(pa Va )

a

0(PaEa )

a

+ div (Pa V a Va ) = -VPa

+ div (pa VaEa ) = -div ( Vapa )

Ра =(Еа " У^/2)ра (У"1) , Еа = СуТа + У^ / 2

В приведенной системе уравнений п, v, T, р, P, д, V, ф, Е, В - концентрация, скорость, температура, плотность, давление, приведенная масса, частота столкновений, потенциал и напряженность электрического поля, индукция магнитного поля. Остальные обозначения общепринятые. Индексы ^ e, a относятся к ионам, электронам и нейтральным частицам.

Математическая модель задачи включает кроме системы уравнений (1) стандартную систему начальных и граничных условий [1-4].

Численная модель задачи

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

временном слое определяется самосогласованное электрическое поле путем решения уравнения Пуассона.

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

В работах Ю.М. Давыдова и его учеников [5,6] обоснована необходимость использования однородных и изотропных вычислительных пространств. В случае тел сложной геометрии предлагается использование метода дробных ячеек. Чтобы не нарушать единообразие вычислений и не применять специальных формул для ячеек на границе расчетной области, предложен метод фиктивных ячеек.

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

Оптимизация вычислительного алгоритма включает:

- оптимизацию размера расчетной области;

- оптимизацию числа крупных частиц (или величин шагов по координатам);

- оптимизацию величины шага по времени;

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

- оптимизацию по числу используемых для расчета членов ряда (например, гармоник ряда Фурье), если используются численно-аналитические методы, связанные с предварительными аналитическими преобразованиями - разложениями в ряды.

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

• г0 = - безразмерный радиус цилиндра;

• ф0 = ф^М^) - безразмерный потенциал тела; (2)

• 8=^/^ - отношение температур ионов и электронов.

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

Рис. 1 Зависимость размера возмущенной зоны вблизи сферического тела от г/Мг и е=Т,/Тс

Здесь Мг = гс = (в0кТ;/(п;оэе2}- масштаб длины.

Как следует из рисунка, размер расчетной области в достаточно широком диапазоне изменения характерных параметров задачи составляет (100го^200го). В интервале изменения параметров задачи, необходимых для практики, А проявила достаточно слабую зависимость от потенциала тела. Более сильная зависимость проявляется от характерного размера тела г0 и отношения температур в.

При наличии направленной скорости в теневой области тела возникает удлиненная возмущенная зона в виде «следа». Следовательно в теневой области размер возмущенной зоны существенно возрастает. Для сокращения необходимых дополнительных ресурсов ЭВМ необходимо учесть следующее:

1. Если целью расчета является детальное исследование всей возмущенной зоны, то «след» должен целиком уместиться в расчетную область. Поскольку размер «следа» со временем нарастает постепенно, увеличивать размер расчетной области можно также постепенно, пропорционально увеличению «следа».

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

3. Если целью расчета является вычисление тока заряженных частиц на тело, необходимо контролировать зависимость тока от размера «следа». Методические расчеты показали, что ток в

<1, (3)

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

Рассмотрим вопрос о выборе шага по времени Дt. Как показано в [8], условие устойчивости Неймана применительно к гиперболическим уравнениям будет выполняться, если использовать известную схему Лакса [7]. Анализ устойчивости схемы, проведенный в [8], приводит к неравенству уД1

дХ

где v - физическая скорость среды, Дt - шаг по времени, Дx - шаг по пространственной координате. Отсюда вытекает ограничение на шаг по времени

* ДХ

Д* < |-1, (4)

У

| тах |

где vmax - максимальная скорость крупных частиц. Это условие для шага по времени является примером условия устойчивости Куранта-Фридрихса-Леви (КФЛ) [9] применительно к гиперболическим уравнениям. Полученное неравенство означает, что в явном методе шаг по времени нужно выбирать меньше наименьшего характерного физического времени задачи, которое в случае уравнений переноса есть время, за которое скорость v приводит к переносу на расстояние Дx. Условие Куранта-Фридрихса-Леви выражает требование, чтобы физическая

Дх

скорость v была меньше сеточной скорости .

Однако многочисленные вычислительные эксперименты [1-4] с пристеночной плазмой показали, что приведенное ограничение на шаг по времени оказывается слишком жестким. Если потенциал стенки достаточно высок, то кривая распределения потенциала имеет максимальный градиент вблизи поверхности, следовательно, вблизи поверхности максимален перенос, связанный с подвижностью. Во многих задачах вблизи поверхности максимален также и градиент концентрации заряженных частиц, поэтому максимален перенос, связанный с процессом диффузии. Из сказанного следует, что именно вблизи поверхности находятся самые высокоскоростные частицы среды, которые и влияют на шаг по времени через условие Куранта-Фридрихса-Леви. Для подтверждения сказанного приведем несколько примеров из конкретных вычислительных экспериментов. На рис. 2,3 рассматривается расчетная сетка и поле скоростей ионов для заряженного цилиндра, помещенного в поперечный поток слабоионизованной столкновительной плазмы. Расчеты проведены при числе Маха М=0.6, электрическом числе Рейнольдса Reэ=25, безразмерном потенциале цилиндра ф0=ф/Мф=-45 и значении безразмерного

радиуса ^=/^=10 (Mф=kTi/e - масштаб потенциала, Мг = гв = (в0кТ;/(п;оэе2}- масштаб длины).

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

Рис. 2. Вычислительная сетка

Рис. 3. Поле скоростей ионов (ro=10;Reэ=25;M=0,6;фo=-45)

Таблица 1

1 2 3 4 5

1 -2,78 -2,20 -1,76 -1,35 -1,05

20 -2,32 -1,88 -1,36 -1,04 -0,71

40 -1,66 -0,99 -0,47 -0,18 0,08

60 -2,17 -1,87 -1,66 -1,45 -1,26

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

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

Приведенный пример и другие вычислительные эксперименты позволяют сделать вывод, что шаг по времени может быть увеличен в 3-5 раз по сравнению с его величиной, вытекающее из условия КФЛ (формула (4)) без нарушения устойчивости и точности решения. Это условие было выведено для линейных гиперболических уравнений и имеет наглядную физическую интерпретацию: крупная частица за один шаг по времени не может продвинуться больше, чем на размер одной Эйлеровой ячейки. Поскольку в задачах электродинамики пристеночной плазмы система уравнений существенно нелинейна, условие устойчивости будет иным и оптимальный шаг по времени определяется из методических расчетов, а условие КФЛ используется как начальное, с которого начинается исследование.

Оптимальные шаги Эйлеровой сетки по радиальной и угловой координатам, а следовательно и оптимальное число узлов Эйлеровой сетки определялись также из методических расчетов: это число должно быть минимальным при сохранении устойчивости и заданной точности решения. Для рассматриваемой в данной работе задачи обтекания заряженного цилиндра потоком слабоионизованной плотной плазмы оптимальное число узлов по радиусу N=120, по углу N9=180.

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

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

Заключение

Задачи, подобные рассмотренному примеру, успешно решаются на ЭВМ среднего класса с частотой одноядерного процессора порядка 1800 Гц. Время расчета одного варианта в зависимости от выбора характерных параметров задачи при указанной оптимизации вычислительного алгоритма колебалось от нескольких часов до нескольких суток.

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

1. Алексеев Б.В., Котельников В.А. Зондовый метод диагностики плазмы. М.: Энергоатомиздат, 1988, 239 с.;

2. Котельников В.А., Гурина Т.А., Демков В.П., Попов Г.А. Математическое моделирование электродинамики летательного аппарата в пристеночной плазме. М.: Изд. Нац. акад. прикл. наук, 1999, 265 с.;

3. Котельников В.А., Ульданов С.В., Котельников М.В. Процессы переноса в пристеночных слоях плазмы. М.: Наука, 2004, 422 с.;

4. Котельников В.А., Гидаспов В.Ю., Котельников М.В. Математическое моделирование обтекания тел потоками бесстолкновительной и столкновительной плазмы. Изд-во Физматлит, 2010, 288 с. Поддержано РФФИ, грант № 08-08-13586 ОФИ-Ц;

5. Давыдов Ю.М., Скотников В.П. Метод крупных частиц: вопросы аппроксимации, схемной вязкости и устойчивости. М.: ВЦ АН СССР, 1978. 71 с.;

6. Давыдов Ю.М., Скотников В.П. Исследование дробных ячеек в методе крупных частиц. М.: ВЦ АН СССР, 1978. 71 с.;

7. Марчук Г.И. Методы вычислительной математики. М.: Наука, 1989, 608 с.;

8. Поттер Д. Вычислительные методы в физике. М.: Мир, 1975, 329 с.;

9. Давыдов Ю.М. Дифференциальные приближения и представления разностных схем. М.: МФТИ, 1981. 131 с.;

10. Котельников М.В. Механика и электродинамика пристеночной плазмы. Дисс. на соиск. ученой степени д.ф.-м.н. Москва, МАИ, 2008г, 271 с.

Сведения об авторах

Котельников Михаил Вадимович, профессор Московского авиационного института (национального исследовательского университета), д.ф.-м.н., доцент. МАИ, Волоколамское ш., 4, Москва, А-80, ГСП-3, 125993; тел.: (499) 158-19-70; e-mail: mvk home@mail.ru

Нгуен Суан Тхау, аспирант Московского авиационного института (национального исследовательского университета), Волоколамское ш., 4, Москва, А-80, ГСП-3, 125993; тел: 8-965-377-20-30; e-mail: nguoikinhbac.ruc.@gmail.com

12.01.2012