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

Применение вариационных принципов для построения разностных схем

9.2. Вариационные схемы для решения задач газовой динамики

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

Рассмотрим следующий пример. В бесконечном объеме (вакууме) находится область \Omega, занятая газом (см. рис. 9.3). Масса газа m = \int\limits_\Omega  {\rho dV} , \rho (x, y, z) — плотность среды.

Движение системы подчиняется следующим уравнениям:

\begin{gather*}
 \rho \frac{{d{\mathbf{W}}}}{dt} + {grad}{P} = 0 \quad \mbox{(движение)}; \\ 
 \frac{{d \rho }}{dt} + \rho {div}{\mathbf{W}} = 0 \quad \mbox{(неразрывность)}; \\ 
 \rho \frac{{d \varepsilon }}{dt} + P{div}{\mathbf{W}} = 0 \quad \mbox{(закон сохранения энергии)}; \\ 
F(\rho , P, \varepsilon ) = 0 \quad \mbox{(уравнение состояния)}. \end{gather*}

Эти уравнения приведены в лагранжевых переменных. Напомним, что полная производная определена как

$ \frac{d}{dt} = \frac{{\partial}}{{\partial}t} + ({\mathbf{W}}, \nabla )  $
. Здесь все обозначения традиционные, через \mathbf{W} обозначена скорость движения среды, а через V — элемент объема. Рассматривается случай невязкого нетеплопроводного газа. Запишем лагранжиан для движения среды, представив ее как систему частиц, а затем устремим число частиц в бесконечность.

$  
L = T - U = \int\limits_\Omega  {\frac{{({\mathbf{W}}, {\mathbf{W}})}}{2} \rho 
dV} - \int\limits_\Omega  {\rho \varepsilon dV} = \int\limits_\Omega {\left({\frac{{({\mathbf{W}}, {\mathbf{W}})}}{2} - \varepsilon }\right)dm}.  $

Движение среды должно доставлять минимум функционалу действия

S = \int\limits_0^{t}{L({\tau})d {\tau}}.

По закону сохранения массы \delta (\rho dV) = \delta (dm) = 0. Найдем \delta S - вариацию S:

\delta S = \int\limits_0^{t}{\delta L({\tau})d {\tau}} = \int\limits_0^{t}{d {\tau}} \int\limits_\Omega  {(({\mathbf{W}}, \delta {\mathbf{W}}) -  \delta \varepsilon )dm}.

Согласно первому началу термодинамики в предположении адиабатичности процесса,

$   \delta \varepsilon = \frac{P}{{\rho ^2}} \delta \rho .  $

Используя кинематическое соотношение

$   \delta {\mathbf{W}} = \delta \left({\frac{{d{\mathbf{r}}}}{dt}}\right) =  \frac{d}{dt}(\delta {\mathbf{r}})  $

получим, что

$   \frac{{\delta \rho }}{{\rho ^2}} = \frac{{\delta \left({\frac{1}{V}}
\right)}}{{\rho \left({\frac{1}{V}}\right)}} = - \frac{{\delta (dV)}}{{\rho dV}} = 
 - \frac{1}{\rho }{div}(\delta {\mathbf{r}}).  $

Тогда

\begin{gather*}
 \delta S = - \int\limits_0^{t}{d {\tau}} \int\limits_\Omega  {\left({\rho  \frac{{d{\mathbf{W}}}}{dt} + {grad}{P}}\right) \delta {\mathbf{r}}dV} + \\ 
 + \left. {\left[{\int\limits_\Omega  {({\mathbf{W}}, \delta {\mathbf{r}}) \rho dV}}\right]}
\right|_0^{t} + \int\limits_0^{t}{d\tau} \oint\oint\limits_{\partial\Omega }{(\delta {\mathbf{r}}, {\mathbf{n}})Pd {\sigma}} = 0. \end{gather*}

Независимые вариации траекторий \delta {\mathbf{r}} полагают равными нулю при \tau  = 0 и \tau  = t. Считая также \left. {(\delta {\mathbf{r}}, {\mathbf{n}})}\right|_{{\partial}\Omega } = 0, получим, что вариация лагранжиана равна нулю на решении системы уравнений Эйлера газовой динамики.


Рис. 9.4.

Рассмотрим дискретный аналог функционала действия. Для этого введем в \Omega сетку с ячейками, пронумерованными по левому нижнему углу. При отображении \Omega на квадрат сетка отображается на равномерную квадратную сетку (см. рис. 9.4). Тогда

