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

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

Том 6, № 2, 2001

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

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

The results of the numerical investigations are presented. Test problems are considered as well as the calculation problem of turbulent trans-sonic flows. Three high-accuracy implicit schemes are used: Runge — Kutta scheme of the third order, Runge — Kutta scheme of the second order with the order increased to the fifth for stationary solutions, and the scheme of the third order with respect to time and the third order with respect to the space directions. All the schemes include artificial diffusion of the same type. The diffusion properties are analysed using Burgers equation by way of example.

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

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

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №99-01-00560.

© В. И. Пинчуков, 2001.

осцилляционными свойствами из класса формул одного порядка для определения каких-либо элементов схем, как правило, потоков через границы разностной ячейки. В настоящее время этот подход ограничен построением явных схем. Подробный обзор ЕМО-методов, а также других методов высоких порядков содержится в [1].

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

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

А — неявная схема Рунге — Кутты [1, с. 72-82, 2] третьего порядка по времени и пространственным переменным со скалярным стабилизирующим оператором;

Б — неявная схема Рунге — Кутты [1, с. 82-95, 3] второго порядка, имеющая пятый порядок на стационарных решениях, также со скалярным стабилизирующим оператором;

В — неявная схема [1, с. 57-59, 4] первого порядка по времени и третьего порядка по пространственным переменным, реализуемая векторными прогонками.

Схемы А и В обладают одинаковыми операторами пространственного дифференцирования и при установлении дают идентичные решения. Как показывают численные эксперименты, частично описываемые ниже, схемы А и Б со скалярными стабилизирующими операторами могут использоваться лишь при числах Куранта, не превышающих 2, схема

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

Искусственная диффузия. Способ построения диффузии проиллюстрируем на примере уравнения Бюргерса

и* + (и2/2)х = 0.

Рассмотрение проведем для дифференциально-разностной схемы

дщ/дг + Д-(/т/2 + ^+1/2) = 0, —го < г < го, (1)

Д = 6 /к, 6 и£ = ±(и± 1 — и).

К ней сводятся явные и неявные схемы разных типов при уменьшении временного шага. Эта схема позволяет исследовать их важные свойства. В уравнении (1) / — энергетически нейтральная компонента потоков, удовлетворяющая для функций с конечным носителем (и£ = 0, если | г |> I) условию

ГО

Е/£+1/2 (и£+1 — и£)/к =0, (2)

— ГО

d — диссипативная их часть, удовлетворяющая для таких функций условию

ГО

У di+l/2(ui+l — и)/к < 0. (3)

—ГО

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

ГО

тТ, ^ о- (4)

г = —ГО

Действительно, умножая разностную схему на щ и суммируя по г, получаем соотношение

ГОГО

д/дЬ ^ («2/2) + ^ «гА (/г+1/2 + ^г+1/2) = 0.

г = —го г = —го

Учитывая тот факт, что мы рассматриваем функции с конечным носителем, для которых оператор разностного дифференцирования А- может быть переброшен на другой сомножитель (с умножением на —1 и заменой индекса на обратный), получаем

ГО ГО ГО

д/дЬ £ («2/2) — £ /¿+1/2А+щ = £ <г.+,/2А+щ < 0.

г = —го г = —го г = —го

Таким образом, последнее неравенство превращается в соотношение (4), так как вычитаемое слева есть нуль. Это соотношение имеет вид априорной оценки с убыванием разностной энергии, которая выполняется для произвольных, в том числе разрывных, начальных условий. Первое условие удовлетворяется, если вычислять /г+1/2, например, по известной формуле второго порядка точности

/г+1/2 = («2 + щ«г+1 + Щ+1)/6-

Поскольку справедлива формула /г+1/2А+щ = (и2 + иг«г+1 + и2+1)(«г+1 — иг)/(6Н) = («3+1 — «3)/(6^), при суммировании по индексу выражения /г+1/2А+щ для функций с конечным носителем получаем нуль. Можно показать, что следующая формула [1, 4], обеспечивающая четвертый порядок аппроксимации конвективного слагаемого уравнения Бюргерса, также удовлетворяет условию энергетической нейтральности (2):

