Конструирование разностных сеток с помощью численного решения обращенных уравнений Бельтрами и уравнений диффузии Текст научной статьи по специальности «Математика»

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

Том 11, часть 2, Специальный выпуск, 2006

КОНСТРУИРОВАНИЕ РАЗНОСТНЫХ СЕТОК С ПОМОЩЬЮ ЧИСЛЕННОГО РЕШЕНИЯ ОБРАЩЕННЫХ УРАВНЕНИЙ БЕЛЬТРАМИ II УРАВНЕНИЙ ДИФФУЗИИ*

В. Д. ЛИСЕЙКИН, И. А. Васева, Ю.В. Лиханова Институт, вычислительных технологий СО РАН, Новосибирск, Россия, e-mail: lvd@ict.nsc.ru, vaseva@ngs.ru, ula@gorodok.net

This paper describes a comprehensive method based on the solution of inverted Beltrami and diffusion equations in monitor metric. Examples of metrics providing generation of adaptive, aligned to vector fields, and balanced numerical grids are demonstrated.

Введение

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

Разработкой методов и компьютерных программ автоматического построения разностных сеток занимаются в основном в США с начала 90-х годов [1] в рамках глобальной национальной инициативы по созданию мощных компьютерных центров и автоматизированных пакетов компьютерных программ численного решения прикладных задач [2]. Их компьютерные программы для конструирования структурных сеток базируются на численном решении обращенных уравнений Пуассона и на фронтальном и методах Делоне для генерирования неструктурных сеток. Уравнения Пуассона обладают существенными

*Работа выполнена при финансовой поддержке Комплексного интеграционного проекта СО РАН (2006, № 1.8) и Российского фонда фундаментальных исследований (грант № 06-01-08009).

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

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

С начала 2000-х годов в Институте вычислительных технологий разрабатывается более перспективный метод конструирования сеток [3]. Он базируется на численном решении обращенных уравнений Бельтрами и более общих уравнений диффузии относительно управляющей метрики. Уравнения Бельтрами позволяют единообразно и независимо от параметризации физической геометрии произвольной размерности строить как структурные, так и неструктурные разностные сетки со свойствами, требуемыми в прикладных задачах. Такие свойства легко реализуются с помощью общей формулы управляющей метрики, описанной в [3]. Некоторые последние результаты, представленные в данной статье, показывают, что этот метод более эффективен и перспективен по технологичности, универсальности и возможности автоматизации по сравнению с методами, базирующимися на численном решении обращенных уравнений Пуассона.

1. Схема метода отображений

Разрабатываемый авторами пакет компьютерных программ построения разностных сеток основывается на технологии метода отображений, который формулируется единообразно для произвольной п-мерной физической геометрии Бхп С 11”+ (кривая, область, поверхность), задаваемой локально параметризацией

где Б” — п-мерная параметрическая область (отрезок при п = 1); х(б) — гладкая вектор-функция ранга п в каждой то чке б € Б ”. Пр и к = 0 поверхнос ть Бхп является областью

х(б) : Б” ^ 1Г+к, х = (х1 х”+к), б = (в1,в”), п > 1,

(1)

Рис. 1. Схема метода отображений для построения сеток с треугольными ячейками.

X” С И”. В этом случае в качестве параметрической области можно рассматривать X”.

Бхп

помощи промежуточного невырожденного преобразования

8(*):5” ^ Б” , * = (£\...,С), (2)

между параметрической областью Б” и соответствующей вычислительной областью 5” (рис. 1). При этом узлы сетки в Бх” определяются отображением эталонной сетки, заданной в 5”, при помощи преобразования

х[б(*)] : 5” ^ Бх” С 1Г+к. (3)

Форма вычислительной области 5” и ячеек эталонной сетки выбирается в зависимости

Бх”

задачи. В частности, эталонная сетка в 5” может быть неструктурной.

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

Основными инструментами для построения разностных сеток в предлагаемом методе отображений являются операторы Бельтрами и диффузии. Эти операторы формулируются

Бх”

управляющей метрикой. Будем обозначать ковариантные (контравариантные) компоненты управляющей метрики в координатах V1,..., V” через д- (д-).

Оператор Бельтрами Дв[ ] для функции f(б), заданной в локальных координатах в1,..., в” физической геометрии Бх” с введенной управляющей метрикой д— определяется следующим образом:

где дй = ). Известно [3], что значение оператора Бельтрами для заданной функ-

ции f (б) не зависит от параметризации физической геометрии Бх”, поэтому оно является инвариантом, называемым вторым дифференциальным параметром Бельтрами функции / (б).

Промежуточное преобразование б(*) между вычислительной и параметрической областями в данной версии метода отображений определяется как обратное к вектор-функции *(б) : Б” ^ 5”, являющейся решением задачи Дирихле:

Дв [*] = °,

(5)

*|с№ = ^(8).

