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

Журнал вычислительной математики и математической физики, 2022, T. 62, № 8, стр. 1386-1401

Анализ законов сгущения сеток в пограничном слое на примере численного решения задачи обтекания пластины вязким газом

А. Н. Кудрявцев 1, 2, *, В. Д. Лисейкин 2, 3, **, А. В. Мухортов 2, ***

1 Институт теоретической и прикладной механики СО РАН
630090 Новосибирск, ул. Институтская, 4/1, Россия

2 Новосибирский государственный университет
630090 Новосибирск, ул. Пирогова, 1, Россия

3 Федеральный исследовательский центр информационных и вычислительных технологий
630090 Новосибирск, пр-т Академика Лаврентьева, 6, Россия

* E-mail: alex@itam.nsc.ru
** E-mail: liseikin.v@gmail.com
*** E-mail: a.mukhortov@g.nsu.ru

Поступила в редакцию 10.10.2021
После доработки 21.01.2022
Принята к публикации 11.03.2022

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

Аннотация

Численно исследуется решение задачи о течении вязкого сжимаемого газа над плоской пластиной, помещенной в сверхзвуковой поток под нулевым углом атаки. Решаются двумерные уравнения Навье–Стокса для различных чисел Рейнольдса с применением адаптивных сеток, сгущающихся в зоне пограничного слоя. Рассматриваются разностные сетки, построенные с помощью координатных преобразований, устраняющих пограничные слои различных типов. В серии вычислительных экспериментов проведен анализ характеристик численных решений (значение и порядок погрешности, значение и порядок скачка решения, время расчета) и сделаны выводы о преимуществах и недостатках, а также допустимости использования каждого закона сгущения сетки в пограничном слое для нахождения численного решения данной задачи. Новизна работы состоит в анализе специальных адаптивных сеток и их использовании для решения задач, имеющих применение в различных областях сверхзвуковой аэродинамики и газодинамики. Библ. 16. Фиг. 14. Табл. 7.

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

ВВЕДЕНИЕ

Конструирование специальных разностных сеток является одним из подходов в создании алгоритмов решения задач с пограничными и внутренними слоями. Такими являются задачи, характеризуемые наличием малого параметра $\varepsilon $ в коэффициентах при старших производных. Наиболее популярной и хорошо изученной является двухточечная задача для обыкновенного дифференциального уравнения:

