Анализ алгоритмов методов конечных элементов и конечного объема на неортогональных сетках при решении уравнений Навье - Стокса Текст научной статьи по специальности «Физика»

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

Том 2, № 4, 1997

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

Э. П. Шурина, Т. В. Войтович Новосибирский государственный технический университет, Россия

e-mail: [email protected]

The finite element method in the Ritz' form (discretization of momentum equation) in combination with the finite volume technique (approximation of continuity equation) are used for incompressible viscous flows in pipes modeling. Analysis of proposed linearization scheme connecting with transition of the momentum equations to a self-adjoined form is presented. One of the possible finite volume form constructed on pipe crossection triangulars and corresponding approximation of continuity equation are investigated.

1. Введение

При моделировании вязких течений несжимаемой жидкости в каналах сложных геометрий возникает задача исследования различных способов дискретизации системы уравнений Навье — Стокса. В данной работе приводятся результаты анализа алгоритма, основанного на методах конечных элементов (МКЭ) и конечного объема (МКО), использующих неструктурированные сетки (триангуляция в случае 2D, разбиение на тетраэдры в случае 3D). Построение сетки, удовлетворяющей ряду условий (отсутствие малых углов, "вырождающихся" треугольников, смежных треугольников со значительно различающимися размерами), в заданной расчетной области представляет собой самостоятельную проблему и здесь рассматриваться не будет.

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

© Э.П. Шурина, Т.В. Войтович, 1997.

Основной сложностью при расчете поля скорости является вычисление поля давления. Для решения этой проблемы, как правило, используют уравнение неразрывности. Известен ряд методов: типа вихрь—функция тока; искусственной сжимаемости; метод Харлоу с расщеплением по времени и т.д. [1-3]. Сравнительный анализ методов позволил отдать предпочтение последовательной во времени коррекции полей скоростей—давления, характерной для алгоритмов типа SIMPLE (SIMPLER) [4], где поле давления вычисляется итерационно из условия выполнения уравнения неразрывности. Для дискретизации уравнения неразрывности был выбран метод конечного объема (бокс-метод) [4], ориентированный на применение ортогональных сеток. При моделировании течений в каналах со сложным поперечным сечением возникает необходимость использования неструктурированных сеток [5]. В работе предлагается способ построения конечного объема на основании треугольных призм, заполняющих пространство между сечениями, и соответствующий способ аппроксимации уравнения неразрывности.

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

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

Математическая модель течения в канале (область П = П и Г), содержащая уравнения Навье — Стокса в приближении пограничного слоя [6], имеет вид: уравнения сохранения импульса

TTdU лгди ВЦ' UТГ + + WTT

дх ду dz

ттдУ лгдУ bY UlT + VlT + Wtt

дх ду dz

TTdW dW BW p U— + V—- + W—-дх ду dz

уравнение неразрывности

д

dp + дх

dp дУ

dp

dz

