1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 63, № 9, 2023

Журнал вычислительной математики и математической физики, 2023, T. 63, № 9, стр. 1554-1563

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

П. В. Чувахов 1, 2, *, И. О. Погорелов 1, 2, **, К. В. Шубин 2, ***

1 Центральный аэрогидродинамический институт им. Н.Е. Жуковского
140180 М.о., Жуковский, ул. Жуковского, 1, Россия

2 Московский физико-технический институт (национальный исследовательский университет)
141701 М.о., Долгопрудный, Институтский пер., 9, Россия

* E-mail: pavel_chuvahov@mail.ru
** E-mail: ilya.pogorelov@phystech.edu
*** E-mail: ottdimile@mail.ru

Поступила в редакцию 27.10.2022
После доработки 26.03.2023
Принята к публикации 28.04.2023

Полный текст (PDF)

Аннотация

Предложен универсальный подход к верификации численного метода решения уравнений Навье–Стокса, который позволяет надежно оценивать его диссипативные свойства. Подход основан на явлении вязкого затухания слабых элементарных возмущений, распространяющихся в однородном потоке. Соответствующее теоретическое решение позволяет определить порядок сходимости численного метода. Библ. 8. Фиг. 5.

Ключевые слова: численная диссипация, численная вязкость, численный метод, уравнения Навье–Стокса, уравнения Эйлера, сходимость.

1. ВВЕДЕНИЕ

Численный метод неизбежно привносит в нестационарное решение дополнительную диссипацию и дисперсию, обусловленную аппроксимацией исходного непрерывного решения на дискретной пространственно-временной сетке. Численная диссипация – или схемная вязкость, или численная вязкость, или диссипативные свойства численного метода и т. д. – зависит от используемого численного метода и обязана уменьшаться по мере уменьшения шага сетки (свойство сходимости численного решения) [1]. При моделировании нестационарных газодинамических явлений численная диссипация проявляется в виде избыточного затухания возмущений, развивающихся в расчетной области. Поскольку это может приводить к искажению исследуемых физических явлений, проводятся исследования диссипативных и дисперсионных свойств конкретных разностных схем (см., например, [2–6]). В качестве примера детального анализа дисперсионных и диссипативных свойств разностных схем для уравнения переноса в линейном случае можно указать монографию [7].

В настоящей работе предлагается универсальный подход к верификации произвольного численного метода решения уравнений газодинамики путем оценки его диссипативных свойств. Рассматривается общий случай полных уравнений Навье–Стокса, а применимость подхода тестируется на примере конкретной реализации метода Годунова. Оценка проводится с помощью методического исследования простейшего процесса – вязкого затухания слабой монохроматической акустической волны или волны энтропии в однородном газовом потоке без газодинамических разрывов. Предполагается, что такая оценка останется справедливой для широкого класса задач.

Для выбранной задачи можно получить точное теоретическое решение.

2. ТЕОРЕТИЧЕСКОЕ РЕШЕНИЕ

Для простоты рассмотрим задачу развития малых возмущений в одномерной постановке, когда вектор направления энергии возмущений направлен вдоль оси $Ox$. После простых алгебраических преобразований уравнения Навье–Стокса, приведенные в [8] в консервативном виде в трехмерной безразмерной постановке, сводятся к виду (добавлено уравнение состояния):

