Geodynamo simulations software package suite based on spectral hereditary models

Full Text

Введение

Задача геодинамо представляет собой задачу о магнитогидродинамической конвекции проводящей вязкой несжимаемой жидкости во вращающейся сферической оболочке с твёрдыми границами (внешнее ядро Земли) [1-3]. Система уравнений геодинамо связывает поля скорости, давления, температуры и магнитной индукции и включает уравнения Навье-Стокса с кориолисовым и лоренцевым членами, уравнение температуропроводности и уравнение индукции.

Для её исследования применяют различные численные методы, в том числе и метод Галеркина. Следует отметить, что, хотя метод Галеркина и является по своему происхождению численным методом решения уравнений, его можно рассматривать при малом числе мод как способ построения малоразмерных моделей [4, 5]. Поля задачи при однородных граничных условиях представляются в виде линейных комбинаций стационарных базисных полей (мод) с зависящими от времени амплитудами. Каждое базисное поле должно удовлетворять граничным условиям, что гарантирует точное выполнение этих условий при любых амплитудах. Геометрически метод Галеркина выполняет тогда проектирование исходных уравнений на подпространства, порождаемые модами разложения.

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

В случае маломодовых моделей, содержащих только крупномасштабные моды, необходимо включать в уравнения динамо в явном виде альфа-эффект. Он описывает генерацию крупномасштабного поля турбулентными мелкомасштабными вихрями. Чтобы избежать неограниченного роста, необходима и обратная связь в виде подавления альфа-эффекта энергией этого поля [8, 9]. Известно, что в турбулентных магнитогидродинамических системах проявляется свойство памяти [10], поэтому подавление α-эффекта в рассматриваемых моделях вводится с учетом этого свойства, а именно, подавление определяется не только актуальным, но и всеми предшествующими значениями энергии. Математически это формализуется в виде функционала от энергии с разностным ядром, и уравнения модели являются интегро-дифференциальными. Для некоторых частных типов ядер возможно сведение уравнений только к дифференциальным. Крупномасштабные приближения неявно предполагают, что влияние мелкомасштабных мод в среднем нулевое. Однако эти моды могут спонтанно синхронизироваться на случайное время, оказывая влияние на генерацию. В модели это введено как аддитивная поправка в интенсивность альфа-эффекта в виде случайного процесса, представляющего собой последовательность прямоугольных импульсов, возникающих в случайные моменты, имеющих случайную длительность и гауссовскую, в среднем нулевую, амплитуду.

Для автоматизированного составления описанных моделей и проведения численного моделирования разработан программный комплекс комбинированных символьно-численных вычислений, описываемый в настоящей работе.

Уравнения эредитарной спектральной модели

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

Жидкое ядро Земли рассматривается как сферическая оболочка вязкой проводящей жидкости с твердыми границами, равномерно вращающейся вокруг оси Oz. В сферической системе координат (r, θ,φ) внутренняя граница ядра (ICB) определяется как $r\mathrm{<mo>=</mo>}{r}_i$, а граница ядро $–$ мантия (CMB) как $r\mathrm{<mo>=</mo>}{r}_o$. Температура на ICB и CMB постоянна и равна $T_i$ и $T_o$ соответственно.

В качестве единицы длины принимается внешний радиус Земли, а в качестве единицы температуры принимается $\delta T\mathrm{<mo>=</mo>}{T}_i-T_o$. Тогда в этих единицах $r_o\mathrm{<mo>=</mo>1}$, $r_i\mathrm{<mo>=</mo>0.35}$ и $T_o\mathrm{<mo>=</mo>}{T}_i+1$.

Поле $T\mathrm{<mo>=</mo>}T\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$ задает отклонение температуры относительно равновесного профиля (без конвекции)

$\begin{array}{ll}{T}_s\mathrm{<mo>=</mo>}\frac{r_i}{1-r_i}\left(\frac{1}{r}-1\right)+T_i-1.\hfill & \hfill \end{array}$

Обезразмеренные уравнения геодинамо в ядре ( $r\le 1$ ) имеют вид:

$\begin{array}{ll}\frac{\partial v}{\partial t}+\left(v\nabla \right)v\mathrm{<mo>=</mo>}-\nabla p+\Delta v-2{\text{E}}^{-1}{e}_z\times v+{\text{RaPr}}^{-1}Tr+\text{rot}B\times B\mathrm{,}\hfill & \hfill \\ \frac{\partial T}{\partial t}+\left(v\cdot \nabla \right)\left(T+T_s\right)\mathrm{<mo>=</mo>}{\text{Pr}}^{-1}\Delta T\mathrm{,}\hfill & \hfill \\ \frac{\partial B}{\partial t}\mathrm{<mo>=</mo>}\text{rot}\left(v\times B\right)+{\text{R}}_{\alpha }\text{rot}\left(\alpha \mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}B\right)+{\text{Pm}}^{-1}\Delta B\mathrm{,}\hfill & \hfill \\ \nabla v\mathrm{<mo>=</mo>0,}\nabla B\mathrm{<mo>=</mo>0,}\hfill & \hfill \end{array}$ (1)

где $\alpha \mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}a\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}\cos \theta$ $–$ заданное осесимметричное поле тензора параметризованного $\alpha$ -эффекта.

Параметры подобия: $\text{E}$ $–$ число Экмана, $\text{Ra}$ $–$ число Релея, $\text{Pr}$ $–$ число Прандтля, $\text{Pm}$ $–$ магнитное число Прандтля, ${\text{R}}_{\alpha }$ $–$ амплитуда $\alpha$ -эффекта.

Среда вне ядра предполагается непроводящей (токи отсутствуют), значит, при $r\mathrm{<mo>></mo>1}$:

$\begin{array}{ll}B\mathrm{<mo>=</mo>}-\nabla \Psi \mathrm{,}\text{где}\Delta \Psi \mathrm{<mo>=</mo>0.}\hfill & \hfill \end{array}$ (2)

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

$\begin{array}{ll}v\left(r\mathrm{<mo>=</mo>}{r}_i\right)\mathrm{<mo>=</mo>}v\left(r\mathrm{<mo>=</mo>1}\right)\mathrm{<mo>=</mo>0,}T\left(r\mathrm{<mo>=</mo>}{r}_i\right)\mathrm{<mo>=</mo>}T\left(r\mathrm{<mo>=</mo>1}\right)\mathrm{<mo>=</mo>0,}\hfill & \hfill \\ B\left(r\mathrm{<mo>=</mo>0}\right)\mathrm{<mo><</mo>}\infty \mathrm{,}B\mathrm{<mo>=</mo>}-\nabla \Psi \text{при}r\mathrm{<mo>=</mo>1,}\Psi \mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo>=</mo>}+\infty \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>0.}\hfill & \hfill \end{array}$ (3)

