Задача о взаимодействии ударной волны с вихрем (SVI problem)

Базовая постановка (Basic formulation)

Задача SVI является многопараметрической задачей. Обычно в качестве объекта моделирования рассматривают совершенный газ с показателем адиабаты γ = 1.4, систему координат связывают с ударной волной (в начальный момент времени t = 0 фронт ударной волны неподвижен и расположен в сечении x = x_S = const), а двумерный вихрь полагают изоэнтропическим.

Основными параметрами задачи SVI являются: M_S = u_x / c – ударно-волновое число Маха, характеризующее интенсивность ударной волны (u_x – скорость потока перед волной, c – скорость звука), V = V(r) – профиль скорости цилиндрического вихря (r – расстояние от центра вихря), M_V = V_m / c – вихревое число Маха, характеризующее интенсивность вихря (V_m – максимальная скорость на профиле). Другими параметрами задачи являются: r₀ – эффективный радиус вихря (определяется точкой на профиле скорости, где V(r₀) = V_m), (x_V, y_V) – координата центра вихря в начальный момент времени.

При решении задачи SVI в качестве расчетной области в системе координат xy обычно используют прямоугольник, поперечный размер которого равен единице. На верхней и нижней границах области задается условие симметрии течения, на левой границе задается сверхзвуковой втекающий поток, а на правой границе – дозвуковой вытекающий поток при заданном давлении (номинальное давление за ударной волной). Расчет задачи проводится до заданного момента времени t = t₁, обеспечивающего полное прохождение вихря через ударную волну и его последующее значимое продвижение в продольном направлении.

Чаще всего задача SVI моделируется в рамках уравнений Эйлера (невязкое приближение) [1–3, 5], однако иногда для ее решения используются уравнения Навье–Стокса (вязкое приближение) [4, 6, 7]. В последнем случае требуется задание числа Рейнольдса Re (далее число Рейнольдса будем определять по параметрам потока перед ударной волной и эффективному диаметру вихря, равному 2r₀) и, строго говоря, числа Прандтля Pr.

Наиболее интересным с вычислительной точки зрения является тестовый случай, когда интенсивный вихрь проходит через сильную ударную волну. В этом случае картина течения имеет сложную структуру с множественными скачками уплотнения и поверхностями разрывов. Однако, если данную задачу решать в невязком приближении на достаточно подробной сетке, то контактные поверхности станут неустойчивыми, и поэтому сеточная сходимость будет отсутствовать. Если задачу решать в вязком приближении, то можно найти значение числа Рейнольдса, при котором сложная ударно-волновая структура течения сохранится, а численное решение при увеличении сеточного разрешения будет оставаться устойчивым.

Для тестирования численных методов предлагается вариант задачи SVI с параметрами γ = 1.4, M_S = 3, M_V = 0.8, Re = 10⁴, Pr = 0.75, x_S = 0, (x_V, y_V) = (–0.5, 0.5), r₀ = 0.075 и профилем скорости вихря, описываемым формулой:

V(r) = V_m ⋅ (r / r₀) ⋅ exp{[1 – (r / r₀)²] / 2},  где r = [(x–x_V)² + (y–y_V)²]^0.5.

В базовой постановке в качестве расчетной области будем использовать прямоугольник [–1, 1] × [0, 1].

Рассмотрим детально процедуру задания полей начальных данных. При описании параметров газа будем применять общепринятые обозначения: u_x и u_y – компоненты вектора скорости газа, ρ – плотность, p – давление. Нижние индексы 1 и 2 припишем, соответственно, параметрам газа перед и за фронтом ударной волны (здесь пока речь идет только о равномерном фоновом потоке перед ударной волной, без учета «вихревого» возмущения). Положив ρ₁ = p₁ = 1, для остальных параметров будем иметь 

u_x1 = M_S (γp₁ / ρ₁)^0.5 = 3⋅(1.4)^0.5,   u_y1 = 0,

u_x2 = u_x1 ⋅ [(γ–1)M_S² + 2] / [(γ+1)M_S²] = 7⋅(1.4)^0.5/9,   u_y2 = 0,

ρ₂ = ρ₁ ⋅ [(γ+1)M_S²] / [(γ–1)M_S² + 2] = 27/7,   p₂ = p₁ ⋅ [2γM_S² – (γ–1)] / (γ+1) = 31/3.

Теперь скорректируем поле начальных данных, добавив цилиндрический вихрь в фоновый поток перед ударной волной; вихрь будем полагать вращающимся по часовой стрелке. Алгоритмически, процедура коррекции состоит из следующих последовательных шагов. Сначала вычисляется максимальная скорость на профиле вихря V_m = M_V (γp₁ / ρ₁)^0.5 = 0.8⋅(1.4)^0.5. Затем в каждой расчетной точке области x < x_S = 0 определяются два локальных параметра:

f(r) = V_m ⋅ exp{[1 – (r / r₀)²] / 2},   g(r) = 1 – f ² ⋅ [(γ–1)ρ₁] / (2γp₁) = 1 – f ²/7.

Окончательно, параметры газа перед ударной волной вычисляются как

u_x = u_x1 + f ⋅ (y–y_V) /  r₀,   u_y = u_y1 – f ⋅ (x–x_V) /  r₀,

ρ = ρ₁ ⋅ g^1/(γ–1) = g^2.5,    p = p₁ ⋅ g^γ/(γ–1) = g^3.5.

Коэффициент динамической вязкости полагается постоянным и равным μ = ρ₁u_x1(2r₀) / Re = 0.45⋅(1.4)^0.5⋅10^–4.

Задача рассчитывается до момента времени t = t₁ = 1.5 / u_x1 = 0.5⋅(1.4)^–0.5.

Замечание относительно начальных (стартовых) ошибок

Такого рода ошибки в англоязычной литературе называют «start-up errors» (см., например, раздел 15.8.4 в книге [10] и раздел 3.3 в статье [9]). Они возникают из-за неподходящего (несогласованного с численной вязкостью) размывания ударной волны в начальных данных. Хотя начальные ошибки носят локальный характер, но могут заметно искажать и без того сложную структуру течения в данной задаче. Если используемый расчетный метод подвержен такого рода ошибкам, то их можно легко устранить. Для этого в процессе расчета следует сделать однократную коррекцию: в момент времени t = t₁/6 (вихрь еще не начал взаимодействовать с ударной волной) восстановить начальные параметры в ячейках, расположенных в области x > x_S + 0.03 = 0.03.