Уравнение адвекции и численные методы его решения

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

, (81)

где U, V, W представляют собой компоненты скорости вдоль трехмерной декартовой системы координат. Это уравнение имеет важное самостоятельное значение, поскольку необходимость в его решении появляется на одном из этапов метода расщепления по физическим процессам. Качество общего решения задачи адвекции - диффузии - реакции сильно зависит от того, каким образом аппроксимированы члены, описывающие адвективный перенос. Ошибки аппроксимации, проявляющиеся в виде паразитной дисперсии и аппроксимационной вязкости, приводят не только к количественному, но и к качественному искажению получаемого решения. Ввиду большой практической значимости, проблема поиска новых направлений кардинального улучшения диссипативных и дисперсионных свойств численных схем для уравнения переноса является чрезвычайно актуальной, а соответствующая задача может быть отнесена к разряду фундаментальных задач вычислительной гидродинамики. Именно важностью приложений объясняется большое внимание исследователей к решению подобного рода уравнений. В последние годы интенсивное изучение проблем механики сплошной среды, прогноза погоды, распространения загрязняющий веществ, динамики океана привело к разработке достаточно универсальных и эффективных методов численного уравнения адвекции (Вабищевич и др., 2000; Марчук, 1974; Марчук, 1982; Рождественский, Яненко, 1978; Роуч, 1980; Самарский, Вабищевич, 1998; Arakawa, 1997; Bryan, 1997). Наряду с этим продолжается и поиск принципиально новых подходов к решению этой задачи (Головизнин, Карабасов, 1998; Головизнин, Самарский, 1998; Головизнин, Самарский, 1998). В данном разделе мы рассмотрим некоторые из таких подходов.

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

, (82)

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

, (83)

где - амплитуда в начальный момент времени, - волновое число, L - длина волны, - мнимая единица. Найдем решение уравнения (82) методом разделения переменных. Представим в виде произведения некоторых двух функций:

. (84)

После подстановки выражения (84) в уравнение (82) и деления результата на получим

Выполнив интегрирование, получим:

или , , где , - постоянные интегрирования. Следовательно, , . Искомая функция -

Используя начальное условие (83), имеем: , , а следовательно, окончательное решение можно представить в следующем виде:

(85)

при условии, что - есть дифференцируемая функция. Это решение описывает процесс начального возмущения вдоль характеристики . Это означает, что на любой прямой . Итак, задача (82) при определяет процесс распространения возмущения в сторону возрастающих значений . Это, хорошо известное положение учитывается и при конструировании разностных схем для уравнения адвекции. Среди наиболее употребительных на практике схем для численного решения одномерного уравнения переноса можно отметить следующие схемы (Белов и др., 1989; Мезингер, Аракава, 1979; Роуч, 1980):

·согласованную, трехуровневую по времени, условно устойчивую схему центральных разностей второго порядка точности. Условие устойчивости для данной схемы соответствует критерию Куранта - Фридрихса - Леви (КФЛ): ;

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

·неявную двухуровенную схему трапеций второго порядка (схема Кранка-Николсона). Эта схема нейтральна и абсолютно устойчива при любых шагах по пространству и времени. Кроме того она свободна от эффекта вычислительной вязкости;

·согласованную двухтактную схему Мацуно (Эйлера с пересчетом), первого порядка точности по времени и второго - по пространству. Эта схема устойчива при выполнении критерия КФЛ. Ей свойственна вычислительная вязкость, причем вязкость максимально воздействует на волны длиной . На волны длиной счетная вязкость не влияет;

·согласованную явную схему четвертого порядка точности по пространству и второго - по времени. Схема условно устойчивая, причем условие устойчивости оказывается более жестким, чем для схем первого и второго порядка по пространству. Для схем четвертого порядка по пространству шаг по времени должен быть примерно на 27% меньше, чем для схем второго порядка (при одних и тех же U и ). В случае применения этой схемы эффект вычислительной вязкости не проявляется;

·явную согласованную схему Лакса-Вендрофа второго порядка точности. Важным достоинством схемы является ее способность подавлять волны длиной , которые нельзя правильно представить на сетке. Тем самым появляется возможность демпфирования коротких волн в процессе численного интегрирования уравнений. Схема условно устойчива. Условие устойчивости соответствует критерию КФЛ. При схема Лакса-Вендрофа диссипативна. По мере увеличения длин волн диссипативность схемы уменьшается.

