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

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

Том 3, № 6, 1998

МЕТОД ЭКВИРАСПРЕДЕЛЕНИЯ ДЛЯ ПОСТРОЕНИЯ АДАПТИВНЫХ СЕТОК *

Г. С. Хлкимзянов, Н. Ю. ШокинА Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: khak@adm.ict.nsc.ru, shokin@ict.nsc.ru

The equidistribution method for constructing motionless regular grids is considered in the present work. The grids are adapted to the curvilinear boundaries of the domain as well as to the control function which is a priori given in the domain. Using the general methodology based on the equdistribution principle, the nonlinear equations of the elliptic type for the definition of the nodes coordinates inside the domain and on its boundaries are obtained.

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

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

Аналогично работе [1] можно выделить основные принципы адаптации сеток.

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

2. Пользователь контролирует как, когда, где и насколько изменяется сетка.

3. Процесс адаптации сетки не генерирует “плохих” ячеек.

4. Ортогональность линий сетки и другие “хорошие” качества не уменьшаются в процессе адаптации.

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

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №97-01-00819.

© Г. С. Хакимзянов, Н.Ю. Шокина, 1998.

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

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

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

1. Одномерный метод эквираспределения

В работе [2] описывается метод эквираспределения для построения одномерных адаптивных сеток, покрывающих отрезок [0; L]. Суть метода эквираспределения заключается в том, чтобы при конструировании сетки добиться постоянства произведения длины ячейки на значение управляющей функции w в центре ячейки. Пусть Xj — координаты узлов сетки, j = 1, ... , N, при этом Xi = 0, Xn = L. Если длину ячейки [xj; xj+i] обозначить через Sj+i/2, то принцип эквираспределения примет вид

w(Xj+i/2)Sj+1/2 = Ch = const, j = 1, ...,N - 1. (1.1)

Таким образом, длины ячеек будут малы там, где w принимает большие значения и, наоборот, узлы сетки будут иметь разрежения в той части П = (0; L), в которой функция w принимает малые значения. Здесь и всюду далее предполагается, что управляющая функция удовлетворяет ограничениям

1 < w(x) < W, W = const < то, x G Й. (1.2)

Принцип эквираспределения (1.1) является разностным аналогом дифференциального уравнения

dx

w(x) — (q) = C, C = const, q G Q = (0; 1), (1.3)

dq

где

x = x(q), q G Q (1.4)

есть искомое взаимно-однозначное отображение Q на Й. Уравнение (1.3) будем называть принципом эквираспределения в дифференциальной форме. На практике при вычислении

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

^ (ю(х) (?) =0- (1-5)

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

1

h

W(Xi+l/2)

Xj+1 Xj

h

— w(xj-i/2)

Xj — Xj-1 h

0,

2,

N- 1

(1.6)

будем называть ЕБ1-уравнением в разностной форме.

Далее вместо (1.1) и (1.3) будем использовать соотношения, в которые входит якобиан

отображения (1.4). Обозначим его через 3 , 3(д) = — (д). Тогда принцип эквираспреде-

ад

ления (1.3) примет вид

w(x(q))J(q) = C, C = const, q G Q.

Если якобиан J в центре ячейки аппроксимировать центральной разностью

xj+i — xj

j+l/2

h

(1.7)

(1.8)

то из (1.7) будет следовать разностный аналог (1.1) принципа эквираспределения, который можно теперь записать как

w(xj+1/2) Jj+1/2 — -h = Ch

const, j — 1,

N- 1.

(1.9)

В одномерном случае принцип эквираспределения и ЕБ1-уравнение эквивалентны в том смысле, что если отображение (1.4) удовлетворяет принципу эквираспределения (1.7), то функция х = х(^) будет и решением ЕБ1-уравнения (1.5), и наоборот, любое решение ЕБ1-уравнения (1.5) удовлетворяет принципу эквираспределения (1.7). Эквивалентность имеет место и на разностном уровне: если координаты узлов сетки удовлетворяют принципу эквираспределения (1.9), то сеточная функция х^- является решением разностных ЕБ1-уравнений (1.6), а решение этих уравнений дает сетку, удовлетворяющую принципу эквираспределения (1.9).

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

2. Двумерный метод эквираспределения

В двумерном случае аналогом длины одномерной ячейки является площадь ячейки двумерной сетки с узлами х^-, покрывающей область П с границей Г. Здесь к = (ж1, ж2),

І = (ІьІ2) — мультииндекс, І = 1, ..., Да. Площадь четырехугольной ячейки с вершинами х^1а-2, х^1+1а2, х^1+1)^2+1, х^1)^2+1 будем обозначать через 5^+1/2. Тогда принцип эквираспределения, аналогичный (1.1), может быть записан в виде

