Опубликован: 25.10.2007 | Доступ: свободный | Студентов: 1275 / 290 | Оценка: 4.40 / 4.36 | Длительность: 21:57:00
Специальности: Математик
Лекция 4:

Введение в методы численного решения уравнений газовой динамики

< Лекция 3 || Лекция 4: 123456 || Лекция 5 >

4.4. Разностная схема И.М. Гельфанда для численного решения одномерной системы уравнений газовой динамики

Система уравнений газодинамики решается в области t \ge 0, \eta \in \left[{0, X}\right], отрезок интегрирования разбивается на интервалы узлами \left\{{\eta_i }\right\}_0^{N}. Все интервалы заполнены газом, что соответствует приближению механики сплошной среды. Величины на интервалах считаются кусочно - постоянными. При численном решении определяются шесть функций: u, p, \varepsilon , \{ T, v, x\} — скорость, давление, удельная внутренняя энергия, температура, удельный объем, эйлерова координата.

Для решения задачи используется система одномерных нестационарных уравнений в частных производных, описывающих поведение газа в лагранжевых переменных ( \eta — лагранжева координата):

\begin{gather*}  
\frac{{dv}}{dt} - \frac{du}{{{\partial}\eta }} = 0, \\ 
 \frac{du}{dt} + \frac{{{\partial}(p + Q)}}{{{\partial}\eta }} = 0, \\ 
 \frac{d}{dt} \left({\varepsilon + \frac{{u^2}}{2}}\right) +  \frac{{{\partial}\left[{(p + Q)u}\right]}}{{{\partial}\eta }} =  \frac{{\partial}}{{{\partial}\eta }} \left[{a({T, v}) \frac{{\partial}t}{{{\partial}\eta }}}\right], \\ 
 \frac{dx}{dt} = u.
  \end{gather*} ( 4.5)

Здесь a(T, v) — заданный коэффициент теплопроводности,

$  Q = \frac{\mu }{v} \left({\frac{{\partial}u}{{{\partial}\eta }} - \left|{\frac{{\partial}u}{{{\partial}\eta }}}\right|}\right) \frac{{\partial}u}{{{\partial}\eta }}
  $

- искусственная вязкость Рихтмайера - Неймана, \mu — коэффициент искусственной вязкости. Очевидно, что $ Q \ne 0 $ если

$ \frac{{\partial}u}{{{\partial}\eta }} < 0,  $
при этом
$ \frac{{{\partial} v}}{{\partial}t} < 0  $
или
$ \frac{{{\partial}{\rho}}}{{{\partial} t}} > 0.  $
Так как плотность со временем увеличивается, происходит сжатие газа. В зонах разрежения, где
$ \frac{{{\partial} u}}{{{\partial} x}} > 0  $
и Q = 0, искусственная вязкость действует только в зонах сжатия.

Коэффициент теплопроводности зависит от температуры и задается как

a(T, v) = T^{\gamma} a(T, v), \quad \gamma > 1, \quad a_1 \le a \le a_2 , \quad a_1 > 0.

Начальные данные для рассматриваемой задачи будут

u(0, \eta ) = u_{0},\  x(0, \eta ) = x_{0},\  v(0, \eta ) = v_{0},\  T(0, \eta ) = T_{0} .

Условия на границах выбираем следующие:

$  u(t, 0) = u_1 (t), \quad a \frac{{\partial}T}{{\partial}\eta }(t, 0) = P_1 (t), \quad T(t, X) = T_2 (t),  $

При численном интегрировании полагаем, что значения функций { u, T, , v, x }известны на n слое и задача состоит в вычислении этих же функций на n + 1 слое. Сеточные функции { u_{m^{n + 1}}, T_{m + 1/2}^{n + 1},  v_{m + 1/2}^{n + 1}, x_{m^{n + 1}} } вычисляются с помощью неявной разностной схемы.

В результате разностная схема записывается как

\begin{gather*}  
 \frac{{{{v}}_{m + 1/2}^{n + 1} - {{v}}_{m + 1/2}^{n}}}{\tau} - \frac{{u_{m + 1}^{n + 1} - u_m^{n + 1}}}{2} = 0, m = 1, \ldots , M,  \\ 
 \frac{u_m^{n + 1} - u_m^{n}}{\tau} + \frac{{(p + Q)}_{m + 1/2}^{n + 1} - {(p + Q)}_{m - 1/2}^{n + 1}}{h_m} = 0 \\ 
{h_m = \eta_{m + 1/2} - \eta_{m - 1/2}, m = 1, \ldots , M - 1, } \\ 
 \frac{1}{\tau} \left[{\varepsilon_{m + 1/2}^{n + 1} + \frac{{(u_m^{n + 1})^2 + 
(u_{m + 1}^{n + 1} )^2}}{4} - \varepsilon_{m + 1/2}^{n} -  \frac{{(u_m^{n} )^2 + (u_{m + 1}^{n} )^2}}{4}}\right] + \\ 
 + \frac{1}{{h_{m + 1/2}}} \left[{(p + Q)_{m + 1}^{n + 1}u_{m + 1}^{n + 1} - (p + Q)_m^{n + 1}u_m^{n + 1}}\right] =  \\ 
 = \frac{1}{{h_{m + 1/2}}} \left[{a_{m + 1} \frac{{T_{m + 3/2}^{n + 1} - T_{m + 1/2}^{n + 1}}}{{h_{m + 1}}} - a_m \frac{{T_{m + 1/2}^{n + 1} - T_{m - 1/2}^{n + 1}}}{{h_m }}}\right],  \\ 
h_{m + 1/2} = \eta_{m + 1} - \eta_m , m = 1, \ldots , M - 1,  \\ 
 \frac{x_m^{n + 1} - x_m^{n}}{\tau} - \frac{u_m^{n + 1} - u_m^{n}}{2} = 0, \quad 
m = 0, \ldots , M.   \end{gather*} ( 4.6)

К этим уравнениям системы добавим разностное выражение для вычисления искусственной вязкости

$  Q_{m + 1/2} = \frac{\mu }{{(hv)_{m + 1/2}}} \left[{(u_{m + 1} - u_m ) - \left| {u_{m + 1} - u_m }\right|}\right](u_{m + 1} - u_m )  $

с фиктивным краевым условием QM + 1/2 = 0 и интерполяционное выражение для pm

$  p_m = \frac{{h_{m - 1/2} \cdot p_{m + 1/2} + h_{m + 1/2} \cdot p_{m - 1/2}}}{{h_{m - 1/2} + h_{m + 1/2}}}.  $

Схема имеет второй порядок аппроксимации по координате, а при весовом коэффициенте, равном 0, 5, и второй порядок по времени, причем все точки спектра лежат на единичной окружности. Подробнее о данной схеме в [14.10].

< Лекция 3 || Лекция 4: 123456 || Лекция 5 >