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

МОДЕЛИРОВАНИЕ ПРОЦЕССОВ

УДК 516.6.531.33

В. В. Кузенов

ПОСТРОЕНИЕ РЕГУЛЯРНЫХ АДАПТИВНЫХ СЕТОК В ПРОСТРАНСТВЕННЫХ ОБЛАСТЯХ С КРИВОЛИНЕЙНЫМИ ГРАНИЦАМИ

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

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

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

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

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

Введем в декартовой системе координат ХУ2 прямоугольный параллелепипед ЛЕГЕБСОИ, для которого непрерывно дифференцируемое отображение в криволинейный параллелепипед (гексаэдр) Л'Е'Г'Е'Б'С'О'И'. При этом прямоугольная сетка, нанесенная на область ЛЕГЕБСОИ, образует гладкую криволинейную сетку в области Л'ЕТ'Е'Б'С'О'И'.

Обозначим через г радиус-вектор в системе координат ХУ2 и введем вектор и = г * — г, характеризующий смещение точек. Здесь г иг * радиусы-векторы точек областей до (г € ЛЕСБ) и после (г* € Л'Е'Г'Е'Б'С'О'И'О) преобразования. Тогда уравнения, определяющие смещения их,иу,иг и описывающие продольную деформацию пластин [1], в декартовой системе координат ХУ2 имеют вид:

ЛД (Лхдих) + (1—£) ± (е, д-Цх) + (1—£) 1 (с,дих) +

дх \ дх ) 2 ду \ ду ) 2 дг \ дг )

+ 2 { дх (Бх ду ) + дх (Ох дг ) } 0

Л д_(АЩ,) + (1 — а) д_ ( диу_) + (1 — а) д_ ( диу_) +

Лу ду\Луду)+ 2 дх\Еу дх) + 2 дг\Су дг ) +

+ ^ {£) + ( (о_£)} = о;

Л.I (л- ^) + ^т (Е.т1) + ^т (с- ^) +

дг \ дг ) 2 дх \ дх ) 2 ду \ ду )

+ ^ {I (Б- £) +1 (О-^)} = 0.

(1)

Здесь принято, что все коэффициенты, входящие в данную систему

уравнений, равны между собой: Ai = Bi = Ci = Di = Gi = F = = (V f, V/)a/2 + e, i G {x, y, z}. Граничные условия, необходимые для решения данной системы уравнений, задаются следующим образом:

Ui\r = ri \d(A'B'F'E'D'C'G'H') — ri\d(ABFEDCGH), i G {x,y, z}, где символ

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

При этом / является функцией, управляющей адаптацией расчетных узлов, коэффициент а задает необходимую степень сгущения узлов расчетной сетки. Коэффициент а G [-1,1] характеризует отношение поперечной деформации к продольной. Коэффициент Xi > 0 увеличивает "диффузию" сеточных линий от границы внутрь расчетной области A'B'F'E'D'C'G'H', обеспечивая нахождение узлов сетки внутри расчетной зоны. Параметр e ~ min (h) выбирается не

hi £ A B' C'D'

равным нулю (e = 0) и порядка шага сетки h, чтобы избежать особенности в тех узлах, где V/ = 0. Эта система уравнений является обобщением метода построения регулярных сеток [1] в случае переменных коэффициентов A,t, B,t ,Ci,Di,Gi, i G {x,y,z}, входящих под производные и трехмерную систему координат.

В результате решения системы уравнений (1) можно построить неравномерную сетку в криволинейном параллелепипеде (гексаэдре) А'Б'F'E'D'C'G'H', в котором поверхности, определяющие форму границы деформированной области, задаются гладкими функциями в системе координат XYZ. При этом расположение узлов расчетной сетки на граничных поверхностях задается до начала решения задачи. Введем в системе координат XYZ равномерную прямоугольную сетку с узлами Xi, yj, zm(i = 1,... ,L; j = 1,..., M; m = 1,..., N) и постоянными пространственными шагами {hx,hy,hz}, которые входят в состав управляющих параметров и выбираются из минимальных значений шагов сетки вдоль граничных поверхностей. Зададим конечно-разностное представление производных на введенной сетке с помощью следующих формул [2]:

д_

дх

dg

дх J i,j,m

ai+1/2,j,m (gi+1,j,m gi,j,m) ai-1/2,j,m (gi,j,m gi—l,j,m)

hX

_ ai,j,m + ai±1,j,m d

Ч±1/2,3,т = 2 • dX

+ O (h

2) •

x

a—

. dy

ij,m

ai+1,j,m (gi+1,j+1,m gi+1,j-1,m) ai-1,j,m (gi-1,j+1,m gi-1,j-1,m)

