12.3.2. Численное решение систем дифференциальных уравнений

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

(136)

с начальными условиями ,.

Численная схема решения системы уравнений (136) согласно методуЭйлерабудет выглядеть:

(137)

где t– интервал дискретизации;– номер интервала;– количество интервалов; с начальными условиями,.

Численная схема решения системы уравнений (136) согласно модифицированному методу Эйлерабудет выглядеть:

(138)

где x_i*,y_i*,t_i*– промежуточные точки, рассчитываемые следующим образом:

,

,

.

Если подставить значения промежуточных точек в формулы (138), получим:

(139)

Численная схема решения системы уравнений (136) согласно методу Эйлера-Кошибудет выглядеть:

(140)

где x_i*,y_i*,t_i*– промежуточные точки, рассчитываемые следующим образом:

,

.

Если подставить значения промежуточных точек в формулы (140), получим:

(141)

Численная схема решения системы уравнений (136) согласно методуРунге-Кутта 4-го порядка будет выглядеть:

(142)

где X1_i,X2_i,X3_i,X4_i,Y1_i,Y2_i,Y3_i,Y4_i– промежуточные точки, рассчитываемые следующим образом:

Пример. Имеется система нелинейных дифференциальных уравнений:

(143)

Задача: составить численные схемы решения системы уравнений (143).

Метод Эйлера.

(144)

Модифицированный метод Эйлера.

(145)

или

(146)

Метод Эйлера-Коши.

(147)

или

(148)

Метод Рунге-Кутта 4-го порядка.

(149)

12.4. Пример выполнения

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

.

Численное интегрирование осуществим методом Эйлера при следующих начальных условиях:

,

,

,

U=1 м/с,l=2 м,l= 0,1 м.

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

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

(150)

где X_A, X_P, X_D– концентрации реагентов (исходных продуктов или продуктов, образующихся в зоне реактора);_A, _P, _D– скорости образования (исчерпывания) веществ, с^-1;l– линейная координата, м;U– линейная скорость подачи веществ в зону реакции, м/с.

Для реакции, приведенной в задании, можно выделить следующие элементарные стадии:

; .

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

,

, (151)

,

где ;.

Модель реактора идеального вытеснения (150) с учетом (151) описывается системой дифференциальных уравнений:

(152)

Начальные условия: ,,.

Проведем аппроксимацию полученной системы дифференциальных уравнений конечно-разностными уравнениями на интервале l[0,L] с шагомl в соответствии с рассмотренными выше методами. Полученная система конечно-разностных уравнений позволит осуществить расчет профиля концентраций по длине реактора (рис. 99).

Численная схема, полученная согласно методу Эйлера:

(153)

где ,, с начальными условиями,,.

Численная схема, полученная согласно модифицированному методу Эйлера:

(154)

или

Численная схема, полученная согласно модифицированному методу Эйлера-Коши:

(155)

или

Численная схема, полученная согласно методу Рунге-Кутта 4-го порядка:

(156)

Рис. 99. Блок-схема алгоритма расчета профиля концентраций