Методы решения моделирующих уравнений

О модели /

Методы решения моделирующих уравнений можно разделить на

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

В первом случае решение ищется в виде

y(r,\,\theta,\,\lambda,\,t)=\sum_n\sum_my_n^m(r)P_n^m(\cos\theta)\exp^{im(\Omega t+\lambda)},\,

или в виде

y(r,\,\theta,\,\lambda,\,t)=\sum_n\sum_m\tilde y_n^m(r,\,t)P_n^m(\cos\theta)\exp^{im\lambda},\,

где P_n^m(\cos\theta) — присоединенные полиномы Лежандра.

Применением метода Галеркина каждое уравнение исходной системы моделируемых уравнений трансформируется в систему дифференциальных уравнений для коэффициентов разложения y_n^m(r) или \tilde y_n^m(r,\,t),\, которые затем уже решаются разностными методами. Спектральные методы в ряде случаев более удобны для анализа и интерпретации результатов, однако для их применения требуется предварительная линеаризация уравнений (или организации громоздкой итерационной процедуры по нелинейности), а для воспроизведения резких пространственных неоднородностей термосферы и ионосферы, связанных с высыпаниями частиц и магнитосферной конвекцией, требуется неприемлемо большое число членов разложения, что делает предпочтительными чисто разностные методы.

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

\large{\large\partial y_\alpha\over\large\partial t}=A_\alpha{\large\partial^2y_\alpha\over\large\partial r^2}+B_\alpha{\large\partial y_\alpha\over \large\partial t}+D_\alpha{\large\partial y_\alpha\over\large\partial\lambda}+E_\alpha y_\alpha+F_\alpha,\,

где коэффициенты A_\alpha\,,B_\alpha\,,C_\alpha зависят в общем случае от неизвестных функций y_\alpha и их производных, т.е. уравнения являются линейными. Коэффициенты уравнений на каждом шаге интегрирования по времени вычисляются по уже найденным значениям неизвестных функций для предыдущего момента времени или предыдущей итерации.

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

\large{\large\partial y^{(1)}\over\large\partial t}=A{\large\partial^2y^{(1)}\over\large\partial r^2}+B{\large\partial y^{(1)}\over\large\partial r}+E^{(1)}y^{(1)}+F^{(1)},\,\hspace{10pt}(5.4)

{{\large\partial y^{(2)}\over\large\partial t}}=C{\large\partial y^{(2)}\over\large\partial\theta}+E^{(2)}y^{(2)}+F^{(2)},\,

{{\large\partial y^{(3)}\over\large\partial t}}=D{\large\partial y^{(3)}\over\large\partial\lambda}+E^{(3)}y^{(2)}+F^{(3)},\,

E^{(1)}+E^{(2)}+E^{(3)}=E,\,

F^{(1)}+F^{(2)}+F^{(3)}=F,\,

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

Одномерное уравнение второго порядка по вертикальной координате вида (5.4) являются наиболее распространенным типом уравнений в задачах ионосферного моделирования. Если переходить к разностному аналогу уравнения (5.4) путем непосредственной замены в нем производных на отношения конечных разностей, то окажется, что в каждой ячейке пространственной сетки не будут выполняться интегральные законы сохранения величин \rho,\,\rho v,\,\rho v^2 (появятся фиктивные источники), т.е. разностная схема будет неконсервативной, что может приводить, вообще говоря, к ее расходимости.

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

a{\large \partial y\over\large \partial t}={\large \partial\over\large \partial r}\left[b{\large \partial y\over\large \partial r }+cy\right]-dy+f,\,\hspace{40pt}(5.8)

Аппроксимация уравнения (5.8) соответствующей консервативной схемой имеет вид:

{\large\bar a_l}{\large\bar y_l-y_l\over\large\tau}={1\over\large h}\left[g_{l+1}{\large\bar y_{l+1}-\bar y_l\large\over h}-g_l{\large\bar y_l-\bar y_{l-1}\over\large h}\right]+{\large c_{l+1}\bar y_{l+1}-c_{l-1}\bar y_{l-1}\over\large h}-d_l\bar y_l+f_l.\hspace{10pt}(5.9)

Здесь используется разностная сетка (рис. 34), где

  • \tau=t_{k-1}-t_k — шаг интегрирования по времени (k=0,1,…, M — номер временного слоя);
  • h=z_{l+1}-z_l — шаг интегрирования по высоте (l=0,1,…,N — номер высотного слоя)

\bar y_l\equiv y_l^{k+1}\equiv y(t_{k+1},r_l),\hspace{10pt}y_l\equiv y_l^{k}\equiv y(t_{k},r_k)\,.\hspace{10pt}(5.10)

Аналогичный смысл имеют обозначения a_l,\,b_1,\,c_1,\,d_1,\,e_1,\,f_1.

{\large g_{l+1}={\large b_{l+1}+b_l\over\large 2},\,\hspace{20pt}g_l={\large b_l+b_{l-1}\over\large 2}}.\hspace{10pt}(5.11)

Тот факт, что в правой части (5.9) стоят \bar y,\, а не y означает, что схема (5.9) неявная. В противном случае явной схемы значения искомой функции \bar y в текущий момент времени явным образом выражались бы через y — значения искомой функции в предыдущий момент времени. Преимущество неявных схем над явными заключается в том, что последние устойчивы лишь при достаточно малом шаге сетки по времени \tau, удовлетворяющем в случае уравнения второго порядка неравенству

\tau\lesssim{h^2\over{b/a}},\,\hspace{20pt}(5.12)

а в случае, когда b=0 и (5.9) становится уравнением первого порядка

\tau\lesssim{h^2\over{c/a}},\,\hspace{20pt}(5.13)

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

Зависящие от y коэффициенты в (5.9) вычисляются по значениям искомых функций, взятым с предшествующего временного слоя или предыдущей итерации. Соответствующее схеме (5.9) итерационное уравнение может быть записано в следующем трехточечном виде:

A_l^s\bar y_{l-1}^{s+1}-C_l^s\bar y_l^{s+1}+B_l^s\bar y_{l+1}^{s+1}=-F_l^s,\,\hspace{10pt}l=1,\,\ldots,\,N-1,\,\hspace{20pt}(5.14)

где

A^s_l={g^s_l-{{1\over 2}{h}c^s_{l-1}}\over h^2},\,\hspace{20pt}

B^s_l={{g^s_{l+1}+{1\over 2}hc^s_{l+1}}\over h^2},\,

C^s_l={g^s_l+g^s_{l+1}\over h^2}+d^s_l+{a^s_l\over\tau},\,

F^s_l=f_l^s+{a^s_ly_l\over \tau},\,\hspace{10pt}(5.15)

s=0,\ 1,\ 2\ \ldots\ — номер итерации.

Система алгебраических уравнений (5.14) линейна относительно \bar{y}^{s+1}_l . За нулевое приближение берутся значения y_l с предыдущего временного слоя. Итерационный процесс заканчивается по заданному числу итераций или заданной точности: \max|\bar y^{s+1}_l-\bar y^s_l|<\varepsilon .

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

\{\begin{array}{c c}B_0\bar{y_1}-C_0\bar{y_0}=-F_0&r=r_0\hspace{20pt}(5.16)\\A_N\bar{y_{N-1}}-C_N\bar{y_N}=-F_N&r=r_n\hspace{20pt}(5.17)\end{array}

и решается одним из вариантов метода прогонки. Решение в этом методе ищется в виде

\bar y^{s+1}_l=\alpha^s_{l+1}\bar y^{s+1}_{l+1}+\beta^s_{l+1}\hspace{10pt}(5.18)  .

Для нахождения прогоночных коэффициентов \alpha и \beta используют рекуррентные формулы

\alpha_{l+1}={B_l\over C-l-A_l\alpha_l},\,\hspace{20pt}(5.19)

\beta_{l+1}={A_l\beta_l+F_l\over C_l-A_l\alpha_l},\,\hspace{30pt}(5.20)

\alpha_1={B_0\over C_0},\,\hspace{25pt}\beta_1={F_0\over C_0}\hspace{25pt}(5.21)

Здесь индекс {s} нумерующий итерации, опущен.

Обратная прогонка — нахождение неизвестных функций {\bar y_l} по формуле (5.18) с помощью рассчитанных предварительно прогоночных коэффициентов начинается со значения {\bar y_N,\,} которое определяется из краевого значения (5.17) и соотношения (5.18) для точки {l=N-1} :

\bar y_{N-1}=\alpha_N\bar y_N+\beta_N\hspace{10pt}(5.22)

Получим

\bar{y_N}=({F_N\over C_N}+{A_N\over C_N}\beta_n)/(1-{A_N\over C_N}\alpha_N)\hspace{20pt}(5.23)

Для обеспечения разрешимости системы (5.13), (5.16), (5.17) и устойчивости метода прогонки необходимо, чтобы коэффициенты системы (5.14) и граничных условий (5.16), (5.17) удовлетворяли неравенствам

A_l>0,\,\hspace{20pt}B_l>0,\,\hspace{20pt}(5.24)

0\leq{B_0\over C_0},\,\hspace{20pt}{A_N\over C_N}\leq1,\,\hspace{20pt}{B_0\over C_0}+{A_N\over C_N}<2\hspace{20pt}(5.25)

Отметим модификации метода прогонки, часто используемые в практике ионосферного моделирования — потоковый вариант и матричная прогонка, а также их комбинации. — потоковый вариант матричной прогонки. Потоковый вариант метода прогонки основан на раздельной разностной аппроксимации одномерного уравнения непрерывности в форме (3,149) и развернутого выражения для вертикального потока ионов {n_iV_{iz},\,} где {V_{iz}} задается выражениями (3,151)-(3,153). Этим обеспечивается большая точность вычисления дивергентного члена уравнения непрерывности и устраняются трудности, связанные с быстрым ростом коэффициента диффузии с высотой. Потоковый вариант особенно удобен при задании в качестве граничных условий значений потоков частиц.

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

Для построения оптимальных разностных сеток, учитывающих рост характерных масштабов явлений с высотой, используют замены независимой пространственной переменной или искомой функции в исходных уравнения с последующей разностной аппроксимацией преобразованных уравнений на равномерной сетке по новой переменной или для новой функции. Например, при интегрировании моделирующих уравнений в дипольных координатах целесообразно от координаты {q,\,} отсчитываемой вдоль силовой линии геомагнитного поля, перейти к новой независимой безразмерной переменной {x} с помощью замены

x={\large{\text sh}(\Gamma q)\over\large {\text sh}(\Gamma q_m)},\,\hspace{20pt}(5.26)

где

  • {q_m} — значение  {q_m} на северном конце геомагнитной силовой линии,
  • \Gamma — численный множитель, определяющий степень сгущения сетки по {q} при равномерном разбиении по {x} .

При интегрировании по вертикали в сферических координатах удобно перейти от координаты {r} к безразмерной переменной {l} с помощью соотношения

\large{\large\partial\over\large\partial r}={1\over\large \Delta r_l}{\large\partial\over\large\partial t},\,\hspace{20pt}(5.27)

где

\Delta{r_l}={r(l)\over N}l\cdot\ln({r_N\over r_0})\hspace{20pt}(5.28)

или

r(l)=r_0\exp^{(\ln{r_N\over r_0})},\,\hspace{20pt}(5.29)

{r_0} и {r_N} — координаты верхней и нижней границы. В этом случае равномерной сетке по {l} соответствует экспоненциально нарастающий шаг по {r} .

При интегрировании уравнений для концентрации нейтральных частиц, экспоненциально убывающих с высотой, целесообразно осуществить замену искомой функции {n\to y=\ln n} .

Выбор закона изменения шага интегрирования по высоте осуществляется обычно таким образом, чтобы в D- и E-областях он составлял 0,5-5 км, в области F1 и в окрестности максимума F2-слоя — 10-20км, в верхней части F2-области — 50-100км. Типичные шаги интегрирования по другим переменным: по широте — 2-10°, по долготе — 15°, по времени — 20 мин для спокойных дневных условий и 2-5 мин для восходно-заходных периодов и возмущенных периодов типа суббурь. Все эти цифры являются ориентировочными. Такие шаги обеспечивают погрешность численного решения в пределах 10%. Окончательный выбор шагов интегрирования в каждой конкретной задаче моделирования осуществляется путем численных экспериментов с последовательным дроблением шагов, проводимых до тех пор, пока решения не перестанут отличаться в пределах заданной погрешности.

Реклама