Алгоритмическое вовлечение дополнительной информации для численного решения уравнений Навье - Стокса в проекте Tristam Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 18, Специальный выпуск, 2013

Алгоритмическое вовлечение дополнительной информации для численного решения уравнений Навье — Стокса в проекте TRISTAM*

В. В. ШАйдуров1,2, А. В. Вяткин1 1 Институт вычислительного моделирования СО РАН, Красноярск, Россия 2 Университет Бейхан, Пекин, КНР e-mail: shaidurov04@mail.ru

Вычислительная математика накопила целую палитру высокоточных численных методов, которые реализуют решение современных математических моделей аэродинамики на высокоэффективных вычислительных системах. В статье схематически описана четырёхстадийная технология для решения уравнений Навье — Стокса вязкого газа, реализуемая в рамках международного проекта ТШБТАМ.

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

Введение

Вычислительная аэродинамика, как правило, дополняет физические эксперименты, а в некоторых случаях опережает их. Успешные физические эксперименты приводят к уточнениям математических моделей, например, системы уравнений Навье — Стокса для вязкого теплопроводного газа и их усреднённой по Рейнольдсу формой. С другой стороны, бурное развитие компьютерных технологий привело к использованию компьютеров на всех стадиях проектирования летательных аппаратов. Однако, несмотря на рост вычислительной производительности, моделирование аэродинамического вязкого течения с большими числами Рейнольдса около сложной поверхности до сих пор обходится довольно дорого и часто не даёт приемлемой точности. Поэтому решение уравнений Навье — Стокса и их модификаций для задач аэродинамики с приемлемой точностью остаётся довольно серьёзной задачей [1-8].

В международном проекте ТШБТАМ мы используем современные вычислительные алгоритмы в рамках четырёхстадийной вычислительной технологии с адаптацией на выбранную цель моделирования.

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

* Работа выполнена при частичной поддержке РФФИ (грант № 14-01-31203).

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

1. Четырёхстадийная вычислительная технология

Разработанная технология решения задачи состоит из четырёх стадий.

Первая стадия называется "Формирование данных и предварительное решение". Её задачами являются формализация геометрии вычислительной области, хранение данных, решение уравнений Навье — Стокса методом конечных элементов на почти равномерной сетке после дискретизации по времени. В нашем подходе мы предпочитаем использовать прямоугольную сетку, равномерную почти везде, за исключением тонкой полосы вдоль поверхности обтекаемого тела. Вычисления на этой стадии должны быть проведены с точностью, достаточной для определения проблемных зон вычислительной области, в которых наблюдаются ударные волны, пограничные слои и другие особенности. На этой стадии необходимо лишь определить проблемные зоны, не пытаясь достигнуть максимальной точности численного решения. Точность приближённого решения будет неоднократно повышаться в последующих стадиях технологии!

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

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

Четвёртая стадия называется "Целевое задание" и состоит из трёх шагов. Целевое задание означает вычисление некоторых важных параметров, таких как лобовое сопро-

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

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

Отметим, что на второй стадии важную роль играют апостериорные оценки точности. В нашем подходе они основаны на прямом вычислении невязки. Для этого используем эрмитовы конечные элементы из классов Н2 и С(1), которые позволяют производить такие вычисления.

Преимущество точного вычисления важных параметров вместо уточнения полного решения во всей области определения уже было продемонстрировано в рамках проекта АБЮМА, когда повышение точности требуемых параметров на несколько порядков достигалось без существенного увеличения вычислительных затрат [1, 10].

2. Некоторые постановки задачи

В качестве иллюстрации рассмотрим типичный пример задачи двумерного моделирования течения газа около тела на сверхзвуковой скорости. Пусть ПС0тр = (0, Н1) х (0, Н2) — прямоугольная вычислительная область (рис. 1) с точками х = (х,у). Её граница Гсотр содержит левый участок втекания газа Г;п, а также правый участок и две другие стороны вытекания газа, обозначенные символом Го^. Внутри вычислительной области размещено твёрдое тело П80па с границей Г80щ. Таким образом, геометрической областью нашей задачи является П = ПС0тр\(П80ш и Г30ца) с границей Г = Гсотр и Г80ш.