w(Xj+i/2)Sj+i/2 — Ch, ja — 1, . . . , Na — 1,

(2.1)

где w(xj+1/2) — значение заданной управляющей функции в центре ячейки, координаты которой определяются как среднее арифметическое соответствующих координат четырех ее вершин.

Совокупность узлов сетки обозначим через Гh. Будем считать, что множество Г является образом равномерной прямоугольной сетки Qh, покрывающей квадрат Q = [0; 1] х [0; 1], при взаимно-однозначном отображении

xa = xa(q1, q2), a = 1, 2 (2.2)

плоскости q1Oq2 на плоскость декартовых координат x1Ox2. Отображение (2.2) является неизвестным. Выведем дифференциальное уравнение для определения этого отображения. Возьмем тождество, которому удовлетворяет произвольное гладкое отображение (2.2):

д i g22 dxa g12 dxa \ + д i gu dxa + gn 5xa \ =0

dq1 V J dq1 J dqV + dq2 V J dq1 + J dq2 j = . ( . )

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

через дав.

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

д12 = 0. (2.4)

И, во-вторых, будем считать, что отображение (2.2) удовлетворяет принципу эквираспределения в дифференциальной форме, аналогичному принципу (1.7) в одномерном случае:

w(x(q))J(q) = C, C = const, q =(q1,q2) G Q, (2.5)

где J — якобиан отображения (2.2). Использование этих условий в тождестве (2.3) приводит к двумерным уравнениям метода эквираспределения или к ЕБ2-уравнениям в дифференциальной форме:

д / dxa \ д / dxa \

д? (wg22 + з? (wg11 = 0- a = 1-2 <2-6)

Итак, в двумерном случае ЕБ2-уравнения следуют из принципа эквираспределения

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

Лемма 2.1. Если отображение (2.2), задаваемое решениями уравнений (2.6), является невырожденным и ортогональным, то для него выполняется принцип эквираспределения (2.5).

Доказательство. Использование равенства g12 = 0 из условия леммы и тождества

д x д _ \ « I /Э i д x

geeдд3-в - g12 ддв = (-1) J^qT, а,в = 1, 2,

в котором суммирования по повторяющимся индексам а и в не производится, приводит к соотношениям

джа . .„+-, тдж3—а джа . тдж3—а

922&? = (-1) '/_&Г’ 911 а? =(-1) }~е^Т’ а = 1,

Подставляя их в уравнения (2.6), получим однородную систему из двух уравнений

ди> 3 джа ди> 3 джа

0, а = 1, 2

дq1 дq2 дq2 дq1

относительно производных дwJ/дqа, определитель которой совпадает с якобианом J. В

силу предположения J = 0 эта система имеет лишь тривиальное решение дwJ/дqa = 0,

a = 1, 2. Отсюда следует, что wJ = const.

Покажем, как следует аппроксимировать уравнение (2.6), чтобы указанное в лемме свойство выполнялось и на разностном уровне. Покроем Q равномерной прямоугольной сеткой с количеством узлов Na и шагом —a в направлении оси Oqa, —a = 1/(Na — 1), a = 1, 2. Пусть qj = (q^, q22), j = (j1,j2), ja = 1, ..., N, qja = (j« — 1)—«. Введем два базисных оператора разностных производных:

^(q1 + h1/2,q2) — ^(q1 — —1/2,q2) (27)

(Dqi ^)(q) = ----------------—------------------, (2.7)

^(q1,q2 + —2/2) — ^(q1,q2 — —2/2) (20)

(Dq2 ^)(q) j . (2.0)

—2

При использовании этих производных для записи разностных операторов необходимо учитывать область определения сеточных функций, так как в некоторых случаях сеточные функции будут относиться не к узлам сетки с целыми индексами, а, например, к центрам ячеек (узлы с двумя полуцелыми индексами) или к центрам сторон ячеек (узлы с одним целым и одним полуцелым индексами). Таким образом, результат действия операторов зависит от того, в каких узлах определена сеточная функция <^. Этот результат будет также зависеть и от того, в какой точке требуется вычислить производную. Если функция ^ определена в целых узлах qj = (q^, q22), то полагаем

(Dq1 ^)ji + 1/2,j2

ji j2

^jl + 1,j2 — ^jl,j2

—1

( П.1 ro) _• 1 /2 _• +

(Dq1 ^)ji,j2

(Dq1 ^)jl -1/2,j2 + (Dq1 ^)jl + 1/2,j2

,j2 =

(Dq1 ^)j1,j2 + 1/2

2

(Dq1 ^)j1 ,j2 + (Dq1 ^)j1,j2 + 1

2

/п \ _ Фд1 ^)л+1/2,^2 + Фд1 +1/2,^2+1

(^д1 ^)л + 1/2,^2 + 1/2 = ----------------------------------2-'

