книги / Течения газа около тупых тел. Метод расчета и анализ течений А. Н. Любимов, В. В. Русанов

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

4.Расчет ударной волны. Расчет формы ударной волны является весьма сущест­

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

Имея в виду сказанное, исследуем структуру разностного уравнения для функции Р и свойства его предельного стационарного решения. Рассмотрим сначала случай плоского или осесимметричного течения, когда зависимость от 0 отсутствует и урав­

Если предположить известной зависимость Х|==1 от ц, I, то из (8.52) получим диффе­ ренциальное уравнение для Р:

(8.53)

которое при 1 -+ оо переходит в стационарное уравнение для предельного решения ^(Л):

Однако рассмотрение (8.53) как независимого уравнения нужно нам здесь лишь для упрощения рассуждений и не имеет принципиального значения.

Из общей разностной схемы, описанной в § 5, вытекает разностная схема для

Для того чтобы (8.55) имело смысл, необходимо определить выражения производных

Если течение несимметрично и к = О соответствует нижней границе области ©

лучим, если в (8.14) положим Р%+® = Р% = Рк и функцию Ф заменим функцией Ф

где Рц%ъ определены формулами (8.56).

Опыт расчетов по схеме (8.55) показал, что решение разностного уравнения (8.57) отличается от точного решения уравнения (8.54) на знакопеременную величину Ск(—1)й, Съ > 0. Величина коэффициентов Скиногда оказывается довольно большой, что приводит к нарушению гладкости волны и к значительным колебаниям газодина­ мических функций во всей области течения.

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

2у = —6/а. Записав для (8.58) разностное уравнение (8.57) с учетом граничного усло­ вия, получим

Нетрудно показать, что при к -> 0 решение разностного уравнения сходится к реше­ нию дифференциального. Однако при каждом конечном к в решении щ присутствует

знакопеременное слагаемое вида ()24* ^ (—1)* О (1 + 2ук)к, величина которо­ го при конечном к может оказаться довольно значительной. Если бы корень был равен —1, то избавиться от этого слагаемого не представляло бы труда с помощью усреднения по индексу к:

(Щ)с? = (м-’к-х + + и*-!)/4- (8-60)

Бели же %2 Ф —1, то усреднение не приводит к полному исчезновению колебаний, а*, лишь несколько гасит их. В случае нелинейного уравнения оно не дает сколько-ни­ будь существенного эффекта, тем более что усреднять приходится не только Р, но и X, что сильно усложняет алгоритм.

Значительно более эффективным является видоизменение разностной схемы та­ ким образом, чтобы корень (в линейном случае) стал равен точно —1. Для этого дос­

таточно включить усреднение в саму разностную схему, а именно в (8.556) вместо Р ^ д) подставлять значения, усредненные по формуле (8.60). Отмечая неусредненные ве­

Для предельного стационарного решения разностных уравнений снова получает­ ся уравнение (8.57) с соответствующими граничными условиями, однако Р%определе­ ны теперь формулами (8.62). Исследуем поведение предельного решения опять на при­ мере уравнения (8.58) с граничным условием и;0 = со. Выражая щ и со^ц через щ,. получим линейную систему:

В отличие от схемы (8.59) усреднение полностью гасит паразитическое решение.

Разностная схема (8.61) имеет стационарное решение вида @2 (—1)к не только в случае уравнения с постоянными коэффициентами, но и в случае линейного урав­ нения вида (8.58) с коэффициентами, зависящими от г\. В самом деле, рааностная схе­ ма выглядит тогда так:

или уравнения на концах интервала [т]0, т)#].

Как было показано, подходящим выбором параметров схемы иногда можно до­ биться, чтобы (32 обратилось в нуль. К сожалению, не удается указать общего способа выбора параметров схемы, обеспечивающего равенство (?2 = 0 для любого линейного уравнения, а тем более для нелинейного уравнения (8.54). Поэтому в общем случае

необходимо проводить усреднение, считая значения У к предварительными, а — окончательными, аппроксимирующими значения функции Р (ц, I) в точках сетки, как и сделано в расчетном алгоритме (см. § 5).

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

нелинейного уравнения Рк вообще нельзя разложить на сумму основного и парази­

тического решений и усреднение лишь означает линейное преобразование Рк в Рк, которые затем объявляются основными. Естественно возникает вопрос о законности такого подхода. Можно ли считать значения Рк основными, исходными для построе­ ния разностной схемы и какой именно схемы?

Наиболее простой и естественный ответ на него состоит в том, что Ркможно считать основными тогда, когда существует однозначное обратное преобразование Рк в Рк.

Всамом деле, тогда достаточно усреднить уравнения (8.61) по к и выразить входящие

вних Рн<через Рк, чтобы получить искомую разностную схему для величин Рк^ ‘\

Существование обратного преобразования Рк в Рк эквивалентно условию одно­

значной разрешимости (8.62) относительно Рк. Чтобы найти это условие, рассмотрим однородную систему, положив Рк = 0, и выясним, когда (8.62) имеет единственное

Отличие от нуля определителя системы (8.65) и есть искомое условие существования обратного преобразования Рк в Рк:

которые удовлетворяют условию (8.66). Рассмотрим сначала симметричный случай, и пусть г\0 = 0. Тогда формула (8.62а), написанная для к = 0 с учетом Р_х =

переходит в ^ 0Н(,)= (^ 1+^ + ^о+(5))/2, т. е. (8.626) при г0 = 1/2. Тогда левая часть (8.66) равна (2г/с — 1)/2 и гк не должно быть равно 1/а- Экспериментальные расчеты пока­ зали, что это действительно так и в предельной форме волны при гк = V* в ряде слу­ чаев возникают нерегулярные колебания, не исчезающие после усреднения. В ре­ зультате численных экспериментов было найдено, что наивыгоднейшим значением является гк — 1.

В несимметричном случае целесообразно выбрать г0 и гк равными единице.

§ 9. Вопросы рсалпсации алгоритма на вычислительной машине

1. Автоматическое управление процессом решения задачи. Алгоритм расчета обтекания тупого тела, арифметическая часть которого описана в § 5, включает в. себя также ряд логических операций контроля и управления. Необходимость автома­ тизации этих операций при решении задачи на ЭВМ очевидна, если иметь в виду вы­ полнение серийных расчетов ряда вариантов. В противном случае расчет каждого ва­ рианта становится самостоятельной сложной задачей, так как требует непрерывного участия квалифицированного специалиста и большого количества календарного вре­ мени, а расчет серии вариантов делается практически невыполнимым. Возможность, автоматизации операций контроля и управления в значительной степени зависит ог возможности их алгоритмизации, в частности от того, насколько четко и полно ука­ заны условия, обеспечивающие корректность процесса решения для класса задач. В этом смысле хорошо обоснованные теоретически методы, к которым принадлежит" метод конечных разностей, имеют преимущество перед менее изученными теоретиче-

значение шага тп. Применение завышенных оценок для А 2 нежелательно,. так как

это приводит к излишнему уменьшению тЛ и увеличению времени счета. Поэтому шаг тп+1 для расчета (п + 1)-го слоя определяется в процессе решения задачи по фор­ муле (6.64), в которой величина Л2 вычисляется на (п — 1)-м слое:

Отнесение Л2 на один слой назад позволяет вычислять Л27~1 одновременно с расчетом 71-го слоя и исключает необходимость специального просмотра тг-го слоя для нахож­ дения максимума А2. Величина у- <1 1 обеспечивает некоторый запас устойчивости и исключает случайное увеличение тп сверх допустимого значения вследствие сме­

щения Л””1 или численных погрешностей. Практически для этого достаточно полагать у- ^ 0.95 -г- 0.96. Эксперимент показывает, что (9.2) в некоторых случаях дает зани­ женные значения тп+1, так что иногда расчет с ух даже несколько большим единицы оказывается сходящимся (см. также стр. 98).

В результате работы блока выбора шага величина тЛ+1 определяется перед нача­ лом счета (п + 1)-го слоя и используется для расчета и остальных величин, зави­ сящих от т.

2) Блок проверки условия на свободной границе. Как показано в § 4, для коррект­ ности решения задачи необходимо, чтобы при ц = Н всегда выполнялось условие (4.34), обеспечивающее пространственный тип этой границы. Нарушение условия мо­ жет возникнуть в процессе решения задачи, например, вследствие расширения дозву­ ковой области. В этом случае необходимо расширить область ©, увеличив Н и пере­ местив границу вниз по течению. В блоке проверки условия (4.34) в процессе счета

ное число, то область © расширяется на один шаг по т], т. е. устанавливается новая граница Н ' = Н + Ат). Величины X и Р экстраполируются на кольцо г\ = Н' ли­ нейно, по последним двум кольцам. Положительное число Д к введено для того, чтобы иметь возможность в окончательном результате отбросить последние кольца, где точ­ ность расчета несколько ниже из-за экстраполяции при счете последнего кольца (см. § 5 п.4), а также чтобы гарантировать выполнение условия (4.34) при случайных ко­ лебаниях функций в процессе установления. В случае необходимости в алгоритм может быть включено сокращение области @, когда разность К — /с * превосходит заданное число.

3) Б лок контроля устойчивости прогонки. Аналогично условию (4.34) в каждой

точке производится вычисление величин и XI, для которых должны быть выполне­ ны неравенства (7.40). На каждом слое вычисляются максимальные и минимальные значения для которых и проверяется (7.40). Так как неравенство (7.40а) практи­ чески никогда не нарушается, то для экономии времени его проверка обычно отклю­ чается. В случае нарушения (7.406) выдается информация об этом и счет может быть либо прекращен, либо продолжен с другой функцией ф (см. § 3), в зависимости от инструкции, введенной в программу.

4) Блок управления процессом установления. Решение задачи методом установ­

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

регулярно проводится сравнение значений Х п и Х п~п°в некоторой норме и проверяется выполнение условия:

где и0> 0 и еуст— заданные величины. Как только (9.4) оказывается выполненным для заданного числа тгх последовательных значений л, т. е. п + 1, п + 2, п + 3, ...

