$$$\def \bm#1{\boldsymbol#1}$

Quasi-Biennial Oscillation

Nonhydrostatic Model

Численные симуляции будем проводить в рамках приближения Буссинеска для несжимаемого уравнения Навье-Стокса. Соответсвующие уравнения принимают вид:

где $B_{bg}(z)=-g{\rho }_{bg}(z)/{\rho }_0$$B_{bg}(z) = -g \rho_{bg}(z)/\rho_0$ – фоновое значение плавучести, связанное с равновесной стратификацией жидкости. Напомним, что в рамках приближения Буссинеска $\rho (\mathit{r},t)={\rho }_0+{\rho }_{bg}(z)+\delta \rho (\mathit{r},t)$$\rho(\bm r,t) = \rho_0 + \rho_{bg} (z) + \delta \rho(\bm r,t)$, причем ${\rho }_0\gg {\rho }_{bg}(z)+\delta \rho$$\rho_0 \gg \rho_{bg}(z) + \delta \rho$. Динамические вариации плотности $\delta \rho (\mathit{r},t)$$\delta \rho (\bm r, t)$ определяют плавучесть $b=-g\delta \rho /{\rho }_0$$b = -g \delta \rho/\rho_0$. В дальнейшем мы будем рассматривать линейно стратифицированную жидкость, и в этом случае $B_{bg}(z)=N^{2}z$$B_{bg}(z) = N^2 z$. Постоянная величина $N$$N$ имеет размерность частоты и называется Brünt-Vaisälä frequency. Окончательно получаем:

где для сокращения записи введено обозначение $D_t={\partial }_t+(\mathit{v}\cdot \mathrm{\nabla })$$D_t = \partial_t + (\bm v \cdot \nabla)$.

Wave Packet

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

где $a(x,z)$$a(x,z)$ – гауссова огибающая, а частота волны должна удовлетворять закону дисперсии:

Начальные условия в нашей задачи должны быть согласованы с полем давления при $t=0$$t=0$. Соответсвующие выражения имеют вид:

Результаты для параметров $N=1,k=8,m=16$$N=1, \; k=8, \; m=16$ показаны ниже:

QBO Setup

Теперь перейдем непосредственно к основному расчету. В направлении оси $X$$X$ граничные условия периодические для полей скорости и плавучести. На верхней границе мы используем stress-free граничные условия для поля скорости и отсутсвие нормального градиента для плавучести:

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

Численные параметры: $L_x=1.517388$$L_x = 1.517388$, $L_z=0.43$$L_z=0.43$, $N=1.56605$$N=1.56605$, $\nu =1.004\cdot {10}^{-6}$$\nu =1.004 \cdot 10^{-6}$, $\kappa =1.5\cdot {10}^{-7}$$\kappa = 1.5 \cdot 10^{-7}$, $ϵ=0.0016$$\epsilon=0.0016$, $\omega =0.43$$\omega = 0.43$, $k=16\pi /L_x$$k = 16 \pi/L_x$.

Расчет проведен до времени $T=50000$$T=50000$, шаг по времени $dt={10}^{-2}$$dt = 10^{-2}$. Сохранение данных проводится с шагом $\mathrm{\Delta }t=10$$\Delta t = 10$. Размер сетки $512\times 512$$512 \times 512$. В вертикальном направлении сетка неоднородная – более мелкая вблизи волнопродуктора, где нужно разрешать вязкий пограничный слой. Характерный размер пограничных слоев для полей скорости и плавучести можно оценить как $\sqrt{\nu /\omega }$$\sqrt{\nu/\omega}$ и $\sqrt{\kappa /\omega }$$\sqrt{\kappa/\omega}$ соответсвенно.

Linear Internal Waves

Для проверки численной схемы было предложено смоделировать распространение инерционной волны в вертикальном направлении. Размер области $L_z$$L_z$ в направлении оси $Z$$Z$ должен быть достаточно большим, чтобы можно было пренебречь отраженной волной. Кроме того, аналитический расчет предполагает выполнение условия $\nu \gg \kappa$$\nu \gg \kappa$. Граничные условия на волнопродукторе аналогичны предыдущей задаче. За пределами пограничных слоев $z\gg \sqrt{\nu /\omega },\sqrt{\kappa /\omega }$$z \gg \sqrt{\nu/\omega}, \; \sqrt{\kappa/\omega}$ для вертикальной компоненты скорости можно найти:

где $m=k\sqrt{N^2/{\omega }^2-1}$$m = k \sqrt{N^2/\omega^2 - 1}$. Численные параметры выберем в соответсвии с расчетом QBO, только уменьшим амплитуду накачки $ϵ={10}^{-5}$$\epsilon = 10^{-5}$ (для выполнения линейного приближения), увеличим вязкость среды $\nu =1.5\cdot {10}^{-6}$$\nu = 1.5 \cdot 10^{-6}$ и размер $L_z=L_x$$L_z=L_x$. Значения остальных параметров: $L_x=1.517388$$L_x = 1.517388$, $N=1.56605$$N=1.56605$, $\kappa =1.5\cdot {10}^{-7}$$\kappa = 1.5 \cdot 10^{-7}$, $\omega =0.43$$\omega = 0.43$, $k=16\pi /L_x$$k = 16 \pi/L_x$.

В этом случае длина распространения волны $l_d\sim 0.31$$l_d \sim 0.31$ и стационарный режим установится за характерное время $\tau \sim l_d/v_g^z\sim 100$$\tau \sim l_d/v_g^z \sim 100$, где групповая скорость $v_g^z={\displaystyle \frac{\omega m}{k^2+m^2}}\sim 3.4\cdot {10}^{-3}$$v_g^z = \dfrac{\omega m }{k^2+m^2}\sim 3.4 \cdot 10^{-3}$.

Расчет проведем на адаптивной сетке $2048\times 2048$$2048 \times 2048$ со сгущением к волнопродуктору. Длительность расчета $50T$$50 T$, где период колебания $T=2\pi /\omega$$T = 2 \pi/\omega$. Данные сохраняем с шагом $\mathrm{\Delta }t=0.1T$$\Delta t = 0.1 T$, шаг интегрирования $dt=5\cdot {10}^{-3}$$dt = 5 \cdot 10^{-3}$.

Для сравнения расчетов с теорией возьмем $t_0=50T$$t_0= 50 T$ и наложим на один график $w(x,z,t_0)$$w(x,z,t_0)$ и теоретическую огибающую. Результаты представлены ниже для значения $\sin (kx)\approx -0.999$$\sin(kx) \approx -0.999$: