Численное моделирование жестких систем с использованием (𝟒, 𝟐)-метода

Аннотация

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

Полный текст

Моделирование процесса термического пиролиза этана для получения этилена является актуальной задачей. При решении таких задач требуется найти решение задачи Коши для систем обыкновенных дифференциальных уравнений, составленных на основе схемы реакций. Большая часть этих систем является жесткой по части фазовых переменных. Выполнение свойств жесткой системы влечет за собой большое различие в величине коэффициентов, что, в свою очередь, не позволяет применить к решению задачи Коши стандартные явные методы. В работах [1; 2] подробно описано построение (𝑚, 𝑘)-методов для решения жестких систем дифференциальных уравнений, основанных на методах типа Розенброка. В данной работе описан (4,2)-метод, примененный к решению задачи Коши для нелинейной системы обыкновенных дифференциальных уравнений, соответствующей математической модели брутто-реакции пиролиза этана.

Рассмотрим брутто-реакцию пиролиза этана [3; 4]:

$\begin{array}{l}{C}_2{H}_6\to C_2{H}_4+H_2,\\ 2C_2{H}_6\to C_2{H}_4+2C{H}_4.\end{array}$

Соответствующие кинетические уравнения имеют вид

$\begin{array}{l}\left[C_2{H}_6\right]\to \left[C_2{H}_4\right]+\left[H_2\right],{\omega }_1=k_1{c}_1,\\ 2\left[C_2{H}_6\right]\to \left[C_2{H}_4\right]+2\left[C{H}_4\right],{\omega }_2=k_2{c}_1^2,\end{array}$

где 𝑐_𝑖 — концентрация -того вещества, 𝑘_𝑗 — константа скорости 𝑗-той прямой стадии, которые вычисляются согласно выражениям Аррениуса при температуре 800 К, таким образом 𝑘₁ = 0.051, 𝑘₂ = 0.7.

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

$\left\{\begin{array}{l}{\dot{c}}_1=-k_1{c}_1-2k_2{c}_1^2\\ {\dot{c}}_2=k_1{c}_1+k_2{c}_1^2\\ {\dot{c}}_3=k_1{c}_1\\ {\dot{c}}_4=2k_2{c}_1^2\end{array}\right.$                                                                             (1)

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

$c_1\left(0\right)=\mathrm{1,}{c}_2\left(0\right)=\mathrm{0,}{c}_3\left(0\right)=\mathrm{0,}{c}_4\left(0\right)=0.$                                                      (2)

Для нахождения изменений количественных характеристик реакции численно решим задачу Коши (1)-(2) (𝑚, 𝑘)-методом и проведем сравнительную характеристику со способом решения методом Эйлера.

Используя (4,2)-метод четвертого порядка точности, разностную схему для каждого компонента вектора решения построенной модели можно представить в виде [1; 2]:

$\begin{array}{c}{y}^{n+1}=y^n+\sum _{i=1}^4{p}_i{q}_i,D_n=E-a\tau {f^{\text{'}}}_n,\\ D_n{q}_1=\tau f\left(y_n\right),\\ D_n{q}_2=q_1,\\ D_n{q}_3=\tau f(y_n+{\beta }_{31}{q}_1+{\beta }_{32}{q}_2)+a_{32}{q}_2,\\ D_n{q}_4=q_3+a_{42}{q}_2,\end{array}$

где 𝜏 – шаг интегрирования решаемой задачи, 𝐸 − единичная матрица размерности (4×4), $y\in R^4,q_i\in R^4,y_0=\left(\mathrm{1,0,0,0}\right)$.

Вещественные коэффициенты, при которых схема обладает устойчивостью и 4-м порядком точности, $a,p_i,a_{j2},{\beta }_{3k}\mathrm{,1}\le i\le \mathrm{4,3}\le j\le \mathrm{4,1}\le q\le 2$, имеют следующие значения [1; 2]:

$\begin{array}{l}a=\mathrm{0.57281606248213,} p_i=\mathrm{1.27836939012447,}\\ p_2=-\mathrm{1.00738680980438,}{p}_3=\mathrm{0.9265539109350,}\\ p_4=-\mathrm{0.33396131834691,}{\beta }_{31}=\mathrm{1.00900469029922,}\\ {\beta }_{32}=-\mathrm{0.25900469029921,}{a}_{32}=-\mathrm{0.49552206416578,}\\ a_{42}=-\mathrm{1.28777648233922.}\end{array}$

Для системы (1) 𝑓(𝑦) 𝑓′_n и имеют вид

$f\left(y\right)=$$\left(\begin{array}{c}-k_1{y}_1-2k_2{y}_1^2\\ k_1{y}_1+2k_2{y}_1^2\\ k_1{y}_1\\ 2k_2{y}_1^2\end{array}\right),$

${f^{\text{'}}}_n=\frac{\partial f\left(y\right)}{\partial y}=$$\left(\begin{array}{cccc}-k_1-4k_2{y}_1& k_1+2k_2{y}_1& k_1{y}_1& 4k_2{y}_1\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right).$

На основе построенной схемы был разработан программный комплекс для решения системы уравнений (1) с начальными данными (2). На рисунке представлены изменения концентраций веществ во времени при температуре 800 К. Для оценки эффективности (4,2)- метода эта же система была решена явным методом Эйлера, разностная схема для которого имеет вид [5]:

Сравнительные результаты испытаний представлены в таблице 1, за опорное было принято решение, полученное с помощью системы Wolfram Mathematica с шагом 𝜏 = 0.1.

$y^{n+1}=y^n+\tau f\left(y^n\right).$

Таблица 1

Результаты работы методов

Число итераций	Время выполнения		Точность, 𝜀
	Явный метод Эйлера	(4,2) − метод	Явный метод Эйлера	(4,2) − метод
10	6 мс	7 мс	0.034	0.001
50	59 мс	49 мс	0.00762	0.00055
93	105 мс	98 мс	0.0043	0.0001

Из таблицы 1 видно, что многостадийный метод имеет преимущество по времени выполнения, точности и корректности по отношению к явному методу.

Рис. 1. Изменение концентраций веществ при 800 К.

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