Краевая задача (5) корректна при любой невырожденной метрике д— А так как значение оператора Бельтрами Дв [*] инвариантно относительно выбора координатной системы в Бх”, сетка, полученная с помощью преобразования (3), где в(£) — это функция, обратная к решению задачи (5), не зависит от выбора параметризации.

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

Бх”

когда вычислительная область 5” выпукла. Более того, при п = 2 справедлива теорема Радо, из которой, в частности, следует, что преобразование *(б), полученное решением задачи (5) с произвольной метрикой, является невырожденным, если область 52 выпукла и отображение границ Б2 и 52, заданное условием Дирихле, взаимно-однозначное. В этих свойствах основная причина того, что в задаче (5) базисная система уравнений Бельт-рами сформулирована для преобразования *(б), обратного к промежуточному отображению б(*). Для нахождения преобразования искомого б(*) численно решаются обращенные уравнения, которые выписаны ниже.

Уравнения в (5) эквивалентны системе уравнений Эйлера — Лагранжа для функцио-

НШТс!

= = (6)

называемого функционалом энергии отображения *(б) : Б” ^ 5” или обобщенным функционалом гладкости разностной сетки.

Выражение функционала (6) подсказывает один из возможных путей определения управляющей метрики для построения разностных сеток, удовлетворяющих какому-либо требованию. Для этого нужно сформулировать контравариантный метрический тензор д-к управляющей метрики в Бх” таким образом, чтобы подынтегральное выражение в (6) описывало меру отклонения разностной сетки от желаемой сетки с требуемым свойством в точке 8 € Б”. Если такой метрический тензор найден, то при минимизации функционала

(6) можно ожидать получения разностной сетки с этим свойством.

Замена ~д* на ад (в) в (5) приводит к более общей краевой задаче

-^-т(т(5)д1к-^г\ = 0, А; = 1,..., п,

дв- V у двк/ >./>>>> ^

*|б№ = ^(б), г = 1, ...,п,

где ю(б) > 0 — весовая функция, усиливающая или уменьшающая влияние метрики в нужных зонах Бх”. Уравнения в (7) будем называть уравнениями диффузии.

Для нахождения искомого преобразования б(*) необходимо обратить уравнения в (5) и (7). В результате задача Дирихле (7) преобразуется в задачу Дирихле для обращенных уравнений диффузии:

д2 д / -А

ы(ъ)дк^——— = ——(т(з)д1г), 1^,к,т = \,...,п,

д£кд£т дв-V у ’ ’ ’ ’ (8)

в^|се™ = ^(£)> г = 1, ...,п,

где -0*(*) — г-я компонента преобразования, обратного к ^(б). Численное решение задачи (8) в точках эталонной сетки, задаваемой в вычислительной области 5”, определяет

Б”

разования х(б) задает сетку в физической геометрии Бх”.

Если в задаче (8) сделать замену ад (в) = у'//*- получим краевую задачу для обращенных уравнений Бельтрами.

3. Управляющая метрика

Для эффективного управления свойствами сеток с помощью уравнений Бельтрами или

Бх”

рик. Эти управляющие метрики определяются через переменные физической задачи, гео-

Бх”

ношению к которым сетка должна адаптироваться.

Бх”

в1,..., в” имеет следующий вид:

д« = Ф)ду + ^(8) , г, ] = 1,..., п , к = 1,..., I, (9)

где ^(б) > 0 — весовая функция; д— = х^ ■ х^ — метрика Бх”, а Рк(б) г = 1,..., п, — компоненты некоторого ковариантного вектора Ек (б) [3].

Функции ^(б) и Рк(б) в (9) должны удовлетворять условию с1е1;(д-) > 0. В частности, ^(д—) > 0 при ^(б) > 0.

Функционал энергии (6), так же как и уравнения Бельтрами и уравнения диффузии, зависит от величины дй = (д—) и контравариантых метрических компонент д— в коор-

динатах в1в”. Поэтому вместо ковариантных метрических компонент д— можно, если это удобно, изначально задавать контравариантные компоненты управляющей метрики д—. Это естественно сделать, например, в виде (9), а именно:

д- = Ф)д—х + ркрк, г, 2 = 1,..., п , к = 1.,..., /, (Ю)

где с(б) > 0, Рк, г = 1,..., п, — компоненты коптравариантного вектора Рк = (Рк,..., Р”), к = 1,..., I. В этом случае величину д& в уравнениях Бельтрами можно определить из соотношения д& = 1/д8, где д8 = йе1(д^-).

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

3.1. Метрика для построения сеток, сгущающихся в зонах больших градиентов функции

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

Бх”

функции f(б). Эта поверхность задается при помощи следующей параметризации:

г(8): Б” ^ #”+к+г , г(8) = [х(8),/ 1(б),...,/ 1 (8)] , (11)

