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

УЧЕНЫЕ ЗАПИСКИ Ц А Г И

Т о м III 197 2 №6

УДК 533.6.011

ПРИМЕНЕНИЕ ПРИНЦИПА МИНИМАЛЬНЫХ ЗНАЧЕНИЙ ПРОИЗВОДНОЙ К ПОСТРОЕНИЮ КОНЕЧНОРАЗНОСТНЫХ СХЕМ ДЛЯ РАСЧЕТА РАЗРЫВНЫХ РЕШЕНИЙ ГАЗОВОЙ ДИНАМИКИ

В. П. Колган,

Развитие вычислительной техники позволило исследовать с помощью численных методов ряд выдвигаемых практикой многомерных газодинамических задач стационарного или нестационарного обтекания тел. При выборе численной схемы расчета при прочих равных условиях всегда предпочтительнее схемы повышенного порядка точности, позволяющие выявлять наиболее тонкие особенности исследуемого течения и экономить время решения задачи на ЭЦВМ. В общем случае решения нестационарных нелинейных уравнений газовой динамики не аналитичны и могут содержать сильные разрывы в виде ударных волн и контактных разрывов. От этого в немалой степени зависит выбор конечноразностной схемы расчета. Если сильные разрывы выделяются при счете, а в остальной части течение считается гладким, то возможно применение конечноразностных схем, близких ко второму порядку точности аппроксимации [1] —[3]. Однако если число разрывов с течением времени увеличивается и получающаяся волновая картина достаточно сложна, то применение указанных схем становится затруднительным. В этом случае предпочтительнее применение схем сквозного счета [4] — [6], обладающих меньшей точностью, но дающих возможность вести расчет сквозь сильные разрывы без выделения поверхностей разрыва. В настоящее время разрабатываются методы сквозного счета повышенной точности [7] — [9].