Общее решение уравнения (2) в виде разложения по сферическим функциям прекрасно известно [11], поэтому решать надо только задачу (1), а выражение для $\Psi$ дает граничные условия для аналогичного (по сферическим функциям) разложения поля внутри ядра [7].

Поля скорости внутри ядра аппроксимируются конечными линейными комбинациями стационарных полей (мод):

$\begin{array}{l}v\left(r\mathrm{,}t\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{l\mathrm{<mo>=</mo>1}}^{L_{\max }}}{\beta }_l\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}{v}_l\left(r\right)\mathrm{,}T\left(r\mathrm{,}t\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{s\mathrm{<mo>=</mo>1}}^{S_{\max }}}{\alpha }_s\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}{T}_s\left(r\right)\mathrm{,}B\left(r\mathrm{,}t\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{p\mathrm{<mo>=</mo>1}}^{P_{\max }}}{\gamma }_p\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}{B}_p\left(r\right)\mathrm{.}\hfill \end{array}$ (4)

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

$\begin{array}{l}\mu v-\nabla p+\Delta v\mathrm{<mo>=</mo>0,}\nabla v\mathrm{<mo>=</mo>0,}\hfill \\ \lambda T+\Delta T\mathrm{<mo>=</mo>0,}\hfill \\ \eta B+\Delta B\mathrm{<mo>=</mo>0,}\nabla B\mathrm{<mo>=</mo>0}\hfill \end{array}$ (5)

с соответствующими граничными условиями. Расчет собственных значений и собственных мод этих задач подробно описан в [7].

Операторы этих задач самосопряжены относительно скалярных произведений:

$\begin{array}{l}{⟨v_1\mathrm{,}{v}_2⟩}_v\mathrm{<mo>=</mo>}{\displaystyle \underset{r_i\le r\le 1}{\iiint }}{v}_1{v}_{2}dV\mathrm{<mo>=</mo>}{\displaystyle {\int }_{r_i}^1}{r}^{2}dr{\displaystyle {\int }_0^{\pi }}\sin \theta d\theta {\displaystyle {\int }_{-\pi }^{\pi }}{v}_1{v}_{2}d\phi \mathrm{,}\hfill \\ {⟨T_1\mathrm{,}{T}_2⟩}_T\mathrm{<mo>=</mo>}{\displaystyle \underset{r_i\le r\le 1}{\iiint }}{T}_1{T}_{2}dV\mathrm{<mo>=</mo>}{\displaystyle {\int }_{r_i}^1}{r}^{2}dr{\displaystyle {\int }_0^{\pi }}\sin \theta d\theta {\displaystyle {\int }_{-\pi }^{\pi }}{T}_1{T}_{2}d\phi \mathrm{,}\hfill \\ {⟨B_1\mathrm{,}{B}_2⟩}_B\mathrm{<mo>=</mo>}{\displaystyle \underset{0\le r\le 1}{\iiint }}{B}_1{B}_{2}dV\mathrm{<mo>=</mo>}{\displaystyle {\int }_0^1}{r}^{2}dr{\displaystyle {\int }_0^{\pi }}\sin \theta d\theta {\displaystyle {\int }_{-\pi }^{\pi }}{B}_1{B}_{2}d\phi \mathrm{,}\hfill \end{array}$ (6)

Системы собственных мод полны и ортогональны. Их можно считать нормированными.

Каждая мода скорости и магнитного поля определяется мультииндексами вида $\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{,}n\mathrm{,}m\mathrm{,}type\mathrm{<mo\; stretchy="false">)</mo>}$, где $k$, $n$ и $m$ соответствуют дискретизации спектра по $r$, $\theta$ и $\phi$ соответственно, а $type$ $–$ бинарный индекс типа моды: тороидальная или полоидальная. Температурные моды определяются мультииндексами $\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{,}n\mathrm{,}m\mathrm{<mo\; stretchy="false">)</mo>}$.

Подстановка разложений (4) в уравнения геодинамо и применение процедуры Галeркина дает динамическую систему для амплитуд:

$\begin{array}{ll}\frac{d{\beta }_l}{dt}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{,}j\mathrm{<mo>=</mo>1}}^{L_{\max }}}{B}_{lij}{\beta }_i{\beta }_j-{\mu }_l{\beta }_l+{\text{E}}^{-1}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^{L_{\max }}}{E}_{li}{\beta }_i+{\text{RaPr}}^{-1}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^{S_{\max }}}{C}_{li}{\alpha }_i+{\displaystyle \sum _{i\mathrm{,}j\mathrm{<mo>=</mo>1}}^{P_{\max }}}{Q}_{lij}{\gamma }_i{\gamma }_j\mathrm{,}\hfill & \hfill \\ l\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{L}_{\max }\mathrm{,}\hfill & \hfill \\ \frac{d{\alpha }_s}{dt}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{,}j\mathrm{<mo>=</mo>1}}^{L_{\max }\mathrm{,}{S}_{\max }}}{F}_{sij}{\beta }_i{\alpha }_j+{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^{L_{\max }}}{H}_{si}{\beta }_i-{\text{Pr}}^{-1}{\lambda }_s{\alpha }_s\mathrm{,}s\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{S}_{\max }\mathrm{,}\hfill & \hfill \\ \frac{d{\gamma }_p}{dt}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{,}j\mathrm{<mo>=</mo>1}}^{L_{\max }\mathrm{,}{P}_{\max }}}{W}_{pij}{\beta }_i{\gamma }_j+{\text{R}}_{\alpha }{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^{P_{\max }}}{W}_{pij}^{\alpha }{\gamma }_i-{\text{Pm}}^{-1}{\eta }_p{\gamma }_p\mathrm{,}p\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{P}_{\max }\mathrm{.}\hfill & \hfill \end{array}$ (7)

В этой системе ${\mu }_l\mathrm{<mo>></mo>0}$, ${\lambda }_s\mathrm{<mo>></mo>0}$ и ${\eta }_p\mathrm{<mo>></mo>0}$ $–$ собственные значения, а прописные буквы $–$ постоянные коэффициенты Галeркина, определяемые скалярными произведениями:

$\begin{array}{ll}{B}_{lij}\mathrm{<mo>=</mo>}-{⟨\left(v_i\nabla \right)v_j\mathrm{,}{v}_l⟩}_v\mathrm{,}{E}_{li}\mathrm{<mo>=</mo>}-2{⟨e_z\times v_i\mathrm{,}{v}_l⟩}_v\mathrm{,}{C}_{li}\mathrm{<mo>=</mo>}{⟨T_{i}r\mathrm{,}{v}_l⟩}_v\mathrm{,}\hfill & \hfill \\ Q_{lij}\mathrm{<mo>=</mo>}{⟨\text{rot}{B}_i\times B_j\mathrm{,}{v}_l⟩}_v\mathrm{,}{F}_{sij}\mathrm{<mo>=</mo>}-{⟨v_i\left(\nabla T_j\right)\mathrm{,}{T}_s⟩}_T\mathrm{,}{H}_{si}\mathrm{<mo>=</mo>}\frac{r_i}{1-r_i}{⟨r^{-2}{v}_i{e}_r\mathrm{,}{T}_s⟩}_T\mathrm{,}\hfill & \hfill \\ W_{pij}\mathrm{<mo>=</mo>}{⟨\text{rot}\left(v_i\times B_j\right)\mathrm{,}{B}_s⟩}_v\mathrm{,}{W}_{pi}^{\alpha }\mathrm{<mo>=</mo>}{⟨\text{rot}\left(\alpha \mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}{B}_i\right)\mathrm{,}{B}_s⟩}_v\mathrm{.}\hfill & \hfill \end{array}$ (8)

Спектральная модель геодинамо $–$ это набор мод и система уравнений (7). Если число мод невелико, говорят о маломодовой модели.

В модель вводится эредитарное подавление $\alpha$ -эффекта энергией поля и случайное возмущение $\alpha$ -эффекта $\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ в виде:

$\begin{array}{ll}{\text{R}}_{\alpha }\mapsto \frac{{\text{R}}_{\alpha }\left(1+\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}\right)}{1+u\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}u\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{1}{T_K}{\displaystyle {\int }_0^t}K\left(\frac{t-\tau }{T_K}\right){\displaystyle \sum _{p\mathrm{<mo>=</mo>1}}^{P_{max}}}{\gamma }_p^2\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>}d\tau \mathrm{,}\hfill & \hfill \end{array}$ (9)

где безразмерная $K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ $–$ ядро функционала подавления со свойствами $K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}\ge 0$, $K\mathrm{<mo\; stretchy="false">(</mo>}+\infty \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>0}$ и ${{\displaystyle \max }}_{t\ge 0}K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1}$. Параметр $T_K$ задает временной масштаб ядра.

Для численного решения система уравнений (7)-(9) дополняется начальными значениями для амплитуд мод:

$\begin{array}{ll}{\beta }_l\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{\beta }_l^0\mathrm{,}{\alpha }_s\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{\alpha }_s^0\mathrm{,}{\gamma }_p\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{\gamma }_p^0\mathrm{,}\hfill & \hfill \\ l\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{L}_{max}\mathrm{,}s\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{S}_{max}\mathrm{,}p\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{P}_{max}\mathrm{.}\hfill & \hfill \end{array}$ (10)

Стохастический процесс $\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ описывает влияние когерентных структур, спонтанно образованных мелкомасштабными модами, не учитываемыми явно в модели. В зависимости от морфологии каждой структуры она может как усиливать, так и ослаблять генерацию поля, а среднее значение этого влияния должно быть равно нулю.

Опишем модельную структуру этого процесса. Предположим, что $k$ -ая когерентная структура ( $k\mathrm{<mo>=</mo>1,2,}\dots$ ) спонтанно образуется в случайный момент времени ${\phi }_k$ и саморазрушается в случайный момент времени ${\theta }_k$. Возможностью одновременного существования двух или более структур пренебрегаем. Пусть также ${\theta }_0\mathrm{<mo>=</mo>0}$.

Тогда случайное время ожидания структуры ${\tau }_k^W\mathrm{<mo>=</mo>}{\phi }_k-{\theta }_{k-1}$, а случайное время ее существования ${\tau }_k^E\mathrm{<mo>=</mo>}{\theta }_k-{\phi }_k$. Интенсивность влияния $k$ -ой структуры на генерацию поля описывается величинами ${\xi }_k$. Все эти три класса величины предполагаются независимыми между собой при различных $k$ и одинаково распределенными в пределах одного класса.

Величины ${\xi }_k$ гауссовские $\mathcal{N}\mathrm{<mo\; stretchy="false">(</mo>0,}{\sigma }^2\mathrm{<mo\; stretchy="false">)</mo>}$, времена ожидания ${\tau }_k^W$ распределены по степенному закону с плотностью

$\begin{array}{ll}{p}_W\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{\nu -1}{c{T}^W}{\left(1+\frac{t}{c{T}^W}\right)}^{-\nu }\mathrm{,}{T}^W\mathrm{<mo>></mo>0,}c\mathrm{<mo>=</mo>}{\left(2^{\frac{1}{\nu -1}}-1\right)}^{-1}\mathrm{,}\nu \mathrm{<mo>></mo>1,}t\ge \mathrm{0,}\hfill & \hfill \end{array}$ (11)

а времена существования ${\tau }_k^E$ распределены по показательному закону с плотностью

$\begin{array}{ll}{p}_E\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{\ln 2}{T^E}\exp \left(-\frac{\ln 2}{T^E}t\right)\mathrm{,}{T}^E\mathrm{<mo>></mo>0,}t\ge 0.\hfill & \hfill \end{array}$ (12)

Параметры $T^W$ и $T^E$ имеют смысл медиан соответствующих случайных величин, принимаемых в качестве характерных времен ожидания и существования. Использование медиан вместо традиционных математических ожиданий в качестве характерных значений связано с тем, что у степенного распределения (11) при $\nu \le 2$ математическое ожидание бесконечно.

Случайных процесс $\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ определяется следующим образом:

$\begin{array}{ll}\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{k\mathrm{<mo>=</mo>1}}^{+\infty }}{\xi }_k\left[H\mathrm{<mo\; stretchy="false">(</mo>}t-{\tau }_k\mathrm{<mo\; stretchy="false">)</mo>}-H\mathrm{<mo\; stretchy="false">(</mo>}t-{\theta }_k\mathrm{<mo\; stretchy="false">)</mo>}\right]\mathrm{,}\hfill & \hfill \end{array}$ (13)

где $H\mathrm{<mo\; stretchy="false">(</mo>}\cdot \mathrm{<mo\; stretchy="false">)</mo>}$ $–$ функция Хевисайда. Он представляет собой последовательность прямоугольных импульсов с амплитудами ${\xi }_k$, возникающих в моменты ${\phi }_k$ и исчезающих в моменты ${\theta }_k$. Между импульсами процесс нулевой.

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

Структура комплекса

Представляемый вычислительный комплекс моделирования геодинамо предназначен для составления уравнений (7)$–$(9), т.е. расчета коэффициентов Галeркина и собственных значений, генерации реализаций шума $\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ и численного решения уравнений с заданными начальными значениями амплитуд.

Структура комплекса представлена на рис. 1.

 

Рис. 1. Структура комплекса.

Figure 1. Structure of the complex.

 

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

В пополняемых файлах параметров мод хранятся собственные значения и числовые параметры мод с идентификаторами в виде мультииндексов.

Исходными данными на входе являются параметры подобия, количества мод $L_{max}$, $S_{max}x$, $P_{max}$, тип и параметры ядра $K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$, параметры процесса $\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$, время моделирования $T_{max}$ и шаг временной сетки $h$.

Общая схема работы следующая:

1. По заданным мультииндексам идет обращение в файлы коэффициентов Галеркина. Если все необходимые коэффциенты в файлах есть, они передаются в модули численного решения. Если каких-то коэффициентов нет, соответствующие мультииндексы передаются на вход модуля расчета коэффициентов.
2. Модуль расчета коэффициентов обращается в файлы параметров мод по заданным мультииндексам. Если параметры всех необходимых мод в файлах есть, производится расчет коэффициентов, они дописываются в файл и передаются в модули численного решения. Если параметров каких-либо мод нет, соответствующие мультииндексы передаются в модуль расчета параметров мод.
3. Модуль расчета параметров проводит вычисления для заданных мультииндексов, пополняет файлы параметров мод и передает параметры в модуль расчета коэффициентов Галеркина.
4. По значениям параметров на входе комплекса и коэффициентам Галеркина как по входным данным работает один из двух модулей численного решения. Выбор одного из двух типов определяется типом ядра подавления. Подробнее это описано ниже в работе.
5. В процессе работы на каждом шаге по времени модуль численного решения обращается к модулю генерации шума за реализацией очередного шумового отсчета.
6. Результатом работы является файл временных рядов решений уравнений модели. Структура этого файла описана ниже.

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

Модули расчет параметров базисных мод

Собственные тороидальные моды скорости и магнитной индукции $v_{knm}^T$ и $B_{knm}^T$, собственные полоидальные моды скорости и магнитной индукции $v_{knm}^P$ и $B_{knm}^P$, собственные моды температуры $T_{knm}$ определяются выражениями

$\begin{array}{ll}{v}_{knm}^T\mathrm{<mo>=</mo>}\text{rot}\left(R_{kn}^T\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}{Y}_n^m\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{,}\phi \mathrm{<mo\; stretchy="false">)</mo>}r\right)\mathrm{,}{B}_{knm}^T\mathrm{<mo>=</mo>}\text{rot}\left(X_{kn}^T\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}{Y}_n^m\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{,}\phi \mathrm{<mo\; stretchy="false">)</mo>}r\right)\mathrm{,}\hfill & \hfill \\ v_{knm}^P\mathrm{<mo>=</mo>}\text{rot}\text{rot}\left(R_{kn}^P\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}{Y}_n^m\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{,}\phi \mathrm{<mo\; stretchy="false">)</mo>}r\right)\mathrm{,}{B}_{knm}^P\mathrm{<mo>=</mo>}\text{rot}\text{rot}\left(X_{kn}^P\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}{Y}_n^m\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{,}\phi \mathrm{<mo\; stretchy="false">)</mo>}r\right)\mathrm{,}\hfill & \hfill \\ T_{knm}\mathrm{<mo>=</mo>}{Z}_{kn}\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}{Y}_n^m\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{,}\phi \mathrm{<mo\; stretchy="false">)</mo>,}\hfill & \hfill \end{array}$ (14)

где $Y_m^m\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{,}\phi \mathrm{<mo\; stretchy="false">)</mo>}$ $–$ сферические гармоники, а функции радиальной переменной $r$ задаются в виде

$\begin{array}{ll}{R}_{kn}^T\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{A}_{kn}^T{j}_n\left(\sqrt{{\mu }_{kn}^T}r\right)+B_{kn}^T{y}_n\left(\sqrt{{\mu }_{kn}^T}r\right)\mathrm{,}\hfill & \hfill \\ X_{kn}^T\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{a}_{kn}^T{j}_n\left(\sqrt{{\eta }_{kn}^T}r\right)+b_{kn}^T{y}_n\left(\sqrt{{\eta }_{kn}^T}r\right)\mathrm{,}\hfill & \hfill \\ R_{kn}^P\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{C}_{kn}^1{j}_n\left(\sqrt{{\mu }_{kn}^P}r\right)+C_{kn}^2{y}_n\left(\sqrt{{\mu }_{kn}^P}r\right)+C_{kn}^3{r}^n+C_{kn}^4{r}^{-n-1}\mathrm{,}\hfill & \hfill \\ X_{kn}^P\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{a}_{kn}^P{j}_n\left(\sqrt{{\eta }_{kn}^P}r\right)+b_{kn}^P{y}_n\left(\sqrt{{\eta }_{kn}^P}r\right)\mathrm{,}\hfill & \hfill \\ Z_{kn}\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{a}_{kn}{j}_n\left(\sqrt{{\lambda }_{kn}}r\right)+b_{kn}{y}_n\left(\sqrt{{\lambda }_{kn}}r\right)\mathrm{.}\hfill & \hfill \end{array}$ (15)

В формулах (14) и (15) индекс $n\mathrm{<mo>=</mo>1,2,}\dots$ у мод скорости и магнитного поля, и $n\mathrm{<mo>=</mo>0,1,2,}\dots$ у мод температуры. Индексы $k\mathrm{<mo>=</mo>0,1,2,}\dots$ и $m\mathrm{<mo>=</mo>}-n\mathrm{,}\dots \mathrm{,}n$. Положительные ${\mu }_{kn}^T$, ${\mu }_{kn}^P$, ${\eta }_{kn}^T$, ${\eta }_{kn}^P$ и ${\lambda }_{kn}$ являются собственными значениями соответствующих мод, а прочие коэффициенты в (15) определяются из соответствующих граничных и нормировочных условий. Эти условия ставятся на функции (15) в результате применения условий (3) к модами типа (14). Уравнения на собственные значения и нормировочные условия приведены в работе [7].

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

Уравнения на собственные значения очень громоздкие, как и нормировочные условия, имеющие интегральную форму. Кроме того, нормировочные интегралы для магнитных мод являются несобственными второго рода. Поэтому модуль разработан в пакете Maple (лицензия 910346) и сочетает в себе символьные и численные вычисления. Технология расчета параметров мод с помощью Maple очень подробно описана в работе [7], поэтому здесь ограничимся общей схемой работы модуля.

На входе модуля: массивы мультииндексов мод скорости, температуры и магнитной индукции, для которых необходимо рассчитать собственные значения и коэффициенты из (15).

Схема работы модуля:

1. Программно формируется выражение радиальной функции моды.
2. Аналитически вычисляется норма с неопределенными собственными значениями и коэффициентами.
3. Программно формируется выражение левой части уравнения на собственные значения.
4. Численно решается уравнение, находятся собственные значения.
5. Численно определяются коэффициенты моды.
6. Собственные значения и коэффициенты подставляются в выражение для нормы и производится пересчет коэффициентов для получения единичной нормы.