$\frac{{\partial \rho }}{{\partial t}} + u\frac{{\partial \rho }}{{\partial x}} + \rho \frac{{\partial u}}{{\partial x}} = 0,$
$\rho \frac{{\partial u}}{{\partial t}} + \rho u\frac{{\partial u}}{{\partial x}} + \frac{{\partial p}}{{\partial x}} - \frac{4}{{3{\text{R}}{{{\text{e}}}_{{\infty ,L}}}}}\frac{\partial }{{\partial x}}\left( {\mu \frac{{\partial u}}{{\partial x}}} \right) = 0,$
$\frac{1}{\gamma }\frac{{\partial \rho T}}{{\partial t}} + \frac{{\partial \rho uT}}{{\partial x}} - \left( {\gamma - 1} \right){\text{M}}_{\infty }^{2}u\frac{{\partial p}}{{\partial x}} - \frac{1}{{{\text{R}}{{{\text{e}}}_{{\infty ,L}}}{\text{Pr}}}}\frac{\partial }{{\partial x}}\left( {\mu \frac{{\partial T}}{{\partial x}}} \right) - \frac{{4\mu }}{{3{\text{R}}{{{\text{e}}}_{{\infty ,L}}}}}\frac{\partial }{{\partial x}}{{\left( {\frac{{\partial u}}{{\partial x}}} \right)}^{2}}\left( {\gamma - 1} \right){\text{M}}_{\infty }^{2} = 0,$
$\gamma {\text{M}}_{\infty }^{2}p = \rho T.$
В этих уравнениях единичная длина, скорость и температура соответствуют характерному размеру задачи $L{\kern 1pt} *$ и параметрам набегающего потока $V_{\infty }^{*}$, $T_{\infty }^{*}$. Звездочка в верхнем индексе обозначает размерную величину. Единичное давление соответствует удвоенному скоростному напору $\rho _{\infty }^{*}V{{_{\infty }^{*}}^{2}}$. Число Маха набегающего потока вычисляется как ${{{\text{M}}}_{\infty }} = V_{\infty }^{*}{\text{/}}a_{\infty }^{*}$, где $a$ – скорость звука, ${\text{R}}{{{\text{e}}}_{{\infty ,~L}}} = \rho _{\infty }^{*}V_{\infty }^{*}L{\kern 1pt} *{\kern 1pt} {\text{/}}\mu _{\infty }^{*}$ – число Рейнольдса, $\Pr = \mu _{\infty }^{*}c_{p}^{*}{\text{/}}\lambda _{\infty }^{*}$ = 0.72 – число Прандтля, $\gamma = c_{p}^{*}{\text{/}}c_{{v}}^{*} = 1.4$ – показатель адиабаты.

Линейное приближение этой системы относительно однородного набегающего потока – $\left( {u,\rho ,p,T} \right) = \left( {1,~1,~1{\text{/}}\gamma \operatorname{M} _{\infty }^{2},~1} \right) + \left( {u{\kern 1pt} ',\rho {\kern 1pt} ',p{\kern 1pt} ',T{\kern 1pt} '} \right)$ – имеет вид

