Thermodynamics & heat transfer: MQL recipes

Recipes

U-value and heat loss of an insulated wall
Carnot limit of a heat pump
Isothermal and adiabatic compression of air
Newton's law of cooling: fit, then predict
Efficiency of an annular fin
Frost depth for a water main (semi-infinite solid)
Aluminium pin fin: heat rate, efficiency, tip temperature
Radiation and convection loss of a steam pipe
Counterflow exchanger: effectiveness-NTU checked by LMTD
Calibrating thermocouples at the ice point

executed on MentDB

U-value and heat loss of an insulated wall

A 12 m2 wall is made of 20 cm brick, 8 cm mineral wool and 1.5 cm plaster. Add the layer resistances e/lambda and the surface resistances, then get the U-value, the heat loss for 20 degC inside / -5 degC outside, and the energy over a 2800 degree-day heating season.

Pipeline

01 · layersSum of e/lambdaunit calc→02 · uU = 1 / R totalmath calc→03 · lossHeat flow U A dTunit calc→04 · seasonEnergy in kWhunit calc

Script

admin

#1. Thermal resistance of each layer R = e/lambda, per m2: brick, mineral wool, plaster;
-> "[R]" (unit calc "0.20 m/(0.72 W/(m*K)) + 0.08 m/(0.035 W/(m*K)) + 0.015 m/(0.50 W/(m*K))" "m^2*K/W");
#2. Add the surface resistances Rsi = 0.13 and Rse = 0.04;
-> "[Rt]" (math calc "[R]+0.13+0.04");
-> "[U]" (math calc "1/[Rt]");
#3. Heat flow through 12 m2 for 20 degC inside, -5 degC outside;
-> "[P]" (unit calc (concat [U] " W/(m^2*K) * 12 m^2 * 25 K") "W");
#4. Energy over a heating season of 2800 degree-days;
concat "R total = " (math round [Rt] 3) " m2K/W, U = " (math round [U] 3) " W/m2K, loss = " (math round [P] 1) " W, season = " (math round (unit calc (concat [U] " W/(m^2*K) * 12 m^2 * 2800 K * 24 h") "kWh") 0) " kWh"

Result

mentdb

R total = 2.763 m2K/W, U = 0.362 W/m2K, loss = 108.6 W, season = 292 kWh

Go further

Find the insulation thickness for U = 0.2 W/m2K with num solve
Compute the temperature at each interface from the series resistances

executed on MentDB

Carnot limit of a heat pump

A heat pump feeds a floor heating at 35 degC from outdoor air at -5 degC. Compute the Carnot COP, a realistic COP at 50 % of Carnot, and the electricity needed for 9000 kWh of heat per year.

Pipeline

01 · kelvinTemperatures in Kunit convert→02 · carnotCOP = Th / (Th - Tc)math calc→03 · energyElectricity per yearmath round

Script

admin

#1. Floor heating at 35 degC, outdoor air at -5 degC, in kelvin;
-> "[Th]" (unit convert 35 "degC" "K");
-> "[Tc]" (unit convert -5 "degC" "K");
#2. Carnot COP of a heat pump Th / Th - Tc, real machine at 50 % of Carnot;
-> "[cop]" (math calc "[Th]/([Th]-[Tc])");
-> "[real]" (math calc "0.5*[cop]");
#3. Electricity for 9000 kWh of heat;
concat "COP Carnot = " (math round [cop] 2) ", real COP = " (math round [real] 2) ", electricity = " (math round (math calc "9000/[real]") 0) " kWh/year instead of 9000 kWh with a resistive heater"

Result

mentdb

COP Carnot = 7.7, real COP = 3.85, electricity = 2337 kWh/year instead of 9000 kWh with a resistive heater

Go further

Tabulate the COP against the outdoor temperature with num table
Compare with radiators at 55 degC

executed on MentDB

Isothermal and adiabatic compression of air

2 L of air at 1 bar and 20 degC are compressed to 0.5 L. Get the amount of gas from pV = nRT, the isothermal work as the integral of p dV (checked against nRT ln 4), and the final temperature if the compression is adiabatic (gamma = 1.4).

Pipeline

01 · statepV and n = pV / RTunit calc→02 · workW = integral of nRT/V dVnum integrate→03 · adiabaticT2 = T1 (V1/V2)^(gamma-1)math calc→04 · celsiusT2 in degCunit convert

Script

admin

#1. 2 L of air at 1 bar and 20 degC;
-> "[nRT]" (unit calc "1 bar * 2 L" "J");
-> "[n]" (unit calc (concat [nRT] " J/(_R*293.15 K)") "mol");
#2. Isothermal compression to 0.5 L: W = integral of nRT/V dV;
-> "[W]" (num integrate (concat [nRT] "/v") "v" 0.5 2);
#3. Adiabatic compression instead, gamma = 1.4: T2 = T1 * V1/V2 ^ gamma-1;
-> "[T2]" (math calc "293.15*(2/0.5)^0.4");
concat "n = " (math round [n] 5) " mol, isothermal work = " (math round [W] 2) " J (nRT ln4 = " (math round (math calc "[nRT]*ln(4)") 2) " J); adiabatic: T2 = " (math round (unit convert [T2] "K" "degC") 1) " degC"

Result

mentdb

n = 0.08206 mol, isothermal work = 277.26 J (nRT ln4 = 277.26 J); adiabatic: T2 = 237.3 degC

Go further

Compute the adiabatic work with num integrate of p1*(V1/V)^1.4
Use a van der Waals gas and compare

executed on MentDB

Newton's law of cooling: fit, then predict

From temperature readings of a coffee in a 20 degC room, fit T = 20 + (T0 - 20) exp(-k t) to get k with its uncertainty, then integrate dT/dt = -k (T - 20) until the coffee reaches 60 degC.

Pipeline

01 · fitNon-linear fit of kstat fit_model→02 · uncertaintyk with its uncertaintymath round_uncertainty→03 · odeIntegrate until T = 60num ode

Script

admin

#1. Coffee cooling measurements in a 20 degC room, t in min;
-> "[t]" "[0, 5, 10, 15, 20, 30]";
-> "[T]" "[90.1, 80.0, 72.1, 64.4, 58.6, 48.3]";
#2. Fit Newton's law T = 20 + T0-20 * exp of -k t;
json load "rx_cool" (stat fit_model "20+(T0-20)*exp(-k*t)" "t" [t] [T] "{\"T0\": 90, \"k\": 0.05}");
-> "[k]" (json select "rx_cool" /parameters/k);
#3. Solve the ODE dT/dt = -k T - 20 until the coffee reaches 60 degC;
json load "rx_ode" (num ode "{\"T\": \"-[k]*(T-20)\"}" "{\"T\": 90}" 0 120 "{\"points\": 3, \"event\": \"T-60\"}");
concat "k = " (math round_uncertainty [k] (json select "rx_cool" /errors/k)) " 1/min, r2 = " (math round (json select "rx_cool" /r2) 4) "; 60 degC reached after " (math round (json select "rx_ode" /event/t) 1) " min"

Result

mentdb

k = 0.03007 ± 0.00024 1/min, r2 = 0.9998; 60 degC reached after 18.6 min

Go further

Plot the data and the fitted curve with plot fit
Add the room temperature as a third fitted parameter

executed on MentDB

Efficiency of an annular fin

An aluminium circular fin on a 25 mm tube: the exact efficiency uses the modified Bessel functions I0, I1, K0, K1. The script compares it with a straight fin of the same length and computes the heat removed by one fin.

Pipeline

01 · mFin parameter mmath calc→02 · etaBessel efficiency formulabessel_i(nu, x)→03 · compareStraight fin tanh(mL)/mLmath calc→04 · heatHeat per finmath calc

Script

admin

