Как решить систему из 8 ODE в Matlab, где для 4 ODE начальное значение находится на высоте z = 0, а для остальных начальное значение находится на z = H - PullRequest
Новая
       77

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

0 голосов
/ 05 октября 2018

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

Я привел их начальные значения (т.е. я знаю начальные значения твердых тел вверху и начальные значения газов внизу).Для решения уравнений баланса тепла и массы я написал 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

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

...