$\frac{{\partial \rho {\kern 1pt} '}}{{\partial t}} + \frac{{\partial u{\kern 1pt} '}}{{\partial x}} + \frac{{\partial \rho {\kern 1pt} '}}{{\partial x}} = 0,\quad \frac{{\partial u{\kern 1pt} '}}{{\partial t}} + \frac{{\partial u{\kern 1pt} '}}{{\partial x}} + \frac{{\partial p{\kern 1pt} '}}{{\partial x}} - \frac{4}{{3{\text{R}}{{{\text{e}}}_{{\infty ,L}}}}}\frac{{{{\partial }^{2}}u{\kern 1pt} '}}{{\partial {{x}^{2}}}} = 0,$
${\text{M}}_{\infty }^{2}\left( {\frac{{\partial p{\kern 1pt} '}}{{\partial t}} + \frac{{\partial p{\kern 1pt} '}}{{\partial x}}} \right) + \frac{{\partial u{\kern 1pt} '}}{{\partial x}} - \frac{1}{{{{{\operatorname{Re} }}_{{\infty ,L}}}\Pr }}\frac{{{{\partial }^{2}}T{\kern 1pt} '}}{{\partial {{x}^{2}}}} = 0,\quad \gamma {\text{M}}_{\infty }^{2}p{\kern 1pt} ' = \rho {\kern 1pt} ' + T{\kern 1pt} '.$

Будем искать решение в виде бегущей волны $\left( {u{\kern 1pt} ',\rho {\kern 1pt} ',p{\kern 1pt} ',T{\kern 1pt} '} \right) = \left( {\hat {u},\hat {\rho },\hat {p},\hat {T}} \right)\exp \left( {i\alpha x - i\omega t} \right)$, где ${{\omega }}$ – действительная частота, а ${{\alpha }}$ – комплексное волновое число. После такой подстановки последняя система дифференциальных уравнений сводится к системе алгебраических уравнений с малым параметром $\epsilon = 1{\text{/R}}{{{\text{e}}}_{{\infty ,L}}} \ll 1$

$\left( {\begin{array}{*{20}{c}} {i\alpha }&{\gamma {\text{M}}_{\infty }^{2}\left( {i\alpha - i\omega } \right)}&{ - i\alpha + i\omega } \\ {i\alpha - \omega t + \frac{4}{3}{{\alpha }^{2}}\epsilon }&{i\alpha }&0 \\ {i\alpha }&{{\text{M}}_{\infty }^{2}\left( {i\alpha - i\omega } \right)}&{\frac{1}{{\Pr }}{{\alpha }^{2}}\epsilon } \end{array}} \right)\left( {\begin{array}{*{20}{c}} {\hat {u}} \\ {\hat {p}} \\ {\hat {T}} \end{array}} \right) \equiv {\mathbf{A}}\left( {\begin{array}{*{20}{c}} {\hat {u}} \\ {\hat {p}} \\ {\hat {T}} \end{array}} \right) = {\mathbf{0}}.$

В главном приближении по $\epsilon $ (при $\epsilon = 0$) условие разрешимости $\det {\mathbf{A}} = 0$ дает известные дисперсионные соотношения для энтропийной и акустической волн: ${{{{\alpha }}}_{0}} = {{\omega }}$ – для энтропийной волны; ${{{{\alpha }}}_{0}}\left( {1 \pm 1{\text{/}}{{{\text{M}}}_{\infty }}} \right) = {{\omega }}$ – для акустических волн (верхний знак соответствует быстрой волне, нижний – медленной). Первое приближение получим, раскладывая искомое собственное значение ${{\alpha }}$ по формуле Тейлора до $o\left( \epsilon \right)$ при $\epsilon \to 0$

${{\alpha }}\left( {{\omega }} \right) \equiv {{{{\alpha }}}_{0}} + {{{{\alpha }}}_{1}} + o\left( \epsilon \right) = {{{{\alpha }}}_{0}} + {{\left( {\frac{{d{{\alpha }}}}{{d\epsilon }}} \right)}_{{\epsilon = 0}}}\epsilon + o\left( \epsilon \right).$

Производную ${{\left( {\frac{{d{{\alpha }}}}{{d\epsilon }}} \right)}_{{\epsilon = 0}}}$ найдем, дифференцируя систему по параметру $\epsilon $:

${{\left( {\frac{{\partial {\mathbf{A}}}}{{\partial \epsilon }}} \right)}_{{\epsilon = 0}}}\left( {\begin{array}{*{20}{c}} {\hat {u}} \\ {\hat {p}} \\ {\hat {T}} \end{array}} \right) = {\mathbf{0}}.$
После ряда алгебраических манипуляций получим

${{\left( {\frac{{d{{\alpha }}}}{{d\epsilon }}} \right)}_{{\epsilon = 0}}} = i\frac{{{{\alpha }}_{0}^{2}}}{{{\text{Pr}}}}\frac{{{\text{M}}_{\infty }^{2}{{{\left( {{{\alpha }} - {{\omega }}} \right)}}^{2}}\left( {4{\text{Pr}} + 3{{\gamma }}} \right) - 3{{\alpha }}_{0}^{2}}}{{9{\text{M}}_{\infty }^{2}{{{\left( {{{\alpha }} - {{\omega }}} \right)}}^{2}} + {{{{\omega }}}^{2}} - {{{\left( {3{{{{\alpha }}}_{0}} - {{\omega }}} \right)}}^{2}}}}.$

Так как ${{{{\alpha }}}_{0}} \in \mathbb{R}$, элементарные волны не затухают в отсутствие вязкости. В линейном приближении вязкость приводит исключительно к демпфированию колебаний. Декремент вязкого затухания для энтропийной волны (${{{{\alpha }}}_{0}} = {{\omega }}$) есть

${{{{\alpha }}}_{1}} = {{{{\alpha }}}_{{1,i}}} = \frac{{{{{{\omega }}}^{2}}}}{{{\text{R}}{{{\text{e}}}_{{\infty ,L}}}{\text{Pr}}}},$
а для акустической волны (${{{{\alpha }}}_{0}}c = {{{{\alpha }}}_{0}}\left( {1 \pm 1{\text{/}}{{{\text{M}}}_{\infty }}} \right) = {{\omega }}$)

(1)
${{{{\alpha }}}_{1}} = {{{{\alpha }}}_{{1,i}}} = \frac{2}{3}\frac{{{{{{\omega }}}^{2}}}}{{{{c}^{3}}{\text{R}}{{{\text{e}}}_{{\infty ,L}}}}}\left( {1 + \frac{3}{4}\frac{{{{\gamma }} - 1}}{{{\text{Pr}}}}} \right).$

Аналогичный результат для акустической волны получен с помощью общих энергетических соображений в классической литературе [9]. В случае нулевой объемной вязкости (гипотеза Стокса) он согласуется с выражением, полученным выше.

Из полученных выражений следует, что в предельном случае ${\text{R}}{{{\text{e}}}_{{\infty ,L}}} \to \infty $ вязкое затухание отсутствует и любое затухание, наблюдаемое в численном моделировании нестационарных течений в рамках уравнений Эйлера, есть проявление диссипативных свойств численного метода. Исследуем этот случай численно, определив декремент численного затухания как

(2)
${{{{\alpha }}}_{{i,{\text{числ}}}}} = {{{{\alpha }}}_{i}} - {{{{\alpha }}}_{{i,{\text{теор}}}}},$
где ${{{{\alpha }}}_{i}}$ – декремент, восстановленный из численного решения, а ${{{{\alpha }}}_{{i,{\text{теор}}}}}$ – его теоретическое значение, которое в невязком случае равно нулю.

3. ЧИСЛЕННЫЙ МЕТОД

Уравнения Навье–Стокса в двумерной постановке интегрировались численно с помощью авторского пакета расчетных программ [8]. В пакете реализован квазимонотонный численный метод конечных объемов второго порядка аппроксимации по времени и пространству. Интерполяция параметров потока на границы ячеек выполняется методом WENO третьего порядка, решение задачи Римана о распаде разрыва осуществляется методом Roe с малой энтропийной коррекцией. Полученная при аппроксимации система нелинейных уравнений решается модифицированным методом Ньютона, на каждой итерации по нелинейности используется метод GMRes для решения линейной системы уравнений. Расширенное описание численного метода дано в [8].

Расчеты ведутся в квазиодномерной постановке. Для этого используется структурированная сетка в прямоугольной области с ${{N}_{x}} \times 5$ равномерно распределенных узлов, величина ${{N}_{x}}$ варьируется. Направление $y$ моделируется формально на ширину шаблона аппроксимации. На левой (входной) границе $x = 0$ накладываются условия I рода ${\mathbf{q}} \equiv \left( {u,{v},p,T} \right) = {{{\mathbf{q}}}_{b}} + {\mathbf{q}}{\kern 1pt} '$, задающие базовый однородный поток и его нестационарное возмущение ${\mathbf{q}}{\kern 1pt} ' = {{\delta }}{\mathbf{q}} \cdot \sin \omega t$ в виде акустической ${{\delta }}{\mathbf{q}} = {{\varepsilon }} \cdot \left( { \pm {{{\text{M}}}_{\infty }},0,1,\left( {{{\gamma }} - 1} \right){\text{M}}_{\infty }^{2}} \right)$ или энтропийной ${{\delta }}{\mathbf{q}} = \left( {0,0,0,\varepsilon } \right)$ волны с различными амплитудами ${{\varepsilon }}$. На правой (выходной) границе накладывается условие линейной экстраполяции всех газодинамических величин изнутри расчетной области: $\partial {\mathbf{q}}{\text{/}}\partial x = {\mathbf{0}}$. На границах по координате $y$ накладывается условие симметрии потока. В качестве начального условия рассматривается однородный поток ${\mathbf{q}} = {{{\mathbf{q}}}_{b}}$.

4. РЕЗУЛЬТАТЫ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ

4.1. Невязкий cлучай

Расчеты, воспроизводящие невязкий случай ${\text{R}}{{{\text{e}}}_{{\infty ,L}}} \to \infty $, выполнены в рамках полных двумерных уравнений Навье–Стокса при числе Рейнольдса ${\text{R}}{{{\text{e}}}_{{\infty ,L}}} = {{10}^{{16}}}$ в квазиодномерной постановке: плоская монохроматическая волна распространяется в однородном потоке в направлении $x$. В базовом случае рассматривается медленная акустическая волна с частотой ${{\omega }} = 131$ и амплитудой давления ${{\delta }}p = {{10}^{{ - 5}}}{{p}_{\infty }}$, число Маха ${{{\text{M}}}_{\infty }} = 6$. Такие параметры характерны для задач восприимчивости сверхзвуковых пограничных слоев к акустическим возмущениям.

Расчетная сетка построена в прямоугольной области $x \in \left[ {0,2} \right] \times \left[ {0,0.4} \right] \ni y$. Для исследования пространственно-временной сходимости результатов рассматривается поле возмущения давления $p{\kern 1pt} '$ (или температуры $T{\kern 1pt} '$ в случае энтропийной волны) в момент времени ${{t}_{0}} \approx 3$. Проводятся расчеты для всевозможных комбинаций ${{N}_{x}} \times {{N}_{t}}$ с различным количеством сеточных узлов вдоль $x$ – ${{N}_{x}} \in $ [501, 595, 707, 841, 1001, 1190, 1415, 1501, 1683, 2001, 2380, 2501, 2830, 3001, 3366, 4003, 4501, 4760, 5661, 6001, 6732, 7501, 9001, 9520, 11322, 13464, 16012, 18001, 19041, 22644, 36001] – и различным суммарным количеством временных шагов – ${{N}_{t}} \in $ [938, 1115, 1326, 1577, 1875, 2230, 2652, 3153, 3750, 4460, 5303, 6307, 7500, 8919, 10607, 12613, 15000, 17838, 21213, 25227, 30000, 35676, 42426, 50454, 60000, 71352, 84853, 100908, 120000]. Такой объем расчетных случаев представляет собой прямоугольную сетку из $31 \times 29$ расчетов. В координатах ${{\log }_{2}}{{N}_{x}},\;{{\log }_{2}}{{N}_{t}}$ эта сетка близка к равномерной с шагом 1/4. Количество узлов на длину волны составляет ${{N}_{{{\lambda }}}} = {{N}_{x}}{{\omega }}\left( {{{x}_{{{\text{max}}}}} - {{x}_{{{\text{min}}}}}} \right){\text{/}}\left( {2{{\pi }}c} \right) = {{N}_{x}}{{\omega /}}\left( {{{\pi }}c} \right)$, где c – фазовая скорость возмущения ($1 - 1{\text{/}}{{{\text{M}}}_{\infty }}$ для медленной акустической волны или $1$ для энтропийной волны). Количество временных шагов на период возмущения составляет ${{N}_{\tau }} = {{N}_{t}}{{t}_{0}}{{\omega /}}\left( {2{{\pi }}} \right) = 1.5{{N}_{t}}{{\omega /\pi }}$. Финальный момент времени ${{t}_{0}}$ может незначительно отличаться из-за равномерности шага по времени в течение всего расчета.

Декремент ${{{{\alpha }}}_{i}}$ определяется из результатов расчета путем аппроксимации локальных максимумов ${{\left| {p{\kern 1pt} '{\kern 1pt} \left( x \right)} \right|}_{{{\text{max}}}}}$ с помощью экспоненциальной зависимости и метода наименьших квадратов. Точным аналитическим решением данной задачи является незатухающая монохроматическая волна, занимающая всю расчетную область:

$p_{{{\text{теор}}}}^{'}\left( {t,x} \right) = {{\varepsilon }}\sin \left( {{{\omega }}\left( {\frac{x}{c} - t} \right)} \right)~.$

Квинтэссенцией проведенного многопараметрического исследования является фиг. 1, на которой изображено поле логарифмического декремента акустической волны в базовом случае для всей сетки расчетов. Сходимость к точному решению ${{{{\alpha }}}_{{i,{\text{теор}}}}} = 0$ наблюдается при одновременном увеличении пространственного и временного разрешения. При увеличении только одного из них величина численной диссипации выходит на некоторый постоянный уровень и более не изменяется. Этот факт отражает теоретическую точность разностных схем, которая математически выражается в виде остаточного члена $o\left( {{{\Delta }}{{x}^{m}} + {{\Delta }}{{t}^{p}}} \right)$ при анализе порядка аппроксимации схемы с помощью формулы Тейлора вблизи некоторой точки $\left( {{{t}_{0}},{{x}_{0}}} \right)$.

Фиг. 1.

Численная диссипация в зависимости от пространственно-временного разрешения. а – для ${{{{\alpha }}}_{{i,{\text{числ}}}}}$; б – для ${{{{\alpha }}}_{{i,{\text{числ}}}}}{{\lambda }}$.

Таким образом, при ${{N}_{{{{\lambda }},{\text{max}}}}}$ можно оценить порядок сходимости численной диссипации возмущений по времени, а при ${{N}_{{\tau ,{\text{max}}}}}$ – по пространству. Результаты такого анализа представлены на фиг. 2: в обоих случаях достигается практически третий порядок сходимости (порядок по пространству 2.98; порядок по времени 2.97). При расчете последнего несколько точек, отмеченных крестиками, отклоняются от прямолинейной зависимости и исключены из рассмотрения. Отклонение объясняется недостаточным пространственным разрешением ${{N}_{\lambda }}$: при ${{N}_{x}} = {{N}_{{x,{\text{max}}}}}$ с ростом временного разрешения ${{N}_{t}}$ в некоторый момент пространственная составляющая ошибки начинает преобладать над временной составляющей. В результате суммарная ошибка выходит на плато, что проявляется в виде поведения, отмеченного крестиками на фиг. 2, и соответствует фиг. 1.

Фиг. 2.

Сходимость ${{{{\alpha }}}_{{i,{\text{числ}}}}}$ к нулю по пространству и по времени. Крестиками отмечены точки, не участвующие в линейной регрессии.

4.2. Сходимость по времени

Используемый численный метод [8] основан на аппроксимации временных и пространственных производных со вторым порядком точности, хотя для реконструкции газодинамических переменных на гранях ячейки используется схема третьего порядка. Тем не менее полученный третий порядок временной сходимости для декремента затухания ${{{{\alpha }}}_{{i,~{\text{числ}}}}}$ не противоречит следствию теоремы Филиппова–Рябенького, в соответствии с которым порядок аппроксимации должен совпадать с порядком сходимости, если схема устойчива.

Рассмотрим одностороннюю аппроксимацию временной производной на временном шаге $n + 1$, представив числитель формулой Тейлора в окрестности ${{t}^{{\left[ {n + 1} \right]}}}$ с учетом обозначений ${{Q}^{{\left[ n \right]}}} \equiv Q\left( {{{t}^{{\left[ {n + 1} \right]}}} - {{\Delta }}t} \right)$ и ${{Q}^{{\left[ {n - 1} \right]}}} \equiv Q\left( {{{t}^{{\left[ {n + 1} \right]}}} - 2{{\Delta }}t} \right)$, где ${{\Delta }}t$ – шаг численного интегрирования по времени:

$\left( {\frac{{\partial Q}}{{\partial t}}} \right)_{h}^{{\left[ {n + 1} \right]}} = \frac{{3{{Q}^{{\left[ {n - 1} \right]}}} - 4{{Q}^{{\left[ n \right]}}} + {{Q}^{{\left[ {n + 1} \right]}}}}}{{2{{\Delta }}t}} = \frac{{\partial {{Q}^{{\left[ {n + 1} \right]}}}}}{{\partial t}} - \frac{{{{\Delta }}{{t}^{2}}}}{3}\frac{{{{\partial }^{3}}{{Q}^{{\left[ {n + 1} \right]}}}}}{{\partial {{t}^{3}}}} + \frac{{{{\Delta }}{{t}^{3}}}}{4}\frac{{{{\partial }^{4}}{{Q}^{{\left[ {n + 1} \right]}}}}}{{\partial {{t}^{4}}}} + o\left( {{{\Delta }}{{t}^{3}}} \right),$
при ${{\Delta }}t \to 0$.

Очевидно, что в общем случае достигается второй порядок аппроксимации. Однако в рассматриваемом случае слабозатухающей монохроматической волны в однородном потоке все нечетные производные поля течения обращаются в ноль в экстремумах гармонического возмущения – вблизи огибающей возмущения порядок аппроксимации повышается до третьего. Поэтому численный декремент затухания, который получен из аппроксимации этой огибающей, сходится к теоретическому значению также с третьим порядком точности по времени ($p = 3$):

$\mathop {\max }\limits_x \left| {{{\alpha }_{i}} - {{\alpha }_{{i,{\text{теор}}}}}} \right| \equiv {{\left\| {{{\alpha }_{i}} - {{\alpha }_{{i,{\text{теор}}}}}} \right\|}_{\infty }} \leqslant C{{\Delta }}{{t}^{p}}.$

Теперь рассмотрим сходимость по времени в классическом определении

(3)
${{\left\| {q - {{q}_{{{\text{теор}}}}}} \right\|}_{\infty }} \equiv \mathop {\max }\limits_x \left| {q - {{q}_{{{\text{теор}}}}}} \right| \leqslant C{{\Delta }}{{t}^{p}}.$

На фиг. 3 представлены примеры распределения возмущений в зависимости от пространственного и от временного сеточного разрешения. В первом случае фиксировано наиболее подробное временное разрешение ${{N}_{\tau }} = {\text{const}}$. Численное решение развивается практически синфазно с точным решением, а амплитуда стремится к амплитуде точного решения с ростом ${{N}_{{{\lambda }}}}$. Во втором случае фиксировано наиболее подробное пространственное разрешение ${{N}_{\lambda }} = {\text{const}}$. Поведение амплитуды волны не изменяется, однако при недостаточном временном разрешении проявляется численная дисперсия, которая приводит к появлению разности фаз между вычисленным и точным решениями в некоторый момент времени. Разность фаз стремится к нулю с ростом ${{N}_{\tau }}$.

Фиг. 3.

Распределение возмущения давления в медленной акустической волне при разных пространственных и временных разрешениях, $t = 3$; а – для ${{N}_{t}} = {{N}_{{t,{\text{max}}}}}$; б – для ${{N}_{x}} = ~{{N}_{{x,{\text{max}}}}}$.

Разность фаз обусловлена тем, что фазовая скорость акустической волны при численном моделировании отличается от истинной на величину $\frac{{{{{\left( {{{\omega \Delta }}t} \right)}}^{2}}}}{3} \sim N_{\tau }^{{ - 2}}$. Это следует из теоретического анализа линеаризованных уравнений Эйлера с учетом аппроксимации временной производной и хорошо согласуется с проведенными методическими расчетами.

Применяя подход разд. 2 к уравнениям Эйлера и полагая процесс распространения малых возмущений адиабатическим, получаем систему дифференциальных уравнений в матричном виде:

Действие операторов дифференцирования на экспоненциальную функцию сводится к умножению экспоненты на коэффициенты $i{{\alpha }}x$ и $--i{{\omega }}t$ соответственно. Исследование получившейся матрицы на собственные значения приводит к указанному ранее дисперсионному соотношению

$\frac{{{\omega }}}{{{\alpha }}} = 1 \pm \frac{1}{{{{M}_{\infty }}}}.$

Рассмотрим действие использованного оператора численного дифференцирования на экспоненциальную функцию:

${{\left( {\frac{\partial }{{\partial x}}} \right)}_{h}}{{e}^{{ - i{{\omega }}t}}} = \frac{{{{e}^{{ - i{{\omega }}t}}}}}{{2\Delta t}}\left( {3{{e}^{{ - 2i{{\omega }}\Delta t}}} - 4{{e}^{{ - i{{\omega }}\Delta t}}} + 1} \right) = - i{{\omega }}{{e}^{{ - i{{\omega }}t}}}\left( {1 + \frac{{{{{{\omega }}}^{2}}\Delta {{t}^{2}}}}{3} + o\left( {{{{{\omega }}}^{2}}\Delta {{t}^{2}}} \right)} \right),\quad {{{{\omega }}}^{2}}\Delta {{t}^{2}} \to 0.$

Подставляя такой оператор в матричное уравнение и пренебрегая остаточным членом $o\left( {{{{{\omega }}}^{2}}\Delta {{t}^{2}}} \right)$, получим дисперсионное соотношение в виде

$\frac{{{\omega }}}{{{\alpha }}}\left( {1 + \frac{{{{{{\omega }}}^{2}}\Delta {{t}^{2}}}}{3}} \right) = 1 \pm \frac{1}{{{{{\text{M}}}_{\infty }}}}.$

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

Фиг. 4.

Изменение волнового числа в численном решении при применении конечно-разностной схемы дифференцирования по времени.

Из-за численной дисперсии абсолютная погрешность численного решения нарастает вниз по потоку и при разности фаз $\pi $ достигает максимума (фиг. 5а) – истинный максимум гармонического возмущения становится минимумом при некотором $x$. Такая ситуация наблюдается только при достаточно грубом временном разрешении, которое далеко от практики. В остальных случаях наибольшая погрешность достигается вблизи правой границы расчетной области. Используя классическое определение сходимости (3), можно вновь оценить порядок временной сходимости численного решения к точному. Он очень близок к двум (фиг. 5б). Следует отметить, что из линейной регрессии исключены результаты расчетов, в которых разность фаз превышала $\pi $ (крестик на фиг. 5б).

Фиг. 5.

Распределение абсолютной погрешности численного решения $|{\kern 1pt} p{\kern 1pt} {\text{'}} - p_{{{\text{теор}}}}^{'}{\kern 1pt} |$ при ${{N}_{x}} = {{N}_{{x,{\text{max}}}}}$ и анализ временной сходимости.

4.3. Инвариантность численной диссипации

Аналогичные расчеты проведены с помощью того же численного метода (см. разд. 3) при числе ${\text{R}}{{{\text{e}}}_{{\infty ,L}}} = {{10}^{6}}$, для которого монохроматическая волна медленно затухает из-за вязкости и теплопроводности. Оказалось, что полученные выше результаты с высокой точностью повторяются для величины ${{{{\alpha }}}_{{i,~{\text{числ}}}}}$ (2), где теоретическое значение вычисляется по формуле (1). Все иллюстрации, в частности, фиг. 1, визуально не изменяются и поэтому не приводятся повторно. Таким образом, численная диссипация является добавочным эффектом, который отделяется от физической картины явления.

Аналогичное методическое исследование проведено при различных числах Маха ${{{\text{M}}}_{\infty }} = 3,\;4.5,\;7.5,\;9$ как для акустических, так и для энтропийных возмущений с различными частотами ${{\omega }} = 65.5,\;131,\;232$ и амплитудами ${{\varepsilon }} = {{10}^{{ - 4}}},\;{{10}^{{ - 3}}}$. Установлено, что численное затухание возмущений на масштабе длины волны, ${{\alpha }_{{i,{\text{числ}}}}}\lambda $, зависит только от количества точек на длину волны ${{N}_{{{\lambda }}}}$ и количества временных шагов на период возмущения ${{N}_{\tau }}$. Используя эту инвариантность (независимость величины ${{\alpha }_{{i,{\text{числ}}}}}\lambda $ от физических параметров задачи), можно оценивать численную диссипацию, масштабируя фиг. 1 на случай произвольных параметров набегающего потока и элементарного возмущения.

5. ОБЩИЕ ЗАМЕЧАНИЯ

Корректное моделирование физического процесса затухания акустической волны требует чрезмерно подробного пространственно-временного разрешения возмущений. Однако в большинстве практически важных случаев не требуется точно воспроизводить вязкое затухание, допуская погрешность результатов нестационарного моделирования. Например, для 40 точек на длину волны и 200 точек на ее период относительная ошибка в расчете декремента ${{{{\alpha }}}_{i}}$ оказывается большой из-за численной диссипации, но абсолютная величина декремента 0.04 остается достаточно малой, и амплитуда волны уменьшается лишь на 4% на единице расчетной области. Таким образом, численная диссипация акустической волны будет приводить к ограниченной прогнозируемой погрешности результатов численного моделирования.

Численный декремент ${{{{\alpha }}}_{{i,{\text{числ}}}}} = 0.1$ подразумевает, что амплитуда элементарного возмущения уменьшится примерно на 10% на единице длины расчетной области, или на 0.4% на масштабе длины волны базового возмущения. То есть, численная вязкость привносит в физическую задачу отрицательное интегральное усиление возмущений с фактором $N = - 0.1$ ($q = {{q}_{0}}{{e}^{{N\left( x \right)}}}$). Базовый расчетный случай является характерным для моделирования процессов ламинарно-турбулентного перехода в сверхзвуковых пограничных слоях в малошумных условиях. В таких задачах N‑факторы достигают характерных значений более пяти на том же масштабе длины. Таким образом, численная диссипация приводит к ограниченной погрешности расчета эволюции возмущений в пределах нескольких процентов, что не может существенно повлиять на основные выводы, сделанные по результатам таких расчетов.

В полной пространственной постановке численная диссипация зависит от сеточного разрешения во всех трех направлениях. По предыдущему опыту расчетов, она проявляется в каждом направлении индивидуально, т.е. практически не зависит от сеточного разрешения в других направлениях. Поэтому анализ численной диссипации, проведенный выше на одномерной задаче, годится и в многомерном случае.

6. ЗАКЛЮЧЕНИЕ

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

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

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

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

Список литературы

  1. Петров И.Б., Лобанов А.И. Лекции по вычислительной математике: учебное пособие М.: Интернет-Университет Информационных Технологий; БИНОМ. Лаборатория знаний, 2006. 523 с.

  2. Рождественский Б.Л., Яненко Н.Н. Системы квазилинейных уравнений и их приложения к газовой динамике. (Изд. 2-е, перераб. и доп.) М.: Наука, 1978. 688 с.

  3. Аристова Е.Н., Астафуров Г.О. Сравнение диссипативно-дисперсионных свойств компактных разностных схем для численного решения уравнения адвекции // Ж. вычисл. матем. и матем. физ. 2021. Т. 61. № 11. С. 1747–1758. https://doi.org/10.31857/S004446692111002

  4. Шестаковская Е.С., Стариков Я.Е., Клиначева Н.Л. Метод исследования диссипативных свойств разностных схем в эйлеровых координатах// Вестн. ЮУрГУ. Сер. Матем. моделирование и программирование. 2021. Т. 14. № 2. С. 108–116.

  5. Rogov B.V. Dispersive and dissipative properties of the fully discrete bicompact schemes of the fourth order of spatial approximation for hyperbolic equations // Applied Numerical Mathematics. 2019. Vol. 139. P. 136–155.

  6. Головизнин В.М., Карабасов С.А., Козубская Т.К., Максимов Н.В. Схема “Кабаре” для численного решения задач аэроакустики: обобщение на линеаризированные уравнения Эйлера в одномерном случае // Ж. вычисл. матем. и матем. физ. 2009. Т. 49. № 12. С. 2265–2280.

  7. Головизнин В.М., Соловьев А.В. Дисперсионные и диссипативные характеристики разностных схем для уравнений в частных производных гиперболического типа. М.: МАКС Пресс, 2018, 198 с.; http://lim.cs.msu.ru/index.php?id=9

  8. Егоров И.В., Новиков А.В. Прямое численное моделирование ламинарно-турбулентного обтекания плоской пластины при гиперзвуковых скоростях потока // Ж. вычисл. матем. и матем. физ. 2016. Т. 56. № 6. С. 1064–1081 (Comput. Math. Math. Phys., 56: 6 (2016), 1064).

  9. Ландау Л.Д., Лифшиц Е.М. Гидродинамика. Теоретическая физика: Т.VI. (3-е изд., перераб.) М.: Наука, 1986. 736 с.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал