navierstokes.online · справка

Постановка, допущения и методы расчёта

Страница описывает обе основные задачи: нестационарный запуск плоского сопла Лаваля и обтекание проволочной сетки в канале. Для каждой приведены входные данные, уравнения, граничные условия, численный метод и формулы для выводимых величин. Обозначения совпадают с заданиями.

Часть I

Сопло Лаваля: запуск и сравнение со стационаром

Плоское (двумерное) сопло «Контур 4». В начальный момент газ во всём сопле покоится. Справа за срезом CD давление понижено до \(P_{03}\), и от этого сопло запускается. Нужно проследить, как в горле возникает и развивается сверхзвуковая зона: сначала у стенок (режим Тейлора), затем через всё горло (режим Мейера), как в работе Крайко.

Геометрия

Стенки — дуги окружностей радиуса \(R\), симметричные относительно оси \(y=0\). Полувысота канала:

\[ \begin{gathered} y_w(x) \;=\; y_m + R - \sqrt{R^2 - x^2},\\ R = 200\ \text{мм},\quad y_m = \tfrac12 h_{кр} = 5\ \text{мм},\quad x\in[-80,\,80]\ \text{мм}. \end{gathered} \]

Горло находится в \(x=0\) и имеет высоту \(h_{кр} = 10\) мм. На концах \(y_w(\pm 80) = 21{,}70\) мм, поэтому геометрическое отношение площадей (на единицу ширины) равно

\[ \frac{A_e}{A_*} = \frac{y_w(80)}{y_m} \approx 4{,}34 . \]

Расчётная область — только само сопло: вход AB при \(x=-80\) мм, выход CD при \(x=+80\) мм, стенки AD и BC.

Входные параметры

ВеличинаЗначениеГде задаётся
Газвоздух, \(\gamma = 1{,}4\), \(R_{г} = 287{,}05\) Дж/(кг·К)константы решателя
Слой ячеек слева от AB (зона 1)\(P_{01} = 1\) атм, \(T_0 = 300\) K, покой, постояннофиксировано
В сопле при \(t=0\) (область 2)\(P_{02} = 1\) атм, \(T_0 = 300\) K, \(\mathbf u = 0\)фиксировано
Слой ячеек справа от CD (зона 3)противодавление \(P_{03}\), \(T_0 = 300\) K, покой, постояннополе \(P_{03}\) (точность \(10^{-6}\) атм) или поле \(n\); кнопки режимов 1…10 и значений \(n\) из табл. 1 Крайко
Сеткаот 180×64 до 640×256 ячеек или своя \(N_x\times N_y\); сгущение к стенкам \(\beta_y\) и к горлу \(\beta_x\)списки «Сетка», «Сгущение»
Выдача результатовнепрерывно или через \(\Delta t^*\); конечное время \(t^*_{кон}\)раздел «Выдача результатов»

Противодавление задаётся от расчёта к расчёту точно: в поле \(P_{03}\) (атм) или в поле \(n = P_{01}/P_{03}\). Если последним изменено \(n\), расчёт берёт \(P_{03} = 1/n\) без округления. Кнопки ниже подставляют значения из таблицы задания; режимы идут с шагом \(\Delta p = 0{,}0472\) атм:

Режим12345678910
\(P_{03}\), атм0,95280,90560,85840,81120,76400,71680,66960,62240,57520,5280
\(n\)1,0501,1041,1651,2331,3091,3951,4931,6071,7391,894

При \(n = 1{,}395\) сопло у Крайко запирается дважды. Для сравнения с табл. 1 статьи есть кнопки \(n = 1{,}35;\ 1{,}38;\ 1{,}385;\ 1{,}39;\ 1{,}395;\ 1{,}40\). Контур у Крайко другой (кривые Безье, \(x\in[-17,4]\,y_m\)), и снаружи сопла у него большая область, поэтому пороги по \(n\) здесь не обязаны совпадать; см. раздел о слое за CD.

Уравнения и допущения

Решаются двумерные нестационарные уравнения Эйлера в консервативной форме:

\[ \frac{\partial \mathbf U}{\partial t} + \frac{\partial \mathbf F}{\partial x} + \frac{\partial \mathbf G}{\partial y} = 0, \] \[ \mathbf U = \begin{pmatrix}\rho\\ \rho u\\ \rho v\\ E\end{pmatrix},\; \mathbf F = \begin{pmatrix}\rho u\\ \rho u^2 + p\\ \rho u v\\ (E+p)\,u\end{pmatrix},\; \mathbf G = \begin{pmatrix}\rho v\\ \rho u v\\ \rho v^2 + p\\ (E+p)\,v\end{pmatrix}, \]

замкнутые уравнением состояния совершенного газа:

\[ E = \frac{p}{\gamma-1} + \frac{\rho\,(u^2+v^2)}{2},\qquad p = \rho R_{г} T,\qquad a = \sqrt{\gamma p/\rho},\qquad M = \frac{|\mathbf u|}{a}. \]

Допущения

  • Невязкий газ. Вязкости и теплопроводности нет, пограничных слоёв нет, стенки скользящие. Расчётное время — десятки микросекунд, за это время вязкий слой у стенки остаётся тонким по сравнению с \(y_m\).
  • Совершенный газ с постоянными \(\gamma\) и \(c_p\).
  • Плоское течение: глубина сопла по \(z\) постоянна, а сопло симметрично относительно оси. Симметрию решатель не навязывает: считается весь канал, а отличие поля от зеркального выводится как «асимметрия M верх/низ» и служит контролем точности.
  • Резервуары бесконечной ёмкости: снаружи AB и CD газ покоится при заданных давлении и 300 K.

Характерный масштаб времени — время пробега звука через полувысоту горла:

\[ t_* = \frac{y_m}{a_0},\qquad a_0 = \sqrt{\gamma R_{г} T_0} = 347{,}2\ \text{м/с},\qquad t_* \approx 14{,}4\ \text{мкс}. \]

На экране время выводится в микросекундах и в долях \(t/t_*\).

Граничные условия

Стенки AD и BC — непротекание

В фиктивной ячейке за стенкой нормальная составляющая скорости отражается:

\[ u_n^{\,g} = -u_n,\qquad u_\tau^{\,g} = u_\tau,\qquad \rho^{g} = \rho,\quad p^{g} = p . \]

В итоге поток массы через стенку равен нулю, а касательная скорость не тормозится (скольжение).

Открытые границы AB и CD — втекание и вытекание

К расчётной области примыкает слой фиктивных ячеек. Его состояние задаётся при \(t=0\) и остаётся постоянным всё время расчёта:

\[ \text{слева от AB:}\;\; (p,\,T,\,\mathbf u) = (P_{01},\,300\ \text{K},\,0),\qquad \text{справа от CD:}\;\; (p,\,T,\,\mathbf u) = (P_{03},\,300\ \text{K},\,0). \]

Поток через грань AB или CD считается тем же решателем задачи Римана HLLC, что и внутри сопла, между постоянным состоянием слоя и ячейкой сопла:

\[ \mathbf F_{AB} = \mathbf F^{HLLC}\big(\mathbf U_1,\;\mathbf U_{1,j}\big),\qquad \mathbf F_{CD} = \mathbf F^{HLLC}\big(\mathbf U_{N_x,j},\;\mathbf U_3\big). \]

Направление течения через границу заранее не задаётся: решение задачи Римана само даёт втекание или вытекание. Волны, приходящие к границе изнутри, уходят в слой, который остаётся неизменным. Так моделируется бесконечный резервуар с постоянными параметрами.

Начальные условия: во всех ячейках сопла (зона 2) \(\mathbf u = 0\), \(p = 1\) атм, \(T = 300\) K. При \(t < 0\) газ в сопле неподвижен, потому что граница CD закрыта мембраной. В момент \(t=0\) мембрана мгновенно убирается. На CD возникает задача о распаде разрыва \(P_{02} \to P_{03}\): внутрь сопла уходит волна разрежения, в резервуар 3 — ударная волна и контактный разрыв. Газ начинает разгоняться к горлу, и сопло запускается.

Задача моделирует запуск реального сопла, присоединённого по AB к бесконечному резервуару с постоянными параметрами и вытекающего через CD в бесконечный резервуар с параметрами из таблицы.

Численный метод

ЭлементРеализация
СеткаКриволинейная, «подогнанная под контур». Узлы по \(x\) сгущены к горлу, \(x_i = 0{,}08\,\sinh\!\big(\beta_x(2i/N_x-1)\big)/\sinh\beta_x\) м, \(\beta_x \in \{0;\,1{,}2;\,2{,}15;\,3\}\), по умолчанию 2,15. Узлы по \(y\) сгущены к стенкам, \(y_{ij} = y_w(x_i)\,\eta_j\), \(\eta_j = \tanh\!\big(\beta_y(2j/N_y-1)\big)/\tanh\beta_y\), \(\beta_y \in \{0;\,1{,}05;\,1{,}6;\,2{,}2;\,3\}\), по умолчанию 1,05; \(\beta = 0\) — равномерная сетка. Вертикальные грани, наклонные грани повторяют стенку.
ПространствоМетод конечных объёмов. Поток через грань — приближённый решатель задачи Римана HLLC. Реконструкция примитивных переменных \((\rho,u,v,p)\) — MUSCL с ограничителем minmod, второй порядок.
ВремяДвухстадийный метод Рунге–Кутты SSP-RK2 (Хойна), второй порядок, сохраняет монотонность.
Шаг по времениНаправленное условие Куранта, число Куранта 0,45: \[\Delta t = \mathrm{CFL}\,\min_{\text{ячейки}}\Big[\frac{|u|+a}{\Delta x} + \frac{|v|+a}{\Delta y}\Big]^{-1}.\] У стенок ячейки тонкие, поэтому \(\Delta t\) определяется наименьшим размером ячейки.
КонтрольНевязка \(\|\rho^{n+1}-\rho^n\|_2/\bar\rho\), асимметрия M между верхней и нижней половинами, нижние пороги \(\rho\) и \(p\) для устойчивости.

