Глава 3. Компьютерное Моделирование и Оптимизация хтп

Оптимальное время проведения химической реакции в аппарате идеального смешения

Рассчитать оптимальное время проведения химической реакции в аппарате идеального смешения, приняв в качестве критерия оптимальности выход целевого продукта P.

Схема реакции: . Порядок обеих стадий реакции – первый. Константы скоростей равны:

Решение

Материальный баланс по компонентам A и P:

При делении уравнений

на расход реагента v получаем:

где  - среднее время пребывания реагентов в реакторе

Выход продукта P выражается:

Необходимое условие

существования экстремума: или

Поскольку  и

Условие экстремума будет иметь вид:

Откуда:  или

Программа: Model 1_ tau ( построение графиков xa от T и xp от T. Графическое решение задачи. Поиск экстремума. Программа расчета оптимального времени пребывания (tau - час) в изотермическом реакторе идеального смешения со стехиометрической схемой реакции A - P –S)

Файл DATA. m Программный код файла DATA.m - задание исходной информации для расчетов

function DATA

global xa0 k1 k2 tau_a tau_b;

Концентрация реагента А в долях в реакции A-P-S (мольные доли)

xa0=1;

Константы скоростей реакций (час^(-1))

k1=0.35; k2=0.13;

Левая граница поиска оптимального времени пребывания в реакторе (час)

tau_a=1;

Правая граница поиска оптимального времени пребывания в реакторе (час)

tau_b=10;

end

Файл model 1_ stat _ tau. m Программный код файла model1_stat_tau.m - расчет выходных концентраций продуктов из реактора

Значения вектора х функции model 1_ stat _ tau зависят от параметра «tau»

function x= model1_stat_tau(tau)

globalxa0 k1 k2;

Обнуляем матрицу a и вектор b

a=zeros(2,2);b=zeros(2,1);

Для решения системы уравнений используем метод обратной матрицы.

Преобразуем систему уравнений:

чтобы получить коэффициенты матрицы:

Задание коэффициентов матрицы А (коэффициенты перед и ) и

вектора В (коэффициенты правой части системы уравнений,

после знака «равно»).

a(1,1)=1+k1*tau;a(1,2)=0;b(1)=xa0;

a(2,1)=tau*k1;a(2,2)=(-1-tau*k2); b(2)=0;

Определение выходных параметров модели. А=х*В. х=А/В. х=А^-1*В

x=inv(a)*b;

Оператор inv (а) - поиск обратной матрицы А

end

Файл GLAV _ maximum _ phi _ p. m Программный код файла GLAV_maximum_phi_p.m - главная управляющая программа

clc;

clear all;

close all;

global t xa xp tau_a tau_b;

i=0;

DATA;

Задаем значение коэффициента tau (время пребывания в реакторе, мин) в интервале от tau_a до tau_b с шагом 0.1

for tau=tau_a:0.1:tau_b

 i=i+1;  

Вызываемфункцию x= model1_stat_tau(tau)

x= model1_stat_tau(tau);

Массив значений коэффициента tau (мин) в интервале от tau _ a до tau _ b с шагом 0.1

t(i)=tau;

Массив значений коэффициента ха (мол.дол) - концентрации продукта А на выходе из реактора в заданном интервале температур от tau _ a до tau _ b

xa(i)=x(1);

Массив значений коэффициента ха (мол.дол) - концентрации продукта P на выходе из реактора в заданном интервале температур от tau _ a до tau _ b

xp(i)=x(2);

end

REPORT;

Файл REPORT. m Программный код файла REPORT.m - отчет о работе программы

global xa0 k1 k2 t xa xp tau_a tau_b;

Выводим на экран необходимый текст.

disp('ПРОГРАММА МОДЕЛИРОВАНИЯ ИЗОТЕРМИЧЕСКОГО РЕАКТОРА ИДЕАЛЬНОГО ПЕРЕМЕШИВАНИЯ ');

disp('          РЕАКЦИЯ: A - P - S');

disp('-------------------------------------------------');

disp('Программавключаетследующиефайлы: GLAV_model1_grafik.m+DATA.m+model1_stat_tau+REPORT.m');

disp('_________________________________________________ ');

disp('ИСХОДНЫЕ ДАННЫЕ ');

disp('-------------------------------------------------');

Оператор num2str(xa0,'%10.2f') переводит числовые значения в строковые.

Для этого необходимо записать оператор num2str, далее в скобках указываем переменную, которую необходимо перевести в строковый формат «xa0», затем ставим запятую «,» либо пробел и в одинарных кавычках через точку указываем число символов (%10)в строковой переменной и количество знаков после запятой(2f).

disp(['1.Концентрация реагента A на входе в реактор (xa0) = ' num2str(xa0,'%10.2f') ' мольные доли']);

disp('--------------------------------------------------');

disp(['2.Константа скорости первой элементарной реакции (k1) = ' num2str(k1,'%10.2f') ' час^(-1)']);

disp(['3.Константа скорости второй элементарной реакции (k2) = ' num2str(k2,'%10.2f') ' час^(-1)']);

disp('---------------------------------------------------');

disp(['4.Левая граница интервала поиска оптимального времени пребывания в реакторе (tau_a) = ' num2str(tau_a,'%10.2f') ' час']);

disp('--------------------------------------------------');

disp(['5. Правая граница интервала поиска оптимального времени пребывания в реакторе(tau_b) = ' num2str(tau_b,'%10.2f') ' час']);

disp('___________________________________________________ ');

disp('---------------------------------------------');

disp ('РЕЗУЛЬТАТЫ РАСЧЕТОВ');

disp('------------------------------------------------');

disp('--------------------------------------------');

disp ('1.Концентрации продуктов на выходе из реактора');

disp('--------------------------------------------------');

disp('_____________________________')

disp(' tau(час)        xa(мол.д.)  xp(мол.д.) '); 

disp('_____________________________')

i=0;

Выводим таблицу данных. Переменные t, xa, xp;

for i=1:length(t)

   tt=t(i);xaa=xa(i);xpp=xp(i);

disp(sprintf('%10.3f\t %10.5f\t %10.5f',tt,xaa,xpp));

end

Построение графиков

hfig(1)=figure;

hfig(2)=figure;

Рисунок 1. График зависимости выхода продукта A от времени пребывания в реакторе.

figure(hfig(1));

Построение графика

plot(t,xa);

Название графика

title('График зависимости выхода продукта A от времени пребывания в реакторе');

Подписи осей графика.

xlabel('tau - час');ylabel('Выходные концентрации: xa - мольные доли');

Наложение сетки на график

gridon;

Рисунок 2. График зависимости выхода продукта P от времени пребывания в реакторе

figure(hfig(2));

Построение графика

plot(t,xp);

Название графика

title('График зависимости выхода продукта P от времени пребывания в реакторе');

Подписи осей графика.

xlabel('tau - час');ylabel('Выходные концентрации: xp - мольные доли');

Наложение сетки на график

grid on;

end

Optim1_tau (Оптимальные значения режимных параметров реактора: времени пребывания, концентрации продукта А на выходе из реактора, концентрации продукта Р на выходе из реактора. Программа расчета оптимального времени пребывания (tau - час) в изотермическом реакторе идеального смешения со стехиометрической схемой реакции A - P –S)

Файл DATA. m Программный код файла DATA.m - задание исходной информации для расчетов

function DATA

global xa0 k1 k2 tau_a tau_b;

Концентрация реагента А в долях в реакции A-P-S (мольные доли)

xa0=1;

Константы скоростей реакций (час^(-1))

k1=0.35; k2=0.13;

Левая граница поиска оптимального времени пребывания в реакторе (час)

tau_a=1;

Правая граница поиска оптимального времени пребывания в реакторе (час)

tau_b=5;

end

Файл model 1_ stat _ tau. m Программный код файла model1_stat_tau.m - расчет выходных концентраций продуктов из реактора

Значения параметра phi _ p функции model 1_ stat _ tau зависят от параметра «tau»

function phi_p= model1_stat_tau(tau)

global xa0 k1 k2 xa xp;

Обнуляем матрицу a и вектор b

a=zeros(2,2);b=zeros(2,1);

Для решения системы уравнений используем метод обратной матрицы.

Преобразуем систему уравнений:

чтобы получить коэффициенты матрицы:

Задание коэффициентов матрицы А (коэффициенты перед и ) и

вектора В (коэффициенты правой части системы уравнений,

после знака «равно»).