На выходе модуля: массивы собственных значений и параметров радиальных функций из (15).

Модуль расчета коэффициентов Галеркина

Коэффициенты Галеркина представляют собой интегралы по объему жидкого ядра от очень громоздких мультипликативных комбинаций базисных мод и операторов векторного анализа в сферических координатах. Сами моды сложным образом выражаются через сферические гармоники (по переменным $\theta$ и $\phi$ ) и сферические функции Бесселя (по радиальной переменной $r$ ). Все это в совокупности приводит к очень сложным подынтегральным выражениям, которые проблематично даже безошибочно ввести вручную в код расчетной программы.

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

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

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

Схема работы модуля:

1. На основе значений мультииндексов и общих формул (14) программно вычисляются алгебраические выражения мод $v_l$, $T_s$, $B_p$ с неопределенными радиальными функциями $R_l\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}$, $Z_s\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}$, $X_p\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}$. Зависимости от $\theta$ и $\phi$ в этих выражениях будут заданы явно.
2. Программно формируются подынтегральные выражения для коэффициентов Галеркина (8).
3. Выполняется аналитическое интегрирование по поверхности сферы $\mathrm{<mo\; stretchy="false">(</mo>}\theta \mathrm{,}\phi \mathrm{<mo\; stretchy="false">)</mo>}$. Результатом этого интегрирования будет либо нуль, либо выражения, зависящие от переменной $r$, неопределенных функций $R_l\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}$, $Z_s\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}$, $X_p\mathrm{<mo\; stretchy="false">(</mo>}r\mathrm{<mo\; stretchy="false">)</mo>}$ и их производных.
4. В случае ненулевого результата подставляются явные выражения для радиальных функций (15) с числовыми значениями параметров мод и выполняется численное интегрирование по радиальной переменной.

На выходе модуля: массивы числовых значений коэффициентов Галeркина.

Необходимо сделать следующее замечание. Каждый коэффициент Галеркина имеет смысл меры взаимодействия мод в процессе, связанном с соответствующим членом уравнений геодинамо. Например $W_{sij}$ $–$ это мера эффективности генерации моды $B_s$ из моды $B_j$ за счет моды скорости $v_i$. Аналитическое интегрирование позволяет определять точно нулевые коэффициенты, что возможность выделять цепочки взаимодействующих мод.

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

Модуль генерации шума

Прежде всего опишем формулы, с помощью которых в модуле генерируются реализации случайных величин, определяющих случайный процесс $\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$. Это три типа величин, которые по определению процесса независимы между собой и, соответственно, могут генерироваться независимо друг от друга: случайное время ожидания $k$ -го прямоугольного импульса ${\tau }_k^W$, случайное время существования этого импульса ${\tau }_k^E$ и случайная амплитуда импульса ${\xi }_k$.

Для генерации случайных времен используется метод обратных функций [12]. Легко получить по схеме этого метода, что величины ${\tau }_k^W$ и ${\tau }_k^E$ при заданных плотностях распределения $p_W\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ (11) и $p_E\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ (12) можно генерировать по формулам

$\begin{array}{ll}{\tau }_k^W\mathrm{<mo>=</mo>}c{T}^W\left(U_1^{\mathrm{1<mo>/</mo><mo\; stretchy="false">(</mo>1}-\nu \mathrm{<mo\; stretchy="false">)</mo>}}-1\right)\text{и}{\tau }_k^E\mathrm{<mo>=</mo>}-\frac{T^E}{\ln 2}\ln U_2\mathrm{,}\hfill & \hfill \end{array}$ (16)

где $U_1$ и $U_2$ равномерно распределены на отрезке $\left[\mathrm{0<mo>;</mo>1}\right]$ и независимы между собой и при различных $k$.

Для генерации случайных амплитуд ${\xi }_k$ используется стандартный метод Бокса-Мюллера [12]:

$\begin{array}{ll}{\xi }_k\mathrm{<mo>=</mo>}\sigma \sqrt{-2\ln U_3}\cos \mathrm{<mo\; stretchy="false">(</mo>2}\pi U_4\mathrm{<mo\; stretchy="false">)</mo>,}\hfill & \hfill \end{array}$ (17)

где $U_3$ и $U_4$ также равномерно распределены на отрезке $\left[\mathrm{0<mo>;</mo>1}\right]$ и независимы между собой и при различных $k$.

Реализующие численное решение уравнений модели (7)$–$(10) модули комплекса выполняют расчет на временной разностной сетке, причем на каждом шаге вычисляются значения решений в очередном временном узле $t_{n+1}$ по их значениям в одном или нескольких предшествующих узлах. Поэтому модуль генерации шума $\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ должен формировать реализацию ${\xi }^{n+1}\mathrm{<mo>=</mo>}\xi \mathrm{<mo\; stretchy="false">(</mo>}{t}_{n+1}\mathrm{<mo\; stretchy="false">)</mo>}$ по уже известной ${\xi }^n\mathrm{<mo>=</mo>}\xi \mathrm{<mo\; stretchy="false">(</mo>}{t}_n\mathrm{<mo\; stretchy="false">)</mo>}$, т.е. по схеме «от текущего $\xi$ в момент $t$, к следующему ${\xi }^{next}$ в момент $t^{next}\mathrm{<mo>=</mo>}t+h$». Именно эта схема и реализуется в модуле генерации шума.

На входе модуль получает следующие параметры и переменные:

• медианное значение $T^W$, порядок $\nu$ и параметр $c\mathrm{<mo>=</mo>}{\left(2^{\frac{1}{\nu -1}}-1\right)}^{-1}$ степенного распределения (11) времени ожидания импульса;
• медианное значение $T^E$ экспоненциального распределения (12) времени существования импульса $T^E$;
• стандартное отклонение $\sigma$ нормального распределения амплитуды импульса;
• булеву переменную текущего состояния процесса pulse (нет импульса $–$ «false», есть импульс $–$ «true»);
• текущее состояние процесса $\xi$;
• ближайшее время переключения switch_time;
• временной шаг $h$ между моментами $t$ и $t^{next}$.

На выходе модуль возвращает:

• булеву переменную состояния процесса pulse в момент $t^{next}$ (возможно, измененную по сравнению со входной);
• следующее состояние процесса ${\xi }^{next}$;
• ближайшее (возможно, измененное по сравнению со входным) время переключения switch_time.

