Метод Годунова

Метод Годунова — реализация схем сквозного счета, с помощью которых можно рассчитывать газодинамические течения с разрывами параметров внутри расчётной области. Эта схема предложена С. К. Годуновым в 1959 г. Метод Годунова — это вариант метода контрольного объёма. Потоки через боковые грани определяются из решения задачи о распаде произвольного разрыва. Поясним на примере.

Пример

Рассмотрим построение численного метода Годунова первого порядка точности на примере решения системы уравнений одномерной нестационарной газовой динамики, записанной в дивергентной форме:

{\begin{cases}{\begin{array}{lll}{\frac {\partial \rho }{\partial t}}+{\frac {\partial \rho u}{\partial x}}&=&0\\{\frac {\partial (\rho u)}{\partial t}}+{\frac {\partial (p+\rho u^{2})}{\partial x}}&=&0\\{\frac {\partial E}{\partial t}}+{\frac {\partial u(E+p)}{\partial x}}&=&0\end{array}}\end{cases}}

Здесь:

$\rho$ — плотность
$p$ — давление
$u$ — скорость
$E$ — полная энергия в единице объёма

Заметим, что:

Первое уравнение выражает закон сохранения массы
Второе уравнение выражает закон сохранения импульса
Третье уравнение выражает закон сохранения энергии
$E=e+{\frac {\rho u^{2}}{2}}$

$e=\rho \varepsilon ,$

$\varepsilon =\varepsilon (p,\rho ),$

$\varepsilon ={\frac {1}{\gamma -1}}\cdot {\frac {p}{\rho }}$ , где $\gamma$ — показатель адиабаты

Дифференциальная форма

Начальная система может быть записана в более компактной форме:

{\frac {\partial q}{\partial t}}+{\frac {\partial f}{\partial x}}=0

где:

$q$ — вектор консервативных переменных
$q=\left({\begin{array}{c}\rho \\\rho \,u\\E\end{array}}\right)$
$f$ — вектор потоков
$f=\left({\begin{array}{c}\rho \,u\\p+\rho \,u^{2}\\u(E+p)\end{array}}\right)$

Интегральная форма

Вместо дифференциальной формы уравнений выведем новую интегральную форму уравнений, более приспособленную для представления слабого решения. Здесь под слабым решением понимается обобщённая функция, определяемая интегральными равенствами, полученными из соответствующих дифференциальных уравнений и начальных условий задачи. Для этого выделим некоторый контрольный объём $\Omega$ и проинтегрируем систему уравнений по этому объёму. Применим обобщённую теорему Стокса к полученному интегралу от дивергенции (при двух независимых переменных это будет теорема Грина, и формула Остроградского-Гаусса в трёхмерном пространстве). При этом введем направление обхода контура против часовой стрелки.

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

\iint \limits _{\Omega }\left({\frac {\partial \rho }{\partial t}}+{\frac {\partial \rho u}{\partial x}}\right)d\,xd\,t=\oint \limits _{\partial \Omega }\rho dx-\rho udt

Для всей системы уравнений

\iint \limits _{\Omega }\left({\frac {\partial q}{\partial t}}+{\frac {\partial f}{\partial x}}\right)d\,xd\,t=0\Leftrightarrow \oint \limits _{\partial \Omega }(qdx-fdt)=0

Записывая систему в развернутом виде:

{\begin{cases}{\begin{array}{lll}\oint \limits _{\partial \Omega }(\rho \;d\,x-\rho u\;d\,t)&=&0\\\oint \limits _{\partial \Omega }(\rho u\;d\,x-(p+\rho u^{2})\;d\,t)&=&0\\\oint \limits _{\partial \Omega }(E\;d\,x-u(p+E)\;d\,t)&=&0\\\end{array}}\end{cases}}

Аппроксимация

Произведен переход от дифференциальной формы записи исходной системы уравнений к интегральной форме. Интегральная форма записывается в виде равенства нулю интегралов по контуру (границе выделенного контрольного объёма) от векторов консервативных переменных и потоков. Контурный интеграл представляем в виде суммы интегралов по участкам (интервалам) 1-2, 2-3, 3-4, 4-1 контрольного объёма на рисунке (которого пока нет) и на каждом участке аппроксимируем интеграл с использованием метода прямоугольников как произведение подынтегрального выражения в центре интервала на длину интервала интегрирования:

{\begin{array}{c}q_{12}\cdot (x_{2}-x_{1})-f_{12}\cdot (t_{2}-t_{1})+\\\qquad q_{23}\cdot (x_{3}-x_{2})-f_{23}\cdot (t_{3}-t_{2})+\\\qquad \qquad q_{34}\cdot (x_{4}-x_{3})-f_{34}\cdot (t_{4}-t_{3})+\\\qquad \qquad \qquad q_{41}\cdot (x_{1}-x_{4})-f_{41}\cdot (t_{1}-t_{4})=0\end{array}}\Rightarrow q_{34}=q_{12}-{\frac {\Delta t}{\Delta x}}(f_{23}-f_{41})

с учётом равенств, справедливых для контрольного объёма, построенного по декартовой расчётной сетке:

$x_{3}-x_{2}=0$
$x_{1}-x_{4}=0$
$t_{2}-t_{1}=0$
$t_{4}-t_{3}=0$

кроме того:

$x_{2}-x_{1}=\Delta x$
$x_{4}-x_{3}=-\Delta x$
$t_{3}-t_{2}=\Delta t$
$t_{1}-t_{4}=-\Delta t$

находим значения вектора консервативных переменных на интервале 3-4, принадлежащем новому слою:

q_{34}=q_{12}-{\frac {\Delta t}{\Delta x}}(f_{23}-f_{41})\qquad \Rightarrow \qquad q_{j}^{n+1}=q_{j}^{n}-{\frac {\Delta t}{\Delta x}}(f_{j+{\frac {1}{2}}}-f_{j-{\frac {1}{2}}})

В данном случае величинами с полуцелыми индексами обозначены потоки сохраняемых величин через границы расчётной ячейки за время или потоки через боковые грани (2-3 и 4-1) контрольного объёма. Если скорость потока направлена в одну сторону с внешней нормалью к боковой грани, то поток отрицательный, то есть вытекает из контрольного объёма и наоборот.

В развернутом виде:

{\begin{cases}{\begin{array}{lll}\rho _{j}^{n+1}&=&\rho _{j}^{n}-{\frac {\Delta t}{\Delta x}}((\rho u)_{j+{\frac {1}{2}}}-(\rho u)_{j-{\frac {1}{2}}})\\(\rho u)_{j}^{n+1}&=&(\rho u)_{j}^{n}-{\frac {\Delta t}{\Delta x}}((p+\rho u^{2})_{j+{\frac {1}{2}}}-(p+\rho u^{2})_{j-{\frac {1}{2}}})\\E_{j}^{n+1}&=&E_{j}^{n}-{\frac {\Delta t}{\Delta x}}((u(p+E))_{j+{\frac {1}{2}}}-(u(p+E))_{j-{\frac {1}{2}}})\end{array}}\end{cases}}

Потоки через боковые грани, $f_{j+{\frac {1}{2}}}$ и $f_{j-{\frac {1}{2}}}$ определяются из решения задачи о распаде произвольного разрыва.

Постановка граничных условий

Особенностью постановки и реализации граничных условий в методах контрольного объёма (в том числе и в методе Годунова) является необходимость задания или расчета потоков через грань контрольного объёма, совпадающую границей расчётной области. Для первой и последней ячеек расчётного слоя надо определить потоки массы, импульса и энергии через грани.

Часто для задания граничных условий вводятся «виртуальные» расчётные ячейки. Для этого слева от первой ячейки и справа от последней ячейки вводится ещё по одной дополнительной ячейке, в каждой из которых задаются такие параметры течения, чтобы при решении задачи Римана на боковой грани моделировались требуемые потоки.

