УДК 551.574.14.001.572

 

 

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

Д.Я. Прессман

Гидрометеорологический научно-исследовательский центр
Российской Федерации, г. Москва, Россия

pressman@mecom.ru

 

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

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

 

 

Coordinate splitting, Riemann’s equations,
grid-characteristic method, alternating direction scheme
and the use of integral conservation laws to correct
the solution of finite-difference approximation of the equation for an ideal baroclinic fluid

D.Ya. Pressman

Hydrometeorological Research Center of Russian Federation, Moscow, Russia

pressman@mecom.ru

 

The system of hyperbolic equations of the dynamical core of the prognostic model for the ideal atmosphere without filtration of any waves is approximated by finite differences. The grid-characteristic method is applied to the Riemann’s equations obtained after the coordinate splitting. The scheme of alternating directions is used, and the solution is corrected using the integral conservation laws for mass, energy and momentum. The flow limited from above by a free material constant-pressure surface propagates over the horizontally periodic relief.

Keywords: finite-difference approximation, ideal baroclinic atmosphere equations, coordinate splitting, Riemann’s equations, grid-characteristic method, alternating direction scheme, solution correction, integral conservation laws

 

 

Введение

При численном решении задачи прогноза погоды мощность современного компьютера позволяет перейти к аппроксимации ‘исходных’ уравнений в надежде устранить искажения, возникающие при фильтрации каких-либо волн. За ‘исходные’ уравнения здесь приняты гиперболические уравнения, описывающие поведение идеальной негидростатической бароклинной среды. Поток, ограниченный сверху свободной материальной поверхностью постоянного давления, движется над заданным рельефом. Делается замена координат, ‘спрямляющая’ подвижную верхнюю и фиксированную нижнюю границы среды. Не учитывается кривизна геоида Земля.

Цель статьи – предложить аппроксимацию прогностической модели с нефильтрованными уравнениями, основанную на максимальном использовании гиперболичности этих уравнений, – свести задачу к ‘переносу’ комбинаций неизвестных функций вдоль характеристик и учесть возможность пересечения характеристик одного семейства. Для этого ‘на дифференциальном уровне’ ‘исходная’ система уравнений расщепляется по координатам [2, 6], после чего в каждой из расщепленных задач делается переход к уравнениям Римана [4]. Уравнения Римана соответствуют характеристикам [4] расщепленных задач (кратной характеристике отвечает равное кратности число уравнений Римана) и имеют вид однородных или неоднородных квазилинейных алгебраических уравнений для производных неизвестных функций вдоль соответствующей характеристики.

На ‘разностном уровне’ применяется аппроксимация, хорошо зарекомендовавшая себя при решении задач газовой динамики. После ‘спуска’ всех функций на одну и ту же прямоугольную равномерную по горизонтальным координатам и равномерную по вертикали (в спрямляющей
системе координат) сетку точек используется сеточно-характерис-тический метод [5] аппроксимации уравнений Римана. Для этого в каждой точке сетки находятся все приходящие в нее к концу временного шага криволинейные характеристики. Положение точек ‘исхода’ этих характеристик, скорости ‘переноса’ вдоль них комбинаций неизвестных функций, величины правых частей и коэффициентов уравнений Римана определяются итерациями, цель которых – получить второй порядок аппроксимации. В ходе итераций величины скоростей переноса, правых частей и коэффициентов каждого из уравнений Римана относятся к середине отрезка соответствующей характеристики, ограниченного точками ее ‘исхода’ и ‘прихода’. На каждой итерации осуществляются упомянутые ‘переносы’ комбинаций неизвестных функций в исследуемую точку сетки, что приводит к системе линейных алгебраических уравнений для определения всех неизвестных функций в этой точке.

В уравнениях Римана для уточнения неизбежных при использовании сеточно-характеристического метода интерполяций функций между точками сетки может быть сделан (и в данной статье он сделан) переход к новым зависимым переменным, за которые будут приняты отклонения от функций  и  фоновой атмосферы. Здесь – отношение теплоемкостей среды при постоянном давлении и постоянном объеме.

Давление  и плотность  фоновой атмосферы получаются осреднением по горизонтальным координатам начальных данных (в принципе, данных какого-либо шага по времени), зависят только от вертикаль-
ной координаты
Z и задаются таблично с максимально возможным разрешением.

В данной статье фон статичен: , где g = 9,8 м/с2
ускорение силы тяжести, и определяется заданием фоновой темпера-
туры .

Для получения второго порядка аппроксимации всей задачи, а не только на каждом этапе расщепления, в алгоритм включена схема переменных направлений [2, 7] при следующей циклической смене координат: (ZXY-ZYX-YZX-YXZ-XYZ-XZY). Каждая тройка отвечает переходу к следующему моменту времени и указывает последовательность координат, относительно которых делается расщепление при этом переходе. Например, тройка ZXY означает, что сначала выделяется Z-координата и решаются соответствующие Z-задаче уравнения Римана, затем –
X-координата и решаются задачи Римана для X-задачи, затем –
Y-координата. Начальными данными для любого из этапов расщепления служат значения неизвестных функций, полученные в конце предыдущего этапа.

Явная сеточно-характеристическая аппроксимация гиперболических уравнений Римана требует для устойчивости счета шаг по времени, определяемый наклоном ‘самых быстрых’ характеристик (выполнение критерия Куранта [3], согласно которому область зависимости разностного решения должна быть не меньше области зависимости решения дифференциальной задачи). Например, в X-задаче должно быть выполнено неравенство , где  – минимальный угол (угол наклона) между X-координатой и характеристикой. Этот наклон меняется по ходу счета вместе с решением.

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

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

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

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

 

1. Дифференциальные уравнения

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

По горизонтальным координатам ставим условие периодичности.

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

В новой системе координат вертикальная скорость . Очевидно, что . Функция  – якобиан матрицы перехода от координат (t,x,y,z) к координатам (.

Легко видеть, что выполнены соотношения:

;    ;

 

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

                  (1)

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

(2)

где  ;

;

;

 – давление;  – плотность; =1,4 – отношение теплоемкости при постоянном давлении к теплоемкости при постоянном объеме сухого воздуха; g – ускорение силы тяжести.

Для полноты к системе (2) следует добавить уравнение (1), а для определения температуры Т – уравнение состояния . Универсальная газовая постоянная 287,05 Дж/кг/К.

 

2. Координатное расщепление

Стандартное координатное расщепление [2, 6] (когда для каждой координаты удается получить отдельную задачу с производными неизвестных функций только по времени и по этой координате) провести не удается, так как в систему уравнений входят произведения производных по  (а именно производных ,  и ) на производные по  и по , входящие в  и .

Поэтому в использованном далее переходе к -задаче считаем  и  функциями-коэффициентами при производных по . В начале каждого временного шага эти коэффициенты нужно найти в точке ‘исхода’ характеристик. Точка ‘прихода’ характеристик – всегда точка разностной сетки. В конце временного шага значения этих коэффициентов в точке ‘прихода’ находятся разностным дифференцированием по x и по y найденного решения. В пределах шага по времени делаются итерации, в ходе которых эти коэффициенты пересчитываются (так же как скорости ‘переноса’ вдоль криволинейных характеристик и положение точки ‘исхода’), заменяясь полусуммой их значений в точках ‘исхода’ и ‘прихода’, что позволяет получить аппроксимацию второго порядка.

 

3. -расщепление

При -расщеплении мы переходим к системе, в которую переменные  и  входят параметрически:

          (3)

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

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

Эти добавки будут компенсироваться соответствующими изменениями правых частей уравнений на этапах X- и Y-расщеплений.

Уравнения Римана получим из первых пяти уравнений системы (3), так как шестое уравнение – обыкновенное дифференциальное уравнение. В матричной записи эти пять уравнений превращаются в

,                                             (4)

где    – вектор-столбец ;

 – вектор-столбец ;

 – единичная матрица;  матрица А:

Очевидны следующие выражения для левых собственных векторов  матрицы  и соответствующих собственных значений  :

;      ;

;  ;

;            ;

;

;

.

Умножая (4) слева на векторы  получим пять уравнений
Римана
, первые два из которых (уравнения (5) и (6)) связывают производные функций вдоль характеристики с наклоном :

;    (5)

;    (6)

                                                                                       (7)

(здесь ), тогда как уравнения (8) и (9) связывают производные функций  вдоль характеристик с наклонами  и :

    +

+        (8)

(здесь ),

          (9)

(здесь ).

Перейдем к разностной аппроксимации полученных уравнений
Римана. Как отмечалось во введении, в данной работе шаг по времени сеточно-характеристической аппроксимации один и тот же внутри каждой тройки циклической смены координат расщепления и определяется наклоном самых ‘быстрых’ характеристик трех задач расщепления.
Задача для -координаты определяет максимальный шаг по времени неравенством , где  
- минимальный угол наклона характеристик -задачи, или эквивалентным неравенством . Аналогичные оценки должны быть получены для задач X- и Y-расщепления. Минимум этих трех оценок используется для определения величина .

Находим точки ‘исхода’ характеристик, приходящих в исследуемую точку сетки в конце указанного шага по времени.

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

Рассмотрим отдельно аппроксимации на границах -задачи.

Так как при  выполнено  и коэффициенты уравнений (5) и (6) не зависят от времени, то из краевого условия  следует, что в -задаче при  функции ,  и W не зависят от времени.

При  снова выполнено , но коэффициенты уравнений (5) и (6) зависят от времени, а из условия = const и четвертого уравнения системы (3) следует:

.                                (10)

Аппроксимация (5), (6) и (10) при  приводит к алгебраическим уравнениям, определяющим ,  и W. При решении этих уравнений на первой итерации величины  и  берутся в начале временного шага, а на следующих – полусуммой значений в начале и конце временного шага. Производные по  в уравнении (10) аппроксимируем односторонне, при этом значения на ближайшем внутреннем -уровне должны уже быть определены.

Перейдем к аппроксимации уравнений Римана для внутренних точек -задачи.

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

Стандартные (фоновые) значения зависят только от координаты Z, и для них желательно иметь таблицы с высоким разрешением. Если фон статичен, то уравнение (7) заменяется на (11):  

                      (11)

(здесь ),

а уравнения (8) и (9) - на (12) и (13):

      ++

   =.    (12)

(здесь ),

+

=,         (13)

(здесь ).

Числитель первого множителя правой части уравнения (11) отрицателен при сверхадиабатическом градиенте температуры фоновой среды.

Для внутренних точек сетки на оси  сеточно-характеристическая аппроксимация уравнений (5) ,(6), (12), (13) приводит к системе четырех линейных алгебраических уравнений, определяющих , , W и  в конце временного шага.

Здесь сеточно-характеристическая аппроксимация требует, прежде всего, определения ‘точки исхода’ кратной характеристики, приходящей со скоростью  в исследуемую точку сетки на оси  (точку прихода) за время . Значения функции  и комбинаций функций , , W, определяемых уравнениями (5) , (6) и (11), ‘переносятся’ из точки исхода вдоль этой характеристики. Затем следует найти 'точки исхода' двух оставшихся характеристик, также приходящих в исследуемую точку на оси  за время , но со скоростями  и . Из этих точек исхода ‘переносятся’ комбинации функций  и , определяемые уравнениями (12) и (13).

Полученные в точке сетки пять перечисленных комбинаций образуют систему алгебраических уравнений для определения неизвестных функций. Если на первой итерации скорости характеристик, коэффициенты и правые части уравнений (5), (6), (12), (13) определялись явно в точках 'исхода', то на последующих итерациях все эти величины определяются полусуммой значений в точках ‘исхода’ и ‘прихода’ характеристик. Таким образом, достигается второй порядок аппроксимации. Конечно, итерации включают аппроксимацию краевых условий.

В проведенных расчетах 5-ти итераций было достаточно для повторения 12-ти из 15-ти десятичных разрядов записи соответствующих чисел в машине.

Существенно, что в (12) и (13) должно быть подставлено = в соответствии с системой (3). Кроме того, из-за изменения в ходе итераций величины  постоянным равномерным уровням -сетки будут соответствовать измененные z-уровни, и, следовательно, параметры фоновой атмосферы на -уровнях модели должны пересчитываться.

Значения  и  находятся из полученных величин  и  и значений  фоновой атмосферы на -уровнях модели.

 

4.  X -расщепление

При -расщеплении мы приходим к системе, в которую переменные  и  входят параметрически:

                                        (14)

Заметим, что искусственная добавка  в правую часть
U-уравнения компенсирует соответствующую добавку в правую часть
U-уравнения на этапе -расщепления. Поэтому величина  при построении -расщепления перенесена в правую часть и не определяет характеристические направления X -задачи, хотя в нее входит производная .

Далее будем поступать аналогично решению -задачи: уравнения Римана получим из первого и четвертого уравнений системы (14), так как остальные уравнения уже имеют форму уравнений Римана.

В матричной записи эти два уравнения превращаются в

,                                        (15)

где  – вектор-столбец  ;  – вектор-столбец;

Е  – единичная матрица;

Матрица .

Левые собственные вектора  матрицы  и соответствующие собственные значения :

 ;           ;

 ;        

Умножая (15) слева на вектора  получим два уравнения Римана X-задачи, связывающие производные функций  вдоль характеристик с наклонами ,  и  :

                                           (16)

(здесь );

                                           (17)

(здесь ).

Снова для уточнения аппроксимации может быть сделан переход
к системе уравнений для отклонений функций  и  от стандартных значений  и фоновой, зависящей только от координаты Z атмосферы. Если фон статичен, то пятое уравнение системы (14) превращается в уравнение (18):

    (18) (здесь ),

а уравнения (16) и (17) – в уравнения (19) и (20):

              (19)

(здесь ),

           (20)

(здесь ).

Числитель первого множителя правой части уравнения (18) отрицателен при сверхадиабатическом градиенте температуры фоновой среды.

Сеточно-характеристическая аппроксимация состоит, прежде всего, в определении ‘точки исхода’ кратной характеристики, приходящей со скоростью  в исследуемую точку сетки (точку прихода) на оси X в конце . Значения функций , W, ,  ‘транспортируются’ в точку прихода в соответствии со вторым, третьим, пятым и шестым уравнениями системы (14). Затем следует найти точки ‘исхода’ двух оставшихся характеристик, также приходящих в исследуемую точку на оси X в конце , но со скоростями  и . Из этих точек исхода ‘переносятся’ комбинации функций  и , определяемые уравнениями (19) и (20).

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

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

В проведенных расчетах 5-ти итераций было достаточно для повторения 12-ти из 15-ти десятичных разрядов записи соответствующих чисел в машине.

Существенно, что в (18), (19) и (20) должно быть подставлено  в соответствии с последним уравнением системы (14). Кроме того, из-за изменения величины  в ходе итераций равномерным уровням -сетки будут соответствовать новые -уровни, и, следовательно, параметры фоновой атмосферы на -уровнях должны пересчитываться.

Значения  и находятся из величин  и  и  значений  фоновой атмосферы на -уровнях модели.

Значения  на нижней границе определены условием обтекания.

 

5. Y-расщепление проводится аналогично X-расщеплению. Поэтому приведем только уравнения системы для отклонений:

;

;

(здесь );

(здесь );

(здесь );

.

Значения  и находятся из величин  и  и таблиц
фоновой атмосферы.

Значения  на нижней границе определены условием обтекания.

6. Коррекция результатов на границах сеточного интервала, над которым пересекаются линии тока

Если линии тока (в задачах, полученных в результате координатного расщепления), выходящие из концов сеточного интервала (i, i+1) на рассматриваемой координате пересекаются над этим интервалом, то на его границах возможна коррекция ранее найденных там величин.

Для этого по формулам статьи [1] определяется скорость  перемещения возникающего в таком интервале контактного разрыва и корректируются значения функций на той границе интервала, к которой этот разрыв приближается. Пусть точка i – такая граница. Тогда в точку i
поступят (полученные с помощью интерполяции) значения функций в точке (i – ), где  – шаг схемы по времени, а
 – шаг сетки вдоль этой оси.

 

7. Коррекция полей с помощью аппроксимации интегральных законов сохранения

Будут использованы следующие, верные для системы (2), интегральные законы сохранения. Потоковые слагаемые не приводятся, так как они пропадают или из-за принятого здесь условия периодичности по горизонтальным координатам или из-за равенства нулю вертикальной скорости  на верхней и нижней границах расчетной области.

Горизонтальные составляющие импульса (сохраняются при плоском рельефе):

                      (21)

                                (22)

Вертикальная составляющая импульса (сохраняется при квазистатическом решении):

   (23)

Масса и энергия среды. Сохраняются и определяются начальными данными задачи: 

                                                                                 (24)

  (25)

‘Энтропийный’ интеграл. Сохраняется на гладких решениях и не убывает в общем случае:

                                                                      (26)

Коррекция решения (через каждые 6 шагов по времени, т. е. в конце цикла смены направлений) проводится следующим образом. Определяются разностное решение задачи, в которой в уравнениях Римана сделан переход к отклонениям от фона (решение №1), и разностное решение
задачи, где такой переход не сделан (решение №2). Поскольку принцип построения разностных схем для получения этих двух решений одинаков, то с точки зрения порядка аппроксимации они неотличимы, и модуль разности этих решений в каждой точке сетки может быть прибавлен (вычтен) без изменения порядка аппроксимации к любому из этих решений.

Коррекцию поля  проведем следующим образом. Найдем интеграл по области от выше указанного модуля разности решений  и обозначим его JJJ. Разность между интегральной массой в конце текущего временного шага для решения №1 и сохраняющейся начальной массой поделим на JJJ. Это отношение (обозначим его М) не зависит от точки сетки и может быть любого знака. В каждой точке сетки из значения  решения №1 вычитается соответствующий этой точке вышеуказанный модуль разности, умноженный на М. Очевидно, что это приводит к выполнению сохранения массы, причем изменение в каждой точке сетки тем меньше, чем ближе друг к другу разностные решения №1 и №2.

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

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

Коррекция составляющих импульса ,  и  аналогична коррекции массы, если предварительно учтены правые части (зависимости от времени импульсных интегралов по области от времени) в уравнениях (21), (22), (23). Эти уравнения интегрируются по времени разностно, так что мы не знаем точных значений интегралов от составляющих импульса в текущий момент времени, и проведенная коррекция носит предварительный характер. Предварительные значения составляющих вектора скорости находятся делением предварительно скорректированных составляющих импульса на ранее полученные значения . Теперь (в текущий момент времени) определим предварительное интегральное значение кинетической энергии.

Из закона сохранения энергии находим ‘точное’, окончательное для текущего момента времени интегральное значение кинетической энергии.  Извлекаем квадратный корень из отношения ‘точного’ значения интегральной кинетической энергии к предварительному значению и в каждой точке сетки умножаем предварительные значения составляющих скорости на этот корень. Так корректируется скорость.

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

Для учета силы Кориолиса применялось физическое расщепление. В конце каждого шага по времени горизонтальные составляющие скорости решения основной задачи (сформулированной как в отклонениях от фона, так и без них) служат начальными данными для системы (27), аппроксимируемой по явной схеме, 

.                                                                                  (27)

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

 

9. Результаты численных экспериментов

Разностная схема сохраняет состояние покоя с точностью, зависящей от стратификации фонового потока. В тестовых расчетах с начальными данными покоящейся среды отклонения от фона при любой стратификации не превосходили 10Е-8. В любых расчетах аппроксимация интегральных законов сохранения выполнялась с точностью не менее, чем до сотых долей процента.

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

Для сокращения времени счета применено OMP-распараллеливание (16 нитей) на 32-процессорном INTEL-XEON. Такое распараллеливание привело к 12-кратному сокращению времени счета.

На приведенных рисунках – результаты расчета задачи (сетка NX=32, NY=128, NZ=20) обтекания симметричного относительно Z-оси, ‘косинусного’ рельефа (высота 2400 м, радиус основания 6 км, центр (i=16; j=11)). Условие периодичности по горизонтальным координатам означает, что поток обтекает равномерную ‘решетку’ таких ‘гор’. Период решетки поперек начального потока равен 32 км (вдоль оси X), а вдоль него (вдоль оси Y) – 128 км. Шаг разностной сетки по горизонтальным координатам равен 1 км. По вертикали шаг сетки равномерный в -системе координат, что соответствует в Z-системе координат ~1 км вне зоны над горой.

Горизонтальная начальная скорость постоянна по всем координатам:  м/с, = 20 м/с, а вертикальная начальная скорость не зависит от вертикали и равна . Верхняя граница в начальный момент горизонтальна и расположена на высоте 20 км. Начальное давление на всех -уровнях равно соответствующему фоновому. Всегда . Нет отклонений давления и плотности от фона в начальный момент. Аппроксимация силы Кориолиса проводится по явной схеме.

Все рисунки с указанным рельефом относятся к 4-часовому сроку прогноза и даны в -системе координат, поток движется слева со скоростью V᷉.

Рис. 1 и 2 относятся к случаю, когда гидростатическая фоновая атмосфера в слое до высоты 2,5 км стратифицирована устойчиво с градиентом температуры (-6 °К/км). Затем до высоты 10,0 км расположен
неустойчивый слой с температурным градиентом (-15 °К/км). Выше
(до высоты 40 км) – устойчивый слой с температурным градиентом (+2 °К/км). Фоновое давление на нижней границе
Z = 0 фона =
100 000
 Па, температура = 273,15 °К. Шаг по времени определяется минимальным наклоном характеристик каждой из расщепленных задач, несколько меняется от итерации к итерации и приблизительно равен
~2,4 c.

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

Влияние силы Кориолиса (рисунки не приводятся) заметнее всего на величине составляющей скорости : ноль меняется на значения в интервале (-6,8 … +2,0), а также на искажении симметрии, особенно там, где конвекция слабее.

На рис. 3 и 4 приведены результаты расчета для случая, когда гидростатическая фоновая атмосфера стратифицирована устойчиво с градиентом температуры 6 °К/км во всей области до высоты 40 км. Фоновое давление на нижней границе = 0 фона =100 000 Па, температура –= 273,15°К. Шаг по времени определяется минимальным наклоном характеристик каждой из расщепленных задач, несколько меняется от итерации к итерации и примерно равен 2,5 с.

 

 

 

 

 

 

 

 

 

а)

 

 

 

 

 

 

 

