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

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

Том 9, № 3, 2004

ПРИМЕНЕНИЕ НЕЯВНЫХ СХЕМ РУНГЕ - КУТТЫ ДЛЯ ЧИСЛЕННОГО ИССЛЕДОВАНИЯ ТЕЧЕНИЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ С КОНТАКТНЫМИ РАЗРЫВАМИ*

В. И. Пинчуков

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

e-mail: pinchvi@ict.nsc.ru

Implicit third order Runge — Kutta schemes are used for solving of two dimensional equations of incompressible fluid flows, written in the vorticity-stream function formulation. The results of numerical calculations are presented to illustrate the efficiency of the numerical method.

Введение

Если течения несжимаемой жидкости являются гладкими в начальный момент времени, они остаются таковыми и в последующие моменты времени. Однако структура их со временем часто усложняется настолько, что любой разностный метод на фиксированной дискретной сетке может катастрофически терять точность. Еще хуже ситуация с численным исследованием течений, содержащих контактные разрывы. Эти течения также могут усложняться, причем использование схем высоких порядков не спасает ситуацию, поскольку точность схем не проявляется на решениях недостаточной гладкости. В этой связи перспективен подход [1, 2] к моделированию контактных разрывов как линий уровня некоторой гладкой функции. В [3] соответствующая система уравнений, а также классическая система в переменных завихренность — функция тока решаются на основе ENO схем высоких порядков. В данной работе в рамках этого подхода анализируется эффективность неявных схем Рунге — Кутты третьего порядка [4-5] для моделирования течений несжимаемой жидкости.

1. Уравнения несжимаемой жидкости в произвольных криволинейных ортогональных координатах

Рассмотрим двумерную систему уравнений несжимаемой жидкости в переменных завихренность w — функция тока ф

Wt + (uw)x + (vw)y = V (Wxx + Wyy ),

* Работа выполнена при поддержке Американского фонда гражданских исследований и развития для независимых государств, образованных на территории бывшего Советского Союза (CRDF) (грант № RM1-2324-NO-02).

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

фхх + Фуу = --, и = фу, V = -фх.

При переходе к произвольным криволинейным ортогональным координатам п, х = х(£, п), У = п) эти уравнения переписываются в виде

+ (и-)« + = ^ {[(хП + уП)-е/1]« + + У2)-п/1]п }>

и = иуп - vyn, V = то« - иу«; (1)

[(4 + УП)^/I]« + [(х2 + у?2)^/I ]п = -I-; (2)

и = (ф«уп - )/7, v = - Ф«уп)/7, 1 = хпу« - УП• (3)

Представляет интерес также система уравнений для описания течений с вихревыми пеленами, которая может быть получена из этой системы в случае нулевой вязкости. Учитывая, что если при V = 0 /(£,п,£) есть решение уравнения (1), то и ш = ш(/) также есть решение этого уравнения, и предполагая, что завихренность локализована на одной лишь линии тока (линии контактного раздела), будем искать завихренность - в виде ш = с$(/ - /о), где / — некоторая гладкая функция, $(/) — дельта-функция, которая понимается как предел при е ^ 0 функции

¿(/) = тах(0,1 - /2/е2)1'5/(0.1875пе). (4)

В результате уравнение (1) заменяется на уравнение

/ + (и/)« + (V/)п = 0, (5)

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

2. Схема для уравнений несжимаемой жидкости

