Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];

(*parallel low*)

m1=2;
m2=1.5;
c1=4186;
c2=c1;
T1=80;
T3=10;
ClearAll[powerParallel];
powerParallel[kS_]:=Module[{eqn1,eqn2,eqn3,eqn4,Δ1,Δ2,solLinear,solFinal},
eqn1=T1*m1*c1+T3*m2*c2==T2*m1*c1+T4*m2*c2;
eqn2=P==m1*c1*(T1-T2);
eqn3=P==kS*((Δ1-Δ2)/(Log[((Δ1)/(Δ2))]));
eqn4=P==kS*((Δ1+Δ2)/2);
Δ1=T1-T3;
Δ2=T2-T4;
solLinear=Solve[{eqn1,eqn2,eqn4}][[1]];
solFinal=Chop@FindRoot[{eqn1,eqn2,eqn3},{P,P/.solLinear},{T2,T2/.solLinear},{T4,T4/.solLinear}];
{P,T2,T4}/.solFinal

];
powerParallel[10000]

ClearAll[powerCounter];
powerCounter[kS_]:=Module[{eqn1,eqn2,eqn3,eqn4,Δ1,Δ2,solLinear,solFinal},

Δ1=T1-T4;
Δ2=T2-T3;
eqn1=T1*m1*c1+T3*m2*c2==T2*m1*c1+T4*m2*c2;
eqn2=P==m1*c1*(T1-T2);
eqn3=P==kS*((Δ1-Δ2)/(Log[((Δ1)/(Δ2))]));
eqn4=P==kS*((Δ1+Δ2)/2);

solLinear=Solve[{eqn1,eqn2,eqn4}][[1]];
solFinal=Chop@Quiet@FindRoot[{eqn1,eqn2,eqn3},{P,P/.solLinear},{T2,T2/.solLinear},{T4,T4/.solLinear}];
{P,T2,T4}/.solFinal

];
powerCounter[10000];

kSmin=100.;
kSmax=100000.;
nkS=100;
q=((kSmax/kSmin)^(1/nkS));
kses=kSmin*(q^(#-1))&/@Range[nkS+1];

dataCounter={#,powerCounter[#]}&/@kses;
dataParallel={#,powerParallel[#]}&/@kses;

pl1=ListLogLinearPlot[dataCounter/.{kS_,{P_,T2_,T4_}}:>{kS,P},AxesLabel->{"kS [W/K]","P [W]"},PlotStyle->{Red,Thick},Joined->True];
pl2=ListLogLinearPlot[dataParallel/.{kS_,{P_,T2_,T4_}}:>{kS,P},AxesLabel->{"kS [W/K]","P [W]"},PlotStyle->{Blue,Thick},Joined->True];
Show[pl2,pl1,PlotRange->All]

mass=100.*100;
c=4186.;
Tin=5;
Tout=90;
energy=mass*c*(Tout-Tin);
Print["energy= ",(energy/(10^9))," GJ"];
time=8*3600;
power=(energy/time);
Print["power= ",(power/1000)," kW"];
Δ=8;
kS=(power/Δ);
Print["kS= ",(kS/1000)," kW/K"]

nTubes=50;
massFlow=(mass/time);
Print["mass flow= ",massFlow," kg/s"]
density=1000;
volumeFlow=(massFlow/density);
Print["volume flow= ",volumeFlow," (m^3)/s"];
Tdef=50;
ν=5.54*(10^(-7));
mm=0.001;
DcharBeer=10mm;
DsmalExt=12mm;
DbigInt=DsmalExt+2*5mm;
Swater=(1/4)*Pi*((DbigInt^2)-(DsmalExt^2));
DcharWater=((4*Swater)/(Pi*(DbigInt+DsmalExt)));