HLLC решает задачу Римана с тремя волнами: скорости \(S_L, S_R\) — оценки крайних акустических волн, \(S_*\) — контактный разрыв:

\[ S_L = \min(u_{nL}-a_L,\;u_{nR}-a_R),\qquad S_R = \max(u_{nL}+a_L,\;u_{nR}+a_R), \] \[ S_* = \frac{p_R - p_L + \rho_L u_{nL}(S_L-u_{nL}) - \rho_R u_{nR}(S_R-u_{nR})}{\rho_L(S_L-u_{nL}) - \rho_R(S_R-u_{nR})}. \]

Одномерная оценка режимов (для сравнения)

Квазиодномерная изэнтропическая теория для \(A_e/A_* = 4{,}34\) даёт три характерных перепада:

\[ \frac{A}{A_*} = \frac1M\Big[\tfrac{2}{\gamma+1}\Big(1+\tfrac{\gamma-1}{2}M^2\Big)\Big]^{\frac{\gamma+1}{2(\gamma-1)}},\qquad \frac{p}{P_0} = \Big(1+\tfrac{\gamma-1}{2}M^2\Big)^{-\frac{\gamma}{\gamma-1}} . \]
Порог\(M_e\)\(p_e/P_0\)\(n = P_0/p_e\)
Первый критический: \(M=1\) в горле, дальше дозвук0,1350,9871,013
Прямой скачок точно на срезе3,03 → скачок0,2753,63
Расчётный режим (без скачков)3,030,02638,2

По одномерной теории все режимы таблицы (\(n = 1{,}05 \ldots 1{,}89\)) должны запирать горло и давать скачок в расширяющейся части. Двумерный расчёт показывает, как это происходит на самом деле. Звуковая линия искривлена: \(M=1\) сначала достигается у стенок, где стенка горла выпукла, и только потом звуковые линии смыкаются на оси. Из-за этой кривизны пороги запирания в 2D отличаются от 1D-оценки.

Что показывает расчёт

В одном решателе два начальных состояния (список «Задача»):

Отдельные страницы «1D» и «2D стационар» в меню «Другие примеры» остаются быстрыми демонстрациями.

Выдача результатов и рисунки

Выдача. «Непрерывно» — картинка обновляется после каждой пачки шагов. «Через интервал \(\Delta t^*\)» — расчёт останавливается точно в моменты \(t^* = \Delta t^*,\,2\Delta t^*,\ldots\): последний шаг перед моментом выдачи укорачивается, \(\Delta t = \min(\Delta t_{CFL},\,t_{выд}-t)\), так что поля относятся ровно к заданному времени. Так же точно расчёт останавливается в \(t^*_{кон}\). В момент выдачи можно:

Звуковые линии запоминаются и в моменты событий запуска (\(t_1, t_2, \ldots\)): решатель прерывает пачку шагов сразу после события. Включив «Пауза на событиях», можно остановиться точно на событии. Кнопка «+ ЗЛ сейчас» запоминает текущую линию. Флажок «Не очищать ЗЛ при новом расчёте» позволяет наложить линии разных сеток или разных \(n\), как на рис. 3 у Крайко.

Рисунок предназначен для анализа и публикации. На нём всегда есть:

По желанию добавляются блок со звуковой линией, сводкой, сеткой и событиями, а также сама расчётная сетка: в лупе видна каждая ячейка, а на всём сопле — каждая \(k\)-я линия, и это подписано.

Звуковая линия: измеряемые величины

Звуковая линия строится по числу Маха в центрах ячеек (метод квадратов с линейной интерполяцией) и доводится до стенки. На рисунке она чёрная с белой обводкой и на всём сопле, и в лупе — одной и той же толщины. У скользящей стенки первые три строки сетки численно загибают ЗЛ к нормали (дефект схемы, а не течения — он не уменьшается при измельчении \(\Delta y_1\)). Поэтому на рисунке эти три строки для M = 1 заменяются линейной экстраполяцией с более глубоких строк, а угол к стенке считается по той же очищенной окрестности. Лупа «у верхней / нижней стенки» центрируется на точке касания.

