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

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

Том 16, № 5, 2011

Применение метода конечных элементов для построения адаптивных сеток*

И. А. Васева, В. Д. Лисейкин Институт вычислительных технологий СО РАН, Новосибирск, Россия

е-таП: vaseva.irina@gmail.com

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

Ключевые слова: адаптивные сетки, метод конечных элементов, уравнения Бельтрами, уравнения диффузии.

1. Метод отображений и управляющая метрика

Метод отображений формулируется для произвольной n-мерной физической геометрии Sхп С Rn+fc, Для построения сеток в двумерных областях и на поверхностях (n = 2, k = 1) физическая геометрпя Sх2 задается при помощи параметризации

x(s) : S2 ^R3, x =(x1,x2,x3), s =(s1,s2), (1)

где S2 — параметрическая область, x(s) — вектор-функция, дважды дифференцируемая в каждой точке s G S2 (рис. 1).

Согласно методу отображений сетка в физической геометрии Sx2 строится при помощи промежуточного невырожденного преобразования

s(|):~2 ^ S2, | =(£\£2) (2)

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

x[s(£)] : S2 ^ Sx2 CR3. (3)

Переменные 2 облает и S2 называются сеточными координатами, а величины s^s2

в (1) — параметрическими координатами. В дальнейшем будем предполагать, что s1, s2 —

S2

*Работа выполнена при финансовой поддержке РФФИ (грант № 10-01-00335-а) и Проекта фундаментальных исследований объединенного ученого совета по математике и механике СО РАН (№ 94).

Обозначим через дХ* (д*Х) ковариантные (контравариантные) элементы физической геометрии Бх2 в координатах в1, в2. Тогда параметризация (1) дает формулу для элементов ковариантного метрического тензора физической геометрии Бх2 в координатах

дХ* = х*; ■ хв,-, г,] = 1, 2. (4)

Для построения адаптивных сеток в физической геометрии Бх2 вводится управляющая метрика, ковариантные (контравариантные) элементы которой обозначаются д-(д^) в параметрических координатах з1,^2 и (д^7) в сеточных координатах 2,

Наиболее общая формулировка управляющей метрики в физической геометрии Бх2 имеет следующий вид |1|:

д-3 = *(8)дХТ + К' фК», г,] = 1, 2, к = 1,...,1, (5)

где г (б) > 0 — весовая функция, Кгк (б) — компоненты некоторого ковариантного вектора ¥к (б).

Здесь и далее в выражениях, имеющих одинаковые индексы и по включающих знаки < + < — ^ и считается, что по этим индексам проведено суммирование.

2. Сеточные уравнения

Для построения сетки в Бх2 необходимо найти значения преобразования б(£) (2) в узлах эталонной сетки, заданной в вычислительной области Е2, Затем полученные значения нужно отобразить из области Б2 в Бх2 при помощи параметризации (1), Для нахождения дискретного преобразования б(£) численно решается задача Дирихле для обращенных уравнений диффузии |1|:

d2 s'

= г, ],1 = 1,2

С1 ' tot t

где

s(0 „ = m (6)

dS3-P ds3-q

ij _ j2 гj _ / 1 \i-\-j-\-p-\-QnPQ_____

a - j - \ i) gs

Pl = (7)

v y w(s) v ys/ 3-m w

i,j,M,m,p,5 = l,

J — якобиан преобразования s(£), w(s) > 0 — весовая функция (можно положить w(s) = 1), : dS2 — dS2 задает отображение границы вычислительной области S2 на границу параметрической области S2,