Для упрощения выкладок введём выражение означающее производную,

которая подобна субстанциональной (лагранжевой) производной:

Ба _ да д(аи) д(аь) БЬ дЬ дх ду

С учётом этого выражения уравнения Навье — Стокса записываются в виде системы четырёх уравнений с неизвестными функциями р, и, ь и е в цилиндре (Ь, х) Е (0, Т) х П :

БР = о; (1)

Рис. 1. Обтекаемое тело в вычислительной области

О(рп) + I дх

дг

ТХХ) Г) Тху 0;

т(рь) д д

ду'

т

дхТху + ду (Р Туу)

т(Ре) + дк + дк = _ р(— + — \ + ф дг дх дх \дх ду /

(2)

(3)

(4)

Здесь р(г,х,у) — плотность; п(г,х,у), ь(г,х,у) — компоненты вектора скорости и = (п,ь); е(г,х,у) —внутренняя энергия; Р(г,х,у) —давление.

Поскольку в четырёх уравнениях фигурирует 5 неизвестных функций, для исключения произвола используется дополнительное алгебраическое соотношение — уравнение состояния идеального газа

Р = (7 - 1)ре, (5)

где константа ^ ^^

1.4 является характеристикой воздуха. Соотношение (5) не является единственно возможным; другие соотношения приведены, например, в [3] для различных сред и при разных ограничениях.

Три функции тхх, тху, туу являются компонентами тензора давления

3

Тх

Тх

ху

Тх

ху

Ту

уу

2 /дп дь\ Тхх = эЧ2дх - ду)

Ту

уу

где

2 /^дь дп

3 V \ ду дх

1

Тх

ху

' дп дь ду дх

^(г,х,у,х) = — ^*(t,x,y,z), Не

а V* — динамический коэффициент вязкости:

V

- 1)тМо)ш еш, 0.76 < ш < 0.9.

Здесь И,е — число Рейнольдса (как правило, достаточно большое: И,е ^ 1); ~ 1 — число Маха; ш — фиксированная константа из заданного интервала. Другие виды соотношения (8) тоже приведены в [3] для различных сред и условий.

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

Чх

7 де Рг ^ дх,

Чу

7 де Рг ^ ду

(9)

Здесь Рг ^ 0.72 — число Прандтля для воздуха.

Функцию диссипации Ф(г,х,у) представим следующим образом:

Ф

V

2/ дп\2 2/ дь 3 \дх) 3 \ду

дь дп\ / дп дь

дх ду дх - ду

Для начала расчётов по времени используем следующие начальные данные:

р(0, х) = ро(х) > 0 при х € П;

и(0, х) = и0(х) при х € П; е(0, х) = е0(х) > 0 при х € П.

:ю)

;ц) '12) :13)

2

2

При иллюстрации применяемых методов будем ориентироваться на моделирование сверхзвукового течения. Поэтому пока ограничимся краевыми условиями на участке втекания газа rin с невозмущённым потоком:

p(t, x) = Pin(t, x) > 0 при (t, x) e (0, T) x r¡n; (14)

u(t, x) = uin(t, x) при (t, x) e (0,T) x rin; (15)

e(t, x) = ein(t, x) > 0 при (t, x) e (0,T) x rin. (16)

Для функции p справедлив закон сохранения массы

= -J p (u ■ n) dr (17)

П' Г'

в любой подобласти П' С П с границей Г'. Здесь и далее выражение n(t, x) означает вектор внешней нормали, а слагаемое u ■ n — скалярное произведение двух векторов. Для упрощения полагаем, что нормальная компонента скорости un = u ■ n является отрицательной на границе втекания rin (это означает приток массы) и неотрицательной на всех остальных частях границы r\rin (это подразумевает отсутствие притока газа через эти участки границы).

Один из простых и эффективных методов аппроксимации уравнения (1) состоит в использовании явных монотонных разностных схем [6-8]. Этот способ имеет ограничение на шаг по времени при большом значении числа Рейнольдса. Вторым законом является закон сохранения момента