Модуль реализован на C++. Псевдокод его основной части представлен в Листинге 1. Предполагается, что функции rnd_ ${\tau }^E\mathrm{<mo\; stretchy="false">(</mo>}{T}^W\mathrm{<mo\; stretchy="false">)</mo>}$, rnd_ ${\tau }^W\mathrm{<mo\; stretchy="false">(</mo>}c\mathrm{,}{T}^W\mathrm{,}\nu \mathrm{<mo\; stretchy="false">)</mo>}$ и rnd_ $\xi \mathrm{<mo\; stretchy="false">(</mo>}\sigma \mathrm{<mo\; stretchy="false">)</mo>}$ генерируют реализации ${\tau }_k^E$, ${\tau }_k^W$ и ${\xi }_k$ по формулам (16) и (17).

 

Листинг 1. Псевдокод генерации значения шума ${\xi }^{next}$ в момент $t^{next}$ по известному значению $\xi$ в момент $t$.

while switch_time ≤ t+h then

if pulse then

ξ^next := 0 

switch_time := switch_time+rnd_τ^W(c,T^W,ν)

pulse := false

else

ξ^next := rnd_ξ(σ)

switch_time := switch_time+rnd_τ^E(T^E)

pulse := true

end if

end do

 

Реализация переключений процесса в цикле связана с тем, что на малом отрезке [t; t+h] возможны, хотя и маловероятны, несколько таких переключений, поскольку величины τ^E и τ^W теоретически могут принимать как угодно близкие к нулю значения. Однако понятно, что в подавляющем большинстве случаев данный цикл на каждом шаге решения по времени будет отрабатывать не более одного раза.

Модуль численного решения уравнений модели для экспоненциального ядра

Модуль реализован на C++ и предназначен для решения уравнений модели (7)$–$(10) для ядер подавления $K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{M}_n{t}^n\exp \mathrm{<mo\; stretchy="false">(</mo>}-t\mathrm{<mo\; stretchy="false">)</mo>,}n\mathrm{<mo>=</mo>0,1,2,}\dots$, имеющих экспоненциальный порядок убывания на бесконечности. Здесь $M_n$ $–$ нормировочный коэффициент, определяемый условием ${{\displaystyle \max }}_{t\ge 0}K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1}$. Такие ядра далее будем называть экспоненциальными (рис. 2).

 

Рис. 2. Экспоненциальные ядра K(t) =M_ntⁿ exp(−t).

[Figure 2. Exponential kernels K(t) =M_ntⁿ exp(−t).]

 

Параметр $n$ определяет задержку отклика $\alpha$ -эффекта на подавление. Действительно, при $n\mathrm{<mo>=</mo>0}$ отклик мгновенный, поскольку $K\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>1}$ является наибольшим значением функции $K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$, и максимальный вклад в интеграл подавления в момент времени $t$ дает значение энергии магнитного поля в этот же момент. Если $n\ne 0$, то $K\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>0}$ и на подавление в момент $t$ влияют только предшествующие значения энергии поля, т.е. имеет место задержка отклика. Эта задержка тем больше, чем больше $n$, поскольку максимум функции $K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ достигается в точке $t\mathrm{<mo>=</mo>}n$. Кроме того, чем больше $n$, тем ближе к нулю значения $K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ в окрестности $t\mathrm{<mo>=</mo>0}$. Все эти рассуждения хорошо иллюстрируются рис. 2, на котором изображены функции $M_n{t}^n\exp \mathrm{<mo\; stretchy="false">(</mo>}-t\mathrm{<mo\; stretchy="false">)</mo>}$ для нескольких $n$.

Несложно показать, подобно работам [13, 15], что в случае ядер данного вида интегральное выражение (9), определяющее $u\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$, оказывается равносильным следующей задаче Коши для функции $u\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$:

$\begin{array}{l}{\left(T_K\frac{d}{dt}+1\right)}^{n+1}u\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}n\mathrm{<mo>!</mo>}{\displaystyle \sum _{p\mathrm{<mo>=</mo>1}}^{P_{max}}}{\gamma }_p^2\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}u\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{u}^{\prime }\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}\dots \mathrm{<mo>=</mo>}{u}^{\mathrm{<mo\; stretchy="false">(</mo>}n\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>0.}\hfill \end{array}$ (18)

Стандартным образом вводя переменные $u_m\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{u}^{\mathrm{<mo\; stretchy="false">(</mo>}m\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$, $m\mathrm{<mo>=</mo>0,1,}\dots \mathrm{,}n$, получаем тогда, что модель (7)$–$(9) с заданными начальными значениями для амплитуд полей (10) становится следующей задачей Коши для нормальной системы уравнений порядка $L_{\max }+S_{\max }+P_{\max }+n+1$ со случайным возмущением $\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$:

$\begin{array}{ll}\frac{d{\beta }_l}{dt}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{,}j\mathrm{<mo>=</mo>1}}^{L_{\max }}}{B}_{lij}{\beta }_i{\beta }_j-{\mu }_l{\beta }_l+{\text{E}}^{-1}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^{L_{\max }}}{E}_{li}{\beta }_i+{\text{RaPr}}^{-1}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^{S_{\max }}}{C}_{li}{\alpha }_i+{\displaystyle \sum _{i\mathrm{,}j\mathrm{<mo>=</mo>1}}^{P_{\max }}}{Q}_{lij}{\gamma }_i{\gamma }_j\mathrm{,}\hfill & \hfill \\ {\beta }_l\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{\beta }_l^0\mathrm{,}l\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{L}_{\max }\mathrm{,}\hfill & \hfill \\ \frac{d{\alpha }_s}{dt}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{,}j\mathrm{<mo>=</mo>1}}^{L_{\max }\mathrm{,}{S}_{\max }}}{F}_{sij}{\beta }_i{\alpha }_j+{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^{L_{\max }}}{H}_{si}{\beta }_i-{\text{Pr}}^{-1}{\lambda }_s{\alpha }_s\mathrm{,}\hfill & \hfill \\ {\alpha }_s\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{\alpha }_s^0\mathrm{,}s\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{S}_{\max }\mathrm{,}\hfill & \hfill \\ \frac{d{\gamma }_p}{dt}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{,}j\mathrm{<mo>=</mo>1}}^{L_{\max }\mathrm{,}{P}_{\max }}}{W}_{pij}{\beta }_i{\gamma }_j+\frac{{\text{R}}_{\alpha }\left(1+\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}\right)}{1+u_0}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^{P_{\max }}}{W}_{pij}^{\alpha }{\gamma }_i-{\text{Pm}}^{-1}{\eta }_p{\gamma }_p\mathrm{,}\hfill & \hfill \\ {\gamma }_p\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{\gamma }_p^0\mathrm{,}p\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{P}_{\max }\mathrm{,}\hfill & \hfill \\ \frac{d{u}_m}{dt}\mathrm{<mo>=</mo>}{u}_{m+1}\mathrm{,}\hfill & \hfill \\ u_m\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>0,}m\mathrm{<mo>=</mo>0,}\dots \mathrm{,}n-\mathrm{1,}\hfill & \hfill \\ \frac{d{u}_n}{dt}\mathrm{<mo>=</mo>}-{\displaystyle \sum _{m\mathrm{<mo>=</mo>0}}^n}{T}_K^{m-n-1}{C}_{n+1}^m{u}_m+\frac{n\mathrm{<mo>!</mo>}}{T_K^{n+1}}{\displaystyle \sum _{p\mathrm{<mo>=</mo>1}}^{P_{max}}}{\gamma }_p^2\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}\hfill & \hfill \\ u_n\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>0,}\hfill & \hfill \end{array}$ (19)

