"""Mars convection at radiator scale (convection module of FUNNY-TM-2026-003). Properties: NIST WebBook (Span-Wagner EOS; Laesecke-Muzny 2017 viscosity; Huber 2016 conductivity) 600 Pa isobar, 216.59-316.59 K, linearly extrapolated below the triple-point temperature (<=16.6 K). """ import math, numpy as np R = 188.92 # J/kgK, CO2 # NIST 600 Pa isobar: T, cp (J/gK), mu (uPa s), k (W/mK) nist = np.array([ [216.59,0.75399,10.951,0.010665],[226.59,0.76527,11.439,0.011343],[236.59,0.77658,11.928,0.012036], [246.59,0.78785,12.415,0.012744],[256.59,0.79905,12.901,0.013465],[266.59,0.81012,13.386,0.014199], [276.59,0.82104,13.870,0.014942],[286.59,0.83178,14.352,0.015696],[296.59,0.84234,14.831,0.016458], [306.59,0.85268,15.308,0.017228],[316.59,0.86282,15.783,0.018005]]) def interp(T, col): x = nist[:,0]; y = nist[:,col] if T < x[0]: s = (y[1]-y[0])/(x[1]-x[0]); return y[0] + s*(T-x[0]) if T > x[-1]: s = (y[-1]-y[-2])/(x[-1]-x[-2]); return y[-1] + s*(T-x[-1]) return float(np.interp(T, x, y)) def props(T, p=600.0): cp = interp(T,1)*1000; mu = interp(T,2)*1e-6; k = interp(T,3) rho = p/(R*T); nu = mu/rho; Pr = mu*cp/k; alpha = k/(rho*cp) return dict(T=T,rho=rho,cp=cp,mu=mu,k=k,nu=nu,Pr=Pr,alpha=alpha) g = 3.71 def natural(L, dT, Tinf=210.0, p=600.0, orient='up'): Tf = Tinf + dT/2; P = props(Tf, p); beta = 1/Tf if orient == 'vert': Lc = L Ra = g*beta*dT*Lc**3/(P['nu']*P['alpha']) Nu = (0.825 + 0.387*Ra**(1/6)/(1+(0.492/P['Pr'])**(9/16))**(8/27))**2 # Churchill-Chu 1975 corr = 'Churchill-Chu (full range)' elif orient == 'up': Lc = L/4.0 # A/P for a square plate Ra = g*beta*dT*Lc**3/(P['nu']*P['alpha']) # Lloyd-Moran 1974 (Incropera eq. 9.30/9.31) if Ra < 1e7: NuLM = 0.54*Ra**0.25 else: NuLM = 0.15*Ra**(1/3) # Raithby-Hollands blended (valid through low Ra, thick-BL correction) Cl = 0.671/(1+(0.492/P['Pr'])**(9/16))**(4/9) NuT = 0.835*Cl*Ra**0.25 Nul = 1.4/math.log(1+1.4/NuT) Nut = 0.14*Ra**(1/3) NuRH = (Nul**10 + Nut**10)**0.1 Nu = NuRH; corr = 'Raithby-Hollands (LM %.2f)' % (NuLM*P['k']/Lc) elif orient == 'down': Lc = L/4.0 Ra = g*beta*dT*Lc**3/(P['nu']*P['alpha']) Nu = 0.52*Ra**0.2; corr = 'Incropera 9.32 (0.52 Ra^1/5)' h = Nu*P['k']/Lc regime = 'laminar' if (Ra < 1e9 if orient=='vert' else Ra < 1e7) else 'turbulent' return dict(Ra=Ra, Gr=Ra/P['Pr'], Nu=Nu, h=h, regime=regime, corr=corr, Lc=Lc) def forced(L, U, dT=60.0, Tinf=210.0, p=600.0, Retr_bypass=6e4): Tf = Tinf + dT/2; P = props(Tf, p) Re = U*L/P['nu']; Pr3 = P['Pr']**(1/3) Nu_lam = 0.664*Re**0.5*Pr3 # Blasius/Pohlhausen, Re<5e5 # bypass transition under high free-stream turbulence (Abu-Ghannam & Shaw: Re_x,tr ~6e4 at Tu 5-10 %) if Re > Retr_bypass: Nu_byp = (0.664*Retr_bypass**0.5 + 0.037*(Re**0.8 - Retr_bypass**0.8))*Pr3 else: Nu_byp = Nu_lam Nu_trip = max(Nu_lam, 0.037*Re**0.8*Pr3) # tripped at leading edge (upper bound) h_lam = Nu_lam*P['k']/L; h_byp = Nu_byp*P['k']/L; h_trip = Nu_trip*P['k']/L # mixed convection with natural (up-facing, n=3 assisting/transverse) nat = natural(L, dT, Tinf, p, 'up') h_mix = (h_lam**3 + nat['h']**3)**(1/3) Gr = g*(1/Tf)*dT*L**3/P['nu']**2 return dict(Re=Re, h_lam=h_lam, h_byp=h_byp, h_trip=h_trip, h_nat=nat['h'], h_mix=h_mix, GrRe2=Gr/Re**2) if __name__ == "__main__": print('== §1 properties') for T in (200,210,220,240,270): for p in (400,600,900): P = props(T,p) if p==600 or T==210: print(f"T={T} p={p}: rho={P['rho']:.5f} cp={P['cp']:.1f} mu={P['mu']*1e6:.3f}e-6 k={P['k']*1e3:.3f}e-3 Pr={P['Pr']:.3f} nu={P['nu']:.3e} alpha={P['alpha']:.3e}") # mean free path P = props(210); lam = P['mu']/600*math.sqrt(math.pi*R*210/2); print('mean free path 210K 600Pa:', lam) print('\n== §2 natural convection (Tinf=210 K, 600 Pa)') for orient in ('up','vert','down'): print('orient', orient) for L in (0.1,0.5,1,3,10): row=[] for dT in (20,60,120): r = natural(L,dT,orient=orient); row.append(f"dT{dT}: Ra={r['Ra']:.2e} {r['regime'][:3]} h={r['h']:.3f}") print(f" L={L}: " + ' | '.join(row) + (' ['+natural(L,60,orient=orient)['corr']+']' if orient=='up' else '')) print('pressure/temperature sensitivity, up, L=1, dT=60:') for p in (400,600,900): for Ti in (200,210,240): print(f" p={p} Tinf={Ti}: h={natural(1,60,Ti,p,'up')['h']:.3f}") print('\n== §3 forced convection (dT=60 K film, 600 Pa, 210 K)') for L in (0.001,0.0025,0.005,0.01,0.05,0.1,0.5,1,3,10): row=[] for U in (0.5,2,4,10,25): r = forced(L,U) row.append(f"U{U}: Re={r['Re']:.1e} lam={r['h_lam']:.2f} mix={r['h_mix']:.2f} byp={r['h_byp']:.2f} Gr/Re2={r['GrRe2']:.2g}") print(f"L={L}: " + ' | '.join(row)) print('\nforced sensitivity L=1 U=4 and L=0.5 U=10 over p,T,dT:') for (L,U) in ((1,4),(0.5,10),(1,10)): for p in (400,600,900): for dT in (20,60,120): r=forced(L,U,dT=dT,p=p) print(f" L={L} U={U} p={p} dT={dT}: lam={r['h_lam']:.2f} mix={r['h_mix']:.2f} byp={r['h_byp']:.2f}") # length scale for 3.5 at 0.5 m/s and 9 at 4 m/s (laminar forced+natural mixed, and pure laminar) def bis(f,a,b): fa=f(a) for _ in range(200): m=(a+b)/2; fm=f(m) if (fm>0)==(fa>0): a,fa=m,fm else: b=m return (a+b)/2 for key in ('h_lam','h_mix'): for (U,ht) in ((0.5,3.5),(4,9.0),(4,3.5),(10,3.5),(10,9.0),(25,3.5),(25,9.0),(0.5,1.5),(4,1.5),(10,1.5),(25,1.5)): f = lambda lnL: forced(math.exp(lnL),U)[key]-ht Lsol = math.exp(bis(f, math.log(1e-4), math.log(20))) print(f"{key}={ht} at U={U}: L={Lsol*1000:.2f} mm") # Bhandari 2013 check: 0.65 m, 15 m/s, 1067 Pa, 196 K fluid Pb = props(196+ (245-196)/2, 1066.6) Pf = props(196,1066.6) Re = 15*0.65/Pf['nu']; h = 0.664*Re**0.5*Pf['Pr']**(1/3)*Pf['k']/0.65 print(f"\nBhandari check (props at 196 K, 8 Torr): Re={Re:.3g} h={h:.2f} (paper: Re 2.9e4, h 1.4; CFD 1.8)") # Soria-Salinas check: rod 1.2 x 5 mm diag ~5.14 mm, 990 Pa, 298 K; Whitaker-type Nu=0.32+0.43Re^0.52 (as tabulated) Ps = props(298, 990); D = math.hypot(1.2e-3,5e-3) for U in (0.8,4,12): Re = U*D/Ps['nu']; Nu = 0.32+0.43*Re**0.52*Ps['Pr']**(1/3); print(f" HABIT rod U={U}: Re={Re:.1f} h={Nu*Ps['k']/D:.1f} (measured band 3-16 over 0.8-12 m/s)") # Bardera 2025 check: MMRTG ~0.64 m, natural, 600 Pa; vertical cylinder approx as vertical plate H~0.66 m for (Ts,Ta) in ((373,223),(423,173),(453,248),(423,273)): r = natural(0.66, Ts-Ta, Ta, 600, 'vert'); print(f" MMRTG-like vertical 0.66 m Ts={Ts} Ta={Ta}: h={r['h']:.3f}") # planetary bulk transfer (Martinez 2014 eq 5, neutral f=1) P = props(210) for z0 in (0.005,0.01,0.015): CH = 0.4**2/math.log(1.6/z0)**2 print(f" bulk: z0={z0}: C_H={CH:.4f}; h_eff/U = {P['rho']*736*CH:.4f}; h_eff at U=4: {P['rho']*736*CH*4:.3f}, U=10: {P['rho']*736*CH*10:.3f}")