Введем обозначения и = ±и(...,хг ± Ахг, ...)^и(..., ,...), А±±г = /Ахг, ДХг = (#Х+ + )/(2Ахг).

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

ь

Не

3 т

_ - _ил« - т2АЯ2« + ^т4^« - К(т//*

X

х

3 Т

_ - _уЛп - т2АЯ2ч + ^т4^ - )П

(^ - <)- + Д = (П« + ппX; (6)

/г

гк

3 т

2 - ^иЛ« - Т2АЯ2? + ^Т4Я4? - К(т//гк)П«

х

х

3т

_ - 2УЛП - Т2АД2п + Ст4Д4п - к(т//гк)ПП

(^(2) „/^А +

- ^гк )3Т +

- ))3Т + (ДГ + Д(11))/2 = (П« + Пп)« + и^);

(7)

/гк [1 - 7Т2Д2« + ^Т4Д4« - (т//гк)П«] X

X [1 - 7т2Д2п + ^Т4Д4п - (т//гк)ПП] (

п! (^Г+( „,Г

+ (т7/гкД2« + ?П« + Т7/гк^п + ?ПП)(Ш* - ^)3/2 + (Д + ЯГ3)/4

(П« + Пп )« + 3^))/4;

(8)

Д( = Л«(/и^)гк + Лп(/уш)гк, Л« = Д0(1 - ¿+¿"/6), Лп = ДП(1 - ¿+¿"/6)

П п П

Цгк = (игкЛпУгк - ^ЛПХгк)//гк, ^ = (^Л«Хгк - МгкЛ«Угк)//гк;

«гк = [(Л«^) (ЛПУгк) - (ЛП^к)(Л«Угк)]//гк, ^ = - [(Лп^гк)(Л«Хк) - (Л«^)(ЛП^)]//гк,

/гк = (ЛПУгк )(Л«Хгк ) - (Л«Ук )(ЛП^ ).

(9)

Здесь и далее разностные операторы, заключенные вместе с разностной функцией в круглые скобки, предполагаются действующими только на функции внутри скобки, т. е. они не действуют на величины правее закрывающей скобки; П«, Пп — операторы диффузии, включающие физическую и искусственную диффузию четвертого порядка:

П«^ = VД"{[(Д0Жг+(/2к)2 + (ДПУг+1/2к)2](^г+1к - ^)//г+1/2к}-

Операторы Д4«, Д2« строятся по формуле

Д4« = /"^"Д+КЦгк )4/гк ]Д"Д+, Д2« = I"1 Д+[и2/гк]Д"

аналогично определяются операторы Д4п, Д2п.

Параметры Л, 7, ^ определены в процессе анализа устойчивости [5] для многомерного диффузионного уравнения переноса с постоянными коэффициентами. В результате этого анализа найдено однопараметрическое семейство, обеспечивающее абсолютую устойчивость схемы:

7 = (р/12 + г)/2, Л = — (р/12 + г)3/4, к = 3/2, ^ = ^3/2,

р = р(р/12 + г)2/16 тах[(р/12 + г)/(4г), 1], г> 0.

Здесь г — свободный параметр; р = 250/243 х т, т — число пространственных переменных, в нашем случае т = 2. В двумерных расчетах течений несжимаемой жидкости применяется набор, полученный с помощью приведенных выше формул при подстановке г = р/24, что минимизирует значение

П = р/16, Л = —3р/32, £ = ^3/2, к = 3/2, ^ = р3/1024,

р = т 250/243 = 2 250/243.

Для замыкания численного метода необходимо трижды определить скорости, а именно (мга,^га), (м(2),^(2)), для вычисления фигурирующих в схеме (6)-(8) величин Д11), Д2). Для этого следует трижды решить уравнение (2) относительно функции тока, подставляя в правую часть величины и(1), и(2). Построим дискретный аналог уравнения (2). Определив формулы четвертого порядка

+ 0(Ап4),

(Zn)ik+i/2 = (9zifc+i - 9zifc + Zik-1 - zifc+2)/(6An) = (Zn)ik+i/2 + O(An )

и используя соотношение

(azn)n = A- [aifc+i/2(zn)ifc+i/2] - A- A+i- aik+i/2(zik+i - zifc)/24 + O(A^4), можно получить разностный аналог уравнения (2)

A- [(ЛпXi+i/2fc)2 + (Лпyi+i/2fc)2](9^i+ifc - 9^ifc + ^i-ifc - ^i+2fc) ? [(x?)i+i/2fc(Лпyi+i/2fc) - (У?)i+i/2fc(ЛпXj+i/2k)]6A£

-A-A+i-[(inXi+i/2k )2 + (i;yi+i/2k )2] +

i+i/2fc

+A- [(Л?Xifc+i/2)2 + (Л?yifc+i/2)2](9^ifc+i - 9^ifc + ^ifc-i - ^ifc+2) V [(X?)ifc+i/2(Л?yifc+i/2) - (У?)ifc+i/2(Л?xifc+i/2)]6An

-A-A+i-[(i0Xjfc+i/2)2 + (i?0yifc+i/2)2] = -Wik /ifc, (10)

241ifc+i/2

где в правой части /ik определяется по формуле четвертого порядка (9). Для описания итерационных методов, использовавшихся для решения этого уравнения, перепишем его схематически:

wt2Vj_2fc + wk-iVi-ifc + W0^ + WfcVi+ik + wk2Vj+2k+

+N-2Vifc_2 + Nifc"1)^ifc_1 + Nfc'Wi + NJW2 = Rk • (11)

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

Wi_2)cik/2 + W^Ctr + Wfvs+1/2 + W^ctl'2 + W2VvS/2+ +N_2)eU + Nik_1)C_1 + Ni(1)VS+i + N2VS+2 + (C+1/2 - С)/T.k = Rik; (12)

Wk_2)c+k/2+w^ct;'2+wfvs+1/2+^ст+wfvi++k/2+