В работах (Головизнин, Карабасов, 1998; Головизнин, Самарский, 1998; Головизнин, Самарский, 1998) предложена принципиально новая базисная линейная разностная схема для одномерного модельного уравнения конвективного переноса с постоянной скоростью, названная схемой "кабаре". Эта схема имеет ряд уникальных особенностей, в связи с чем ее дальнейшее изучение представляет как теоретический, так и прикладной интерес. Главной особенностью схемы "кабаре" является то, что она оказывается точной при двух различных числах Куранта, лежащих в области ее устойчивости. Это означает, что в этих случаях данная схема не содержит также и аппроксимационной дисперсии, что обеспечивает наличие у нее достаточно хороших транспортных свойств. Так, на нерегулярной расчетной сетке схема без существенного искажения переносит нетривиальные начальные распределения на расстояния, на порядок больше, чем это позволяют классические линейные схемы, такие, например, как схема Лакса-Вендрофа. Другая особенность схемы заключается в том, что она имеет второй порядок аппроксимации на нерегулярных расчетных сетках и исключает распространение возмущений вверх по потоку. Кроме того, схема "кабаре" является консервативной (дивергентной), и для нее автоматически выполняется закон сохранения квадрата переносимой функции, который является аналогом закона сохранения энергии. Поэтому можно говорить о том, что эта схема обладает свойством полной консервативности. Путем тестовых расчетов было показано, что при наличии значительных пространственных вариаций у переносимого профиля решение теряет свою монотонность. Поэтому в схему "кабаре", как и в другие малодиссипативные линейные схемы, необходимо включать алгоритмы нелинейной коррекции, монотонизирующие получаемые решения (Головизнин, Карабасов, 1998). Есть все основания полагать, что распространение данной схемы на многомерный случай может также привести к получению качественно новых результатов.

Если скорость - переменная по пространству и во времени, то нахождение задачи (82) в аналитическом виде, в общем случае, невозможно. Именно в таких ситуациях и оказывается необходимым применение численных методов, основанных на разностных аппроксимациях. Следует только иметь виду, что даже при использовании неявных диссипативных разностных схем возможно нарушение счетной устойчивости (Белов и др., 1989). Особенно это проявляется в нелинейных задачах. Это объясняется следующим (Мезингер, Аракава, 1979). Если разложить решение разностной задачи и коэффициент в ряд Фурье, то произведения рядов Фурье приведут к гармоникам как более длинным, чем взаимодействующие, так и более коротким. В результате такого процесса в ряде случаев может произойти перекачка "энергии" от ошибок округления из длинных волн в более короткие, и процесс вычислений окажется неустойчивым, даже несмотря на то, что данная разностная схема с постоянным коэффициентом будет счетно устойчивой. Она также иногда появляется и при решении линейных задач с переменными коэффициентами. Поэтому построение разностных схем для нелинейных уравнений с переменными коэффициентами, устойчивых в отношении к любым возмущениям, является чрезвычайно актуальной задачей. В большинстве случаев подавление счетной неустойчивости возможно с помощью диссипативных разностных схем, отвечающих определенному выбору коэффициента счетной вязкости. Однако такие схемы, как правило, оказываются схемами первого порядка аппроксимации либо по , либо по , либо и по , и по .

Рассмотрим теперь задачу адвекции на плоскости :

в , , (86)

на , (87)

в при t = 0. (88)

Здесь , . Будем предполагать, что компоненты вектора скорости U, V в каждый момент времени удовлетворяют уравнению неразрывности

. (89)

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

в , , (90)

на , (91)

в при t = 0. (92)

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

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

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

. (93)

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

(94)

не удовлетворяет условию (93). В этом случае мы получаем следующие соотношения:

, .

Таким образом, операторы в форме (94) не могут быть взяты в качестве элементарных для построения схемы последовательного расщепления.

Выберем теперь более сложную форму представления операторов :

. (95)

Каждый из таких операторов теперь удовлетворяет условию (93), а сумма этих операторов в точности равна . В самом деле, составим сумму

Здесь мы воспользовались уравнением неразрывности (89).

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

(96)

Как было отмечено выше, если , то при достаточной гладкости решения задачи и элементов матриц система разностных уравнений (96) абсолютно устойчива и эта схема аппроксимирует исходную задачу со вторым порядком точности по . Кроме того, данный прием снимает весьма сильное требование коммутативности операторов. Это - важный и поучительный пример того, насколько существенно учитывать математические особенности решаемой задачи и сущность описываемых процессов. Формальное, прямолинейное применение методов расщепления может скомпрометировать саму идею.

Выписывая схему (96) мы еще ничего не сказали о путях рациональной аппроксимации оператора конвективного переноса по пространственным переменным . Удобным методом аппроксимации с сохранением аддитивных свойств оператора и его основных качественных особенностей является метод покоординатной аппроксимации (Марчук, 1980). Предположим, что поле скоростей достаточно гладкое и рассмотрим уравнение адвекции в дивергентной (консервативной) форме:

в , , (97)

на , (98)

в при t = 0. (99)

За основу построения разностной схемы возьмем оператор конвективного переноса С, который определяется выражением

, (100)

а разностный аналог этого соотношения рассмотрим в виде

(101)

Разностное выражение (101) аппроксимирует исходный оператор (100) со вторым порядком относительно пространственных шагов и сетки при достаточной гладкости решения и поля скоростей . Однако выражение (101) имеет существенный недостаток, поскольку в такой форме оператор нарушает свою кососимметрическую структуру, то есть теперь . А это означает, что внешне привычная аппроксимация оказывается неудовлетворительной для построения вычислительного алгоритма решения поставленной задачи. Для того, чтобы соотношение выполнялось, необходимо осуществить аппроксимацию оператора (100) в следующем виде (Марчук, Саркисян, 1988):

, (102)

где , , , . Выражение (102), так же как и (101), аппроксимирует исходный оператор (100) со вторым порядком относительно пространственных шагов сетки и . Таким образом построенный оператор сохраняет важное свойство кососимметричности и, более того, каждый из операторов и , определяемых выражениями

, (103.1)

(103.2)

также удовлетворяет условиям , (Марчук, Саркисян, 1988).

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

, (104)

при t = 0, (105)

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

(106)

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

Рассмотрим теперь другой вариант построения конечно-разностной схемы для двумерного уравнения адвекции (86) - (88) на интервале времени, равном шагу по времени от до - схему предиктор-корректор (Белов и др., 1989; Моделирование процессов переноса и трансформации вещества в море, 1979). Здесь уравнение (86) расщепляют на два уравнения так, чтобы в каждом из них присутствовали производные только по одной пространственной переменной:

, . (107)

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

(108)

(109)

Эти конечно-разностные схемы применяются последовательно так, что решение, полученное на предыдущем дробном шаге, используется в качестве начального условия для следующего дробного шага. В совокупности схемы (108), (109) обеспечивают получение решения на момент времени . Поскольку в этих схемах применяются направленные разности по времени и по пространству, то аппроксимация ими уравнения (86) имеет первый порядок точности. Схемы (108), (109) называются схемами предиктора.

Далее, на интервале времени от до , уравнение (86) решается с помощью схемы центральных разностей

, (110)

которая аппроксимирует уравнение (86) уже со вторым порядком точности. Эта схема называется схемой корректора.

Рассмотрим более подробно решение системы уравнений (107). Поскольку оба уравнения этой системы однотипны, рассмотрим решение только первого уравнения

. (111)

Для простоты разобьем шаг по времени на два равных интервала: от до и от до . На первом интервале (дробном шаге) от до используется схема предиктора

(112)

Двухточечные схемы (112) формально можно записать в виде единой трехточечной схемы, объединяющей случаи и :

, (113)

где , , , . Схема (113), записанная для точек i = 1, 2, ..., I-1, представляет собой систему алгебраических уравнений, в которых содержатся неизвестные функции в трех соседних точках: , , . Матрица коэффициентов этой системы уравнений - трехдиагональная. В точках i = 0, i = I ставятся соответствующие задаче граничные условия (например, , ). Система уравнений (113) с граничными условиями является замкнутой. Для решения этой системы применяется метод разностной факторизации (прогонки). Запишем систему (113) в следующем виде:

(114)

Для того, чтобы обеспечить точное соблюдение граничных условий, полагаем , , , , , . Выразим из первого уравнения . Получим:

, (115)

где , . Подставляя (115) во второе уравнение, найдем :

, (116)

где , . Повторяя далее процесс подстановки, для некоторого промежуточного значения i получим:

, (117)

где , . Для i = I-1 из (117) получим :

, (118)

где , . Наконец, подставляя (118) в последнее уравнение системы (114), получаем:

. (119)

Поскольку , , , то из (119) следует, что граничное условие при в точности выполняется. Таким образом, видно, что метод прогонки реализуется с помощью двух циклов вычислений, которые называются прогонками вперед и назад. На прогонке вперед (индекс i возрастает, пробегая значения i = 1, 2, ..., I-1) рассчитываются вспомогательные величины , , , . На прогонке назад (индекс i убывает, пробегая значения i = I-1, I-2, ..., 1) по формуле определяются искомые величины, причем .

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

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

Покажем, что аппроксимация уравнения (111) конечно-разностным уравнением (112) приводит к появлению эффекта счетной вязкости (Белов и др., 1989). Для этого разложим решение уравнения (111) в ряды Тейлора в окрестности точки , , :

Подставляя данные разложения в конечно-разностное уравнение (113), получаем:

, (120)

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

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

. (121)

Таким образом, применение схем предиктора и корректора обеспечивает получение решения в момент времени . Аппроксимация уравнения (111) схемами предиктора и корректора имеет второй порядок точности.

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

. (122)

Из конечно-разностных схем (121) и (122) исключим неизвестную величину с дробным временным индексом . Для этого запишем схему (122) для точек и :

Вычитая из первого соотношения второе, получим:

. (123)

Схему (121) запишем для точки : . Отсюда находим

. (124)

Используя выражение (124), уравнение (123) можно представить в следующем виде:

. (125)

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

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

Для приближенного решения начально-краевых задач для уравнения переноса широко используются также разностные схемы, для которых выполнен принцип максимума - монотонные разностные схемы. Для уравнений переноса легко строятся безусловно монотонные разностные схемы первого порядка по пространству с использованием направленных разностных производных (Самарский, 1983). Хорошо известно (Годунов, 1959), что в классе однородных линейных схем отсутствуют схемы с порядком аппроксимации, по пространству выше первого. При численном решении задач механики сплошной среды строятся различные классы нелинейных монотонных схем, которые в основной части расчетной области имеют повышенный порядок - второй и выше. В настоящее время большое внимание в прикладных исследованиях по механике сплошной среды уделяется разностных схемам типа TVD(Total Variation Diminishing). Безусловная монотонность в этих и им подобных схемах достигается за счет нелинейного ограничения потоков в исходной схеме с направленными разностями второго порядка. Эти исследования были инициированы в работах (Колган, 1972; Harten, 1983) и продолжены во многих других (Hartenetal., 1986; Harten, Osher,1987; Hartenetal., 1987). В работе (Самарский, Вабищевич, 1998) отмечается возможность построения монотонных однородных нелинейных разностных схем для уравнения переноса на основе регуляризации (возмущения) схемы второго порядка аппроксимации. В качестве производящей схемы выбираются схемы со стандартными аппроксимациями первой производной второго порядка - центральными разностями и направленными разностями второго порядка, а также их комбинациями.

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

. (126)

Уравнение (126) дополняется начальным условием

. (127)

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

, (128)

где - нормаль к границе расчетной области. При ограничениях (128) нет необходимости формулировать какие-либо граничные условия для однозначного определения решения задачи (126), (127). Для поставленной задачи справедлив принцип максимума: , если . Кроме того, верна априорная оценка в : , где . Естественно ориентироваться на построение таких разностных схем, для которых эти важнейшие свойства были бы выполнены.

Рассмотрим теперь разностную схему с направленными разностями. После дискретизации по пространственным переменным от (126) придем к дифференциально-операторному уравнению

. (129)

В соответствии с (127) положим

. (130)

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

(131)

При построении схем с направленными разностями учитывается знак скорости в соответствующем направлении. Будем использовать следующие обозначения:

С учетом (128) положим ,

(132)

Разностная схема (131), (132) устойчива в при следующих ограничениях (условие Куранта) на шаг по времени:

. (133)

При этом для решения имеет место оценка , где . Разностная схема (131) с центрально-разностными аппроксимациями конвективных слагаемых

(134)

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

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

(135)

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

(136)

В этом случае значение соответствует использованию центрально-разностных аппроксимаций (134), - аппроксимаций направленными разностями второго порядка. Среди аппроксимаций второго порядка (136) отдельного упоминания заслуживает случай и - аппроксимации третьего порядка.

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

В нелинейной регуляризованной схеме используются аппроксимации конвективного переноса, которые имеют форму аппроксимации направленными разностями со специальными регуляризующими множителями:

. (137)

Принимая во внимание (136), для регуляризующих множителей в работе (Вабищевич и др., 2000) получено следующее представление

, (138.1)

. (138.2)

Здесь и этот коэффициент выступает в качестве параметра регуляризации. В предельном случае мы приходим к нелинейной аппроксимации (136) второго или третьего (при ) порядка. Монотонность схемы (131), (137), (138.1) - (138.2) (неотрицательность регуляризующих коэффициентов ) обеспечивается выбором коэффициента . Принимая во внимание то, что для , приходим к следующим ограничениям на временной шаг: , где (см. (133)) - минимальный шаг для схемы с направленными разностями первого порядка. Тем самым при минимально допустимом значении параметра предельный шаг по времени в нелинейной монотонной схеме в два раза меньше по сравнению с обычной схемой (131), (132). Отметим, что нет простых соображений в пользу того или иного однозначного выбора параметра . На тестовых расчетах для одномерной задачи показано, что параметр влияет на перемещение профиля с большей или меньшей скоростью по сравнению с точным решением тестовой задачи. Поэтому его связывают с дисперсионными характеристиками разностной схемы. Установлено, что наиболее оптимальным является значение . Важным достоинством рассматриваемых схем является их простота и прозрачность базовых математических конструкций.