A ( d£\ + dUL\

дх \ dx J ду \ dy J

d ( dV

dx dx

d ( dW\ d

A ( dV\

дУ V дУ J dW

дх дх ) + ду \ dy

+ pFi, + PF2, + pF3;

дд - (pU) + - (py) + ^ (pw )

(i) (2) (3)

(4)

где и = (и, V, Ш) —векторное поле скорости, р — давление, р — плотность, —

поле внешних массовых сил, ^ — молекулярная вязкость. На стенках канала заданы условия прилипания:

U |г = V |г = W |г = 0.

(5)

На начальном сечении канала заданы распределения компонент вектора скорости и давления:

и и = ио, V и = Ус, Ш и = Шо, ри = Ро. (6)

p

p

0

3. Эквивалентная вариационная постановка в форме Рит-ца для уравнений движения

Пусть ось канала разбита точками zi,...,ZNZ и соответствующие z = Zk, k = 1,...,NZ, сечения обозначены как Ei,..., Enz (расчетная область схематично представлена на рис. 1).

Рис. 1.

Предположим, что триангуляция сечения Ек— задана и переносится с сохранением числа треугольников и их подобия на сечение Ек.

Построим дискретную модель уравнения движения для компоненты скорости и (в предположении отсутствия внешних массовых сил):

Р

TTdü лгди тт-rdU

UTT + VJT + W-^-

dx dy dz

dp dx

А ( dü\ dL ( d_U_\

dx \ dx ) dy \ dy )

(7)

Применение чисто неявной разностной схемы по £ для аппроксимации слагаемого

тхди ик — ик-1 Л

р№ —— & рм---, Агк = хк — хк-1 дает уравнение

дх Агк

д

3U\

dx у dx )

А (dü\

dy V dy )

+Р

TTdük TrdUk TTr Uk U—--h V—--h W-

dx

dy

Azk

= dp+pWU k-i

dx Azk

(8)

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

Переносом слагаемых, соответствующих конвективным членам, в правую часть (8) получим (в операторной форме)

L(U) = fu(U, V, W,p),

(9)

d f dü\ d ( 3U\ pW

L(U) = + a • U, a = t—•

dx \ dx j dy \ dy) Azk

fu (U,V,W,p) = - dx - Р

u— + V9-U

dx dy

+ pW■

U

k-i

Azk

(10)

Аналогично уравнения сохранения импульса по у, г могут быть приведены к виду

Ь(У) = ¡и(и, V, Ш,р), Ь(Ш) = ¡щ(и, V, Ш,р). (11)

Самосопряженность оператора Ь(и) позволяет выписать эквивалентную вариационную постановку в форме Ритца (сильная форма), согласно которой решение, например, исходной дифференциальной задачи (9) с краевыми условиями (5) на сечении ^к эквивалентно минимизации квадратичного функционала вида

1и(у) = (ь(у),у) - 2(^,и)

( ди

"Ш + »[ву

ду\'

¿П + у ау2(1П - 2 у ¡иу<Ю.

£к £к

(12)

на множестве допустимых функций и £ Н1; здесь а и ¡и определены выражениями (10). Функционалы, соответствующие У и Ш, аналогичны. Нижний индекс и в обозначении полученного функционала 1и(и) означает, что он вычисляется с использованием параметров а и ¡и, зависящих не от решения (и,У,Ш) задачи (9), (11), а от некоторого вектора и = (и' ,У', Ш'), каждая компонента которого есть функция из Н1.

Развязывая нелинейные по коэффициенту и правой части уравнения (9), (11) относительно неизвестных функций и, У, Ш [7], будем учитывать нелинейность по конвекции при помощи итераций (теоремы 1, 2 [8]). На 5-й итерации по нелинейности в этом случае будет производиться минимизация функционалов

1ив-! (у) = (ь(у),у) - 20и ,и), (у) = (ь(у),у) - 2(Л,и), ^(у) = (ь(у),у) - 2(^,и),

2

в которых правые части ¡и, ¡у, ¡щ и коэффициент а содержат компоненты вектора скорости, вычисленные при минимизации функционалов 1и-2 (и), 1у3-2 (и), (и). Если построенный таким образом процесс сходится в первой норме, то он сходится к истинному решению нелинейной задачи (9), (11) на сечении ^^.

4. Дискретная модель уравнений сохранения импульса

Для дискретизации трех линейных уравнений вида

Ь(ик'3-1) = ¡и (ик'3-1,У к'3-1,Ш к'3-1,рк),

1(ук,8-1 ) = ¡и (ик,8-\, у к,8-1,Ш к,3-1,рк ),

Ь(Шк'3-1) = ¡Щ (ик'3-1,Ук'3-1 к'3-1,рк) (13)

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

конечномерных подпространств Бсодержащихся в НЕ, элементы которых представимы в виде линейной комбинации финитных базисных функций фг:

= ^ О-гФг = ^Ф-

г=1

Так как процедура коррекции полей скоростей—давления имеет итерационный характер, на итерации I этой процедуры давление находится как рк'1 = рк'1-1 + ДрМ-1 [6]. Тогда дискретным аналогом уравнений (13) на сечении Ек являются матричные уравнения вида

Аи = Фи - С1Арк'1-1,

АУ = Фу - С2Арк1-\

АШ = Фш - - С3Дрм-1

(14)

Векторы Фи, Фу, Фш и матрица массы, входящая в А, пересчитываются на каждой итерации по нелинейности. В формировании матриц 01, 02, аппроксимирующих операции дифференцирования ——,——, в этом случае участвуют интегралы вида [ -Ф— [ -Ф— дх ду J дх J ду

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

Матрица СЛАУ А и матрицы для учета поправки к давлению Ог имеют разреженную структуру, что делает возможным использовать удобную схему хранения — разреженный строчный формат. Для упаковки всех несимметричных матриц Ог можно использовать один портрет.

Так как осевая скорость может быстро изменяться на конечном элементе, например в пристеночных областях, принималось, что коэффициент а = рШкув-1/Агк представлен своим линейным интерполянтом на каждом из треугольников; аналогичное представление было использовано для давления.

5. Дискретная модель уравнения неразрывности

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

Метод конечного объема традиционно строится с использованием ортогональных сеток, поскольку в этом случае значительно упрощаются выражения для потоков искомых величин через грани конечных объемов и их легче аппроксимировать [4]. Так как для дискретизации сечений в методе конечных элементов были выбраны неортогональные сетки, естественным оказывается построение конечных объемов с использованием триангуляции сечений.

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

Определим вид конечного объема для рассматриваемой задачи. Пусть триангуляция сечения Ек—1 переносится с сохранением числа треугольников и их подобия на сечение Ек. Воспользуемся тем, что в результате этой процедуры пространство Пк, заключенное между сечениями Ек-1 и Ек, оказывается разбитым на N призм. Конечный объем, соответствующий узлу г, может быть получен переносом по г отрезков медиан прилежащих к узлу г треугольников (рис. 2, а). Если г — внутренняя точка триангуляции, будем строить "неполный" конечный объем (рис. 2, б). Примем следующие обозначения: Пг — конечный объем, соответствующий узлу г; (г,],т) — один из треугольников;

Е^к-1, Ек — основания конечного объема (грани, принадлежащие сечениям Ек_1, Ек); Ет — грань конечного объема, основаниями которой являются отрезки медиан треугольников номер п, выходящих из узла т;

Б'т — отрезок медианы треугольника номер п, выходящей из узла т.

Рис. 2. Способ построения контрольного объема. Построение дискретной модели уравнения неразрывности методом конеч-

ного объема.

Используем метод конечного объема для аппроксимации уравнения неразрывности

д д д

дХ (Ри) + - (РУ) + - ) = „. (15)

Фиксируем некоторую внутреннюю точку триангуляции г. Проинтегрируем уравнение (15) по соответствующему конечному объему П с использованием формулы Остроградского — Гаусса:

J (т = J гг(р[1)(а,

П а

здесь и — векторное поле, П — произвольный объем, а — его поверхность, п — орт внешней нормали к поверхности. Балансовое соотношение, формализующее закон сохранения массы для конечного объема ПДсм. рис. 2, а), будет иметь вид:

J рШ(а - J рШ(а + ^ ^ J р(ипх + Упу)(а = 0, (16)

Як-1 п пЬ ЯПь

где п — номера треугольников, встречающихся в узле г, пЬ — номера узлов треугольника номер п, смежных с узлом г, п = (пх, пу, 0) — внешняя нормаль к боковой грани ЕПь. Два первых слагаемых в (16) представляют собой потоки вектора скорости через основания конечного объема Е^-1, Е^, а третье — поток через боковые грани.

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

Рассмотрим некоторый треугольник (г,],т), содержащий узел г (рис. 3).

Рис. 3.

Вклад треугольника (г,],т) в выражение для потока через общую боковую поверхность П содержит потоки через две боковые грани Е'ПЬ, пЬ £ {],т}. Выражение для каждого из этих потоков в предположении несжимаемости жидкости может быть преобразовано следующим образом:

^ППь = Р / (и(х,У,г)пх + У(х,у,г)пу)(а = р / (и(х,у,г)пх + У(х,у,г)пу)(а =

пЬ

пЬ

Zk

xnb

УпЬ

Zk

р f <— f V(x,y(x),z)dx + J U(x(y),y, z)dy > dz = p f {—Iv(z) + IJ(z)} dz. (17)

Zk-1

Ус

Zk-1

Здесь и далее использованы следующие обозначения для отрезка медианы Б^,: (хпЬ1 упЬ) — координаты основания медианы Б^; (хс,ус) — координаты центроида треугольника п; ¡ПЬ, 1П — длины проекций медианы БП, на оси х, у соответственно; уравнение медианы: у = кпЬх + ЬпЬ. Используя по г квадратурную формулу, получим аппроксимацию потока Щ,:

Zk

.1щь = рУ {-1у (г) + 1и(г)} dz и

Zk-1

и р {ш1у(гк) + ш!ц(гк) - (1 - ш)1у (гк-1) + (1 - ш)1ц(гк-1)} (Агк),

(18)

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

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

Vi av 1 xi yi

Vj = Bn bv , Bn = 1 x2 y2

Vm . °v . 1 x3 y3 _

К = (Вп )-1,

получим выражение для 1у (номер сечения временно опустим):

хпЬ

1у = [ ау + Ьух + Су [кпЬх + ЬпЬ] dx = {(ау + суЬпЬ)1пЬ + (Ьу + сукпЬ)0.5 Р)2} Sign|Bn|:

= {(kiiV + kuVj + ki2Vm + bn k;nVi + bnbk32V + bnbk33Vm)C + + k2lVl + k22V+ +k2aVm + knbhiV + knbk32Vj + knbk33Vm)0.5 (lf)2}Sign|Bn| = dV?bv, где v = \Vi,Vj, Vm]T — локальный вектор,

n.n dV

¡пЬк1и + 0.5 (¡пЬ) к2и + (упЬ1ХпЬ - хпЬ1пЬ + 0.5/пЬС)кз^| Sign|Bп| Аналогично можно получить

3

v=1

Ij = djn и,

U =[Ui,Uj ,Um]

T

n.n

dJ

ffkiv + 0.5 (if)2 k3v + (xnbin — ynbinb + 0.5/nbC)k2v) Sign|Bn|

3

v=i

(19)

Таким образом, на каждом из сечений Ек—, Ек определены интегралы, входящие в (18) (для локальных векторов и и и восстановим номера сечений):

4 = (>ь ик, /ки-1 = (^ьик-1,

/к = (^пЬУк, /к-1 = 4ПЬик-1. (20)

В результате аппроксимация потока через боковую грань с учетом (20) имеет вид:

!Пь = ! р(ипх + Упу)йа «

Яп

пЬ

р {-ийупЬик + ийипЬик - (1 - и)йупЬик-1 + (1 - и)йипЬик-1} (Д^), (21)

7 _ Г • Ч п.пЬ п.пЬ

где пЬ £ \3,т\; здесь векторы аи , йу содержат вклады в элементы г-х строк соответствующих матриц дискретного аналога уравнения (15).

Переход к генерации матриц дискретного аналога по конечным элементам.

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

Рис. 4.

Любой треугольник является основанием частей трех конечных объемов (рис. 4). Законы сохранения массы, записанные в интегральной форме, для конечных объемов г,], т дадут коэффициенты в строках г,], т соответственно. Таким образом, на каждом из сечений Ек-1, Ек вклад треугольного конечного элемента (г],т) в матрицы дискретных аналогов может быть представлен в виде таких локальных матриц коэффициентов , Пу *, П^) , использование которых с (21) позволяет записать

П(е) = Пи =

а £

и

и

ат

и

П(е) = Пу =

а й{

V

V

йт

V

здесь

d'U = J

Е

il,nb lU

nbE{i,j,m}\{l}

J =

dlV = J

E

l,nb V

nbE{i,j,m}\{l}

pjAzk для сечения ,

p(1 — j)Azk для сечения Ek-1.

l E {i,j, m}

Итак, получены формулы для локальных матриц D^, D^ вклада элемента (i,j,m) в балансовые соотношения, матрица будет определена ниже.

Учет ориентации поверхности. Коэффициент Sign|Bn| в выражениях для векторов d'Unb, dVnb возникает следующим образом. Пусть для конечного объема в плоскости любого из сечений задано положительное направление обхода (рис. 5). Тогда при расчете вкладов треугольника (i,j, m) поток через отрезок медианы f p(Unx + Vny)ds, nb E {j, m} ,

nb

nb E {j, m} входит в балансовые соотношения со знаком "+", если точки i, nb, {j, m}\{nb} образуют левый поворот, или " —" — в противоположном случае, т.е. необходимо вычислить знак определителя вида

det

1 Xi yi

1 xnb ynb

1 x{j,m}\{nb} y{j,m}\{nb}

Рис. 5.

В случае использования треугольных сеток оказались возможными только два случая (рис. 6); здесь использована локальная нумерация узлов.

Рис. 6. а, б — точки 1-3 образуют соответственно левый и правый поворот.

Интегрирование производится от центра масс до основания медианы, тогда интегралы по граням конечных объемов вносятся в локальные матрицы с чередующимися знаками, а именно (фиксируем узлы по порядку, оставшиеся узлы обходим по возрастанию):

a) sn +; on S3 —; b) n s2 s3n +;

a) sn +; b) sn +; s3n

a) Sn +; n s2 —; b) n s1 —; s2n +

