I'm trying to solve a 2nd order related ODE system with five dependent variables.
The code is:
Nc = 8; Nn = Nc-4; Pt = 0.495; R = 0.08206; Ru = 8.314; Cebo = 0.992; Cwo = 5.76; Pebo = 0.073;
Pwo = 0.422; yebo = 0.1428; ywo = 0.8571; T0 = 893.15;
yeb[z_] := (Ceb[z] R T[z])/(Pt 1000); yw = ywo; yst[z_] := Abs[(1 - Ceb[z] - ywo)/2];
c2[z_] := (yw Pt 1000)/(R T[z]); c3[z_] := (yst[z] Pt 1000)/(R T[z]); c4[z_] := c3[z];
Mmeb = 106.17; Mmw = 18; Mmst = 104.15; Mmh2 = 2; Mmbe = 78.11; Mmet = 28.05; Mmto = 92.14; Mmme= 16.04;
M[1] = Mmeb; M[2] = Mmw; M[3] = Mmst; M[4] = Mmh2; M[5] = Mmbe; M[6] = Mmet; M[7] = Mmto; M[8] = Mmme;
j = 0;
Do[i = 1; j = j + 1;
M[i_,j_] := 2*(1/M[i] + 1/M[j])^-1, Nc];
Teb = 409.15; Tw = 373.15; Tst= 418.15; Th2 = 20.27; Tbe = 353.25; Tet = 169.5; Tto = 384.00; Tme= 111.65;
Tb[1] = Teb; Tb[2] = Tw; Tb[3] = Tst; Tb[4] = Th2; Tb[5] = Tbe; Tb[6] = Tet; Tb[7] = Tto;Tb[8] = Tme;
Veb = 139.24; Vw = 18.789; Vst= 131.27; Vh2 = 28.16; Vbe = 96.017; Vet = 49.29; Vto = 118.29; Vme= 35.64;
Vm[1] = Veb; Vm[2] = Vw; Vm[3] = Vst; Vm[4] =Vh2; Vm[5] = Vbe; Vm[6] = Vet; Vm[7] = Vto; Vm[8] = Vme;
Deb = 0.58; Dw = 1.8546; Dst = 0.13; Dh2 = 0; Dbe = 0; Dett = 0; Dto = 0.36; Dme = 0;
Dp[1] = Deb; Dp[2] = Dw; Dp[3] = Dst; Dp[4] = Dh2; Dp[5] = Dbe; Dp[6] = Det; Dp[7] = Dto; Dp[8] = Dme;
i = 0;
Do[i = i + 1;
DD[i_] = (1940*Dp[i]^2)/(Vm[i]*Tb[i]); SS[i_] = ((1.585*Vm[i])/(1 + 1.3*DD[i]^2))^(1/3); EK[i_] = 1.18 (1 + 1.3*DD[i]^2)*Tb[i];
CO[i_,z_] := (1.16145*(T[z]/EK[i])^-0.14874) + 0.52487*(Exp[-0.77320 (T[z]/EK[i])]) + 2.16178*(Exp[-2.43787 (T[z]/EK[i])]);
Mu[i_,z_] := (2.6709 10^-6)*Sqrt[M[i]*T[z]]/((SS[i]^2)*CO[i,z]),Nc];
y[1,z_] := yeb[z]; y[2,z_] := yw; y[3,z_] := yst[z]; y[4,z_] := yst[z];
i = 0; j = 0;
Do[i = i + 1;
PHI[i,i,z] = 1;
Do[j = j + 1;
PHI[i_,j_,z_] := (1/Sqrt[8]*(1 + M[i]/M[j])^-0.5)*(1 + ((Mu[i,z]/Mu[j,z])^0.5)*(M[j]/M[i])^0.25)^2, Nc],Nc];
Mum[z_] := \!\(\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(4\)]\((\*FractionBox[\(y[i, z]*Mu[i, z]\), \(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(4\)]\(\(y[j,z]\)\(\ \)\(PHI[i, j, z]\)\(\ \)\)\)])\)\);
Daeb := 8*15.9 + 10*2.31; Daw := 4.54 + 2*2.31; Dast := 8*15.9 + 8*2.31; Dah2 := 2*2.31;
Da[1] = Daeb; Da[2] = Daw; Da[3] = Dast; Da[4] = Dah2;
i = 1; j = 0;
Do[j = j + 1;
Dc[i_,j_,z_] := (0.00143*T[z]^1.75)/(Pt M[i,j]*(Da[i]^(1/3) + Da[j]^(1/3))^2), Nc];
aa = 1; bb = 0; cc = 1; dd = 1;
ZZ= (cc + dd - aa - bb)/aa;
Dm[z_] := 1/(1/(1 + ZZ*y[1,z])*((1/Dc[1,2,z])*(y[2,z] - bb/aa y[1,z]) + ((1/Dc[1,3,z])*y[3,z] + cc/aa y[1,z]) + ((1/Dc[1,4,z])*(y[4,z] + dd/aa y[1,z]))));
n = 3;
Cp[1,z_] := (5 + n)*(1/2)*((8.314*T[z]^0)/M[1]); Cp[2,z_] := 143.05 - 183.54*(T[z]/100)^0.25 + 82.751*(T[z]/100)^0.5 - 3.6989*(T[z]/100)^1; Cp[3,z_] := (5 + n)*(1/2)*(8.314 /M[3])*T[z]^0 ; Cp[4,z_] := 56.505 - 702.24 (T[z]/100)^-0.75 + 1165.0 (T[z]/100)^-1 - 560.70 (T[z]/100)^-1.5;
Cpm[z_] := \!\(\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(4\)]\ \(Cp[i, z]\ y[i,z]\)\);
i = 0;
j = 0;
Do[i = i + 1;
k[i_,z_] := ((8.3127*10^-2)*Sqrt[(T[z]/M[i])])/((SS[i]^2)*CO[i,z]) + 1.32*(Cp[i,z] - (5/2)*8.314/M[i])*((2.6709*10^-6)*Sqrt[M[i]*T[z]])/((SS[i]^2)*CO[i,z]),Nn];
km[z_] := \!\(\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(4\)]\((\*FractionBox[\(y[i, z]\ k[i, z]\), \(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(4\)]\(\(y[j,z]\)\(\ \)\(PHI[i, j, z]\)\(\ \)\)\)])\)\);
CDeb = -102.22; CDst = -104.56; CDh2 = -117.95;
Ea[1] = 175.38; Ea[2] = 296.29; Ea[3] = 474.76; Ea[4] = 213.78;
Aeb = (1.014*10^-5)/0.9869; Ast = (2.678*10^-5)/0.9869; Ah2 = (4.519*10^-7)/0.9869;
A[1] = 5.594*10^9; A[2] = 1.060*10^15; A[3] = 1.246*10^26; A[4] = 8.024*10^10;
i = 0;
Do[i = i + 1; kc[i_,z_] := A[i]*Exp[-Ea[i]/(Ru/1000 T[z])],Nn];
keb[z_] := Aeb*Exp[-CDeb/(Ru/1000 T[z])]; kst[z_] := Ast*Exp[-CDst/(Ru/1000 T[z])]; kh2[z_] := Ah2*Exp[-CDh2/(Ru/1000 T[z])];
KEQ[z_] := (y[3,z]*Pt*y[4,z]*Pt)/(y[1,z]*Pt);
r1[z_] := (kc[1,z]*keb[z]*(y[1,z]*Pt - (y[3,z]*Pt*y[4,z]*Pt)/(KEQ[z]*Pt)))/(1 + keb[z]*y[1,z]*Pt + kh2[z]*y[4,z]*Pt + kst[z]*y[3,z]*Pt)^2;
r2[z_]:= (kc[2,z]*keb[z]*y[1,z]*Pt)/(1 +keb[z]*y[1,z]*Pt + kh2[z]*y[4,z]*Pt + kst[z]*y[3,z]*Pt)^2;
r3[z_] := (kc[3,z]*keb[z]*y[1,z]*Pt*kh2[z]*y[4,z]*Pt)/(1 + keb[z]*y[1,z]*Pt + kh2[z]*y[4,z]*Pt + kst[z]*y[3,z]*Pt)^2;
CDH1 = 124.83; CDH2=101.50; CDH3=-65.06;
u=((0.005 + 0.00585)/2)*3600; q = (35.52 + 31.8)/2; L=6;
eqn1 = Dm[z]*Ceb''[z] - u*Ceb'[z] == 1442*r1[z];
eqn2 = km[z]*T''[z] - q*Cpm[z]*T'[z] == 1442*(r1[z]*(-CDH1) + r2[z]*(-CDH2) + r3[z]*(-CDH3));
initConds={T[0] == 893.15, T'[L] == 0, Ceb[0] == 0.992, Ceb'[L] == 0 };
eqns=Join[{eqn1,eqn2}];
NDSolve[{eqns,initConds},{Ceb[z],T[z]},{z,0,6,0.01}]
I'm using to use NDSolve
as shown below, but the code did not solve the system
Assuming that Mathematica can solve my system, how could I write code using NDSolve
or any other method?