3
$\begingroup$

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?

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ C and D are predefined functions in Mathematica and can't be used as your own functions without grief. Perhaps rename those CapitalC and CapitalD. Next, NDSolve expects every variable and function you aren't solving for to have a predefined, preferably numeric value, so u,a,r1,b,r2,r3,d don't have values that the reader can use to try in his own worksheet and see if his potential answer actually works. Perhaps this will give you enough correct the first wave of things in your code. $\endgroup$
    – Bill
    Dec 30, 2018 at 19:13
  • $\begingroup$ Thanks Bill! I would write the corresponding functions of r1,r2,r2 but the problem is that those functions are related to other in "waterfall" method. If so, I should write all codec here and I think it's not the way to ask. Those functions can be whatever only if they depends of T[z]. For the rest, "u" is equal to 0.005, "q" is equal to 32, "a"=1442, "b"=124, "c"=100 and "d"=-65. $\endgroup$
    – atavistic mephit
    Dec 30, 2018 at 19:32

3 Answers 3

Reset to default
1
$\begingroup$

I changed the code a bit, but the final system of equations does not differ from that obtained by the author's code. To solve the system we use the method of the false transient.

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;

Table[
    M1[i, j] = 2*(1/M[i] + 1/M[j])^-1, {i, 1, Nc}, {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;

Do[
    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]);, {i, 1, Nc}];
y[1, z_] := yeb[z]; y[2, z_] := yw; y[3, z_] := yst[z]; 
y[4, z_] := yst[z];

Do[
    PHI[i, i][z_] := 1;
    Do[
      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;, {i, 1, Nc}], {j, 1, 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;

Table[
   Dc[j] = (0.00143*T[z]^1.75)/(Pt *
      M1[1, j]*(Da[1]^(1/3) + Da[j]^(1/3))^2), {j, 1, 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[2])*(y[2, z] - bb/aa y[1, z]) + ((1/Dc[3])*
          y[3, z] + 
         cc/aa y[1, z]) + ((1/Dc[4])*(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[
    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]), {i, 1, 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; 

Table[ kc[i] = A[i]*Exp[-Ea[i]/(Ru/1000 T[z])], {i, 1, 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]*
     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]*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]*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; t0 = 2;
eqns = { Dm[z]*Ceb''[z] - u*Ceb'[z] - 1442*r1[z], 
    km[z]*T''[z] - q*Cpm[z]*T'[z] - 
     1442*(r1[z]*(-CDH1) + r2[z]*(-CDH2) + r3[z]*(-CDH3))} /. {T[z] ->
      U[t, z], Ceb[z] -> W[t, z], T'[z] -> D[U[t, z], z], 
    T''[z] -> D[U[t, z], z, z], Ceb'[z] -> D[W[t, z], z], 
    Ceb''[z] -> D[W[t, z], z, z]};
bc = {U[t, 0] == 893.15, Derivative[0, 1][U][t, L] == 0, 
   W[t, 0] == 0.992, Derivative[0, 1][W][t, L] == 0 };
ic = {U[0, z] == 893.15, W[0, z] == 0.992}; L = 6; t0 = 2;
sol = NDSolve[{eqns == {D[W[t, z], t], D[U[t, z], t]}, bc, ic}, {U, 
    W}, {t, 0, t0}, {z, 0, L}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 80, "MaxPoints" -> 100, 
       "DifferenceOrder" -> 2}}, MaxSteps -> 10^6];
{Plot3D[Evaluate[U[t, z] /. sol], {t, 0, t0}, {z, 0, L}], 
 Plot3D[Evaluate[W[t, z] /. sol], {t, 0, t0}, {z, 0, L}],
 Plot[Evaluate[U[t0, z] /. sol], {z, 0, L}, AxesLabel -> {"z", "T"}], 
 Plot[Evaluate[W[t0, z] /. sol], {z, 0, L}, 
  AxesLabel -> {"z", "Ceb"}], 
 Plot[Evaluate[U[t0, z] /. sol], {z, 0.95*L, L}, 
  AxesLabel -> {"z", "T"}], 
 Plot[Evaluate[W[t0, z] /. sol], {z, 0.95*L, L}, 
  AxesLabel -> {"z", "Ceb"}]}

fig1

Improve this answer
$\endgroup$
0
$\begingroup$

Compare this character by character with what you have. There is a missing ] in your code that I guessed where to insert, there appear to be differences in capitalization and naming, etc. Maybe this will help you get started. Please run this code as is in a fresh empty notebook and verify if it works for you. Then make small changes, one at a time, trying to bring this closer to your actual problem and verify that it still works after each tiny change, finding and fixing the source of any errors that it displays before making the next tiny change. If you can provide corrections then I will try to edit those into this

L=9;u=0.005;a=1442;b=124;c=100;d=65;p=5;q=32;
CapitalD[z_]:=9;r2[z_]:=11;r3[z_]:=13;
k[z_]:=17;kc1[z_]:=19;keb[z_]:=23;
keq[z_]:=29;kh2[z_]:=31;kst[z_]:=37;
y1[z_]:=41;y3[z_]:=43;y4[z_]:=47;
r1[z_]:=(kc1[z]*keb[z]*(y1[z]*p-(y3[z]*p+y4[z]*p)/keq[z]*p))/
  (1+keb[z]*y1[z*p]+kh2[z]*y4[z]*p+kst[z]*y3[z]*p)^2;
eqn1 = CapitalD[z]*CapitalC''[z]-u*CapitalC'[z]==a*r1[z];
eqn2 = k[z]*T''[z]-q*T'[z]==a*(r1[z]*b+r2[z]*c+r3[z]*CapitalD[z]);
initConds={T[0]==893.15, T'[L]==0,CapitalC[0]==0.992,CapitalC'[L]==0};
s = NDSolve[{eqn1, eqn2, initConds}, {T[z], CapitalC[z]}, {z, 0, 9}][[1]];
Plot[Evaluate[CapitalC[z]/.s],{z,0,9}]
Plot[Evaluate[T[z]/.s],{z,0,9}]
Improve this answer
$\endgroup$
10
  • $\begingroup$ Okay, will rewrite and post the code as soon as possible. $\endgroup$
    – atavistic mephit
    Dec 31, 2018 at 1:33
  • $\begingroup$ Bill! I used your code in my worksheet and well, it plot something. At the moment this plot has no sense, but at least it plot something. By the way, in your code are lot of constants that should be variables as r2[z], r3[z], y1[z] etc.. $\endgroup$
    – atavistic mephit
    Dec 31, 2018 at 9:45
  • $\begingroup$ Okay, using your code changing all constants to his respective functions, the program give me the next output: NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}. $\endgroup$
    – atavistic mephit
    Dec 31, 2018 at 10:42
  • $\begingroup$ as NDSolve is not able to solve this system, Do u know any other method capable of it? $\endgroup$
    – atavistic mephit
    Dec 31, 2018 at 15:16
  • $\begingroup$ okay, srry then! Didn't want to bother you. Thanks for the help anyways! $\endgroup$
    – atavistic mephit
    Dec 31, 2018 at 17:41