Это свойство позволяет компактно записать схему вычисления локальных матриц с учетом ориентации с использованием умножения на знак определителя det Bn, а так как во все выражения для локальных матриц входит величина (det Bn) 1 (поскольку K = (Bn)-1), то

Sign det Bn 1

можно использовать --- = —--, и единственной дополнительной операцией,

det Bn |det Bn|

необходимой для учета ориентации поверхности, будет взятие модуля.

Для завершения построения аппроксимации уравнения неразрывности определим локальную матрицу dWW . Аппроксимация потоков через основания конечного объема Q имеет вид:

У pWda « pWjmes Elz, l e{k,k - 1} ,

si

а компоненты локальной матрицы (отрезки медиан делят треугольник на части равной площади) могут быть представлены соотношениями

D

vk

а/ 0,

3pWimes (i,j,m), v = k,

v = k.

Процесс вычисления локальных матриц элемента (i,j,m) не зависит от наличия в нем узлов, принадлежащих границе Г.

Суммируя по всем конечным элементам расширенные локальные матрицы, получим на каждом из сечений (I — номер сечения):

DU = Е D^, Dlv = £ DDlw = £ D^1. (e)

(e)

(e)

В результате аппроксимация уравнения неразрывности в пространстве между сечениями Ек-1, Ек будет иметь вид:

DkUk + Dk Uk-1 + Dk Vk + Dk Vk-1 + Dk Wk - Dk Wk-1 = 0.

(22)

Подобно работе [6], используем совместно аппроксимации уравнений движения и уравнения неразрывности для построения итерационной процедуры коррекции полей скоростей— давления. Исключение рассчитываемых на сечении Ек значений векторов скорости в узлах ик, Ук, Шк из уравнения (22) и дискретных аналогов уравнений движения (14) дает СЛАУ относительно поправки к давлению:

ПДр = Б,

П = DU A-1Gi + Dkv A-1G2 + DW A-1GS,

Б = Dkv А-1Фи + DV А-1Фу + DW А-1Фщ + Dk-1Uk-1 + Dkv-1Vk-1 + D^1W

)yA-1'

k -1

k-1T Tk-1

■¡k-U rk-1

ïk-1TT rk-1

Матрица П системы невырождена, но плохо обусловлена из-за формального исключения давления из числа неизвестных предположением pk,l,m = pk,i,m-1 + Apk,l,m-1, и задача вычисления поля давления в случае течения несжимаемой жидкости в канале становится некорректной. Для решения СЛАУ (23) использовалась регуляризация [6], т.е. находилось решение Ape, доставляющее минимум функционалу

Q(Ap) = ||ПАр — Б||2 + е ||Ap||2.