Для величин ниже в окне \(|x| < 5\,y_m\) на каждой строке сетки ищется точка, где M впервые переходит через 1 при движении вниз по потоку: это левая (входная) ветвь звуковой линии.

События запуска (как табл. 1 Крайко)

ОбозначениеЧто фиксируется
\(t_1\)M = 1 впервые на стенке в окне горла — зарождение сверхзвуковых зон (режим Тейлора).
\(y<0{,}1\)Головка ЗЛ впервые ниже \(0{,}1\,y_m\) — момент, для которого у Крайко сравниваются сетки на рис. 3.
\(t_2,\ y_{m2}\)Головка ЗЛ дошла до оси (\(y_{m2}=0\), ЗЛ делится на левую и правую ветви, сопло заперто — режим Мейера) или, не дойдя, пошла назад (\(y_{m2}>0\) — наименьшая ордината головки; поворот засчитывается, если головка поднялась на 0,02 \(y_m\) и держится выше 2 \(t_*\), или если СЗЗ исчезли совсем).
СЗЗ исчезли / сноваСверхзвуковых зон нет нигде дольше 1 \(t_*\) (течение снова всюду дозвуковое) и момент их повторного появления. Бывает около критического \(P_{03}\), когда первая волна разгона лишь кратко доводит M до 1 у стенки.
\(t_3,\ y_{m3}\)Первое запирание прекратилось: на оси в окне горла нет M ≥ 1 не меньше 1 \(t_*\); \(y_{m3}\) — ордината головки в этот момент.
\(t_4,\ y_{m4}\)Повторное запирание (\(y_{m4}=0\)); в журнале указан наибольший отход головки от оси между \(t_3\) и \(t_4\).
\(t_\infty,\ y_\infty,\ \Delta x_\infty\)СЗЗ перестала меняться: за последние 25 \(t_*\) головка ЗЛ изменилась не больше чем на 0,01 \(y_m\), длина СЗЗ на оси — на 0,02 \(y_m\), площадь СЗЗ в окне горла — на 2 %, запертость не менялась. \(t_\infty\) — начало этого интервала. Если потом СЗЗ снова заметно меняется, в журнале появляется «установление нарушено».

Строку сводки текущего расчёта можно добавить в «серию». Так от расчёта к расчёту с разными \(n\) или сетками собирается своя таблица 1, её можно сохранить в CSV. Серия хранится в браузере.

Сетка и её паспорт

Для публикации на панели и в паспорте расчёта приводятся:

На сетке 240×96 при \(\beta_y = 1{,}05\) у стенки в горле \(\Delta y_1 \approx 0{,}011\,y_m\), и в слое \(0{,}1\,y_m\) лежат 8 строк. При \(\beta_y = 2{,}2\) — \(\Delta y_1 \approx 0{,}0024\,y_m\) и 18 строк, на 480×192 — \(0{,}0012\,y_m\) и 36 строк. Шаг по времени падает примерно пропорционально \(\Delta y_1\).

Сохранение результатов

Слой за CD и сравнение с Крайко

Эффективный перепад меньше заданного

Слой неподвижного газа с давлением \(P_{03}\) примыкает прямо к срезу CD. Когда газ вытекает со скоростью \(u_e\), решение задачи Римана на грани CD даёт давление выше \(P_{03}\) примерно на акустический импеданс: \(p_e \approx P_{03} + \rho_3 a_3 u_e\). Пример: при \(n = 1{,}395\) (\(P_{03} = 0{,}7168\) атм) на сетке 180×64 после установления \(u_e \approx 40\) м/с и \(p_e \approx 0{,}838\) атм. Значит, фактический перепад на сопле около 1,19, а не 1,395, и горло не запирается. У Крайко такие же «примыкающие» слои стоят на внешней границе большой области за срезом (рис. 1 статьи), и струя успевает расшириться до \(P_{03}\), поэтому его пороги \(n\) нельзя переносить на эту постановку напрямую. Запирание здесь наступает при больших \(n\): например, при \(n = 1{,}894\) (режим 10) — \(t_1 \approx 59{,}6\), \(t_2 \approx 60{,}1\), \(t_\infty \approx 133\) на сетке 180×64. Этот режим близок к критическому, и ранняя стадия зависит от сетки: около \(t_* \approx 38\)–42 первая волна разгона доводит M у стенки почти до 1. На 180×64 порог не достигается, на 240×96 СЗЗ появляются и сразу исчезают (\(t_1 \approx 42{,}0\)), на 360×128 (\(\beta_y = 2{,}2\)) успевает сложиться эпизод Тейлора (\(t_1 \approx 37{,}8\), поворот головки при \(y_{m2} \approx 0{,}26\)). Запирание наступает на всех трёх сетках в \(t_* \approx 60{,}1\). Угол ЗЛ со стенкой после установления: 62,4° (180×64, \(\beta_y = 1{,}05\)) и 61,8° (\(\beta_y = 2{,}2\)).

