🔬
Solver Purpose & Physical Scope
Calculate turbulent Prandtl number (Prt) profiles and eddy diffusivity ratios for heat transfer closure.
📂 Discipline: Cfd
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 344 times
📄 Source File:
turbulent_prandtl.f90
📁 calcul/CFD /
turbulent_prandtl.f90
program turbulent_prandtl
implicit none
double precision :: nu_t, alpha_t_in, Pr_t_in
double precision :: rho, mu, cp_val, k_cond, U_mean, dTdy, dudy
double precision :: nu, alpha_mol, Pr_mol
double precision :: Pr_t, alpha_t, nut_ratio
double precision :: Pr_t_KC, Pr_t_JR, Pr_t_K94, Pr_t_RA, Pr_t_Ceb, Pr_t_const
double precision :: q_lam, q_turb, q_total, k_eff, keff_ratio
double precision :: ra_factor
double precision :: x, arg
integer :: i, n_pts, ios
double precision :: nr, p1, p2, p3, p4, p5, p6, ker, ar
read(*,*,iostat=ios) nu_t; if(ios/=0)then;write(*,*)'ERROR: Invalid nu_t.';stop;end if
read(*,*,iostat=ios) alpha_t_in; if(ios/=0)then;write(*,*)'ERROR: Invalid alpha_t.';stop;end if
read(*,*,iostat=ios) Pr_t_in; if(ios/=0)then;write(*,*)'ERROR: Invalid Pr_t.';stop;end if
read(*,*,iostat=ios) rho; if(ios/=0)then;write(*,*)'ERROR: Invalid rho.';stop;end if
read(*,*,iostat=ios) mu; if(ios/=0)then;write(*,*)'ERROR: Invalid mu.';stop;end if
read(*,*,iostat=ios) cp_val; if(ios/=0)then;write(*,*)'ERROR: Invalid cp.';stop;end if
read(*,*,iostat=ios) k_cond; if(ios/=0)then;write(*,*)'ERROR: Invalid k_cond.';stop;end if
read(*,*,iostat=ios) U_mean; if(ios/=0)then;write(*,*)'ERROR: Invalid U_mean.';stop;end if
read(*,*,iostat=ios) dTdy; if(ios/=0)then;write(*,*)'ERROR: Invalid dTdy.';stop;end if
read(*,*,iostat=ios) dudy; if(ios/=0)then;write(*,*)'ERROR: Invalid dudy.';stop;end if
if(nu_t<=0d0)then;write(*,*)'ERROR: nu_t must be > 0.';stop;end if
if(rho<=0d0)then;write(*,*)'ERROR: rho must be > 0.';stop;end if
if(mu<=0d0)then;write(*,*)'ERROR: mu must be > 0.';stop;end if
if(cp_val<=0d0)then;write(*,*)'ERROR: cp must be > 0.';stop;end if
if(k_cond<=0d0)then;write(*,*)'ERROR: k_cond must be > 0.';stop;end if
nu = mu / rho
alpha_mol = k_cond / (rho * cp_val)
Pr_mol = nu / alpha_mol
nut_ratio = nu_t / nu
! Determine Pr_t and alpha_t
if (Pr_t_in > 0d0) then
Pr_t = Pr_t_in
alpha_t = nu_t / Pr_t
else if (alpha_t_in > 0d0) then
alpha_t = alpha_t_in
Pr_t = nu_t / alpha_t
else
Pr_t = 0.85d0
alpha_t = nu_t / Pr_t
end if
! 6 correlations for Pr_t
! 1. Kays & Crawford (simplified)
x = nut_ratio
if (x > 0.1d0) then
arg = 5.165d0 / x
if (arg > 500d0) arg = 500d0
Pr_t_KC = 1d0 / (0.5882d0 + 0.228d0*x - 0.0441d0*x*x*(1d0-exp(-arg)))
if (Pr_t_KC < 0.5d0) Pr_t_KC = 0.85d0
if (Pr_t_KC > 3.0d0) Pr_t_KC = 0.85d0
else
Pr_t_KC = 0.85d0
end if
! 2. Jischa & Rieke (1979): Pr_t = 0.85 + 0.015/Pr
Pr_t_JR = 0.85d0 + 0.015d0 / Pr_mol
! 3. Kays (1994): Pr_t = 0.85 + 0.7/(Pr*(nu_t/nu))
if (Pr_mol * nut_ratio > 1d-10) then
Pr_t_K94 = 0.85d0 + 0.7d0 / (Pr_mol * nut_ratio)
else
Pr_t_K94 = 0.85d0
end if
! 4. Reynolds analogy: Pr_t = 1.0
Pr_t_RA = 1.0d0
! 5. Cebeci: near wall ~1.0, outer ~0.9 (average)
Pr_t_Ceb = 0.9d0
! 6. Constant = 0.85 (standard CFD default)
Pr_t_const = 0.85d0
! Heat transfer
q_lam = -k_cond * dTdy
q_turb = -rho * cp_val * alpha_t * dTdy
q_total = q_lam + q_turb
k_eff = k_cond + rho * cp_val * alpha_t
keff_ratio = k_eff / k_cond
! Reynolds analogy factor
if (Pr_t > 0d0) then
ra_factor = 1d0 / Pr_t
else
ra_factor = 1d0
end if
! Output
write(*,'(A)') '============================================================'
write(*,'(A)') ' TURBULENT PRANDTL NUMBER (Pr_t) CALCULATOR'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUT CONDITIONS ----------------------------------------'
write(*,'(A,ES14.6,A)') ' nu_t (eddy viscosity) = ', nu_t, ' m2/s'
write(*,'(A,ES14.6,A)') ' Density (rho) = ', rho, ' kg/m3'
write(*,'(A,ES14.6,A)') ' Dyn. Viscosity (mu) = ', mu, ' Pa.s'
write(*,'(A,F12.4,A)') ' Spec. Heat (cp) = ', cp_val, ' J/kg.K'
write(*,'(A,F12.6,A)') ' Conductivity (k) = ', k_cond, ' W/m.K'
write(*,'(A,F12.4,A)') ' U_mean = ', U_mean, ' m/s'
write(*,'(A,F12.4,A)') ' dT/dy = ', dTdy, ' K/m'
write(*,'(A,F12.4,A)') ' du/dy = ', dudy, ' 1/s'
write(*,*)
write(*,'(A)') '--- MOLECULAR PROPERTIES ------------------------------------'
write(*,'(A,ES14.6,A)') ' Kinematic Visc (nu) = ', nu, ' m2/s'
write(*,'(A,ES14.6,A)') ' Thermal Diffus (alpha) = ', alpha_mol, ' m2/s'
write(*,'(A,F12.4)') ' Prandtl Number (Pr) = ', Pr_mol
write(*,'(A,F14.4)') ' nu_t / nu = ', nut_ratio
write(*,*)
write(*,'(A)') '--- TURBULENT PRANDTL NUMBER --------------------------------'
write(*,'(A,F12.6)') ' Pr_t (used) = ', Pr_t
write(*,'(A,ES14.6,A)') ' alpha_t = nu_t/Pr_t = ', alpha_t, ' m2/s'
write(*,'(A,F12.4)') ' alpha_t / alpha = ', alpha_t / alpha_mol
write(*,*)
write(*,'(A)') '--- Pr_t CORRELATIONS (6 models) ----------------------------'
write(*,'(A,F12.6)') ' Kays-Crawford = ', Pr_t_KC
write(*,'(A,F12.6)') ' Jischa-Rieke = ', Pr_t_JR
write(*,'(A,F12.6)') ' Kays (1994) = ', Pr_t_K94
write(*,'(A,F12.6)') ' Reynolds analogy = ', Pr_t_RA
write(*,'(A,F12.6)') ' Cebeci = ', Pr_t_Ceb
write(*,'(A,F12.6)') ' Constant (default) = ', Pr_t_const
write(*,*)
write(*,'(A)') '--- HEAT TRANSFER -------------------------------------------'
write(*,'(A,ES14.6,A)') ' q_laminar = ', q_lam, ' W/m2'
write(*,'(A,ES14.6,A)') ' q_turbulent = ', q_turb, ' W/m2'
write(*,'(A,ES14.6,A)') ' q_total = ', q_total, ' W/m2'
write(*,'(A,ES14.6,A)') ' k_eff (effective cond) = ', k_eff, ' W/m.K'
write(*,'(A,F14.4)') ' k_eff / k = ', keff_ratio
write(*,*)
write(*,'(A)') '--- REYNOLDS ANALOGY ----------------------------------------'
write(*,'(A,F12.6)') ' 2*St/Cf = 1/Pr_t = ', ra_factor
if (abs(ra_factor - 1d0) < 0.05d0) then
write(*,'(A)') ' Status: FULL REYNOLDS ANALOGY (Pr_t ~ 1)'
else
write(*,'(A)') ' Status: MODIFIED REYNOLDS ANALOGY'
end if
write(*,*)
! Profile sweep
write(*,'(A)') '--- PROFILE vs nu_t/nu --------------------------------------'
write(*,'(A)') ' nut/nu Pr_t(KC) Pr_t(JR) Pr_t(K94) keff/k at/a'
write(*,'(A)') ' ------------------------------------------------------------------'
n_pts = 40
do i = 1, n_pts
nr = 1d0 + (999d0) * dble(i-1) / dble(n_pts-1)
! KC
if (nr > 0.1d0) then
arg = 5.165d0 / nr
if (arg > 500d0) arg = 500d0
p1 = 1d0 / (0.5882d0 + 0.228d0*nr - 0.0441d0*nr*nr*(1d0-exp(-arg)))
if (p1 < 0.5d0 .or. p1 > 3d0) p1 = 0.85d0
else
p1 = 0.85d0
end if
! JR
p2 = 0.85d0 + 0.015d0 / Pr_mol
! K94
if (Pr_mol * nr > 1d-10) then
p3 = 0.85d0 + 0.7d0 / (Pr_mol * nr)
else
p3 = 0.85d0
end if
! k_eff/k
ker = 1d0 + Pr_mol * nr / p1
! alpha_t / alpha
ar = nr / p1
write(*,'(F10.2,2X,F10.4,2X,F10.4,2X,F10.4,2X,F10.2,2X,F10.2)') &
nr, p1, p2, p3, ker, ar
end do
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
write(*,'(A)') ' Pr_t = nu_t / alpha_t = epsilon_m / epsilon_H'
write(*,'(A)') ' Jischa-Rieke: Pr_t = 0.85 + 0.015/Pr'
write(*,'(A)') ' Kays (1994): Pr_t = 0.85 + 0.7/(Pr*nu_t/nu)'
write(*,'(A)') ' k_eff = k * (1 + Pr*nu_t/(Pr_t*nu))'
write(*,'(A)') ' q_turb = -rho*cp*alpha_t*dT/dy'
write(*,'(A)') ' 2*St/Cf = 1/Pr_t (Reynolds analogy)'
write(*,'(A)') '============================================================'
end program turbulent_prandtl
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 turbulent_prandtl.f90 -o turbulent_prandtl
2. Execution with input.txt redirection:
turbulent_prandtl < input.txt
📄 Sample input.txt File Structure
Sample Data:
5e-3 0 0 1.225 1.789e-5 1006 0.0262 30 100 500
Parameter Description:
nut_i\nat_i\nprt_i\n[kg/m]\n[Pas]\ncp [J/kgK]\nk_i\nU_i\ndT_i\ndu_i