Q(Ap)

В этом случае из необходимого условия экстремума ——-— = 0 получается СЛАУ для

д Ap

определения регуляризированного решения Ap£:

(Пт П + еЕ )Ap£ = Пт Б. (24)

Выбор параметра регуляризации е может оказать значительное влияние на скорость сходимости итераций, но не может изменить поле давления, к которому сходится итерационный процесс. Для алгоритмов типа SIMPLE показано [4], что поскольку в сходящемся решении достигаются такие поля давления, при которых соответствующие поля скорости удовлетворяют уравнению неразрывности, то детали построения уравнения для Ap не влияют на корректность сходящегося решения, изменяя лишь скорость сходимости метода. Для ускорения сходимости первая поправка к давлению Apk'° на каждом из сечений вычислялась из условия постоянства расхода через поперечное сечение канала, как это было предложено в [6]; при этом для расчета расхода использовалась линейная интерполяция продольной составляющей на каждом из конечных элементов.

Окончательно алгоритм моделирования течения в канале с использованием сочетания методов конечных элементов и конечного объема может быть записан в виде: Ввод параметров алгоритма и определение начальных условий (Uo,Vo,Wo,po) триангуляция первого из сечений цикл по сечениям: do k = 1,NZ

генерация матриц Gi,G2,G3, и вектора численного интегрирования Rk итерации по усреднению поля давления:

do l = 1, max L

начальное приближение: Apk,1,° = Apk'1-1 цикл по коррекции полей скоростей-давления:

do m = 1, max M

итерации по нелинейности уравений Навье — Стокса:

do it = 1, max NL

решение уравнений движения при фиксированном поле давления pk'l,m-1: генерация и решение СЛАУ

AU k's = Фи — G1Apk'l'm-1 AV k's = Ф[ — G2Apk'l'm-1 AW k's = Ф№ — G3Apk'l'm-1 if||Uk's — Uk's-1\\ <£U и \\Vk's — Vk's-1\\ <£V и

|| Wk's — Wk,s-1 \\ < £U-переход к вычислению Apk'l'm

enddo

вычисление поправки к давлению:

if[первая итерация по уточнению поля давления на сечении] then нахождение Apk'0 из расходного соотношения

else генерация матриц дискретного аналога уравнения неразрывности Du,Dy,Dw, СЛАУ ПApke'l'm _ B,

л k,l,m

нахождение Ape регуляризациеи:

(Пт П + eE)Apk'l'm = ПТ B. if[||nTB|| < £Р]

then переход к следующему сечению

else коррекция поля давления:

pk,l,m _pk,l,m-1 + Apk'l'm

enddo

pk,i _ epk,i-1 + (i _ e)pk'l'm

enddo enddo

Так как метод конечного объема был построен по аналогии с методом конечных элементов, соответствующие матрицы, входящие в дискретный аналог уравнения неразрывности, также собираются поэлементно, а их портрет полностью совпадает с портретом несимметричных матриц, используемых для генерации векторов правых частей в методе Ритца, — G1, G2, G3. Таким образом, для хранения матриц как метода конечных элементов, так и метода конечного объема возможно использование одних и тех же массивов указателей для ненулевых элементов матриц, что позволяет уменьшить требования к оперативной памяти. При выборе метода хранения симметричной разреженной матрицы метода конечных элементов A предпочтение было отдано способу представления RR(U)U [9], при котором диагональ хранится отдельно от матрицы, так как данная матрица не имеет нулевых диагональных элементов, а такие матричные операции, как неполное разложение Холесского, содержат частые обращения к диагональным элементам.

При выборе метода решения СЛАУ проводилось сравнение методов сопряженных градиентов, сопряженных градиентов с предобусловливанием неполной факторизацией Хо-лесского с извлечением квадратного корня и предобусловленного метода сопряженных градиентов с неполной факторизацией Холесского без извлечения корня [10]. Последние были выбраны, поскольку реализация алгоритма показала плохую обусловленность матрицы СЛАУ относительно поправки к давлению. Сравнение методов по скорости, точности и требованиям к оперативной памяти показало, что предпочтение можно отдать методу сопряженных градиентов с предобусловливанием при помощи разложения Холесского без извлечения корня в сочетании со схемой хранения RR(U)U.

6. Результаты расчета и особенности алгоритма

В качестве тестовой рассматривалась задача расчета поля течения на входном участке гладкого прямоугольного канала с различным отношением сторон и равномерным распределением потока на входе. При моделировании варьировались числа Рейнольдса, при этом исследовалось влияние параметров алгоритма на скорость сходимости. На первом из сечений задавались нулевые поперечные составляющие скорости U0 = V0 = 0 и удовлетворяющая расходному соотношению J pWda = q продольная составляющая; принималось,

что

W1 = W0 = const во внутренних точках триангуляции, W 1|г = 0. Результатами расчета являлись компоненты вектора скорости U,V,W, значение давления p в каждой точке триангуляции на каждом из сечений Ek, массив, характеризующий падение давления от

Ртах " Ргшп = 9.40997е - 007

Ртах - Ргшп = 9.18962е - 009

[/•Ю-9, м/с

-и-1

0.5

о

-0.5

В-и-■

О 0.01 ж, м 0.02

Ртах - Ршт = 6.31154в - 008

Ртах " Ртт = 4.76104в - 009

Рис. 7. Профили составляющих векторов скорости , и и поля давлений (масштаб вывода различен) на разных сечениях канала. Ие = 19.

Ш, м/с

У- 10"6, м/с

Рис. 8. Профили составляющих векторов скорости Ие = 190, б — Ие = 380.

б

Ш, м/с

и- 10"6, м/с

У- 10"6, м/с

и, V на разных сечениях канала. а —

Рис. 9. Профили составляющих векторов скорости W, U и изменение вида поперечной составляющей U на разных сечениях канала; масштаб вывода U(x, y) различен. Re = 800.

сечения к сечению (значения р(0, 0, гк)), значение нормы ||ПТВ|| на итерациях по уточнению полей скоростей-давления.

На рис. 7-10 приведены результаты расчета в канале с поперечным сечением х = 2 см, у = 1 см. Числа Рейнольдса входного потока варьировались от 19 до 800; на графиках представлены профили составляющих вектора скорости !],У,Ш в сечениях плоскостью у = 0.5, х = 1. Расчеты для И,е = 19 проводились с наиболее подробным шагом по г (Дг = 1 • 10-4); при И,е = 190, 380, 800 шаг по г фиксировался и составлял Дг = 1 • 10-2.

Для получаемых в результате расчета профилей поперечной составляющей скорости (см. рис. 2, 3) на первых сечениях канала характерно наличие максимумов вблизи стенок. По мере развития течение отходит от стенки и максимум продольной составляющей перемещается к центру канала, что наблюдается на рис. 7-9.

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

Рассчитываемое поле давления, как правило, имеет максимум вблизи стенок (рис. 10) (приводимый на рисунках параметр ртах — рт;п характеризует перепад давления на сечении). Как правило, уменьшение параметра ртах — ртт от сечения к сечению таково, что при графическом выводе полей давления на двух соседних сечениях с одинаковым масштабом последний из графиков будет казаться плоскостью. Вследствие этого масштаб вывода значений давления для соседних сечений различен.

Рис. 10. Характерный вид поля давления на двух последовательных сечениях. Ие = 190.

Поскольку в систему уравнений Навье — Стокса входит градиент давления, то оно имеет относительное, с точностью до постоянной, значение, и интересен лишь перепад давления. Поэтому его необходимо фиксировать в одной из точек триангуляции. При расчетах принималось, что после вычисления первой поправки к Дрк'° из расходного соотношения компонента последующих поправок, соответствующая точке (0, 0,гк), равна нулю

(учет этого условия производился модификацией первой строки матрицы ПтП + еЕ СЛАУ относительно поправки Дре).

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

Было исследовано влияние параметра регуляризации е на скорость сходимости итераций по уточнению давления. При этом чрезмерное увеличение параметра регуляризации е не влияет на качественное поведение поправки к давлению, но значительно уменьшает каждую из компонент вектора Др£, что очень замедляет сходимость (рис. 11, а). На первых итерациях в этом случае могут наблюдаться осцилляции нормы ||ПТ В|| , малость которой является условием выхода.

Число итераций

Рис. 11. Влияние параметра регуляризации на скорость сходимости алгоритма.

Зависимость числа итераций по уточнению полей скорости—давления от параметра регуляризации всегда имеет четко выраженный экстремум. Так, на примере расчета И,е = 380 (рис. 11, б), характеризующего сходимость на втором сечении, можно заметить первоначальное снижение числа итераций с 20 до 3 при уменьшении е, 1 • 10-3 до 5 • 10-5. Уменьшению е до 1.75 • 10-5 соответствует увеличение числа итераций до четырех.

Дальнейшее уменьшение е приводит к неопределенной СЛАУ ПтП + еЕ, для которой неполная факторизация Холесского невыполнима, и получить численное решение становится невозможным.

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

несколько завышенного значения е.

Было также проведено исследование предложенной схемы линеаризации уравнений движения, включающей пересчет матрицы массы и вектора правой части. В любом случае при моделировании поля течения с использованием маршевого прохода по области наиболее трудоемким является расчет входного участка, где течение еще не установилось. Так, второму сечению соответствовало максимальное количество итераций по нелинейности уравнений Навье — Стокса (от 10 до 16 итераций для канала 2 х 1 см), в то время как уже на третьем сечении оно составляло, как правило, 4-5 итераций и с установлением уменьшалось до 1-2 итераций по нелинейности. С изменением формы канала это соотношение может изменяться; так, при расчете "щелевидного" канала 4 х 1 см расчет первого профиля потребовал 26 итераций, но уже через несколько сечений их число уменьшилось до двух.

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

Цшк'8 — Шм-1|| < 1 • 10-18, ||им — ик'8-1Ц + — Vк'®-11| < 1 • 10-24. (25)

Как видно из рис. 12, трех—четырех итераций достаточно для достижения практически полного совпадения графиков на соседних итерациях по нелинейности. Этот факт позволяет не доводить итерационный процесс до полной сходимости по условию (25), так как поле давления, входящее с градиентом в правую часть уравнения, приблизительно, и уменьшить время счета на итерациях по совместной коррекции О — р. Компонента Шг, удовлетворяющая расходному соотношению, также изображена на графике и используется для ускорения сходимости в дальнейших итерациях.

Рис. 12. Сходимость компоненты Ш на итерациях по нелинейности на сечении £2 и скорость , удовлетворяющая расходному соотношению.

Число итераций по нелинейности уменьшается по мере установления режима течения, и на каждом из сечений — с увеличением номера итерации по уточнению полей скорости-давления. При этом на итерациях по нелинейности нормы ||Ок,в — Ок'в-1|, к,в — Vк,в-11|, ||Шк,в — Шк,в-11| всегда монотонно стремятся к нулю.

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

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

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

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

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

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

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

[1] Роуч П. Вычислительная гидродинамика. Мир, М., 1980.

[2] Белоцерковский О. М. Численное моделирование в механике сплошных сред. Наука, М., 1984.

[3] Harlow F. H., Welch J. E. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Phys. Fluids, 8, №12, 1965, 2182-2189.

[4] Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. Энергоатомиздат, М., 1984.

[5] Ильин В. П., Туракулов А. А. Об интегро-балансных аппроксимациях трехмерных краевых задач. Новосибирск, 1993 (Препр. ВЦ СО РАН, №986).

[6] Захаров В. П., Соловейчик Ю.Г., Шурина Э. П. Моделирование гидродинамических процессов и теплопереноса во вращающихся криволинейных каналах. В "Численные методы и программное обеспечение", Изд-во ОВМ АН СССР, М., 1990, 42-54.

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

[8] Шурина Э. П., Соловейчик Ю.Г. Решение краевых задач в составных областях: Уч. пособие НЭТИ. Новосибирск, 1986.

[9] Писсанецки С. Технология разреженных матриц. Мир, М., 1988.

[10] Ортега Дж., Рейнболдт В. Введение в параллельные и векторные методы решения линейных систем. Мир, М., 1991.

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