a(1,1)=1+k1*tau;a(1,2)=0;b(1)=xa0;

a(2,1)=tau*k1;a(2,2)=(-1-tau*k2); b(2)=0;

Определение выходных параметров модели. А=х*В. х=А/В. х=А^-1*В

Оператор inv (а) - поиск обратной матрицы А

x=inv(a)*b; xa=x(1);xp=x(2);

Вычисление критерия оптимальности phi _ p (целевой функции);

В дальнейшем будет использоваться оператор «fminbnd», который находит минимум функции. Для того чтобы найти максимум функции с помощью оператора «fminbnd», мы ставим знак «-» перед x (2)/ xa 0

phi_p=-x(2)/xa0;

end

Файл GLAV _ maximum _ phi _ p. m Программный код файла GLAV_maximum_phi_p.m - главная управляющая программа

clc;

clear all;

close all;

global tau_opt phi_p_max tau_a tau_b;

DATA;

disp('Информация об итерационном процессе расчета');

Информация об итерационном процессе расчета

iteration=optimset('Display','iter');

Реализация алгоритма поиска максимума функции одной переменной

[tau_opt,phi_p]=fminbnd('model1_stat_tau',tau_a,tau_b,iteration);

Знак «-» перед phi _ p необходим для корректного вывода данных экстремума функции. Оператор «fminbnd» находит минимум функции со знаком «-». Нам нужен максимум функции со знаком «+».

phi_p_max=-phi_p;

REPORT;

Файл REPORT Программный код файла REPORT.m - отчет о работе программы

function REPORT

global xa0 k1 k2 tau_opt phi_p_max xa xp tau_a tau_b;

Выводим на экран необходимый текст.

disp('ПРОГРАММА ОПРЕДЕЛЕНИЯ ОПТИМАЛЬНОГО ВРЕМЕНИ ПРЕБЫВАНИЯ');

disp('В ИЗОТЕРМИЧЕСКОМ РЕАКТОРE ИДЕАЛЬНОГО ПЕРЕМЕШИВАНИЯ ');

disp('                  РЕАКЦИЯ: A - P - S');

disp('-------------------------------------------------');

disp('Программа: GLAV_maximum_phi_p.m+DATA.m+model1_stat_tau.m+REPORT.m');

disp('___________________________________________________ ');

disp('ИСХОДНЫЕ ДАННЫЕ ');

disp('-------------------------------------------------');

Оператор num2str(xa0,'%10.2f') переводит числовые значения в строковые.

Для этого необходимо записать оператор num2str, далее в скобках указываем переменную, которую необходимо перевести в строковый формат «xa0», затем ставим запятую «,» либо пробел и в одинарных кавычках через точку указываем число символов (%10)в строковой переменной и количество знаков после запятой(2f).

disp(['1.Концентрация реагента A на входе в реактор (xa0) = ' num2str(xa0,'%10.2f') ' мольные доли']);

disp('-------------------------------------------');

disp(['2.Константа скорости первой элементарной реакции (k1) = ' num2str(k1,'%10.2f') ' час^(-1)']);

disp('---------------------------------------------');

disp(['3.Константа скорости второй элементарной реакции (k2) = ' num2str(k2,'%10.2f') ' час^(-1)']);

disp('-----------------------------------------------');

disp(['4.Левая граница интервала поиска оптимального времени пребывания в реакторе (tau_a) = ' num2str(tau_a,'%10.2f') ' час']);

disp('-----------------------------------------------');

disp(['5. Правая граница интервала поиска оптимального времени пребывания в реакторе(tau_b) = ' num2str(tau_b,'%10.2f') ' час']);

disp('____________________________________________________ ');

disp('---------------------------------------------');

disp ('РЕЗУЛЬТАТЫ РАСЧЕТОВ');

disp('------------------------------------------------');

disp('-----------------------------------------------------');

disp ('1.Оптимальные значения режимных параметров реактора');

disp('-----------------------------------------------');

disp('____________________________________________________')

disp(' tau_opt(час) phi_p_max xa(мол.д.) xp(мол.д.) '); 

disp('________________________________________________');

Выводим таблицу данных. Переменные tau _ opt, phi _ p _ max, xa, xp;

disp(sprintf('%10.3f\t %10.5f\t %10.5f\t %10.5f',tau_opt,phi_p_max,xa,xp));

end