Журнал вычислительной математики и математической физики, 2019, T. 59, № 9, стр. 1570-1580

ВВЕДЕНИЕ

В сложных прикладных задачах часто встречаются жесткие системы дифференциальных уравнений, в которых характерные скорости изменения отдельных компонент очень сильно отличаются друг от друга. Для численного решения таких задач необходимо применять методы с повышенной степенью устойчивости. В работах [1], [2] предложены и рассмотрены мульти-неявные методы со второй производной (mISD, m-Implicit with Second Derivative). В классических многошаговых методах используется информация в нескольких предшествующих точках. В противоположность этому в мульти-неявном методе используют несколько последующих точек, для каждой из которых записывается отдельное разностное уравнение. Совокупность этих разностных уравнений рассматривается как единый разностный метод. Итерационный метод Ньютона, обычно используемый для разрешения неявных схем, требует вычисления матриц производных правой части, которые непосредственно связаны со вторыми производными решения. Использование этой информации в конструкции самой разностной схемы позволяет существенно повысить порядок аппроксимации метода без больших дополнительных вычислительных затрат. Построены и протестированы симметричные А-устойчивые методы 6-го, 8-го и 10-го порядка аппроксимации. В работах [3], [4] построены и исследованы расширенные семейства 2ISD и 3ISD методов, среди которых обнаружены L-устойчивые методы. Тестирование на простых примерах показало, что эти методы обладают повышенной степенью устойчивости и предпочтительны для расчетов неравновесных процессов. Однако для успешного решения сложных нелинейных задач одной L-устойчивости может быть недостаточно, что связано с недостатками метода Ньютона из-за его чувствительности к начальному приближению. На локальных участках решения при относительно большом шаге интегрирования метод Ньютона может не сходиться или сходиться к побочному нефизическому решению.

Устранить этот недостаток позволяют специальные алгоритмы автоматического контроля точности и изменения шага интегрирования в процессе решения задачи. Такие алгоритмы известны как для явных (см. [5]–[8]), так и для неявных (см. [9], [10]) методов Рунге–Кутты. Идея состоит в использовании так называемых встроенных методов, т.е. двух методов с одинаковыми стадиями, но разного порядка аппроксимации из-за различий в заключительной строке коэффициентов таблицы Бутчера. Расчет проводится по методу более высокого порядка, а контроль точности и шага интегрирования осуществляется по формулам метода меньшего порядка.

В данной работе рассмотрен аналогичный подход на базе mISD-схем высокой точности. Описан алгоритм автоматического выбора шага интегрирования для 3ISD-схемы 8-го порядка с контролем точности по 2ISD-схеме 6-го порядка. Продемонстрирована эффективность алгоритма выбора шага на примере задачи неравновесной кинетики химических реакций водородно-кислородной смеси, включающей 6 компонент (H₂, O₂, H, O, OH, H₂O).

1. МУЛЬТИ-НЕЯВНЫЕ МЕТОДЫ СО ВТОРОЙ ПРОИЗВОДНОЙ

Рассматриваем задачу Коши для системы обыкновенных дифференциальных уравнений

где u(t) и f(u) являются векторными функциями размерности p. Введем равномерную сетку {t_n = τn, n = 0, 1, 2…} с постоянным шагом τ > 0. Функции, определенные на сетке и приближенные решения обозначим как v_n = v(t_n), f_n = f(v_n).

Основная идея мульти-неявных методов со второй производной для задачи (1.1) заключается в совместном использовании нескольких разностных уравнений, включающих в неявном виде несколько искомых значений численного решения, а также функций и производных правой части в моменты времени t_n, t_n + 1, … t_n + m:

где производные в правой части выражаются через матрицу Якоби J:

При известном v_n = v(t_n) разностная схема (1.2) представляет собой систему взаимосвязанных нелинейных разностных уравнений относительно неизвестных v_n + 1, v_n + 2, …, v_n + m в m последующих точках t_n + 1, t_n + 2, …, t_n + m. Такую схему будем называть мульти-неявной.

Для нахождения коэффициентов схемы (1.2) запишем ошибку аппроксимации, используя свойство решений системы (1.1) f = u', f' = u'':

Разложения в (1.3) по степеням τ с центром в точке t_n позволяют выписать систему линейных уравнений на коэффициенты, которые обеспечивают для всех k равенство Полученные коэффициенты удобно записывать в матричном виде. Для m = 1 получим одношаговую разностную схему Энрайта 4-го порядка (см. [11]) Для m = 2 получим дважды неявную схему (2ISD) 6-го порядка из работы [1] Для m = 3 получим трижды неявную схему (3ISD) 8-го порядка, которая рассмотрена в работе [2]

В работе [3] приведены коэффициенты для схемы 10 -го порядка аппроксимации. Там же показано, что все mISD-схемы максимального порядка (т.е. порядка 2m + 2) являются А-устойчивыми, причем область устойчивости в точности совпадает с левой полуплоскостью комплексной плоскости. Это свойство следует из свойств центральной симметрии для коэффициентов a_ki и антисимметрии для коэффициентов b_ki, а именно

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

Нахождение разностного решения m-неявной схемы из системы (1.2) осуществляется итерационным методом Ньютона, на каждой итерации которого нужно решать систему линейных уравнений размерности в m раз больше, чем размерность исходной системы (1.1). Вычислительная сложность, т.е. количество операций, при этом увеличивается в m³ раз. Однако в прикладных задачах основную трудоемкость составляют вычисления матриц Якоби J системы (1.1), количество которых в m-неявной схеме пропорционально лишь первой степени m.

2. РАСШИРЕННЫЕ СЕМЕЙСТВА МУЛЬТИ-НЕЯВНЫХ МЕТОДОВ

Способ расширения мульти-неявных методов продемонстрируем на примере трижды неявного метода 8-го порядка (1.4). Запишем систему разностных уравнений (1.2) в виде, в котором левая часть каждого уравнения содержит лишь одно неизвестное значение

Коэффициенты в правой части системы (2.1) будем искать из условия, чтобы лишь последнее уравнение обеспечивало максимальный 8-й порядок аппроксимации, а порядок аппроксимации первых двух снизим на 1. Порядок аппроксимации последнего разностного уравнения в (2.1) не изменится, если значения v_n + 1 и v_n + 2, используемые в его правой части, будут вычисляться лишь с 7-м порядком. Этот подход позволяет ввести в первые два уравнения два свободных параметра с сохранением 8-го порядка точности для v_n + 3 в последнем уравнении. Матрицы коэффициентов схемы 8 -го порядка для (2.1) с двумя свободными параметрами α и β выглядят следующим образом:

В работе [3] проведено исследование такого двухпараметрического семейства трижды неявных схем. В частности показано, что при $\alpha = 1{\text{/}}54$ и $\beta = - 1{\text{/}}216$ получается L₂-устойчивая схема 8 -го порядка аппроксимации, а при $\alpha = 1{\text{/}}540$ и $\beta = 1{\text{/}}1080$ получается А-устойчивая схема, которая на линейных задачах достигает 10-го порядка аппроксимации. Теоретические выводы в [3] подтверждены серией численных тестов.

3. ВЫБОР ШАГА ИНТЕГРИРОВАНИЯ С ЗАДАННОЙ ПОГРЕШНОСТЬЮ АППРОКСИМАЦИИ

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

Опишем процедуру автоматического контроля шага интегрирования применительно к методу 8-го порядка (2.1). Пусть v = {v_n, v_n + 1, v_n + 2, v_n + 3} – решение этой схемы L₈(τ, v_n, v_n + 1, v_n + 2, v_n + 3) = 0 на шаге от t_n до t_n + 3τ. Контроль шага τ будем осуществлять по погрешности аппроксимации второго разностного уравнения схемы 6 -го порядка. Это уравнение имеет симметричную форму

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

будет включать в себя величину шага τ как дополнительное неизвестное. Можно решать совместную систему (3.2) итерационным методом Ньютона с начальным значением τ с предыдущего шага. Но учитывая, что контролирующая разностная схема симметрична и на два порядка отличается от расчетной, можно построить экономичный способ определения τ на основе простой итерации. Различие в два порядка позволяет рассматривать решение более точной схемы в качестве точного решения дифференциального уравнения по отношению к менее точной схеме. С другой стороны, ошибка аппроксимации контролирующей симметричной схемы (3.1) на решении дифференциального уравнения содержит только четные степени по τ, т.е. Тогда итерационный процесс определения τ из условия при малых τ сходится со скоростью, близкой к квадратичной. Итерационное уточнение шага происходит следующим образом. Пусть есть предварительное значение τ_p. С этим шагом выполняется схема 8 -го порядка. На ее решении вычисляется ошибка аппроксимации схемы 6 -го порядка и уточненное значение τ из условия (3.4):

Для целенаправленного контроля точности желательно связать задаваемую ошибку аппроксимации δ с погрешностью численного решения. Для явных схем произведение ошибки аппроксимации на шаг τ представляет собой норму погрешности численного решения на одном шаге интегрирования при старте с точного решения. Следовательно, результат подстановки численного решения 8-го порядка в схему 6 -го порядка характеризует локальную погрешность контролирующей схемы 6 -го порядка. Для оценки желаемой погрешности ε численного решения расчетной схемы в конце участка интегрирования Δ можно использовать сумму норм локальных погрешностей контролирующей схемы, что приводит к связи

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

Таким образом, алгоритм автоматического контроля и изменения шага интегрирования состоит из следующих этапов.

1. От текущего v_n с использованием величины предыдущего шага τ_p находим v_n + 1, v_n + 2, v_n + 3 по схеме 8 -го порядка (3ISD).

2. Подстановкой этого решения в схему 6 -го порядка (2ISD) вычисляем ее ошибку аппроксимации S и уточненное значение шага τ по формуле (3.5).

3. Если τ отличается от τ_p более, чем на 1%, то повторяем шаг интегрирования с новым значением τ, в противном случае переходим к следующему шагу.

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

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

4. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ КИНЕТИКИ ВЗРЫВА ГРЕМУЧЕЙ СМЕСИ

Для проверки работоспособности выше описанных алгоритмов рассмотрим задачу поведения взрывчатой газовой смеси (H₂ + O₂) в замкнутом объеме с подвижным поршнем. С помощью поршня осуществляется цикл сжатия-расширения газовой смеси по заданному закону. При начальных условиях реакции в смеси не происходят. На этапе сжатия повышается температура смеси, в результате чего возбуждается механизм химических реакций с выделением тепла. Задача относится к классу жестких, так как скорости химических реакций в процессе сжатия-расширения изменяются на несколько порядков.

Масса смеси в замкнутом объеме постоянна, поэтому сжатие-расширение по заданному закону будем описывать заданным поведением плотности смеси во времени ρ(t). Для описания реакций использовалась модель химической кинетики из работы [12]. Смесь состоит из 6 компонент (H₂, O₂, OH, H₂O, H, O), каждая из которых является совершенным газом с постоянными теплоемкостями. В табл. 1 приведены энтальпии образования H_i, молекулярные веса κ_i, и отношения теплоемкостей γ_i, участвующих в реакциях газообразных компонент.

Химическая кинетика состоит из 6 реакций, включающих в себя один из основных цепных механизмов окисления, при этом в реакциях диссоциации учитывается различная способность реагентов выступать в роли дополнительной частицы. Формулы и константы скоростей прямых и обратных реакций $c_{r}^{ + }$, $c_{r}^{ - }$, r = 1, 2, …, 16 приведены в табл. 2. Считаем, что скорости реакций подчиняются обобщенной аррениусовской зависимости от температуры T в виде с = AT ^mexp(–E/RT). Газовая постоянная R = 1.98647 кал/(моль K).

Реагирующая среда рассматривается как однородная по пространству смесь из шести идеальных совершенных газов. Мы пренебрегаем скоростью движения газа внутри объема, что оправдано в случае, если скорость поршня гораздо меньше скорости звука. Каждая компонента имеет собственную плотность ρ_i. Температура всех компонент считается одинаковой и равна температуре смеси T. Таким образом, состояние смеси в каждый момент времени характеризуют 7 параметров: молярные концентрации α₁, α₂, …, α₆, и температура T, где ${{\alpha }_{i}} = {{\rho }_{i}}{\text{/}}({\rho }{{\kappa }_{i}})$. Система дифференциальных уравнений для этих величин будет иметь вид

Молекулярный вес κ и отношение теплоемкостей γ смеси, присутствующие в последнем уравнении, связаны с молярными концентрациями α_i известными соотношениями

Стехиометрические коэффициенты $\nu _{{ir}}^{ + }$ и $\nu _{{ir}}^{ - }$ обозначают количество атомов или молекул i-го компонента, вступающих в r-ю прямую или обратную реакцию. Полученная система уравнений описывает поведение смеси (H₂ + O₂) при заданном изменении плотности смеси ρ(t).

5. РЕЗУЛЬТАТЫ РАСЧЕТОВ

Начальные условия смеси удобнее задавать числовыми долями компонент n_i, начальным давлением p₀ и температурой T₀. Для смеси стехиометрического состава n_i = (2, 1, 0, 0, 0, 0). Через них находим молярные концентрации и плотность смеси:

Рассматриваем процесс сжатия-расширения, график плотности смеси для которого изображен на фиг. 1a. На интервале [0, t_a] плотность увеличивается от ρ₀ до ρ_max. Далее на интервале [t_a, t_b] уменьшается от ρ_max до ρ_min, после чего на интервале [t_b, t_end] плотность постоянна. Для задания переходных процессов используем функцию с непрерывными старшими производными на всем промежутке времени:

здесь θ(x) – гладкая монотонная переходная функция [0, 1] → [0, 1] с нулевыми производными на краях. Алгоритм расчета шага интегрирования использует свойства ошибки аппроксимации, определяемой на гладком решении. Чтобы обеспечить гладкость функции ρ(t) в данной работе использовалась Такой вид переходной функции обеспечивает непрерывность всех производных для ρ(t) в точках стыка t = t_a и t = t_b.

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

На фиг. 1б изображены зависимости концентраций компонент с течением времени. Отчетливо видна зона задержки, в течение которой происходит накопление отсутствующих компонент, необходимых для возбуждения механизма цепной реакции. В районе t = 5 мкс наблюдается взрывообразное изменение концентраций за время менее 1 мкс. Очевидно, что из соображений точности расчет на промежутке взрыва нужно проводить с малым шагом по времени. Значительные изменения концентраций наблюдаются также на завершающей стадии расширения смеси, когда активизируются реакции рекомбинации и образования H₂O. Этот промежуток также требует малого шага интегрирования. Между этими участками, а также после момента остановки поршня (t > t_b) концентрации изменяются медленно, и шаг по времени может быть увеличен.

Фигура 2 иллюстрирует изменение температуры смеси. На периоде индукции температура немного растет из-за сжатия смеси. После достижения 1300 K рост резко усиливается и за очень короткий интервал температура превосходит 3500 K. Взрывной рост температуры объясняется большим и быстрым выделением тепла во время основной фазы реакций. На стадии расширения температура снижается, несмотря на то, что здесь также происходят экзотермические реакции. После остановки поршня температура изменяется незначительно, смесь на этой стадии медленно приходит в равновесное состояние.

Приведенные результаты были получены методами 8-го и 6-го порядков с использованием процедуры автоматического выбора шага интегрирования, изложенной в разд. 3, с заданной постоянной ε = 10^–8.

Опишем некоторые детали реализации применяемых методов. Компоненты вектора решения $x = \{ {{\alpha }_{1}},{{\alpha }_{2}},\; \ldots ,\;{{\alpha }_{6}},T\} $ имеют разные единицы измерения. Компоненты вектора относительного изменения δx вычислялись через абсолютные Δx как

При вычислении локальной погрешности в формуле (3.5) использовалась евклидова норма вектора δx. Эта же норма использовалась при оценке точности численного решения, а также в критерии остановки метода Ньютона на шаге интегрирования. Итерации в методе Ньютона завершались, когда отношение нормы добавки в очередной итерации к норме суммы добавок в предыдущих итерациях становилось менее 10^–11. В качестве начальных приближений для v_n + 1, v_n + 2, v_n + 3 в методе Ньютона использовалось v_n, т.е. какая-либо экстраполяция для начальных приближений не применялась. Матрица производных для метода Ньютона вычислялась по аналитическим формулам, однако производные от якобианов, входящих в разностную схему, не учитывались. Первые интегралы системы уравнений (4.1), связанные с сохранением количества атомов водорода и кислорода, в расчетах не использовались.

На фиг. 3 изображены зависимости шага интегрирования от времени в процессе численного решения с заданной погрешностью ε = 10^–8 двумя методами. Кривая 1 относится к расчету методом 8-го порядка с контролем шага по формулам 6-го порядка (пара “8-6”), кривая 2 изображает аналогичную зависимость для пары “6-4” (расчет по методу 6-го порядка с контролем по методу 4-го порядка). В обоих случаях для расчета и контроля использованы А-устойчивые схемы. Здесь под шагом τ подразумевается полный шаг мульти-неявной схемы, выраженный в микросекундах.

Для расчета с использованием пары “8-6” потребовалось 129 шагов по времени при 5 итерациях метода Ньютона в среднем на каждый шаг. Для пары “6-4” при той же точности потребовалось 1408 шагов по времени при 4 итерациях метода Ньютона в среднем на каждый шаг. Несмотря на то что метод 8-го порядка в среднем потребовал приблизительно в 2.5 раза больше операций на 1 шаг интегрирования по сравнению с методом 6-го порядка, его суммарные затраты машинного времени оказались в 4 с лишним раза меньше, чем для метода 6-го порядка. Наибольший вклад в различие по количеству шагов двух методов дает зона активных реакций, которая отмечена буквой R. Пара меньшего порядка “6-4” в этой зоне потребовала величину шага τ в ∼20 раз меньше, чем пара “8-6”. Также видно, что уменьшение шага интегрирования наблюдается в зонах сопряжения функции ρ(t) при t = 15 мкс и t = 30 мкс, что связано с большой величиной и переменой знака экстремумов старших производных у переходной функции (5.1). Этим же объясняются локальные всплески шага интегрирования для пары “6-4” (кривая 2). Для пары “8-6” они проявляются в меньшей степени в виде немонотонного поведения τ на заключительном этапе расширения (20 < t < 30).

Реальная оценка погрешности численного решения осуществлялась путем сравнения с наиболее точным численным решением. В качестве такого эталона использовалось решение L₂-устойчивого метода 8-го порядка с постоянным шагом τ = 2 × 10^–4 мкс, то есть на равномерной сетке с 225 000 узлами. Коэффициенты метода приведены в матрицах формулы (2.2). Кривая 1 на фиг. 4 представляет распределение погрешности (т.е. нормы разности с эталонным решением) численного решения по схеме 8 -го порядка с переменным шагом, соответствующим кривой 1 на фиг. 3. Жирные точки на кривой соответствуют узлам расчетной сетки. Видно, что прогнозируемая точность ε = 10⁻⁸ не выдерживается лишь в зоне активных реакций. Кривая 2 изображает погрешность аналогичного расчета с переменным значением ε, величина которого только на участке индукции (t < 4.5 мкс) была уменьшена в 20 раз. В этом расчете потребовалось лишь на 12 шагов больше, но погрешность численного решения в зоне реакций и далее уменьшилась на 1.5–2 порядка. Обнаруженный эффект говорит о том, что точность расчета в зоне индукции, в которой происходит накопление изначально отсутствующих компонент смеси, существенно влияет на точность всего расчета.

Связь параметра ε с количеством N узлов расчетной сетки хорошо описывается соотношением Nε^1/p = const, где p – порядок аппроксимации контролирующей схемы. Это приводит к тому, что с ростом ε преимущество пары “8-6” над парой “6-4” в соотношении количества узлов уменьшается. С учетом более высоких вычислительных затрат на каждом шаге суммарное преимущество пары “8-6” над парой “6-4” исчезает при ε > 10^–4. При таких ε выгоднее использовать менее точную, но более быструю разностную схему. На фиг. 5 для ε = 10^–2, 10^–3 и 10^–4 приведены погрешности численного решения пары “6-4”, то есть метода 6-го порядка с контролем шага по методу 4-го порядка. Их коэффициенты указаны в разд. 1. Как и в предыдущем случае, в этих вариантах для параметра ε на участке индукции (t < 4.5 мкс) применялся коэффициент 1/20. В данных вариантах прогнозируемая точность выдерживается на всем промежутке расчета. Указанное в подписи к фигуре количество шагов N, которое потребовалось для каждого варианта, хорошо подчиняется зависимости Nε^1/4 = const.

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

ЗАКЛЮЧЕНИЕ

Основной целью данной работы являлось исследование эффективности алгоритма автоматического выбора шага интегрирования для mISD-схем высокой точности применительно к сложным математическим моделям с многокомпонентной химической кинетикой. Тестирование проводилось на задаче о взрыве водородно-кислородной смеси в замкнутом пространстве переменного объема. Показано, что при одинаковых требованиях к точности численного решения методы более высокого порядка позволяют в несколько раз сократить затраты машинного времени. Это оправдывает применение mISD-схем для приложений, требующих высокую точность расчета.