Формула (2.7) для сеточных функций, определенных в центрах ячеек Ч^+1/2 =

(91+1/2’ С+1/2 = С + ^«/2, принимает следующий вид:

^^1+1/2,^2+1/2 - <£л —1/2 ,л + 1/2

q1 г )j1,j2+1/2

—1

/п ч Фд1 ¿2-1/2 + Фд1 ^)іі,І2+1/2

^)Л;І2 =--------------------------------------------^-’

/п ч _ (^д1 ^)Іі,І2+1/2 + Фд1 ^)І1+1,72+1/2

(^д1 ^)л+1/2,і2+1/2 =---------------------------------------------~,

(Дд1 ^)^1 ,^2 + (^д1 ^)І1 + 1,І2

'д1 ^)І1+1/2,І2 = 2

(Дд1 ^)^1 + 1/2,_

Если сеточная функция ^ определена в центрах Ч^1+1/2,^2 = (9)1+1/2’ ?|2) горизонтальных сторон ячеек, то под разностной производной (2.7) понимаются следующие выраже-

ния:

< ТЛ ,Л _ ^71+1/2,72 ^І1-1/2,І2

(Дд1 ^)і = ^

(Дд1 ^л+1/2 (Дд1 ^л ¿2 + 1/2 (Дд1 ^)І1 + 1/2,І2 + 1/2 =

= (Дд1 ^л ,л + (Дд1 ^Л + 1,72 І2 = 2

= (Дд1 У)^1 ,^2 + (Дд1 У)Л,72+1

= 2 ’

(Дд1 ^)І1+1/2,і2 + (Дд1 ^л+1/2,72+1

2

Аналогично, с использованием при необходимости осредненных значений <^, вычисляется разностная производная 0,1 ^ для сеточных функций, определенных в центрах ^2+1/2 = (?! ,9|2+1/2) вертикальных сторон ячеек, а также производная 0,2<^. Разностные уравнения для координат узлов получаются интегро-интерполяционным методом. Для этого уравнение (2.6) заменяется интегральным соотношением

джа джа

^922^д2 - и>0п—^д1 = 0’ (2.9)

с

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

(911)7+1/2 = (-0,1 ж1 Цд ж1 + Цд ж2 0,1 ж2)^+1/2, (2.10)

(912)7+1/2 = (0,1 ж10,2 ж1 + 0,1 ж2 0,2 ж2) 7+1 /2, (2.11)

(922)7+1/2 = (0,2 ж10,2 ж1 + 0,2 ж2 0,2 ж2) 7+1 /2, (2.12)

3+1/2 = (0,1 ж10,2ж2 - 0,2ж10,1 ж2) ,+1/2 . (2.13)

В результате будем иметь следующие разностные ЕБ2-уравнения:

(0,1 (^922^^1 Xа) + 0,2 (идп0,2 жа))^. = 0, (2.14)

при этом следует учитывать, что сеточные комплексы Дд3-вжа, а, в = 1, 2 определены

также в центрах ячеек.

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

Отметим, что если якобиан отображения вычисляется по формуле (2.13), то

Sj+1/2 = Jj+1/2h1h2; (2-15)

поэтому принцип эквираспределения в разностной форме (2.1) можно переписать в виде, аналогичном (1.9):

Ch

w(Xj+i/2)J?-+l/2 = 7—^ = Ch = const, ja = 1, . . . , N - 1. (2.16)

hlh2

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

(g12) j+l/2 = 0. (2.17)

Сетку с выпуклыми ячейками будем называть адаптивной, если для каждой ячейки

выполняется принцип эквираспределения в разностной форме (2.16).

Выше было отмечено, что для ортогональных двумерных отображений (2.2) принцип эквираспределения в дифференциальной форме (2.5) и ЕБ2-уравнение (2.6) эквивалентны. На разностном уровне удается доказать, что сеточные функции, определяющие криволинейную сетку, удовлетворяющую принципу эквираспределения в разностной форме (2.16), являются решениями разностных ЕБ2-уравнений (2.14).

Лемма 2.2. Координаты узлов любой квазиортогональной адаптивной сетки с выпуклыми ячейками удовлетворяют уравнению (2.14).

Доказательство. Используя операторы (2.7), (2.8), легко убедиться, что для сеточных функций справедливы формулы

(gee xa - g12Dqe x“).+1/2 = (-1)a+e-1 (JD,e Х3-а)^ + 1/2 , (2.18)

а, в = 1, 2, суммирование по индексам а и в не производится. Следовательно, с учетом условия квазиортогональности (2.17), левую часть уравнения (2.14) можно переписать в следующем виде:

1 {((Дд )Dq2 xa)j1j2 + 1/2 + ((Dq! )Dq2 ^ *^-1/2 ~

- ((Dq2 )Dq1 xa)j1+1/2;i2 - ((Dq2 )Дд ) j1 - 1/2,j2 } . (2.19)

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

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

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

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

Поэтому будем использовать следующий метод определения узлов на границе, в котором вместе с двумерным методом эквираспределения внутри области используется одномерный метод эквираспределения на границе. Для определенности возьмем нижнюю сторону д2 = 0 квадрата Q и соответствующую ей при отображении (2.2) часть границы дП. Пусть эта часть границы задана в параметрическом виде

и имеет длину Ь. Если через в = в(д) обозначить длину дуги границы, то для в(д) получим следующую формулу:

Чтобы узлы на границе имели сгущение в соответствии с управляющей функцией и>, потребуем выполнения на рассматриваемом участке границы аналога уравнения (1.5) для длины дуги:

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

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

3. Трехмерный метод эквираспределения

Теперь рассмотрим принцип и уравнения метода эквираспределения для построения сетки внутри трехмерной односвязной области П с неподвижной границей Г. Пусть

(2.20)

(2.21)

(2.22)

где ад(в) = ^(«(д1)) = и>(/ 1(д(д1)), /2(д(д1))), а параметры д1 и д связаны соотношением (2.21), в левой части которого следует использовать величину «(д1). Уравнение (2.22) дополняется краевыми условиями

з(0) = 0, з(1) = Ь,

Xа = Xа(д1, д2, д3), а = 1, 2, 3

(3.1)

есть взаимно-однозначное, невырожденное, достаточно гладкое отображение единичного куба Q на область О (рис. 1). Для такого произвольного отображения выполняется тождество, аналогичное (2.3) и вытекающее из записи уравнения Лапласа в новых коор-

динатах qe:

9 Jg7/3=°, а = 1 2, 3- (3-2)

dqY \ dqe

В отличие от (2.3) здесь для сокращения записи мы использовали контравариантные компоненты #7в метрического тензора, а также тензорное правило суммирования по повторяющимся верхнему и нижнему индексам 7 и в. Это правило будет, если не оговорено противное, использоваться и в дальнейшем. Для ортогонального отображения (3.1) вы-

Рис. 1. Физическая область О и вычислительная область Q. полняются условия

#12 = 0, #13 = 0, #23 = 0, (3-3)

и тождества (3.2) принимают вид уравнений, не содержащих смешанных производных от функций ха:

д (#22#33 dxa \ + _д_ (gllg33 dxa \ + _д_ (#11 #22 dxa \ =0 (34)

dq1 V J dqV + dq2 V J dq2 У + dq3\ J dq3 / ' ( - )

Если потребовать, чтобы отображение (3.1) подчинялось принципу эквираспределения в дифференциальной форме

w(x(q))J(q) = C, C = const, q = (q1 ,q2,q3) G Q, (3.5)

где J — якобиан отображения (3.1), а x = (x1, x2, x3), то функции xa = xa(q 1,q2, q3)

будут решениями трехмерных уравнений метода эквираспределения в дифференциальной

форме (ЕБ3-уравнений):

d / dxa \ d f dxa \ d f dxa \

L22 Sn #33, j + _q. L11 #22W_ = 0- (3.6)

Также как и в двумерном случае, справедливо утверждение, обратное к сформулированному выше.

Лемма 3.1. Если отображение (3.1), задаваемое 'решениями уравнений (3.6), является невырожденным и ортогональным, то оно будет и адаптивным в смысле выполнения равенства (3.5).

Доказательство. Из условия леммы следует, что дав = 0 при а = в, поэтому функции ^ = ха удовлетворяют уравнению

д дха

д дхт /и) = 0.

д«т д«в

Дифференцируя его и учитывая тождество (3.2), получим систему уравнений относительно величин д (Зи^/д«7:

дха д3и .

£ = 0, а = 1, 2,3. (3.7)

д«в д^ ^ •

Определитель этой системы равен 33. В силу невырожденности отображения отсюда будет следовать, что однородная система (3.7) имеет лишь тривиальное решение д(3и)/дд7 = 0, 7 = 1, 2, 3. Таким образом, выполняется равенство (3.5).

Покроем единичный куб 5 прямоугольной равномерной сеткой с числом узлов Да в направлении оси 0«а. Узлы сетки будем обозначать через qj = (д]1, «|2, «33), 3 — мультииндекс, 3 = (31, 32, 3з), За = 1,..., Да, 3 = (За - 1)^а, = 1/(Да - 1). Совокупность

всех узлов qj обозначим <5^, внутренних — 5^ и граничных — 7^, т. е. <5^ и 7^.

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

Яп + 1/2,.72,.73 = (д11 + 1/2, 0-32 , 3 ) , qjl ,¿2 +1/2,73 (3, «32+1/2 ^ «33),

qл ,32,3з+1/2 (Qjl, Qj2, д3з+1/2),

где да«+1/2 = С + ^аА

Для узлов, совпадающих с центрами граней, полагаем

^?1+1/2,32+ 1/2,3з = /2+1/2, д3з), qjl,j2+1/2,33+1/2 = (4, q|з+l/2),

^1 + 1/2,32,33+1/2 = (031+1/2, 3 , 03з+1/2)-

Наконец, центр ячейки имеет координаты

Ч?+1/2 = qjl+ 1/2,32+1/2,3з+1/2 = (о11+1/2, о22+1/2, °3з+1/2)-

Введем базисные операторы разностных производных:

(Д,,.^) = ^ + ^1/2-92.93)~ ^ - У2.«2, (3.8)

(В,2.^) = + У2^ - ^ОД«2 - ), (3.9)

^2

= ^(q1■q2■q3 + ^3/2) - ^(q1■q2■q3 - ^з/2) (310)

(Д93 т . (3.10)

Д3

При вычислении производных (3.8)-(3.10) необходимо, так же как и при вычислении разностных производных (2.7), (2.8), учитывать область определения сеточной функции ^ и при необходимости использовать усреднение. Например, при вычислении разностных производных от функций жа, определенных в целых узлах q7•, полагаем

Ап Ж71 + 1,72 + 1/2,73 + 1/2 Ж7Ъ72 + 1/2,73 + 1/2 /О 1 1 А

№ Ж 77 + 1/2 = ---------------------^-------------------’ С3-11;

где ж“ 72+1/2 73+1/2 — координаты центра грани с номером ^'1.

Ячейка криволинейной сетки (рис. 2), покрывающей физическую область П, является шестигранником с неплоскими, вообще говоря, гранями. Хотя грань является неплоским четырехугольником, тем не менее отрезки, соединяющие середины противоположных сторон грани, лежат в одной плоскости и пересекаются в некоторой точке, которую и будем называть центром грани. Например, для грани ^ ее центром будет точка х71,72+1/2,7-3+1/2, каждая координата которой равна среднему арифметическому соответствующих координат четырех вершин этой грани.

Рис. 2. Ячейка криволинейной трехмерной сетки, средние линии ячейки и граней, центры граней и центр ячейки.

Указанные выше отрезки будем называть средними линиями граней. Средними линиями ячейки будем называть отрезки, соединяющие центры противоположных граней. Легко показать, что средние линии выпуклой ячейки пересекаются в одной точке. Эту точку будем называть центром ячейки и обозначать х^+1/2.

Объем ячейки будем вычислять по формуле

£¿+1/2 = (хп + 1,72 + 1/2,73 + 1/2 - хл ,72 + 1/2,73 + 1/2) ‘ [ (Х71 + 1/2,72 + 1,73 + 1/2 -

-ХЛ + 1/2,72,73 + 1/2) Х (Х71 + 1/2,72 + 1/2,73 + 1 - Х71 + 1/2,72 + 1/2,7з)]

Если при вычислении якобиана J отображения (3.1) использовать базисные операторы (3.8)-(3.10), то получим равенство, аналогичное (2.15):

Sj+1/2 = Jj+1/2h1h2 h3, (3.13)

поэтому аналогом (3.5) на разностном уровне будет требование постоянства произведения объемов ячеек на управляющую функцию w:

w(Xj+i/2)Sj+i/2 = Ch, ja = 1, ...,Na - 1. (3.14)

Это равенство будем называть трехмерным принципом эквираспределения в разностной

форме. Очевидно, что в силу (3.13) этот принцип можно записать и в ином виде:

Ch

w(Xj+l/2)Jj+1/2 = , , , = Ch = const, ja = 1, ...,Na - 1. (3.15)

h1h2h3

Уравнение (3.6) аппроксимируется со вторым порядком следующим конечно-разностным уравнением:

(Л1Ха + Л2^а + Лзха)j = 0, q7 Є Qh, (3.16)

(A1Xa)j = Dql (g22933WDqlxa) . , (Л2Жа)7 = Dq2 (g11033WDq2Xа) . ,

где

(Л1.ха)7 = ^„1 (о22о33?/;^„1 „у- у-~ 7

(АэЖа)7 = Д„3 (^11022^Д„3ха)7 ,

при этом предполагается, что сеточные функции определены в целых узлах q7•, а

компоненты метрического тензора ¿ад и управляющая функция — в центрах ячеек

q7'+l/2• Уравнения (3.16) будем называть уравнениями метода эквираспределения или ЕБ3-уравнениями в разностной форме.

Для решения системы (3.16) использовался метод стабилизирующей поправки:

дП+1/3 _

д^ = А1^га+1/3 + Л2^га + Аэ^га, (3.17)

т

м«+2/^ дп+1/3

ДД------------= Л2 (дп+2/3 - дп), (3.18)

т

/д-^1 _ ,дп+2/3

дд----------= Лз (дп+1 - дп), (3.19)

т

здесь д = ха, п — номер итерации, т — итерационный параметр, подбираемый экспериментально.

Каждый из шагов (3.17) - (3.19) реализовывался скалярными прогонками, при этом предполагалось, что координаты узлов на границе Г = дП уже известны.

4. Построение сетки на пространственных поверхностях

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

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

Методику вывода поясним на примере получения уравнений для части Г bot граничной поверхности Г, являющейся образом нижней грани q3 = 0 куба Q. Предполагается, что эта часть границы задана в следующем параметрическом виде

x

= /°V,P2), 0 < < 1, а = 1, 2, 3, в = 1, 2.

(4.1)

Относительно преобразования (3.1) делаются следующие предположения:

— координатные линии этого отображения ортогональны на Г^;

— кривизна координатных линий третьего семейства в точках их пересечения с равна нулю;

— координатные поверхности q3 = const “параллельны” поверхности Tbot в некоторой ее окрестности.

Отметим, что второе предположение использовалось ранее в работе [9] при получении уравнений для поверхностной расчетной сетки из уравнений [8] для внутренних узлов.

Будем использовать следующую математическую форму этих предположений:

012 = 013 = g23 = 0, q = (q\ q2, 0)

d Эз 1 Ô033,

dq3 2033 dq3

033 = const, q

^ q = (q ,q , 0)

(q1, q2,0).

(4.2)

(4.3)

(4.4)

Через Эа обозначен вектор, касательный к координатной линии ^а, то есть Эа = дr/дqa, где г — радиус-вектор точки.

Лемма 4.1. Пусть отображение (3.1) удовлетворяет условиям (4-2) - (4-4)- Тогда

д0

3а

dq3

Доказательство. Из равенства

0, а = 1, 2.

(4.5)

q3 =0

дЭ3

dq

— Г7 Э

3 = Г 33Эт

при учете (4.3) будет следовать

Г133 = 0, Г233 = 0.

Ввиду ортогональности (4.2) системы координат получим

Га____ааГ = 1 1 [-g3a + dg3a -g33 ^_______2

Г33 =g Г33а = gaa'На? + a?-sq^l■ а = 1-2

Использование здесь условия (4.4) приводит к равенству (4.5).

Лемма 4.2. Пусть выполнены условия (4-2) - (4-4) и уравнение (3.6) выполняется вплоть до границы Tbot. Тогда

wJ(q\q2,0) = const, 0 < qa < 1, a = 1, 2. (4.7)

Доказательство. В силу (4.2) и полученного равенства (4.5) уравнение (3.6) можно записать на Tbot в следующем виде:

д / -та \

т—в ( wJ2geY —— ) =0, a =1, 2, 3. (4.8)

dqe V OqYj

Для любого в =1, 2, 3 имеем равенство

я7—та = —та = ,m)m -¿_^= —qe (4 9)

g dqY -T™ -T™ -qY -T™ ™ -X« , ( . )

поэтому уравнение (4.8) можно переписать в виде

-qe — — / —qe \

Jt—ттв(wJ)+ wJ^-e ( J^—) = 0, a = 1, 2, 3. (4.10)

-xa dqe^ ’ dqe \ dxay v ’

Пусть u = (и1,«2,«3), причем uY = £7a, тогда ve = uYdqe/dxY = dqe/дта. Отсюда получается, что

. 1 д / т dqe

dlV U = J -

3 дqв \ дж°

В силу того, что и = 0, второе слагаемое в (4.10) обращается в нуль, и тогда равенство (4.7) будет следовать из того, что величины д(и>3)/дqв являются решением однородной системы (4.10), определитель которой не равен нулю.

Доказанное равенство (4.8) означает, что отображение (3.1) удовлетворяет принципу эквираспределения в дифференциальной форме на граничной поверхности Г^-

При выводе уравнений для вычисления координат узлов на поверхности будем использовать тот факт, что величины д (и>3)/дqa = 0 можно считать решением однородной системы

дра ди> 3 дра ди>3

-q2 -q1 -q1 -q2

0, a = 1,2, (4.11)

определитель которой 3 = не равен нулю. Пусть дар — ковариантные ком-

поненты метрического тензора взаимно-однозначного невырожденного преобразования

ра = р“^1 ^2), а = 1, 2 (4.12)

с якобианом 3 > 0, т. е.

др1 др1 др2 др2

д“в = драдрв + драдрв' а'в =1'2 (4Л3)

Тогда для отображения (4.12) справедливы формулы

где

Jg/Î 1 dq1 + ^ dq2 = (-1Г+в dq3-e , а’в = 1 2, (414)

J11 = gJ2, g22 = ^, J12 = - ^. (4.15)

J2 J2 J2

Переписав уравнения (4.11) в дивергентной форме

d ( Tdpa\ д ( Tdpa

wJ^^- — wJ

dq1 \ dq2/ dq2 \ dq1

и использовав в них тождества (4.14) при дополнительном предположении об ортогональности системы координат (4.12), получим следующие уравнения для функций р“^1,?2), а = 1, 2:

д ( ----------д22 др“\ . д ( ------------дра

д^Г ( {\/022033 — д^Г ) + д^2 ( {^022033 — д^2

|î (w^g22g33JJ2§pr) + dq2 (wV0^ydpi) = °. (416)

При выводе данного уравнения использовалось также условие (4.4).

При численном решении (4.16) вместо компонент gaa используются ковариантные компоненты 'д'ав метрического тензора поверхности Г bot. Последние связаны с величинами дад соотношениями :

dp7 ôp£

g«e = дтвdqadqe, a,e,Y,e = 1,2, (4.17)

при этом

dr7 dr7 dpa dpe

g«e = ^7,7, a,e, = 1, 2 (4.18)

(здесь по индексу 7 = 1, 2, 3 производится суммирование).

Конечно-разностное уравнение для вычисления координат узлов х^^д было получено аппроксимацией уравнения (4.16) на 9-точечном шаблоне с помощью интегро-интерполя-ционного метода. По виду это разностное уравнение совпадает с (2.14):

(А?1 (^022^1 ра) + ^2 (^011^2 ра))7ь72д = 0, (4Л9)

где сеточная функция гй = ^^/д22д33/— относится к центрам ячеек. Решаются эти системы нелинейных разностных уравнений итерационным методом продольно-поперечных прогонок.

Отметим, что если поверхность ГЬо* является плоской, например, лежит на плоскости х3 = х3, и в качестве параметров ра взяты декартовы координаты (а = 1, 2), то 0«в = ^«в, 0«в = 0«в, — = 10 и уравнения (4.16) переходят в уравнения (2.14) ЕБ2-метода построения сеток в плоских областях.

Метод построения адаптивных сеток на основе уравнений (4.16) будем называть ЕББ-методом для поверхностей.

°

5. Построение сетки на пространственных кривых

Получим уравнения метода эквираспределения для вычисления координат узлов на пространственной кривой, являющейся одним из ребер криволинейного шестигранника П. Пусть для определенности это ребро (далее обозначается L) принадлежит двум поверхностям, являющимся образами при отображении (3.1) граней q3 = 0 и q2 = 0 единичного куба Q.

Исходя из общего принципа метода эквираспределения, будем строить сетку на рассматриваемом ребре так, чтобы длины отрезков между двумя соседними узлами сетки были обратно пропорциональны значениям управляющей функции в центрах этих отрезков. Причем для лучшего согласования положения узлов внутри области и на кривой L желательно использовать одну и ту же управляющую функцию w. Однако в таком случае w будет отвечать за степень концентрации узлов по разным признакам: внутри области значения w регулируют объем ячеек, а на L — длину одного из ребер шестигранной ячейки. В общем случае использование одной и той же управляющей функции может приводить, например, к тому нежелательному эффекту, что в приграничном слое, примыкающем к L, объемы ячеек при движении вдоль L будут уменьшаться, а длины ребер этих ячеек, лежащих на L, будут увеличиваться. Чтобы избежать таких последствий, приходится накладывать на отображение (3.1) дополнительные ограничения.

Итак, пусть функции ха, определяющие отображение (3.1), являются решениями ED3-уравнений (3.6). Сформулируем дополнительные предположения о поведении функций в окрестности кривой L, при выполнении которых эти функции удовлетворяют принципу эквираспределения на L:

л/дйw = const, q = (q1,0, 0). (5.1)

Так как L является линией пересечения двух боковых граней области П, на которых уравнения EDS-метода получены при условиях типа (4.2)-(4.4), эти же условия можно использовать и при выводе уравнений для сетки на L. В действительности же для их

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

используются следующие ограничения на преобразование (3.1):

д„в = 0, а = в, x е L, (5.2)

дд0

W = Г33Э3- fp2 = Г22Э2• X е L (5-3)

= 0, а = 2,3, в = 1,2,3, в = а, x е L, (5.4)

dqe

означающие, что

— координатные линии этого отображения ортогональны на L;

— кривизна координатных линий второго и третьего семейства в точках их пересечения с L равна нулю;

— координатные поверхности q2 = const и q3 = const “параллельны” кривой L в некоторой ее окрестности.

Лемма 5.1. Пусть выполнены условия (5.2) - (5.4) и уравнение (3.6) выполняется в П вплоть до граничного ребра L. Тогда на L будет справедливо равенство (5.1).

Доказательство. В силу условия ортогональности (5.2) и предположения (5.3) имеем

где индексы а и в принимают значения, указанные в условии (5.4). Учет этого условия в

(5.5) приводит к равенствам

д0ав = 0, а = 2, 3, в = 1, 2, 3, в = а, (5.6)

дд« ^ '

которые выполняются в каждой точке кривой С.

В силу условия ортогональности (5.2) полученные равенства приводят к следующим соотношениям на кривой С:

д 0а в

0ав = 0, = 0, а = 2,3, в = 1,2,3, в = а, (5.7)

поэтому уравнение (3.6) можно переписать на С в виде (4.8) или (4.10), откуда и получим соотношение (5.1), учитывая, что на С компоненты 022 и 033 постоянны.

Так как С является пересечением двух поверхностей д2 = 0 и д3 = 0, то можно считать, что их параметрические уравнения (4.1) совпадают на С и имеют в точках С вид

= Г(р), 0 < р < 1, а = 1, 2, 3. (5.8)

Это параметрическое уравнение кривой С используется в уравнениях метода эквираспределения (в ЕБСЗ-уравнениях)

d Л2 + (d/2V + (/ 3У dp , 0

дд1 1 V V др / \ др ) \ др ) дд1

которые получаются дифференцированием равенства (5.1) по переменной д1. Для плоской кривой с естественной параметризацией ЕБСЗ-уравнение переходит в уравнение (2.22). Если С является прямолинейным отрезком, то из (5.9) следует ЕБ1-уравнение (1.5).

Для аппроксимации (5.9) будем использовать 3-точечное разностное уравнение

А?1 (ЙШ?1 р)Л = 0, = 2, ...,N1 - 1, (5.10)

где ______________________________

w = w(f 1(p), / 2(p), / 3(p))^/(D„/1)2 + (Dp/2 )2 + (Dp/3)2,

при этом предполагается, что сеточная функция w определена в полуцелых узлах pj1+1/2 неравномерной сетки р^х, а функция p(q) — в целых узлах равномерной сетки qjx.

Нелинейное уравнение (5.10) решается итерационным методом. В качестве начального приближения для функции р = p(q1) бралась линейная функция. При вычислении величины р на (v + 1)-й итерации использовались значения wj1+1/2 с v-й итерации. Если итерационный процесс сходится, то для предельных значений координат узлов на L будет выполняться равенство

wji+1/2,1,1 |Xji+i,i,i - Xj-1,1,11 = const, ji = 1,...,Ni - 1. (5.11)

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

Рис. 3. Сетки после 100 итераций с управляющими функциями: а — w = 1; б — w = 1 + 10(z + 1); в — w = 1; г — w = 1 + 10(z + 1).

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

На рис. 3 приведены примеры трехмерных сеток, построенных с помощью метода эквираспределения .

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

[1] Patel M. K., Periclequs K. A., Baldwin S. The development of a structured mesh

grid adaption technique for resolving shock discontinuities in upwind Navier — Stokes codes. Int. J. Numer. Methods in Fluids, 20, 1995, 1179-1197.

[2] BOOR C. Good approximation by splines with variable knots. II. Lecture Notes in

Mathematics, 363, 1974, 12-20.

[3] DWYER H. A., Kee R. J., SANDERS B. R. An adaptive grid method for problems in fluid

mechanics and heat transfer. AIAA J., 18, 1980, 1205-1212.

[4] DWYER H.A. Grid adaption for problems in fluid dynamics. Ibid., 22, No. 12, 1984, 1705-1712.

[5] SHYY W. Computation of complex fluid flows using an adaptive grid method. Int. J. Numer. Methods in Fluids, 8, 1988, 475-489.

[6] CHANG P.Y., Shyy W. Adaptive frid computation of three-dimensional natural convection in horizontal high-pressure mercury lamps. Ibid., 12, 1991, 143-160.

[7] HuANG W., SLOAN D. M. A simple adaptive grid method on two dimesions. SIAM J. Sci. Comput., 15, No. 4, 1994, 776-797.

[8] THOMPSON J.F., Warsi Z.U.A., Mastin C.W. Numerical grid generation, foundations and applications. Amsterdam: Nort-Holland, 1985, 483.

[9] TAKAGI T., Miki K., Chen B.C. J., SHA W. T., Numerical Generation of Boundary-Fitted Curvilinear Coordinate Systems for Arbitrarily Curved Surfaces. J. of Comp. Physics, 58, 1985, 67-79.

Поступила в редакцию 19 августа 1998 г.