djpu = -J ((u ■ n)pu + Pn -Sn) dr, (18)

П' Г'

справедливый для любой подобласти П' С П с границей Г'.

Оба закона сохранения (17) и (18), как правило, используются для построения дискретной аппроксимации методом конечных объёмов в своей исходной формулировке. Что касается внутренней энергии, то закон сохранения для неё совмещает e с u:

p (e + (u2 + v2)/2) dn =

П'

= - y p (e + (u2 + v2)/2) (u ■ n) dr - J (Pu - Su + q) ■ n dr (19)

Г' Г'

для любой подобласти П' С П с границей Г'. Первое подынтегральное выражение со слагаемым e соответствует взвешенной норме в пространстве Li (в силу положительности e). Однако два следующих слагаемых соответствуют взвешенной норме в пространстве L2. Таким образом, разные части одного равенства соответствуют разным нормам, что приводит к усложнению исследования свойств решения задачи. Более того, функциональное пространство L1 не является гильбертовым. В связи с этим функцию e целесообразно заменить на другую искомую функцию

П = e1/2.

(20)

С физической точки зрения внутренняя энергия e является неотрицательной функцией, поэтому n(t, x) — вещественная функция. Таким образом, заменяя в (5) e на п2, получим эквивалентное уравнение

D(pn2) , dqx dqy п f du | dv

+ + = --+--+ф. (21)

БЬ дх ду \дх ду)

После замены е на п2 закон сохранения внутренней энергии перепишется в следующем виде:

|/Р (П2 + (и2 + ь2)/2) СП =

П'

= _ ^ Р (п2 + (и2 + V2)/2) (и ■ п) СГ _ ^ (Ри _ Зи + я) ■ п сСГ. (22)

Г' Г'

Введём ещё одну функцию:

а = р1/2.

Вновь с физической точки зрения плотность р является неотрицательной функцией, поэтому о(Ь, х) — вещественная функция. Таким образом, закон сохранения (22) может быть сформулирован как в терминах ¿2-нормы для вектор-функции (ап, аи, аь), так и в терминах -нормы для другой вектор-функции с компонентами (ре, ри2, рь2). Мы будем использовать первую формулировку. Для этого оставим уравнения (1)-(3) без изменения, а уравнение (21) перепишем в следующем эквивалентном виде:

Б(рп) 1 дцх 1 дцу Р (ди ду\ Ф . ,

+ = _ + + , (23)

БЬ 2п дх 2п ду 2п \дх ду) 2п

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

Р/п = (7 _ 1)рП.

