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