Если в формулах (2) положить w(s) = -/(f, то уравнения (6) становятся эквивалентными обращенным уравнениям Бельтрами (gs — определитель матрицы gj).

Для нахождения численного решения краевая задача (6) заменяется на нестационарную краевую задачу относительно функций s'(£,t), l = 1, 2:

= 941 _ Pi dt д£Ю£з

s'(*,t) = ^(0, £ e dS2, t> 0,

0) = s0(D, | e S2, (8)

^ j,M = 1,2,

где s0(£) — l-я компонента начального преобразования

so(£) : S2 - S2, so(|) = [s0(|),s2(|)],

задаваемого пользователем,

В случае, когда левая часть системы (6) является эллиптическим оператором, решение задачи (8) сходится к решению (6) при t — то [2].

3. Метод конечных элементов

Рассмотрим метод конечных элементов применительно к задаче (8), Согласно методу Фаэдо—Галеркина для нестационарных задач [3, 4] введем в вычислительной области с2 Ж-мерное подпространство пробных функций V\ где N — число внутренних узлов эталонной сетки в с2. Необходимо найти функцию м(£, £), принадлежащую при каждом £ > 0 подпространству Vй и удовлетворяющую при всех ун Е Vй уравнению

которое является дискретизацией слабой формы для уравнений (8), а решение п € Vн уравнения (9) есть приближение к решению в1, I = 1, 2, уравнений (8),

Чтобы записать полученную задачу в операторной форме, выберем в пространстве пробных функций Vй базис ...,

(£ ) = Г ^ вСЛИ к = p, ,

) | о, если к = р, ^ '

где <^р($к) — значение функции в к-м узле эталонной сетки = (£¿,£2). В данной работе в качестве базисных функций использовались финитные кусочно-линейные функции.

Тогда можно разложить функцию п(£,Ь) по базису

п(£,г) = пг + Пк <^к, к =1,...,Ж,

где пг — значения фун кции п(£, Ь) на границе вычислительной области дЕ2, Для задания граничных значений пг введем на дЕ2 дополнительный набор функций фд аналогично (10), В этом случае имеем

пг = ПТфд, д = N +1,...,Жг, (11)

где Nг — количество внутренних и граничных узлов сетки.

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

^ I <рк<рт<1£ = -ик I

з2 з2 з2

к,т = 1,..., N, (12)

или в матричном виде

ди

М— = Ки-¥, (13)

дЬ к '

где М — матрица масс, К — матрица жесткости, Е — вектор нагрузки определяются по формулам

Мтк = Рк <Рт<1 £,

Кт = / Рфтй(14)

2

2

С учетом граничных условий уравнение (13) записывается в следующем виде:

М

V

Е