/г+1/2 = («2 + «г«г+1 + «2+1)/6 — [(«г + «г+1)(^г+3/2 — ^г—1/2)/2+

+ (1/2 — С)^г+1/2(^г+3/2 + ^г—1/2) + С (^г+3/2 + ^г2—1/2)]/12, ^г+1/2 = ^+«г

(С — произвольная константа).

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

^г+1/2 = в (2 — а+1 — а)$+иг — Ь5+аг5 5+щ. (5)

Теорема 1 [1, 4]. Если коэффициенты соотношения (5) удовлетворяют условиям

0 < аг < 1, в > 0, Ь > 0,

то для функций с конечным носителем (щ = 0, если | г |> I) обеспечивается свойство (3).

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

ЕШ [1].

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

Щ = «г + (^г+1/2 — ^г—1/2)т/^- (6)

Используя обозначения гг+1/2 = $+иг, р = Ьт/Л,, д = вг/^ и проводя несложные преобразования, получим

£г*+1/2 = гг+1/2[1 — 4д + (аг + аг+1)(2д — 3Р)] + гг+3/2[д(2 — аг+1 — аг+2) + Раг+2 + 3Раг+1] +

+гг— 1/2[д(2 — аг — аг—1) + раг—1 + Зраг] + ¿г+5/2 (—раг+2) + гг—3/2 (— раг— 1) - (7)

Для сохранения монотонности достаточно, чтобы при гг+1/2 > 0 (< 0), —го < г < го имело бы место £г*+1/2 > 0 (< 0). Вычислим коэффициенты аг, например, по формуле

аг = тах[0, ш1п(1, Згг+^/гг—1/2 — 1-5, Згг—1/2М+1/2 — 1-5)]. (8)

Эти коэффициенты удовлетворяют условию 0 < аг < 1. Если величины гг+1/2, гг—1/2 имеют разные знаки или соотношение их модулей превышает 2, то аг = 0. Отсюда, если все гг+1/2 одного знака, то, как несложно показать, выполняются неравенства

2гг+1/2аг — гг— 1/2аг > 0, 2гг—1/2аг — ¿г+^аг > 0, если ^+1/2 > 0, —го < г < го, (9)

2гг+1/2аг — ¿г—1/2аг < 0, 2гг—1/2аг — ¿г+1/2аг < 0, если ^+1/2 < 0 — го < г < го. (10)

Этих неравенств достаточно для доказательства нужной нам оценки для разностной функции гг?+1/2 при выборе определенных значений параметров р, д. Действительно, рассмотрим, например, случай гг+1/2 > 0. Тождественные преобразования позволяют привести формулу (7) к виду

гг*+1/2 = гг+1/2[1 — 4д + Зр(аг + аг+1)] + ¿г+3/2[2д — (д + р)аг+2] +

+гг—1/2[2д — (д + р)аг—1] + [д — 3р]{(2гг+1/2 — ¿г—1/2)аг + (2гг+1/2 — гг+3/2)аг+1)}+

+ (2гг+3/2 — ¿г+5/2)Раг+2 + (2гг—1/2 — ¿г—3/2)Раг—1.

Если выбрать значения параметров д, р такими, чтобы удовлетворялись условия д < 1/4, д > Зр > 0, то выражение в каждой квадратной скобке неотрицательно, а следовательно, первые три слагаемых этой формулы есть комбинация с неотрицательными коэффициентами величин г с разными индексами. Следующие три слагаемых неотрицательны ввиду неравенства (9). Отсюда г*+1/2 также неотрицательно. Таким образом, верна Теорема 2. Если коэффициенты формулы (5) удовлетворяют условиям

вт/й < 1/4, в > ЗЬ > 0

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

