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

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

Том 5, № 5, 2000

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

В. И. ПААСОНЕН Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: paas@net.ict.nsc.ru

A system of Navier-Stokes equations with respect to variable stream function and velocity vortex is considered. The aim of this work is to construct a compact scheme which is defined on a pattern with the dimensions 3 x 3 on a rectangular non-uniform grid. A third-order approximation scheme has been constructed. In the particular case of a uniform grid this scheme has the fourth order of approximation. The stability of the linearized scheme has been investigated and stability criterion has been obtained.

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

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

С другой стороны, при использовании схем повышенной точности, в том числе компактных схем, обычно наблюдается улучшение результатов расчетов (см., например, [1, 3-5]), даже при крутых профилях решения. Однако на равномерных сетках в области больших градиентов и для схем повышенного порядка аппроксимации также возможны значительные локальные погрешности [?, ?].

© В. И. Паасонен, 2000.

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

1. Постановка задачи

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

и = фу, V = —фх, Ш = иу — Ух, известным путем исключают давление и приходят к ф — ш — системе

+ (фу ш)х — (фхш)у = Мшхх + Шуу) + Р,

фхх + фуу = Ш. (1.1)

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

Пусть г — одна из независимых переменных х или у. Зафиксируем точку неравномерной сетки по переменной г и обозначим через к+ и к_ шаги сетки справа и слева от этой точки, а через Д+ и Д_ — разделенные разности “вперед” и “назад” первого порядка. Введем также обозначения для разности, суммы и произведения соседних шагов к+ и к_:

d = к+ — к_, в = к+ + к_, р = к+к_. (1.2)

Тогда для разностных аппроксимаций

Д= к_Д+ + к+Д_, л = Д+ — Д_ (1.3)

в ’ в/2 У 1

операторов одно- и двукратного дифференцирования по г справедливы разложения

Дп = п' + -и)'" + 0(к3),

6

Ли = и" + dw''' + Ъп'ш + 0(к3), Ъ = , (1.4)

3 12

где к — максимальный шаг.

Кроме того, ниже для выражения (дпх)г нам понадобится разностная аппроксимация

