Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
Метод Годунова. Для расчета разрывных решений широко используется метод, основанный на решении задачи о распаде разрыва. Пусть начальные данные являются кусочно-постоянными на некоторой сетке {хт}"=0, Т.е. и, Р, е(0, х) = {и°т+1/2, рв+1/2, е°т+1/2} при jc Є (Jcfn, хт+1). Оказывается, эта задача имеет точное явное решение. Оно строится так: в каждой точке Jcm нужно решить задачу о распаде разрыва независимо от всех остальных разрывов. Такое решение можно использовать до того момента времени, когда при каком-то значении т правая волна, образовавшаяся от разрыва в точке Jcfn, встретится с левой волной, идущей от распада в точке Jcm+1.
Перейдем к описанию схемы. Основными счетными величинами являются ит +1/2, р"+]/2, е" +1/2, где п — номер временного слоя. Величина unm+il2, например, представляет собой приближенное значение функции u(in, (jcm + jcm+, )/2). Пусть начальные данные
ит+1/2> Pm+1/2> ет+1/2 превращены в кусочно-постоянные на сетке {jcJ функции. Построим точное решение уравнений газовой динамики и определим малое (порядка min (jcm+1 — Xm)) время т, в течение которого это явное решение существует.
Один шаг численного интегрирования соответствует продвижению по t на х. Для того чтобы сделать второй шаг, нужно и при t = х иметь кусочно-постоянное решение. Точное решение, конечно, таковым не является. Поэтому оно аппроксимируется кусочно-постоянным с сохранением в пределах каждого интервала основных физических величин: массы, импульса и энергии. Например,
xM+1
(*m + I _ ^m)Pm+1/2= \ Р(Т> *) d*•
. *»
Аналогично вычисляются
Pm+1/2 ит+1/2> Рт+1/2Іет+1/2 +
Однако и, р, е при Z = X являются сложными функциями JC, вычислять интегралы от них трудно. Это препятствие обходится следующим образом.
Рассмотрим ячейку (хт, *т+() х (0, t) и запишем уравнения в интегральной форме (3), взяв ячейку в качестве ?2. В результате
ОДНОМЕРНЫЕ УРАВНЕНИЯ ГАЗОВОЙ ДИНАМИКИ
299
получим соотношение
Хт+\ ' X X-n+l X
J Л[0, х] dx — J Q[t, xm+l] dt - J Л[х, x]dx+\ Q[t, xm) dt = 0,
0 (17)
где R[t, x] = R(u(t, x), p(t, x), e(t, x)). B (17) третий интеграл есть то, что нам требуется. Первый интеграл вычисляется элементарно, так как Л[0, х] = R(u°m+1/2, р® +1/2, +|/2) в силу постоянст-
ва начальных данных на (xm, Jcm+,).
Так же легко вычисляются второй и четвертый интегралы. В силу атомодельности решения задачи о распаде произвольного разрыва на ЛИНИИ X = Xm (т.е. I = const) все функции постоянны. Обозначим их Pm, ет. Тогда второй интеграл есть Q(um, рт, ет)х.
Формально схему Годунова можно записать в виде
^С+\/2~ Rm+l/2 , Qm+\~~Qm __ q
Т Xm+l Хт
Следует иметь в виду, что «поток» Qm вычисляется решением задачи о распаде разрыва. Она сводится к решению системы нелинейных уравнений. Это относительно «дорогая» операция (ведь она проводится при всех т, п). Значительные усилия прилагаются к тому, чтобы снизить ее трудоемкость. В частности, используется то, что в большинстве узлов (п, т) величины слева и справа от хт мало отличаются друг от друга. Разработанная в середине пятидесятых годов схема до сих пор применяется в расчетах; при этом она, разумеется, обобщается и совершенствуется.
Расчет контактного разрыва. Проблемы расчета течения, содержащего контактный разрыв, рассмотрим, используя известное нам точное решение уравнений газовой динамики типа «чистый контактный разрыв». В этом случае из всех уравнений газовой динамики нетривиально только одно: pt + Upx = 0 (здесь M = Const). Будем решать его методом сеток.
Построим равномерную сетку с шагом h по х и t по t. Узлы сетки jcm = mh, tn = пх. Приближенное решение ищем в виде сеточной
функции р". Используем простейшую явную схему (предполагая и > 0):
(РГ1 - Pm)/* + “(Р2, - Pm-l)M = 0. (18)
Как известно, эта схема устойчива при условии Куранта их/Л < 1 (см. § 12). (Заметим, что такая схема весьма популярна при расчете задач в эйлеровых координатах: во все уравнения в этом случае
300
ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
[Ч. II
входит характерный оператор «субстанциальной производной» д/д1 + и д/дх.) При и < 0 используется аппроксимация с шаблоном, ориентированным в противоположную сторону (против потока):
(Pm+1 - Pm)/X + “(Pm + l - Рт)/А = °-
Когда решается полная система уравнений газовой динамики,
функция u(t, х) может менять знак. Соответственно, и разностные
формулы строятся в точках (п, т) в зависимости от знака и“. Что же получается
при расчете контактного разрыва по такой схеме? Происходит неприятное явление. В расчете контактный разрыв размывается, его «ширина» растет с ростом времени. На рис. 31 показана характерная картина расчета: точное решение («ступенька») и приближенное.
Это явление особенно неприятно в тех задачах, в которых контактный разрыв разделяет среды с разными уравнениями состояния. Нужно знать достаточно точно границу между разными газами, чтобы пользоваться в данной точке нужным уравнением состояния. На рис. 31 представлены результаты, полученные в вычислительном эксперименте. Ho для того чтобы бороться с «размазыванием контактного разрыва», нужно иметь хоть какую-то теорию. Этим мы сейчас и займемся. Заметим, что излагаемый ниже аппарат имеет достаточно общее значение, не ограничивающееся только задачей о контактном разрыве. Он может применяться при построении разностных схем, лучших в том или ином отношении (смотря по тому, что нам нужно в задаче).
