Solution for the fluid dynamics problems by the discontinuous Galerkin method in the curvolinear coordinate system

Abstract

The development of a method for solving the one-dimensional fluid dynamics equation by the discontinuous Galerkin method is considered. Test calculations have been done for the problem of a strong point explosion.

Full Text

Решение задач газовой динамики разрывным методом Галёркина в криволинейной системе координат¹

Введение. Численные методы повышенного порядка точности для решения задач механики сплошной среды на сегодняшних день развиваются многими научными коллективами [1–4]. Наиболее перспективным методом является метод Галёркина с разрывными базисными функциями [5]. Он сочетает в себе качества конечно-объемных и конечно-элементных методов, что позволяет получать с достаточной точностью решать задачи, решения которых характеризуются высокими значениями градиентов. В данной работе построена методика решения одномерных задач газовой динамики, обладающих сферической симметрией. Выполнены тестовые расчеты для задачи о сильном точечном взрыве [6; 7].

Разрывный метод Галёркина. Определим разрывный метод Галёркина [5] для одномерной неоднородной системы уравнений Эйлера.

$\frac{\partial q}{\partial t}+\frac{1}{r^n}\frac{\partial \left(r^{n}f\right)}{\partial r}+h=\mathrm{0,}\text{ (1)}$

где

$\text{q=}\left(\begin{array}{c}\rho \\ \rho u\\ \rho E\end{array}\right);\text{ (2)}$

$\text{f=}\left(\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uH\end{array}\right);\text{ (3)}$

$\text{h=}\frac{n}{r}\left(\begin{array}{c}0\\ -p\\ 0\end{array}\right)\text{. (4)}$

В приведенных выше уравнениях ρ – плотность жидкости, u – скорость жидкости, p – статическое давление, E и H – полная энергия и полная энтальпия на единицу массы.

$\text{E=}\frac{p}{\rho \left(\gamma -1\right)}+\frac{u^2}{2},\text{ (5)}$

$\text{H=E+}\frac{p}{\rho }\text{. (6)}$

Для нахождения давления используем уравнение состояния идеального газа

$\text{p=}\left(\gamma -1\right)\rho \epsilon ,\text{ (7)}$

где $\gamma$ – показатель адиабаты.

В уравнениях (1) – (4) n = 0 для декартовых координат, n = 1 для цилиндрических координат и n = 2 для сферических координат.

Для применения разрывного метода Галеркина введем равномерную сетку

$0<r_{\frac{1}{2}}<\mathrm{...}<r_{i-\frac{1}{2}}<r_{i+\frac{1}{2}}<\mathrm{...}<r_{N-\frac{1}{2}}<1.\text{ (8)}$

На каждом отрезке определим базис

${\phi }_{0k}\left(t\right)=\mathrm{1,}\text{ (9)}$

${\phi }_{1k}\left(t\right)=\frac{r-r_c}{\Delta r},\text{ (10)}$

${\phi }_{2k}\left(t\right)={\left(\frac{r-r_c}{\Delta r}\right)}^2,\text{ (11)}$

$\Delta r=r_{i+\frac{1}{2}}-r_{i-\frac{1}{2}}.\text{ (12)}$

Приближенное решение 𝑞ℎ в каждой ячейке $\left[r_{i-\frac{1}{2}},r_{i+\frac{1}{2}}\right]$ представляется в виде разложения по базисным функциям (9) – (12).

$q_h\left(r,t\right)=\sum _{i=0}^2{q}_{ik}\left(t\right){\phi }_{ik}\left(r\right),\text{ (13)}$

Разрывный метод Галеркина получается путем умножения уравнения Эйлера (1) на базисную функцию ${\phi }_k$ и интегрирование по некоторому объему (K).

$\underset{K}{\int }\frac{\partial q}{\partial t}{\phi }_k\left(r\right)dK+\underset{K}{\int }\frac{1}{r^n}\frac{\partial \left(r^{n}f\right)}{\partial r}{\phi }_k\left(r\right)dK+\underset{K}{\int }h{\phi }_k\left(r\right)dK=\mathrm{0,}\text{ (14)}$

где 𝑑𝐾= 𝑟^𝑛𝑑𝑟 – для одномерного потока.

Для одного пространственного измерения уравнение (14) принимает вид

$\begin{array}{l}\underset{r_{i-\frac{1}{2}}}{\overset{r_{i+\frac{1}{2}}}{\int }}\frac{\partial q}{\partial t}{\phi }_k\left(r\right)r^{n}dr+{\left[f{r}^n{\phi }_k\left(r\right)\right]}_{\frac{i+1}{2}}-{\left[f{r}^n{\phi }_k\left(r\right)\right]}_{\frac{i-1}{2}}-\\ -\underset{r_{i-\frac{1}{2}}}{\overset{r_{i+\frac{1}{2}}}{\int }}f\frac{d{\phi }_k\left(r\right)}{dr}{r}^{n}dr+\underset{r_{i-\frac{1}{2}}}{\overset{r_{i+\frac{1}{2}}}{\int }}h{\phi }_k\left(r\right)r^{n}dr=0.\text{ (15)}\end{array}$

Функция f не определена в точках $r_{i-\frac{1}{2}}$ и $r_{i+\frac{1}{2}}$ – решение может быть разрывным на гранях элементов, что приведет к неоднозначности. Эта проблема преодолевается с помощью соответствующего выбора численных потоков на гранях, зависящих от предельных значений функции q_ℎ слева и справа от точек $r_{i-\frac{1}{2}}$ и $r_{i+\frac{1}{2}}$.