б)

 

 

 

 

 

 

 

в)

 

Рис. 1. Изолинии центральных (Y, )-сечений; координата Y (слева направо) меняется от 0 до 132км; величина Z(Y) каждого из трех рисунков меняется от высоты горы до ~20км: а) поле температуры: Tmin=104°К в левом экстремуме, затем чередование 7 ложбин с 6 экстремумами с T=120°К; T≈272°К на нижней границе за горой; T≈128°К на верхней границе; изолинии через 16°К; б) составляющая скорости V᷉: V᷉min≈ -10м/с; V᷉max≈ +15м/с; после двух экстремумов +15м/с и -10м/c чередование 12 экстремумов со значениями +10м/с и -10м/с, изолинии через 5м/с; в) составляющая скорости W: Wmin≈ 
-24м/с;
Wmax≈ +11м/с, после двух экстремумов -24м/с и +11м/с чередование 12 экстремумов со значениями -14м/с +11мc;  изолинии через 5м/с.

Fig. 1. Contour lines of central (Y, )-cross sections; Y-coordinate (to the right) changes from 0 up to 132km; the value Z(Y) of every three figures changes from mountain height up to ~20km: a) temperature field: Tmin=104°К in the left extremum; after that alternation of 7 hollows and 6 extrema with T=120°К; T≈272°К at the bottom after the mountain; T≈128°К at the upper boundary; contour intervals are every 16°К; б) V-component of the wind velocity: Vmin≈
-10m/s; V
max≈+15m/s; after two extrema +15м/с and -10м/c - alternation of 12 extrema with the values +10m/s and -10m/s;  contour intervals are every 5m/s; в) W-component of the wind velocity: Wmin≈-24m/s; Wmax≈+11m/s; after two extrema -24m/s and +11m/s – alternation of 12 extrema with the values -14m/s and +11m/s; contour intervals are every 5m/s.

 

 

 

 

 

 

 

 