#1. Aluminium annular fin: k = 200 W/m.K, h = 40 W/m2.K, t = 1 mm, r1 = 12.5 mm, r2 = 30 mm;
-> "[k]" 200; -> "[h]" 40; -> "[t]" 0.001; -> "[r1]" 0.0125;
-> "[r2c]" (math calc "0.030+[t]/2");
-> "[m]" (math calc "sqrt(2*[h]/([k]*[t]))");
#2. Fin efficiency with the modified Bessel functions I and K;
-> "[eta]" (math calc "2*[r1]/([m]*([r2c]^2-[r1]^2))*(bessel_k(1,[m]*[r1])*bessel_i(1,[m]*[r2c])-bessel_i(1,[m]*[r1])*bessel_k(1,[m]*[r2c]))/(bessel_i(0,[m]*[r1])*bessel_k(1,[m]*[r2c])+bessel_k(0,[m]*[r1])*bessel_i(1,[m]*[r2c]))");
#3. Straight fin of the same length, for comparison;
-> "[eta_s]" (math calc "tanh([m]*([r2c]-[r1]))/([m]*([r2c]-[r1]))");
#4. Heat from one fin, base 100 degC, air 25 degC;
-> "[q]" (math calc "[eta]*[h]*2*pi*([r2c]^2-[r1]^2)*(100-25)");
concat "m = " (math round [m] 2) " 1/m, annular fin efficiency = " (math round [eta] 4) " (straight fin " (math round [eta_s] 4) "), heat per fin = " (math round [q] 2) " W"

Result

mentdb

m = 20 1/m, annular fin efficiency = 0.937 (straight fin 0.9589), heat per fin = 13.67 W

Go further

Find the outer radius that gives an efficiency of 0.9 with num solve.
bessel_k(nu, x) is also the kernel of the steady heat loss of a buried pipe.

executed on MentDB

Frost depth for a water main (semi-infinite solid)

Soil at 20 C (alpha = 0.138E-6 m^2/s) sees its surface drop to -15 C for 60 days. With the semi-infinite solid solution (T - Ts)/(Ti - Ts) = erf(x / (2 sqrt(alpha t))), find the depth reached by 0 C, the minimum burial depth of a water main, and the temperature at 1 m.

Pipeline

01 · time60 days in secondsunit convert→02 · depthInverse error functionmath erf_inv→03 · profileTemperature at 1 mmath erf→04 · roundDepth to the next 10 cmmath ceil

Script

admin

-> "[a]" 0.138E-6;
-> "[t]" (unit convert 60 "d" "s");
-> "[Ti]" 20;
-> "[Ts]" -15;
#1. Diffusion length 2 sqrt of alpha t;
-> "[L]" (math calc "2*sqrt([a]*[t])");
#2. Depth where T = 0;
-> "[x0]" (math calc (concat "[L]*" (math erf_inv (math calc "(0-[Ts])/([Ti]-[Ts])"))));
#3. Temperature at 1 m;
-> "[T1]" (math calc (concat "[Ts]+([Ti]-[Ts])*" (math erf (math calc "1/[L]"))));
concat "2 sqrt(alpha t) = " (math round [L] 3) " m; 0 C reached at " (math round [x0] 3) " m, bury the main at least " (math calc (concat (math ceil (math calc "[x0]*10")) "/10")) " m; T at 1 m = " (math round [T1] 2) " C"

Result

mentdb

2 sqrt(alpha t) = 1.692 m; 0 C reached at 0.677 m, bury the main at least 0.7 m; T at 1 m = 5.89 C

Go further

The penetration depth grows like sqrt(t): 4 times longer gives only twice deeper.
With a convective surface condition use the erfc-exp solution.

executed on MentDB

Aluminium pin fin: heat rate, efficiency, tip temperature

A pin fin D = 5 mm, L = 50 mm, k = 200 W/m.K, in air h = 25 W/m^2.K, base 80 K above the air. m = sqrt(4h/(kD)); with the corrected length Lc = L + D/4: q = sqrt(hPkA) theta_b tanh(m Lc), efficiency tanh(m Lc)/(m Lc), and the adiabatic-tip excess temperature theta_b/cosh(m L).

Pipeline

