How to code this iterative process?

Question

I am facing difficulty to realize this double iterative process.

The equations in question are

The flow chart for the iterative process is given as

The different parameters are defined as

alphan=1.72*10^(-4);
alphap=2.037*10^(-4);
L=1.3*10^(-3);
A=2.08*10^(-6);
kp=1.265;
kn=1.011;
sigmap=1.314e-5;
sigman=1.119e-5;
alphapn=alphap-alphan;
Rpn=L/(A)*(sigmap+sigman);
RL=1;
Kpn=(A/L)*(kp+kn);
cf=4205;
cc=4153;
hf=80;
hc=1000;
Tfin=773;
Tcin=353;
mf=20;
mc=20;

The equations are

qh[i, j] = cf*mf*(Tf[i, j] - Tf[i + 1, j])/ny

qh[i, j] = hf*Sf[i, j]*(Tfav[i, j] - Th[i + 1, j])

qh[i, j] = alphapn*I1*Th[i, j] + Kpn*(Th[i, j] - TL[i, j]) - 0.5*I1^2*Rpn

qL[i, j] = cc*mc*(Tc[i + 1, j] - Tc[i + 1, j])/ny

qL[i, j] = hc*Sc[i, j]*(TL[i, j] - Tcav[i, j])

qL[i, j] = alphapn*I1*TL[i, j] + Kpn*(Th[i, j] - TL[i, j]) + 0.5*I1^2*Rpn

I1 = Sum[alphapn*(Th[i, j] - TL[i, j]), {i, 1, nx}, {j, 1, nx}]/(nx*ny*Rpn + RL)

P = Sum[alphapn*(qh[i, j] - qL[i, j]), {i, 1, nx}, {j, 1, nx}]

eta = 100*P/Sum[qh[i, j], {i, 1, nx}, {j, 1, nx}]

Tf[1, j] = Tfin

Tc[1, j] = Tcin

Tfav[i, j] = (Tf[i, j] + Tf[i + 1, j])/2

Tcav[i, j] = (Tc[i, j] + Tc[i + 1, j])/2

Answer 1

We can organize computation in 3 steps. First, we calculate transition matrix tij for every temperature required as follows

Tfav[i, j] = (Tf[i, j] + Tf[i + 1, j])/2;
Tcav[i, j] = (Tc[i, j] + Tc[i + 1, j])/2;

qh1 = cf*mf*(Tf[i, j] - Tf[i + 1, j])/ny;
qh2 = hf*Sf[i, j]*(Tfav[i, j] - Th[i, j]);
qh3 = alphapn*I1*Th[i, j] + Kpn*(Th[i, j] - TL[i, j]) - 0.5*I1^2*Rpn;
qL1 = cc*mc*(Tc[i + 1, j] - Tc[i, j])/ny;
qL2 = hc*Sc[i, j]*(TL[i, j] - Tcav[i, j]);
qL3 = alphapn*I1*TL[i, j] + Kpn*(Th[i, j] - TL[i, j]) + 0.5*I1^2*Rpn;

sol = Solve[{qh1 == qh2, qh2 == qh3, qL1 == qL2, 
    qL2 == qL3}, {Tf[i + 1, j], Tc[i + 1, j], Th[i, j], TL[i, j]}][[1]]

tij = {Tf[i + 1, j], Tc[i + 1, j], Th[i, j], TL[i, j]} /. sol

On the second step we use tij to compute temperature with I1=0. We don't know how surface areas Sc[i,j], Sf[i,j] can be defined, but in this code its are given constants

alphan = -1.72*10^(-4);
alphap = 2.037*10^(-4);
L = 1.3*10^(-3);
A = 2.08*10^(-6);
kp = 1.265;
kn = 1.011;
sigmap = 1.314 10^-5;
sigman = 1.119 10^-5;
alphapn = alphap - alphan;
Rpn = L/(A)*(sigmap + sigman);
RL = 1;
Kpn = (A/L)*(kp + kn);
cf = 4205;
cc = 4153;
hf = 80;
hc = 1000;
Tfin = 773;
Tcin = 353;
mf = 20;
mc = 200; nx = 20; ny = 10; a = 0.9; b = 0.49;

