Из выражения (16.137) видно, что LS имеет форму Ь{ф1г}).. Дискретизацию этого выражения можно провести по формуле

Lr,{фL) = {ФiU2{l+l - - /-1/2(/ - (16.142)

где /+1/2 = 0.5(/+ /+i) и т. д. Индексы п и / в формулах (16.141)) и (16.142) определяют точки сетки в направлении I и г\ соответственно.

Для построения эффективного маршевого алгоритма решения системы (16.141) удобно перейти к вектору зависимых переменных q и провести линеаризацию относительно известного на поверхности i решения. Разложение в ряд Тейлора дает

F4 АА?\ G G + BAq\ § +1 s4MAq \

-2 dF dG dS

(16.143>

Знак над символами означает, что компоненты q определяются на п-м слое, в то время как метрические коэффициенты - на {п+ 1)-м.

При линеаризации Fsi в подслое следует заметить, что Fsi = Fsi (Ч Psi) при w < а (1 + 8). Тогда

FsY = Fi +

L J

(16.144)

L dPsi J

F, + AsiA?4f2, где FAp,{0, 1, 0, uf{ljjy\ a Ap заменяется на Лр, после экстраполяции из области вверх по потоку. Если

u>a{l+es), то Asi = A и Fp = 0.

После подстановки (16.143) и (16.144) в (16.141) получается следующая линейная система уравнений относительно AqS

Asi4-

(LB-LM) А? =

= (i) RHS + - - AF? + (16.145)

Дополнительные члены AAsi и AF? возникают из-за необходимости проведения линеаризации на слоях п и п-{- \. Это обеспечивает консервативность, которая необходима для правильного

сквозного счета скачков (см. гл. 14 и [Schiff, Steger, 1980]). Для подавления высокочастотных осцилляции добавлен диссипатив-ный член четвертого порядка 3)qJ (п. 18.5.1), который имеет вид

= еД1 {ГУ (VA) (ЧГ (16.146)

где для устойчивости 8е < 1/8 и

(VA) ? = ?-2 - Ях + б9/ - + ?/+2- (16.147)

Область вне подслоя также описывается системой уравнений (16.145), в которой только следует заменить Asi на А и положить Fp = 0.

При использовании центральных разностей для аппроксимации Lyi система (16.145) есть (4 X 4)-блочно-трехдиагональная при 2 у JMAX - 1 на каждом слое п -f 1 вниз по потоку. При определении правой части (16.145) используются лишь уже известные значения, расположенные по \ на л-м и {п-1)-м слоях. Граничные условия на твердой поверхности (/=1) и на удаленной (/ = JMAX) позволяют получить соответствующие значения AqJ и AQ/max- решения (16.145) применим метод решения блочно-трехдиагональных систем, описанный в п. 6.2.5.

Шифф и Стегер [Schiff, Steger, 1980] использовали описанный выше алгоритм для расчета двумерных ламинарных и турбулентных течений у параболического профиля. Полученное решение хорошо согласуется с решением, определенным методом установления [Steger, 1978], которое также использовалось для определения начальных данных при х/с = 0.1 (с -длина хорды профиля).

Однако, если шаг маршевой переменной (Л = Ах) сделать слишком маленьким (Ах/с 0.005), маршевый алгоритм расходится. Данное свойство связано со слабой эллиптичностью, обусловленной вязкими членами в (16.38). Расходимость можно устранить, если применить повторяющиеся (итерационные) проходы в маршевом направлении, при этом для psi следует использовать некоторую взвешенную комбинацию значений с предыдущих итераций, т. е. необходимо построить сходящиеся глобальные итерации. Шифф и Стегер получили, что достаточно точное решение получается в результате трех-четырех итераций. Необходимо подчеркнуть, что для достаточно больших чисел Рейнольдса корректное решение получается в результате одного прохода на достаточно грубой сетке, не приводящей к расходимости решения.

( ди . dv\ . , др . V др

где 0) = 1 - (7 - 1)(o/y. Для проведения анализа Фурье (см. (п. 16.1.2)) производные от плотности и температуры заменены

Основываясь на том же алгоритме, Шифф и Стегер [Schiff, Steger, 1980] получили решение трехмерной задачи. На рис. 16.19 и 16.20 представлены распределения давления и профили скорости на сферически затупленном цилиндре, расположенном под углом 5° к набегающему потоку. Маршевое решение получено на интервале 3.07 x/Rn 40.0, где Rn - радиус затупления. На подветренной и наветренной сторонах для давления и продольной скорости имеется хорошее совпадение с расчетами методом установления [Pulliam, Steger, 1980] и экспериментальными данными [Hsieh, 1976]. Шифф и Стегер отмечают, что в методе установления в диапазоне 9.0 ; ?лг 14.0 используется слишком грубая сетка, в результате чего теряется точность решения.

Другой способ расчета дозвукового слоя у поверхности заключается в домножении градиента давления в уравнении импульса в направлении маршевой переменной на некоторый параметр (О, в результате чего устраняется экспоненциальное нарастание решения при одном маршевом проходе. Возможный выбор (О определяется уравнением (16.153). Если и> а, градиент давления учитывается обычным образом. Такой подход использовался для расчета сверхзвуковых вязких течений у крыльев дельтообразной формы в работах [Vigneron et al., 1978; Tannehill et al., 1982] и для расчета сверхзвуковых вязких течений у конусов под углом атаки в работе [Rackich et al, 1982].

Метод Виньерона будет здесь продемонстрирован применительно к укороченным уравнениям Навье - Стокса, описывающим сверхзвуковое вязкое течение, т. е. к уравнениям (16.12),

(16.13), (16.15) и (16.17), в которых лишьр* = р7р(/, и поэтому скорость звука определяется выражением == р/р, а не (16.33). Уравнения могут быть записаны в виде