0
$\begingroup$

Assuming that previously all the needed functions and constant are defined, you can try

kc1[z_] := 1 + z
keb[z_] := 1 + z
y1[z_] := 1 + z
y3[z_] := 1 + z
y4[z_] := 1 + z
keq[z_] := 1 + z
kh2[z_] := 1 + z
kst[z_] := 1 + z
d[z_] := 1 + z
k[z_] := 1 + z
r1[z_] := (kc1[z]*
 keb[z]*(y1[z]*p - (y3[z]*p + y4[z]*p)/keq[z]*p))/(1 + 
  keb[z]*y1[z]*p + kh2[z]*y4[z]*p + kst[z]*y3[z]*p)^2;
r2[z_] := r1[z]
r3[z_] := r1[z]
eqn1 = d[z]*CC''[z] - u*CC'[z] == a*r1[z];
eqn2 = k[z]*T''[z] - q*T'[z] == a*(r1[z]*b + r2[z]*c + r3[z]*d0);
initConds = {T[0] == 893.15, T'[L] == 0, CC[0] == 0.992, CC'[L] == 0};
parms = {p -> 1, q -> 1, a -> 1, b -> 1, c -> 1, d0 -> 1, u -> 1};
L = 10;
zmax = 9;
delz = 0.01;
eqns = {eqn1, eqn2} /. parms;
DE = Join[eqns, initConds];
s = NDSolve[DE, {T, CC}, {z, 0, zmax, delz}]
Plot[Evaluate[{T[z], CC[z]} /. s], {z, 0, zmax}]
Improve this answer
$\endgroup$
16
  • $\begingroup$ Hello Cesareo, I also used your code before and the result is the same. NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}. $\endgroup$
    – atavistic mephit
    Dec 31, 2018 at 10:49
  • $\begingroup$ @atavisticmephit To avoid misunderstandings in this edit, I defined the declared functions and associated values to the parameters used. I hope this helps. $\endgroup$
    – Cesareo
    Dec 31, 2018 at 11:21
  • $\begingroup$ @atavisticmephit This code as it is, runs without problems at my computer with MATHEMATICA v 11.0 $\endgroup$
    – Cesareo
    Dec 31, 2018 at 11:27
  • $\begingroup$ I going to try it, let's see if it works but I don't understant those parameters. Why all of them are 1? $\endgroup$
    – atavistic mephit
    Dec 31, 2018 at 11:28
  • $\begingroup$ Well, using it i have a solution. The program draws a plot, with no sense, but a plot. There are things that I don't understant as {p -> 1, 33.66 -> 1, a -> 1, b -> 1, c -> 1, d0 -> 1, 19.53 -> 1} ( output message). $\endgroup$
    – atavistic mephit
    Dec 31, 2018 at 11:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.