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

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

Том 11, № 6, 2006

АНАЛИЗ ВЫЧИСЛИТЕЛЬНЫХ СХЕМ МЕТОДОВ КОНЕЧНЫХ ЭЛЕМЕНТОВ И КОНЕЧНЫХ РАЗНОСТЕЙ ДЛЯ МОДЕЛИРОВАНИЯ ТЕЧЕНИЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ

А. В. Говыш

Новосибирский государственный технический университет, Россия

e-mail: agobysh@mail.ru

Н. Ю. ШокинА High Performance Computing Center Stuttgart, Germany e-mail: shokina@hlrs.de

The paper is devoted to the numerical modelling of two- and three-dimensional incompressible fluid flows. Finite element methods (the Uzawa algorithm, a mixed finite element method, a stabilized finite element method) and the finite difference method are employed for the numerical modelling of fluid flows. The results of some numerical experiments demonstrating numerical convergence are presented.

Введение

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

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

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

Несмотря на достоинства стандартного МКЭ при решении задач гидродинамики, такие как консервативность, абсолютная устойчивость, возможность конструирования конечноэлементных аппроксимаций на неструктурированных сетках, остаются актуальными следующие проблемы: выбор интерполяционных полиномов для вектора скорости и давления,

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

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

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

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

Математическая постановка задачи об установившемся протекании идеальной жидкости через односвязную область П в плоскости декартовых координат x1 Ox2 заключается в определении вектора скорости u и давления р, удовлетворяющих в П уравнениям неразрывности и движения:

V ■ u = 0, (1)

(u ■ V)u + Vp = f, x = (x1,x2) Є П, (2)

с заданным вектором массового расхода u на входе Г1 , условиями непротекания на твердых стенках и с заданной нормальной составляющей u ■ n вектора потока массы на выходе Г2. Плотность жидкости положена равной І.

Постановка задачи в переменных функция тока ф — функция вихря ш, которые вводятся как

дф дф _ дщ ди2 ,оЧ

и1 = , и2 = — , ш = rot u = —^~2 + ^ТТ, (3)

дx2 дx1 дx2 дx1

имеет вид

Аф = — ш; (4)

u ■ Vu = rot f. (5)

При этом граничные значения для ф однозначно определяются из условий (3) с точностью

й й , дф до аддитивной постоянной ф0, ——

дп

находится из этих же условий, шL можно определить Г і 111

с помощью условий для ф при численном решении задачи [7].

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

2. Конечно-элементные вычислительные схемы для решения задачи в естественных переменных

Конечно-элементное решение уравнений (1), (2) ищется в пространствах Н1(П) := (Н0(П))2 для и и Ь2 (П) для р. Для вывода слабой формулировки задачи умножим два уравнения в (2) на произвольные У\ £ Н0(П), у2 £ Н^П), уравнение (1) умножим на д £ Ь2(П).

Вариационная постановка в форме Галёркина имеет следующий вид. Найти и £ Н1(П) и р £ Ь2(П) такие, что Vv £ Н^(П), д £ Ь2(П),

J([и • V]u)vdП — Jp(div и)^П = J f • vdП, п п п

! (V ■ и)д^П = 0.

п

Выполним триангуляцию расчетной области. Пусть заданы пространства конечных элементов И и Рн для аппроксимации и и р. Конечно-элементное решение определяется из следующих соотношений. Найти иь € Иь и рн € Рн такие, что € И^, € Р^,