а)

 

 

 

 

 

 

 

б)

 

Рис. 2. Изолинии (Y,X)- поверхностей, Y-координата (слева направо) меняется от 0 до 132км, X-координата (снизу вверх) меняется от 0 до 32км:
a) (Y,X)-сечение на среднем -уровне составляющей скорости W: Wmin≈ 
-24м/с; Wmax≈ +14м/с; левый экстремум -24м/c, затем чередование 9 вертикальных полос со значениями +10м/с и -14м/с, затем экстремум -16м/с и снова чередование 3 полос; изолинии через 10м/с; б) высота верхней границы потока: Z
top меняется в интервале (19790м, 20170м); очевидны три области: левый экстремум 19790м, затем полосы постепенного подъема к 20170м и полосы спуска к 19940м; изолинии через 30м.

Fig. 2. Contour lines of (Y,X)-cross sections, Y-coordinate (to the right) changes from 0 up to 132km, X-coordinate (upwards) – from 0 up to 32km: а) (Y,X)-cross section at the middle -level of W component of the wind velocity: Wmin≈ 
-24m/s; Wmax≈ +14m/s; left extremum -24m/s, after that alternation of 9 zones with values +10m/s and -14m/s, after that the extremum -16m/s and again the alternation of 3 zones; contour intervals are every 10m/s; b) the height of the flow upper boundary: Ztop belongs to interval (19790m, 20170m); discernable three regions: left extremum 19790, after that the zone of gradually increase up to 20170m and the zone of gradually decrease up to 19940m; contour intervals are every 30m.                           

 

 

Согласно рис. 4a, заметные конвективные ячейки существуют только над горой и на небольшом расстоянии за ней.

Влияние силы Кориолиса (рисунки не приводятся) на течение при устойчивом фоне заметно сильнее, чем при неустойчивом, и состоит в нарушении симметрии потока, во вращении по часовой стрелке оси конвективных ячеек над и за горой и в смещении возмущений верхней границы среды. 

Часто качество разрабатываемого алгоритма проверяется расчетом двумерного (X,Z)-обтекания рельефа, профиль которого был предложен в [8]: , h = 250 м, = 5000 м, L = 4000 м. Так, в статье [9] приведены результаты расчетов, в которых поток ограничен сверху горизонтальной стенкой, так что для подавления отраженных от этой стенки волн верхняя 10 км часть расчетной области занята абсорбирующим слоем. На боковых границах ставились паллиативы условий открытых границ. Рис. 5 в статье [9] дублируется рис. 7a в Z-системе координат этой статьи.

 

 

 

 

 

 

 

 

а)

 

 

 

 

 

 

 

 

б)

 

 

 

 

 

 

в)

 

Рис. 3. Изолинии центральных (Y, )-сечений; координата Y (слева направо) меняется от 0 до 132км; величина Z(Y) каждого из трех рисунков меняется от высоты горы до ~20км: а) поле температуры: Tmin =156°К на верхней границе; Tmax ~270°К на нижней границе вне горы; изолинии через 12°К; б) составляющая скорости V᷉: V᷉min≈ -7м/с на нижней границе перед горой; V᷉max≈+23м/с над горой в экстремуме верхней части рисунка; два экстремума в нижней части рисунка слева направо: 20м/с и 8м/с; изолинии через 2м/с; в) составляющая скорости W: Wmin≈-9.0м/с; Wmax≈+3.0м/с; экстремумы слева направо: +1.5 м/с, -7м/с, +3 м/с: изолинии через 2м/с.  

Fig. 3. Contour lines of central (Y, )-cross sections; Y-coordinate (to the right) changes from 0 up to 132km; the value Z(Y) of every three figures changes from mountain height up to ~20km: a) temperature field: T≈270°К at the bottom after the mountain; T≈156°К at the upper boundary; contour intervals are every 20°К; б) V-component of the wind velocity: Vmin≈ -7m/s at the bottom before the mountain; Vmax≈+25m/s = extremum in the upper fig. part exactly above the mountain; two extrema in the lower part of fig. from the left 20m/s and 8m/s; contour intervals are every 2m/s; в) W-component of the wind velocity; Wmin≈-9.0m/s; Wmax≈+3.0m/s; extrema from the left +1.5m/s, -7m/s, +3m/s; contour intervals are every 2m/s.

 

 

 

 

 

 

 

а)

 

 

 

 

 

 

б)

Рис. 4. Изолинии (Y,X)-поверхностей, Y-координата (слева направо) меняется от 0 до 132км, X-координата (снизу вверх) каждого из двух рисунков – от 0 до 32км: а) (X,Y)-сечение на среднем -уровне составляющей скорости W: Wmin≈-1.9м/с; Wmax≈+1.7м/с; экстремумы слева направо: +1.2м/с, -1.9м/с,  +1.7м/с; изолинии через 0.5м/с; б) высота верхней границы потока: Ztop меняется в интервале (20000м, 20020м); над горой три экстремума слева направо: 20015м, 20000м, 20020м; изолинии через 5м.

Fig. 4. Contour lines of (Y,X)-cross sections, Y-coordinate (to the right) from 0 up to 132km, and X-coordinate (upwards) of every picture– from 0 up to 32km: а) (Y,X)-cross section at the middle -level of W component of the wind velocity: Wmin≈-1.9m/s; Wmax≈+1.7m/s; extrema above the mountain to the right: +1.2m/s, -1.9m/s,  +1.7m/s;  contour intervals are every 0.5m/s; б) height of the flow upper boundary: Ztop belongs to interval (20000m, 20020m); extrema to the right above the mountain: 20015m, 20000m, 20020m; contour intervals are every 5m.

 

Рис. 5. Поле W(X,Z) из [8] при обтекании холма через 5ч, изолинии через 0.05м/с, значения ΔX и ΔZ в работе не приведены, но судя по изображению профиля холма, они порядка 100м, , средняя скорость U потока  10м/с. В центрах модуль скорости ~0.4м/с.

Fig. 5. Vertical velocity profiles, W(X,Z), for flow over a  hill after 5 h from [9],  contour intervals are every 0,05·m/s, values of  ΔX and ΔZ  approximately 100m, , background velocity 10m/s.  Module of W in the centers ~0.4m/s.

На рис. 6 (в -системе координат) приведены изолинии вертикальной скорости при расчете на 5 ч по предложенному здесь алгоритму двумерного (X,)-обтекания (сетка 256×80 при разрешении 250 м по координате X и при почти таком же, слабо зависящем от X, разрешении по Z) этого же рельефа, но с принятыми в данной работе ‘чистыми’ граничными условиями: сверху (примерно на высоте 20 км) – свободная материальная поверхностью постоянного давления, по горизонтальной координате – периодичность. Параметры фонового потока взяты из [9]: начальная горизонтальная скорость 10 м/с, температурная стратификация  переменна и отвечает постоянной частоте Брента-Вяйсяля – 0,01с-1. Течение довольно быстро устанавливается, так что картина потока уже через 3 ч весьма близка к картине потока через 5 часов. 

 

 

Рис. 6. Поле W(X,) при обтекании холма через 5ч, изолинии через 0.10м/с, ΔX =250м ~ ΔZ, Δt ~0.7c, высота рисунка ~20км, длина 62.5км, начальная скорость U᷉ потока 10м/с. Wmax = 1.8м/с, Wmin = -1.85м/с. Max и Min величины достигаются на склонах холма. В центрах меандров максимум модуля W ~0.35м/с.

Fig. 6. Vertical velocity profiles, W(X,) for flow over a hill after 5h, contour
intervals are every
0.10m/s, values of
ΔX and ΔZ are approximately 250m,
Δt ~0.7s, initial U᷉(X,Z) = 10m/s. Module of W in the centers ~0.4m/s: Wmax=1.8m/s, Wmin = -1.85m/s. Max и Min values are at the hill boundary. At the centers max module of W ~0.35m/s.

 

Очевидно, что различие в постановке краевых условий может привести к заметному различию картины течений в описанных экспериментах. Так, количество центров на рис. 5 в нижнем 12 км слое такое же, как во всем ~20 км слое на рис. 6, При этом значения максимумов и минимумов вертикальной скорости в этих центрах довольно близки.

 

Заключение

Цель статьи – применить хорошо известные методы приближенного решения задач газовой динамики для аппроксимации уравнений идеальной атмосферы. Не менее важная цель – получить не уменьшающий порядок аппроксимации способ коррекции разностного решения с помощью аппроксимаций интегральных законов сохранения. Мотивация разработки состоит в максимальном использовании гиперболичности системы уравнений: сведении задачи к ‘переносу’ определенных комбинаций неизвестных функций вдоль характеристик в задачах Римана и в учете возможности пересечения характеристик одного семейства.

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

Удается преодолеть трудности, связанные с невозможностью провести классическое расщепление из-за ‘сильной’ нелинейности уравнений, возникающей при наличии рельефа.

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

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

Автор благодарен В.А. Гордину за указание использовать схему переменных направлений, Ю.В. Алферову за пакет программ расчерчивания изолиний, В.Д. Жупанову за советы по ОМР распараллеливанию, Ю.А. Степанову за представленную возможность проводить длительные расчеты на ЭВМ.

 

Список литературы

1. Сафронов А.В., Фомин Ю.В. Метод численного решения уравнений газодинамики с помощью соотношений на разрывах // Труды МФТИ. 2010. Том 2, № 2. С. 137-148.

2. Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука – Сибирское отделение, 1967. 197 с.

3. Годунов С.К., Рябенький В.С. Разностные схемы. М.: Наука, 1973. 400 с.

4. Уизем Дж. Линейные и нелинейные волны. М.: Мир, 1977. 624 с.

5. Roe P.L. Approximate Riemann Solvers, Parameter Vectors and Difference Scheme // J. Comput. Phys. 1981. Vol. 43. P. 357-372.

6. Peaceman D.W., Rachford H.H. The numerical solution of parabolic and elliptic differential equations // J. Soc. Ind. Appl. Math. 1955. Vol. 3. Р. 28-41.

7. Douglas J.Ir. Alternating direction methods for three space variables // Num. Math. 1962. Vol. 4, is. 1. P. 41-63.

8. Schär C., Leuenberger D., Fuhrer O., Lüthi D., Girard C. A new terrain-following vertical coordinate formulation for atmospheric prediction models // Mon. Wea. Rev. 2002. Vol. 130, no. 10. P. 2459-2480.

9. Melvin T., Dubal M., Wood N., Staniforth A., Zerroukat M. An inherently mass-conserving iterative semi-implicit semi-Lagrangian discretization of the non-hydrostatic vertical-slice equations // QJRMS. Apr 2010. Vol. 136 (648). P. 799-814. DOI:10.1002/qj.603.

 

References

1. Safronov A.V., Fomin Y.V. Metod chislennogo resheniya uravnenii gazodinamiki s pomosch'yu sootnoshenii na razryvah [Method of the numerical solution of a gas dynamics equation by correlations on ruptures]. Trudy MFTI [Proceedings of MIPT], 2010, vol 2, no. 2, pp. 137-148 [in Russ].

2. Yanenko N.N. Metod drobnyh shagov resheniya mnogomernyh zadach matematicheskoi fiziki. Novosibirsk, Nauka publ. Sibirskoe otdelenie, 1967, 197 p. [in Russ].

3. Godunov S.K., Ryaben'kii V.S. Raznostnye skhemy. Moscow, Nauka publ., 1973, 400 p. [in Russ].

4. Whitham G.B. Lineinye i nelineinye volny [Linear and Nonlinear Waves]. Moscow, Mir publ., 1977, 624 p. [in Russ]. Whitham G.B. Linear and Nonlinear Waves. New York, John Wiley & Sons publ., 1974, 635 p.

5. Roe P.L. Approximate Riemann Solvers, Parameter Vectors and Difference Scheme. J. Comput. Phys., 1981, vol. 43, pp. 357-372.

6. Peaceman D.W., Rachford H.H. The numerical solution of parabolic and elliptic differential equations. J. Soc. Ind. Appl. Math., 1955, vol. 3, pp. 28-41.

7. Douglas J. Ir. Alternating direction methods for three space variables. Num. Math., 1962, vol. 4, no. 1, pp. 41-63.

8. Schär C., Leuenberger D., Fuhrer O., Lüthi D., Girard C. A new terrain-following vertical coordinate formulation for atmospheric prediction models. Mon. Wea. Rev., 2002, vol. 130, no. 10, pp. 2459-2480.

9. Melvin T., Dubal M., Wood N., Staniforth A., Zerroukat M. An inherently mass-conserving iterative semi-implicit semi-Lagrangian discretization of the non-hydrostatic vertical-slice equations. QJRMS, 2010, vol. 136, no. 648, pp. 799-814. DOI:10.1002/qj.603.

 

Поступила в редакцию 08.07.2019 г.

Received by the editor 08.07.2019.