Теперь опишем идею конформного метода конечных элементов на временном слое Ьк = кт с одновременной аппроксимацией производной Б(-)/БЬ по времени. Для этого на слое Ьк рассмотрим одну из базисных функций (р(х,у) (рис. 2) и найдём её восполнение 6(Ь,х,у) между временными слоями (рис. 3) как решение уравнения

дв д6 д6 г дЬ + идх + ьду = 0 Ь е [Ьк-1,Ьк].

Она обладает следующим свойством:

вБМ = БСрив! = вБМ + ^ (24)

БЬ БЬ БЬ ' дЬ v у

Поэтому умножение на эту функцию приводит уравнение (2) к виду

¥ + 6» (Р _ тхх) _ вд|утх, = 0. (25)

Рис. 2. Базисная функция на слое tk с носи- Рис. 3. Восполнение функции между слоями телями и 1 и tk

Проведём интегрирование этого уравнения по носителю V функции в и используем формулу Гаусса—Остроградского:

{^(ври) д ч _ д , 7ТЛ

' - + в— (Р - Тхх) - в — ТХу) ¿V =

II О / \ ХХ /

У \ д£ дх ду

= / Р^ ^^ + / рив + / (дх (Р - Тхх) - дУтх^ в ^ = 0.

ш д V

После деления на т и замены последнего интеграла на интервале , tk] квадратурной формулой правого треугольника со значением в точке £ = tk получается приближённое равенство

[ /ри дР дтхх дтх^\ . [ ри

Дт + дх - дтх- - * «у тв •

ш д

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

Аналогичные уравнения выводятся для V и п.

3. Граничные условия на границе вытекания

С математической точки зрения уравнения (2)-(4) или (2), (3), (23) принадлежат к параболическому типу уравнений второго порядка. Чтобы постановка задачи была корректной, необходимо задать граничные условия в каждый момент времени £ Е (0, Т) на всех границах вытекания. Эти условия основываются на различных предположениях и до сих пор вызывают бурную дискуссию [8]. Наряду с другими подходами к реализации граничных условий будем формировать постановки задач с помощью одного из трёх следующих методов.

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

(0 на Гоиъ

1 на расстоянии > е от Гои^ достаточно гладкая в ПСОтр.

С её помощью переопределим функцию вязкости

КЬ,х,у) = в(х,у) V*(Ь,х,у).

Ке

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

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

Б(ри) = 0 БЬ

Б(рь) 0 при (Ь, х) е [0,Т] х Гоиъ

Dt

D(pn)

0

Dt

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

Более приемлемую форму граничных условий для субзвуковых течений реализует подход "do nothing", который возник при численном решении уравнений Навье — Стокса для вязкой несжимаемой жидкости [12]. Его идея основывается на законах сохранения (18) и (19). Чтобы не создавать дополнительного искусственного импульса в методе конечных элементов от границы вычислительной области, необходимо скорректировать два последних слагаемых в интеграле по границе из (18):

Pn - Sn = Poutn при (t, x) e [0, T] x rout. (26)

С физической точки зрения это означает, что давление на границе вытекания становится постоянным, а трение незначительное.

Подобные рассуждения для (19) (сокращение дополнительной искусственной энергии в методе конечных элементов от границы вычислительной области) с учётом (26) приводят к граничному условию

q ■ n = 0 при (t, x) e [0, T] x rout. (27)

Заключение

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

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

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

[1] ADIGMA — A European initiative on the development of adaptive higher-order variational methods for aerospace applications // Notes on Numerical Fluid Mechanics and Multidisciplinary Design. Springer, 2010. Vol. 113.

[2] Vos J.B., Rizzi A., Darracq D., Hirschel E.H. Navier — Stokes solvers in European aircraft design // Progress in Aerospace Sci. 2002. Vol. 38. P. 601-697.

[3] The Handbook of Fluid Dynamics. CRC Press LLC & Springer, 1998.

[4] Bosnyakov S., Kursakov I., Lysenkov A. et al. Computational tools for supporting the testing of civil aircraft configurations in wind tunnels // Progress in Aerospace Sci. 2008. Vol. 44. P. 67-120.

[5] Oberkampf W.L., Trucano T.G. Verification and validation in computational fluid dynamics // Progress in Aerospace Sci. 2002. Vol. 38. P. 209-272.

[6] Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer, 2009.

[7] Titarev V.A., Toro E.F., Comput J. Finite-volume WENO schemes for three-dimensional conservation laws // J. of Comput. Phys. 2004. Vol. 201. P. 238-260.

[8] Lodato G., Domingo P., Vervisch L. Three-dimensional boundary conditions for direct and large-eddy simulation of compressible viscous flows // J. of Comput. Phys. 2008. Vol. 227. P. 5105-5143.

[9] Pironneau O. On the transport-diffusion algorithm and its applications to the Navier — Stokes equations // Numerische Mathematik 1982. Vol. 38. P. 309-332.

[10] Bangerth W., Rannacher R. Adaptive Finite Element Methods for Differential Equations. Berlin, Birkhauser, 2003.

[11] Shaydurov V., Liu T., Zheng Z. Four-stage computational technology with adaptive numerical methods for computational aerodynamics // American Institute of Physics: Conf. Proc. 2012. Vol. 1484. P. 42-48.

[12] Rannacher R. Methods for numerical flow simulation // Hemodynamical Flows: Aspects of Modeling, Analysis and Simulation Oberwolfach: Lecture Notes of Oberwolfach Seminar / G.P. Galdi, R. Rannacher et al. (eds.), Birkhauser, Basel, 2007.

Поступила в 'редакцию 29 ноября 2013 г.