CC BY
62
Вычислительные технологии
Том 5, № 1, 2000
ОБ ИНТЕРАКТИВНОМ КОМПЛЕКСЕ ПРОГРАММ ПОСТРОЕНИЯ ДВУМЕРНЫХ СТРУКТУРНЫХ
СЕТОК*
В. Д. ЛисЕйкин, Ю. И. Молородов, Г. С. Хлкимзянов Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: liseikin@net.ict.nsc.ru, yumo@net.ict.nsc.ru e-mail: khak@adm.ict.nsc.ru
The paper is devoted to the description of interactive software for the generation of curvilinear structure grids covering two-dimensional one-connected domains of a complicated form. The algorithms for the generation of boundary grids, the initial grids for iterative processes and final grids of elliptic methods are given. The criteria are described which allow one to estimate the quality of grids generated by different methods and automatically to select the best of them. Some examples of grids constructed making use of the described software are given.
В настоящее время при численном решении задач математической физики широко используются алгоритмы расчета на неравномерных структурных или неструктурных сетках. Довольно широкий круг алгоритмов основан на применении конечно-разностных методов расчета на структурных сетках, адаптирующихся как к границам области, так и к особенностям отыскиваемого решения. Такие сетки, как правило, связаны с зонами больших градиентов решения, поэтому их сгущение в нужных подобластях может приводить к повышению точности численного решения по сравнению с результатами, получаемыми на неадаптивных сетках того же размера. Разнообразные методы построения адаптивных сеток приведены в монографиях [1, 4, 14, 22, 25], обзорах [5-7, 23, 26], электронных публикациях [28].
В комплексе программ, о котором пойдет речь в настоящей статье, реализованы методы построения двумерных адаптивных сеток, основанные на итерационных процессах решения эллиптических уравнений для координат узлов. В частности, используются уравнения, полученные из принципа эквираспределения [10, 12, 16] и не относящиеся к типу "обращенных" эллиптических уравнений, которые составляют ядро широко известной программы EAGLE [27] и ее различных модификаций (см., например, [21]).
Для того чтобы начать итерационный процесс получения численного решения рассматриваемых эллиптических уравнений, необходимо задать "качественное" начальное приближение. В статье уделяется много внимания построению такого приближения, так как
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, гранты №97-01-00819, 96-01-01705.
( В. Д. Лисейкин, Ю.И. Молородов, Г. С. Хакимзянов, 2000.
скорость сходимости итераций, да и сама сходимость процесса сильно зависят от начального приближения. Для построения начальной сетки (начального приближения) в комплексе программ реализованы алгоритмы двух алгебраических методов: метода построения сетки в односвязных областях звездного типа с границей, составленной из конечного числа отрезков прямых (далее — метод STAR) и метода трансфинитной интерполяции (далее — метод TFI). При использовании классического метода TFI для областей сложной формы не всегда удается получить начальную сетку со сгущениями в требуемых подобластях и без самопересекающихся ячеек. На этот случай в комплексе предусмотрена возможность построения начальной сетки модифицированным методом TFI с введением в интерактивном режиме дополнительных внутренних сечений [18], применение которых позволяет повысить эффективность регулирования сгущений сетки и избавиться от самопересекающихся ячеек.
Обращает на себя внимание тот факт, что в большей части работ, посвященных построению расчетных сеток, даны только описания алгоритмов, сопровождаемые изображением полученных сеток. Поэтому о качестве сеток в этих работах можно судить лишь визуально. Между тем "красивая" сетка может оказаться подходящей для одних задач и совершенно непригодной для решения в той же самой области других задач. О качестве сетки можно судить лишь в строгой связи с конкретной задачей, для решения которой она используется. С другой стороны, накопленный опыт использования криволинейных сеток в расчетах задач математической физики говорит о том, что существуют некоторые объективные характеристики сеток (ортогональность, локальная равномерность, невытянутость ячеек и т. д.), по которым можно априорно оценить качество сетки и ее пригодность для расчетов. Таким образом, возникает необходимость в автоматическом определении таких характеристик. Некоторые подходы для анализа качества сеток и сравнения программ их построения приведены в [9]. В рассматриваемом комплексе программ имеется блок, в котором с помощью критериев качества производится оценка качества как начальных, так и окончательных сеток (сеток, полученных в результате итерационного решения эллиптических уравнений).
В работе приводятся также описание программы, некоторые особенности ее использования, демонстрационные примеры и результаты решения на адаптивной сетке смешанной задачи для уравнения Лапласа.
1. Расстановка граничных узлов
Предполагается, что область П, лежащая в плоскости декартовых координат x1Ox2, топологически эквивалентна единичному квадрату Q из плоскости q1Oq2. Пусть взаимно однозначное соответствие между точками этих областей задается отображением
где х =(ж1,ж2), q = (д1,д2), х^) = (х1^1, д2), х2^1, д2)). Тогда границу Г = дО области О можно представить в виде объединения четырех ее частей: Г = Г^ У Г^ и и Гюр> являющихся образами левой, нижней, правой и верхней сторон квадрата Q соответственно.
Отображение (1.1) заранее неизвестно, поэтому разбиение границы Г на четыре последовательные части должно указываться пользователем комплекса. Для описания каждой из частей он может использовать либо параметрическое задание
x = x(q),
(1.1)
x
= fe(q), 0 < q < 1, в = 1, 2,
(1.2)
либо табличное в виде совокупности точек, называемых опорными. Чтобы при табличном задании криволинейная граница имела гладкую форму, необходимо задавать достаточно много опорных точек. Алгоритм построения сетки на границе существенно использует предположение о том, что граница состоит из прямолинейных отрезков, соединяющих соседние опорные точки. Поэтому даже при аналитическом задании некоторой части границы программа выстраивает достаточно плотное множество опорных точек и расставляет граничные узлы по полученной ломаной.
В комплексе программ для расстановки узлов на каждой из четырех частей границы используется одномерный метод эквираспределения [10, 12]. Вначале узлы на границе расставляются равномерно по длине дуги. Затем в ходе итерационного процесса их положение изменяется так, чтобы в итоге удовлетворялся бы принцип эквираспределения: произведение длины дуги между соседними узлами на значение управляющей функции т в середине этой дуги должно быть величиной, постоянной для всех ячеек одномерной сетки, покрывающей рассматриваемую часть границы. Таким образом, сетка будет более густой там, где т принимает большие значения. Пользователь имеет возможность задавать разные или одинаковые управляющие функции для разных частей границы, а также использовать для расстановки узлов на границе такую же управляющую функцию ^(ж1,^2), как и внутри области П (см. соответствующий пример в работе [8]).
2. Начальное приближение
Выше было указано, что в комплексе программ для построения начальных сеток используются два алгебраических метода: STAR и TFI. Пусть область П требуется покрыть сеткой с количеством узлов N х N2. Предположим, что на границе сетка уже построена, при этом на частях T^t и rtop находится Ni узлов, а на Г[е£4 и rright — N2 узлов. Так как граничные узлы входят составной частью в множество узлов двумерной сетки, то для нумерации граничных узлов также будем использовать двойные индексы. Например, узлы, лежащие на rieft, будем обозначать через x1j2, j2 = 1, ... , N2.
При построении начального приближения в звездной области методом STAR вначале вычисляются координаты xp, x^ полюса звездности хр [3], а затем строятся координатные линии первого семейства. Они являются ломаными, задаваемыми различным количеством опорных точек. Одна из ломаных проходит через xp и узлы сетки x1;[n2/2] и xn1,[n2/2], располагающиеся "посередине" на границах Tieft и rright соответственно. Эта ломаная определяется тремя опорными точками, целиком лежит внутри П и делит ее на две части.
Алгоритмы построения других ломаных в обеих частях аналогичны, поэтому рассмотрим построение координатных линий первого семейства в "нижней" части П, т. е. линий с номерами 1 < j < [N2/2]. Опорные точки для j-й линии лежат на отрезках, соединяющих узлы границы Xj £ rieft, Xji,i £ Tbot и XNj £ rright (ji = 1, ... , Ni, j2 = 2, ... , 1j) с полюсом Xp. Всего на j-й координатной линии будет находиться N1 + 2(j — 1) опорных точек. Их координаты вычисляются по формулам
o , Lleft(x1,j2, Xp) /• • w \ • • 1 о
Xj_j2+1 = X1,j2 + [Na/2] — j2 (j — j2)(XP — X1>-b )' J2 = J' J — 1 2
o . L^ot^h,^ Xp) /• -t\/ \ -t tvt
X0i+j_1 = xji>1 + [N2/2] — 1 (j— 1)(xP — J1 =1..., N1,
o Lright (xNi ,j2, xp) /• • w \ D
xN1+j-2+h = xNi,j2 + [Na/2] - j2 (j - j2)(xp - xNi,j2) j2 = ^ . . . , j
где Lieft, Lbot, Lright — длины отрезков, соединяющих соответствующий узел границы с полюсом звездности x^. Далее через опорные точки проводится ломаная, и вдоль нее с помощью одномерного метода эквираспределения расставляются узлы.
Предлагаемый метод STAR построения начальной сетки является более простым по сравнению с методом окаймляющих линий [3] и позволяет строить сетки в произвольной звездной области. В качестве примера на рис. 1, 2 приведены начальные сетки, построенные методом STAR. На рис. 2 показан результат работы STAR для невыпуклой звездной области сложной формы, взятой из статьи [3].
При построении начальной сетки методом TFI в каждом из координатных направлений q1, q2 задаются сечения вычислительной области q^ (/1 = 1, ... , Li) и q22 (/2 = 1, ... , L2), где Lp > 2, причем первые (/^ = 1) и последние (/^ = Lp) сечения совпадают с боковыми границами qe = 0 и соответственно qe =1 области Q (в =1, 2). Затем на каждом сечении определяется вектор-функция A^q), представляющая значение функции x(q) на этом сечении:
All (qii, q2) = x(qii, q2), A* (q1,^ )= x(q1,qi2 ), /p =1,...,L^, в = 1, 2.
На пересечении сечений q1 = q^ и q2 = q22 функции Aii и A^ должны быть согласованы:
Ali (qi\,qi22) = A22 (qi\, q* ), /в = 1, . . . , Le, в = 1, 2.
Формулы трансфинитной интерполяции реализуются в два этапа. На первом этапе интерполируются значения координат узлов сетки с левой и правой границ, а также значения во внутренних сечениях вектор-функций A^:
Li
1 ^а,1^
F(q\q2) = J] < (q1)A1i (qj, q2).
11=1
На втором этапе интерполируются значения координат узлов сетки с нижней и верхней границ:
V2
х(д1, д2) = Г(д1, д2) + £ а|; (д2) (А*, - Г)(д1, ^ ),
¿2 = 1
где ), в = 1, 2 гладкие скалярные функции, зависящие от одной переменной и
удовлетворяющие условиям
а? (<£ ) = 4, =1,-..,Ьв.
Здесь — символ Кронекера.
Функции а^ ) называют переходными функциями. Они служат для переноса значений вектор-функции А^ с границ и внутренних сечений расчетной области вовнутрь. В качестве переходных функций используются полиномы Лагранжа, сплайны, кусочно-линейные функции [25]. Полиномы Лагранжа и сплайны высокого порядка при большом числе сечений приводят к осцилляциям в процессе построения сетки и выходу ее узлов за границы области О. Кусочно-линейные функции обладают другим недостатком, а именно в зонах их изломов получаются сетки с резким изменением шага.
Рис. 1. Сетка, построенная методом STAR. Рис. 2. Сетка для области из работы [3], w = 1 + 10x(1 — x) + y3. построенная методом STAR (w = 1).
В работе предлагается схема построения переходных функций а^ (дв), свободных от недостатков, присущих полиномам Лагранжа, сплайнам высокого порядка и кусочно-линейным функциям. Эта конструкция использует произвольную положительную функцию р(£) : [0, то) ^ [0, то), удовлетворяющую следующим ограничениям:
р(0) = 0, р(1) = 1.
Построение переходных функций осуществляется в два этапа. Вначале для сечений , I = 1, ..., Ьр определяются две серии функций (дв) и ргв(дв), при этом для I = 1
Если 1 < / < Lf, то
pf (qe) = p(1 — qe), 0 < qe < 1.
pf (qe)
p
ql
q в
q в
0 < qe < qf-1; qf-1 < qf < 1.
Аналогично
Если 1 < / < Le, то
pB) = p(9e), 0 < qe < 1.
pB (qв)
0,
Окончательно переходные функции ae (qe) определяются по формуле
ql +1 — ql
в
в
1 > qв > qf+1, 0 < qe < qf+1.
of (qf) = pf (qв)pf (qf), / = 1, ..., Le.
(2.1)
Отметим, что если базисная функция р(£) гладкая и р/(0) = 0, то переходные функции получаются всегда гладкими.
0
в
q
На рис. 3, а представлены гладкие переходные функции для
Ьв = 4, дв = 0, чв = 0.333, дв = 0.666, дв = 1,
представляющие комбинацию гладких отображений на основе базисной функции ^(¿) = 8т2(0.5п£). Для сравнения на рис. 3, б представлены переходные функции, построенные на основе полиномов Лагранжа
ьв в — в
ав(чв) = П , 1 = I. (2.2)
7=1 ч - Ъ
В случае, когда учитываются только границы физической области и не вводятся дополнительные внутренние сечения, мы имеем классическую трансфинитную интерполяцию [19].
На рис. 4, а показана Swan-область [22], граница которой состоит из двух прямолинейных отрезков Г[е£4 и ГЬо4 и двух парабол (ж = 1 + 2у — 2у2, 0 < у < 1) и Г^ор (у = 1 — 3ж + 3ж2, 0 < ж < 1). На рис. 4, б представлена сетка, полученная по формулам классической трансфинитной интерполяции. Эта сетка имеет самопересекающиеся ячейки и потому является неудовлетворительной для расчетов.
Дополнительный контроль за распределением внутренних точек может быть получен путем введения между границами промежуточных сечений. Из рис. 5 видно, что введение всего одного дополнительного внутреннего сечения 1 = 21 позволяет получить сетку без самопересекающихся ячеек.
3. Уравнения для вычисления координат узлов
После того, как начальная сетка построена, производится ее "улучшение". Для этого используются итерационные методы решения уравнений, которым должны удовлетворять координаты узлов сеток. В настоящем комплексе программ окончательные сетки строятся на основе двумерного метода эквираспределения (ЕБ2-метода) и ИЬ-метода из работы
[24].
Основы ЕБ2-метода изложены в статье [16], в работах [10, 12] уравнения этого метода получены непосредственно из двумерного варианта принципа эквираспределения. Этот метод является развитием идеи [15, 24] о возможности построения сеток путем решения в вычислительной области Q = (0,1) х (0,1) уравнений эллиптического типа для декартовых координат узлов жв (в = 1, 2). ЕБ2-уравнения имеют следующий вид:
д / джв\ д / джв\ , 1 2ч ^ , .
ддт 1^22 д^) + дф) = 0 а =(д,д) € Q, (3Л)
где gee — ковариантные компоненты метрического тензора искомого преобразования координат (1.1), w(x) > 1 — управляющая функция, x = (ж1, ж2) G П. Уравнения (3.1) получены из уравнений для конструирования ортогональных сеток, приведенных в [24], и условия постоянства произведения wJ, где J — якобиан отображения (1.1). На разностном уровне это требование сводится к свойству wS ~ const, где S — площадь произвольной ячейки криволинейной сетки. Следовательно, задавая функцию w меняющейся в области П, мы должны получить неравномерную сетку со сгущениями и разрежениями узлов в
Рис. 3. Переходные функции ав(дв). а — гладкие переходные функции (2.1), б — полиномы Лагранжа (2.2).
Рис. 4. 8'теап-область (а) и сетка, построенная классическим методом трансфинитной
интерполяции (б).
Рис. 5. Вычислительная область с дополнительным внутренним сечением (а) и сетка,
построенная ТЕТ-методом (б).
подобластях соответственно больших и малых значений функции и>(х). На рис. 6, 7 показаны примеры окончательных сеток, построенных на основе ЕБ2-метода с использованием начальных приближений, изображенных на рис. 1 и 2 соответственно.
Рис. 6. Сетка, построенная ЕБ2-методом (т = 1 + 10х(1 - х) + у3).
Рис. 7. Сетка для области из работы [3], построенная ЕБ2-методом (т = 1).
ЯЬ-уравнения имеют следующий вид:
д / дхв\ д /1 дхв \
V д^У + \1 = ' ( ' )
где управляющая функция А(х) является отношением длин смежных сторон ячейки. В работах [17, 20] отмечаются сложности, связанные с построением удовлетворительных сеток на основе уравнения (3.2). В основном они вызваны тем, что функция А должна выбираться специальным образом до начала итераций либо подправляться в ходе итераций, чтобы полученная сетка была близка к ортогональной. В работе [2] отмечается, что вопрос выбора функции А тесно связан с вопросом о количестве краевых условий на функции хв.
Опыт работы с комплексом программ показал, что ЯЬ-метод очень чувствителен к выбору функции А(х). Во многих тестовых примерах с невыпуклыми областями сетку удавалось построить ЯЬ-методом только после кропотливой работы по подгонке функции А. ЕЭ2-метод является в этом смысле более универсальным, он строит сетки практически для любой функции т(х) > 1.
4. Критерии оценки качества сетки
Рассмотрим теперь критерии, по которым можно судить о качестве сетки. Авторы считают, что окончательную оценку сетке можно дать лишь после решения на ней предметной задачи. Однако количественная информация о сетке (различные критерии качества) может оказаться полезной для исследователя при предварительной оценке сетки, позволяя ему отбросить заведомо непригодные из них до решения основной задачи. В настоящей работе описываются несколько критериев, в том числе два критерия из [9].
Рассмотрим четыре узла сетки х^¿2, Xj1+1^2, Xj1+1^2+1, Xj1 ¿2+1 и, следуя [9], присвоим им номера 1, 2, 3, 4. Ячейка своими диагоналями разбивается на две пары треугольников. Вычислим площади этих треугольников, приписывая каждой из них положительное значение, если обход треугольника в сторону возрастания номеров вершин совершается против часовой стрелки, и отрицательное значение в противном случае. Величина /(!\/2
полагается равной минимальному из четырех отношений площадей этих треугольников к половине площади рассматриваемой ячейки, которую будем нумеровать мультииндек-сом ] + 1/2 = (^1 + 1/2,^2 + 1/2). Значения /(+)1/2 могут лежать в интервале (-то; 1], при
этом для выпуклой ячейки будем иметь 0 < /(+>1/2 < 1, для ячеек, вырождающихся в
треугольник, — /+1/2 = 0, для невыпуклых — /+1/2 < 0. Наибольшее значение /(+)1/2 = 1 достигается только в том случае, когда ячейка является параллелограммом. С помощью этого подхода определяются и самопересекающиеся ячейки [9].
Второй критерий связан с оценкой ортогональности координатных линий сетки. Для
г(2)
каждой ячейки вычисляется величина /.,-+1/2, равная минимальному значению синусов
(2
внутренних углов. Функция /)+1/2 может принимать значения из интервала [-1; 1]. Наибольшее значение /?(+)1/2 = 1 достигается в том случае, когда ячейка является прямоугольником [9].
Одним из основных требований к сеткам, рассматриваемым в настоящей работе, является адаптация сетки к априорно заданной управляющей функции и>. Причем адаптация здесь понимается в смысле адаптации площадей ячеек Б-+1/2:
(™%ц/2 = С, (4.1)
при этом постоянная вычисляется по формуле
М1-1,М2-1 I М1-1,М2-1
С = .7= Б+"»/ - 1 —V <4-2)
1, 2=1 1, 2=1
Это требование учитывалось, например, при выводе уравнения (3.1) ЕБ2-метода. Поэтому было бы интересно проверить выполнение соотношения (4.1) для окончательной сетки, а также близость (в смысле выполнения (4.1)) начальной сетки к адаптивной. Количественную оценку можно дать с помощью третьего критерия:
Г(3) -тщ/(^)-7+1/2 _С_1 (43)
/+1/2 = т1^ С , (^Б)^/. (4.3)
Для данной ячейки величина /7(+)1/2 характеризует меру отклонения ее площади от предписанного соотношением (4.1) значения Б. Функция /7(+)1/2 может принимать значения из интервала [0; 1], при этом малое значение /7(+)1/2 свидетельствует о сильном несоответствии площади ячейки закону адаптации, а при выполнении (4.1) получается значение
/(3) = 1
/-7 + 1/2 = 1
Величина /7(+)1/2 служит для оценки вытянутости ячеек. Ее значение равно отношению длины наименьшей к длине наибольшей из сторон ячейки, поэтому 0 < /7(+)1/2 < 1.
Значение /7(+)1/2 = 0 соответствует случаю вырождения ячейки в треугольник или отрезок. Для обычного ромба, а также невыпуклого вырожденного в уголок ромба будем
иметь /7(+)1/2 = 1. В последнем случае первые два критерия покажут "плохие" значения: /(1) < 0, /(2) < 0.
/7+1/2 < 0, /7+1/2 < 0.
После того, как значения критериев для каждой ячейки вычислены, определяются их
Лк) „ «
значения для сетки в целом, при этом они полагаются равными / +1/2 в самой "худшей"
ячейке [9]:
f(k) = . min 1 fj+1/2, в = 1,2, k =1,2,3,4. (4.4)
Зв = 1,...,Ne-1 JT '
Отметим, что окончательная сетка, полученная после итераций, должна быть "лучше" начальной и критерии f(k) должны уловить это улучшение. Вместе с тем может случиться, что значения одного или нескольких критериев для начальной и окончательной сеток одинаковы. Это, однако, не всегда означает, что окончательная сетка не стала лучше начальной. Может случиться, что "худшая" ячейка начальной сетки не улучшила свои характеристики, но тем не менее в других ячейках значения критериев стали большими. Поэтому имеет смысл кроме критериев по худшей ячейке ввести в рассмотрение осред-ненные критерии, характеризующие качество сетки в среднем:
Ni-1,N2-1 /
fOk) = £ fj+W (N1 - 1)(N - 1), k =1,2,3,4. (4.5)
31,32 = 1 /
5. Описание программы
Описанные выше алгоритмы реализованы в виде программного комплекса, написанного в основном на алгоритмическом языке Watcom Fortran 77/32 Optimizing Compiler Version 10.5. Некоторые модули (работа с манипулятором "мышь", построение гладких кривых по опорным точкам и т.д.) реализованы на алгоритмическом языке WATCOM С32 Optimizing Compiler.
Основная часть программного комплекса состоит из программ построения начального приближения и программ построения окончательной сетки с помощью итерационных методов решения эллиптических уравнений. При запуске комплекса появляется приглашение к интерактивной работе по описанию частей rieft, T^ot, rrjght и rt0p границы области П, для которой будет строиться сетка. Если область уже была описана ранее, то достаточно указать соответствующий файл. Для новой области такой файл с описанием границы будет создаваться в ходе интерактивной работы с комплексом. Если какая-то из указанных частей границы является отрезком прямой, то для ее описания достаточно ввести координаты начала и конца этого отрезка. Для границы, заданной таблично, пользователь вводит координаты опорных точек. Для аналитического задания некоторых частей границы пользователь получает доступ к строго определенным Fortran-подпрограммам, в которые ему необходимо внести изменения в соответствии с конкретными параметрическими зависимостями (1.2). Работа этой части программы заканчивается созданием и сохранением нового файла опорных точек границы и отрисовкой изображения области. При неудовлетворительном результате можно вернуться к началу и исправить описания некоторых или всех частей границы.
Далее работают программы построения начальных сеток. По умолчанию управляющие функции для сетки на границе и внутри области имеют постоянное значение, равное единице. Если необходимо строить сетку со сгущениями, то предварительно исправляются управляющие функции в соответствующих Fortran-подпрограммах, доступных для пользователя. Для каждого метода полученные сетки изображаются с помощью графических процедур компилятора Watcom Fortran 77/32 и пакета SMOG [13]. Кроме этого, на экран выводятся значения критериев качества сеток.
Первой работает программа, в которой реализован метод STAR. Сначала исходная область проверяется на выпуклость. Если она не выпукла, то проводится дополнительная
проверка на звездность. Если область является выпуклой или звездной, то определяется полюс звездности и строится начальная сетка. В случае успешной работы программы STAR координаты узлов начальной сетки записываются в файл.
Далее происходит переход к интерполяционному TFI-методу. Для построенной сетки проверяется, не будет ли какой-либо отрезок, соединяющий два соседних узла на одной и той же внутренней координатной линии, пересекаться со звеньями ломаных, которыми представляются границы. В случае пересечения вырабатывается соответствующий признак. Если некоторые внутренние координатные линии пересекают границу, то часть узлов начальной сетки может оказаться за пределами области. Такое начальное приближение считается плохим. Тем не менее, как было показано ранее [8, 10], ЕВ2-метод может работать и с таким начальным приближением, поэтому оно также сохраняется в файле.
Неудовлетворительную сетку, построенную с помощью классического TFI-метода, можно подправить, введя несколько внутренних сечений. Построение сечения начинается с указания курсором "мыши" первой опорной точки в том месте границы, где будет начинаться сечение. Далее с помощью "мыши" устанавливаются опорные точки сечения внутри области и последняя точка на противоположной границе. Результатом работы этой части программы будет изображение опорных точек сечения, соединенных отрезками прямых, и гладкой кривой, построенной по этим точкам с помощью B-сплайна.
При необходимости можно изменить форму сглаженной кривой. Для этого следует с помощью "мыши" переместить некоторые из опорных точек, что приведет к изменению этой кривой. Окончание работы по построению данного сечения происходит при нажатии левой клавиши "мыши" в произвольной точке рабочего поля, не совпадающей с опорными точками. Координаты опорных точек B-сплайнов всех сечений записываются в файл данных и используются в модифицированном TFI-методе.
Далее работает блок сравнения начальных сеток. Здесь вначале отбираются сетки с условием f(1) > 0, т. е. сетки с выпуклыми ячейками. Для всех таких сеток f(2) > 0. В полученном множестве сеток выбирается подмножество, состоящее из сетки S с максимальным значением f(2) и сеток, для которых значения осредненных критериев f2) больше осредненного критерия f2) сетки S. Аналогичное сравнение проводится по критериям f(3) и f3) и в последнюю очередь — по f(4) и f4). Если в результате такого отбора последнее подмножество все еще состоит более чем из одной сетки, то производится окончательный отбор путем последовательного сравнения значений f(2), f(3) и f(4).
Если для всех сеток значение f(1) < 0, то для дальнейшего отбора оставляются сетки без самопересекающихся ячеек и без пересечений внутренних координатных линий с границами. Если такое множество не пусто, то, в отличие от предыдущего, сравнение начинается с критерия f(1): оставляются сетка S с максимальным значением f(1) и сетки, у которых средние характеристики fo превышают среднюю характеристику fo сетки S. Далее отбор идет как для сеток с выпуклыми ячейками.
Наконец, в случае, когда все начальные сетки имеют самопересекающиеся ячейки или внутренние координатные линии, пересекающиеся с границами, отбор идет путем последовательного сравнения средних характеристик fifc), k = 1, 2, 3, 4.
В любом случае будет выбрана одна из начальных сеток. Ее координаты будут прочитаны из файла, а на экране появится изображение выбранной автоматически начальной сетки, "лучшей" в описанном выше смысле. Отметим, что пользователь имеет возможность выбрать по своему усмотрению другую из построенных начальных сеток.
Окончательная сетка строится на основе итерационного решения эллиптических уравнений, рассмотренных в разд. 3. При этом на каждой из итераций вычисляются критерии
Рис. 8. Начальная (а) и окончательная (б) сетки для Swan-области (и> = 1).
качества сетки и если эти критерии перестают изменяться от итерации к итерации, то итерационный процесс прекращается. По окончании итерационного процесса на экране отображаются полученная окончательная сетка и итоговые критерии качества.
На рис. 8, б показана окончательная сетка, полученная из начальной сетки, изображенной на рис. 8, а. Окончательная сетка стала более оргональной, на что указывает критерий ортогональности: 0.22 в начале против 0.33 в конце итераций.
6. Пример расчета на различных сетках
В качестве области О в тестовой задаче брался криволинейный четырехугольник, ограниченный слева и справа прямыми х = 0 и х = 1, снизу — прямой у = — 1, сверху — кривой у = п(х) = —аеов(2пх). В этой области решалась смешанная краевая задача для уравнения Лапласа, при этом на верхней части границы задавалось условие Дирихле, а на других частях — условие Неймана. В качестве точного решения задачи была выбрана функция
Рех(х,у) = —— еовЬ[2п(у + 1)] еов(2пх). еовп(2п)
Уравнение Лапласа для функции ^ и краевые условия были переписаны в новых координатах д1,д2. Полученное уравнение аппроксимировалось со вторым порядком схемой типа "косой крест" [11] на девятиточечном шаблоне.
Сначала были проведены расчеты на криволинейных сетках простой структуры, адаптирующихся только к границам области (рис. 9, а). Одно семейство координатных линий состояло из отрезков прямых, соединяющих нижнюю и верхнюю части границы. Каждый из этих отрезков делился на одно и то же число элементарных отрезков равной длины. Кривые, соединяющие точки деления, образовывали второе семейство координатных линий. Таким образом, для построения сетки вначале использовался метод вертикальной интерполяции [8], основанный на "спрямляющем" преобразовании координат вида
х = д1, у = —1 + д2(п(х) + 1).
На таких "квазиравномерных" сетках были проведены расчеты с последовательным удвоением числа узлов, которые показали, что разностная схема [11] имеет второй порядок точности.
Использованная сетка адаптируется к границе области, но не учитывает особенностей решения и соответственно имеет низкое значение критерия адаптивности (0.06). Решение
Рис. 9. Квазиравномерная (а) и адаптивная (б) сетки.
рассматриваемой тестовой задачи очень быстро меняется в верхней части области О, поэтому сетка должна иметь здесь мелкие ячейки. На рис. 9, б показана адаптивная сетка, построенная БВ2-методом с использованием управляющей функции вида
ы(х,у) = а/1 + ао(х,у)\ + а1\У^ех(х,у)\.
Построенная сетка лучше адаптируется к решению (значение критерия адаптивности равно 0.29), чем квазиравномерная сетка, но имеет худшие значения других критериев, в особенности критерия вытянутости. Тем не менее даже на такой неудовлетворительной адаптивной сетке точность численного решения повышается. Так, при а0 = 0.5, а1 =2 точность повышается практически в три раза по сравнению с квазиравномерной сеткой с таким же количеством узлов. Еще большего увеличения точности численного решения можно добиться подбором соотношения между числами N1 и N2 (N1 > N2), при котором значение критерия вытянутости адаптивной сетки лишь незначительно уступает этому же критерию квазиравномерной сетки с тем же количеством узлов N1 х N2.
7. Заключение
Создан комплекс программ, позволяющий в интерактивном режиме общения с пользователем конструировать адаптивные сетки в двумерных односвязных областях. В комплексе реализованы алгоритмы описания области, задания ее внутренних сечений, несколько алгебраических и дифференциальных методов построения сеток и алгоритм оценки качества получающихся сеток.
Программный комплекс имеет модульную структуру и вследствие открытой архитектуры допускает простое расширение за счет добавления в него новых модулей с реализациями других алгебраических и дифференциальных методов построения двумерных структурных сеток.
Список литературы
[1] Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. 1—2, Мир, М., 1990.
[2] ВольцингЕр Н.Е., Клеванный К. А., ПЕлиновский Е. Н. Длинноволновая динамика прибрежной зоны. Гидрометеоиздат, Л., 1989.
[3] ГАСИЛОВА И. А. Алгоритм автоматического построения начального приближения криволинейной сетки для областей звездного типа. Вопросы атомной науки и техники. Сер. Матем. моделир. физ. процессов, вып. 3, 1994, 33-40.
[4] ГОДУНОВ С. К. И ДР. Численное решение многомерных задач газовой динамики. Наука, М., 1976.
[5] ЛИСЕЙКИН В. Д. Технология конструирования трехмерных сеток для задач аэрогазодинамики. Вопросы атомной науки и техники. Сер. Матем. моделир. физ. процессов, вып. 3, 1991, 31-45.
[6] ЛИСЕЙКИН В. Д. Обзор методов построения структурных адаптивных сеток. Журн. вычисл. матем. и матем. физ., 36, № 1, 1996, 3-41.
[7] МЕТОДЫ расчета обтекания элементов летательных аппаратов при трансзвуковых скоростях Ч. II. Методы расчета сеток. ОНТИ ЦАГИ, М., 1989.
[8] МОЛОРОДОВ Ю. И., ХАКИМЗЯНОВ Г. С. Построение и оценка качества регулярных сеток для двумерных областей. Вопросы атомной науки и техники. Сер. Матем. моделир. физ. процессов, вып. 1, 1998, 19-27.
[9] ПРОКОПОВ Г. П. Об организации сравнения алгоритмов и программ построения регулярных двумерных разностных сеток. Там же, вып. 3, 1989, 98-108.
[10] ХАКИМЗЯНОВ Г. С., ШОКИНА Н. Ю. О методе эквираспределения для построения двумерных адаптивных сеток. В "Вычисл. технологии". ИВТ СО РАН, Новосибирск, 4, №13, 1995, 271-282.
[11] ХАКИМЗЯНОВ Г. С., ЧУБАРОВА Э. В. О некоторых способах аппроксимации уравнения для потенциала. Вычисл. технологии. 2, №5, 1997, 82-90.
[12] ХАКИМЗЯНОВ Г. С., ШОКИНА Н. Ю. Метод эквираспределения для построения адаптивных сеток. Там же, 3, №6, 1998, 63-81.
[13] MATSOKIN A. M., DEBELOV V. A., SlROTIN V. G., UPOLNIKOV S.A. Multi-Purpose Computer Graphics System SMOG-85. In "Comput. and Graphics". Pergamon Press, NY, 12, No. 3/4, 1988, 441-456.
[14] ФЛЕТЧЕР К. Вычислительные методы в динамике жидкостей. 2, Мир, М., 1991.
[15] AMSDEN A. A., Hirt C. W. A simple scheme for generating general curvilinear grids. J. Comp. Phys, 11, 1973, 348-359.
[16] CHRISTOV C. I. Orthogonal coordinate meshes with managable jacobian. In "Numerical Grid Génération; Appl. Math, and Comp.", 10/11, 1982, 885-894.
[17] DURAISWAMI R., PROSPERETTI A. Ortogonal mapping in two dimensions. J. Comp. Phys. 98, 1992, 254-268.
[18] ERIKSSON L. E. Practical three-dimensional mesh generation using transfinite interpolation. SIAM J. Sci. and Statist. Comput., 6, No. 3, 1985, 712-742.
[19] GORDON W. J., Thiel L.C. Transfinite mappings and their applications to grid generation. In "Numerical Grid Generation; Appl. Math. and Comp.", 2/3, 1982, 171192.
[20] KANG I.S., Leal L.G. Ortogonal grid generation in a 2D domain via the boundary integral technique. J. Comp. Phys., 102, 1992, 78-87.
[21] Kim H.J. AND Thompson J.F. Three-dimensional adaptive grid generation on a composite-block grid. AIAA J., 28, No. 3, 1990, 470-477.
[22] Knupp P., Steinberg S. Fundamentals of grid generation. CRC Press, 1994.
[23] Liseikin V. D. Survey of grid generation technology. Advanced Mathematics: Comp. and Appl. Proc. of AMCA-95, Eds. A. S. Alexeev, N. S. Bakhvalov. Novosibirsk, 1995, 511-517.
[24] Ryskin G., Leal L.G. Orthogonal mapping. J. Comp. Phys., 50, 1983, 71-100.
[25] Thompson J.F., Warsi Z.U.A., Mastin C.W. Numerical grid generation, foundations and applications. North-Holland, Amsterdam, 1985.
[26] THOMPSON J.F. A survey of dynamically adaptive grids in the numerical solution of partial differential equations. Appl. Numer. Math., 1, 1985, 3-28.
[27] THOMPSON J.F. A composite grid generation code for general 3-D regions — the EAGLE code. AIAA J., 26, 1988.
[28] Thompson J.F., Warsi Z.U.A., Mastin C.W. Numerical grid generation: Foundations and Applications. Electronic Edition. WWW.erc.msstate.edu, 1997.
nocmyuuAa e peda^uw 25 deKaSpa 1998 г.
