How to code this iterative process?

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
• What difficulty are you facing specifically? What have you tried so far? Commented Jan 9, 2021 at 2:44
• @MarcoB I have no clue from where to start...otherwise I would have posted my try..
– zhk
Commented Jan 9, 2021 at 2:45
• You will need to provide more information on the problem you are trying to solve, at the very least, with some context for the “instructions” you posted. Also translations of the formulae in MMA format would be nice. It still seems unlikely that somebody will take it upon themselves to do your task for you from scratch... Commented Jan 9, 2021 at 4:19
• The biggest problem with your post, I think, is that you do not say what do you want to find (unknowns) and what are the starting values of the iterative process. Commented Jan 14, 2021 at 21:14
• @zhk Could you refer paper or book where you taken this model? Commented Jan 15, 2021 at 13:24

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}
• First of thank you for your efforts. How to plot P vs Tfin, P vs mf and P vs Sf? This will show whether we are getting the correct results or not?
– zhk
Commented Jan 17, 2021 at 2:15
• @zhk To get result as in the paper we should use same input data. But I didn't see definition of Sc[i,j], Sf[i,j] in the paper and in your code. Also your question is not about P vs Tfin, P vs mf, or P vs Sf, but about code of iterative process. Could you add this extension to your post with some picture what exactly you asking about? Commented Jan 17, 2021 at 10:36
• I agree with you about the incomplete info from the paper. But by plotting P vs Tfin, P vs mf, or P vs Sf or eta vs Tfin or eta vs mf will enable us to interpret the results to the best of our understanding.
– zhk
Commented Jan 22, 2021 at 3:10
• @zhk Ok! Have you any version of Sc, Sf? Commented Jan 22, 2021 at 14:00
• When can take a fix value for it say Sc=Sf=0.5.
– zhk
Commented Jan 23, 2021 at 2:00