Л%(;) = 2 (д(г + к+) + д(г) Д+ — д(г> + д2(г — к_>Д^ „(г), (1.5)

1 Здесь рассматривается только дивергентная форма уравнения для вихря. Для недивергентной формы построение компактной схемы на неравномерной сетке проводится аналогично.

для которой справедливо разложение

И

Лд 'ш(г) = (д'ш'У +— (їди}"1 + Зд'и)" + 3д"и)') + 0(Л,2), 6

(1.6)

а также для выражения (дгих)г — аппроксимация вида

(1.7)

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

2. Аппроксимация уравнения Пуассона

Очевидно, что простейшая схема для второго уравнения системы (1.1)

имеет на неравномерной сетке всего лишь первый порядок аппроксимации. Следуя принципу построения компактных схем [?, ?], усовершенствуем ее, применив к слагаемым схемы одномерные осредняющие операторы Бх, Бу по ортогональным направлениям:

Здесь и в дальнейшем Е означает тождественный оператор, а знак “«” означает равенство с точностью до членов 0(к3). В качестве оператора Б осреднения правой части можно взять произведение БхБу или, например, выражение

Аппроксимация уравнения Пуассона схемой (2.1) с третьим порядком следует из равенств

д2 д2 ( д2 д2 \

Лх ~ Бхдх2, Лу ~ Буду2, БуЛх + БхЛу ~ Б (дХ2 + ду2) ,

справедливых ввиду специального выбора коэффициентов осредняющих операторов и в силу разложения (1.4). Нетрудно заметить, что при равномерной сетке все члены разложения нечетного порядка обращаются в нуль, и (2.1) превращается в классическую схему Микеладзе [?] четвертого порядка точности.

Установим знаки коэффициентов схемы (2.1). Для этого сначала запишем операторы Б* и Л* в индексах. Первый из них

(Лх + Лу )ф = и

(Бу Лх + БхЛу )ф = Би,

(2.1)

Бги = ах иі-і + а0и + а+иі+1

имеет коэффициенты

1 5

а± = (^2 + р ± 24г), а0 = - + . (2.2)

6^ Лі* 6 6р

Коэффициенты второго оператора

2 2

л± = тг~, л0 = -;г• (2.3)

З,г Л±.г р,г

Очевидно, что коэффициенты оператора БуЛж + £жЛу определяются через коэффициенты (2.2)-(2.3). Их выражения удобно представить, в соответствии с размером шаблона, в виде матрицы размерности 3 х 3 (ось х направлена вправо, а ось у — вверх):

а+ Лж + ах Л+ а+ Лж + аХЛ+ а+ Л+ + а+ Л+

Л^» I а^ Лу ау Лх I а^ Лу Л^ I ^аЛу

аУЛж + аж Лу аУЛХ + аХЛУ аУЛ+ + а+\ • (2.4)

Так как коэффициенты Л0 отрицательны, а а0 — положительны, то центральный элемент

матрицы (2.4) отрицателен. Сумма всех элементов матрицы равна нулю. Следовательно,

если потребовать неотрицательность всех остальных элементов (2.4), то для схемы (2.1)

будет справедлив принцип максимума [?].

Угловые элементы матрицы (2.4) в силу (1.2) и (2.2)-(2.3) с точностью до различных

~ 2

положительных множителей приводятся к виду

Рж + Ру + ^ж^ж Зж + 0у ^у Зу ,

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

|Л+ж — ^Хж! + !^+у — ЛХу1 < Л+жЛ-ж + Л+у Л-у • (2-5)

Заметим, что даже замена (2.5) более грубым условием

¡Л+г - Л—1< Л+гЛ—, £ = X, у (2.6)

дает вполне приемлемое ограничение на коэффициент изменения шага = Л-г/Л+г. При этом имеет место следующее ограничение на коэффициент :

^5 - 1 ^5+ 1

0.62 и ------ < аг < -------и 1.62.

2 - Уг - 2

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

2Интерес представляют только знаки элементов матрицы (2.4).

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

З; + Рж ± ^уЗу - Ру, З2 + Ру ± 4Зж - Рж.

Требование их неотрицательности в силу (1.2) эквивалентно симметричной системе неравенств

Л+жЛ-ж + |Л+ж - Л-ж! < Л+у + ЛХу + 3Л+уЛ-у,

Л+уЛ-у + |Л+у - Л-у| < Л+ж + Л-ж + 3Л+жЛ-ж, (2.7)

которая означает недопустимость резкого различия между средними местными шагами по различным переменным. В предельном случае равномерной сетки по обоим направлениям условия (2.7) совпадают с известным ограничением для схемы Микеладзе

— < 'т- < -5

V 5 Лу

Таким образом, условия устойчивости схемы (2.1) для уравнения Пуассона, рассматриваемого отдельно, а не в составе системы (1.1), задаются неравенствами (2.5), (2.7). Первое из них регулирует отношение соседних шагов сетки по каждой переменной, а второе ограничивает степень отклонения сетки от квадратной.

3. Аппроксимация уравнения для вихря

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

М^ЖЖ + ^УУ) = /> / = (фУ — (фж^0у + ^ — (3-1)

Тогда схема третьего порядка аппроксимации для него с точностью до множителя ^ совпадает с (2.1). Запишем ее в виде

МБУ ЛЖ + БжЛу)^ = ЬЛ (3-2)

причем оператор Ь ^ Б возьмем на этот раз в дифференциальной форме

„ -ж д -У д д2 -Ж-У д2 д2 Ь = Е + + Ьж— + ж У „ + Ьу

3 дх 3 ду ж дх2 9 дхду у 5у2’

Покажем, что правую часть (3.2) можно аппроксимировать с погрешностью О (Л3) разностными выражениями на нашем компактном шаблоне. Имеем

Ь/ и Б(^ - ^) + ^0 + ф1 + ^2, (3.3)

где — суммы членов О(Лг):

^° = (ифу )ж — (ифж)у, ф1 = -3- ^Ж + "3Г ^У >

02 = 6 0° + О0 + 6 0°

° 6ж°хх + 9 °жу + 6У °уу •

Представим члены нулевого порядка в виде суммы 0° ^ Л° + 0° главной части в разностной форме и погрешности второго порядка:

Я° = Дж(^Дуф) — Ду(иДжф),

0° "66” ^(иф”)УУУ (ифУУУ)^) ”66” ^(ифУ)жжж (ифжжж)у) • (3-4)

Затем выполним аналогичное разложение О1 ^ ^ + 01 в отношении слагаемых первого

порядка

Лі — "зр ГЛж(иДуф) — ДуЛжф^ —зу ( Лу(шДжф) — Джлуф

^2 / з \

°1 = "9 ( (^фж)ууу — (ифууу)ж — 2(^уФ”)”” ) —

/ з \

—9” I (иф”)жжж — (^фжжж)у— 2(^хфж)ж^ • (3-5)

Из равенств (3.4)-(3.5) следует, что 0°, 01 можно заменить разностными выражениями Я°, Л1, включив в остаточный член 02 дополнительные члены погрешности 0°, 02. Тогда равенство (3.3) преобразуется к виду

¿/ ~ Б(и — Р) + ^° + ^1 + ^2, (3.6)

где Я2 = 02 + 0° + 02 — сумма всех членов второго порядка малости.

Заметим, что если пренебречь слагаемым Д2 и подставить затем (3.6) в качестве правой части в (3.2), то уже получается схема второго порядка аппроксимации. Однако, как будет показано ниже, порядок можно повысить еще на единицу. Препятствием к этому служит наличие в выражении Д2 третьих производных по одноименным переменным от искомых функций, для непосредственной аппроксимации которых необходимо минимум четыре точки по одной координате. Для его преодоления необходимо представить Д2 в иной форме.

Для этого сначала, дифференцируя произведения, выделим из Д2 слагаемые с производными от функции ф порядков 3 + 1 и 1 + 3 (т. е. третьего по одной переменной и

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

р* ^ р* + ^ /р* . ^

а* = 7^ + 77* , 6* = ———, с* = - — + —

12 36 12 ’ V 6 9

при подобных членах равна нулю.

Затем все слагаемые разобьем на три группы Л — Р1 + Р2 + Р3. К первой и второй отнесем слагаемые с третьими производными (порядков 3 + 0 и 0 + 3) от функций ш и ф соответственно

Р1 ау фжшууу ажфу Р2 ажфжжжшу ау фууу

(3.7)

а к третьей — все остальные слагаемые

р3 — 3ау (фхушу )у 3аж (фжушж)ж +

+ ~ (фж----------- (фу^у )жу + —— ( (фу^)жжу — (фж^)жуу ) . (3-8)

+ Ьу (2фуу ^у + фу ^уу )ж — 6ж(2фжж^ж + фЖ^ЖЖ)у +

Г2 Г2 ^ ^

иж( / ) у ( / ) , и*цу

6 (фж^ж)жу — 6 (фу^у)жу + д

Слагаемые (3.8) не содержат тройного дифференцирования по одноименной переменной и могут быть непосредственно аппроксимированы на компактном шаблоне. В первых двух группах (3.7) третьи производные выражаются через производные низшего порядка по каждой переменной из однократно продифференцированных уравнений системы (1.1). При этом Р, удобно представить в виде суммы Р, = Р,0 + Р,1, выделив члены, содержащие (^ — Р). В результате получим

р0 = ^ (Ш, — р )у — ^ (Ш, — р )„,

1 Оуфж о \ «жфу о \

Р1 = “ (^ — ^^жж)у-------- (^ — ^уу )ж,

Р2 (ож Оу )^ж^у + Оу фжжу Ожфжуу ^у • (3'9)

Аппроксимируя после этих преобразований выражения (3.8)-(3.9) с использованием разностных операторов (1.3), (1.5), (1.7), получим с погрешностью третьего порядка

Р0 = о,ф д — Р),

р, = О,.^ / ЛД,„ ф — ЛД, + ^ Ду \ — / ЛД, * ^ — ЛД, „ ф — ^ Д \ ,

Р2 = («ж — Оу )Дж^Ду^ + Оу Дж^Л^Ду ф — «жДу ^Лу Джф,

Рз = 3оу Г-у Джф — 3ажГ^ж Ду ф+

+ЬуДж (2ЛуфДу^ -+ ДуфЛу^) — ЬжДу(2ЛжфДж^ + ДжфЛж^)+

Л2 Г2 Г Г

+у Ду (С ф) — ^ Дж(г-у ф) + (ЛжЛ-ф — Лу л-ф).

Теперь осталось заменить в (3.6) слагаемое Р2 суммой полученных выражений, а затем подставить результат в правую часть схемы (3.2). Объединяя при этом члены, содержащие (^ — Р), получим схему третьего порядка аппроксимации без дискретизации по времени:

В — Р) = ^(БуЛж + Лу — (Р0 + R 1 + Р12 + Р2 + Р3), (3.10)

на компактном шаблоне, где оператор

В = Б — ф Дж + «у^ Ду

с эквивалентной погрешностью может быть заменен факторизованным:

В = ВжВу = (Е + ^жДж + ЬЖЛЖ)(Р + ^у Ду + Ьу Лу),

Гж ОжДуф Гу Оу Джф

= —----------, 5у = — + ,

3 ^ 3 ^

обеспечивая расщепление разностной задачи на одномерные. Выбор конкретного способа дискретизации (3.10) по времени дает определенный вариант компактной схемы, не

обязательно двухслойной. На равномерной сетке погрешность схемы составляет величину O(h4). Очевидно, (3.10) включает в себя также и стационарный случай.

Заметим, что коэффициенты схемы при диссипативных членах имеют порядок O(h-2), а при конвективных — O(h-1). Поэтому в стационарном случае в предположении постоянства коэффициентов и при малых шагах сетки знаки коэффициентов схемы (3.10) определяются знаками коэффициентов схемы (2.1) для уравнения Пуассона. Следовательно, критерии устойчивости схемы (3.10) при достаточно малых h близки к ограничениям (2.5),

(2.7). Однако при конечных шагах или при коэффициенте вязкости ^, сравнимом с h, это утверждение становится неверным, т. е. принцип максимума нарушается, если в дополнение к чисто сеточным условиям (2.5), (2.7) не ограничить сверху также и разностное число Рейнольдса.

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

[1] ВАЛИУЛЛИН А.Н. Схемы повышенной точности для задач математической физики. Изд. НГУ, Новосибирск, 1973.

[2] КАЛИТКИН Н.Н. Численные методы. Наука, М., 1978.

[3] МИКЕЛАДЗЕ Ш.Е. О численном интегрировании уравнений эллиптического и параболического типов. Известия АН СССР. Сер. матем., 5, №1, 1941, 57-74.

[4] ПААСОНЕН В.И. Об одном методе построения высокоточных разностных схем и их применении в механике жидкости. Численные методы механики сплошной среды, 7, №6, 1976, 111-126.

[5] Толстых А.И. Компактные разностные схемы и их применение в задачах аэрогидродинамики. Наука, М., 1990.

Поступила в редакцию 21 февраля 2000 г.