CC BY
35
УЧЕНЫЕ ЗАПИСКИ Ц А Г И
Т о м IV
19 73
№ 2
УДК 533.6.011
О РАЗНОСТНЫХ СХЕМАХ ПОВЫШЕННОЙ ТОЧНОСТИ ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ НЕКОТОРЫХ ЗАДАЧ
АЭРОДИНАМИКИ
Описывается метод построения неявных консервативных разностных схем повышенной точности, пригодных для расчета течений невязкого и вязкого газов. Рассматриваются вопросы аппроксимации и устойчивости, а также способы решения разностных уравнений. В качестве примера приводятся результаты расчета течений идеального газа за отошедшей ударной волной.
1. При численном исследовании различных классов течений, как правило, применяются схемы либо первого, либо второго порядка точности на гладких решениях исходных уравнений. Однако часто бывает желательным использовать схемы более высокого порядка точности. Такая необходимость возникает, прежде всего, при расчете течений вязкого газа, когда погрешность аппроксимации инерционных членов может стать сравнимой с членами, содержащими коэффициент вязкости. Применение более точных схем может оказаться эффективным и при расчете течений идеального газа, при этом желаемая точность может быть достигнута при сравнительно крупных шагах расчетной сетки.
В статье [1] указаны способы построения консервативных разностных схем до третьего порядка аппроксимации по всем переменным, некоторые из них использовались для численного решения уравнений Навье — Стокса сжимаемого газа. Ниже рассматриваются вопросы, связанные с аппроксимацией и устойчивостью этих схем и, в качестве примера, приводятся результаты расчетов течения идеального газа за отошедшей ударной волной перед затупленным телом.
2. Чтобы построить схемы третьего порядка аппроксимации для простейшего уравнения
А. И. Толстых
введем на сетке xi = ih операторы Л(+ и ЛИ' по формулам: ■*7 + 1
f (х) с/х ---= /г /’,■ Fi+1 — ^2 +
+ 0(/г4) = /г(Л^)/7). L+0(ft4),
‘ + 2
f F(x)dx = h(^ Fi+1 + 4- ^ - it ^-i) + 0 (*4) =
^hiA^F) ,+0(/г4). г+ У
Разностные схемы для уравнения (1) запишем в виде
(О
*7+1
Аф/«
/г
+ (Л<4^) ,=о, * 2
h
а !!>/,
/г
+ (л!!>5) ,=о,
' +т
+ (л(:,§-); 1=о,
где
(<).
(За)
(36)
(Зв)
Подставив в (За) — (Зв) достаточно гладкие функции f н g и разложив их в ряд Тейлора, например, в окрестности точек х , и хг, получим для погрешности аппроксимации уравнений (1)
/+~2
соответственно выражения
г+Т
/г d / df ,
2 dx \ d.* 1
+ -Ti£(lr + s) +0(лз)*
В силу (1) все выражения в скобках и производные от этих выражений исчезают, так что схемы (За) — (Зв) аппроксимируют уравнение (1) с третьим порядком точности на гладких решениях последнего. Отметим, что схемы (36) и (Зв) выделены из совокупности схем (За) для выяснения характера аппроксимации в узле хс.
Для уравнения
df(u, х,у) дх
> (и, х, у)
ду
+ -ь (и, х, у) = 0
(4)
на сетке х1 = й, У] = у7 можно получить аппроксимацию вида О (/г3 + /*), заменив производную ду/ду при х = х1 каким-либо разностным отношением (дер/ду) с порядком точности & и приме-
нив операторы Д^ и Л^ соответственно к и
Схемы (3) в этом случае запишутся в виде
+{А^т/ду) +ф))і+_ = о,
Щрі- + {А(і)((д<?Іду)+'ї))і+г = О,
1±Т ' /
а погрешности аппроксимации, например, в случае схемы (36) будут иметь вид
В качестве простейших выражений для (д<?!ду) можно использовать двух- или трехточечные формулы
Такая „неравномерная" аппроксимация по всем переменным имеет смысл в тех случаях, когда в одном из направлений желательно иметь высокий порядок точности, например, вдоль нормали к телу в случае вязких течений, а для других направлений это не столь существенно.
Для аппроксимации уравнения (4) с третьим порядком по всем переменным используем операторы Л^, Д-р и Л^?, Д^\ соответствующие координатам х я у; тогда одну из совокупности полученных схем можно записать, например, в виде
Разлагая достаточно гладкие функции / и ? в ряд Тейлора
которое для /, ? и ф, удовлетворяющих уравнению (4), дает члены порядка О (Л3 4- Р).
Схемы третьего порядка точности по всем переменным, построенные на основании идей метода Рунге — Кутта, рассматривались в работах [2] и [3]. Большое количество узлов шаблона для каждой пространственной переменной, по-видимому, усложняет применение их для решения краевых задач.
Отметим, что схемы (3), (5) и (7) можно получить другим путем, минуя соответствующие дифференциальные уравнения. Для этого достаточно составить баланс потоков через границы элементарных ячеек, применяя квадратурные формулы (2).
к д ТШ
(66)
(6а)
в окрестности точки получим для погрешности аппро-
ксимации выражение вида
Консервативность описанных выше схем позволяет, по-видимому, применять их и для описания разрывных решений исходных уравнений (например, для расчета течений идеального газа с ударными волнами).
Для аппроксимации уравнения второго порядка, записанного в дивергентном виде (например, в задачах вязкого обтекания), можно применить описанный выше способ построения разностных схем, либо сведя это уравнение к системе уравнений первого порядка, либо считая выражения для потоков через соответствующие границы ячеек зависимыми от первых производных искомых функций.
В частности, для простейшего уравнения
(8,
гг ди
полагая /=/—£—, аппроксимируя производную со вторым порядком точности и используя схемы (3), получим аппроксимацию на решении (8) вида О (шах (ек2, Л3)), что при малых е, по крайней мере в области слабых градиентов функции и, дает практически высокий порядок точности.
3. При применении схем (5) — (7) для решения краевых задач возникает вопрос о выборе оператора А+ или А-. Такой выбор можно сделать на основании следующих соображений: во-первых, используемые узлы сетки не должны выходить за границы рассматриваемой области, и, во-вторых, схема должна быть устойчивой по крайней мере в линейном приближении для задачи Коши.
Для уточнения последнего требования исследуем устойчивость схем (5) — (7) применительно к простейшему линейному уравнению
ди ди , /пч
Ж = а Тх > а ^ соп^‘ (9>
На сетке £т = ягт, хп = пк рассмотрим сначала аппроксимацию (5), (6а), которая, например, при а>0 запишется в виде
т+1 _ т т+Х __ т „т+г „т ит+х — мт+1
5 “л+1 ип+1 , 8 «п ип 1 п—1 —ип-1 „ ип+1 ип мт
12 х ~г 12 1 12 х к ' ^
Разыскивая собственную функцию оператора перехода от т-го слоя к (яг + 1)-му в виде ип = еіап, получим для модуля собственного числа X выражение
1Х!2= (1— а>1*е№*)2 +(1ш И7а)2 > О1)
где & здесь и в дальнейшем — это минус или плюс; г—х/Н; , ]Х^+ — функции комплексного аргумента, соответствующие операторам А- и А(+:
еы — 1 <>-1 .
= -=------------------------------ , Ц7+ =
5 І* _1_ 8 1 -и ’ “ 5 8 и 1 91*
------- Р і — ______________________ _______ Р *а .1 I ------ Л1а ________ - - pH ®
12 е 12 12 12 12 12 е
Легко усмотреть, что Ие1^_<0 и Не1Р + >0, так что для оператора Л*"*, обеспечивающего выполнение неравенства а Ие Ц7Й<0, схема (5), (6а) является абсолютно устойчивой. ;
В случае аппроксимации (5), (66), а также схемы третьего порядка типа (7) выражения для собственного числа X приобретают соответственно вид
5(Х) —Хаг ИГ* = 0, 5 (X) =
2+21-
(12а)
X-1-/?(Х)аг 1^ = 0, Я(Х) =
12
х+.
12
1
12 X
(126)
Очевидно, что для выполнения соотношения |Х|< 1 также необходимо удовлетворить условию аЯеМ/лА<0. Для проверки достаточности последнего были проведены вычисления значений X при различных значениях аргумента а(0<а-<2тг) и параметра \а\г. На фиг. 1 показано изменение величины А = тах|Х| с изменением
0<а<2п
\а\г для схем (5), (6а) и (7). Видно, что при аКей^А<;0 схема (5), (6а) абсолютно устойчива, а схема (7) — условно устойчива (устойчива при ]а(г^1).
В случае уравнения второго порядка вида
і
✓ И
Область неустойки <{120)
\\\ 1,\\' ЛЧ' XV Ч\\' \\\
{12а)
ди
д(
— а
ди
дх
д2 и
(13)
\а\ г/7і
схемы (5) порядка с
Фиг. 1
(7) для производных первого надлежащей аппроксимацией производной д2 и/дх* на (т + 1)-м временном слое при аИеи^СО образуются абсолютно устойчивые неявные схемы.
В случае применения разностных схем (5) — (7) к уравнениям газовой динамики оценку устойчивости, необходимую для выбора нужных индексов используемых операторов, можно получить из разностных уравнений для соответствующей „недивергентной“ системы с „замороженными" коэффициентами. В качестве простого примера рассмотрим схемы (5), (6 6) для одномерных нестационарных уравнений газовой динамики, записанных в переменных Эйлера в виде
др_ дt
I дРы дри
' дх ’ ді
д(?и2 + р)
дх
0.
до
к и2
Т + ~2
дх
+
дри ( к + ~2~ дх
(14)
= 0, р —
г- 1
где и, р, р, к, — соответственно скорость, плотность, давление, энтальпия и показатель адиабаты.
Пусть в некоторой области выполняется условие \и\~>с = = — 1)А, тогда, вводя на сетке Ьт — 1т, хп = пЬ оператор Л*л)
таким образом, чтобы «ИеН^^О, получим для собственного числа X оператора перехода к последующему временному слою уравнение вида
[5 (X) + иЖк г\ {[5 (X) -\-uWk г]« - с* \РІ г2} = 0,
где г = т/8, а 5 (X) — функция, определяемая в формуле (12а). Это уравнение распадается на ^ри уравнения вида (12а) с коэффициентами а ——и и а = — (и + с), каждое из которых, согласно предыдущим результатам, определяет значения 'к, не превосходящие по модулю единицы.
В областях, где ]и|<с, при аппроксимации первого и третьего уравнений системы (14) используем оператор А[п), для которого Не1^>-0, а при аппроксимации второго уравнения изменим знак в выбранном индексе на противоположный. Тогда, замечая, что 1ш Т&+ = 1ш V?-, получим для определения величин ). уравнения
Рассматривая второе из этих уравнений как квадратное относительно 5 (А), придем к двум уравнениям вида (12 а) с чисто мнимыми значениями коэффициента а, обеспечивающими выполнение условия |Х|< 1. Таким образом, схема (5), (6 6) для рассмотренного способа аппроксимации системы (14) является абсолютно устойчивой в приближении „замороженных11 коэффициентов.
При решении систем разностных уравнений в случаях, когда происходит переход от одного слоя к другому, можно воспользоваться методом прогонки. Рассмотрим для простоты случай смешанной задачи для уравнения (9) в полосе 0^х-<1:
Разностные уравнения при применении схемы (5), (6 а) (аИеи^_<.0) можно записать в виде
где через (1п обозначены члены, содержащие значения функций на т-м слое. Второе граничное условие для (14) можно получить из системы двух уравнений, записанных при п= 1 и 2 с учетом изменения индекса оператора Ак на противоположный. После исключения величины исГ+1 получим соотношение
где через обозначены члены, зависящие от й2 и йА. Используя выражения для прогоночных коэффициентов, можно показать, что уравнение (15) с условием (16) позволяют реализовать устойчивую прогонку.
Отметим, что возможен также итерационный способ решения системы (15), при котором значение и„-\ берется из предыдущей итерации и уравнения сводятся к двухточечным; можно показать, что собственные числа соответствующего оператора перехода меньше единицы и условие сходимости выполнено. Применение итерационного процесса имеет смысл в случае нелинейных уравнений.
8{1)±иЖкг = О,
Я2 (X) + 2 иг (1т \Рк) 5 (X), 4- | \ * х* (с2 - и2) = 0.
и {х, 0) = 9(*), «(1, *) = ф(0 (а>0).
(15)
а2 г2 4- аг-\
(16)
3
4. Рассмотрим коротко вопрос о применении схем описанного типа к расчету стационарного осесимметричного течения совершенного газа между затупленным телом и отошедшей ударной волной. Поверхность ударной волны будем считать поверхностью разрыва, на которой выполняются соотношения Гюгонио. Выбрав систему координат ял, связанную с контуром тела (координаты в и п направлены соответственно вдоль поверхности и по нормали к ней) и введя новую переменную 1 — п1г, где е — расстояние отхода волны, запишем уравнения движения и неразрывности в „дивергентном" виде (4), понимая под / и ср следующие векторы:
Здесь и, V — касательная и нормальная составляющие скорости, V — некоторая комбинация и и V, возникающая из-за введения переменной Н — коэффициент Лямэ. Роль координат х и у в (4) играют соответствено X и «. В качестве уравнений, замыкающих систему, будем использовать интеграл Бернулли, уравнение состояния и геометрическое соотношение для скачка уплотнения, выражающее связь между величиной е и углом наклона волны к телу. В полосе 0<£<;1 [точка (0,0) расположена
в критической точке тела] на сетке «(==/Д«, = (/ = 0, 1
/' = 0, 1, ..., N1 введем аппроксимации (5), (6 6)* с /г = Д£ и /=== Дя.
Для решения системы разностных уравнений использовался итерационный метод с релаксацией, эквивалентный некоторой неявной схеме метода установления. Под итерацией будем понимать расчет всего поля течения по неявной схеме, аналогичной схеме (10) для эволюционного уравнения (9); при этом роль временной переменной Ь играет координата в (5 =-= 0, Дз, ..., УИДв). Окончательный переход к следующей (т + 1)-й итерации для любой функции Т7 осуществлялся по формуле
уравнений, со — релаксационный параметр. Для различных расчетных вариантов применялась как „одновременная", так и „последовательная" релаксация; в последнем случае значения функций релаксировались сразу же после того, как они определялись во всех точках луча 5 = 5г.
Последовательное вычисление искомых функций при переходе с одного луча на другой можно реализовать различными способами. Использовался следующий прием: общая система разностных уравнений разбивалась на подсистемы, соответствующие каждому дифференциальному уравнению; в каждой из этих подсистем некоторые газодинамические функции считались известными, т. е. взятыми либо из предыдущей итерации, либо из найденных в течение текущей итерации решений других подсистем. Например,
* Производные др/дв аппроксимировались с помощью центральных разностей.
/ = (рVII, рот/ + р, рг>).
в разностных уравнениях, соответствующих уравнению движения в проекции на ось s, считались неизвестными значения скорости ии и известными все остальные функции ру, vtj, hц. Найденные таким образом значения во всем поле после релаксации использовались в остальных уравнениях и т. д.
Специфика решения разностных уравнений вдоль лучей s = const состояла в следующем. Значения функций, умноженные на коэффициент 1/12, предполагались известными из предыдущей „внутренней" итерации для каждого луча, так что трехточечные разностные уравнения сводились к двухточечным. Это позволяло вести „бегущий счет“ от волны или от тела с последующими внутренними итерациями. Практически иногда достаточно было вообще не итерировать значения функций на луче, а величины, перенесенные в правую часть уравнений, брать из предыдущей „внешней" итерации. ’
Разностные уравнения, аппроксимировавшие уравнение движения в проекции на ось s, рассматривались как квадратные относительно значений скорости utj и решались последовательно от волны к телу.
Уравнения неразрывности и движения в проекции на ось п считались линейными соответственно относительно скорости V и плотности р. Значения v,j определялись последовательно от тела к волне, значения р^ — от волны к телу; величина отхода е находилась из разностной аппроксимации уравнения неразрывности.
Значения функций, вычисленные на i-м луче, использовались в разностных уравнениях для (г+1)-го луча. На луче s = 0 использовались условия симметрии и производные (du/ds)s=о, вычисленные из предыдущей итерации; при S = S[ всегда, без потери точности вследствие условий симметрии, применялись двухточечные разностные аппроксимации производных по s; на луче s = sM эти производные записывались с помощью левых разностных отношений.
Выбор операторов А+’ или А-\ если он оказывал влияние на устойчивость схемы в линейном приближении, осуществлялся так же, как и в случае уравнения (9), с той лишь разницей, что вместо (9) рассматривались „недивергентные“ уравнения с „замороженными" коэффициентами для соответствующих независимых уравнений или систем.
Выбор индекса & оператора А{Р в известной степени аналогичен выбору направления несимметричных разностей в зависимости от знаков скоростей в несимметричных разностных схемах; при этом ввиду нелинейности уравнений окончательное заключение об устойчивости используемой аппроксимации можно получить только в результате численных экспериментов.
Расчеты проводились при значениях релаксационного параметра Фиг. 2
— И
• сетка M=10-,N=J0 х » M=20-,N=10
<в = 0,1; число итераций для получения трех совпадающих значащих цифр в течение последних 50 итераций составляло 200 — 300.
В качестве примера на фиг. 2 изображены положения ударной волны и звуковой линии при числе М«, = 3 и 10 и показателе адиабаты ? == 1,4, вычисленных на различных разностных сетках (5лг=Ь®)! там же приведены данные работы [4]. На фиг. 3 показано распределение давления вдоль контура тела в сравнении с результатами [4].
Как видно из фиг. 2 и 3, полученные данные достаточно близки к общепринятым результатам и мало изменяются при увеличении
шага Дя. Для контроля за точностью расчетов можно использовать значение энтропии 5 = /?/рт на нулевой линии тока, поскольку уравнение сохранения энтропии нигде не использовалось в исходной системе.
На фиг. 3 приведена функция 5 (в), изменяющаяся в интервале 0, не более чем на 4%.
ЛИТЕРАТУРА
1. Толстых А. И. О методе численного решения уравнений Навье — Стокса сжимаемого газа для широкого диапазона чисел Рейнольдса. Доклады АН СССР, т. 210, № 1, 1973.
2. Русанов В. В. Разностные схемы третьего порядка точности для сквозного счета разрывных решений. Доклады АН СССР, 180, № 6, 1303 (1968).
3. В и г s t е i n S. Z., М i г i п A. Third order difference methods for hyperbolic equations. Труды секции по численным методам второго международного коллоквиума по газодинамике взрыва и реагирующих систем. Т. 1. М., ВЦ АН СССР, 1971.
4. Белоцерковский О. М. Расчет обтекания осесимметричных тел с отошедшей ударной волной (расчетные формулы и таблицы полей течения). М., ВЦ АН СССР, 1961.
Рукопись поступила 9jVI 1972 г.