где $C_{n+1}^m$ $–$ биномиальные коэффициенты.

Модуль строит численное решение задачи Коши (19) с помощью схемы типа «предиктор$–$корректор» Адамса-Башфорта-Моултона 4-го порядка. Для краткости запишем эту задачу в следующем общем виде:

$\begin{array}{l}\frac{dx}{dt}\mathrm{<mo>=</mo>}f\left(x\mathrm{,}\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}a\right)\mathrm{,}x\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{x}^0\mathrm{,}\hfill \end{array}$ (20)

где фазовый вектор $x$ состоит из амплитуд ${\alpha }_s\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$, ${\beta }_l\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$, ${\gamma }_p\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ и переменных $u_m\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$, а вектор параметров $a$ включает в себя все параметры подобия, все параметры процесса $\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ и параметры ядра, все коэффициенты Галeркина.

Пусть $f^k\mathrm{<mo>=</mo>}f\left(x^k\mathrm{,}{\xi }^k\mathrm{,}a\right)$. Тогда реализованная в модуле разностная схема Адамса-Башфорта-Моултона 4-го порядка записывается в виде [14]:

$\begin{array}{l}{x}^{pr}\mathrm{<mo>=</mo>}{x}^n+\frac{h}{24}\left(55f^n-59f^{n-1}+37f^{n-2}-9f^{n-3}\right)\mathrm{,}\hfill \\ x^{n+1}\mathrm{<mo>=</mo>}{x}^n+\frac{h}{24}\left(9f\left(x^{pr}\mathrm{,}{\xi }^{n+1}\mathrm{,}a\right)+19f^n-5f^{n-1}+f^{n-2}\right)\mathrm{.}\hfill \end{array}$ (21)

Поскольку данная схема является 4-шаговой, необходимо до ее применения определить $x^1$, $x^2$ и $x^3$. Они вычисляются в модуле по стандартной явной схеме Рунге-Кутта 4-го порядка, причем начальное значение шума ${\xi }^0\mathrm{<mo>=</mo>0}$.

На входе модуль получает следующие параметры из текстового файла:

• параметры $L_{max}$, $S_{max}$, $P_{max}$, параметры подобия $\text{E}$, $\text{Ra}$; $\text{Pr}$, $\text{Pm}$, ${\text{R}}_{\alpha }$; ненулевые коэффициенты Галеркина (все прочие по умолчанию полагаются нулевыми);
• параметры ядра подавления $T_K$, $n$;
• параметры $T^E$, $T^W$, $\nu$, $\sigma$ процесса $\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$;
• массивы начальных значений ${\alpha }_s^0$, ${\beta }_l^0$, ${\gamma }_p^0$;
• общее время моделирования $T_{max}$ и шаг по времени $h$.

На выходе модуль формирует текстовый файл, каждая строка которого содержит:

• отсчет времени $t_n$;
• отсчеты амплитуд мод скорости ${\beta }_l^n$, $l\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{L}_{max}$;
• отсчеты амплитуд мод температуры ${\alpha }_s^n$, $s\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{S}_{max}$;
• отсчеты амплитуд магнитных мод ${\gamma }_p^n$, $p\mathrm{<mo>=</mo>1,}\dots \mathrm{,}{P}_{max}$;
• отсчет шума ${\xi }^n$;
• отсчет интегрального члена $u^n$;
• отсчет переменной интенсивности $\alpha$ -эффекта ${\text{R}}_{\alpha }\left(1+{\xi }^n\right)\mathrm{<mo>/</mo>}\left(1+u^n\right)$.

Модуль численного решения уравнений модели для степенного ядра

Модуль реализован на C++ и предназначен для решения уравнений модели (7)-(9) с начальными значениями (10) при ядрах подавления

$\begin{array}{l}K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{M_{\rho \mathrm{,}\psi }{t}^{\rho }}{{\mathrm{<mo\; stretchy="false">(</mo>1}+t\mathrm{<mo\; stretchy="false">)</mo>}}^{\rho +\psi }}\mathrm{,}\rho \ge \mathrm{0,}\psi \mathrm{<mo>></mo>0,}\hfill \end{array}$ (22)

где $M_{\rho \mathrm{,}\psi }$ $–$ нормировочный коэффициент, определяемый условием ${{\displaystyle \max }}_{t\ge 0}K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>1}$. Ясно, что параметр $\psi$ определяет порядок асимптотики ядра, поскольку $K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}\sim \mathrm{1<mo>/</mo>}{t}^{\psi }$ при $t\to +\infty$. Параметр же $\rho$ определяет задержку отклика $\alpha$ -эффекта на подавление. Это следует из рассуждений, полностью аналогичных рассуждениям о роли параметра $n$ в экспоненциальных ядрах, а максимум функции $K\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ достигается в точке $t\mathrm{<mo>=</mo>}\rho \mathrm{<mo>/</mo>}\psi$. Иллюстрирующие примеры ядер приведены на рис. 3.

 

Рис. 3. Степенные ядра K(t) =M_ρ,ψt^ρ/(1+t)^ρ+ψ: вверху ψ = 0.5, внизу ψ = 2.

[Figure 3. Power kernels K(t) =M_ρ,ψt^ρ/(1+t)^ρ+ψ: ψ = 0.5 top, ψ = 2 bottom.]

 

Для описания реализованной в модуле разностной схемы запишем уравнения (7)-(9) с условиями (10) в следующем общем виде

$\begin{array}{l}\frac{dx}{dt}\mathrm{<mo>=</mo>}f\left(x\mathrm{,}u\mathrm{,}\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}a\right)\mathrm{,}u\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\int }_0^t}K\left(\frac{t-\tau }{T_K}\right)q\left(x\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>,}a\right)d\tau \mathrm{,}x\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>}{x}^0\mathrm{,}\hfill \end{array}$ (23)

где фазовый вектор $x$ состоит из амплитуд ${\alpha }_s\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$, ${\beta }_l\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$, ${\gamma }_p\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$, а вектор параметров $a$ включает в себя все управляющие параметры и коэффициенты Галeркина.

