Как решить систему из 8 ОДУ в Matlab, где для 4 ОДУ начальное значение находится на высоте z = 0, а для остальных начальное значение находится на z = H

У меня есть реактор, в который твердое вещество загружается сверху, а газ - снизу.

Я предоставил их начальные значения (т.е. я знаю начальное значение твердых объектов вверху и начальные значения газов внизу). Для решения уравнений баланса тепла и массы я написал 8 ОДУ.

Но как решить их одновременно?
Некоторое начальное значение ODE находится на z=0, а для некоторого другого начального значения z оно находится на z = H;, где H - высота печи.

function f = odefun( Z,Y )

  % FeO + CO -> Fe + CO2 ---(1)
  % C + 1/2 O2 -> CO ---(2)
  % C + O2 -> CO2  ---(3)
  % C + CO2 -> 2CO ---(4)
 
  Fg    = Y(1);     %volumetic flow rate of gas Nm3s-1 //initialize
  Cfeo  = Y(2);     % initial FeO concentration mol/m3
  Cc    = Y(3);     % initial C concentration mol/m3
  yco   = Y(4);     % mole fraction of CO initially
  yco2  = Y(5);     % mole fraction of CO2 initially
  yo2   = Y(6);     % mole fraction of O2 initially
  Tg    = Y(7);     % temp of gas K
  Ts    = Y(8);     % temp of solid K
 
  Fs = 0.044;      %Mass flow rate of solid m3/s //initialize and constant

 P = 5;         %total pressure atm
 R = 8.314;     %gas constant
 M = 104;       %molecular mass //constant
 keqm1 = 3.33;     %equilibrium constant for rx 1
 keqm2 = 3 * 10^10;     %equilibrium constant for rx 2
 keqm3 = 4.75 * 10^20;     %equilibrium constant for rx 3
 
 Az = 4;        %area of the bed m2 //constant
 e = 0.4;       %voidage //constant
 O = 0.6;       %Avg shape factor of solid particle //constant
 dp = 0.02;     %diameter of particles // constant m
 Cpg = 293000;    %specific heat capacity of gas (JKg-1K-1) //constant
 Cps = 440000;    %specific heat capacity of solid (JKg-1K-1) //assume constant
 
 mu = 0.005;        %viscosity of gas Kg m-1 s-1 //assume constant
 
 a = 10;        %surfce area of DRI per unit volume of bed
 Ac = 20;       % surfce area of coal per unit volume
 hgs = 200;      %heat transfer coefficient b/w solid and gas J m-2 s-1 K-1

 H1 = -16000;     %enthalpy of reaction 1 KJ mol-1
 H2 = -110520;    %enthalpy of reaction 2 KJ mol-1
 H3 = -393509;    %enthalpy of reaction 3 KJ mol-1
 H4 = -282989;      %enthalpy of reaction 4 KJ mol-1
 
 
 C0 = 28/22.4;    % Kg/m3
 C1 = 44/22.4;
 C2 = 32/22.4;
 
 Dg = (C0 * yco + C1 * yco2 + C2 * yo2); %Density of Gas Kg m-3
 Db = (1-e)*(Cfeo*72 + (Y(2) - Cfeo)*56 + Cc*12 )*1000 ; %Density of Bed Kg m-3
 Re = (Dg * Fg * e * dp)/(mu*Az) ; %Reynolds number
 
 f1 = (150 * (1-e)^2 * mu)/(e^2 * Dg * dp^2); %friction factor parameter
 f2 = (1.75*(1-e))/(e^3 * Dg * dp);
 
 dFgdZ = -(Fg/(P * Az^2))*(f1*Az*Fg + f2 * Fg^2); %volumetric flow rate of gas Nm3 s-1
 
 kp1 = 2.5*exp(-8857.6/Ts);  % chemical rxn rate constant s-1
 kp2 = 6.52*10^(5) * exp(-22000/Ts) * Ts^(0.5);
 kp3 = kp2;
 kp4 = 4*10^12 * exp(-40000/Ts);
 
 Vg = yco + yco2; %mass fraction of volatile evolve 
 kf = 2.0*Re^(-0.336)* Vg; %mass transfer coefficient m/s
 
 Po2 = yo2; % partial pressure atm
 Pco = yco;
 Pco2 = yco2;
 Pco2eqm = 0.2;  % Change  
 
 Pre2 = (Po2 - (Pco^(2)/keqm2^(2))); % driving force atm
 Pre3 = (Po2 - (Pco2^(2)/keqm3^(2)));
 Pre4 = (Pco2 - Pco2eqm);
 
 R1 = a*(Fs/Y(1))^(1/3)*kp1*(yco-yco2/keqm1)*P/(R*Ts);  %reaction rate constant m-3s-1
 R2 = (Ac/((6/(dp*kp2))+(1/kf))) * (9.869*10^(-6)*Pre2)/(M*R*Ts);
 R3 = (Ac/((6/(dp*kp3))+(1/kf))) * (9.869*10^(-6)*Pre3)/(M*R*Ts);
 R4 = (Ac/((6/(dp*kp4))+(1/kf))) * (9.869*10^(-6)*Pre4)/(M*R*Ts);


dCfeodZ = (R1*Az*(1-e))/Fs; %for FeO Cfeo Kg/s
dCcdZ = ((R2+R3+R4)*Az*(1-e))/Fs; %for C

dycodZ  = 22.4*Az*(R2*(yco/2 - 1)+(yco - 2)+R1)/Fg; %for CO
dyco2dZ = 22.4*Az*(yco2*(R2/2 + R4)+R4-R1-R3)/Fg;  %for CO2
dyo2dZ  = 22.4*Az*(R2/2 + R3 + yo2*(R2/2 + R4))/Fg; %for O2


%Heat transfer equations

S = (6*(1-e)*Az)/(dp*O) ; %surface area of the solid

B = (6*(1-e)*hgs)/(dp*O);
A = 16*R1+12*(R2+R3+R4);
E = H1*R1+ H2*R2+ H3*R3+ H4*R4;

dTgdZ = (22.4/(Fg*Cpg))*(Az*Cpg*Tg*(R2/2 + R4)+ S*(Tg - Ts)*hgs );
dTsdZ = (Az*(1-e)/(Fs*Cps*Db))*(A*Cps*Ts + B*(Tg - Ts) - E);



f = [dFgdZ;dCfeodZ;dCcdZ;dycodZ;dyco2dZ;dyo2dZ;dTgdZ;dTsdZ];


end

Это функция с ОДУ и начальными значениями Fg (z = 0), Cfeo (z = H), Cc (z = H) yco (z = 0), yco2 (z = 0), yo2 (z = 0 ), Tg (z = 0), Ts (z = H) любезно помогли мне решить эту проблему.

Это краевая задача (BVP), для которой требуется решатель BVP. Вы найдете это ключевое слово в документации по Matlab с примерами приложений.

Lutz Lehmann 05.10.2018 21:51

Если я приму это как BVP, то мне потребуются 2 начальных условия (т. Е. При z = 0 и при z = H) для всех ODE. Но для некоторых у меня есть при z = 0, а для других - при z = H. Можете прислать материалы. Это срочно.

Aman 07.10.2018 10:19

Нет, вы получаете столько же граничных условий, сколько состояние имеет размеры. То, что граничные условия могут быть построены из полного состояния в 0 и H, не противоречит первому утверждению.

Lutz Lehmann 07.10.2018 10:45

Можете ли вы дать мне решение или код для этой проблемы.

Aman 22.10.2018 06:40
Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
2
4
53
0

Другие вопросы по теме