4hx hy

+ O (hX, hy).

f

Теперь можно привести конечно-разностное операторное представление исходной задачи Л11% 3 = 0:

, . . -г+1/2,з,т (их,г+1,з,гп — — -г—1/2,з,гп (их,г^.т — их,г—1^,т)

)х = Лх — |-

пх

(1 — а) \ + 1/2(их,-1,з + 1,гп — их,г,з,т) — —г,]-1/2,т {их,г,],т — их,г,з~1,т) ]

+ ^ I ч Г

(1 — а) ( -г,з,т+1/2 (их,г,з,т + 1 — ^х,г,^,гш) — -г,д,-т — 1/2 ("^х,г,^,гш — их,1,д,т — 1) Л

~Ц / +

(1 + а) Г —г+1,0,т (иу,г+1,з + 1,т и'у,г + 1,у — 1,т) —г—1,з,т (иу,г—1,у + 1,т и'у,г—1,^ — 1,т)

I \ ' / J - I у ¿М- I I _I -■- > J /_- -■- ? " ^ \ Л)" -■- I ^ ; " у_»)' -■- > J -■- ? " у / I I

2 " '

+ (1 + a) f Fi+l,j,m (Uz,i+l ,j,m+l Uz,i+l,j,m — l) Fi—l,j,m (Uz,i— l,j,m + l Uz,i-l,j,m-l)

2 l 4hxhz

/ ЛТТ \ _ \ Fi,j+l/2,m (Uy,i,j+l,m — Uy,i,j,m) — Fi,j — l/2,m (Uy,i,j,m — Uy,i,j — l,m)

(AUi,j,m )y — Лу — T

hy

(1 — a) ( Fi+l/2,j,m (Uy,i + l,j,m — Uy,i,j,m) — Fi—l/2,j,m (Uy,i,j,m — Uy,i—l,j,m) 1

+ hx ; +

(1 — a) f Fi,j,m + l/2 (Uy,i,j,m+l — Uy,i,j,m) — Fi,j,m — l/2 (Uy,i,j,m — Uy,i,j,m — l) 1

+ _4 / +

+ (1 + a) f Fi,j + l,m (Ux,i+l,j+l,m Ux,i—l,j + l,m ) Fi,j — l,m (Ux,i+l ,j — l ,m Ux,i — l,j — l,m ) ^ +

2 l 4hyhx

+ (1 + a) f Fi,j + l,m (Uz,i,j+l,m + l Uz,i,j+l,m — l ) Fi,j — l,m (Uz,i,j — l,m + l Uz,i,j — l,m — l)

2 \ 4hyhz

/ ЛТТ \ Л Fi,j,m + l/2 (Uz,i,j,m+l — Uz,i,j,m) — Fi,j,m — l/2 (Uz,i,j,m — Uz,i,j,m — l)

(AUi,j,m )z — h22 +

+ (1 — a) ( Fi+l/2,j,m (Uz,i + l,j,m — Uz,i,j,m) — Fi—l/2,j,m (Uz,i,j.m — Uz,i—l ,j,m) "I +

+ l 4 J+

+ (1 — a) \ Fi,j + l/2,m (Uz,i,j+l,m — Uz,i,j,m) — Fi,j — l/2,m (Uz,i,j,m — Uz,i,j — l,m) 1 +

+ 1 4 ) +

(1 + a) / Fi,j,m + l (Ux,i+l,j,m + l Ux,i—l,j,m+l) Fi,j,m — l (Ux,i+l,j,m — l Ux,i—l,j,m — l~) .

2 \ 4hzhx

+ (1 + a) / Fi,j,m + l (Uy,i,j + l,m + l Uy,i,j — l,m + l ) Fi,j,m — l (Uy,i,j+l,m — l Uy,i,j — l,m — l)

2

где г = 2,...,Ь — 1; ] = 2,...,М — 1; т = 2,. — 1.

Для решения задачи Ли = 0 используем метод установления [3]. При этом шаг по "времени" т найдем с помощью итерационного метода вариационного типа. Для этого определим вектор невязок

К = (Лх,Лу ): к к

^х,г,],гп = (Ли%,],т )х —

= ,Ьк • ,Ьк • • | ^0,0, ^ ,

%,3,т I х,г,],по "у , %,],по ,г,],т) )

г = 2,...,Ь — 1; 2 = 2, ...,М — 1; т = 2,...,И — 1, где к — номер итерации.

Введем скалярное произведение следующим образом [3]:

Ь-1 М-1 N-1

i=2 з=2 т=2

Ь-1 М-1 N-1

+ ^ ^ ^ ^ ^ ^ aytitjt'mbytitjtmhxhy+