(1)
$\begin{gathered} - {{[\varepsilon + d(y)]}^{\alpha }}u{\text{''}} + a(y,u)u{\text{'}} + f(y,u) = 0,\quad \alpha > 0,\quad d(y)\; \geqslant \;0,\quad 0 < y < 1, \\ u(0) = {{A}_{0}},\quad u(1) = {{A}_{1}}, \\ \end{gathered} $
моделирующая качественное поведение решений различных многомерных задач относительно координаты $y$, трансверсальной к пограничному слою (см. [1]). Задача (1) для $d(y) = 0$ аналитически исследована в [2], [3] и [4], где было показано, что ее решение может иметь экспоненциальные пограничные слои и степенные II рода внутренние слои, а в [1] продемонстрирoвано, что для более произвольных функций $d(y)$, $a(y,u)$ и $f(y,u)$ решения могут иметь другие слои, в часности, степенные I рода, логарифмические и смешанные пограничные и внутренние слои. Для экспоненциального пограничного слоя возле $y = 0$ справедливы следующие оценки производных решения:
(2)
$\left| {{{u}^{{(p)}}}(y,\varepsilon )} \right|\;\leqslant \;M\left[ {{{\varepsilon }^{{ - pk}}}\exp \left( { - by{\text{/}}{{\varepsilon }^{k}}} \right) + 1} \right],\quad b > 0,\quad 1\;\leqslant \;p\;\leqslant \;n,\quad 0\;\leqslant \;y\;\leqslant \;m,$
для степенного первого типа –
(3)
$\left| {{{u}^{{(p)}}}(y,\varepsilon )} \right|\;\leqslant \;M\left[ {{{\varepsilon }^{{bk}}}{\text{/}}{{{({{\varepsilon }^{k}} + y)}}^{{b + p}}} + 1} \right],\quad b > 0,\quad 1\;\leqslant \;p\;\leqslant \;n,\quad 0\;\leqslant \;y\;\leqslant \;m,$
для степенного второго типа –
(4)
$\left| {{{u}^{{(p)}}}(y,\varepsilon )} \right|\;\leqslant \;M\left[ {{{{\left( {{{\varepsilon }^{k}} + y} \right)}}^{{b - p}}} + 1} \right],\quad 1 > b > 0,\quad 1\;\leqslant \;p\;\leqslant \;n,\quad 0\;\leqslant \;y\;\leqslant \;m,$
для логарифмического пограничного слоя –
(5)
$\left| {{{u}^{{(p)}}}(y,\varepsilon )} \right|\;\leqslant \;M\left[ {{{1 + {{{\left( {{{\varepsilon }^{k}} + y} \right)}}^{{ - p}}}} \mathord{\left/ {\vphantom {{1 + {{{\left( {{{\varepsilon }^{k}} + y} \right)}}^{{ - p}}}} {\ln (1{\text{/}}\varepsilon )}}} \right. \kern-0em} {\ln (1{\text{/}}\varepsilon )}}} \right],\quad 1\;\leqslant \;p\;\leqslant \;n,\quad 0\;\leqslant \;y\;\leqslant \;m,$
где $M$, $b$ и $m$ – некоторые положительные константы, не зависящие от $\varepsilon $, а $k > 0$ – масштаб пограничного слоя, т.е. первая производная в точке $y = 0$ оценивается величиной $m{\text{/}}{{\varepsilon }^{k}}$. В частности, для задач газовой динамики с граничным условием прилипания масштаб пограничного слоя равен $1{\text{/}}2$ в предположении $\varepsilon = M{\text{/Re}}$.

Конечно, для реальной задачи пограничный слой может быть и нового типа, который еще предстоит открыть, однако, выражения (2)–(5) показывают многообразие структур пограничных слоев.

Очевидным, но не единственным, представителем функции пограничного слоя является для: экспоненциального типа –

(6)
$u(y,\varepsilon ) = \exp \left( { - y{\text{/}}{{\varepsilon }^{k}}} \right),\quad y > 0,$
степенного первого рода –
(7)
$u(y,\varepsilon ) = {{{{\varepsilon }^{k}}} \mathord{\left/ {\vphantom {{{{\varepsilon }^{k}}} {\left( {{{\varepsilon }^{k}} + y} \right)}}} \right. \kern-0em} {\left( {{{\varepsilon }^{k}} + y} \right)}},\quad y > 0,$
степенного второго рода –
(8)
$u(y,\varepsilon ) = {{\left( {{{\varepsilon }^{k}} + y} \right)}^{b}},\quad 0 < b < 1,\quad y > 0,$
логарифмического типа –

(9)
$u(y,\varepsilon ) = {{\ln \left( {{{\varepsilon }^{k}} + y} \right)} \mathord{\left/ {\vphantom {{\ln \left( {{{\varepsilon }^{k}} + y} \right)} {\ln (\varepsilon )}}} \right. \kern-0em} {\ln (\varepsilon )}},\quad y > 0.$

Слоем для $n$-й производной функции является узкий интервал, в котором значения данной производной не ограничены при $\varepsilon $, стремящемся к нулю, а вне слоя они равномерно ограничены (см. [1, с. 17]). Длина этого интервала будет шириной слоя. Так, для функции (6) с экспоненциальным пограничным слоем его ширина равна $(kn{\text{/}}b)\varepsilon \ln (1{\text{/}}\varepsilon )$; для функции (7) ширина равна $m{{\varepsilon }^{{k/(n + 1)}}}$; для функции (8) ширина равна $m$; для функции (9) ширина равна $m{\text{/}}{{(\ln (1{\text{/}}\varepsilon ))}^{{1/n}}}$, где $m$ – константа, не зависящая от $\varepsilon $.

Высокоэффективными для решения задачи (1) с экспоненциальными слоями являются сетки Бахвалова (см. [2]) и Шишкина (см. [5]). Однако для других слоев вида (3)–(5) их сетки не столь эффективны, так как эти слои несравнимо шире экспоненциальных (см. [6]). Сетки для решения задачи (1) с такими слоями описаны в [7].

Перспективными для численного решения задачи (1) являются алгоритмы, использующие координатные преобразования $y(\xi ,\varepsilon ):[0,\;1] \to [0,\;1]$, устраняющие сингулярности высокого порядка решений $u(y,\varepsilon )$, т.е. производные относительно зависимой переменной ${{u}_{1}}(\xi ,\varepsilon ) = u[y(\xi ,\varepsilon ),\varepsilon ]$ по $\xi $ должны быть равномерно ограничены по $\varepsilon $:

(10)
$\left| {\frac{{{{d}^{j}}}}{{d{{\xi }^{j}}}}u\left[ {y(\xi ,\varepsilon ),\varepsilon } \right]} \right|\;\leqslant \;M,\quad j\;\leqslant \;n,\quad 0\;\leqslant \;\xi \;\leqslant \;1,$
где константа $M$ не зависит от параметра $\varepsilon ,$ число $n$ зависит от порядка аппроксимации задачи (чем больше порядок, тем больше $n$). Координатные преобразования $y(\xi ,\varepsilon ):[0,\;1] \to [0,\;1]$, устраняющие сингулярности, порождают сгущающиеся в слоях сетки ${{y}_{i}},$ $i = 0,\;1,\; \ldots ,\;N$, ${{y}_{0}} = 0$, ${{y}_{N}} = 1$ с помощью формулы ${{y}_{i}} = y(i{\text{/}}N,\varepsilon )$. На таких сетках разности значений численного решения в соседних узлах (как в слое, так и вне слоя) при увеличении числа узлов уменьшаются равномерно по $\varepsilon $ в такой же пропорции, т.е. при утроении числа узлов скачок решения уменьшается втрое. Принцип построения таких преобразований для $n\;\leqslant \;3$ был сформулирован в [1] для решения задачи (1) со специальными функциями $d(y)$, $a(y,u)$ и $f(y,u)$ на простой схеме первого порядка с направленными разностями. Была доказана равномерная по малому параметру сходимость первого порядка численного решения к точному. Однако для гарантии большего порядка сходимости для схем высокого порядка эти преобразования не подходят. Причина в том, что сетки, полученные с помощью этих координатных преобразований, не обеспечивают сгущение узлов в слоях для высших производных, входящих в погрешность аппроксимации задачи схемами высокого порядка, так как эти слои шире слоев для производных меньшего порядка. Далее, схемы высокого порядка, в отличие от схемы с направленными разностями, не обладают, как правило, свойством обратной монотонности, которая позволяет установить связь между ошибкой решения и погрешностью аппроксимации и, таким образом, доказать равномерную по малому параметру сходимость. Для таких схем равномерную сходимость можно обосновать только численными расчетами (см. [8]). В этих случаях координатные преобразования для построения сеток, обеспечивающих равномерную сходимость высокого порядка, должны удовлетворять (10) для больших $n$. Такие преобразования могут быть использованы для формулировки метрик в методах построения адаптивных многомерных сеток как структурных (см. [9]), так и неструктурных (см. [10]).

Оценки (10) могут гарантировать равномерную сходимость для преобразованной задачи относительно зависимой переменной ${{u}_{1}}(\xi ,\varepsilon ) = u\left[ {y(\xi ,\varepsilon ),\varepsilon } \right]$ на равномерной сетке, так как погрешность аппроксимации будет равномерно ограничена для $n\; \geqslant \;p + 2$, где $p$ – порядок схемы. Возможно, равномерная сходимость высокого порядка будет и для задачи в физической переменной $y$. Более того, численное решение может быть равномерно проинтерполировано на весь интервал $[0,\;1]$.

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

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

В настоящей статье на примере данной задачи рассматривается стратегия построения эффективных расчетных сеток в подобных случаях. Она основана на использовании для построения сеток известных координатных преобразований, устраняющих особенности (2)–(5). Можно надеяться, что такой подход позволит строить более эффективные сетки и повысить точность решения и при численном моделировании обтекания тел сложной формы, часто встречающихся в различных практических приложениях.

1. ФОРМУЛИРОВКА ЗАДАЧИ И ЧИСЛЕННЫЙ МЕТОД

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

(11)
$\frac{{\partial Q}}{{\partial t}} + \frac{{\partial F}}{{\partial x}} + \frac{{\partial G}}{{\partial y}} = \frac{M}{{{\text{Re}}}}\left( {\frac{{\partial {{F}^{v}}}}{{\partial x}} + \frac{{\partial {{G}^{v}}}}{{\partial y}}} \right).$
Здесь векторы консервативных переменных $Q$, невязких потоков $F$, $G$ и вязких потоков ${{F}^{v}}$, ${{G}^{v}}$ равны соответственно
$Q = \left( {\begin{array}{*{20}{c}} \rho \\ {\rho u} \\ {\rho v} \\ e \end{array}} \right),\quad F = \left( {\begin{array}{*{20}{c}} {\rho u} \\ {\rho {{u}^{2}} + p} \\ {\rho uv} \\ {(e + p)u} \end{array}} \right),\quad G = \left( {\begin{array}{*{20}{c}} {\rho v} \\ {\rho uv} \\ {\rho {{v}^{2}} + p} \\ {(e + p)v} \end{array}} \right),$
(12)
${{F}^{v}} = \left( {\begin{array}{*{20}{c}} 0 \\ {{{\tau }_{{xx}}}} \\ {{{\tau }_{{xy}}}} \\ {u{{\tau }_{{xx}}} + v{{\tau }_{{xy}}} + \varkappa \frac{{\partial T}}{{\partial x}}} \end{array}} \right),\quad {{G}^{v}} = \left( {\begin{array}{*{20}{c}} 0 \\ {{{\tau }_{{xy}}}} \\ {{{\tau }_{{yy}}}} \\ {u{{\tau }_{{xy}}} + v{{\tau }_{{yy}}} + \varkappa \frac{{\partial T}}{{\partial y}}} \end{array}} \right),$
${{\tau }_{{xx}}} = \mu \left( {\frac{4}{3}\frac{{\partial u}}{{\partial x}} - \frac{2}{3}\frac{{\partial v}}{{\partial y}}} \right),\quad {{\tau }_{{xy}}} = \mu \left( {\frac{{\partial u}}{{\partial y}} + \frac{{\partial v}}{{\partial x}}} \right),\quad {{\tau }_{{yy}}} = \mu \left( {\frac{4}{3}\frac{{\partial v}}{{\partial y}} - \frac{2}{3}\frac{{\partial u}}{{\partial x}}} \right).$
Система уравнений замыкается уравнением состояния совершенного газа

(13)
$p = \frac{{\rho T}}{\gamma }.$

В (11)–(13) $u$, $v$ – компоненты вектора скорости вдоль осей $x$ и $y$ соответственно; $\rho $ – плотность; $p$ – давление; $e = p{\text{/}}(\gamma - 1) + \rho {{\left( {{{u}^{2}} + {{v}^{2}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{u}^{2}} + {{v}^{2}}} \right)} 2}} \right. \kern-0em} 2}$ – полная энергия на единицу объема; $\tau $ – тензор вязких напряжений; $\varkappa = \mu {\text{/}}[(\gamma - 1){\text{Pr}}]$ – коэффициент теплопроводности; $T$ – температура; $\gamma $ – отношение  удельных  теплоемкостей;  $\mu $ –  коэффициент  динамической вязкости (зависимость вязкости от температуры аппроксимировалась по формуле Сазерленда $\mu = {{T}^{{1.5}}}(1 + C{\text{/}}T_{\infty }^{ * }){\text{/}}(T + C{\text{/}}T_{\infty }^{ * })$, $C$ – постоянная Сазерленда); ${{\operatorname{Re} }_{L}} = \rho _{\infty }^{ * }U_{\infty }^{ * }{{L}^{ * }}{\text{/}}\mu _{\infty }^{ * }$ – число Рейнольдса, вычисляемое по параметрам набегающего потока (индекс “$\infty $”) и длине пластины $L{\text{*}}$; ${\text{M}}$, ${{U}_{\infty }}$ – число Маха и скорость набегающего потока; ${\text{Pr}}$ – число Прандтля; индекс “$ * $” соответствует размерным величинам.

