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

kJ=1000;
EnergyOfOneBeer=(1850/2)kJ;
m=85;
g=9.81;
Solve[EnergyOfOneBeer==h*m*g]

Quiet@Remove["Global`*"];
cH2O=4186;
ΔT=1;
g=9.81;
m=1;
Solve[m*cH2O*ΔT==h*m*g]

Quiet@Remove["Global`*"];
sol=DSolve[{T''[x]==0,T[0]==T1,T[d]==T2},T[x],x][[1]];
temperature=T[x]/.sol;
Print["T[x]= ",temperature]
q=-λ*Grad[temperature,{x,y,z}]
qx=q[[1]]
Print["qx= ",qx]
P=S*qx;
Print["P= ",P]

(Sqrt[1000.])

Quiet@Remove["Global`*"];
(*convection, wire of 400kV overhead lines *)
Ts=70;
Ta=0;
cm=0.01;
d=5cm;
Dchar=d;
Tdef=((Ts+Ta)/2);
Print["Tdef= ",Tdef," °C"]
λ=29.99*(10^(-3));

(*https://www.engineersedge.com/heat_transfer/nusselt_number_13856.htm*)

v=20;
ν=16.45*(10^(-6));(*konematic viscosity*)
Rey=((v*Dchar)/ν);
Print["Re= ",Rey," -"];
Pr=0.7;
Nu=0.027*(Rey^0.805)*(Pr^(1/3));
Print["Nu= ",Nu," -"];
αForced=Nu*(λ/Dchar);
Print["αForced= ",αForced," W.(m^(-2)).(K^(-1))"]

NuNew[re_]:=0.3+((0.62*(re^0.5)*(Pr^(1/3)))/(((1+((((0.4/Pr)))^(2/3))))^(1/4)))*(((1+((((re/282000)))^(5/8))))^(4/5));
Print["NuNew= ",NuNew[Rey]," -"]

length=1;
S=Pi*d*length;
Pconv=S*αForced*(Ts-Ta);
Print["Power of convection from 1 meter of wire= ",Pconv," W"]

ϵ=0.8;
σ=5.67*(10^(-8));
Prad=S*ϵ*σ*((((Ts+273))^4)-(((Ta+273))^4));
Print["Power of radiation from 1 meter of wire= ",Prad" W"]

PtotalForced=Pconv+Prad;
Print["Total power from 1 meter of wire= ",PtotalForced" W"]

ρAl20=2.65*(10^(-8));
αR=0.00429;
ρAl=ρAl20*(1+αR*(Ts-20));
ScrossSection=((Pi*(d^2))/4);
R=((ρAl*length)/ScrossSection);
Print["R1m= ",(10^6)R," μΩ"]
electricalLoses=R*(current^2);
Solve[electricalLoses==PtotalForced]

(*start of free convection*)
Series[(1/(1+β*ΔT)),{ΔT,0,3}]

β=(1/(Tdef+273.));
g=9.81;
Gr=((β*(Ts-Ta)*g*(Dchar^3))/(ν^2));
Print["Gr= ",Gr," -"];
Print["Gr*Pr= ",Gr*Pr," -"]
NuFree=0.48*(((Gr*Pr))^0.25);
Print["NuFree= ",NuFree," -"]

αFree=NuFree*(λ/Dchar);
Print["αFree= ",αFree," W.(m^(-2)).(K^(-1))"]

PconvFree=S*αFree*(Ts-Ta);
Print["Power of free convection from 1 meter of wire= ",PconvFree," W"]

PtotalFree=PconvFree+Prad;
Print["Power of radiation from 1 meter of wire= ",Prad" W"]
Print["Total power from 1 meter of wire free= ",PtotalFree" W"]

Solve[electricalLoses==PtotalFree]

Quiet@Remove["Global`*"];
(*heater in a room*)
Ts=50;
Ta=20;
a=0.5;
b=1;
Dchar=a;
Tdef=((Ts+Ta)/2);
Print["Tdef= ",Tdef," °C"];
λ=26.99*(10^(-3));
ν=16.45*(10^(-6));(*konematic viscosity*)
β=(1/(Tdef+273.));
g=9.81;
Gr=((β*(Ts-Ta)*g*(Dchar^3))/(ν^2));
Pr=0.7;
Print["Gr= ",Gr," -"];
Print["Gr*Pr= ",Gr*Pr," -"]

Ra=Gr*Pr;
Nu=0.68+((0.67*(Ra^0.25))/(((1+((((0.492/Pr)))^(9/16))))^(4/9)));
α=Nu*(λ/Dchar);
NuMIT=(((0.825+((0.387*(Ra^(1/6)))/(((1+((((0.492/Pr)))^(9/16))))^(8/27)))))^2);
αMIT=NuMIT*(λ/Dchar);
{α,αMIT}

NuOld=0.135*(((Gr*Pr))^(1/3));
αOld=NuOld*(λ/Dchar);
{α,αMIT,αOld}

S=2*a*b;
Pconv=αMIT*S*(Ts-Ta);
Print["Power of convection= ",Pconv," W"]

ϵ=0.8;
σ=5.67*(10^(-8));
Prad=S*ϵ*σ*((((Ts+273))^4)-(((Ta+273))^4));
Print["Power of radiation from the heater= ",Prad" W"]

Print["Total power from the heater= ",Prad+Pconv," W"]

Quiet@Remove["Global`*"];
(*heat exchangers*)
Plot[(1-((2*(ξ-1))/((ξ+1)*Log[ξ])))*100,{ξ,0.5,2},AxesLabel->{"ξ[-]","δ[%]"},GridLines->Full,PlotStyle->{Red,Thickness[0.006]}];
Limit[((Δ1-Δ2)/(Log[((Δ1)/(Δ2))])),Δ2->Δ1]

(*parallel or counter flow?*)
volume=2.5*5*11;
Print["Volume= ",volume," (m^3)"];
Tinside=20;
Toutside=-12;
cAir=1000;
density=(349./(Tinside+273));
Print["ρAir= ",density," (kg/(m^3))"]
timeOfAirExchange=3*3600;
massOfAir=volume*density;
Print["mAir= ",massOfAir," kg"]
energy=massOfAir*cAir*(Tinside-Toutside);
Print["energy for one exchange= ",(energy/(10^6))," MJ"]
power=(energy/timeOfAirExchange);
Print["power for ventilation= ",power," W"]
mDot=(massOfAir/timeOfAirExchange);
Print["mDot= ",mDot," kg/s"]

(*counter flow*)
T1c=Tinside;
T3c=Toutside;
powerC[kS_]:=Module[{eqns,Δ1,Δ2,P,T4,T2},
Δ1=T1c-T4;
Δ2=T2-T3c;
eqns={P==kS*((Δ2+Δ1)/2),P==mDot*cAir*(T1c-T2),
mDot*cAir*T1c+mDot*cAir*T3c==mDot*cAir*T2+mDot*cAir*T4
};
P/.Solve[eqns][[1]]
];
(*parallel flow*)
powerP[kS_]:=Module[{eqnsLinear,Δ1,Δ2,P,T4,T2,solLin,forFindRoot,eqns},
Δ1=T1c-T3c;
Δ2=T2-T4;
eqnsLinear={P==kS*((Δ2+Δ1)/2),P==mDot*cAir*(T1c-T2),
mDot*cAir*T1c+mDot*cAir*T3c==mDot*cAir*T2+mDot*cAir*T4
};
solLin=Solve[eqnsLinear][[1]];
forFindRoot=solLin/.Rule->List;
eqns={P==kS*((Δ1-Δ2)/(Log[((Δ1)/(Δ2))])),P==mDot*cAir*(T1c-T2),
mDot*cAir*T1c+mDot*cAir*T3c==mDot*cAir*T2+mDot*cAir*T4
};
Quiet[Chop[P/.FindRoot[eqns,forFindRoot]]]
];

dataC=Table[{kS,powerC[kS]},{kS,0.1,100,1}];
dataP=Table[{kS,powerP[kS]},{kS,0.1,100,1}];
plC=ListPlot[dataC,AxesLabel->{"kS [W/K]","P [W]"},Joined->True,PlotStyle->{Red,Thickness[0.006]},GridLines->Full];
plP=ListPlot[dataP,AxesLabel->{"kS [W/K]","P [W]"},Joined->True,PlotStyle->{Blue,Thickness[0.006],Dashed},GridLines->Full,PlotRange->All];
Show[plC,plP,PlotRange->{Automatic,{0,power}}]