J([иЛ ■ V]uh)vhdП - ^рЛ(&у и^П = J I ■ ^П,

п п п

J(V ■ иЛ)д^П = 0.

п

Для того чтобы перейти от дискретной вариационной постановки к системе алгебраических уравнений, перенумеруем все базисные функции из И^: и^, г = 1,... ,пи (сначала для первого компонента вектора скорости, потом — для второго), далее занумеруем все базисные функции из Рн: ^, г = 1,... , пр. Получим систему уравнений вида

С -Я \ (и\ = (I

Ят о ) \р) (,0

где и € И^; р € Рн; С — конвективная матрица; Я — матрица, соответствующая дивергентному оператору; I — вектор правой части.

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

Смешанный метод конечных элементов. Начальные значения вектора скорости и0 и давления р0 заданы, е > 0, п > 0, компоненты вектора скорости и давления на (п+1)-й итерации вычисляются по алгоритму:

с -д \ /ига+1\ ( і

Q1 єМр ) \eMpp

(6)

где матрица Мр — матрица массы; матрицы С, д и вектор і определены выше; є — произвольный параметр; п — шаг итерации.

Алгоритм Удзавы. Начальные значения вектора скорости и0 и давления р0 заданы, є > 0, п > 0, компоненты вектора скорости и давления на (п + 1)-й итерации вычисляются по следующему алгоритму

с +1 ямр-1ят) ига+1 = I + Ярп, (7)

рП+1 = рп - 1 м-1 ятип+1. (8)

е 1

Здесь матрицы Мр, С, Я и вектор I определены выше, е — произвольный параметр, п — шаг итерации.

Для существования и единственности решения и^ и рн необходимо выполнение условия согласованности на пространствах И и Рн (условие Ладыженской — Бабушки — Брецци (ЛББ)) [5, 6].

Стабилизированные методы конечных элементов (БИРО/РБРО) позволяют использовать пространства, для которых условие ЛББ не выполняется [4]:

С —(Я + К8ИРС) А (иА = ((М + СтиРС)/А (9)

(Я + КР8РС)Т ^РвРО )\р) V ЯР8РС/ /

где матрицы КдиРо, Крдро , ^рдро , Сзирс, Ярзрс — матрицы стабилизации; М — матрица массы; матрицы С, Я определены выше.

Для аппроксимации конвективных членов рассмотрены противопотоковые схемы с расчетом неизвестных в узлах симплекса и серединах его сторон [8].

3. Итерационный конечно-разностный метод для решения задачи в переменных функция тока — функция вихря

С помощью взаимно-однозначного невырожденного отображения единичного квадрата Я, лежащего в плоскости координат с1 Ос2, на область П уравнения (4), (5) записываются в новых независимых переменных са [7]. Полученные уравнения аппроксимируются на криволинейной сетке, которая строится методом эквираспределения [9].

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

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

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

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

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

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

4. Результаты расчетов

4.1. Сравнение конечно-элементных вычислительных схем

Сравнительный анализ конечно-элементных вычислительных схем учета конвективных членов в уравнениях (1), (2) проводился на задаче конвективно-диффузионного типа с гладким решением на последовательности вложенных пространственных сеток

—div (vTgrad T) + u ■ grad T = f, П = [0,1]2, (10)

с известным полем скоростей

u = (x2y + y3, —y2x — x3)T. (11)

Точное решение имеет вид T = x3 + y3. Краевые условия и правая часть задачи задаются в соответствии с этими функциями. В расчетной области задано регулярное симплици-альное разбиение. Исследуются интерполяционные свойства следующих вычислительных схем: метода конечных элементов, в котором используется противопотоковая схема с расчетом неизвестных в вершинах симплексов (схема 1) [8], стандартного метода конечных элементов (МКЭ), метода конечных элементов, в котором используется противопотоковая схема с расчетом неизвестных в серединах сторон симплексов (схема 2) [8], стабилизированного метода конечных элементов Петрова — Галёркина (SUPG) [4], метода конечных элементов/конечных объемов (МКЭ/МКО). В табл. 1 и 2 приведены значения евклидовой

Таблица 1. Значения гт, полученные на сетке 41 х 41

VT SUPG МКЭ/ МКО Схема 1 МКЭ Схема 2

1 7.32E-5 2.09E-5 4.16E-5 1.44E-4 4.62E-5

1.0E-1 2.00E-4 2.06E-4 3.77E-4 1.42E-3 3.89E-4

1.0E-2 3.55E-5 1.05E-3 1.39E-3 7.56E-3 1.15E-3

1.0E-3 6.36E-6 1.54E-3 1.69E-3 1.11E-2 1.36E-3

Таблица 2. Значения гт, полученные на сетке 81 х 81

VT SUPG МКЭ/ МКО Схема 1 МКЭ Схема 2

1 2.14E-5 1.05E-5 1.38E-5 7.41E-5 7.21E-5

1.0E-1 6.70E-5 1.04E-4 1.29E-4 7.31E-4 7.17E-4

1.0E-2 1.96E-5 5.39E-4 5.12E-4 3.88E-3 5.75E-4

1.0E-3 3.03E-6 7.80E-4 6.46E-4 5.73E-3 5.65E-4

Таблица 3. Исследование интерполяционных свойств метода ЯиРС/РЯРС на уравнениях Эйлера

N U — Uh UX Ux U — Uh Uy Uy Uz — Uh ■e 1 n

3115 3.42E-3 3.44E-3 1.17E-3 3.23E-2 77

7041 2.03E-3 2.03E-3 8.48E-4 2.52E-2 130

9108 1.76E-3 1.76E-3 7.93E-4 2.46E-2 260

нормы вектора погрешности Гу = ||Т — Th|| для различных значений vy. При доминировании конвективных членов МКЭ дает наихудший результат. Метод SUPG показывает хорошие результаты при уменьшении значений vy.

Численно решаются уравнения Эйлера в области П G R3, задача имеет точное решение для компонент вектора скорости

sin (пх) ■ cos (ny) cos (nz)

U = I cos (пх) sin (ny) cos (nz)

— 2 cos (пх) cos (ny) sin (nz)

и давления p = 3п cos (пх) cos (ny) cos (nz).

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

4.2. Сравнение схем методов конечных элементов и конечных разностей

Проведем сравнительный анализ описанных вычислительных схем методов конечных элементов и конечных разностей для задачи о плоском существенно вихревом течении несжимаемой жидкости в канале с разворотом потока на 270° с известным точным решением.

Граница канала состоит из прямолинейных отрезков Г и Г2:

Г1 = 1(x,y)

3п 3п п

------< х<-------------, y = —

8 “ “ 16 У 2

п 11п ^ ^ 7п к

х = 0, -----------< y < — > ,

16 8

Г2 = |(х,у)

и криволинейных участков Г0 и Г0, заданных в параметрическом виде

х = /1(q^ y = /2(q), 0 < q < 1

где

eW \ I 3 Л f2( \ -I cos х0

/ (q) = х0 cos п -q — 1 , / (q) = arcsin

2 / \cosf 1(q)

х0 = 3п/8 для Г0 и х0 = 3п/16 для Г0;.

В качестве точного решения задачи (1)-(3) для данной области П возьмем функции

Uexact = Ul (х, y) = cos х cos y, Vexact = U2 (х, y) = sin х sin y,

Рис. 1. Криволинейная сетка, построенная при Рис. 2. Конечно-элементная сетка (466 узлов и

управляющей функции ад = 1 + а^ехао^, раз- 796 треугольников).

мерность сетки 31 х 11.

Тогда

Рехас! = р(х, у) = Ро 1 (сОв 2х — СОв 2у).

ФехасЬ = Ф(х, у) = СОв х від у, ^ехас! = ^(х, у) = 2 СОв X вІП у.

Сетка строится методом эквираспределения [9] и адаптируется к решению путем задания управляющей функции ,ш в виде ,ш = 1 + а|^ехасЬ | (рис. 1), конечно-элементная сетка представлена на рис. 2.

В табл. 4 показана точность решения задачи всеми описанными методами, где Рг, і = 0,1, 2, — полиномы нулевого, первого и второго порядка для аппроксимации компонент вектора скорости и давления.

Для моделирования течения в криволинейном канале с разворотом потока на 180° [10, 11] рассмотрим уравнения Эйлера в области П С И3 с краевыми условиями вида

и|

хЄГ і

^(х) и ■ п|хЄГо = 0, и ■ П|ХЄГ2 = ^2 (х) ,

где п — внешняя нормаль к границе области Г = Г и Го и Г2; Г1 — входное отверстие в области П; VI • п < 0 на Г1, Г0 — непроницаемая часть границы, Г2 — выходное отверстие; У2 > 0 на Г2 (VI = (0,1), У2 = 1).

Таблица 4. Исследование методов моделирования течения несжимаемой жидкости в канале с разворотом потока на 270°

Метод 1 «и 1 «1

Метод конечных разностей 1.95Е-2 1.96Е-2

ВиРО/РЯГО (Р1 - Р1) 2.11Е-2 2.08Е-2

Смешанный МКЭ (Р2 — Р:) 2.25Е-2 2.24Е-2

Смешанный МКЭ (Р1 — Р0) 2.29Е-2 2.31Е-2

Таблица 5. Значения дискретной дивергенции (ББ1У), полученные при расчете трехмерного течения в криволинейном канале

N Стабилизированный МКЭ, Смешанный МКЭ,

ББ1У ББ1У

1188 0.183Е-3 0.198Е-3

3078 0.105Е-3 0.119Е-3

Рис. 3. Значения компоненты вектора скорости иу в правом изгибе канала: 1 — результаты, полученные методом конечных разностей, 2 — смешанным МКЭ, 3 — стабилизированным МКЭ.

Проведено моделирование на измельчающихся сетках стабилизированным методом БИРО/РБРО, смешанным методом конечных элементов и методом конечных разностей с использованием переменных вектор-потенциал — вектор завихренности. Вектор скорости направлен перпендикулярно к входной границе Г! канала. В табл. 5 приведены значения дискретной дивергенции. При расчете стабилизированным методом конечных элементов значения дискретной дивергенции меньше по сравнению со значениями, полученными смешанным МКЭ. На рис. 3 даны результаты расчета (значения компоненты вектора скорости Иу в правом изгибе канала), полученные при моделировании течения стабилизированным МКЭ, смешанным МКЭ и методом конечных разностей [10]. Численные значения компонент скорости, полученные этими методами, отличаются друг от друга на величину 10-3.

Заключение

Анализ вычислительных схем моделирования течений жидкости позволяет сделать следующие выводы.

Сравнение данных, полученных при моделировании двумерных уравнений конвекции-диффузии с доминированием конвективного переноса, говорит о преимуществе метода БИРО и рассмотренных противопотоковых схем метода конечных элементов над методами МКЭ/МКО, стандартным МКЭ.

При моделировании поля течения в 2D смешанным МКЭ, алгоритмом Удзавы на базисных функциях Тейлора —Худа, удовлетворяющих условию ЛББ, более точное решение достигается на квадратичных базисных функциях для скорости и линейных — для давления (P2 — P:). Стабилизированные методы, в которых интерполяционные пространства скорости и давления одного порядка (линейные), позволяют более точно вычислять давление по сравнению с рассмотренными методами на базисных функциях Тейлора —Худа. Использование функций P2 — P1 приводит к увеличению размерности систем линейных алгебраических уравнений (СЛАУ) по сравнению с размерностью СЛАУ стабилизированных методов на величину, равную числу узлов, лежащих на серединах сторон конечных элементов расчетной области.

Проведено моделирование течения в криволинейных каналах в двумерном и трехмерном случаях стабилизированным МКЭ, смешанным МКЭ и методом конечных разностей. Результаты расчетов, полученные рассмотренными методами, отличаются друг от друга на величину 10-3.

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

[1] Warsi Z.U.A. Fluid dynamics. Theoretical and Computational Approaches. FL: Taylor & Francis, CRC Press, 2005.

[2] Liu C. Ho. Numerical solutions of three-dimensional Navier — Stokes equations by a velocity — vorticity method // Intern. J. Numer. Meth. Fluids. 2001. N 35. P. 533-557.

[3] Годунов С.К., РявЕнький В.С. Разностные схемы. М.: Наука, 1977.

[4] Tezduyar T. Stabilization parameters and local length scales in SUPG and PSPG formulations // Fifth World Congress on Computational Mechanics, July 7-12, 2002, Vienna, Austria. www.mems.rice.edu/tezduyar.

[5] Brezzi F., Fortin M. Mixed and Hybrid Finite Element Methods. N.Y.: Springer-Verlag, 1991.

[6] Темам Р. Уравнения Навье — Стокса. Теория и численный анализ. М.: Мир, 1981.

[7] ШокинА Н.Ю. Численное моделирование на адаптивных сетках двумерных установившихся течений жидкости и газа // Вычисл. технологии. 1998. T. 3, № 3. С. 85-93.

[8] Harig J. A 3-D Finite Element Upwind Approximation of the Stationary Navier — Stokes Equations. Univ. Heidelberg. 1992. Preprint 92-09.

[9] Khakimzyanoy G.S., Shokina N.Yu. Equidistribution method for the construction of adaptive grids // Russ. J. Numer. Anal. Math. Model. 1999. Vol. 14, N 4. P. 339-358.

[10] Хакимзянов Г.С., ШокинА Н.Ю., Кюстер У. Численное моделирование трехмерных течений идеальной жидкости с использованием векторного потенциала и вектора вихря // Вычисл. технологии. 2004. Т. 9. (Совмест. вып. по матер. Междунар. конф. “Вычислительные технологии и математическое моделирование в науке, технике и образовании”). Вест. КазНУ. 2004. № 3(42). Ч. IV. С. 168-178.

[11] Говыш А.В., ШокинА Н.Ю., Шурина Э.П. Анализ вычислительных алгоритмов для уравнений Навье — Стокса, основанных на методе конечных элементов // Тр. Междунар. конф. по вычислительной математике, МКВМ-2004 / Под ред. Г.А. Михайлова, В.П. Ильина, Ю.М. Лаевского. Новосибирск: Изд-во ИВМиМГ СО РАН. 2004. Ч. 1. С. 467-472.

Поступила в редакцию 31 октября 2006 г.