$  L_{h} = \sum\limits_{i, j}{m_{ij} \left({ \frac{1}{2} < u^2 + v^2  >_{ij} - \varepsilon_{ij}}\right)},   $

где m_{ij} = \rho _{ij}V_{ij} — масса ячейки ( \rho _{ij} — ее плотность) < u2 + v2 > — усредненная по ячейке плотность кинетической энергии

$  < u^2 + v^2  >_{ij} = \frac{1}{4} \sum\limits_{l, k = 
0}^1 {(u_{i + l, j + k}^2 + v_{i + l, j + k}^2 )}.  $

Условие адиабатичности течения имеет следствием для дискретной системы выражение m_{ij} d\varepsilon _{ij} = - P_{ij} dV_{ij}.

Кинематические соотношения (уравнения движения для узлов сетки) есть

$   \frac{dx_{ij}}{dt} = u_{ij},  \frac{dy_{ij}}{dt} = v_{ij}.  $

Объем каждой ячейки вычисляется в предположении, что границы ячейки — отрезки прямых:

$  V_{ij} = \frac{1}{2}[(x_{i + 1j} - x_{ij + 1})(y_{i + 1j + 1} - y_{ij} ) - (x_{i + 1j + 1} - x_{ij} )(y_{i + 1j} - y_{ij + 1} )].  $

Теперь необходимо получить соответствующую систему дифференциально - разностных соотношений, аппроксимирующих уравнения Эйлера.

Запишем функционал действия

S_{h} = \int\limits_0^{t}{L_{h} ({\tau})d {\tau}}

и найдем его вариацию:

\begin{gather*}
0 = \delta S_{h} = \int\limits_0^{t}{\delta \left({\sum\limits_{\omega_{h}}{m_{ij} \left({\frac{{< u^2 + v^2  >_{ij}}}{2} - \varepsilon_{ij}}\right)}}\right)} d {\tau} =  \\ 
 = \sum\limits_{\omega_{h}}{\left\{{\delta \left({\frac{{m_{ij}}}{8} \sum\limits_{l, k = 0}^1 {(u_{i + l, j + k}^2 + v_{i + l, j + k}^2 )}}\right) - \delta (m_{ij} \varepsilon_{ij} )}\right\}} = \\ 
 = \sum\limits_{\omega_{h}}{\frac{1}{8} \delta (m_{ij} \sum\limits_{l, k = 0}^1 {u_{i + l, j + k}^2} )} +  \sum\limits_{\omega_{h}}{\frac{1}{8} \delta (m_{ij} \sum\limits_{l, k = 0}^1 {v_{i + l, j + k}^2} )} +  \sum\limits_{\omega_{h}}{P_{ij} \delta V_{ij}}. 
\end{gather*}

Первое слагаемое в последнем равенстве можно переписать как

$  
 \sum\limits_{\omega_{h}}{\delta (u_{ij}^2 \cdot \frac{1}{8} \sum\limits_{l, k = 0}^1 {m_{i - k, j - l}} )} = M_{ij} u_{ij} \delta u_{ij} = M_{ij} \frac{{du_{ij}}}{dt} \delta x_{ij},  $

где

$  M_{ij} = \frac{1}{4}(m_{i - 1j} + m_{i - 1j - 1} + m_{ij} + m_{ij - 1} ).  $

После преобразований получаются следующие соотношения:

\begin{gather*}
M_{ij} \frac{{du_{ij}}}{dt} =  \sum\limits_{l, k = 0}^1 {P_{i - l, j - k} \frac{{{\partial}V_{i - l, j - k}}}{{{\partial}x_{ij}}}} ; \\ 
M_{ij} \frac{{dv_{ij}}}{dt} =  \sum\limits_{l, k = 0}^1 {P_{i - l, j - k} \frac{{{\partial}V_{i - l, j - k}}}{{\partial y_{ij}}}} ; \end{gather*}

т.е. возникает система дифференциально - разностных уравнений движения.

Уравнение адиабатичности трактуется, как уравнение для определения энергии:

\begin{multline*}
M_{ij} \frac{{d \varepsilon_{ij}}}{dt} =  - P_{ij} \frac{{dV_{ij}}}{dt} = \\ 
 =  - P_{ij} \sum\limits_{l, k = 0}^1 {\left({\frac{{\partial V_{ij}}}{{{\partial}x_{i + l, j + k}}}
u_{i + l, j + k} + \frac{{{\partial}V_{ij}}}{{{\partial}y_{i + l, j + k}}}y_{i + l, j + k}}\right)} ;
 \end{multline*}
  $

Система дифференциально - разностных соотношений замыкается уравнением состояния.

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