CC BY
44
Вычислительные технологии
Том 3, № 3, 1998
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ НА АДАПТИВНЫХ СЕТКАХ ДВУМЕРНЫХ УСТАНОВИВШИХСЯ ТЕЧЕНИЙ ЖИДКОСТИ И ГАЗА *
Н. Ю. ШокинА
Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: shokin@adm.ict.nsc.ru
Numerical modeling of two-dimensional steady-state flows of fluid and gas using adaptive grids. A finite-difference method of calculating of plane steady-state flows of ideal incompressible fluid and gas in the channels of complicated geometry with inlet, outlet and impermeable walls is presented. For numerical solving of the problem new dependent variables — stream function ф and vorticity и are introduced. The problems of constructing a curvilinear grids, obtaining of finite-difference approximations on these grids, calculating of vorticity on the inlet are discussed. The results of calculating of test problem about the steady-state non-isoenergetical vortex gas flow in curvilinear channel.
Введение
В настоящее время моделирование течений газа успешно осуществляется на основе уравнений Навье — Стокса [1-4]. Тем не менее имеется достаточно широкий круг задач, которые удовлетворительно решаются и в рамках модели идеального газа, так как эта модель хорошо передает некоторые интегральные и локальные характеристики потока [5]. Кроме того, результаты расчетов на основе уравнений Эйлера можно использовать в качестве начального приближения для итерационных методов решения стационарных уравнений Навье — Стокса, а также для получения предварительной картины течения и оценки учета вязкости газа. В силу сказанного разработка надежных и эффективных алгоритмов расчета течений идеального газа остается актуальной задачей.
В представленной работе рассматриваются дозвуковые установившиеся плоские течения идеальной сжимаемой жидкости в каналах сложной формы. Для численного решения задачи вводятся новые зависимые переменные — функция тока ф и завихренность ш. Известно [3], что при расчете течений вязкой несжимаемой жидкости такой подход дает два уравнения эллиптического типа относительно ф и ш и требует задания ф и вычисления ш на всей границе области. В нашем случае дозвукового течения идеального газа уравнение для функции тока также будет иметь эллиптический тип, а для вихря получается уравнение первого порядка гиперболического типа, поэтому при численном решении
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №97-01-00819.
© Н. Ю. Шокина, 1998.
потребуются значения вихря не на всей границе, а только на входных участках, через которые траектории входят в область течения.
Алгоритм расчета течения газа с использованием функций ф и ш изложен в [6] при предположении, что геометрия канала не является слишком сложной. Такое ограничение было связано с использованными в [6] прямоугольными сетками. В работе [7] этот алгоритм обобщен для криволинейных сеток, согласованных с границами канала.
В настоящей работе, в отличие от [7], расчеты ведутся на криволинейных сетках, адаптирующихся не только к границам криволинейного канала, но и к некоторой априорно заданной функции, которая выбирается так, чтобы происходило измельчение ячеек сетки в предполагаемых подобластях больших значений градиента решения и, наоборот, укрупнение ячеек там, где решение слабо меняется. Для построения таких сеток используются уравнения двумерного варианта метода эквираспределения. Кроме того, в работе сравниваются некоторые методы расчета вихря на входе в область.
1. Математическая модель
Пусть П — односвязная область в плоскости декартовых координат х10х2. Математическая постановка задачи об установившемся протекании газа через область П заключается в определении вектора скорости V, давления р, плотности р, полной энергии Н, удовлетворяющих в П уравнениям неразрывности, движения и энергии:
V- рV = 0,
(рТ? ^^ + Vp = /, х = (х1,х2) € П, V ■ руН = 0, и краевым условиям на границе Г = дП :
рУ
хеГх
и1 (х), V ■ п
0, рУ ■ и
хеГо
^г(х),
хеГ2
Н
хеГх
Н (х),
(1.1)
(1.2) (1.3)
(1.4)
(1.5)
где V = (м1,м2), иа(а = 1, 2) — компоненты вектора скорости в направлении осей 0ха, дд
V = ( т^-т, тг-^т ), п — внешняя нормаль к Г, Г = Г1 и Г0 и Г2, Г1 — вход в П (V ■ и < 0),
дх1 дх2_
Г0 — непроницаемая часть границы, Г2 — выход из П (и2 > 0), Г1 П Г2 = 0.
Газ предполагается совершенным, поэтому для полной энергии имеем следующее выражение:
V 2
Н
тр
+
р(7 - 1) 2
(1.6)
где 7 — показатель адиабаты, V = |у |.
В дополнение к условиям (1.4), (1.5) будем считать, что в некоторой точке М0 € П известно давление:
р(Мо)= ро. (1.7)
Для численного решения задачи введем новые зависимые переменные — функцию тока ф и функцию вихря и:
рщ =
дф
дх2 дх1
Тогда для ф и ш получим уравнения
д /1 дф
- дф ш = rotV = - — + —
дх2 дх1
дхг \ р дхг '
i=i
.V2
V • руш = -rot(PV-^ + rot/, при этом граничные значения для ф однозначно определяются из условий (1.4).
(1.8)
(1.9) (1.10)
2. Переход к криволинейным координатам
Пусть
ха = ха(д\д2) (2.1)
есть взаимно-однозначное невырожденное отображение квадрата (, лежащего в плоскости координат д10д2, на область П. Ради простоты будем предполагать, что участки входа и выхода являются связными кусками границы, причем Г1, Г2 — это образы при отображении (2.1) соответственно левой 71 и правой 72 стороны квадрата а непроницаемая часть Го состоит из двух участков Г0 и Г0, являющихся образами соответственно нижней 70 и верхней 70 сторон (. Область П и ее прообраз ( показаны на рис. 1.
Рис. 1. Схема физической области Q и вычислительной области Q.
Уравнения (1.9), (1.10) и (1.3), записанные в новых независимых переменных да, выглядят следующим образом:
д (, дф дф\ д ( дф дф\ . . тгт к11^П + + тг^ кпт-т + = -ии, (2.2)
дд1 \ дд1 дд2) дд2 \ дд1 дд2)
А( А(У1Х\ -/ / (23)
д^1 ^д^Н + дд2 Удц1 V - дд2 + Ъд1, ( ^
дд
дф (р3у1Н) + дф (р^Н) = 0, (2.4)
где
, _ д22 , _ , _ д12 , _ 911
к11 = -7, к12 = к21 =--7 , к22 = -7,
р. р.! р.
9ав — компоненты метрического тензора, а, в =1, 2,
911 = х?1 + , 912 = V V + Уд1 Уд2 , 922 = х^2 + ^2 ,
.] — якобиан преобразования (2.1), .] = хд1 уд2 -хд2уд1, уа — контравариантные компоненты скорости, связанные с ф соотношениями
^ = дф ^ =__дф (25)
р.1 дд2 р.1 дд1) '
/а — ковариантные компоненты вектора /.
При вычислении давления будут использоваться уравнения движения в форме Громе-ки — Лэмба, записанные в криволинейных координатах:
т 2 , д V V дР ,
+ р—т + ^ = Ь,
д V2 др
р^ш + + дф = /2 (2.6) 3. Разностные уравнения
Разностные уравнения для сеточных функций ф, ш, Н получаются интегро-интерполяци-онным методом путем аппроксимации интегральных соотношений
/ (кп+к12ъдф) ^ - {к21 ^+к22^ф) ^ =- / / Jшdg1dg2, (3.1)
С Ба
£ рJv1шdq2 - рJv2шdq1 = - £ ^рд^^Т - dq2 + £ (^р-^^ - (к1, (3.2)
С С С
£ рJv1Hdq2 - = 0, (3.3)
С
где Бс — прямоугольник с контуром С, изображенным штриховой линией на рис. 2 и 3.
Аппроксимация соотношения (3.1) такая же, как и в случае несжимаемой жидкости [13], при этом разностное уравнение для ф получается девятиточечным, а функции ф, р, 9ав предполагаются определенными в целочисленных узлах сетки.
Предполагается, что сеточные функции ш и Н определены в центрах ячеек сетки, покрывающей вычислительную область Q. При аппроксимации (3.2) учитываются знаки контравариантных компонент скорости и ш берется либо из центра ячейки, охватываемой
Рис. 2. Шаблон разностного уравнения для функции тока.
Рис. 3. Шаблон разностных уравнений для вихря и полной энергии.
контуром С, либо из центра соседней ячейки. При этом для вычисления величин рЗп0, используются разностные аналоги соотношений (2.5). Так, при аппроксимации интеграла в левой части (3.2) по стороне ВС полагается
(Е) = Г ы0, если у1(Е) > 0, Ы(Е) = \ ыз, если у1(Е) < 0,
и^л ф(с) - ф(в)
pJv1(E) =
В результате получается разностное уравнение, аппроксимирующее (2.3) с первым порядком. Шаблон этого уравнения может состоять из одного, трех или четырех узлов, т.е. является переменным. Аналогично [13] можно показать, что при отсутствии замкнутых линий тока и точек покоя газа величину ы в любом узле можно вычислить через значения вихря на входе в область. Для этого можно использовать безитерационный метод — обобщенный метод бегущего счета [13].
Соотношения (3.2) и (3.3) отличаются лишь правыми частями, поэтому разностные уравнения для вычисления Н отличаются от уравнений для вихря только тем, что являются однородными:
4
Нк = 0, (3.4)
к=0
где
= Г -ф(Я) + ф(А), ф(Б) >ф(А), = Г ф(В) - ф(А), ф(А) >ф(В), I 0, ф(Б) < ф(А), в2 \ 0, ф(А) < ф(В),
= Г ф(С) - ф(В), ф(С) <ф(В), = Г ф(Б) - ф(С), ф(С) >ф(Б), I 0, ф(С) > ф(В), в4 \ 0, ф(С) < ф(Б),
в0 = -(в1 + в2 + вз + в4),
при этом
вк < 0, к = 1, 2, 3,4; в0 > 0.
Давление определяется на основе следующего выражения, вытекающего из уравнений (2.6): 2
Р(Яз)= Р0 + J (^1 + -
7(® ,<Ц)
/ д V 2\
+ (/2-р^1ш -р—Ф^) (3.5)
Будем считать, что сеточная функция р определена в целых узлах. Тогда в качестве кривой 7, соединяющей точку д^, в которой требуется определить давление, с точкой д0, в которой давление задано изначально в (1.7), будем брать ломаную, звенья которой параллельны осям 0да и проходят по сторонам ячеек. Сравнивая (3.5) и (3.2) видим, что подин-тегральные выражения в этих интегральных соотношениях одинаковы, поэтому давление не будет зависеть от пути интегрирования 7 в случае, если интегралы в (3.5) аппроксимировать точно по тем же формулам, которые применялись при выводе разностного уравнения для вихря. Тем самым для вычисления давления используется метод согласованной аппроксимации, для случая прямоугольных сеток рассмотренный в работе [9]. Для течений вязкой несжимаемой жидкости он известен как маршевый метод расчета давления [8].
Плотность вычисляется на основе (1.6), при этом производится переинтерполяция полной энергии в целые узлы.
4. Итерационный процесс
Криволинейная сетка строилась до начала итерационного процесса и в ходе итераций не менялась. Координаты узлов рассчитывались на основе двумерного аналога метода эквираспределения [12]:
д ( дха\ д ( дха\
ид? г922и?) + ^911 д^) =°, а = 1,2, (4.1)
где — управляющая функция. При = 1 получается криволинейная "квазиравномерная" сетка. При переменном сетка измельчается в тех подобластях, где принимает большие значения. На границе области П для расстановки узлов использовался одномерный метод эквираспределения с той же управляющей функцией ,ш.
При аппроксимации (4.1) применялся интегро-интерполяционный метод. Полученная система разностных уравнений решалась итерационным методом переменных направлений.
Опишем порядок вычислений в итерационном процессе решения разностной задачи. Пусть фп, шп+1/2, Нп+1/2, рп, рп — п-е итерационное приближение. Процесс получения и + 1-го приближения разбивается на следующие этапы.
I. Расчет течения идеальной несжимаемой жидкости, которое будет начальным приближением для расчета течения газа.
1. Методом последовательной верхней релаксации решается система разностных уравнений для ф™+1.
2. Вычисляются граничные значения вихря на входе Г1 . Для этого используется выражение для вихря в криволинейной системе координат
1 ( дvl дv2^ _
ш = ! (-1? +1^) ■ <4-2>
Ковариантные компоненты скорости вычисляются по заданному на входе расходу (1.4). При вычислении касательной производной дv1/дq2 используются компоненты вектора и1
и значения плотности рп. При вычислении нормальной производной дь2/дд1 используется величина v2, определенная в центре приграничной ячейки по значениям ф™+1. Разностный аналог формулы (4.2) применяется для определения предварительного значения Со, а окончательные значения вихря на входе получаются релаксированием Со и граничных значений вихря с предыдущей итерации:
сп+1
5ШС + (1 - 5Ш)сп , (4.3)
Г1
где 0 < 5Ш — параметр релаксации [10], 5Ш = 0(к1).
3. Обобщенным методом бегущего счета решается система уравнений относитель-
п+1
но вихря с^+1/2 во внутренних узлах.
Итерации ф - с для жидкости проводились до тех пор, пока ип+1 и сп различались на величину, большую некоторого заданного еш. II. Расчет течения газа.
4. Обобщенным методом бегущего счета вычисляется полная энергия Нп_+11/2.
5. Методом согласованной аппроксимации находятся значения давления рп+1.
6. Определяется плотность рп+1, при этом на основе формулы (1.6) получаются предварительные значения р, а в качестве окончательных значений плотности берутся их значения, интерполированные по р и рп:
рп+1 = V +(1 - ¿рК+1> (4.4)
где параметр 5Р > 0 линейно возрастает от нуля до единицы по ходу итерационного процесса.
Итерации для газа продолжаются до сходимости значений плотности. Тестовые расчеты проводились для задачи с известным точным решением, приведенным в [14]. В этой работе решение выписано в неявном виде, поэтому приведем все формулы, по которым вычисляется точное решение. Оно определяется через функцию Г(х), являющуюся решением следующего обыкновенного дифференциального уравнения:
Г2 (1 - в г-Г Г7"1) Гх =--т-г^--, (4.5)
Х 7 (1 + 1) 1
в1 у + ) Г7-1 - 1 2(7 - 1)
где 7 — показатель адиабаты, в — некоторая постоянная, определенная ниже. Уравнение (7.1) дополняется краевым условием
Г (хо) = Го. (4.6)
В работе [14] указано, что в зависимости от выбора Г0 могут получиться решения, соответствующие дозвуковым или сверхзвуковым течениям. Точное решение вычисляется по формулам:
А +
и = уаV = -уа+1А^,
р = у2Ь ГА~а, р = С2Г7, Н = у2аАа(^ сз + ^ А), (4.7)
г
1
где A(x) = 2c3F2 ^ 1 — в—F1 , a — произвольная постоянная, Cj > 0 (j = 1, 2, 3), а постоянная в вычисляется по Cj как
C2
в
С1С3
Алгоритм получения точного решения для тестирования в нашей работе имел следующий вид. Сначала методом Рунге—Кутты 4-го порядка численно находилась функция Г на интервале [х0,х1], затем по полученным значениям Г вычислялись линии тока по формуле
У = ь • А-1/2,
где Ь — произвольная константа. Два значения Ь1 и Ь2 дают две линии тока, ограничивающие тестовую область. Для данного тестового примера принимались следующие значения параметров:
7 = 1.4, в = 0.3, Г0 = 0.86, х0 = 0, х1 = 0.8, С1 = 1, сз = 1, Ь1 = 0.05, Ь2 = 0.2.
Найденная область приведена на рис. 1. Затем в данной области строилась сетка, на которой вычислялось точное решение по формулам (4.7).
Полученные результаты расчетов течения газа показывают первый порядок точности описанного алгоритма.
Автор выражает глубокую благодарность Гаязу Салимовичу Хакимзянову за постоянное внимание к работе.
Список литературы
[1] КОВЕНЯ В. М., ЯНЕНКО Н. Н. Метод расщепления в задачах газовой динамики. Наука, Сиб. отд-ние, Новосибирск, 1981.
[2] ФЛЕТЧЕР К. Вычислительные методы в динамике жидкостей. Т. 2. Мир, М., 1991.
[3] РОУЧ П. Вычислительная гидродинамика. Мир, М., 1980.
[4] Андерсон Д., ТАННЕХИЛЛ Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. Т. 1-2, Мир, М., 1990.
[5] ПЕйРЕ Р., ТеЙЛОР Т.Д. Вычислительные методы в задачах механики жидкости. Гидрометеоиздат, Л., 1986.
[6] ХАКИМЗЯНОВ Г. С., ЯУШЕВ И. К. Итерационный метод расчета двумерных дозвуковых установившихся внутренних течений идеальной сжимаемой жидкости. Препринт/АН СССР. Сиб. отд-ние. Ин-т теорет. и прикл. механики, №4-87, Новосибирск, 1987.
[7] Shokin Yu. I., Khakimzyanov G. S. Numerical calculation of stationary subsonic gas dynamic problems. Proc. of the Eighth GAMM-Conf. on Numerical Methods in Fluid Mechanics. Notes on Numerical Fluid Mechanics. Vol.29, Veiweg, Braunschweig, 1990, 513-521.
[8] RICHARDS C. W., Crâne C. M. Pressure marching schemes that work. Int. J. Numeric. Meth. Eng., 15, 1980, 599-610.
[9] Хлкимзянов Г. С., ЯуШЕВ И. К. О расчете давления в двумерных стационарных задачах динамики идеальной жидкости. Журн. вычислит. математ. и мат. физики, 24, №10, 1984, 1557-1564.
[10] ТАРУНИН Е. Л. О выборе аппроксимационной формулы для вихря скорости на твердой границе при решении задач динамики вязкой жидкости. В "Численные методы механики сплошной среды", 9, №7, 1978, 97-111.
[11] CHRISTOV C. I. Orthogonal coordinate meshes with managable jacobian. Numerical Grid Generat., Appl. Math, and Comput., 10/11, 1982, 885-894.
[12] ХАКИМЗЯНОВ Г. С., ШОКИНА Н. Ю. О методе эквираспределения для построения двумерных адаптивных сеток. В "Вычислительные технологии". ИВТ СО РАН, Новосибирск, 4, №13, 1995, 271-282.
[13] ХАКИМЗЯНОВ Г. С., ШОКИНА Н. Ю. Численное моделирование на адаптивных сетках двумерных установившихся внутриканаловых течений идеальной жидкости. Там же, 283-294.
[14] Андреев В. К., Капцов О. В., Пухначев В. В., Родионов А. А. Применение теории групповых методов в гидродинамике. Наука, Сиб. отд-ние, Новосибирск, 1994.
Поступила в редакцию 12 декабря 1997 г.