г=2 j=2 m=2

Ь-1 M-1 N-1

+ ^ ^ ^ ^ ^ ^ az,i,j,mbz>i>j>mhxhy .

г=2 j=2 m=2

Минимизируем значение невязки Я^^^ = Аи.^^^ — используя модифицированный вариант итерационного метода вариационного типа — метод минимальных невязок [4]. В этом случае итерации следует проводить по формулам [4]: = Щз^ + стЯ^^, с « 0,9,

т = .

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

В настоящей работе для построения регулярных адаптивных сеток, близких к ортогональным, вместо системы уравнений (1) использовали модифицированную в соответствии с работой [5] систему уравнений:

Ах ± (^) + £ (^) + (Ц£) £ (^) +

дх \ дх ) 2 ду \ ду ) 2 дх \ дх )

+ IТ (Р1Т) + Т (Р1Г » +

2 [ дх \ ду ) дх \ дх

+ к jd( 1922 933 + JL f /вззди ЗиЛ + d_ f Ig iig22 9UX opt ' dx \У gii dx J dy \y 922 dy J dz \y 933 dz

9 (vdUy\ + (1 -v)d (FdUy) + (!-v)Л (FdUy ) +

dy \ dy J 2 dx \ dx J 2 dz \ dz

+ 2 { dy ( dx ) + dy ( dz ) } +

+ kopt Ü (^'-U-) +1 (^^) +1 (2)

dx yJ 9ii dx ) dy \y 922 dy J dz \V 933 dz

X1 F

z

cUz z

+ (1 - v) ± (FdU±

2 dx \ dx

+ (1 - v) ± fF9U±

2 dy \ dy

+

(1+ v)

d fD dUx dz V z dx

+ SL (GzdUy

dz \ z dy

+ +

2

+ k í — ( 1922933 диЛ + д( 1933.gn BU А + д( 1911922 BUA \ =0 op{ dx \\¡ gii dx ) dy \\¡ 922 dy J dz \\¡ 933 Bz )) '

где коэффициент krjpt регулирует интенсивность ортогональности сеточных линий.

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

3

9ik =

^ да''- dqk

a=l

3

к? г ? I 1, г — 1 1 2 3

^ 9гк9 ? — 0? — | 0 г — 1 ' д — x, д — У, д — ^ к=1 ^ а контравариантные компоненты метрического тензора находятся с помощью формулы

9гк — Агк /9,

где 9 — фундаментальный определитель det \\(9^к ||; Агк — алгебраическое дополнение элемента 9^к в этом определителе.

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

Два первых примера показывают возможности генерации адаптивных сеток внутри расчетной области вследствие регулирования расположения точек сетки на границе заданной области с помощью одномерных численных [6] и аналитических алгебраических методов [7].

В случае численной перестройки сетки вдоль одной из координатных линий или вдоль криволинейной границы расчетной области (или ее части) применяется принцип равномерного распределения весовой функции ы. В этой численной методике входными параметрами, заданными пользователем, являются максимальный (Ахтах) и минимальный (Ахтщ) шаги сетки и некоторая положительная (дополнительно монотонизированная, т.е. имеющая один минимум и максимум) управляющая функция f. Весовая функция ы отличается от приведенной в работе [6] и задается следующим обра-[1 + АГ]1'В , В — 0 — (/_/.)/(/ _

ч О п ' зтш// \Jmax зтт)->

зом:ы — • 1, В — 0

А — (Ахтах/Ахт1п)В — 1. Для определения шагов сетки Ах^ мо-

^ Г ^х

гут служить уравнения — ы—

= 0, x(0) = 0, x(1) = L, £ е [0,1],

/( N Л

или Ах^ — Ь — , где Ь — длина границы области в фи-

/V тч

зическом пространстве. В случае А — 0 узлы сетки расположены

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

равномерно. Постоянную Б находят из условия равенства минимального расчетного шага сетки min Axi заданному значению Axmin.

i

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

Пример 1. В рассматриваемом примере демонстрируются возможности предложенного метода, в котором в качестве управляющей функции / используют аналитическое алгебраическое преобразование. Это преобразование производит растяжение физического пространства в определенных, предположительно заданных областях: в пограничном слое, слоях смешения, на фронте ударных волн и в зоне контактных разрывов. В данной ситуации сетка построена с помощью одномерных численных [6] и аналитических алгебраических методов [7] в зоне критической части сопла, адаптирована к пограничному слою и представлена на рис. 1, а-г.

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

Для случая алгебраического метода адаптации управляющая функция имеет вид [8]:

f (х) = ^- + (1 - р) (: - ШМ) •

П* = (п - па)/(пв - па), П g [пв, па], ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. "Машиностроение". 2008. № 1 9

Рис. 2. Регулярная расчетная сетка, построенная с применением объемной управляющей функции f (х, у) = ехр ^—А (х — у)2^ с помощью разработанной численной методики (а) и численной адаптации с использованием ортогонализи-рующей процедуры (б)

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

Пример 2. В этом примере в прямоугольной области физического пространства расчетная сетка адаптируется к функции вида

/ (х,у) = ехр (—А (х — у)2) , А = 20.

Этот пример, как отмечается в работе [6], моделирует достаточно сложную ситуацию для адаптации регулярной четырехугольной сетки: управляющая функция / имеет градиенты, образующие с границей расчетной области угол, близкий к 45°, что может приводить к сильной скошенности адаптированной сетки.

Результаты расчетов, представленные на рис. 2, а, основаны на системе уравнений (1).

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

Работа выполнена при частичной поддержке программ фундаментальных исследований РАН и грантов РФФИ №05-01-00780 и №07-01-00133.

СПИСОК ЛИТЕРАТУРЫ

1. Игнатьев А. А. Построение регулярных сеток с помощью механической аналогии // Математическое моделирование. - 2000. - Т. 12, № 2. - С. 101-105.

2. Р о у ч П. Вычислительная гидродинамика. - М.: Мир, 1980. - 616 с.

3. Самарский А. А.,Николаев Е. С. Методы решения сеточных уравнений. - М.: Наука, 1978. - 590 с.

4. Альшина Е. А., Болтнев А. А., К а ч е р О. А. Эмпирическое улучшение простейших градиентных методов // Математическое моделирование. -2005. - Т. 17, № 6. -C. 43-57.

5. L i s e i k i n V. D. Grid generation methods. - Berlin: Springer, 1999.

6. Накахаси К., Дейуэрт Дж. C. Автоматический метод построения адаптирующихся сеток и его применение в задачах обтекания профиля // Аэрокосмическая техника. - 1987. - № 12. - C. 10-18.

7. С у р ж и к о в С. Т. Аналитические методы построения конечно-разностных сеток для расчета аэротермодинамики спускаемых аппаратов // Вестник МГТУ им. Н.Э. Баумана. Сер. "Машиностроение". - 2004. - № 2. - C. 24-50.

8. R o b e r t s G. O. Lecture Notes in Physics. Vol. 8. - Berlin, Heidelberg: Springer, 1971. - P. 171-176.

Статья поступила в редакцию 10.03.2007 Виктор Витальевич Кузенов окончил МГУ им. М.В. Ломоносова в 1983 г. Канд. техн. наук, доцент, ст. научн. сотр. лаборатории "Радиационная газодинамика" ИПМех РАН, доцент кафедры "Теплофизика" МГТУ им. Н.Э. Баумана. Автор более 50 научных работ в области теплофизики, радиационной газовой динамики и плазмодина-мики.

V.V. Kuzenov graduated from the Lomonosov Moscow State University in 1983. Ph. D. (Eng.), senior researcher of "Radiation Gas-dynamics" laboratory of Institute for Problems in Mechanics of RAS, assoc. professor of "Thermal Physics" department of the Bauman Moscow State Technical University. Author of more than 50 publications in the field of thermal physics, radiation gas- and plasma-dynamics.

ЖУРНАЛ "ВЕСТНИК МОСКОВСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА имени Н.Э. БАУМАНА"

В журнале публикуются наиболее значимые результаты фундаментальных и прикладных исследований и совместных разработок, выполненных в МГТУ им. Н.Э. Баумана и других научных и промышленных организациях.

Журнал "Вестник МГТУ им. Н.Э. Баумана" в соответствии с постановлением Высшей аттестационной комиссии Федерального агентства по образованию Российской Федерации включен в перечень периодических и научно-технических изданий, в которых рекомендуется публикация основных результатов диссертаций на соискание ученой степени доктора наук.

Подписку на журнал "Вестник МГТУ им. Н.Э. Баумана" можно оформить через агентство "Роспечать".

Подписывайтесь и публикуйтесь!

Подписка по каталогу "Газеты, журналы" агентства "Роспечать"

Индекс Наименование серии Объем выпуска Подписная цена (руб.)

Полугодие 3 мес. 6 мес.

72781 "Машиностроение" 2 250 500

72783 "Приборостроение" 2 250 500

79982 "Естественные науки" 2 250 500

Адрес редакции журнала "Вестник МГТУ им. Н.Э. Баумана": 105005, Москва, 2-я Бауманская ул., д. 5.

Тел.: (495) 263-62-60; 263-60-45. Факс: (495) 261-45-97. E-mail: [email protected]