01 · geometryPerimeter and section with piconstant math_pi→02 · mFin parametermath sqrt→03 · qHeat rate and efficiencymath tanh→04 · tipTip temperaturemath cosh

Script

admin

-> "[pi]" (string regexp_extract (constant math_pi) "[-0-9.]+E?[-0-9]*" 0);
-> "[D]" 0.005;
-> "[L]" 0.05;
-> "[k]" 200;
-> "[h]" 25;
-> "[tb]" 80;
#1. Geometry and fin parameter;
-> "[P]" (math calc "[pi]*[D]");
-> "[A]" (math calc "[pi]*[D]^2/4");
-> "[m]" (math sqrt (math calc "[h]*[P]/([k]*[A])"));
-> "[mLc]" (math calc "[m]*([L]+[D]/4)");
#2. Heat rate and efficiency;
-> "[th]" (math tanh [mLc]);
-> "[q]" (math calc "sqrt([h]*[P]*[k]*[A])*[tb]*[th]");
concat "m = " (math round [m] 3) " 1/m, q = " (math round [q] 3) " W, efficiency = " (math round (math calc "[th]/[mLc]") 4) ", effectiveness = " (math round (math calc "[q]/([h]*[A]*[tb])") 1) "; tip excess " (math round (math calc (concat "[tb]/" (math cosh (math calc "[m]*[L]")))) 2) " K"

Result

mentdb

m = 10 1/m, q = 1.482 W, efficiency = 0.9208, effectiveness = 37.8; tip excess 70.95 K

Go further

An effectiveness well above 2 justifies the fin.
For an array: overall surface efficiency 1 - N Af/At (1 - eta_f).

executed on MentDB

Radiation and convection loss of a steam pipe

An uninsulated steam pipe, D = 70 mm, surface 200 C, emissivity 0.8, in a room at 25 C with h = 15 W/m^2.K. Heat loss per metre by radiation (Stefan-Boltzmann) and by convection, and the linearised radiation coefficient h_r = eps sigma (Ts + Tsur)(Ts^2 + Tsur^2).

Pipeline

01 · offsetCelsius offset from the unit listunit list→02 · sigmaStefan-Boltzmann constantconstant physics_stefan_boltzmann→03 · radT^4 termsmath power→04 · totalLoss per metremath calc

Script

admin

#1. Kelvin offset of degC read in the unit table;
-> "[off]" (string regexp_extract (unit list "temperature") "degC.*?offset_si\": ([0-9.]+)" 1);
-> "[Ts]" (math calc "200+[off]");
-> "[Tr]" (math calc "25+[off]");
-> "[sg]" (string regexp_extract (constant physics_stefan_boltzmann) "[-0-9.]+E?[-0-9]*" 0);
-> "[A]" (math calc "pi*0.07");
#2. Radiation and convection per metre;
-> "[qr]" (math calc (concat "0.8*[sg]*[A]*(" (math power [Ts] 4) "-" (math power [Tr] 4) ")"));
-> "[qc]" (math calc "15*[A]*([Ts]-[Tr])");
#3. Equivalent radiation coefficient;
-> "[hr]" (math calc "0.8*[sg]*([Ts]+[Tr])*([Ts]^2+[Tr]^2)");
concat "offset " [off] " K, sigma = " [sg] "; radiation " (math round [qr] 0) " W/m, convection " (math round [qc] 0) " W/m, total " (math round (math calc "[qr]+[qc]") 0) " W/m; h_r = " (math round [hr] 2) " W/m2.K vs h = 15"

Result

mentdb

offset 273.15 K, sigma = 5.6703744191844314E-8; radiation 421 W/m, convection 577 W/m, total 998 W/m; h_r = 10.94 W/m2.K vs h = 15

Go further

At 200 C radiation is as large as natural convection: never neglect it.
Add an insulation layer and solve the surface temperature with num solve.

executed on MentDB

Counterflow exchanger: effectiveness-NTU checked by LMTD