+4_2)«й1 + + NIVs+i + Nik2)^;2+-1 + («+1 - c+1/2)/nk = R* • (13)

В расчетах выявилась необходимость использования более эффективного метода. В качестве такового был применен метод неполной факторизации, который является обобщением метода [6], развитого для двумерной системы разностных уравнений с тремя точками по каждому координатному направлению (т.е. с пятиточечным шаблоном). Отметим, шаблон системы уравнений (11) имеет пять точек по каждой переменной, всего в шаблоне 13 точек. Первый шаг метода по форме совпадает с первым шагом метода переменных направлений (лишь индекс m + 1/2 заменен на m +1)

W_24m+k + w^Vm+k1 + Wk0) ^Г1 + WlVm+k + wk24m+k+ +4_2)^т_2 + + + Nk2Vm+2 + (С+1 - фп = Rk. (14)

Однако мы не определяем из этой системы уравнений переменные "0S+1. Вместо этого из нее выводятся соотношения, связывающие по три переменных. На каждой координатной линии (n = const) разностные уравнения образуют замкнутую пятидиагональную систему уравнений. Исключая с помощью метода Гаусса, например, нижние диагонали, получаем систему соотношений

•С+1 + Dk c+k1 + Dk C+k = Cik, (15)

где правая часть Сц включает также слагаемые с искомыми переменными фЦТ с предыдущего итерационного слоя. Далее рассматривается нефакторизованная система

W_2Vr+k + W^Vr+k1 + Wk0) С+1 + w^Clk1 + WL2VS£+ +Nk_2Vm+! + N_1Vm+1 + ^Vms + Nk2Vm+2 + wr1 - №v^ = Rk, (16)

в которой переменные ^™+k, ^m+k исключаются с помощью соотношений (15). В результате мы получаем преобразованную систему уравнений, которая может быть решена с помощью пятиточечных прогонок по линиям £ = const последовательно при убывающем индексе i от максимального значения до минимального. Далее этот процесс повторяется с изменением ролей переменных и направления счета (т. е. вместо исключения нижних

диагоналей в системе (14) исключаются верхние диагонали). В итоге расчетный цикл включает четыре шага.

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

т"1 = W? /6.

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

3. Результаты расчетов

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

Wk"2) = -0.1, W,^1) = -0.4, wk = -0.4, Wk2) = -0.1, Nfc"2) = -0.2,

N"1) = -0.8, Nfc) = -0.8, = -0.2, Wik0) = -W"2) - w(-1) - wd - wk2) + N"2) - Nik-1) - N? - Nik2) + d.

Параметр d определяет диагональное доминирование системы уравнений. При приближении к границам расчетной области некоторые узлы шаблона начинают выходить за границы. В этом случае соответствующие коэффициенты системы уравнений полагаются равными нулю. В качестве точного решения выбрана дискретная функция = (i - 1)(/ - i) - (k - 1)(K - k), I и K — количество узлов по каждой из переменных. Правая часть системы вычисляется в результате подстановки этой функции в уравнения (11).

В табл. 1 приведены данные о сходимости для метода переменных направлений (12), (13), в табл. 2 — для метода неполной факторизации (14)-(16). В таблицах отражены максимальное отклонение текущего итерационного решения от точного (в норме C) и невязка, получаемая подстановкой текущего решения в систему (11) (также в норме C).

