Расчёт концентраций для случая сложного равновесия
- В этой теме 2 ответа, 1 участник, последнее обновление 15 лет, 1 месяц назад сделано
Zheka.
Задача:
A+A <=> B; k1
B+A <=> C; k2
C+A <=> D; k3
Эти реакции характеризуются константами равновесия k1, k2, k3.
Равновесные концентрации задаются в виде c(A), c(B), c(C), c(D).
Начальные условия c(A0)=A0, c(B0)=0, c(C0)=0, c(D0)=0.
Расписав константы, имеем:
k1 = c(B) / c(A)^2;
k2 = c(C) / c(A)*c(B);
k3 = c(D) / c(A)*c(C);
Выразив из этих констант соответствующие концентрации
через cA, имеем:
cB = k1 * c(A)^2;
cC = k1*k2 * c(A)^3;
cD = k1*k2*k3 * c(A)^4;
Материальный баланс: c(A) + 2*c(B) + 3*c(C) + 4*c(D) = c(A0);
Получаем уравнение 4-й степени:
4*k1*k2*k3 *c(A)^4 + 3*k1*k2 *c(A)^3 + 2*k1 *c(A)^2 +c(A) — c(A0) = 0;
Выразив количество превратившегося вещества A через x, имеем:
c(A) = c(A0) — x
Подставив в выражение материального баланса, имеем:
4*k1*k2*k3 *(c(A0) — x)^4 + 3*k1*k2 *(c(A0) — x)^3 + 2*k1 *(c(A0) — x)^2 — x = 0;
Как бы его теперь правильно и с достаточной точностью решить?
Произвёл попытку решить эту задачу аналитически и численно в Maple:
### A+A<=>B (k1); B+A<=>C (k2); C+A<=>D (k3)
cA:=cA0-x; cB:=k1*cA^2; cC:=k1*k2*cA^3; cD:=k1*k2*k3*cA^4;
### Задаём константы
Conditions:={cA0=1e-40, k1=2.11E+27, k2=7.93E+13, k3=3.48E+10};
### Задаём уравнение материального баланса
MaterialBalance := cA+2*cB+3*cC+4*cD=cA0;
### Выводим его в красивой форме на экран
MaterialBalanceEquation:= collect(expand(lhs(MaterialBalance)-rhs(MaterialBalance)),x)=0;
### Решаем выражение аналитически относительно x
AnlSol:=solve(MaterialBalance,x);
### Решаем частный случай уравнения с подставленными константами
NumSol:=subs(Conditions, AnlSol);
### Вычисляем концентрации
ResultcA:=evalf(subs(Conditions, (cA0 — NumSol)) );
ResultcB:=evalf(subs(Conditions, k1*ResultcA^2));
ResultcC:=evalf(subs(Conditions, k1*k2*ResultcA^3));
ResultcD:=evalf(subs(Conditions, k1*k2*k3*ResultcA^4));
### Расчёт проверочных параметров
SummC:=ResultcA+2*ResultcB+3*ResultcC+4*ResultcD;
Const1:=ResultcB/(ResultcA^2);
Const2:=ResultcC/(ResultcB*ResultcA);
Const3:=ResultcD/(ResultcC*ResultcA);
После всех перестановок получается уравнение:
4*k1*k2*k3*x^4 + (-3*k1*k2-16*k1*k2*k3*cA0)*x^3 +
(2*k1+9*k1*k2*cA0+24*k1*k2*k3*cA0^2)*x^2 +
(-4*k1*cA0-9*k1*k2*cA0^2-1-16*k1*k2*k3*cA0^3)*x +
3*k1*k2*cA0^3+2*k1*cA0^2+4*k1*k2*k3*cA0^4 = 0
Задавая задачу в таком виде получается странное решение.
При нормальных начальных концентрациях c(A0) решение получается нулевым (а нужно посчитать для c(A0)=2.00 и 0.1 ).
При сверхнизких (1.00e-40) всё вроде сходится.
Обходной путь найден. Нужно лишь находить из материального баланса напрямую c(A),
эта переменная содержит малое значение, поэтому не округляется, в отличие от x.
### A+A<=>B (k1); B+A<=>C (k2); C+A<=>D (k3)
cB:=k1*cA^2; cC:=k1*k2*cA^3; cD:=k1*k2*k3*cA^4;
### Задаём константы и начальные условия
Conditions:={cA0=2.0, k1=2.11E+27, k2=7.93E+13, k3=3.48E+10};
### Задаём уравнение материального баланса
MaterialBalance := cA+2*cB+3*cC+4*cD=cA0;
### Выводим его в красивой форме на экран
MaterialBalanceEquation:= collect(expand(lhs(MaterialBalance)-rhs(MaterialBalance)),cA)=0;
### Решаем выражение аналитически относительно cA
AnlSol:=solve(MaterialBalance,cA);
### Решаем частный случай уравнения с подставленными константами
NumSol:=subs(Conditions, AnlSol);
### Вычисляем концентрации
ResultcA:=evalf(subs(Conditions, NumSol) );
ResultcB:=evalf(subs(Conditions, k1*ResultcA^2));
ResultcC:=evalf(subs(Conditions, k1*k2*ResultcA^3));
ResultcD:=evalf(subs(Conditions, k1*k2*k3*ResultcA^4));
### Расчёт проверочных параметров
SummC:=ResultcA+2*ResultcB+3*ResultcC+4*ResultcD;
Const1:=ResultcB/(ResultcA^2);
Const2:=ResultcC/(ResultcB*ResultcA);
Const3:=ResultcD/(ResultcC*ResultcA);
Вопрос нахождения равновесных концентраций закрыт. 8)
Для ответа в этой теме необходимо авторизоваться.