Настоящая статья посвящена рассмотрению некоторых вопросов построения конечноразностных схем для расчета разрывных решений газовой динамики с помощью применения принципа минимальных значений производной. Построена разностная схема для модельного уравнения и(-\- их = 0, обладающая вторым порядком точности по координате и первым порядком точности по времени. Исследованы свойства полученной схемы, выведены условия устойчивости счета и условия сохранения монотонности функций. Результаты исследований положены в основу при составлении конечноразностной схемы для решения одномерных нестационарных уравнений газовой динамики. Приведены примеры расчетов распада произвольного разрыва, результаты сравниваются с расчетами, проведенными другими методами.

1. Рассмотрим некоторые общие свойства конечноразностных схем различных типов. Для этого изучим модельное уравнение

и( + их = 0. (1.1)

Это уравнение имеет решение вида

м=/(х-0. (1.2)

где /(х) = и (х, 0) — начальное значение функции при £ = 0. Если решение ищется на конечном отрезке, то граничное условие нужно задавать только на левом конце отрезка, а значения на правом конце отрезка определяются из свойств уравнения (1.1).

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

(|) ийх исН — 0. (1.3)

В большинстве случаев решение задач газовой динамики при использовании конечноразностных схем сводится к следующей процедуре расчета одного шага но времени. Вся область предполагается разбитой на ячейки, которым I риписаны заданные осред-ненные значения массы, импульса и энергии в момент времени t = t0. Для модельного уравнения (1. ) рассмотрим случай равномерной сетки с ша ом Дл- /г. Обозначим значения функции и в центральных точках л =• хп = пк через ип:

и„ — и (пк, (п). (1.4)

Координатам границ ячеек будут соответствовать точки

хп+1/2= + ' к = (х„ + хп+1). (1.5)

Далее строится некоторая аппроксимация поля течения в промежутках между центральными точками ячеек. Аппроксимация служит отправным пунктом дл 1 последующего приближения нахождения решения Че|в1 некоторый промежуток времени X и определения в центральных точках новых значений функции при £ = ^0 + х'

ип = и'пк, + х). (1.6)

Точность аппроксимации определяется порядком точности схемы. Зависимость от координаты для параметров течения между центральными точками ячеек может быть различных типов: кусочнопостоянная [5], линейная, параб лическая и т. д. Имеются численные схемы [10], в которых непрерывное поле течения аппроксимируется дискретными движущимися частицами.

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

— обмен через границы ячеек осуществляется путем переноса дискретными частицами, обладающими своими значениями массы, импульса и энергии;

— течение на границе между ячейками находится в виде ряда Тейлора по времени; при этом временные производные выражаются через известные из аппроксимации пространственные производные (этот случай предполагает в окрестности границы непрерывное гладкое течение);

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

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

<о + х

ДЛ+1/2= ^ u(xn+V2, t)dt, (1.7)

А)

где и(хп+ц2, t) определяется видом аппроксимации и механизма обмена; Дя+1/2 обозначает величину обмена параметром и между п и (п + 1)-й ячейками. Записывая для каждой ячейки интегральные законы сохранения (1.3) в течение времени т, получаем конечноразностную схему того или другого вида. При этом необходимы еще дополнительные ограничения, которые следует наложить на шаг по времени т для обеспечения устойчивости счета.

Отметим справедливость следующего утверждения; любая конечноразностная схема с постоянными коэффициентами для уравнения (1.1) вида

N N

ип= 2 akun-k\ 2 Я*=1'> ак = const (1.8)

k = -N k=-N

может быть представлена в виде

Un = U„ + Fn-l/2 — Fn + 1/2, (1-9)

где функция Fn+i/2 характеризует механизм обмена между га и (га + 1)-й ячейками:

N

Fn+\ti= 2 Ckun+k> (1-Ю)

k=\-N

причем если коэффициенты ак обладают свойством четности

ak = a-k (Л=1, 2, . . . , TV), (1.11)

то коэффициенты ck обладают свойством нечетности:

ci+i + c_t = 0 (г = 0, 1, 2, . . . , TV — 1). (1.12)

И наоборот, при нечетных ak

ak = — a-k {k = 1, 2, . . . , N) (1.13)

коэффициенты ck обладают свойством четности:

ci+i = c-i (i = 0, 1, 2 , . . . , N~ 1). (1.14)

В самом деле, приравнивая (1.8) и (1.9), имеем

N N N

ак itti+k ип ~j~ ск ип^.к—1 ск iiti+k •

k=~N k=l-N k = l—N

N N— 1

— + 2 C* + l &Л + А-

A=l— Л/ k = -N

(1.15)

Это равенство должно быть справедливым при любых наборах ип, т. е. сумма коэффициентов при каждом ип должна обращаться в нуль. Отсюда получаем:

а_дг == С\-„ ;

а\-п = — С\—N + С2—Ы ;

а-1 = — С—1 с0;

а0 = 1 — с0 сй

fli = —- Ci -j- Со',

a.v-i = — слг—i + Слг; = — Cn.

(1.16)

Соотношения (1.16) дают однозначную связь между ак и ск. Отсюда же следуют соотношения о нечетности коэффициентов ск при четных ак вида (1.12): из первого и последнего уравнений (1.16) следует соотношение из (1.12) при i=N—1, из второго и предпоследнего уравнений выводим другое соотношение из (1.12) при i = N — 2 и т. д. до i = 0. Доказательство соотношений о четности ск при нечетных ак проводится аналогично.

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

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

ип ■■

(1.17)

Отсюда имеем выражение для функции Рп+)/2, характеризующей обмен между ячейками:

Га+1П = ип-^- 0-18)

Эта схема для модельного уравнения (1.1) является аналогом схемы, предложенной С. К. Годуновым [5]. Схема устойчива при т//г<1 и является среди всех трехточечных схем с постоянными коэффициентами первого и второго порядка точности наилучшей в том смысле, что сохраняет монотонность функций.

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

2. Примем следующую аппроксимацию всего поля для функции и через значения ип в центральных точках ячеек: будем считать, что внутри каждой ячейки функции распределены по линейному закону, который выбирается с учетом значений функций в центральных точках близлежащих соседних ячеек. Рассмотрим поле внутри п-й ячейки: для *л-1/2<Х-*:л+1/2

Величина градиента для функции и внутри ячейки с точностью до членов 0('<8) не влияет на суммарное значение интеграла по координате в каждой ячейке, входящего в законы сохранения (1.3):

Для построения аппроксимации второго порядка точности по координате величину /г„ можно определять различными способами через значения ип-ь ип и ип+1 . В частности, можно определить величину кп как односторонние производные, выражаемые через значения и в двух точках: .

Из двух значений градиента выбираем минимальное по модулю значение, которое и обеспечит минимальную особенность для градиента функции и внутри п-й ячейки:

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

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

и (х) = ип + /г„ {х - хп).

(2.1)

*/1+1/2

(2.2)

■*0—1/2

^лп|ав при | Д„ | Д*—1 [;

&ллев в противном случае

или

(2.4)

Д„ при |Д„|< Дл_11;

Д„_1 в противном случае.

со второго порядка до первого. Отметим, что в численных методах установления, когда ищется решение при оо, порядок точности метода по времени несуществен.

Для модельного уравнения при таком подходе к рассмотрению распада разрыва функция и принимает на границе п-й и (я + 1)-й ячеек постоянное значение при t^>t0, равное левому предельному значению функции в точке разрыва при t = t0. Допускаемая при этом погрешность равна О (-с).

Применяем к каждой ячейке интегральные законы сохранения в течение времени т. Получаем следующую конечноразностную схему, состоящую из четырех вариантов:

1) и — 11п (Ил_1 Мл+1) — 11п 2Д 1 ~Ь ^я)>

! Дя-21 > | Дл-11 I;

X X

2) и" = и„ +-2Д (— М«-2 + 4 «л-т — 3«л) = И„+2Д-(Д/1-2-ЗДя-1),

| Дя-2 I I 1 I I Д„ I;

3) и» = и„ + -^(ип_1 —ия)==и„ —-^Дв-ь [(2.5)

| Дл-21 1 Дл-11 I |;

4) ип = “л + ^ (— “«-2 + 3 1 ~ ип — “л+О =

= м" ~2К ^п~2 - ^ Дя~1—Дл^’

| Д л—2 I < I Дл-1 I > | |.

Приведенная конечноразностная схема является пятиточечной с переменными коэффициентами. Ошибка аппроксимации исходного дифференциального уравнения (1.1) для гладких решений составляет

£ = О (/г2) + О (х),

т. е. схема обладает вторым порядком точности по координате и первым порядком точности по времени.

3. Одним из требований, предъявляемых к разностным схемам, является их устойчивость. Другое требование [5] заключается в сохранении качественного характера поведения функций, или, другими словами, требование сохранения монотонности функций при счете.

Выведем условия, налагаемые на шаг по времени т, при которых схема (2.5) удовлетворяет этим требованиям. Так как схема имеет переменные коэффициенты, то обычные методы исследования (метод Фурье) не пригодны в данном случае, поэтому доказательство проведем путем исследования неравенств.

Прежде всего докажем следующее свойство схемы (2.5): для всех 1/2 выполняются неравенства

при Д„_1>0 и„-1 <«"<«„; (3.1)

в противном случае и„_1 >-н" > и„.

Доказательство приведем для случая Дл_1>0 (случай Дл_1<0 доказывается аналогично).

Пусть реализуется первый вариант схемы (2.5). Тогда имеем, с одной стороны,

ип — ип = — -~ (Д„_1 + Д„) < 0 при | Д„_21 > | Д/1—1 ] > | Д„ |, но, с другой стороны, получаем

чп — Ип—1 — Д/1—1 ■ 2^ [Д/i—i ~Ь Дл] ^ 1 ^ Д/7—1 ^ О-

Последнее неравенство справедливо при

Пусть теперь реализуется второй случай схемы (2.5). Имеем: ип — ип — ~ (Дл-2 —3 Дл_0 <0 при | Дл_2| < | Дя-11 < | Дл|.

Но ип ограничено и снизу:

ип - Un-x = Дя_1 + -^(Дя-2 - з ДЯ—i) > (1 - Х-) д»-1 > 0.

Последнее неравенство предполагает, что выполняется условие 1/2.

Рассмотрим теперь третий случай схемы (2.5). Имеем: и" —и„= —тДл_i<0 при |ДЛ-8|>|ДЯ_1|<|ДЯ|.

Но при этом выполняется также условие

ип ил_1 = Дл_1 Дл—1 0,

что справедливо при т//г<;1.

Остался не рассмотренным четвертый вариант схемы (2.5). Имеем:

ип — Мл_! = Дл_1 + 2^ [А/г—2 — 2 дл_1 — Д„] < 0 при | Дл_2 К | Дл_11 > | Д„ |, но, с другой стороны, получаем:

ип — Мл-i — Дл—1 + 2^" (Дл-2 — 2 ДЛ_1 — Дл)> Дл—1 ^ 1 — 2 j .

Последнее неравенство справедливо при т/Л<1/2. Свойства (3.1) полностью доказано.

Пользуясь результатами доказанного утверждения, легко показать, что схема устойчива.

Докажем, что если шах |м„|<Л4, то шах | ип\ < шах] ип |<М.

п п п

В самом деле, используя (3.1), получаем:

| ип К шах [ | ип-11, | ип П < шах | ип |,

П

откуда имеем:

max | и” |< мах | ип | < М,

П П

что и требовалось доказать.

Покажем, что схема (2.5) сохраняет монотонность функций при х//г<;1/2. Рассмотрим случай, когда все Дл^-0.{ Тогда из (3.1) получаем:

ил_1<и"<мл; м„<ип+1 <«л+ь

Отсюда сразу получаем, что и"<мл+1.

Аналогично доказывается свойство сохранения монотонности при всех Дл<0. Более детальным исследованием можно показать, что схема (2.5) сохраняет монотонность функций при -с/й-<2/3.

4. Рассмотрим нелинейные уравнения одномерного нестационарного течения идеального газа:

д? д(ри) ді дх ’

д_

ді

д (> и) ді

р р и2

+ (р + рц2) = 0;

д

дх

: 0.

(4.1)

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

Из решения задачи о распаде разрыва находим значения в начальный момент £ = £0 + 0 для величин потоков массы, импульса и энергии, которыми обмениваются ячейки.

После применения интегральных законов сохранения, получаем конечноразностную схему для нахождения новых значений рп, рл, ип в момент времени і = і0 Ь т:

Р” = Рп + [(Я^)/г—1/2 •

(Я^Ол-н/г ];

(ри)« = (р и)п + Т[{Р + 1/2 - (Р + Шг)й+1/2 ];

и

р , р«2У ( р , р«2

Т^т + іт =1т=гт + -т

+ т

и(^+№

т— і

тр . ки2 \

ЗТ + -2-

л + 1/2

л—1/2

(4.2)

Величины Р, Я, С/ соответствуют давлению, плотности и скорости на границе ячеек при £=£0-|-0. Очевидно, что схема (4.2) обладает первым порядком точности по времени и вторым по координате. Можно ожидать, что схема (4.2) будет устойчива при

1

Ф<

2\с\

(4.3)

где с — максимальная скорость распространения возмущений, образовавшихся в результате распадов разрывов. Условие (4.3) взято по аналогии с условием устойчивости для модельного уравнения.

Физически оно означает, что возмущения с границ ячеек за время т не должны достигать центральных точек ячеек. Практические расчеты показали, что условие (4.3) обеспечивает устойчивость счета.

Для примера приведем данные численного расчета задачи о распаде разрыва в газе. Начальные условия задавались в следующем виде при £ = 0:

при *<0, р — 2, р = 2, « = 0; 1

при х > 0, р = 1, р ===== 1, и = 0. |

При £>0 в результате распада разрыва образуются три волны: волна разрежения, движущаяся влево по более плотному газу; контактный разрыв, движущийся вправо и отделяющий плотный

газ от менее плотного, и ударная волна. Результаты расчетов

представлены на фиг. 1 и 2.

—точное решение методом

8 о О •1 р итераций • распределение дадления по координате о расчет по метода С.К-. Го дохода для эалероош координат Г//1

°\

е

\ <• о « в

•-§ ® о О

О

с

1,0 4 о -А-*] 2 8

———------———-----——ШЛ——-------——------- I» «18 81

-го -70 О !0 х/»

Фиг. 1

- точнее решение методом итерации •распределение плотности по координате о расчет по метода С.Х. Году-

в а ° О • • р

о

о координат УП

* г»

• О

1 9 .

** ® п 1

0 с

¥ ( >°о -1а-+

—————————\ьИл——————— и.18<1—

-20 -М О }0 х/Ь

Фиг. 2

Следует отметить более точное описание течения в окрестности разрывов, даваемое схемой (4.2). Приведенные данные соответствуют числу шагов N = 37. При дальнейшем увеличении числа шагов качественный характер поведения решений не меняется.

1. Lax P. D., Wendroff B. Systems of conservation laws. Comm. Pure and Appl. Math., v. 13. No 2, 1960.

2. Бабенко К. И., Воскресенский Г. П. Численный метод расчета пространственного обтекания тел сверхзвуковым потоком газа. „Журн. вычислит, матем. и матем. физ.“, т. 1, Jsit 6, 1961.

3. А л а л ы к и н Г. Б., Годунов С. К., Киреева И. Л., Плинер Л. А. Решение одномерных задач газовой динамики в подвижных сетках. М , „Наука", J970.

4. Lax P. D. Weak solutions of nonlinear hyperbolic equations and their numerical computations. Coinm. Pure and Appl Math., v. 7, No 1, 1954.

5. Годунов С. К. Разностный метод численного расчета разрывных решений уравнений гидродинамики. Матем. сб., т. 3, № 47, 1959.

6. Белоцерковский О. М, Давыдов Ю. М. Нестационарный метод „крупных частиц" для решения задач газовой динамики. „Журн. вычислит, матем. и матем. физ.* т. 11, Н 1, 1971.

7. Русанов В. В. Разностные схемы третьего порядка точности для сквозного счета разрывных решений. Докл. АН СССР, 1968, т. 180, № 6.

8. Бала кин В. Б. О методах типа Рунге — Кутта для газовой динамики. „Журн. вычислит, матем. и матем. физ.“, т. 10, Mi; 6, 1970.

9. Косарев В. И. О расчете сверхзвуковых установившихся течений газа с внутренними скачками уплотнения. „Журн. вычислит, матем. и матем. физ.“, т. 11, № 5, 1971.

10. Harlow F. Н. Hydrodynamic problems involving large fluid distortions. Journ. Assoc. Comput. Machinery, 1957, 4, No 2.

11. Годунов С. К., Забродил А. В,, Прокопов Г. П. Разностная схема для двумерных нестационарных задач газовой динамики и расчет обтекания с отошедшей ударной волной. „Журн. вычислит, матем. и матем. физ.‘, т. 1, № 6, 1961.

Рукопись поступила 10/JII 1972