Типы граничных процедур

Все предположения производятся относительно левой границы

Неподвижная жесткая стенка

Главное условие — отсутствие перетекания потока массы газа через границу, что соответствует условию нулевой скорости потока на данной грани $U=0$ В виртуальной ячейке тогда нужно задать следующие параметры течения:

{\begin{cases}{\begin{array}{lcl}p_{w}&=&p_{1}\\\rho _{w}&=&\rho _{1}\\u_{w}&=&-u_{1}\\\end{array}}\end{cases}}

«w» — параметры в виртуальной ячейке
«1» — параметры в первой ячейке

Получаемые в задаче распада разрыва параметры течения на боковой грани реализуют нулевой поток массы через эту грань.

Резервуар неограниченной ёмкости

Этому случаю математически соответствует задание на грани значение давления ${\tilde {P}}$ . Скорость втекания можно определить по формуле

{\tilde {U}}=u_{1}+{\frac {{\tilde {P}}-p_{1}}{c_{1}}}

При этом:

если ${\tilde {P}}>p_{1}$ , то

$c_{1}={\frac {\sqrt {\rho _{1}\cdot ((\gamma +1){\tilde {P}}+(\gamma +1)p_{1})}}{2}}$
если ${\tilde {P}}<p_{1}$ , то

$c_{1}={\frac {a_{1}\rho _{1}(\gamma -1)\left(1-{\frac {P}{p_{1}}}\right)}{2\gamma \left(1-\left({\frac {P}{p_{1}}}\right)^{\frac {\gamma -1}{2\gamma }}\right)}}$

Втекающий сверхзвуковой поток

Пусть верхнее подчеркивание обозначает параметры сверхзвукового потока, тогда, если ${\bar {U}}>{\bar {c}}={\sqrt {\frac {\gamma {\bar {P}}}{\bar {\rho }}}}$ , то

{\begin{cases}{\begin{array}{lcl}p_{w}&=&{\bar {P}}\\\rho _{w}&=&{\bar {\rho }}\\u_{w}&=&{\bar {U}}\\\end{array}}\end{cases}}

Вытекающий сверхзвуковой поток

При этом в виртуальной ячейке задаются следующие параметры течения:

{\begin{cases}{\begin{array}{lcl}p_{w}&=&p_{1}\\\rho _{w}&=&\rho _{1}\\u_{w}&=&u_{1}\\\end{array}}\end{cases}}

Выбор параметров сетки

Шаг расчётной сетки по временной координате в методе Годунова можно определить из критерия устойчивости Куранта — Фридрихса — Леви. Применительно к рассматриваемой схеме это условие формулируется следующим образом:

Волны, возникающие в задаче распада произвольного разрыва в точке $j+{\frac {1}{2}}$ , не должны за время $\Delta t$ достигать боковых граней $j+{\frac {3}{2}}$ и $j-{\frac {1}{2}}$ и искажать автомодельное решение.

Реализация этого принципа приводит к следующим соотношениям:

\Delta t_{j}=r{\frac {\Delta x_{j}}{\max \left(|D_{j+{\frac {1}{2}}}^{L}|,|D_{j-{\frac {1}{2}}}^{L}|\right)}}

где

$D_{j+{\frac {1}{2}}}^{L}$ — значение скорости самой левой волны в распаде разрыва;
$D_{j-{\frac {1}{2}}}^{R}$ — значение скорости самой правой волны в распаде разрыва;

В итоге мы берем:

\Delta t=\min _{j}\Delta t_{j}

Литература

Численное решение многомерных задач газовой динамики. Альбом / редактор Годунов С. К. . — М.: Наука, 1976. — 400 с. — 6500 экз.
Самарский А.А., Попов Ю.П. Разностные методы решения задач газовой динамики. — М.: Наука, 1992. — 2470 экз.