CC BY
29
УЧЕНЫЕ ЗАПИСКИ Ц АГ И
Том XIX 19 88 Мб
УДК 533.6.071.4 532.525
ПРИМЕНЕНИЕ КОМПАКТНЫХ СХЕМ ТРЕТЬЕГО-ЧЕТВЕРТОГО ПОРЯДКА ДЛЯ РАСЧЕТА ТЕЧЕНИЯ ГАЗА В СОПЛАХ С БОЛЬШИМИ СВЕРХЗВУКОВЫМИ ЧИСЛАМИ М НА ОСНОВЕ УПРОЩЕННЫХ УРАВНЕНИЙ НАВЬЕ —СТОКСА
А. П. Быркин, Т. А. Тимофеева, А. И. Толстых
Реализован численный алгоритм решения упрощенных уравнений Навье —■ Стокса для сжимаемого газа, основанный на применении компактных схем третьего-четвертого порядка аппроксимации.
С использованием данного алгоритма проведено численное моделирование течения вязкого газа в соплах аэродинамических труб. Получено удовлетворительное согласие результатов расчета с экспериментальными данными.
Для эффективного численного исследования внутренних течений вязкого газа существенным является рациональный выбор численной модели, т. е. исходной постановки задачи и алгоритма ее решения. В этой связи для расчета стационарных течений при наличии сверхзвукового ядра перспективным оказывается использование упрощенных уравнений Навье—Стокса, допускающих применение маршевых алгоритмов с переходом от одного сечения, перпендикулярного оси (или плоскости) симметрии, к другому [1—-3]. Существенную роль при построении численной модели играет выбор аппроксимации исходных уравнений, от которого зависит не только качество получаемых решений, но и в некоторых случаях сама возможность их получения.
В данной работе используется схема, основанная на компактных аппроксимациях третьего порядка [4—6]. Их применение позволяет вести устойчивый счет без введения каких-либо искусственных диссипативных механизмов. Получаемые численные данные при этом отличаются достаточно высокой точностью и практически не искажаются схемной немонотонностью.
Отметим, что на основе упрощенных уравнений Навье—Стокса решению обратной задачи о нахождении контура сверхзвукового сопла посвящены работы [7, 8], решению прямей задачи — работа [9].
1. Не приводя полную систему уравнений Навье—Стокса, упрощенную систему запишем для интересующих нас случаев течения газа в каналах и соплах, обладающих плоскостью или осью симметрии. При этом перейдем от декартовой (или цилиндрической) системы ко-
ординат х, у (ось х направлена вдоль оси канала) к системе координат 1=х, г\ = у/ут{х), где Ую(х) —заданная форма контура сопла. Сохраняя составляющие скоростей и и V вдоль осей х и у, используемые упрощенные уравнения Навье—Стокса получим из полных уравнений при отбрасывании в них вязких членов со вторыми производными по 1 и смешанными производными по | и г). Полученные таким образом уравнения содержат полностью все члены уравнений Эйлера и вязкие члены С Производными ПО Т1 и имеют вид
dF , до
д£ d-q
Н
Re
Т,
(1)
где F, G, Н, Т—векторы-столбцы:
Г=УУг
р и р(®
р иг+р р(®
puv ) G=j/V Р (v
(, , ф + V2 \ р»(4+ 2 ) 0 P(v Л
//= V 0 , т= т2
Р п
0 Т<
■у. W)
■y'w)h +
и2 + V3
В уравнениях (1) все величины безразмерные; р — плотность, р— давление, Н — энтальпия, р, — динамическая вязкость, Рг,— число Прандтля, Яе = р1и1уго х/уц (индекс 1 отвечает значениям размерных величин на оси сопла в критическом сечении, индекс да — на стенке); значение v = 0 соответствует плоскому, л>=1—осесимметричному случаю; уш—производная контура ую(х)\ Т1 = 0, составляющие Г2, Т3, Т.± могут быть получены, например, так же, как и в работе [10].
Кроме того, газ предполагается совершенным:
х — 1 и
р=——рА
(х — показатель адиабаты) и считается заданной зависимость |л(/г). Граничные условия для системы (1) следующие:
р = р(т,), к = и(т)), г> = г>(т)), /г = А(т)) при 1 — 0,
— -4^- = -4^- = 0, v = 0
дт\ df| av) *
u — v = 0, h = hw (£)
при v) = 0,
ПрИ 7|=1.
(2)
Для численного решения системы (1) с граничными условиями (2) будем применять итерационный процесс, при котором в течение каждой итерации используется маршевый алгоритм, а аппроксимация продольной составляющей градиента давления строится с учетом найденных величин давления вниз по потоку от рассматриваемого сечения |=const. Применение такого подхода, часто называемого методом глобальных итераций, сводится к тому, что при Мж<1 + е (Мх—число М в расчетной точке, определенное по продольной составляющей ско-
рости, є>0) член д/д^УфУ р) в уравнении импульсов в направлении
І, ответственный за передачу воз-мущений вверх по потоку, представляется в разностном виде следующим образом:
где I, / — номера узловых точек расчетной области в направлениях | ит), п — номер текущей итерации.
В схеме глобальных итераций необходимо, кроме того, задание граничного условия для давления на выходе из расчетной области.
При дискретизации уравнений (1) предварительно осуществлялось сгущение узлов в пристенной области, достигавшееся преобразованием переменной г) и введением разностной сетки с постоянными шагами. Поясним конструкцию разностного алгоритма для исходной системы (1), записанной в: виде
где вектор U=(p, и, v, Л), L = -^j, Q = -щ- — якобиевы матрицы.
Для гиперболической части системы (4) при Мх > 1 на сетке шй: т]j = jh, Л = const введем компактные „трехточечные“ операторы ВГ1 и Сп [5], основанные на диагонализации матрицы Z.-1 Q:
Матрица М имеет вид /,-1, где столбцами матрицы 5
являются собственные векторы матрицы £-1 <2, а О —диагональная матрица, образованная величинами signXг, где Хг — собственные значения матрицы I-1 С?.
Вводя на сетке шг: ?, = « (т — шаг по I) простейшую дискретизацию производной дР)д%, запишем разностный алгоритм для системы (1) в виде
где 0,5 < о < 1 — весовой множитель, Тн—сеточная аппроксимация вектора Т.
В данном случае для функций / = и, V, к использовалась аппроксимация вторых производных вида
приводящая к погрешности схемы (5) порядка 0(й34-А/Ие) в сверхзвуковой области.
— (УшУ'Р)
д%
і. /
п
(Ув. УvР)"+1, і — (Уч. У" Р)" j AS
(3)
Re
т,
(4)
dF
дО
Bnfj = (A9 - 0,25Д Сп = 0,5k'1 [&0 — (ffi+l/aA+ - M;-i72A _)\f},
A0fj = + -f-fj + ~-//+1, До fj=fJ+1 - /у_,,
Д -/;• = fj - /,•_! . Д+/у = //+, -/; , ^/=Fl/2 = (Лїу + Жт)/2.
+ С, [oG, + (1 - a) G,-i] = Я, [а ( Д, + -і- Ті) +
(5)
Элементы матрицы М предполагаются вычисленными по известным значениям функций при |=|г-1. Граничные условия для схемы (5) основаны на аппроксимации условий (2) и одном дополнительном условии для значения плотности (или давления) на стенке. Варианты последнего подробно рассмотрены в работе [6].
При больших числах Re члены порядка 0(h) в этой схеме можно рассматривать как малую поправку к коэффициенту вязкости ц,. Как показали методические расчеты, они не влияют на выходную точность метода.
Используя результаты работы [5], можно показать, что в приближении «замороженных» коэффициентов схема (5) является абсолютно устойчивой. На практике она обеспечивает возможность счета при больших значениях чисел Куранта.
В дозвуковой области (точнее, при Мх<1 + е) применялась «центрированная» дискретизация первых производных по г] с погрешностью порядка 0(/i4). Она получается при использовании операторов Вч и С,,, в которых формально положено М= 0, и локальная погрешность аппроксимации при этом имеет вид 0(/i4+A2/Re).
Использование центрированной («симметричной») аппроксимации может приводить к схемным осцилляциям решения (в данном случае плотности). Оказалось, что выбором весового множителя а (изменяющегося от 0,5 при Мж>1 + е до 0,8 вблизи стенки) эти осцилляции можно устранить. Отметим, что такой прием применялся в работе '[3] при использовании центрированной схемы второго порядка.
Реализацию схемы (5) можно свести к решению системы разностных уравнений с блок-трехдиагональной матрицей
Aj Uи+ Bj Ut. i + Cj Uit /+1 = Dj, (6)
где Aj, Bj, — матрицы 4X4, Ut. /— вектор искомого решения, Dj — вектор правых частей. Структура матриц А, В, С и вектора D различна для сверхзвуковой и дозвуковой областей течения: в первом случае в них входит матрица М, во втором М = 0, но учитывается представление (3) для величины dld%(ywy>р). Нелинейные относительно компонент U члены в схеме (5) линеаризируются по методу Ньютона. С учетом граничных условий на оси и стенке сопла система (6) решается методом векторной прогонки.
2. По описанному выше алгоритму проведены методические тестовые расчеты.
В качестве первой тестовой задачи рассмотрена задача о сверхзвуковом радиальном течении невязкого газа (и=1,4; Re = oo) в коническом (v=l) сопле, jy^,=0,25. Значение числа М в начальном сечении на оси было задано равным Moi = 3. Расчет, проведенный до значения Мо=20 по схеме третьего порядка с равномерным шагом по г) при N = 40 (N — число узлов в поперечном направлении) и неравномерным шагом по т, показал, что имеет место практически полное совпадение точного и численного решений. Вместе с тем попытка применения симметричной схемы четвертого порядка точности не увенчалась успехом из-за возникающей большой схемной немонотонности.
Другой тестовой задачей была задача об автомодельном течении вязкого газа с радиальными линиями тока в плоском (v = 0) расширяющемся канале (рис. 1) при числе Мо=3 и числе Re = 500 (х=1,4; ц = const; Рг = 0,71; yw =0,23295) [11].
у
X
Точное решение уравнений Навье—Стокса [11], являющееся также точным решением системы (1), обладает тем свойством, что продольная и поперечная скорости и температура (энтальпия) зависят только от г), а плотность, зависящая от т], изменяется по гиперболическому закону от | (р = р(ті)Д). Для такого течения во всем поле от оси до стенки выполняется равенство д/дЪ,{уюр) =0.
По этой причине оказалось достаточным выполнить лишь одну глобальную итерацию, е = 0,1. Расчеты проводились на неравномерной сетке по г) (Л/'=40), полученной в результате введения новой поперечной координаты
С=1-1п Ї1 ~ ^ а =10 (7)
^ 1п(1+а) » “ Ш V'
с последующим построением равномерной сетки по переменной £.
На рис. 1 представлены результаты сопоставления точного и численного решений при х = 3 (х= \ отвечает начальному сечению), указывающие на их хорошее соответствие.
Благоприятные свойства схемы, выявленные при решении тестовых задач, создали предпосылки для проведения расчетов течения газа в конкретных соплах.
Так, было осуществлено численное моделирование течения вязкого газа в профилированном осесимметричном сверхзвуковом сопле М = 6 [12]. Координаты контура сопла у№(х) получены в результате ре-
шения обратной задачи, исходя из требования, чтобы в выходном сечении число М однородного ядра потока равнялось М = 6 (г/ю1 = 2,5 мм, рабочий газ—воздух, давление торможения ро=5-103 Па, температура торможения 70=450 К, Г„,==290 К, Ре=1098). В качестве поправки идеального контура сопла у(х) на влияние вязкости принималась толщина вытеснения пограничного слоя б*. Пограничный слой на стенке предполагался ламинарным.
На рис. 2 представлены контуры сопла у{х), уи,(х) и характеристический ромб. Видно, что поправка на вязкость в большей части области течения соизмерима с радиусом невязкого контура.
Проведенное в работе [12] экспериментальное исследование описанного сопла показало, что в нем реализуется течение, близкое к расчетному.
При решении прямой задачи для найденного контура сопла ую(х) на основе схемы (5) принималось, что в критическом сечении пограничный слой отсутствует и имеет место однородный поток газа с числом М= 1,01. Расчеты осуществлялись при ц = кп (« = 0,76), Рг = 0,71. Значение плотности на стенке находилось с использованием экстраполяционного условия.
Кроме начальной (нулевой) итерации, когда при Мж<1 + е (є = 0,1) полагалось, что д/ді (ут у" р) =0, были сделаны две итерации. Итерации проводились с граничным условием для р на выходе из сопла д2р/д^2 = 0. Результаты расчетов для первой и второй итерации оказались практически совпадающими между собой и мало отличающимися от результатові расчета для нулевой итерации.
• На рис. 2 приведены расчетные данные в виде распределения чисел М по оси — М0(х) и в поперечном сечении — М(у) при 'х=60,85, сопоставленные с данными эксперимента [12]. Из этого сопоставления можно сделать вывод об удовлетворительном их согласии. В качестве
Сопло М=6
Рис. 2
3—Учены'; записки .Vі 6
33
примера на рисунке приведены профили газодинамических величин и, V, И, и р ПО Г) (я = 64,85, сплошные линии), по которым можно судить о внутренней структуре течения. В частности, видно, что имеет место практически равномерное сверхзвуковое ядро течения. На рис. 2, кроме того, показаны в масштабе профили и (у) в поперечных сечениях сопла; они характеризуют размеры невязкого ядра течения.
Была решена также задача о течении газа в рассматриваемом сопле в случае заданного давления рс на его срезе; оно полагалось равным 0,3 -10-3 и составляло менее половины того значения, которое было получено с использованием «мягкого» граничного условия д2р/д^=0. Для получения решения с заданной точностью в этом случае потребовалось около десяти глобальных итераций. На рис. 2 приведены продольные распределения давления на стенке рт(х) в концевой части сопла, соответствующие обоим граничным условиям для р: <32р/д|2 = 0 (сплошная линия) и задание давления рс на срезе сопла (штрихпунктирная линия). Там же представлены также распределения газодинамических величин по г]. Из рис. 2 следует, что влияние передачи возмущений вверх по потоку, несмотря на толстый пограничный слой, проявилось в основном на характеристиках течения в пограничном слое. Число М и однородность в ядре течения на выходе из сопла практически не изменились. Поскольку для практики основной интерес представляет структура ядра течения, все дальнейшие расчеты течений в соплах проводились с граничным условием для давления д2р/д%2=0.
В качестве следующего примера применения разработанного маршевого алгоритма приведем результаты расчета течения воздуха в осесимметричном профилированном сопле, предназначенном для получения в рабочей части ударной трубы потока с числом М= 16 [13].
о эксперимент
Контуры у (х), уго(х) и характеристический ромб представлены на рис. 3. Контур Угс(х) определялся так же, как и в предыдущем случае; исходные данные при этом были следующие: давление ро=50- 105 Па, температура 7'0= 1650 К, температура стенки сопла ввиду малой продолжительности режима работы установки (несколько сотых секунды) принималась равной комнатной (7’№ = 290 К}, значение ут 1=1,2 мм. При расчете невязкого сверхзвукового контура сопла у(х) воздух рассматривался как двухатомный газ (х=1,4), что соответствует замороженным колебательным степеням свободы молекул газа. Последнее на основании данных работы 4] близко к действительности.
Ввиду сравнительно небольших значений чисел Рейнольдса Ие*, вычисленных по параметрам на внешней границе пограничного слоя и текущей координате х, а также больших сверхзвуковых значений числа М, можно предполагать ламинарный режим рассматриваемого течения.
При проведении расчетов течения вязкого газа воздух также рассматривался как двухатомный газ, зависимость |я(Г) задавалась по формуле Сезерленда, Рг = 0,71; ро=45-105 Па, То = 1650 К, Ую 1=1,2 мм (Ие = 0,95 • 105), Г„ = 290 К.
На рис. 3 полученная в результате расчета зависимость М0(я) и профиль М (у) в выходном сечении сопла сопоставляются с результатами эксперимента рЗ]. Экспериментальные значения чисел М при этом определялись на основании измерений давления р'0 приемником полного давления, как для совершенного газа с показателем х=1,4. Видно, что имеет место удовлетворительное согласие расчетных и экспериментальных данных. Реализованное на выходе значение числа М невязкого ядра течения равно примерно 15,5, что примерно на 3% ниже проектного значения М=16.
Отметим, что в отличие от результатов расчетов [9] течения газа в этом сопле при идентичных граничных условиях распределение М0(л:), представленное на рис. 3, после стадии разгона является почти горизонтальной линией. Размер невязкого ядра течения при этом практически полностью совпадает с экспериментальным. Профили н(т]), у(г]), /1(14) и р(т)), приведенные на рис. 3, соответствуют сечению в окрестности среза сопла.
На рис. 4 приведены результаты расчета течения гелия (х = 5/3) в коническом гиперзвуковом сопле, сопоставленные с данными эксперимента, приведенными в работе [15]. Расчетные параметры торможения газа были следующими: р0 = 8- 105 Па, 7о=300 К. Радиус критического сечения сопла ут 1=1,36 мм, полуугол раствора стенок сопла 0«) = 60, а зависимость динамической вязкости гелия от температуры имеет вид (х ■—■ 7’0’647; число Рг принято равным 0,68.
Для заданных параметров торможения 1?е = 0,585 • 105; при этом число Кеж, вычисленное по параметрам течения в невязком ядре на выходе (на основании данных экспериментов) и длине сопла, примерно равно 1,4-106. При таком числе 1?еж можно предполагать ламинарный режим течения газа в пограничном слое. В данном случае, как и при расчетах [15], проведенных в предположении разделения потока на пограничный слой и невязкое ядро, температура стенки Тю принималась близкой к теплоизолированной (Г№~0,83 Т0).
Контур сопла Уго{х) в окрестности критического сечения представляет собой дугу окружности радиуса Я=10 • у™ и сопрягающуюся с прямой линией в точке, соответствующей значению 0Ц,=6°.
Ионическое сопло
Рис. 4
На рис. 4 приведены профили газодинамических величин по т] на выходе из сопла; в распределении р(ц) пр'И этом виден небольшой излом, возникший, видимо, вследствие резкой перестройки структуры течения в поперечном направлении. Из представленных на рис. 4 результатов расчета и эксперимента видно, что между ними имеется удовлетворительное согласие.
На рис. 5 представлены результаты расчета течения воздуха в осесимметричном профилированном сопле с числом М=18. Контур данного сопла уц,(х) получен по описанной выше методике. Он должен приближенно обеспечить при расчетных условиях (у№ 1=1 мм, Ро = = 1 • 107 Па, То=2000 К, Тю=373 К, Ие= 1,56• 105) в выходном сечении однородное невязкое ядро течения с числом М=18.
Пограничный слой предполагался ламинарным и рассчитывался по методу [17]. При этом расчеты пограничного слоя и невязкого контура у(х) проводились для значения х= 1,4, с тем чтобы исключить влияние эффектов неравновесности и основное внимание уделить влиянию эффектов вязкости.
Расчет контура у(х), проведенный методом характеристик [16] при заданном распределении М0(х) на разгонном участке, выполнен
В. П. Верховским. Разгонный участок включает конический участок с полууглом раствора 0,0= 12°. Контур у(х), характеристический ромб и распределения Мо(х), М(дс), соответствующие оси и стенке сопла без учета вязкости, показаны на рисунке штриховыми линиями. После определения контура ут(х) (обратная задача) при помощи алгоритма (5) была решена прямая задача.
Интересно отметить, что контрольный расчет по схеме (5) течения невязкого газа в сопле с контуром у(х) привел к практически пол-
Рис. 5
ному совпадению распределений чисел М(х) и М0(х), соответствующих прямой и обратной задачам.
Расчеты течения газа в. сопле на основе упрощенных уравнений Навье—Стокса проводились при тех же параметрах торможения, температурных условиях на стенке и теплофизических свойствах, что и расчет пограничного слоя. Значение Re=l,56 - 105, соответствующее указанным выше условиям, назовем Rep. Как в данном, так и в рассмотренных выше случаях значения функций у'т (х) и yw(x) в каждой точке по х определялись на основе сплайновой аппроксимации функции y'w (х), заданной не менее, чем в 100 табличных точках.
Из приведенного на рис. 5 распределения М.0(х) для значения Re = Rep видно, что оно сравнительно близко к распределению Nlo(x), отвечающему обратной задаче (штриховая линия), хотя на основной части сопла реализуемое число М не достигает проектного значения 18 примерно на 3%. Получаемый при этом продольный градиент dM.0/dx незначителен. На рисунке приведены, кроме того, профили и(у) в разных сечениях сопла, профили величин и, v, h, р по rj и профиль М (у) вблизи выходного сечения, характеризующий размеры невязкого ядра потока на выходе из сопла.
С целью выявления нерасчетных условий на характеристики течения в рассматриваемом сопле численные данные были получены также при Re = ReP/2 и Re = 2Rep; они представлены на рис. 5. Можно видеть, что равномерность реализуемого ядра течения на выходе из сопла заметно ухудшается.
На основе проведенных расчетов можно заключить, что использование при проектировании профилированных сопл с большими сверхзвуковыми числами М в качестве поправки на влияние вязкости толщины вытеснения пограничного слоя обеспечивает вполне приемлемые результаты даже в случае толстых пограничных слоев с поперечными размерами, сравнимыми с размерами невязкого ядра течения. Этот факт представляется удивительным.
Все представленные выше результаты получены на неравномерной по г] сетке с параметром а в преобразовании (7), меняющимся от 10 до 300, и числом узлов, равным 65; величина е = 0,1. Шаг в продольном направлении изменялся от 0,05 в окрестности критического сечения до 12,5 в концевой области сопла. При этих параметрах сетки расход газа сохранялся по длине сопла с точностью до долей процента. Следует отметить, что с такой же точностью сохранялась в невязком ядре течения и энтропийная функция pjр*, хотя никакого выделения невязкой части течения в процессе решения не производилось.
Высокая точность полученных решений, отсутствие проблемы подавления схемных осцилляций и надежность алгоритма создают благоприятные предпосылки для его применения и в других задачах о внутренних течениях газа.
ЛИТЕРАТУРА
1. Lin Т. С., Rubin S. G. A numerual for supersonic viscous flow over a slender reentry vehicle. —A1AA Paper, 1979, (N 205.
2. К о вен я В. М., Черный С. Г. Решение упрощенных уравнений вязкого газа маршевым методом. — Числ. методы мех. сплошной среды,
1979, т. 10, № 1.
3. Вой но вич П. А., Фу рее нк о А. А. Расчет струйных и внутренних течений вязкого газа. — ЛФТИ им. А. Ф. Иоффе АН СССР, Ленинград, Препринт № 860, 1983.
4. Толстых А. И. О разностных схемах повышенной точности для численного решения некоторых задач аэродинамики. — Ученые записки ЦАГИ, 1973, т. 4, № 2.
5. Толстых А. И. О неявных схемах повышенной точности для систем уравнений. — Ж- вычисл. матем. и матем. физ., 1981, т. 21, № 2.
6. Савельев А. Д., Толстых А. И. Алгоритмы расчета течений вязкого газа, основанные на компактных аппроксимациях третьего порядка.— Ж- вычисл. матем. и матем. физ., 1987, т. 27, № 11.
7. М и х а й л о в В. В. Метод расчета сверхзвуковых сопл с учетом влияния вязкости. — Изв. АН СССР МЖГ, 1969, № 1.
8. Денисенко О. В. Метод расчета сверхзвуковых сопл при сильном влиянии вязкости. — Ученые записки ЦАГИ, 1982, т. 13, № 4.
9. М у ч н а я М. И. Исследование течений в гиперзвуковых соплах в рамках упрощенных уравнений Иавье — Стокса. — Изв. АН СССР, МЖГ, 1986, № 6.
10. Асланов Т. Д., Бы р кин А. П., Щенников В. В. Численный расчет внутренних течений вязкого газа с использованием уравнений Навье—Стокса. — Ученые записки ЦАГИ, 1981, т. 12, № 3.
11. Быркин А. П. Об одном точном решении уравнений Навье — Стокса для сжимаемого газа.—■ ПММ, 1969, т. 33, № 1.
12. К У Д р я в ц е в а Л. И., М е ж и р о в И. И., П о н о м а р е в С. П., Якушева В. Л. Экспериментальное исследование осесимметричных профилированных сопл при малых числах Ие. — Ученые записки ЦАГИ, 1973. т. 4, № 3.
13. Межиров И. И. Исследование течений в гиперзвуковых соплах аэродинамических труб. — Труды ЦАГИ, вып. 2119, 1981.
14. К о м а р о в В. Н., Полянский О. Ю. Методика определения параметров неравновесного потока воздуха в гиперзвуковых аэродинамических установках по данным измерения Г0, р0 и р0 (7’0<2000К). — Ученые записки ЦАГИ, 1987, т. 18, № 5.
15. Б ы р к и н А. П., Межиров И. И. О расчете течения газа в ги-перзвуковом сопле с учетом влияния вязкости (прямая задача). — Ученые записки ЦАГИ, 1971, т. 2, № 1.
16. Кацкова О. Н., Наумова И. Н., Шмыглевский Ю. Д.. Ш у л и ш н и н а И. П. Опыт расчета плоских и осесимметричных сверхзвуковых течений методом характеристик. — М.: ВЦ АН СССР, 1961.
17. Быркин А. П., Щенников В. В. Об одном численном методе расчета ламинарного пограничного слоя. Ж. вычисл. матем. и матем. физ., 1970, т. 10, № 1.
Рукопись поступила 11/]/ 1987