Расчетная область для численного решения представляет собой прямоугольник размером $1 \times 0.32$. Для системы уравнений (11)–(13) ставятся следующие граничные и начальные условия. На поверхности пластины требуется выполнение условия прилипания: ${{\left. u \right|}_{\Gamma }} = 0$. Кроме того, требуется задать граничное условие для температуры: на поверхности пластины было наложено адиабатическое температурное условие, т.е. тепловой поток через поверхность пластины предполагается равным нулю, что означает равенство нулю производной по нормали к поверхности: ${{\left. {(\partial T{\text{/}}\partial n)} \right|}_{\Gamma }} = 0$. На входной границе значения всех переменных задаются равными определенным значениям, соответствующим параметрам входного потока. Выходная граница полагается свободной.

Если часть границы расчетной области лежит на оси $y = 0$, то продольная компонента $u(x,y)$ вектора скорости $u = (u,v)$ решения задачи для уравнений (11)–(13) является функцией пограничного слоя в окрестности этой части границы относительно переменной $y$. В общем случае структура пограничного слоя для данной задачи неизвестна. Однако некоторую информацию о возможной структуре можно получить из анализа модельной задачи с малым параметром $\varepsilon $:

(14)
$\begin{gathered} - \varepsilon u{\text{''}} + a(y)u{\text{'}} + f(y,u,\varepsilon ) = 0,\quad a(0) = 0,\quad 0 < y < 1, \\ u(0) = 0,\quad u(1) = {{A}_{1}}, \\ \end{gathered} $
являющейся частным случаем задачи (1). Уравнение в (14) получается из стационарного уравнения $({{u}_{t}} = 0)$ для компоненты $u$, которое может быть записано в виде
(15)
$ - \frac{M}{{{\text{Re}}}}\frac{{{{\partial }^{2}}u}}{{\partial {{y}^{2}}}} + \rho v\frac{{\partial u}}{{\partial y}} = g\left( {u,\frac{{\partial u}}{{\partial x}},\frac{{{{\partial }^{2}}u}}{{\partial {{x}^{2}}}},\frac{{\partial v}}{{\partial x}},\frac{{{{\partial }^{2}}v}}{{\partial x\partial y}},x,y} \right),$
если предположить, что $\varepsilon = M{\text{/Re}}$, $x$ – параметр, правая часть равномерно ограничена по $\varepsilon $, $a(y)$ соответствует $\rho v(x,y)$, а $f(y,u,\varepsilon )$ соответствует правой части уравнения (15).

В [1] было показано, что если ${{f}_{u}}(y,u,\varepsilon ) > 0$, то в зависимости от значения производной функции $a(y)$ в точке $y = 0$ решение задачи (14) может иметь экспоненциальный пограничный слой, если $a{\text{'}}(0) = 0$; степенной слой первого типа, если $a{\text{'}}(0) > 0$; смешанный (комбинация экспоненциального и степенного слоя второго типа), если $a{\text{'}}(0) < 0$.

Оценки производных функции $u(x,y)$ по переменной $y$ в пограничном слое для обтекания пластины вязким несжимаемым газом можно было бы получить из автомодельного уравнения $f{\kern 1pt} {\text{'''}} + ff{\kern 1pt} {\text{''}} = 0,\quad u = f{\kern 1pt} {\text{'}}\partial \xi {\text{/}}\partial y,\quad \xi = 1{\text{/}}2y\sqrt {{{u}_{\infty }}{\text{/}}\mu x} $, если было бы известно его аналитическое решение. Однако такое решение неизвестно, а из численных расчетов непонятно, как получить оценки производных.

В данной работе исследуется гиперзвуковое течение газа (число Маха ${\text{M}}$ = 6, 10). Уравнения решались для различных чисел Рейнольдса $\operatorname{Re} $ (10 000, 100 000, 1 000 000). Число Прандтля $\Pr $ предполагается равным 0.72.

Расчеты численного решения производились разработанным в ИТПМ СО РАН (см. [11]) расчетным кодом CFS. При вычислении невязких потоков использовалась схема MP5 (Monotonicity–Preserving, 5th–order), имеющая 5-й порядок на гладких решениях, вязкие члены вычислялись с помощью центральных разностей 4-го порядка на компактном шаблоне (см. [11]). Уравнения интегрировались по времени с помощью схемы Рунге–Кутты 2-го порядка. Расчеты производились до установления стационарного решения (бралось значение времени установления стационарного решения ${{T}_{{{\text{setsol}}}}} = $ 1), шаг $\Delta t$ интегрирования по времени выбирается из условия Куранта–Фридрихса–Леви

$\Delta t\;\leqslant \;\frac{{CFL}}{{\mathop {\max \left[ {\frac{{\left| {{{u}_{{ij}}}} \right| + {{a}_{{ij}}}}}{{\Delta {{x}_{i}}}} + \frac{{\left| {{{v}_{{ij}}}} \right| + {{a}_{{ij}}}}}{{\Delta {{y}_{j}}}} + \frac{{2\gamma \mu M}}{{\rho \operatorname{Re} \Pr }}\left( {\frac{1}{{{{{(\Delta {{x}_{i}})}}^{2}}}} + \frac{1}{{{{{(\Delta {{y}_{j}})}}^{2}}}}} \right)} \right]}\limits_{i,j} }},$
где $a$ – скорость звука, $CFL$ – число Куранта–Фридрихса–Леви, предполагаемое равным 0.4 в данной задаче.

Данный программный код создан для расчета на многопроцессорных ЭВМ и был распараллелен с помощью библиотеки MPI (Message Passage Interface) путем геометрической декомпозиции расчетной области на подобласти, каждая из которых присваивается определенному вычислительному ядру процессора. В расчетах использовалось восемь шестиядерных процессоров Intel Xeon X5670 с тактовой частотой 2932 МГц вычислительного комплекса Информационно-вычислительного центра НГУ (см. [12]).

1.1. Преобразования, генерирующие адаптивные сетки

Для построения сетки вблизи границы $y = 0$ берется бесконечно дифференцируемая функция $\phi (\xi ,\varepsilon )$, $0\;\leqslant \;\xi \;\leqslant \;{{\xi }_{0}}$, а для значений ${{\xi }_{0}}\;\leqslant \;\xi \;\leqslant \;1$ эта функция доопределяется многочленом с соблюдением гладкости класса ${{C}^{l}}[0,\;1]$:

(16)
$y = \Phi (\xi ,\varepsilon ,l,{{\xi }_{0}}) = \left\{ \begin{gathered} {{c}_{1}}\phi (\xi ,\varepsilon ),\quad 0\;\leqslant \;\xi \;\leqslant \;{{\xi }_{0}}, \hfill \\ {{c}_{1}}[\phi ({{\xi }_{0}},\varepsilon ) + \phi {\text{'}}({{\xi }_{0}},\varepsilon )(\xi - {{\xi }_{0}}) + \frac{1}{2}\phi {\text{''}}{\kern 1pt} ({{\xi }_{0}},\varepsilon )({{\xi }_{0}})(\xi - {{\xi }_{0}}{{)}^{2}} + \; \ldots \; \hfill \\ \ldots \; + \;\frac{1}{{l!}}{{\phi }^{{(l)}}}({{\xi }_{0}},\varepsilon )({{\xi }_{0}})(\xi - {{\xi }_{0}}{{)}^{l}} + {{c}_{0}}{{(\xi - {{\xi }_{0}})}^{{l + 1}}}],\quad {{\xi }_{0}}\;\leqslant \;\xi \;\leqslant \;1, \hfill \\ \end{gathered} \right.$
где $l\;\leqslant \;n$, ${{c}_{0}} > 0{\kern 1pt} $, и, например, ${{c}_{0}} = 1$; константа ${{c}_{1}} > 0$ выбирается из условия $\Phi (1,\varepsilon ,l,{{\xi }_{0}}) = 1$. Параметр ${{\xi }_{0}}$ в (16) означает долю длины отрезка, отображаемую в предполагаемую зону пограничного слоя, границей которого условно считается точка, в которой производная решения по физической координате равна константе, не зависящей от $\varepsilon $. В частности, при равномерной сетке по $\xi $ число ${{\xi }_{0}}$ означает долю от общего числа шагов сетки, попадающих в пограничный слой.

1.1.1. Базисные локальные преобразования. Для построения адаптивных сеток, сгущающихся в слое возле границы $y = 0$, в качестве функции $\Phi (1,{\kern 1pt} \varepsilon ,{\kern 1pt} l,{\kern 1pt} {{\xi }_{0}})$ в (16) используются локальные преобразования ${{y}_{i}}(\xi ,\varepsilon ,k,a)$, $i = 1,\;2,\;3,\;4,\;5,$ имеющие следующий вид:

(17)
${{y}_{1}}(\xi ,\varepsilon ,k,a) = - \frac{{{{\varepsilon }^{k}}}}{a}\ln (1 - d\xi ),\quad k > 0,\quad a > 0,\quad 0\;\leqslant \;\xi \;\leqslant \;{{\xi }_{0}},$
(18)
${{y}_{2}}(\xi ,\varepsilon ,k,a) = {{\varepsilon }^{k}}\left[ {{{{(1 - d\xi )}}^{{ - 1/a}}} - 1} \right],\quad k > 0,\quad a > 0,\quad 0\;\leqslant \;\xi \;\leqslant \;{{\xi }_{0}},$
(19)
${{y}_{3}}(\xi ,\varepsilon ,k,a) = {{\left( {{{\varepsilon }^{{ka}}} + d\xi } \right)}^{{1/a}}} - {{\varepsilon }^{k}},\quad k > 0,\quad 1{\text{/}}n > a > 0,\quad 0\;\leqslant \;\xi \;\leqslant \;{{\xi }_{0}},$
(20)
${{y}_{4}}(\xi ,\varepsilon ,k,a) = {{\varepsilon }^{k}}\left[ {{{{(1 + {{\varepsilon }^{{ - k}}})}}^{{a\xi }}} - 1} \right],\quad k > 0,\quad a > 0,\quad 0\;\leqslant \;\xi \;\leqslant \;{{\xi }_{0}},$
(21)
${{y}_{5}}(\xi ,\varepsilon ,k,a) = 2\sigma \xi {\kern 1pt} ,\quad \sigma = \min \left\{ {0.5,(n{\text{/}}a){{\varepsilon }^{k}}\ln N} \right\},\quad 0\;\leqslant \;\xi \;\leqslant \;{{\xi }_{0}}.$

В качестве малого параметра $\varepsilon $ рассматривается $M{\text{/Re}}$.

Параметр $k$ отвечает за масштаб и ширину пограничного слоя, а параметр $a$ влияет на ширину пограничного слоя.

Преобразование ${{y}_{1}}(\xi ,\varepsilon ,k,a)$ при $k = 1$ было предложено в [2], а ${{y}_{5}}(\xi ,\varepsilon ,k,a)$ – в [5] для генерации сеток в экспоненциальных слоях, тогда как ${{y}_{i}}(\xi ,\varepsilon ,k,a)$, $i = 2,\;3,\;4$, были предложены в [13], [14], а ${{y}_{2}}(\xi ,\varepsilon ,k,a)$ при $a = 1$, $k = 1{\text{/}}2$ в [15]. Преобразование ${{y}_{2}}(\xi ,\varepsilon ,k,a)$ применяется для построения сеток в экспоненциальных и степенных слоях первого рода, ${{y}_{3}}(\xi ,\varepsilon ,k,a)$ – в степенных слоях II рода и ${{y}_{4}}(\xi ,\varepsilon ,k,a)$ – в логарифмических слоях.

1.1.2. Глобальные модифицированные преобразования, устраняющие пограничные слои. Локальные преобразования (17)–(21) используются для построения глобальных преобразований, устраняющих пограничные слои (2)–(4), с помощью их масштабирования и склейки с полиномиальной функцией.

Преобразования для экспоненциальных пограничных слоев. Модифицированное преобразование Г.И. Шишкина класса ${{C}^{2}}[0,\;1]$, действующее из $[0,\;1]$ на $[0,\;1]$, полученное из функции (21), имеет вид

(22)
${{y}_{{{\text{Sh}}}}}(\xi ,\varepsilon ,k,a) = \left\{ \begin{gathered} 2\sigma \xi ,\quad 0\;\leqslant \;\xi \;\leqslant \;1{\text{/}}2, \hfill \\ \sigma + 2\sigma (\xi - 1{\text{/}}2) + \omega {{(\xi - 1{\text{/}}2)}^{3}},\quad 1{\text{/}}2\;\leqslant \;\xi \;\leqslant \;1, \hfill \\ \end{gathered} \right.$
где $\sigma = \min \left\{ {0.5,(n{\text{/}}a){{\varepsilon }^{k}}\ln N} \right\}$, число ${{\xi }_{0}} = 1{\text{/}}2$ и $\omega $ выбирается из соотношения ${{y}_{{{\text{Sh}}}}}(1,\varepsilon ,k,a) = 1$, т.е. $\omega = 8(1 - 2\sigma )$.

Преобразование класса ${{C}^{2}}[0,\;1]$, действующее из $[0,\;1]$ на $[0,\;1]$, полученное из локального преобразования Н.С. Бахвалова (17) путем гладкого продолжения за пределом интервала пограничного слоя, устраняющее экспоненциальные особенности масштаба $k$, имеет вид

(23)
${{y}_{{{\text{Bakh}}}}}(\xi ,\varepsilon ,k,a) = \left\{ \begin{gathered} {{c}_{1}}\left[ { - \frac{{{{\varepsilon }^{k}}}}{a}\ln \left( {1 - d\xi } \right)} \right],\quad 0\;\leqslant \;\xi \;\leqslant \;{{\xi }_{0}}, \hfill \\ {{c}_{1}}\left[ {\frac{{k{{\varepsilon }^{k}}}}{{na}}} \right.\ln \left( {{{\varepsilon }^{{ - 1}}}} \right) + \frac{{d{{\varepsilon }^{{k(1 - 1/n)}}}}}{a}(\xi - {{\xi }_{0}}) + \left. {\frac{{{{d}^{2}}{{\varepsilon }^{{k(1 - 2/n)}}}}}{{2a}}{{{(\xi - {{\xi }_{0}})}}^{2}} + {{c}_{0}}{{{(\xi - {{\xi }_{0}})}}^{3}}} \right], \hfill \\ {{\xi }_{0}}\;\leqslant \;\xi \;\leqslant \;1, \hfill \\ \end{gathered} \right.$
где $d = {{\left( {1 - {{\varepsilon }^{{k/n}}}} \right)} \mathord{\left/ {\vphantom {{\left( {1 - {{\varepsilon }^{{k/n}}}} \right)} {{{\xi }_{0}}}}} \right. \kern-0em} {{{\xi }_{0}}}}$, $a > 0$, ${{c}_{0}} > 0$ и константа ${{c}_{1}}$ выбирается из соотношения ${{y}_{{{\text{Bakh}}}}}(1,\varepsilon ,k,a) = 1$, т.е.
$\frac{1}{{{{c}_{1}}}} = \left[ {\frac{{k{{\varepsilon }^{k}}}}{{na}}\ln \left( {{{\varepsilon }^{{ - 1}}}} \right) + \frac{{d{{\varepsilon }^{{k(1 - 1/n)}}}}}{a}(1 - {{\xi }_{0}}) + \frac{{{{d}^{2}}{{\varepsilon }^{{k(1 - 2/n)}}}}}{{2a}}{{{(1 - {{\xi }_{0}})}}^{2}} + {{c}_{0}}{{{(1 - {{\xi }_{0}})}}^{3}}} \right].$
Константа $a$ влияет на ширину пограничного слоя, а константа ${{\xi }_{0}}$ задает количество точек сетки в пограничном слое. В [7] было доказано, что это преобразование устраняет экспоненциальную особенность (2) до порядка $n$ при $0 < a\;\leqslant \;b{\text{/}}{{n}^{2}}$.

Преобразования для экспоненциальных и степенных пограничных слоев. Преобразование класса ${{C}^{2}}[0,\;1]$, действующее из $[0,\;1]$ на $[0,\;1]$, предложенное В.Д. Лисейкиным в [7] и устраняющее в пограничном слое экспоненциальные и степенные I рода особенности масштаба $k$, имеет вид

(24)
${{y}_{{{\text{Lis1}}}}}(\xi ,\varepsilon ,k,a) = \left\{ \begin{gathered} {{c}_{1}}{{\varepsilon }^{k}}\left[ {{{{(1 - d\xi )}}^{{ - 1/a}}} - 1} \right],\quad 0\;\leqslant \;\xi \;\leqslant \;{{\xi }_{0}}, \hfill \\ {{c}_{1}}\left[ {{{\varepsilon }^{{k\beta n}}} - {{\varepsilon }^{k}} + d\frac{1}{a}{{\varepsilon }^{{k\beta (n - 1)}}}(\xi - {{\xi }_{0}}) + \frac{1}{2}{{d}^{2}}\frac{1}{a}\left( {\frac{1}{a} + 1} \right){{\varepsilon }^{{k\beta (n - 2)}}}{{{(\xi - {{\xi }_{0}})}}^{2}} + {{c}_{0}}{{{(\xi - {{\xi }_{0}})}}^{3}}} \right] \hfill \\ {{\xi }_{0}}\;\leqslant \;\xi \;\leqslant \;1, \hfill \\ \end{gathered} \right.,$
где $d = {{\left( {1 - {{\varepsilon }^{{k\beta }}}} \right)} \mathord{\left/ {\vphantom {{\left( {1 - {{\varepsilon }^{{k\beta }}}} \right)} {{{\xi }_{0}}}}} \right. \kern-0em} {{{\xi }_{0}}}}$, при этом $\beta = a{\text{/}}(1 + na)$, $a > 0$, ${{c}_{0}} > 0$ и константа ${{c}_{1}} > 0$ выбирается из условия ${{y}_{{{\text{Lis1}}}}}(1,\varepsilon ,k,a) = 1$, т.е.

$\frac{1}{{{{c}_{1}}}} = \left[ {{{\varepsilon }^{{k\beta n}}} - {{\varepsilon }^{k}} + d\frac{1}{a}{{\varepsilon }^{{k\beta (n - 1)}}}(1 - {{\xi }_{0}}) + \frac{1}{2}{{d}^{2}}\frac{1}{a}\left( {\frac{1}{a} + 1} \right){{\varepsilon }^{{k\beta (n - 2)}}}{{{(1 - {{\xi }_{0}})}}^{2}} + {{c}_{0}}{{{(1 - {{\xi }_{0}})}}^{3}}} \right].$

В [7] было доказано, что это преобразование устраняет экспоненциальную особенность (2) до порядка $n$ при произвольном $a > 0$, а степенную особенность I рода (3) при $0 < a\;\leqslant \;b{\text{/}}{{n}^{2}}$.

Для устранения степенной особенности II рода (4) используется масштабированное преобразование (19):

(25)
${{y}_{3}}(\xi ,\varepsilon ,k,a) = \frac{{{{{\left( {{{\varepsilon }^{{{{k}_{1}}a}}} + d\xi } \right)}}^{{1/a}}} - {{\varepsilon }^{{{{k}_{1}}}}}}}{{{{{\left( {{{\varepsilon }^{{{{k}_{1}}a}}} + d} \right)}}^{{1/a}}} - {{\varepsilon }^{{{{k}_{1}}}}}}},\quad {{k}_{1}} = \frac{k}{{1 - a}},\quad k > 0,$
для которого в [7] было доказано, что это преобразование устраняет до порядка $n$ степенную особенность II рода при $0 < a\;\leqslant \;\min (b{\text{/}}n,1{\text{/}}n).$

Преобразование для логарифмических особенностей. Преобразование класса ${{C}^{2}}[0,\;1]$, действующее из $[0,\;1]$ на $[0,\;1]$, предложенное в [7] для логарифмических особенностей масштаба $k$ (5), имеет вид

(26)
${{y}_{{{\text{Lis2}}}}}(\xi ,\varepsilon ,k) = \left\{ \begin{gathered} {{c}_{1}}{{\varepsilon }^{k}}\left[ {{{{\left( {1 + \frac{1}{{{{\varepsilon }^{k}}\ln ({{\varepsilon }^{{ - k}}})}}} \right)}}^{{\xi /{{\xi }_{0}}}}} - 1} \right],\quad 0\;\leqslant \;\xi \;\leqslant \;{{\xi }_{0}}, \hfill \\ {{c}_{1}}\left[ {{{{\ln }}^{{ - 1}}}\left( {{{\varepsilon }^{{ - k}}}} \right) + 2\left( {{{\varepsilon }^{k}} + \mathop {\ln }\nolimits^{ - 1} \left( {{{\varepsilon }^{{ - k}}}} \right)} \right)\ln \left( {1 + \frac{1}{{{{\varepsilon }^{k}}\ln \left( {{{\varepsilon }^{{ - k}}}} \right)}}} \right)(\xi - {{\xi }_{0}}) + {{c}_{0}}{{{(\xi - {{\xi }_{0}})}}^{2}}} \right] \hfill \\ {{\xi }_{0}}\;\leqslant \;\xi \;\leqslant 1, \hfill \\ \end{gathered} \right.,$
где ${{c}_{0}} > 0$ и константа ${{c}_{1}} > 0$ выбирается из условия ${{y}_{{{\text{Lis2}}}}}(1,\varepsilon ,k) = 1$, т.е.

$\frac{1}{{{{c}_{1}}}} = \left[ {\mathop {\ln }\nolimits^{ - 1} \left( {{{\varepsilon }^{{ - k}}}} \right) + 2\left( {{{\varepsilon }^{k}} + \mathop {\ln }\nolimits^{ - 1} \left( {{{\varepsilon }^{{ - k}}}} \right)} \right)\ln \left( {1 + \frac{1}{{{{\varepsilon }^{k}}\ln \left( {{{\varepsilon }^{{ - k}}}} \right)}}} \right)(1 - {{\xi }_{0}}) + {{c}_{0}}{{{(1 - {{\xi }_{0}})}}^{2}}} \right].$

Преобразование (20) (при не очень малом $\varepsilon $) применяется также в качестве глобального для логарифмических особенностей масштаба $k$:

(27)
${{y}_{4}}(\xi ,\varepsilon ,k,a) = \frac{{{{{\left( {1 + {{\varepsilon }^{{ - k}}}} \right)}}^{{a\xi }}} - 1}}{{{{{\left( {1 + {{\varepsilon }^{{ - k}}}} \right)}}^{a}} - 1}},\quad k > 0,\quad a > 0.$

При этом преобразование (27) имеет отношение к известному преобразованию Эриксона (см. [16]) $y(\xi ,{\kern 1pt} d) = {{\left( {{{e}^{{d\xi }}} - 1} \right)} \mathord{\left/ {\vphantom {{\left( {{{e}^{{d\xi }}} - 1} \right)} {\left( {{{e}^{d}} - 1} \right)}}} \right. \kern-0em} {\left( {{{e}^{d}} - 1} \right)}}$, которое не имеет явную зависимость от $\varepsilon $, но при $d = \ln \left( {1 + {{\varepsilon }^{{ - k}}}} \right)$ совпадает с преобразованием (27) для $a = 1$.

1.2. Исследуемые характеристики численных решений

Расчетная область численного решения представляет собой прямоугольник размером $1 \times 0.32$, нижняя сторона которого совпадает с поверхностью пластины. Расчеты численных решений задачи выполнялись на прямоугольной сетке, состоящей из ${{N}_{x}}_{y}$ ячеек. Сетка по продольной координате $x$ равномерная с фиксированным числом узлов ${{N}_{x}} = $ 192; при построении сетки в нормальном к пластине направлении $y$ использовались исследуемые адаптивные сетки.

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

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

(28)
${{\delta }_{t}} = \mathop {\max }\limits_{i = 1}^{{{N}_{t}}} \left| {u_{{3i}}^{{{{N}_{{t + 1}}}}} - u_{i}^{{{{N}_{t}}}}} \right|,\quad t = 1,\;2,$
где ${{u}^{{{{N}_{t}}}}}$ – приближенное решение на сетке с ${{N}_{t}}$ узлами.

По оценкам погрешности вычислялась оценка реально наблюдаемого порядка точности

(29)
${{p}_{1}} = \mathop {\log }\nolimits_3 ({{\delta }_{1}}{\text{/}}{{\delta }_{2}}).$

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

(30)
$d{{u}_{t}} = \mathop {\max }\limits_{i = 1}^{{{{\hat {N}}}_{t}}} \left| {u_{{i + 1}}^{{{{N}_{t}}}} - u_{i}^{{{{N}_{t}}}}} \right|,\quad t = 1,\;2,\;3,$
где ${{\hat {N}}_{t}}$ – количество узлов сетки ниже положения ударной волны.

По полученным значениям скачка численного решения оценивался порядок скачка решения

(31)
$p_{2}^{t} = \mathop {\log }\nolimits_3 (d{{u}_{t}}{\text{/}}d{{u}_{{t + 1}}}),\quad t = 1,\;2.$

Заметим, что если решение задачи не имеет пограничных и внутренних слоев, то для численного решения этой задачи с использованием схемы порядка $p$ на равномерной сетке значения ${{p}_{1}}$ близки к $p$, а ${{p}_{2}}$ близки к 1.

2. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ

2.1. Выбор параметров и рассмотренные случаи

Во всех рассматриваемых случаях предполагается, что параметр $k$, отвечающий за масштаб пограничного слоя, для данного типа задач равен 0.5 (см. [7]). Параметр $n$, характеризующий порядок схемы, предполагается равным 4.

Для каждого преобразования подбирались оптимальные значения неизвестных параметров. Идея алгоритма поиска заключается в следующем: значения параметров варьировались в интервале их допустимых значений и подбирались так, чтобы их влияние на время расчетов и значение погрешности численного решения было минимальным при изменении малого параметра $\varepsilon = {\text{M/Re}}$, т.е. при изменении чисел M и ${\text{Re}}$.

Для всех предложенных преобразований были произведены расчеты численных решений при параметре Маха ${\text{M}} = 6$ для двух различных чисел Рейнольдса ${\text{Re}}$ (10 000, 100 000).

Для каждого численного эксперимента помимо характеристик (28)–(31) приводится время расчета в часах ${{T}_{{{\text{clock}}}}}$.

2.2. Расчеты с использованием глобальных модифицированных преобразований

Были произведены расчеты с использованием преобразования (25) ${{y}_{3}}(\xi ,\varepsilon ,a,k)$, где $d = 1$ и значение параметра $a$ равно 0.65. Полученные результаты предвставлены в табл. 1 и на фиг. 1, 2. Можно сделать следующие выводы: использование данного преобразования возможно только при больших (>100 000) значениях числа Рейнольдса ${\text{Re}}$, иначе время расчетов ${{T}_{{{\text{clock}}}}}$ превышает 400 ч. Кроме того, большая часть узлов попадает в зону вне пограничного слоя.

Таблица 1.  

Значения погрешности и скачка решения с использованием преобразования (25)

${{N}_{{y,t}}}$ M = 6, Re = 10 000 M = 6, Re = 100 000
${{T}_{{{\text{clock}}}}}$ ${{\delta }_{t}}$ ${{p}_{1}}$ $d{{u}_{t}}$ $p_{2}^{t}$ ${{T}_{{{\text{clock}}}}}$ ${{\delta }_{t}}$ ${{p}_{1}}$ $d{{u}_{t}}$ $p_{2}^{t}$
51 1 – – 0.297 – 1 – – 0.568 –
153 6 0.017 – 0.099 1 4 0.052 – 0.192 0.987
459 420 0.003 1.58 0.033 1 72 0.009 1.6 0.064 1

Примечание. В табл. 1–7 ${{N}_{{y,t}}}$ – количество узлов в поперечном направлении $y$, ${{T}_{{{\text{clock}}}}}$ – время расчета в часах.

Фиг. 1.

Профили продольной скорости и значения в узлах сетки с использованием преобразования (25): (а) – при Re = 10 000, (б) – при Re = 100 000 (M = 6).

Фиг. 2.

Расчетные поля продольной скорости с использованием преобразования (25): (а) – при Re = 10 000, (б) – при Re = 100 000) (M = 6).

Так же были произведены расчеты с использованием преобразования (27) ${{y}_{4}}(\xi ,\varepsilon ,a,k)$ при значении $a = 1$. Полученные результаты представлены в табл. 2 и на фиг. 3, 4. Можно сделать следующие выводы: использование данного преобразования невозможно из-за слишком большого времени расчетов ${{T}_{{{\text{clock}}}}}$, требующегося для достижения установившегося решения.

Таблица 2.  

Значения погрешности и скачка решения с использованием преобразования (27)

${{N}_{{y,t}}}$ M = 6, Re = 10 000 M = 6, Re = 100 000
${{T}_{{{\text{clock}}}}}$ ${{\delta }_{t}}$ ${{p}_{1}}$ $d{{u}_{t}}$ $p_{2}^{t}$ ${{T}_{{{\text{clock}}}}}$ ${{\delta }_{t}}$ ${{p}_{1}}$ $d{{u}_{t}}$ $p_{2}^{t}$
51 1 – – 0.152 – 2 – – 0.25 –
153 15 0.003 – 0.056 0.909 18 0.001 – 0.093 0.9
459 336 0.002 0.3669 0.019 0.984 373 0.001 0 0.033 0.943
Фиг. 3.

Профили продольной скорости и значения в узлах сетки с использованием преобразования (27): (а) – при Re = 10 000, (б) – при Re = 100 000 (M = 6).

Фиг. 4.

Расчетные поля продольной скорости с использованием преобразования (27): (а) – при Re = 10 000, (б) – при Re = 100 000 (M = 6).

Были произведены расчеты с использованием преобразования (22) ${{y}_{{{\text{Sh}}}}}(\xi ,\varepsilon ,a,k)$ при значении параметра $a = 2$. Полученные результаты представлены в табл. 3 и на фиг. 5, 6. Можно сделать следующие выводы: использование данного преобразования возможно при любых значениях числа Рейнольдса ${\text{Re}}$, время расчетов ${{T}_{{{\text{clock}}}}}$ оказалось одним из наименьших среди всех рассмотренных законов сгущения. Недостаток данного преобразования – высокое значение погрешности, которое не уменьшается при увеличении значения числа ${{N}_{y}}$. Кроме того, сгущение узлов осуществляется только на части пограничного слоя, что свидетельствует о том, что пограничный слой не экспоненциального типа.

Таблица 3.  

Значения погрешности и скачка решения с использованием преобразования (22)

${{N}_{{y,t}}}$ M = 6, Re = 10 000 M = 6, Re = 100 000
${{T}_{{{\text{clock}}}}}$ ${{\delta }_{t}}$ ${{p}_{1}}$ $d{{u}_{t}}$ $p_{2}^{t}$ ${{T}_{{{\text{clock}}}}}$ ${{\delta }_{t}}$ ${{p}_{1}}$ $d{{u}_{t}}$ $p_{2}^{t}$
51 1 – – 0.0628 – 1 – – 0.0468 –
153 1 0.406 – 0.0269 0.773 1 0.305 – 0.02 0.774
459 12 0.411 –0.011 0.011 0.813 20 0.319 –0.041 0.0081 0.823
Фиг. 5.

Профили продольной скорости и значения в узлах сетки с использованием преобразования (22): (а) – при Re = 10 000, (б) – при Re = 100 000 (M = 6).

Фиг. 6.

Расчетные поля продольной скорости с использованием преобразования (22): (а) – при ${\text{Re}} = $ 10 000, (б) – при ${\text{Re}} = $ 100 000 (M = 6).

Расчеты с использованием преобразования (23) ${{y}_{{{\text{Bakh}}}}}(\xi ,\varepsilon ,k,a)$ оказались безуспешными, так как ни варьирование параметра $a$, отвечающего за ширину пограничного слоя, ни варьирование параметра ${{\xi }_{0}}$, отвечающего за долю длины отрезка, отображаемую в предполагаемую зону пограничного слоя, не дало таких значений параметров, при которых расчеты численных решений производились бы за “разумное” время – данное преобразование чрезмерно сильно сгущает расчетную сетку вблизи точки $\xi = 0$, вследствие чего шаг по времени слишком мал. Данное преобразование не подходит для расчета численных решений задач с таким типом слоев.

Были произведены расчеты с использованием преобразования (26) ${{y}_{{{\text{Lis2}}}}}(\xi ,\varepsilon ,k)$ при значении параметра ${{\xi }_{0}} = 0.5$. Полученные результаты представлены в табл. 4 и на фиг. 7, 8. Можно сделать следующие выводы: использование данного преобразования возможно при любых значениях числа Рейнольдса ${\text{Re}}$, но время расчетов ${{T}_{{{\text{clock}}}}}$ превышает 160 ч при значении числа ${{N}_{y}} = $ 459. Кроме того, сгущение узлов происходит на части пограничного слоя.

Таблица 4.  

Значения погрешности и скачка решения с использованием преобразования (26)

${{N}_{{y,t}}}$ M = 6, Re = 10 000 M = 6, Re = 100 000
${{T}_{{{\text{clock}}}}}$ ${{\delta }_{t}}$ ${{p}_{1}}$ $d{{u}_{t}}$ $p_{2}^{t}$ ${{T}_{{{\text{clock}}}}}$ ${{\delta }_{t}}$ ${{p}_{1}}$ $d{{u}_{t}}$ $p_{2}^{t}$
51 1 – – 0.36 – 1 – – 0.76 –
153 8 0.009 – 0.133 0.906 9 0.056 – 0.237 1.06
459 204 0.001 2 0.046 0.966 168 0.003 2.66 0.079 1
Фиг. 7.

Профили продольной скорости и значения в узлах сетки с использованием преобразования (26): (а) – при Re = 10 000, (б) – при Re = 100 000 (M = 6).

Фиг. 8.

Расчетные поля продольной скорости с использованием преобразования (26): (а) – при ${\text{Re}} = $ 10 000, (б) – при Re = 100 000 (M = 6).

Были произведены расчеты с использованием преобразования (24) ${{y}_{{{\text{Lis1}}}}}(\xi ,\varepsilon ,a,k)$ при значении параметров ${{\xi }_{0}}$,  равном 0.8 и $a = 2$. Полученные результаты представлены в табл. 5 и на фиг. 9, 10. Можно сделать следующие выводы: использование данного преобразования возможно при любых значениях числа Рейнольдса ${\text{Re}}$, время расчетов ${{T}_{{{\text{clock}}}}}$ оказалось самым наименьшим среди всех рассмотренных законов сгущения. Так же отметим малое значение погрешности, которое стремительно уменьшается при увеличении значения числа ${{N}_{y}}$.

Таблица 5.  

Значения погрешности и скачка решения с использованием преобразования (24)

${{N}_{{y,t}}}$ M = 6, Re = 10 000 M = 6, Re = 100 000
${{T}_{{{\text{clock}}}}}$ ${{\delta }_{t}}$ ${{p}_{1}}$ $d{{u}_{t}}$ $p_{2}^{t}$ ${{T}_{{{\text{clock}}}}}$ ${{\delta }_{t}}$ ${{p}_{1}}$ $d{{u}_{t}}$ $p_{2}^{t}$
51 1 – – 0.15 – 1 – – 0.213 –
153 4 0.005 – 0.052 0.964 3 0.017 – 0.074 0.962
459 18 0.002 0.834 0.018 0.966 12 0.003 1.58 0.025 0.988
Фиг. 9.

Профили продольной скорости и значения в узлах сетки с использованием преобразования (24): (а) – при ${\text{Re}} = $ 10 000, (б) – при Re = 100 000 (M = 6).

Фиг. 10.

Расчетные поля продольной скорости с использованием преобразования (24): (а) – при Re = 10 000, (б) – при Re = 100 000 (M = 6).

2.3. Расчеты для других чисел ${\text{M}}$ и ${\text{Re}}$ с использованием преобразования для экспоненциальных и степенных особенностей

Проанализировав результаты численных экспериментов с использованием каждого из рассматриваемых законов сгущения для числа Маха ${\text{M}} = 6$ и чисел Рейнольдса 10 000 и 100 000, было выявлено, что преобразование (24) ${{y}_{{{\text{Lis1}}}}}(\xi ,\varepsilon ,a,k)$ отличается от всех других глобальных преобразований не только значительно меньшим временем расчетов, но и лучшими значениями исследуемых численных характеристик, а также лучшим распределением узлов в зоне пограничного слоя (так как узлы не только наиболее распределены в зоне пограничного слоя, но и также более равномерно) по сравнению с другими сетками. Поэтому были произведены расчеты с использованием данного преобразования при тех же значениях параметров ${{\xi }_{0}}$ и $b$, но при других значениях чисел Маха M и Рейнольдса ${\text{Re}}$, чтобы проанализировать поведение закона сгущения сеток при различных числах M и ${\text{Re}}$.

Произведены расчеты с использованием преобразования (24) ${{y}_{{{\text{Lis1}}}}}(\xi ,\varepsilon ,a,k)$ при значении числа Маха ${\text{M}} = 6$ и числа Рейнольдса ${\text{Re}} = 1\,000\,000$. Полученные результаты представлены в табл. 6 и на фиг. 11, 12. Можно сделать следующие выводы: время расчетов ${{T}_{{{\text{clock}}}}}$ не превысило 10 ч, а значение погрешности также осталось невелико.

Таблица 6.  

Значения погрешности и скачка решения с использованием преобразования (24)

${{N}_{{y,t}}}$ M = 6, Re = 1 000 000
${{T}_{{{\text{clock}}}}}$ ${{\delta }_{t}}$ ${{p}_{1}}$ $d{{u}_{t}}$ $p_{2}^{t}$
51 1 – – 0.283 –
153 2 0.064 – 0.099 0.956
459 10 0.005 2.32 0.034 0.973
Фиг. 11.

Профили продольной скорости и значения в узлах сетки с использованием преобразования (24) при Re = 1 000 000 и M = 6.

Фиг. 12.

Расчетное поле продольной скорости с использованием преобразования (24) при Re = 1 000 000 и M = 6.

Также были произведены расчеты с использованием преобразования (24) ${{y}_{{{\text{Lis1}}}}}(\xi ,\varepsilon ,b,k)$ при значении чисел Маха ${\text{M}} = 10$ и 15 и числа ${\text{Re}} = 100\,000$. Полученные результаты представлены в табл. 7 и на фиг. 13, 14. Можно сделать следующие выводы: время расчетов ${{T}_{{{\text{clock}}}}}$ значительно не увеличилось после увеличения числа Маха ${\text{M}}$, и значение погрешности так же осталось невелико.

Таблица 7.  

Значения погрешности и скачка решения с использованием преобразования (24)

${{N}_{{y,t}}}$ M = 10, Re = 100 000 M = 15, Re = 100 000
${{T}_{{{\text{clock}}}}}$ ${{\delta }_{t}}$ ${{p}_{1}}$ $d{{u}_{t}}$ $p_{2}^{t}$ ${{T}_{{{\text{clock}}}}}$ ${{\delta }_{t}}$ ${{p}_{1}}$ $d{{u}_{t}}$ $p_{2}^{t}$
51 1 – – 0.313 – 1 – – 0.348 –
153 2 0.018 – 0.108 0.969 3 0.022 – 0.12 0.969
459 18 0.004 1.37 0.037 0.975 57 0.008 0.921 0.041 0.978
Фиг. 13.

Профили продольной скорости и значения в узлах сетки с использованием преобразования (24): (а) – при M = 10, (б) – при M = 15 (Re = 100 000).

Фиг. 14.

Расчетные поля продольной скорости с использованием преобразования (24): (а) – при M = 10, (б) – при ${\text{M}} = $ 15 (Re = 100 000).

ЗАКЛЮЧЕНИЕ

В ходе исследования были проанализированы базисные преобразования для генерации адаптивных сеток и преобразования Н.С. Бахвалова, В.Д. Лисейкина и Г.И. Шишкина с оптимальными значениями параметров для построения сеток в пограничных слоях. В серии численных экспериментов для задачи сверхзвукового обтекания плоской пластины вязким газом были исследованы оценки погрешности и порядка точности численных решений, оценки значений скачка решения в близлежащих узлах и порядок скачка решения, а также рассмотрено время расчетов с использованием данных сеток. Показано, что преобразование (24), сконструированное для численных расчетов задач с экспоненциальными и степенными слоями, является наиболее эффективным для исследования сверхзвукового обтекания плоской пластины вязким газом и может использоваться при численных расчетах задач такого типа для достижения качественных результатов с наименьшими затратами времени.

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

  1. Liseikin V.D. Layer Resolving Grids and Transformations for Singular Perturbation Problems. VSP, Utrecht, 2001.

  2. Бахвалов Н.С. Об оптимизации методов решения краевых задач при наличии пограничного слоя // Ж. вычисл. матем. и матем. физ. 1969. Т. 9. № 4. С. 841–859.

  3. Berger A.E., Han H., Kellog R.B. A priori estimates of a numerical method for a turning point problem // Math. Comput. 1984. V. 42. № 166. P. 465–492.

  4. Farrell P.A. A uniformly convergent difference scheme for turning point problems. In Miller J.J.H. ed. Boundary and Interior Layers, Computational and Asymptotic Methods, Boole Press, Dublin, 1980. P. 270–274.

  5. Шишкин Г.И. Разностная схема для сингулярно возмущенного уравнения параболического типа с разрывным начальным условием // Ж. вычисл. матем. и матем. физ. 1988. Т. 28. № 11. С. 1649–1662.

  6. Liseikin V.D., Karasuljic S. Numerical analysis of grid-clustering rules for problems with power of the first type boundary layers // Comput. Technolog. 2020. V. 25. № 1. P. 49–65.

  7. Liseikin V.D. Grid generation for problems with boundary and interior layers // Novosibirsk: NGU, 2018.

  8. Liseikin V.D., Paasonen V.I. Compact difference schemes and layer-resolving grids for numerical modeling of problems with boundary and interior layers // Numer. Anal. and Appl. 2019. V. 12. № 1. P. 1–17.

  9. Liseikin V.D. Grid Generation Methods. Springer, third ed. Berlin, 2017.

  10. Белокрыс-Федотов А.И., Гаранжа В.А., Кудрявцева Л.Н. Построение сеток Делоне в неявных областях с обострением ребер // Ж. вычисл. матем. и матем. физ. 2016. Т. 56. № 11. С. 1931–1948.

  11. Кудрявцев А.Н., Поплавская Т.В., Хотяновский Д.В. Применение схем высокого порядка точности при моделировании нестационарных сверхзвуковых течений // Матем. моделирование. 2007. Т. 19. № 7. С. 39–55.

  12. Информационно-вычислительный центр Новосибирского государственного университета. [Электронный ресурс]. URL: http://nusc.nsu.ru/ (дата обращения: 25.01.2021).

  13. Лисейкин В.Д. О численном решении уравнений со степенным пограничным слоем // Ж. вычисл. матем. и матем. физ. 1986. Т. 26. № 12. С. 1813–1820.

  14. Лисейкин В.Д. О численном решении сингулярно возмущенных уравнений с точками поворота // Ж. вычисл. матем. и матем. физ. 1984. Т. 24. № 12. С. 1812–1818.

  15. Vulanovic R. Mesh construction for numerical solution of a type of singular perturbation problems // Numer. Meth. Approx. Theory. 1984. P. 137–142.

  16. Eriksson L.E. Generation of boundary-conforming grids around wing-body configurations using transfinite interpolation // AIAA J. 1982. V. 20. P. 1313–1320.

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

Инструменты

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

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

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