Таблица1

Данные о сходимости для метода переменных направлений

d = 1 d = 1 d = 0.1 d = 0.1 d = 0.01 d = 0.01 d = 0 d=0

погрешн. невязка погрешн. невязка погрешн. невязка погрешн. невязка

4.36e3 6.03e3 4.36e3 2.46e3 4.36e3 2.10e3 4.36e3 2.06e3

2.02e2 2.28e2 2.60e3 3.77e2 3.84e2 2.33e2 4.02e3 2.51e2

9.98e0 1.08e1 1.61e3 2.27e2 3.53e2 8.37e1 3.87e3 8.08e1

4.96e-1 5.07e-1 1.01e3 1.39e2 3.26e3 7.18e1 3.74e3 4.28e1

2.48e-2 2.53e-2 6.39e2 8.70e1 3.03e3 6.33e1 3.65e3 2.71e1

1.28e-3 1.28e-3 4.06e2 5.46e1 2.82e3 5.71e1 3.55e3 2.42e1

1.64e-4 6.45e-5 2.58e2 3.45e1 2.64e3 5.18e1 3.48e3 2.17e1

1.36e-4 3.28e-6 1.65e2 2.19e1 2.47e3 4.76e1 3.40e3 2.01e1

1.35e-4 1.67e-7 1.05e2 1.39e1 2.32e3 4.38e1 3.34e3 1.85e1

1.35e-4 8.50e-9 6.73e1 8.87e0 2.17e3 4.06e1 3.27e3 1.74e1

Таблица2

Данные о сходимости для метода неполной факторизации

й = 1 й = 1 й = 0.1 й = 0.1 й = 0.01 й = 0.01 й = 0 й=0

погрешн. невязка погрешн. невязка погрешн. невязка погрешн. невязка

4.36е3 6.03е3 4.36е3 2.46е3 4.36е3 2.10е3 4.36е3 2.06е3

4.83е0 5.04е0 8.72е2 1.11е2 3.07е3 5.96е1 3.65е3 2.85е1

5.69е-3 5.86е-3 1.83е2 2.25е1 2.37е3 4.10е1 3.35е3 1.73е1

1.52е-4 2.03е-4 4.09е1 4.98е0 1.87е3 3.08е1 3.12е3 1.37е1

1.56е-4 2.09е-3 9.01е0 1.08е0 1.48е3 2.36е1 2.94е3 1.16е1

1.56е-4 2.09е-3 2.04е0 2.44е-1 1.18е3 1.85е1 2.79е3 1.01е1

1.56е-4 2.09е-3 4.59е-1 5.42е-2 9.46е2 1.46е1 2.65е3 9.05е0

1.56е-4 2.09е-3 1.04е-1 1.22е-2 7.60е2 1.16е1 2.52е3 8.22е0

1.56е-4 2.09е-3 2.33е-2 2.70е-3 6.11е2 9.24е0 2.41е3 7.56е0

1.56е-4 2.09е-3 5.07е-3 6.76е-4 4.93е2 7.40е0 2.30е3 7.01е0

В строках приводятся данные через каждые четыре итерационных шага, что соответствует одному итерационному циклу метода неполной факторизации.

Как видно, метод неполной факторизации сходится быстрее. Особенно быстро уменьшаются погрешность и невязка на первых итерациях. Естественно, при увеличении параметра й, определяющего диагональное преобладание, скорость сходимости существенно возрастает. В расчетах течений несжимаемой жидкости применяется метод неполной факторизации. Спецификой решения уравнений (11) является близость начальных данных, которые берутся с предыдущего временного слоя, и сошедшегося решения данной конечно-разностной системы. Тем не менее приходится делать сто и более итерационных циклов метода (14)—(16).

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

ш = е(1-г2)/2б/(2п)(г2 - 2), - = у - х + е(1-г2)/2е/(2п), б = 5, г2 = х2 + у2,

и =1 - в(1-г2)/2уе/(2п), V =1 + е(1-г2)/2хб/(2п), -5 < х < 5, -5 < у < 5.

Точное решение заключается в смещении исходных распределений в соответствии со средними значениями компонент скорости и = 1, V = 1. Численное интегрирование проводилось на сетке 81 х 81 до момента времени £ = 100, когда исходный вихрь переместился на 10 размеров расчетной области по каждой переменной. В процессе численного интегрирования расчетная область периодически сдвигалась вслед за перемещением вихря.

На рис. 1 приведены результаты моделирования конвекции вихря схемой (6)-(8) для чисел Куранта 1.1 и 1.7. Изображены профили завихренности плотности в сечении у = 0.

С целью изучения возможностей численного моделирования течений несжимаемой жидкости с контактными разрывами произведено решение уравнений (2), (3), (5), в которых принято т = 0.7£(/). Дельта-функция моделируется формулой (4), б = 12Ду. Расчетная область — квадрат [-1 < < 1, -1 < у < 1], в котором верхняя и нижняя границы соответствуют твердым стенкам (— =1 и — = -1 соответственно), слева жидкость втекает, справа вытекает. Если в качестве начальных данных задать плоскопараллельное сдвиговое течение, то оно является точным стационарным решением уравнений несжимаемой жидкости. Однако задается течение с возмущением /0 = у + 0.05в1п(пх), которое приводит к развитию совершенно другого течения, а именно к сворачиванию вихревой пелены.

J_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_|_|

-4 -2 0 2 4 X

Рис. 1. Задача о пассивной конвекции вихря, завихренность в сечении у = 0.

Рис. 2. Сворачивание вихревой пелены: Ь = 1.74, число Куранта 1.1, переменные ф - f; слева — сетка 89 X 87, справа — сетка 133 X 131.

На рис. 2, 3 приведены линии уровня / = ±е, где е примерно равно 0.05. Используются сетки 89 X 87 и 133 X 131. Моменту времени £ = 1.74 соответствует рис. 2, моменту времени £ = 2.32 — рис. 3. Для сравнения эта задача решена также в классической постановке в переменных ф - ш. На рис. 3, справа, приведены изолинии завихренности ш. Как видим, даже предельно простое начальное течение развивается в весьма сложную структуру.

Итак, предложенные схемы являются работоспособным инструментом исследования течений несжимаемой жидкости. Следует отметить, по удельному расходу арифметических операций они в 11 раз превосходят аналогичные схемы Рунге — Кутты третьего порядка для решения задач газовой динамики. Это объясняется трудоемкостью итерационного решения уравнения для функции тока, которое производится трижды на каждом временном шаге. Однако явные схемы высоких порядков (например, типа ЕМО) также требуют неоднократного решения этого уравнения, в то же время они нуждаются также

Рис. 3. Сворачивание вихревой пелены: I = 2.32, число Куранта 1.1, переменные ф — /; слева — сетка 89 х 87, посередине — сетка 133 х 131, справа — сетка 133 х 131; переменные ф — ш.

и в большем количестве временных шагов.

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

[1] OsHER S., SETHIAN J. Fronts propogating with curvature-dependent speed: algorithms based on Hamilton-Yacoby formulation // J. of Comput. Phys. 1988. Vol. 79. P. 12-49.

[2] Sethian J. Level set methods: Evolving interfaces in geometry, fluid dynamics, computer visions, and material Science // Cambridge Monographs on Appl. and Comput. Math. N.Y.: Cambridge Univ. Pres, 1996.

[3] Harabetian E., Osher S., Shu C.-W. An Eulerian approach for vortex motion using a level set regularization procedure // J. of Comput. Phys. 1996. Vol. 127. P. 15-26.

[4] Пинчуков В.И., Чи-ВАнг Шу. Численные методы высоких порядков для задач аэрогидродинамики. Новосибирск: Изд-во Сиб. отд-ния РАН, 2000. 232 с.

[5] Пинчуков В.И. Сравнение неявных схем Рунге — Кутты третьего порядка // Вычисл. технологии. 2002. Т. 7, № 5. С. 44-57.

[6] Зверев В.Г. Об одном итерационном алгоритме решения разностных эллиптических уравнений // Вычисл. технологии. 1999. Т. 4, № 1. С. 55-65.

Поступила в редакцию 29 января 2004 г-