Итак, диссипативный поток (5) обеспечивает как сохранение монотонности, так и выполнение свойства диссипации энергии. Аналогичные адаптивные процедуры с диссипацией энергии можно построить на основе диффузионных операторов более высокого порядка, например пятого:

^г+1/2 = в(2 — аг — аг+1)5+щ + Ь5+5 агаг+15+5 5+щг, (11)

где величины аг снова вычисляются по формуле (8). Подставив эту формулу в уравнение (6) и подействовав на него оператором 5+, получаем

¿г*+1/2 = ¿г+1/2 + д[(2 — аг—1 — аг)гг—1/2 — 2(2 — аг — аг+1)гг+1/2+

+ (2 — аг+1 — аг+2)гг+3/2] + р[аг+2аг+3 (гг+7/2 — 2гг+5/2 + ¿г+3/2) —

—4аг+1аг+2(гг+5/2 — 2гг+3/2 + гг+1/2) + 6аг аг+1(гг+3/2 — 2гг+1/2 + ¿г—1/2) —

—4аг—1 аг (гг+1/2 — 2^г—1/2 + ¿г—3/2) + аг—2аг—1(^г—1/2 — 2^г—3/2 + ¿г—5/2^

где р = Ьт/й, д = вт/й. Дальнейшие преобразования позволяют найти следующую форму этого соотношения:

г*+1/2 = ¿г+1/2[1 — 4д — 4раг+1 аг+2 — 12рагаг+1 — 4раг—1аг аг+1/2] +

+^г—1/2[д(2 — аг—1) — 2раг—2аг—1 + 6рагаг+1] + ¿г+3/2[д(2 — аг+2) —

—2раг+2аг+3 + 6рагаг+1] + раг+2аг+3 {^г+7/2 — (0.5 + 1.5)^г+5/2 + 3^г+3/2}+ +4раг+1аг+2(2^г+3/2 — ¿г+5/2) + д(2^г+1/2 — ¿г—1/2)аг + д(2^г+1/2 — ¿г+3/2)аг+1+

+4раг—1аг(2гг—1/2 — ¿г—3/2) + раг—2аг— 1{ 3^г—1/2 — (1.5 + 0.5)£г—3/2 + ¿г—5/2 }.

Если наложить на коэффициенты р, д ограничения 0 < р < д/2, 5р + д < 1/4, то, как легко показать, выражение в каждой из квадратных скобок неотрицательно. Таким образом, первые три слагаемых образуют неотрицательное выражение при неотрицательных значениях ¿г+1/2. Остальные слагаемые неотрицательны ввиду неравенства (9). Таким образом, верна

Теорема 3. Если коэффициенты формулы (11) удовлетворяют условиям

(в + 5Ь)т/й < 1/4, в > 2Ь > 0

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

Естественно, теорема 1 справедлива также для формулы (11).

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

При решении уравнений сжимаемого газа по схемам А и В использовалось обобщение формулы (5), по схеме Б — обобщение формулы (11). При этом в каждом уравнении

параметры а вычисляются по давлению, т. е. искусственная диффузия носит скалярный

характер, что позволяет удобно представлять ее в стабилизирующем операторе. Используется формула более общего вида, нежели (8):

аг = шах[0, шт(адг+, дг—, 1) — с]/(1 — с), (12)

2+ = [(5—Рг)(5+Рг) + е2)]/[(5+Рг)2 + 10—9],

£- = [(5+Рг)(5— Рг) + ер]/[(5—Рг)2 + 10—9],

константы а и с позволяют регулировать свойства диффузии, параметр ер предназначен для предотвращения понижения порядка диффузии на гладких экстремумах или участках почти постоянного решения. По результатам пробных расчетов выбраны следующие значения эмпирических величин: с = 0.4, а = 1.3, ер = 2(Ртах — Рт1П). В формуле (5) использовались параметры в = 0.5, Ь = 1/8, в формуле (11) — в = 0.5, Ь = 1/32.

Тестовые расчеты. Продемонстрируем работоспособность схем с рассмотренной искусственной диффузией на задачах идеального газа. Важными характеристиками схем является глубина и скорость сходимости к стационарным решениям. Для иллюстрации этих качеств рассмотрим результаты расчета двумерной задачи распада разрыва. Прежде всего представим расчеты по схеме Рунге — Кутты третьего порядка, хорошей сходимости которой трудно добиться ввиду ее многослойного характера и скалярного стабилизирующего оператора. Пусть в начальный момент заданы два параллельных полубесконечных потока с параметрами Р =1, р =1, и = 2.4(1.4)1/2 при у > 0 и Р = 0.25, р = 0.5, и = 4(1.4)1/2 при у < 0. Отметим, что при этом ударная волна распространяется вниз, далее имеет место контактный разрыв и выше — веер волн разрежения. Рассматривается область 0 < х < 1, —0.3 — 0.3х < у < 0.3 + 0.3х. Используется сетка из 45 х 60 узлов. На рис. 1 приведена динамика сходимости для схем А и В при разных числах Куранта, показана эволюция невязки плотности Л = шах^ | рП+^/рПс — 1 I /т. Видно, что скорость сходимости для схемы В значительно выше.

Моделирование нестационарных разрывных течений представляет трудную задачу для существующих неявных схем. Схемы А и Б имеют скалярные консервативные стабилизирующие операторы, что позволяет использовать их для этой цели. Рассмотрим результаты решения задачи о взаимодействии вихря с ударной волной. Используется расчетная область [0,1] х [0,1]. Стационарный скачок расположен посередине области нормально к оси х. Параметры слева от скачка принимают значения (р,щ,^,Р) = (1,1.1 ^/7, 0, 1), 7 —

Рис. 1. Динамика сходимости для схем А и В (а и б соответственно).

показатель адиабаты. Небольшой вихрь находится в потоке левее скачка с центром в точке (xc, yc) = (0.25, 0.5). Вихрь задается возмущением скорости (u, v) и температуры (T = Р/р) по отношению к среднему потоку:

u = 1 + e“(1-r2/r2 ^(у — yc)/re, v = — e“(1-r2/r2 ^(x — xc)/re,

(Y _ 1)e2«(1-r2/r2)

T = P/p =1 — (Y 1)e------------£2, P/pY = 1,

7 4«y 7

где r = \/(x — xc)2 + (y — yc)2; £ имеет смысл интенсивности вихря, а управляет скоростью убывания возмущения; re — критический радиус, внутри которого вихрь имеет максимальную интенсивность. В расчетах выбрано £ = 0.3, re = 0.05 и а = 0.204. Используется сетка 251 х 100, которая равномерна по у, но сгущается по x возле скачка. Результаты расчетов (изолинии давления) для разных моментов времени показаны на рис. 2.

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

s = P/pY = 1, T = P/p =1 — £2 Y—1 e1-r2, £ = 5, r2 = x2 + y2,

8Yn2

u = 1 — e(1-r2)/2y/(2n)£, v = 1 + e(1-r2)/2x/(2n)£, —5 < x < 5, —5 < у < 5.

Точное решение заключается в смещении исходных распределений в соответствии со средними значениями компонент скорости u = 1, v = 1. На рис. 3 приведены результаты численного моделирования этой задачи на сетке 81 х 81 для момента времени t = 100, когда

а г = 0.05

Рис. 2. Взаимодействие вихря и ударной волны: сетка 251 х 100, 1.2 < К < 1.4; слева — схема третьего порядка, справа — второго порядка (пятого для стационарных решений).

Рис. 3. Задача о пассивной конвекции вихря. Сечение по центру вихря, плотность, Ь = 100, сетка 81 х 81.

исходный вихрь переместился на 10 размеров расчетной области по каждой переменной. В процессе расчета исходная область периодически сдвигалась вслед за перемещением вихря. На рис. 3 изображены профили плотности в сечении x = 0. В [1, гл. 17] приведены аналогичные данные для схем второго — пятого порядков. Погрешность схем второго порядка значительно больше, расчеты по схемам четвертого и пятого порядков [1] позволяют получить очень хорошее решение.

Сверхкритическое обтекание крыла. Рассматриваются турбулентные течения возле профиля NASA 0012 под углом атаки и возле симметричного профиля, составленного из дуг окружности. На стенке наряду с условиями прилипания предполагается выполненным соотношение Tw = . Остальные граничные условия вполне традиционны и детально

описаны в [5]. Используется алгебраическая модель турбулентности Себечи — Смита:

pt = min(pt1, pt2), pt1 = pl2 | du/dz |, l = 0.4z[1 — exp(-z+/26)],

z+ = z(p | du/dz | /p)0'5, p*2 = 0.0168pUe^a/[1 + 5.5(z/5a)6],

где z — расстояние до стенки; 5a — толщина вытеснения; ue — скорость на внешней границе пограничного слоя; p — ламинарный коэффициент вязкости.

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

На рис. 5 в виде изолиний плотности приведены результаты расчета турбулентного обтекания 15 %-го профиля из дуг окружностей по схемам Б и В. Следует отметить, что при фиксированной на стенке температуре плотность в пограничном слое непостоянна, поэтому картина изолиний плотности отражает пограничный слой. В данном случае его толщина близка толщине линии. Можно также видеть контактный скачок, соответствующий отрыву пограничного слоя от поверхности крыла в точке падения на нее скачка.

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

Рис. 4. Профиль МЛЯЛ 0012 в идеальном газе при угле атаки 1.25°, Мто = 0.8: а — изолинии плотности (221 х 81); б — сетка 111 х 41.

а

Рис. 5. Турбулентное обтекание 15 %-го профиля. Сетка 111 х 81. Мто = 0.9: а и б — схемы третьего и пятого порядков соответственно.

расчетных данных на разных сетках и для разных методов — пятого порядка Б и третьего порядка В. Сравнение рисунков говорит о разумной точности расчетов.

Рис. 6. Турбулентное обтекание 15 %-го профиля. Линии тока. Мто = 0.9. Сетки 56 х 41 (а) и 111 х 81 (б). Схема третьего порядка (сверху); схема пятого порядка (снизу).

На рис. 7 изображены изолинии плотности для турбулентного обтекания 20 %-го профиля из дуг окружности. Структура течения аналогична структуре течения для 15 %-го профиля.

Резюмируя представленные результаты, можно констатировать неплохое разрешение

Рис. 7. Турбулентное обтекание 20%-го профиля. Изолинии плотности. Мте = 0.9. Сетка 111 х 81. Схема третьего порядка.

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

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

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

[2] Пинчуков В. И. Абсолютно устойчивые схемы Рунге — Кутты третьего порядка // Журн. вычисл. математики и мат. физики. 1999. Т. 39, №11. С. 1855-1868.

[3] Пинчуков В. И. Компактная схема шестого порядка для решения уравнений Эйлера в криволинейных координатах // Журн. вычисл. математики и мат. физики. 1998. Т. 38, №10. С. 1717-1720.

[4] Пинчуков В. И. Энтропийные и квазиэнтропийные условия для схем высоких порядков аппроксимации // Вычисл. технологии. 1997. Т. 2, №6. С. 71-87.

[5] КАРАМЫШЕв В. Б., Пинчуков В. И. Турбулентные до- и сверхкритические течения возле крыловых профилей // Изв. СО АН СССР. Сер. техн. наук. 1987. №18, вып. 5.

С. 20-24.

[6] Пинчуков В. И. Об одном методе построения адаптивных разностных сеток в задачах аэродинамики // Моделирование в механике. 1988. Т. 2(19), №2. С. 133-141.

Поступила в редакцию 25 октября 2000 г. в переработанном виде — 18 декабря 2000 г.