..., п + ^1, процесс установления считается законченным и значение Х п*пгпринимает­ ся за окончательное значение X .

Норма разности определяется следующим образом:

|Г , - Г -"»||= тах {Л^ ;})

(9.5)

N 0 = |и00|~1т а х |# к >,[,

1сг1

где [Х^ ^ — *-я компонента вектора Х ^ х и а б '1 — максимальное представимое в

машине число. Опыт проведения расчетов показал, что практически достаточно брать п0 = 1, уменьшая соответственно е уСт, и проводить проверку условия (9.4) на каж­ дом шаге.

При решении задач методом установления существенным является то, насколько начальные данные отличаются от искомого решения. Точнее, пусть и 3°к — набор параметров, определяющих граничные условия для начального вектора Х° и искомо­

•близки между собой. Однако понятно, что практически такая ситуация далеко не

дить постепенно, вводя зависимость текущих значений 2Рот времени, например пола­

гая З5 (0 = ё (I) 3йн + 11 — 8 (*)]^к, где ё (*н) = 1 и # (*) = О при I > *к.

Режим перехода от 3йн к Зьк естественно назвать режимом вывода на заданные значения 3&к, в процессе которого получается только весьма грубое приближение к искомому решению, после чего производится окончательное установление. Выбор функции $ (г) можно осуществлять различными способами, например полагая в (I) =

тельно крупной сетке содержится несколько тысяч расчетных точек. Так как каждой точке соответствует пять значений компонент вектора X (и, г;, ш, р, р), то общее коли­ чество чисел, описывающих течение в трехмерной области, может достигать нескольких десятков тысяч. Очевидно, что исследование и обозримое представление этой инфор­ мации невозможно без специального алгоритма обработки. Хорошо организованная система обработки дает возможность свести до минимума потери информации, опе­ ративно проводить техническую и методическую отладку алгоритма на ЭВМ, авто­ матически управлять процессом решения и контролировать его точность и, наконец, «с максимальной полнотой и эффективностью изучать физические свойства многомер­

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

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

1) Блок интерполяции. Результат численного решения задачи обтекания пред­ ставляет собой совокупность значений Х т^1 и Р ^гв точках дискретной сетки. Для возможности анализа течения необходимо прежде всего иметь однозначный алгоритм восстановления значений этих функций во всех точках области @. Эта задача ре­ шается в блоке интерполяции, который является основным блоком программы обра­ ботки. Вычисление значений X или Р в любой точке (^, ц, 0) области © производится

ся лишь номера точек сетки, по которым производится интерполяция. Интерполяция

Рвыполняется аналогично. Алгоритм одномерной интерполяции реализован в блоке

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

2)Б лок вычисления координат. В этом блоке производится вычисление цилиндри­

может быть различен и в некоторых случаях может включать в себя несколько ите­ рационных процессов. Следует, однако, отметить, что арифметическая сложность ал­

по заданным аргументам. Число процедур никак не фиксировано и может быть лю­ бым, в зависимости от задания. Совместно с блоками 1) и 2) блок 3) позволяет вычислятьлюбую газодинамическую величину в любой точке течения.

4) Блок вычисления координат изоповерхностей функции. Пусть / = / (X) — некоторая заданная функция газодинамических величин и /* — некоторое ее значе­ ние. Совокупность точек (^, т], 0), в которых / (X) = /*, будем называть изоповерхностыо функции /. Исследование изоповерхностей и их сечений различными плоскос­ тями (изолинии) существенно облегчает и упрощает качественное и количественноеизучение течений. В рассматриваемом блоке проводится вычисление координат точек пересечения заданных изоповерхностей с координатными линиями г\ и 0. Для этого-

чения определяются грубо, путем обратной линейной интерполяции по ^ на нуле­ вые значения величины / (ХтМ) — /* при фиксированных к и I, а затем уточняются с помощью итерационного процесса типа метода хорд.

Блок 4) состоит в основном из логических и управляющих операторов, использую­ щих в качестве арифметических подпрограмм первые три блока.

5)Блок перехода к новой системе. Предполагается, что переход от (|, т), 0) к какой-нибудь другой системе координат (дг, д2, д3) описывается аналитическими фор­ мулами. Требуется вычислить значения X, Р в точках равномерной сетки в координа­

тах дъ #2, #3. В блоке 5) вычисляются координаты (!*, ц, 0) этих точек, после чего ра­ ботает блок 1). В частном случае координатами д19д21 д3 могут быть цилиндрические или декартовы координаты.

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

мощью интерполяционного многочлена второй или третьей степени (от параметра з).

7) Блок интегрирования системы обыкновенных уравнений. Система обыкновенных уравнений вида 62/63= /^ дг/дз= /2, ду/дз = /3, где Д — функции 2, г,фи X , интегри­ руется методом Рунге — Кутта или методом Эйлера второго порядка. Этот блок применяется для построения линий тока, характеристик и других кривых, опреде­ ляемых дифференциальными уравнениями. Для вычисления значений правых час­

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

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