Oil (Ch = 2000 W/K, 120 C) cools against water (Cc = 4000 W/K, 20 C) in a counterflow exchanger of UA = 3000 W/K. Effectiveness eps = (1 - e^(-NTU(1-Cr)))/(1 - Cr e^(-NTU(1-Cr))) computed with expm1, outlet temperatures, and the LMTD check UA = q / LMTD.

Pipeline

01 · ntuNTU and capacity ratiomath calc→02 · epsEffectiveness with expm1math expm1→03 · lmtdLog mean temperature differencemath log→04 · inverseNTU back from eps with log1pmath log1p

Script

admin

-> "[Cmin]" 2000;
-> "[Cr]" 0.5;
-> "[NTU]" (math calc "3000/[Cmin]");
#1. Effectiveness, -expm1 = 1 - exp;
-> "[em]" (math expm1 (math calc "-[NTU]*(1-[Cr])"));
-> "[eps]" (math calc "-[em]/(1-[Cr]*([em]+1))");
#2. Heat duty and outlets;
-> "[q]" (math calc "[eps]*[Cmin]*(120-20)");
-> "[Tho]" (math calc "120-[q]/2000");
-> "[Tco]" (math calc "20+[q]/4000");
#3. LMTD and back to UA;
-> "[d1]" (math calc "120-[Tco]");
-> "[d2]" (math calc "[Tho]-20");
-> "[lmtd]" (math calc (concat "([d1]-[d2])/" (math log (math calc "[d1]/[d2]"))));
#4. NTU from eps: ln of eps-1 over eps.Cr-1, over Cr-1;
-> "[ntu2]" (math calc (concat (math log1p (math calc "([eps]-1)/([eps]*[Cr]-1)-1")) "/([Cr]-1)"));
concat "NTU = " [NTU] ", eps = " (math round [eps] 4) ", q = " (math round [q] 0) " W, oil out " (math round [Tho] 2) " C, water out " (math round [Tco] 2) " C; LMTD = " (math round [lmtd] 3) " K, q/LMTD = " (math round (math calc "[q]/[lmtd]") 1) " W/K; NTU back = " (math round [ntu2] 6)

Result

mentdb

NTU = 1.5, eps = 0.6908, q = 138157 W, oil out 50.92 C, water out 54.54 C; LMTD = 46.052 K, q/LMTD = 3000 W/K; NTU back = 1.5

Go further

For Cr = 1 the formula is 0/0: use eps = NTU/(1 + NTU).
Size the exchanger for a target outlet: solve for NTU then UA.

executed on MentDB

Calibrating thermocouples at the ice point

Two type-K thermocouples read an ice bath (0 C). The reference indicator has a known noise of 0.10 K. A z-test tells whether sensor A has an offset, the population variance gives its noise, and a Mann-Whitney test compares A and B without assuming normality.

Pipeline

01 · zOffset test, known sigmastat ztest→02 · noisePopulation variancestat pvariance→03 · compareNon-parametric A vs Bstat mann_whitney→04 · offsetCorrection to applymath avg

Script

admin

-> "[A]" "[0.18, 0.05, 0.22, 0.11, 0.15, 0.09, 0.20, 0.13]";
-> "[B]" "[-0.04, 0.02, -0.08, 0.05, -0.01, 0.03, -0.06]";
#1. Offset of A against 0 C;
-> "[z]" (stat ztest [A] 0 0.10);
#2. Noise of A;
-> "[sd]" (math sqrt (stat pvariance [A]));
#3. A vs B, two-sided;
-> "[mw]" (stat mann_whitney [A] [B]);
concat "z-test A: " [z] "; correction " (math round (math calc (concat "-" (math avg [A]))) 4) " K, noise " (math round [sd] 4) " K; Mann-Whitney A vs B: " [mw]

Result

mentdb

z-test A: {"z": 3.995153313703993, "p_value": 6.465240129782352E-5}; correction -0.1413 K, noise 0.054 K; Mann-Whitney A vs B: {"u": 55.5, "p_value": 0.0017633794447228074, "method": "asymptotic"}

Go further