Часть II

Обтекание проволочной сетки в канале

Канал квадратного сечения 100×100 мм. Поперёк него под углом \(\alpha\) к оси \(x\) натянута сетка из проволоки. Нужно найти, насколько падает полное давление, и получить профили скорости за сеткой в сечении \(x = 525\) мм, в 100 мм за концом сетки.

Геометрия

\(\alpha\)\(a\), мм\(L\), мм\(S_c\), м²Прутков в срезе
20°150,25274,750,0028654
40°305,80119,200,0019129
60°367,3057,700,0016421
80°407,4017,600,0015419

Таблица 1 задания. \(S_c\) — площадь полотна сетки в половине канала.

Входные параметры и режимы

ВеличинаЗначение
Вход (\(x=0\))полное давление \(P^*_о = 101325\) Па, полная температура \(T^* = 288\) K, направление — вдоль \(x\)
Выходмассовый расход \(G\) по таблице 2
Стенки и пруткиприлипание \(\mathbf u = 0\), адиабатические: \(\partial T/\partial n = 0\)
Газвоздух, \(\gamma = 1{,}4\), \(R_г = 287\), вязкость по Сазерленду
Критическая скорость\(a_{кр} = \sqrt{\frac{2\gamma}{\gamma+1}R_г T^*} = 310{,}55\) м/с
\(\lambda_о\)\(\rho_о\), кг/м³\(G\), кг/с\(M_о\)\(p_о\), Па\(q_о\), Па\(\mathrm{Re}_d\)
0,21,2060,3740,18398 9802 325≈3 800
0,31,1810,5500,27696 1045 122≈5 500
0,41,1460,7120,37092 1798 839≈7 200
0,51,1020,8560,46687 30213 285≈8 600
0,61,0500,9780,56581 59518 229≈9 800

\(\rho_о\) и \(G\) — таблица 2 задания. Остальные столбцы вычислены по формулам ниже. \(\mathrm{Re}_d = \rho_о u_о d/\mu\) рассчитано по набегающему потоку; в просветах сетки скорость выше. Обязательны режимы \(\lambda_о = 0{,}2\ldots0{,}5\) при всех четырёх углах, всего 16 расчётов.

Газодинамические функции

Через приведённую скорость \(\lambda = V/a_{кр}\) записываются все изэнтропические соотношения:

\[ \tau(\lambda) = \frac{T}{T^*} = 1 - \frac{\gamma-1}{\gamma+1}\lambda^2,\qquad \pi(\lambda) = \frac{p}{P^*} = \tau^{\frac{\gamma}{\gamma-1}}, \] \[ \varepsilon(\lambda) = \frac{\rho}{\rho^*} = \tau^{\frac{1}{\gamma-1}},\qquad M^2 = \frac{2}{\gamma+1}\,\frac{\lambda^2}{\tau(\lambda)} . \]

Отсюда состояние на входе и расход из таблицы 2:

\[ \rho^* = \frac{P^*_о}{R_г T^*} = 1{,}2257\ \text{кг/м}^3, \] \[ \rho_о = \rho^*\varepsilon(\lambda_о),\qquad u_о = a_{кр}\lambda_о,\qquad G = \rho_о\,u_о\,S_о . \]

Например, при \(\lambda_о = 0{,}3\): \(\varepsilon = 0{,}9629\), \(\rho_о = 1{,}180\), \(G = 1{,}180\cdot 93{,}17\cdot 0{,}005 = 0{,}550\) кг/с, что совпадает с таблицей.

Искомые величины

В каждой точке выходного сечения \(x = 525\) мм по статическому давлению \(p\), плотности \(\rho\) и скорости \(V\) вычисляются

\[ T = \frac{p}{\rho R_г},\qquad M^2 = \frac{V^2}{\gamma R_г T}, \] \[ P^* = p\Big(1+\tfrac{\gamma-1}{2}M^2\Big)^{\frac{\gamma}{\gamma-1}},\qquad \sigma(y) = \frac{P^*(y)}{P^*_о},\qquad \lambda_2(y) = \frac{V(y)}{a_{кр}} . \]

Средние по высоте:

\[ \bar\sigma_{100} = \frac{1}{100}\int_0^{100}\sigma\,dy,\qquad \bar\sigma_{50} = \frac{1}{50}\int_{25}^{75}\sigma\,dy, \]

и так же \(\bar\lambda_{100}\) и \(\bar\lambda_{50}\). Интегралы берутся по формуле трапеций на 801 точке; профиль снимается на 520 точках, сгущённых у стенок.

Коэффициент сопротивления сетки:

\[ C_x = \frac{F_x}{S_c\,q},\qquad q = \frac{\rho_о u_о^2}{2},\qquad u_о = a_{кр}\lambda_о , \]

где \(F_x\) — сила давления и трения, действующая на прутки вдоль \(x\) и отнесённая к половине канала, то есть к глубине 50 мм.

Эталон: RANS — осреднённые по Рейнольдсу уравнения Навье–Стокса

Мгновенные величины раскладываются на среднее и пульсацию: \(\phi = \bar\phi + \phi'\). Для сжимаемого газа скорость и температура осредняются по Фавру, \(\tilde u = \overline{\rho u}/\bar\rho\). Осреднённые уравнения:

\[ \begin{aligned} &\frac{\partial\bar\rho}{\partial t} + \nabla\!\cdot(\bar\rho\tilde{\mathbf u}) = 0,\\ &\frac{\partial(\bar\rho\tilde{\mathbf u})}{\partial t} + \nabla\!\cdot(\bar\rho\tilde{\mathbf u}\otimes\tilde{\mathbf u}) = -\nabla\bar p + \nabla\!\cdot\big[(\mu+\mu_t)\,2\mathbf S^{d}\big] - \tfrac23\nabla(\bar\rho k),\\ &\frac{\partial(\bar\rho\tilde E)}{\partial t} + \nabla\!\cdot\big[(\bar\rho\tilde E+\bar p)\tilde{\mathbf u}\big] = \nabla\!\cdot\Big[(\mu+\mu_t)\,2\mathbf S^{d}\!\cdot\tilde{\mathbf u} + c_p\Big(\frac{\mu}{\Pr}+\frac{\mu_t}{\Pr_t}\Big)\nabla\tilde T\Big], \end{aligned} \]

где \(\mathbf S^{d} = \tfrac12(\nabla\tilde{\mathbf u}+\nabla\tilde{\mathbf u}^{T}) - \tfrac13(\nabla\!\cdot\tilde{\mathbf u})\mathbf I\), \(\Pr = 0{,}72\), \(\Pr_t = 0{,}9\). Молекулярная вязкость задаётся законом Сазерленда:

\[ \mu(T) = \mu_{0}\Big(\frac{T}{T_{0}}\Big)^{3/2}\frac{T_{0}+S}{T+S}, \]

где \(\mu_0 = 1{,}716\cdot10^{-5}\) Па·с, \(T_0 = 273{,}15\) K, \(S = 110{,}4\) K.

Напряжения Рейнольдса \(-\overline{\rho u_i'u_j'}\) по гипотезе Буссинеска выражаются через турбулентную вязкость \(\mu_t\). Её даёт модель турбулентности.

Модель турбулентности k-ω SST (Ментер)

\[ \begin{aligned} &\frac{\partial(\rho k)}{\partial t} + \nabla\!\cdot(\rho\mathbf u k) = P_k - \beta^*\rho\,\omega k + \nabla\!\cdot\big[(\mu+\sigma_k\mu_t)\nabla k\big],\\ &\frac{\partial(\rho\omega)}{\partial t} + \nabla\!\cdot(\rho\mathbf u\,\omega) = \frac{\gamma_\omega\,\rho}{\mu_t}P_k - \beta\rho\,\omega^2 + \nabla\!\cdot\big[(\mu+\sigma_\omega\mu_t)\nabla\omega\big] + 2(1-F_1)\frac{\rho\,\sigma_{\omega2}}{\omega}\nabla k\cdot\nabla\omega,\\ &\mu_t = \frac{\rho\,a_1 k}{\max(a_1\omega,\;\Omega F_2)},\qquad a_1 = 0{,}31 . \end{aligned} \]

Функция смешения \(F_1\) переключает модель между k-ω у стенок и k-ε вдали от них. Ограничитель в \(\mu_t\) не даёт завышать турбулентную вязкость в областях торможения. На входе задана интенсивность турбулентности 1 % и \(\mu_t/\mu = 10\). У стенки \(k = 0\), \(\omega = 60\nu/(\beta_1 y_1^2)\).

RANS: граничные условия

ГраницаУсловие
Вход \(x=0\)Полные параметры \(P^*_о = 101325\) Па, \(T^* = 288\) K, направление \((1,0)\). Скорость на входе не задаётся, её определяет течение.
ВыходСтатическое давление \(p_{вых}\). Оно подбирается внешним итерационным циклом (метод секущих) так, чтобы расход совпал с таблицей 2: \[p_{вых}^{(k+1)} = p_{вых}^{(k)} + \frac{G_{цель} - G^{(k)}}{\partial G/\partial p}\bigg|_{\text{секущая}},\qquad |\Delta p| \le 0{,}25\,q .\] Критерий: \(|G/G_{цель}-1| < 0{,}3\,\%\), дисбаланс вход/выход < 0,5 %. В компрессорном SU2 нет простого условия «массовый расход на выходе», поэтому он задаётся так.
Стенки канала, пруткиПрилипание, тепловой поток \(q_w = 0\).

Удлинённый канал

Расчётная область продолжена до \(x = 725\) мм на грубой сетке. Если поставить выходную границу прямо в \(x = 525\) мм, через неё проходят следы от прутков, в них возникает обратный ток, и условие постоянного давления становится неустойчивым. Профили \(\sigma\) и \(\lambda_2\) снимаются внутри области, в сечении \(x = 525\) мм, интерполяцией по ячейкам. Расход при этом не меняется: течение стационарное, и уравнение неразрывности даёт одно и то же значение в любом сечении.

Второе семейство прутков: эквивалентный диаметр

В сетке два семейства прутков. Прутки первого семейства идут вдоль \(z\), и в плоском срезе \(x\)–\(y\) они видны как круги. Прутки второго семейства лежат в плоскости полотна поперёк первых, и в двумерный расчёт они не попадают. Если считать только первое семейство с реальным \(d = 0{,}9\) мм, живое сечение получается \((5{,}4-0{,}9)/5{,}4 = 0{,}833\) вместо 0,694, и потери сильно занижены.

Поэтому в плоском расчёте прутки берутся эквивалентного диаметра: при том же шаге 5,4 мм живое сечение равно живому сечению реальной сетки:

\[ \frac{t - d_э}{t} = \Big(\frac{t-d}{t}\Big)^2 \] \[ d_э = t\left[1 - \Big(1-\frac dt\Big)^2\right] = 5{,}4\left[1-\Big(\frac{4{,}5}{5{,}4}\Big)^2\right] \approx 1{,}65\ \text{мм}. \]

Проверка на режиме α = 80°, λо = 0,3

Модель\(\bar\sigma_{100}\)\(\bar\sigma_{50}\)\(C_x\)
RANS, \(d = 0{,}9\) мм (только первое семейство)0,9890,9900,67
RANS, \(d_э = 1{,}65\) мм (эквивалентное живое сечение)0,9700,9711,74
LBM в браузере (второе семейство — распределённое сопротивление)0,9740,978—

По словам преподавателя, значение LBM совпало с экспериментом. С эквивалентным диаметром RANS даёт такое же падение полного давления. Сайт показывает этот вариант.

RANS: сетка и численный метод

ЭлементРеализация
Сетка (gmsh)Неструктурированная гибридная. Вокруг каждого прутка и вдоль стенок — пограничный слой из четырёхугольников: первая ячейка 4 мкм, рост 1,2, 12 слоёв, \(y^+\lesssim 1\). На окружности прутка 48 узлов. Сгущение в полосе вдоль полотна и в следе до \(x \approx 535\) мм, треугольники вне пограничных слоёв. Дополнительное сгущение у стенок (0,4 мм) и в четырёх углах канала: там начинается пограничный слой, и без сгущения в нескольких ячейках возникает ложная сверхзвуковая область. Всего от 180 тыс. (80°) до 345 тыс. (20°) ячеек.
РешательSU2 8.3, сжимаемый RANS, неявная схема Эйлера по псевдовремени, линейный решатель FGMRES с ILU-предобусловливателем, MPI на 8 ядрах.
КонвекцияСхема Роэ. Сначала первый порядок (800 + 400 итераций) для начального поля, затем второй порядок: реконструкция MUSCL, ограничитель Venkatakrishnan–Wang. Для k и ω — противопоточная скалярная схема.
Число КурантаАдаптивное, от 5 до 40.
СходимостьСопротивление, расход и полные давления выходят на постоянные значения. Среднеквадратичная невязка плотности падает примерно до \(10^{-4}\) и дальше не меняется: остаются небольшие колебания в следах прутков. Величины осредняются по последним 600 итерациям.

Почему решение стационарное. При \(\mathrm{Re}_d \sim 10^3\ldots10^4\) за отдельными прутками физически идёт срыв вихрей. Но RANS описывает осреднённое течение: пульсации входят в турбулентную вязкость, поэтому следы получаются гладкими и одинаковыми у всех прутков. Нестационарный вариант (URANS) тоже проверялся. На этой постановке внутренние итерации шага по времени не сходились, баланс массы нарушался, а расчёт шёл заметно дольше. Поэтому эталоном выбрано стационарное решение второго порядка.

Запирание сетки при α = 80°, λо = 0,5

При заданных \(P^*_о\) и \(T^*\) расход через сечение пропорционален площади и функции расхода \(q(\lambda)\). Её максимум \(q = 1\) достигается при \(\lambda = 1\):

\[ G = m\,\frac{P^*}{\sqrt{T^*}}\,A\,q(\lambda),\qquad q(\lambda) = \lambda\Big(\frac{\gamma+1}{2}\Big)^{\frac{1}{\gamma-1}}\Big(1-\frac{\gamma-1}{\gamma+1}\lambda^2\Big)^{\frac{1}{\gamma-1}},\qquad m = \sqrt{\frac{\gamma}{R_г}\Big(\frac{2}{\gamma+1}\Big)^{\frac{\gamma+1}{\gamma-1}}}. \]

Просветы между прутками составляют долю \(f_α\) высоты канала. Без потерь полного давления до просветов неразрывность требует

\[ q(\lambda_{просв}) = \frac{q(\lambda_о)}{f_α},\qquad q(0{,}5) = 0{,}709,\qquad f_{80^\circ} \approx 0{,}705\;\Longrightarrow\; q(\lambda_{просв}) \approx 1{,}01 > 1 . \]

Табличный расход через сетку пройти не может, даже без учёта сужения струй за прутками. Струи в просветах становятся звуковыми, расход ограничивается, а избыточный перепад давления срабатывается в скачках за сеткой. Расчёт показывает именно это: проходит ≈ 95,5 % табличного расхода, \(\bar\sigma_{100} = 0{,}73\), \(\bar\lambda_{100} \approx 0{,}98\). На наклонных сетках полотно длиннее: \(f_{60^\circ} \approx 0{,}80\), и сетка ещё пропускает расход, хотя течение в просветах уже близко к звуковому (\(\bar\sigma_{100} = 0{,}86\)). Реальная сетка имеет то же живое сечение 0,694, поэтому запирание при α = 80° физическое, а не следствие эквивалентного диаметра. На сайте этот режим помечен звёздочкой.

Быстрая оценка в браузере: решёточный метод Больцмана

Интерактивный режим «LBM · браузер» решает уравнение Больцмана на решётке D2Q9:

\[ f_i(\mathbf x + \mathbf c_i\Delta t,\;t+\Delta t) = f_i - \Omega_i^{TRT}(f - f^{eq}), \] \[ f_i^{eq} = w_i\rho\Big[1 + \frac{\mathbf c_i\!\cdot\mathbf u}{c_s^2} + \frac{(\mathbf c_i\!\cdot\mathbf u)^2}{2c_s^4} - \frac{u^2}{2c_s^2}\Big],\quad c_s = \tfrac{1}{\sqrt3}. \]

Где применим LBM

При \(\lambda_о = 0{,}3\) (\(M \approx 0{,}28\)) влияние сжимаемости мало, и LBM согласуется с RANS. При \(\lambda_о \ge 0{,}4\) в просветах сетки \(M\) доходит до 0,5–0,6, коэффициент потерь уже зависит от числа Маха, и несжимаемая оценка расходится с правильным ответом. Кроме того, Re в браузерном расчёте ниже физического, иначе решётка неустойчива. Эталонный результат — RANS.

Ограничения и допущения — сводка

ДопущениеСледствие
Плоская задача (срез \(x\)–\(y\))Второе семейство прутков учтено эквивалентным диаметром, а не геометрически. Трёхмерные эффекты у стенок \(z\) не учитываются.
Модель турбулентности k-ω SSTСтандартный выбор для отрывных течений у стенок. Отрыв с цилиндра RANS описывает приближённо, и абсолютное значение \(C_x\) зависит от модели сильнее, чем \(\sigma\).
Стационарное решениеДаёт осреднённые профили. Мгновенная вихревая дорожка не разрешается.
Расчётная сетка coarseПограничный слой у прутков разрешён до стенки, \(y^+ \lesssim 1\). Проверка на сетке medium при α = 20°, λо = 0,3: \(\bar\sigma_{100}\) = 0,9533 против 0,9549 на coarse, \(C_x\) = 1,495 против 1,451. Результат от сетки практически не зависит.
Потери при α = 20° больше, чем при 80°Так получается и с реальными прутками \(d = 0{,}9\) мм: \(\bar\sigma_{100}\) = 0,973 при 20° против 0,989 при 80°. Значит, это не следствие эквивалентного диаметра. \(C_x\), отнесённый к площади полотна, при наклоне падает (1,74 → 1,45), но полотно при 20° почти вдвое больше (\(S_c\) = 0,00286 против 0,00154 м²). Полная сила и потеря полного давления \(\Delta P^* \approx C_x q S_c/S_о\) поэтому растут.
Выходное давление подбирается итерационноРасход совпадает с таблицей 2 с точностью ≈0,5 %.

Программы: браузерные решатели — solver-worker.js (Лаваль) и mesh-worker.js (LBM). RANS — rans/mesh.py (сетка gmsh), rans/case.cfg.tmpl (настройки SU2), rans/run_case.py (расчёт и подбор расхода), rans/post.py (профили, средние, \(C_x\)).