tf = ConstantArray[Tfin, {nx, ny}];
tc = ConstantArray[Tcin, {nx, ny}]; sf = 
 ConstantArray[a/nx b/ny, {nx, ny}]; sc = sf;

Do[th[i, j] = 
  tij[[3]] /. {i -> i, j -> j, Tc[i, j] -> tc[[i, j]], 
    Tf[i, j] -> tf[[i, j]], Sf[i, j] -> sf[[i, j]], 
    Sc[i, j] -> sc[[i, j]], I1 -> 0}; 
 tL[i, j] = 
  tij[[4]] /. {i -> i, j -> j, Tc[i, j] -> tc[[i, j]], 
    Tf[i, j] -> tf[[i, j]], Sf[i, j] -> sf[[i, j]], 
    Sc[i, j] -> sc[[i, j]], I1 -> 0};, {i, nx}, {j, ny}]

I0 = Sum[alphapn*(th[i, j] - tL[i, j]), {i, 1, nx}, {j, 1, 
    ny}]/(nx*ny*Rpn + RL)

Here we have out 7.63885. On the last step we organize iterations up to state where $I$ converges:

ii[0] = I0; Do[
 Do[Do[tf[[i + 1, j]] = 
     tij[[1]] /. {i -> i, j -> j, Tc[i, j] -> tc[[i, j]], 
       Tf[i, j] -> tf[[i, j]], Sf[i, j] -> sf[[i, j]], 
       Sc[i, j] -> sc[[i, j]], I1 -> ii[k - 1]}; 
    tc[[i + 1, j]] = 
     tij[[2]] /. {i -> i, j -> j, Tc[i, j] -> tc[[i, j]], 
       Tf[i, j] -> tf[[i, j]], Sf[i, j] -> sf[[i, j]], 
       Sc[i, j] -> sc[[i, j]], I1 -> ii[k - 1]};, {i, 1, 
     nx - 1}];, {j, ny}];
 Do[th[i, j] = 
   tij[[3]] /. {i -> i, j -> j, Tc[i, j] -> tc[[i, j]], 
     Tf[i, j] -> tf[[i, j]], Sf[i, j] -> sf[[i, j]], 
     Sc[i, j] -> sc[[i, j]], I1 -> ii[k - 1]}; 
  tL[i, j] = 
   tij[[4]] /. {i -> i, j -> j, Tc[i, j] -> tc[[i, j]], 
     Tf[i, j] -> tf[[i, j]], Sf[i, j] -> sf[[i, j]], 
     Sc[i, j] -> sc[[i, j]], I1 -> ii[k - 1]};, {i, nx}, {j, ny}];
 ii[k] = Sum[
    alphapn*(th[i, j] - tL[i, j]), {i, 1, nx}, {j, 1, 
     ny}]/(nx*ny*Rpn + RL);, {k, 1, 10}] 

We can plot $I-I0$ on every step to check convergence

ListLinePlot[Table[ii[k] - I0, {k, 0, 10}], PlotRange -> All]

Finally we calculate

P = 
  Sum[alphapn*(qh3 - qL3) /. {i -> i, j -> j, Th[i, j] -> th[i, j], 
     TL[i, j] -> tL[i, j], Sf[i, j] -> sf[[i, j]], 
     Sc[i, j] -> sc[[i, j]], I1 -> I0}, {i, 1, nx}, {j, 1, ny}];
eta = 100*
   P/Sum[qh3 /. {i -> i, j -> j, Th[i, j] -> th[i, j], 
       TL[i, j] -> tL[i, j], Sf[i, j] -> sf[[i, j]], 
       Sc[i, j] -> sc[[i, j]], I1 -> I0}, {i, 1, nx}, {j, 1, ny}];

{P, eta}

Out[]= {0.0208876, 0.00331435}