Будем использовать дискретные потоки Русанова-Лакса-Фридрихса.

${\widehat{f}}_{i+\frac{1}{2}}=\widehat{f}\left(q_{i+\frac{1}{2}}^L,q_{i+\frac{1}{2}}^R\right)=\frac{1}{2}\left(f\left(q_{i+\frac{1}{2}}^R\right)+f\left(q_{i+\frac{1}{2}}^L\right)-\alpha \left(q_{i+\frac{1}{2}}^R-q_{i+\frac{1}{2}}^L\right)\right),\text{ (16)}$

$\alpha =\max \left\{\left|q_{i+\frac{1}{2}}^L\right|+c^L,\left|q_{i+\frac{1}{2}}^R\right|+c^R\right\},\text{ c=}\sqrt{\frac{\gamma p}{\rho }}.\text{ (17)}$

Используя приближение для функции 𝑞_ℎ, получим систему дифференциальных уравнений

$\frac{dq}{dt}=-M^{-1}\left[\left(\widehat{f}{r}^n\left(\begin{array}{c}{\phi }_0\\ {\phi }_1\\ {\phi }_2\end{array}\right)\right)\left|{}_{r_{i+\frac{1}{2}}}^{r_{i-\frac{1}{2}}}\right.-\underset{r_{i-\frac{1}{2}}}{\overset{r_{i+\frac{1}{2}}}{\int }}f{r}^n\frac{\partial }{\partial r}\left(\begin{array}{c}{\phi }_0\\ {\phi }_1\\ {\phi }_2\end{array}\right)dr+\underset{r_{i-\frac{1}{2}}}{\overset{r_{i+\frac{1}{2}}}{\int }}h{r}^n\left(\begin{array}{c}{\phi }_0\\ {\phi }_1\\ {\phi }_2\end{array}\right)dr\right].\text{ (18)}$

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

Согласно [5] будем обозначать действие оператора лимитирования на функцию q следующим образом: ΛП_ℎq.Одним из простых подходов ограничения служит лимитер Кокбурна.

Запишем разложение решения по линейному базису в ячейке:

$q=q_0+q_1\frac{r-r_c}{\Delta r},\Delta r\text{= }{r}_{i+\frac{1}{2}}-r_{i+\frac{1}{2}}.\text{ (19)}$

Для функции (19) действие оператора ΛПℎ𝑞 запишем как

$Л{П}_{h}q(r,t)=q_0+{\tilde{q}}_1\frac{r-r_c}{\Delta r},\Delta r\text{= }{r}_{i+\frac{1}{2}}-r_{i-\frac{1}{2}}.\text{ (20) }$

Значение функции ${\tilde{q}}_1$ вычисляется как

${\tilde{q}}_1=2\ast \text{minmod}\left[q\left(r_{i+\frac{1}{2}}\right)-q_{0i},\alpha \left({\overline{q}}_{i+\frac{1}{2}}-q_{0i}\right),\alpha \left(q_{0i}-{\overline{q}}_{i-\frac{1}{2}}\right)\right],\text{ (21)}$

где $q_{0i}$ - среднее значение на интервале $r_{i-\frac{1}{2}}$, $r_{i+\frac{1}{2}}$, а ${\overline{q}}_{i+\frac{1}{2}}$ и ${\overline{q}}_{i-\frac{1}{2}}$ равны

${\overline{q}}_{i+\frac{1}{2}}=\frac{q_{0i+1}+q_{0i}}{2},{\overline{q}}_{i-\frac{1}{2}}=\frac{q_{0i-1}+q_{0i}}{2}\text{. (22) }$

Если порядок полинома больше 1, то

$q^l(r,t)=q_0^l+q_1^l\frac{r-r_c}{\Delta r},\Delta r=r_{i+\frac{1}{2}}-r_{i-\frac{1}{2}},\text{ (23)}$

где $q_0^l=q_0+\frac{q_2}{12}$ среднее интегральное значение на интервале $r_{i+\frac{1}{2}}$, $r_{i-\frac{1}{2}}$ для полинома степени равной 2.

После применения лимитера

$Л{П}_h{q}^l(r,t)={\tilde{q}}_o^l+{\tilde{q}}_1^l\frac{r-r_c}{\Delta r},\Delta r\text{= }{r}_{i+\frac{1}{2}}-r_{i-\frac{1}{2}}.\text{ (24) }$

Постановка задачи Седова о сильном точечном взрыве [6]. Поместим небольшое количество безразмерной энергии $\epsilon$ = 1 в небольшую область радиуса dr в центре сетки, заполненной неподвижной средой плотности ${\rho }_0$, с давлением $p_0$. Безразмерное давление внутри этого объема определяется:

$p_0=\frac{3\left(\gamma -1\right)\epsilon }{\left(\nu +1\right)\pi d{r}^{\nu }},x<xc,\text{ (25)}$

$p_0={10}^{-5},x>xc,\text{ (26) }$

где v = 2 для цилиндрических координат, v = 3 для сферических координат.

Выбираем dr в 3,5 раза больше, чем шаг сетки. Плотность устанавливается равной ${\rho }_0$=1 по всей сетке, скорость изначально в покое u₀ = 0.

На рисунке 1 представлены результаты расчетов, масштабированные относительно значений на фронте ударной волны. Результаты расчетов хорошо согласуются с данными из классических источников [6; 7].

Рис. 1. Распределения плотности, скорости и давления, нормированные относительно значений на фронте ударной волны.

¹ Работа выполнена при финансовой поддержке РФФИ и Правительства Республики Мордовия (проект № 18-41-130001 р_а)