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

Материал из Википедии — свободной энциклопедии
Перейти к: навигация, поиск

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


Пример[править | править вики-текст]

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


\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}

Здесь:

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

Дифференциальная форма[править | править вики-текст]

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


\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\,x d\,t
 = \oint\limits_{\partial \Omega} \rho dx -  \rho u d t

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

\iint\limits_{\Omega} \left( \frac{\partial q}{\partial t} + \frac{\partial f}{\partial x} \right)d\,x d\,t = 0
\Leftrightarrow \oint\limits_{\partial \Omega} (q dx -  f d t)  = 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 - (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^{L}_{j+\frac{1}{2}}|, |D^{L}_{j-\frac{1}{2}}| \right)}

где

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

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

\Delta t = \min_j \Delta t_j


Литература[править | править вики-текст]

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

Ссылки[править | править вики-текст]