где х(б) — параметризация физической геометрии Бх”; f (б) = (/ 1(б), ..., /1 (б)) — вектор-функция, для которой требуется построить разностную сетку, сгущающуюся в зонах, где эта функция имеет большие градиенты. Поверхность, задаваемая параметризацией (11), обозначается через Бг” и называется мониторной поверхностью над Бх”. Управляющей метрикой д— обеспечивающей сгущение узлов сетки в зонах больших градиентов, является метрика мониторной поверхности Бг”. Ковариантные метрические компоненты д-мониторной поверхности Бг” в координатах в1,..., в” вычисляются по формуле

д— = г5; ■ г^' = д- + /к/к, г,2 = 1,...,п, к = 1,...,/. (12)

' шт.

лйр

Рис. 2. Примеры разностных сеток, адаптированных к градиентам функции.

Эта управляющая метрика имеет вид (9) при я (в) = 1, ^ = /*, к = 1,...,1. Для метрики (12) в случае, когда Ї(в) является скалярной функцией /(в), контравариантные компоненты метрики и якобиан д8 вычисляются следующим образом:

1

к д/ д/

<5^

9

О

двк <^те

і,І,1,Р = 1,...,п,

Узхдз1 дзк

(13)

і, 1,1, Р = 1,...,п.

На рис. 2 представлены примеры сеток, построенных для метрики (12) с помощью решения краевой задачи (8) при ги(в) = л/д®.

3.2. Метрика для построения сеток, адаптирующихся к значениям функции

Для построения сеток на физической геометрии Бх” с узлами, сгущающимися в зонах больших значений какой-либо функции ^(б), например погрешности решения, мера отклонения сетки от требуемой может быть выражена в виде

<т(8) = ФМ(8 1,к,т = \,...,п,

где Ф[^] > 0 положительный оператор, такой что Ф[^](б) принимает большие (малые) значения там, где ^(б) мало (велико). В этом случае контравариантные компоненты управляющей метрики имеют следующий вид:

д- = ф[и](8)д— г,2' = 1,...,п.

Для построения таких сеток решаются уравнения диффузии (7) для ,ш (б) = 1.

(14)

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

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

На рис. 3 представлены примеры сеток, построенных с помощью решения краевой задачи (8) при w (б) = 1 для метрики (14).

3.3. Метрика для построения сеток, согласованных с векторным полем

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

Использование такой метрики эффективно для построения сеток, согласованных с магнитным полем, при численном решении задач с сильной анизотропией, в частности, задач теории плазмы [4]. Контравариантные элементы требуемой управляющей метрики в случае, когда Бхп = Бп, можно задавать следующим образом [5]:

д- = е(8)<— + РкР-, г, 3 = 1,... ,п , к = 1,... ,1. (15)

На рис. 4 представлен пример сетки, согласованной с магнитным полем, построенной с помощью численного решения задачи (8) в осесимметричном сечении камеры токамака.

3.4. Метрика для построения сбалансированных сеток

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

д— = е1 (8)да( + е2(з)д— + еэ(з)д—, г, 3 = 1,..., П (16)

Рис. 4. Линии магнитного поля (слева) и сетка, согласованная с магнитным полем в сечении камеры токамака.

Рис. 5. Сбалансированные сетки.

где 6j(s) > 0, i = 1, 2, 3, — весовые функции, отвечающие за вклад контравариантных элементов gal ? gads и gadv, задаваемых формулами (15), (13) и (14) соответственно. Аналогично выписывается ковариантный метрический тензор gij для построения сбалансированных разностных сеток.

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

gj (s) = (1 - a)ga + a((1 - e )gadg + ^dv) > = 1,...,n. (17)

На рис. 5 представлены примеры двумерных сбалансированных сеток с треугольными ячейками, построенных при помощи численного решения задачи (8) в метрике (17). Слева на рисунке изображена сетка, адаптирующаяся к значениям функции ^i(s) и согласо-

^ / д0 д0 ч

ванная с векторным полем В = — —, ——- , где w(s) — модельная функция. В центре

V ds2 ds1/

показана сетка, согласованная с тем же самым векторным полем и адаптирующаяся к градиентам функции f2 [^2(s)j. Справа на рисунке демонстрируется сетка, согласованная с этим же полем и адаптирующаяся к значениям одной функции и градиентам другой.

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

[1] Thompson J.F., Soni В.К., Thompson J.F. ет al (eds). A reflection on grid generation in the 90s: trends, needs influences // Numerical Grid Generation in CFD. Mississippi State Univ., 1996. Vol. 1. P. 1029-1110.

[2] Raveche H.J., Lawrie D.H., Despain A.M. A National Computing Initiative. Philadelphia: SIAM, 1987.

[3] Liseikin V.D. A computational differential geometry approach to grid generation. Berlin: Springer-Verlag, 2004.

[4] Glasser A.H., Tang X.Z. The SEL macroscopic modeling code // Comp. Phys. Comm. 2004. Vol. 164. P. 237-243.

[5] Методы римановой геометрии в задачах построения разностных сеток / Ю.И. Шокин, В.Д. Лисейкин, А.С. Лебедев и др. Новосибирск: Наука, 2005.

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