🔬
Solver Purpose & Physical Scope
Estimate turbulent kinetic energy (k) and dissipation rate (epsilon) boundary conditions for solver setup.
📂 Discipline: Cfd
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 490 times
📄 Source File:
turbulent_ke.f90
📁 calcul/CFD /
turbulent_ke.f90
program turbulent_ke
implicit none
double precision, parameter :: C_mu = 0.09d0
double precision :: up, vp, wp, U_mean, delta, rho, mu, eps_in, L_in
double precision :: nu, k_val, rho_k, u_rms
double precision :: TI, TI_pct
double precision :: eps_val, L_val, lambda_val
double precision :: eta, tau_eta, u_eta
double precision :: nu_t_ke, nu_t_kw, omega_val, nut_ratio
double precision :: Re_t, Re_lam
double precision :: b11, b22, b33, II_inv, III_inv
double precision :: P_est
logical :: isotropic_flag
integer :: i, n_pts, iostat_val
double precision :: ti_cur, k_c, eps_c, nut_c, ret_c, eta_c, urms_c
! ---- read --------------------------------------------------------
read(*,*,iostat=iostat_val) up; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid u_prime.';stop;end if
read(*,*,iostat=iostat_val) vp; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid v_prime.';stop;end if
read(*,*,iostat=iostat_val) wp; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid w_prime.';stop;end if
read(*,*,iostat=iostat_val) U_mean; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid U_mean.';stop;end if
read(*,*,iostat=iostat_val) delta; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid delta.';stop;end if
read(*,*,iostat=iostat_val) rho; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid rho.';stop;end if
read(*,*,iostat=iostat_val) mu; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid mu.';stop;end if
read(*,*,iostat=iostat_val) eps_in; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid epsilon.';stop;end if
read(*,*,iostat=iostat_val) L_in; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid L_scale.';stop;end if
if(up<=0d0.and.vp<=0d0.and.wp<=0d0)then;write(*,*)'ERROR: At least one fluctuation must be > 0.';stop;end if
if(U_mean<=0d0)then;write(*,*)'ERROR: U_mean 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(delta<=0d0)then;write(*,*)'ERROR: delta must be > 0.';stop;end if
nu = mu / rho
! ---- TKE ---------------------------------------------------------
isotropic_flag = .false.
if (vp <= 0d0 .and. wp <= 0d0) then
isotropic_flag = .true.
vp = up; wp = up
end if
k_val = 0.5d0 * (up**2 + vp**2 + wp**2)
rho_k = rho * k_val
u_rms = sqrt(2.0d0 * k_val / 3.0d0)
TI = u_rms / U_mean
TI_pct = TI * 100.0d0
! ---- Length scale & dissipation ----------------------------------
if (L_in > 0d0) then
L_val = L_in
else
L_val = 0.09d0 * delta
end if
if (eps_in > 0d0) then
eps_val = eps_in
else
eps_val = C_mu**0.75d0 * k_val**1.5d0 / L_val
end if
! ---- Kolmogorov scales -------------------------------------------
eta = (nu**3 / eps_val)**0.25d0
tau_eta = sqrt(nu / eps_val)
u_eta = (nu * eps_val)**0.25d0
! ---- Taylor microscale -------------------------------------------
if (eps_val > 0d0) then
lambda_val = sqrt(10.0d0 * nu * k_val / eps_val)
else
lambda_val = 0d0
end if
! ---- Turbulent viscosity -----------------------------------------
if (eps_val > 0d0) then
nu_t_ke = C_mu * k_val**2 / eps_val
omega_val = eps_val / (C_mu * k_val)
nu_t_kw = k_val / omega_val
else
nu_t_ke = 0d0; omega_val = 0d0; nu_t_kw = 0d0
end if
nut_ratio = nu_t_ke / nu
! ---- Reynolds numbers --------------------------------------------
if (eps_val > 0d0 .and. nu > 0d0) then
Re_t = k_val**2 / (nu * eps_val)
else
Re_t = 0d0
end if
Re_lam = u_rms * lambda_val / nu
! ---- Anisotropy --------------------------------------------------
if (k_val > 0d0) then
b11 = up**2 / (2.0d0*k_val) - 1.0d0/3.0d0
b22 = vp**2 / (2.0d0*k_val) - 1.0d0/3.0d0
b33 = wp**2 / (2.0d0*k_val) - 1.0d0/3.0d0
else
b11=0d0; b22=0d0; b33=0d0
end if
II_inv = 2.0d0*(b11**2 + b22**2 + b33**2)
III_inv = 3.0d0*b11*b22*b33
! ---- Production estimate -----------------------------------------
P_est = nu_t_ke * (U_mean / delta)**2
! ---- output ------------------------------------------------------
write(*,'(A)') '============================================================'
write(*,'(A)') ' TURBULENT KINETIC ENERGY (k) CALCULATOR'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUT FLUCTUATIONS --------------------------------------'
write(*,'(A,F12.6,A)') ' u_prime (streamwise) = ', up, ' m/s'
write(*,'(A,F12.6,A)') ' v_prime (cross-stream) = ', vp, ' m/s'
write(*,'(A,F12.6,A)') ' w_prime (spanwise) = ', wp, ' m/s'
if (isotropic_flag) then
write(*,'(A)') ' ** Isotropic assumed (v_prime=w_prime=u_prime) **'
end if
write(*,'(A,F12.4,A)') ' U_mean = ', U_mean, ' m/s'
write(*,'(A,ES14.6,A)') ' Density (rho) = ', rho, ' kg/m3'
write(*,'(A,ES14.6,A)') ' Dynamic Viscosity (mu) = ', mu, ' Pa.s'
write(*,'(A,ES14.6,A)') ' Kinematic Visc (nu) = ', nu, ' m2/s'
write(*,'(A,ES14.6,A)') ' BL Thickness (delta) = ', delta, ' m'
write(*,*)
write(*,'(A)') '--- TURBULENT KINETIC ENERGY --------------------------------'
write(*,'(A,ES14.6,A)') ' k = 0.5*(u2+v2+w2) = ', k_val, ' m2/s2'
write(*,'(A,ES14.6,A)') ' rho*k (per volume) = ', rho_k, ' Pa'
write(*,'(A,F12.6,A)') ' u_rms = sqrt(2k/3) = ', u_rms, ' m/s'
write(*,*)
write(*,'(A)') '--- TURBULENCE INTENSITY ------------------------------------'
write(*,'(A,F12.6)') ' TI = u_rms / U_mean = ', TI
write(*,'(A,F12.4,A)') ' TI (percent) = ', TI_pct, ' %'
write(*,*)
write(*,'(A)') '--- DISSIPATION AND SCALES ----------------------------------'
write(*,'(A,ES14.6,A)') ' epsilon (dissipation) = ', eps_val, ' m2/s3'
write(*,'(A,ES14.6,A)') ' Integral Scale L = ', L_val, ' m'
write(*,'(A,ES14.6,A)') ' Taylor Microscale lam = ', lambda_val, ' m'
write(*,'(A,ES14.6,A)') ' Kolmogorov eta = ', eta, ' m'
write(*,'(A,ES14.6,A)') ' Kolmogorov tau_eta = ', tau_eta, ' s'
write(*,'(A,ES14.6,A)') ' Kolmogorov u_eta = ', u_eta, ' m/s'
write(*,*)
write(*,'(A)') '--- TURBULENT VISCOSITY -------------------------------------'
write(*,'(A,ES14.6,A)') ' nu_t (k-epsilon model) = ', nu_t_ke, ' m2/s'
write(*,'(A,F14.4)') ' nu_t / nu = ', nut_ratio
write(*,'(A,ES14.6,A)') ' omega (specific dissip) = ', omega_val, ' 1/s'
write(*,'(A,ES14.6,A)') ' nu_t (k-omega model) = ', nu_t_kw, ' m2/s'
write(*,*)
write(*,'(A)') '--- REYNOLDS NUMBERS ----------------------------------------'
write(*,'(A,ES14.6)') ' Re_t = k^2/(nu*eps) = ', Re_t
write(*,'(A,ES14.6)') ' Re_lambda = ', Re_lam
write(*,*)
write(*,'(A)') '--- ANISOTROPY (Lumley) -------------------------------------'
write(*,'(A,F12.6)') ' b11 (streamwise) = ', b11
write(*,'(A,F12.6)') ' b22 (cross-stream) = ', b22
write(*,'(A,F12.6)') ' b33 (spanwise) = ', b33
write(*,'(A,ES14.6)') ' II invariant = ', II_inv
write(*,'(A,ES14.6)') ' III invariant = ', III_inv
if (abs(b11) < 0.01d0 .and. abs(b22) < 0.01d0 .and. abs(b33) < 0.01d0) then
write(*,'(A)') ' Status: NEARLY ISOTROPIC'
else
write(*,'(A)') ' Status: ANISOTROPIC'
end if
write(*,*)
write(*,'(A)') '--- PRODUCTION ESTIMATE -------------------------------------'
write(*,'(A,ES14.6,A)') ' P ~ nu_t*(U/delta)^2 = ', P_est, ' m2/s3'
write(*,'(A,F12.4)') ' P / epsilon = ', P_est / eps_val
write(*,*)
! ---- profile sweep -----------------------------------------------
write(*,'(A)') '--- PROFILE vs TURBULENCE INTENSITY -------------------------'
write(*,'(A)') ' TI% k epsilon nu_t/nu Re_t eta'
write(*,'(A)') ' ------------------------------------------------------------------'
n_pts = 40
do i = 1, n_pts
ti_cur = dble(i) * 0.75d0 ! TI% from 0.75 to 30
urms_c = ti_cur / 100.0d0 * U_mean
k_c = 1.5d0 * urms_c**2
eps_c = C_mu**0.75d0 * k_c**1.5d0 / L_val
nut_c = C_mu * k_c**2 / eps_c / nu
ret_c = k_c**2 / (nu * eps_c)
eta_c = (nu**3 / eps_c)**0.25d0
write(*,'(F8.2,2X,ES12.4,2X,ES12.4,2X,F12.2,2X,ES12.4,2X,ES12.4)') &
ti_cur, k_c, eps_c, nut_c, ret_c, eta_c
end do
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
write(*,'(A)') ' k = 0.5*(u''^2 + v''^2 + w''^2)'
write(*,'(A)') ' TI = sqrt(2k/3) / U_mean'
write(*,'(A)') ' nu_t = C_mu * k^2 / epsilon (k-epsilon)'
write(*,'(A)') ' eta = (nu^3/epsilon)^(1/4) (Kolmogorov)'
write(*,'(A)') ' lambda = sqrt(10*nu*k/epsilon) (Taylor)'
write(*,'(A)') ' Re_t = k^2 / (nu*epsilon)'
write(*,'(A)') ' b_ij = <u_i u_j>/(2k) - d_ij/3 (anisotropy)'
write(*,'(A)') '============================================================'
end program turbulent_ke
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 turbulent_ke.f90 -o turbulent_ke
2. Execution with input.txt redirection:
turbulent_ke < input.txt
📄 Sample input.txt File Structure
Sample Data:
1.5 1.2 1.0 30 0.05 1.225 1.789e-5 0 0
Parameter Description:
up_init\nvp_init\nwp_init\nUmean_init\nδ (BL thickness) [m]\nÏ [kg/m³]\nμ [Pa·s]\nε [m²/s³]\nε [m²/s³]