Модуль численно решает задачу (23) с помощью комбинированной схемы типа «предиктор$–$корректор». Для дифференциальной части задачи в основе схемы лежит метод Коши-Эйлера, а для интегральной $–$ метод трапеций.

Отметим прежде всего, что для каждого момента времени $t_n$ должны быть вычислены значения ядра $K\left(\frac{nh}{T_K}\right)$ и записаны в одномерный массив $K^n$. Размер массива определяется отношением времени моделирования к шагу $h$.

Начальное значение фазового вектора $x^0$ известно, начальные значения $u^0\mathrm{<mo>=</mo>0}$, $q^0\mathrm{<mo>=</mo>}q\left(x^0\mathrm{,}a\right)$, ${\xi }^0\mathrm{<mo>=</mo>0}$.

Рассмотрим теперь основной расчетный этап $—$ переход от момента $t_n$ к моменту $t_{n+1}$. Будем считать, что для момента времени $t_n$ известны $x^n$, ${\xi }^n$ и $u^n$, а также записаны на всех предыдущих шагах вычислений элементы массива $q^k\mathrm{<mo>=</mo>}q\left(x^k\mathrm{,}a\right)$ для $k\le n$.

Шаги этапа следующие:

1. Вычисляем по формуле трапеций интеграл ${\displaystyle {\int }_0^{t_n}}K\left(\frac{t_{n+1}-\tau }{T_K}\right)q\left(x\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>,}a\right)d\tau$:
   $\begin{array}{l}I\mathrm{<mo>=</mo>}\left\{\begin{array}{c}\mathrm{0,}\text{еслиn=0;}\\ \frac{h}{2}{K}^{n+1}{q}^0+h{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^{n-1}}{K}^{n+1-i}{q}^i+\frac{h}{2}{K}^1{q}^n\mathrm{,}\text{еслиn}\ne \text{0}\text{.}\end{array}\right.\hfill \end{array}$
2. Вычисляем предикторы:
   $\begin{array}{l}{x}^{pr}\mathrm{<mo>=</mo>}{x}^n+hf\left(x^n\mathrm{,}{u}^n\mathrm{,}{\xi }^n\mathrm{,}a\right)\mathrm{,}{u}^{pr}\mathrm{<mo>=</mo>}I+\frac{h}{2}{K}_1{q}^n+\frac{h}{2}{K}_{0}q\left(x^{pr}\mathrm{,}a\right)\mathrm{.}\hfill \end{array}$
3. C помощью модуля генерации шума формируем ${\xi }^{n+1}\mathrm{<mo>=</mo>}{\xi }^{next}$ по известному ${\xi }^n$
4. Вычисляем корректоры:
   $\begin{array}{l}{x}^{n+1}\mathrm{<mo>=</mo>}{x}^n+\frac{h}{2}\left[f\left(x^n\mathrm{,}{u}^n\mathrm{,}{\xi }^n\mathrm{,}a\right)+f\left(x^{pr}\mathrm{,}{z}^{pr}\mathrm{,}{\xi }^{n+1}\mathrm{,}a\right)\right]\mathrm{,}\hfill \\ q^{n+1}\mathrm{<mo>=</mo>}q\left(x^{n+1}\mathrm{,}a\right)\mathrm{,}\hfill \\ u^{n+1}\mathrm{<mo>=</mo>}I+\frac{h}{2}\left[K^1{q}^n+K^0{q}^{n+1}\right]\mathrm{.}\hfill \end{array}$
5. Переходим к следующему этапу.

Видно, что формула для $x^{pr}$ $–$ это шаг явного метода Эйлера, а формула для $u^{pr}$ представляет собой формулу трапеций для интеграла ${\displaystyle {\int }_0^{t_{n+1}}}K\left(\frac{t_{n+1}-\tau }{T_K}\right)q\left(x\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>,}a\right)d\tau$, где в последнем слагаемом вместо неизвестного $x^{n+1}$ использован $x^{pr}$. Формула для $x^{n+1}$ $–$ это шаг неявного метода Эйлера с использованием предикторов вместо неизвестных величин. А формула для $u^{n+1}$ $–$ это обычная формула трапеций для интеграла, выражающего $u\mathrm{<mo\; stretchy="false">(</mo>}{t}_{n+1}\mathrm{<mo\; stretchy="false">)</mo>}$. Поэтому и получается, что используемая расчетная схема комбинирует метод Коши-Эйлера для дифференциальной части системы (23) и формулы трапеций для интегральной части этой системы.

На входе модуль получает следующие параметры из текстового файла:

• параметры $L_{max}$, $S_{max}$, $P_{max}$, параметры подобия $\text{E}$, $\text{Ra}$; $\text{Pr}$, $\text{Pm}$, ${\text{R}}_{\alpha }$; ненулевые коэффициенты Галеркина (все прочие по умолчанию полагаются нулевыми);
• параметры ядра подавления $T_K$, $\rho$, $\psi$;
• параметры $T^E$, $T^W$, $\nu$, $\sigma$ процесса $\xi \mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$;
• массивы начальных значений ${\alpha }_s^0$, ${\beta }_l^0$, ${\gamma }_p^0$;
• общее время моделирования $T_{max}$ и шаг по времени $h$.

На выходе модуль формирует текстовый файл, структура которого полностью идентична структуре выходного файла предыдущего модуля.

Заключение

Авторами разработан комплекс программ для моделирования геодинамо в рамках класса спектральных моделей с эредитарным подавлением турбулентного генератора магнитного поля (α-эффекта).

В основе моделей лежат построения галеркинских аппроксимаций путем разложения полей задачи по собственным модам свободных затуханий. Поэтому первая группа задач, решаемых комплексом, $—$ это расчет параметров базисных мод и расчет коэффициентов Галеркина. Эти задачи решаются двумя модулями комплекса за счет комбинированных численно-аналитических вычислений.

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

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

About the authors

Gleb M. Vodinchar

Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS

Author for correspondence.
Email: gvodinchar@ikir.ru
ORCID iD: 0000-0002-5516-1931

PhD (Phys & Math), Associate Professor, Leading Researcher

Russian Federation, 684034, Paratynka, Kamchatka, Mirnaya str., 7

Evgeny A. Kazakov

Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS

Email: gvodinchar@ikir.ru
ORCID iD: 0000-0001-7235-4148

Junior Researcher

Russian Federation, 684034, Paratynka, Kamchatka, Mirnaya str., 7

Liybov K. Feshchenko

Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS

Email: gvodinchar@ikir.ru
ORCID iD: 0000-0001-5970-7316

PhD (Phys & Math), Researcher

Russian Federation, 684034, Paratynka, Kamchatka, Mirnaya str., 7

References

Supplementary files