( дщ \

т

дим

~дГ

/

и

N+1

\ им

\

К

(

Кг

Е

/

и1

UN

\

и

N+1

\ ^ /

где Кг = К при условии замены ^ на фк при к > N (11).

( Р \

N

0

V /

(15)

4. Численный алгоритм

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

К

суммы К = К11 + К12 + К22, где К11, К12 и К22 соответствуют слагаемым исходного

<92 <92 >

уравнения (9), содержащим производные

К12

К11

_д

дс

д^1' 2 5£25£2

<эе

, А именно:

/ 11 \ 0^рт

5/12 \дЧ>тА£

„ 22 /5/22

5е

5е2

(16)

Тогда для системы уравнений (15) можно записать аналог схемы стабилизирующей поправки [5]:

ип+1/2 _ ЛЛп

М-

— = К11ип+1'2 + К12ип + К22ип - Е,

т

М

ип+1 ип+1/2

К22ип+^ К22ип,

т

или в более компактном виде

ь1ип+1/2 = д1, ь2ип+1 = д2,

где

Ь1 = ^М-К1\ Ц1= (^М + К12 + К22УП-¥, 1} = \М- к22, д2 = -Мип+1/2 - к22ип,

тт

11п = I ип ип ип ип

и — | и1 , . . . , иN, UN+1, . . . , иNг

(17)

(18) (19)

г

На каждом шаге по времени системы линейных уравнений (18), (19) решаются при помощи метода сопряженных градиентов. Каждая система решается дважды дня

u

(в1)" = (s%..., (sV

l

1, 2,

где

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

1

max — 1<fc<N т

(s1)" - (s1)n+1 + (s2)n - (s2)n+1

< e

(20)

для некоторого достаточно малого е.

На каждом шаге по времени необходимо заново собирать матрицы Ь1, Ь2 и векторы Q1, Q2, причем Ь1, Ь2, Q1 можно собирать одновременно, а в Q2 необходимо добавить

слагаемое — Ми'1+1^2 уже после решения системы (18),

т

5. Вычислительные аспекты

5.1. Триангуляция вычислительной области

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

Итак, пусть в вычислительной области S2 задана триангуляция т1 либо т2 (см, рис, 2), Ячейка сетки T1 при дискретизации т1 представляет собой равнобедренный прямоугольный треугольник с длиной катета h, катеты треугольника направлены вдоль координатных линий £1 = const, £2 = const. Ячей ка T2 для дискретиз ации т2 является равно-

h

натному направлению £2 = const,

Рис. 2. Триангуляция Т1 (слева) и Т2 (справа)

5.2. Локальные матрицы масс, жесткости и вектор нагрузок

т1

Для каждой ячейки Т1 триангуляции т1 необходимо задать матрицы М1ос, К1ос и вектор Е1ос, Аналогично индексам т и к глобальных матриц в (14) обозначим индексы локальных матриц через то = (0,1, 2) и ко = (0,1, 2). Таким образом, то, к0 — локальные номера вершин треугольника Т1 (см, рис, 2), а также номера строк и столбцов локальных матриц.

Согласно формулам (10), (14), (16) и учитывая, что

К^ос _ К 111°с + К 121°с + К221ос

т1

жесткости:

211

1ос то

м1ос =

24

121 1 1 2

К

111ос

К

121ос

К

221ос

а22(1ко)

-1 1 0

1 -1 0

0 0 0

-2 1 1

1 0 -1

1 -1 0

-1 0 1

0 0 0

1 0 1

Здесь предполагается, что а4 (^ко) — значен не а4 [в(£)] из (2) в узле эталонной сетки

ко к

Р г(£то) — значен не Р г[в(£)] в узл е т.

5.3. Локальные матрицы масс, жесткости и вектор нагрузок

т2

Аналогично случаю ячейки Т1 для каждой ячейки Т2 триангуля ции т2 (см, рис, 2) по формулам (10), (14), (16) и учитывая, что

К

1ос

К

111ос

К

121ос

К

221ос

задаются вектор нагрузок и локальные матрицы масс и жесткости:

^1ос = 4-/г2Рг(|то), М1ос

то

12

48

2 / 2 11 1 2 1 1 1 2

К

111ос

л/3 и —« (€к0)

-1 1 0 1 -1 0 0 0 0

1

2

г.-12/ос N

А =

-1 0 1 0 1 -1 1 -1 0

11

К

221ос

1 -1 -1

Здесь предполагается, что т0, ко _ локальные номера вершин треугольника Т2, аг(^ко) — значение аг^(£)] го (2) в узле с локальным номером к0, который соответствует глобальному номеру к. Аналогично Р1($то) — значение Р1[б(^)] в узле т.

Пример вычисления локальной матрицы жесткости для триангуляции т2. Согласно формулам (14) и определению базисных функций (10) для задания КО^ необходимо вычислить величины

то

д

У д£г

Т2

— т-о, к0 = 0,1, 2, /../ 1.2.

Дня упрощения дальнейшего вывода введем обозначения

* = У = Ц\ Ак0 = а^[ с1 = ^к

и предположим, что ячейка Т2 расположена одной вершиной в начале координат, как показано на рис. 3.

В новых обозначениях распишем вычисление следующих интегралов:

/1

д_

дх

Ако (х, у) СхСу, /2

Т2

д_

дУ

Ако(х,у) йхйу

с1у

Ь,/2

дА

ко

дх

Сх +

Т2

ь-у/у/3

I

Н/2

дА

ко

дх

Сх

Ако{к/2, у) - А^у/у/З, у) + Ако{Н - у/у/З, у) - Ако{Н/2, у)] Су

По формуле трапеций

/1 « -С 2

-Ако(Ь/2, С) + Ако(Ь/2, С) - Ако (0,0) + Ако (Ь, 0)

у/Ък

Ако(Ь, 0) - Ако(0,0)

к х

а

1

о

а

о

4

0

Рис. 3. Базисная ячейка Т2

Следовательно,

л/Зк

4

/1 « <

Ь/2 у/Зх

а4(£ко), если ко = 0,

(£ко), если к0 = 1,

если к0 = 2, н -\/зж+2 а

(21)

/2 = /

4 0,

дА

ко

дУ

+ / ^ж

дА

ко

/1/2

Л/2

дУ

Аж(ж, л/зж) - Ад,о(ж,0)

^ж +

Ако(ж, -л/Зж + л/3/г) - Ако(ж, 0) ¿у.

//2

По формуле трапеций

г к

-/ 9 ~

Считаем, что Ако{к/2,0) « - Ако(0, 0) + Ако(М)

Ако(к/2,^) - Ако(к/2,0) .

Следовательно,

к ■■

—¿а%3(£ко), если /с0 = 0,1, к4

2аЧ(^о)> если ко = 2-

(22)

Используя формулы (21), (22) и определение базисных функций (10) можно получить значения К»^ для триангуляции т2, Локальная матрица жесткости для т1 вычисляется аналогично на треугольнике вида Т1,

5.4. Вычисление производных в коэффициентах уравнений

Значения коэффициентов уравнений а4 [б(^)] и Р1[б(^)] в (2) зависят от производных д/д^1, д/д£2, Поскольку эти значения требуются только во внутренних узлах сетки, они могут быть вычислены при помощи центральных разностей с использованием значений с предыдущего временного слоя. При этом, следуя логике метода конечных элементов, для вычисления производных в узле к необходимо хранить номера соседних узлов Рь т2, р2 (см, рис, 2), т1

правлений £ \ £2, а для триангуляции т2 производные по направлению £2 будем вычислять, используя формулу производной по направлению (см, рис, 2):

д/ д/ д/ . ~д1=д1тс°8а+д^ со8(тг/2 - а).

В данном случае а = п/3. Следовательно,

<9£2 у/з\ д1 д£М

д/ _, /Р2 /т2 д/ _, У»

I »1

д/

2к

д£1

2к

о

о

о

/

то

6. Примеры расчетов

На рис, 4-6 представлены сбалансированные адаптивные сетки, построенные с помощью описанного метода с использованием управляющей метрики для адаптации к значениям некоторой функции [1]

gj = F(s)<j, i,j = 1, 2, (23)

где gj _ контравариантные компоненты управляющей метрики, позволяющей полу-

F(s)

иия. Данная метрика является частным случаем формулы (5), Под сбалансированной понимается сетка, удовлетворяющая двум и более требованиям адаптации [1].

Численные эксперименты показали, что выбор шага по времени т = 0.01 в методе (17) для данных примеров является оптимальным, поскольку его увеличение приво-

т

чивает число итераций, необходимых для сходимости метода. Итерационный процесс продолжается до тех пор, пока левая часть в условии сходимости (20) не достигнет величины порядка 10-5,

На рис, 4, о показана сетка, сгущающаяся в окрестности трех точек с координатами (xk, yk) k = 1, 2, 3, Управляющая функция в формуле (23) вычислялась следующим образом:

F(s) = 2 min рк(s),

k=1,2,3

где pfc(s) = \J(s1 — Xk)2 + (s2 — ук)2 — расстояние до точки (Хк,Ук)■

Форма вычислительной области S2 совпадает с формой S2, эталонная сетка в S2 задается в виде триангуляции ti Количество итераций 500,

1ШГЛ1

улт

'f'A'A'

Шщш

ШI

шщШ

ми—

-..шА'ЛГЛ'Л_____.

WßBBSXSK

-Sssfu.-Л^вюЩв

ШШ8Ш

Ш

шшш

яявяШ

тташШШ

| WMTAViWi'tVi

i 'лшш'Л'Л'Л

П&Ш*/ШШММЖ1№ММММММ№

а

Рис. 4. Пример адаптивной сетки с триангуляцией Ti (а) и Т2 (б)

На рис, 4, б представлена адаптивная сетка в области с границами

в1 = г, в2 = 0.05яп(3пг), в1 = (1.5 - г)г3, в2 = 0.7г, в1 = 1 - 0.5г2, в2 = 0.7г, 0 < г < 1.

Управляющая метрика (23) задавалась с использованием функции

0.4 2

F (s)

0.4 + f i (s) + f2(s) + /a(s)

где

ffc (s) = exp ((s)2/0.01

k = 1, 2, 3;

^1(s) = s2 - 0.05sin(3ns1),

<^2(s) = Vis1 - 0.5)2 + (s2 - 0.5)2 - 0.001,

^s(s) = s2 - 0.25 - 0.7(s1 - 0.5)2.

Такое задание управляющей метрики позволяет получить сетку, сгущающуюся в окрестности кривых ^(s), <^2(s), <^3(s). Для сгущения сетки на границе области использовался алгоритм, основанный на решении одномерных обращенных уравнений Бельтрами [1]. Вычислительная область S2 имеет форму равностороннего треугольника, в котором задана триангуляция т2. Количество итераций 1000,

На рис, 5 изображены адаптивные сетки, сгущающиеся вблизи разрезов (слева), и соответствующие вычислительные области (справа), Форма каждого разреза представляет собой аэродинамический профиль крыла NACA,

i- S * с i? ^ w

ó í» N % ^ ^ tf1 ^ •• ■ ■'■

.< o N N ^i?^¿ч

' í * i?«?* с ^ í', ¡i •:

J J y í 1 I WWlw

¿' >■'' ;'''' '' i" - . ik>\\<' >

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

Разрезы на рне, 5, о, слева были заданы аналитически по формуле

5 = у0 ± 5£с

0.2969(х/с) - 0.126(х/с) - 0.3537(х/с)2 + 0.2843(х/с)3 - 0.1015(х/с)4

где х = 51—х0, (х0, у0) — начальная точка разреза, с — длина разреза, £ — максимальная толщина разреза относительно его длины (см, http://www.aerospaceweb.org). Управляющая метрика (23) задавалась при помощи функции

0.7

т \о.7 + ш + ш + ш, Д (в) = ехр((в)2/0.01) , к =1, 2, 3,

где ^(в) — функция, задающая линию к-го разреза. Количество итераций 200,

У сетки, представленной на рис, 5, б, слева, разрезы имеют форму профилей крыла МП 43, МП 60, МП 91, заданных дискретно по данным с сайта http://www.mh-aerotools.de/airfoils/. Пусть N — общее число точек, задающих разрезы (х&, к = 1,..., N. — их координаты. Тогда управляющую метрику (23) для построения сетки на рис. 5, б, слева можно задать с использованием следующих функций:

Р (в)

0.7

,0.7 + / (в),

/(в) = ехр( -^(в)2/°.01), ^(в) = тт (в)

V / К=1,...,Мс

где рк{в) = а/(в1 — Хк)2 + (в2 — Цк)2 — расстояние до к-й точки разрезов. Количество 400

На рис. 6 приведены сетки, построенные на поверхности сферы. Поверхность получена склейкой двух полусфер, каждая из которых задавалась параметризацией круга радиусом К = 1/3 с центром в точке хо = 0.5, уо = л/З/6. Вычислительная область

Рис. 6. Примеры сеток, построенных на сфере, с единичной метрикой (а) и с адаптацией к заданным кривым (б)

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

х1 = х0 + «(в1 - х0), Х2 = У0 + а(в2 - У0), х3 = ±(а - 1) Л,,

где

к = к2 + (к1-к2)г/Н, г=Л/(^-хо)2 + (з2-уо)2,

к2 + ^/к4 - (к2 + г2)(к2 - К2)

а =---.

к2 + г2

Сфера па рис, 6, а задавалась для = Л,2 = Д, на рис. 6, б для = 0.5Д, Л,2 = 1.5Д,

Для построения сетки, представленной на рис, 6, о, использовалась единичная метрика, т. е, в формуле (23) Р(в) = 1, В этом случае обращенные уравнения диффузии (6)

100

Сетка на рис, 6, б построена при

п.) ■ °'4

,0.4 + /1 (s) + /2(s)

где

= к= 1,2,

0.005 0.005 +

<^(s) = s2 - 0.2sin(3ns1) - 0.3, p2(s) = s2 - 0.2 sin(3ns1 + n) - 0.3.

300

одномерного алгоритма [1].

Авторы выражают благодарность д-pv физ.-мат. наук В.П. Ильину за консультации по вопросам метода конечных элементов.

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

[1] Технология построения разностных сеток / В.Д. Лисейкин, Ю.И. Шокин, И.А. Васева, Ю.В. Лиханова. Новосибирск: Наука, 2009.

[2] Калиткин H.H. Численные методы. М.: Наука, 1978.

[3] Стренг Г., Фикс Дж. Теория метода конечных элементов. М.: Мир, 1977.

[4] Зенкевич О., Морган А. Конечные элементы и аппроксимация. М.: Мир, 1986.

[5] Яненко H.H. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967.

Поступила в редакцию 7 сентября 2010 г., с доработки — 23 мая 2011 г.