🔬
Solver Purpose & Physical Scope
Compute Kolmogorov length, time, and velocity scales for direct numerical simulation (DNS) mesh grids.
📂 Discipline: Cfd
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 384 times
📄 Source File:
kolmogorov_scales.f90
📁 calcul/CFD /
kolmogorov_scales.f90
program kolmogorov_scales
implicit none
double precision, parameter :: C_mu = 0.09d0, C_K = 1.5d0
double precision :: eps_in, nu_in, rho, mu, k_tke, U_mean, L_dom
double precision :: nu, eps, u_rms
double precision :: eta, tau_eta, u_eta
double precision :: lam_g, lam_iso, Re_lam
double precision :: L_int, T_int, u_int
double precision :: L_eta, L_lam, lam_eta
double precision :: Re_L, Re_dom
double precision :: N_grid, N_time, E_L, E_eta
integer :: i, n_pts, ios
double precision :: ec, etc, ttc, lmc, Lc, rlc, nc
read(*,*,iostat=ios) eps_in; if(ios/=0)then;write(*,*)'ERROR: Invalid epsilon.';stop;end if
read(*,*,iostat=ios) nu_in; if(ios/=0)then;write(*,*)'ERROR: Invalid nu.';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) k_tke; if(ios/=0)then;write(*,*)'ERROR: Invalid k.';stop;end if
read(*,*,iostat=ios) U_mean; if(ios/=0)then;write(*,*)'ERROR: Invalid U_mean.';stop;end if
read(*,*,iostat=ios) L_dom; if(ios/=0)then;write(*,*)'ERROR: Invalid L_domain.';stop;end if
if(eps_in<=0d0)then;write(*,*)'ERROR: epsilon must be > 0.';stop;end if
if(k_tke<=0d0)then;write(*,*)'ERROR: k must be > 0.';stop;end if
if (nu_in > 0d0) then
nu = nu_in
else
if(rho<=0d0.or.mu<=0d0)then;write(*,*)'ERROR: Need nu>0 or rho,mu>0.';stop;end if
nu = mu / rho
end if
eps = eps_in
! Kolmogorov scales
eta = (nu**3 / eps)**0.25d0
tau_eta = sqrt(nu / eps)
u_eta = (nu * eps)**0.25d0
! RMS velocity
u_rms = sqrt(2d0 * k_tke / 3d0)
! Taylor microscale
lam_g = sqrt(10d0 * nu * k_tke / eps)
lam_iso = sqrt(15d0 * nu * u_rms**2 / eps)
Re_lam = u_rms * lam_g / nu
! Integral scales
L_int = k_tke**1.5d0 / eps
T_int = k_tke / eps
u_int = u_rms
! Scale ratios
L_eta = L_int / eta
L_lam = L_int / lam_g
lam_eta = lam_g / eta
! Reynolds numbers
Re_L = u_rms * L_int / nu
if (L_dom > 0d0 .and. U_mean > 0d0) then
Re_dom = U_mean * L_dom / nu
else
Re_dom = 0d0
end if
! DNS estimation
N_grid = L_eta**3
N_time = T_int / tau_eta
! Energy
E_L = u_rms**2
E_eta = sqrt(nu * eps)
! Output
write(*,'(A)') '============================================================'
write(*,'(A)') ' KOLMOGOROV SCALES CALCULATOR'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUT CONDITIONS ----------------------------------------'
write(*,'(A,ES14.6,A)') ' epsilon (dissipation) = ', eps, ' m2/s3'
write(*,'(A,ES14.6,A)') ' Kinematic Visc (nu) = ', nu, ' m2/s'
write(*,'(A,ES14.6,A)') ' k (TKE) = ', k_tke, ' m2/s2'
write(*,'(A,F12.4,A)') ' u_rms = sqrt(2k/3) = ', u_rms, ' m/s'
write(*,'(A,F12.4,A)') ' U_mean = ', U_mean, ' m/s'
write(*,'(A,F12.4,A)') ' L_domain = ', L_dom, ' m'
write(*,*)
write(*,'(A)') '--- KOLMOGOROV SCALES (smallest) ----------------------------'
write(*,'(A,ES14.6,A)') ' eta (length) = ', eta, ' m'
write(*,'(A,ES14.6,A)') ' tau_eta (time) = ', tau_eta, ' s'
write(*,'(A,ES14.6,A)') ' u_eta (velocity) = ', u_eta, ' m/s'
write(*,'(A,F12.4)') ' Re_eta (= 1 by def.) = ', u_eta * eta / nu
write(*,*)
write(*,'(A)') '--- TAYLOR MICROSCALE ---------------------------------------'
write(*,'(A,ES14.6,A)') ' lambda (general) = ', lam_g, ' m'
write(*,'(A,ES14.6,A)') ' lambda (isotropic) = ', lam_iso, ' m'
write(*,'(A,F14.4)') ' Re_lambda = ', Re_lam
write(*,*)
write(*,'(A)') '--- INTEGRAL SCALES (largest) -------------------------------'
write(*,'(A,ES14.6,A)') ' L_int = k^3/2 / eps = ', L_int, ' m'
write(*,'(A,ES14.6,A)') ' T_int = k / eps = ', T_int, ' s'
write(*,'(A,F12.6,A)') ' u_int ~ u_rms = ', u_int, ' m/s'
write(*,*)
write(*,'(A)') '--- SCALE RATIOS --------------------------------------------'
write(*,'(A,ES14.6)') ' L / eta = ', L_eta
write(*,'(A,F14.4)') ' L / lambda = ', L_lam
write(*,'(A,F14.4)') ' lambda / eta = ', lam_eta
write(*,'(A,ES14.6)') ' T_int / tau_eta = ', T_int / tau_eta
write(*,'(A,F14.4)') ' u_int / u_eta = ', u_int / u_eta
write(*,*)
write(*,'(A)') '--- REYNOLDS NUMBERS ----------------------------------------'
write(*,'(A,ES14.6)') ' Re_L = u_rms*L/nu = ', Re_L
write(*,'(A,F14.4)') ' Re_lambda = ', Re_lam
write(*,'(A,ES14.6)') ' Re_domain = ', Re_dom
write(*,*)
write(*,'(A)') '--- DNS ESTIMATION ------------------------------------------'
write(*,'(A,ES14.6)') ' N_grid ~ (L/eta)^3 = ', N_grid
write(*,'(A,F14.2)') ' (L/eta)^(1/3) per dim = ', L_eta
write(*,'(A,ES14.6)') ' N_timesteps ~ T/tau_eta = ', T_int / tau_eta
write(*,'(A,ES14.6)') ' Total cost ~ N*Nt = ', N_grid * T_int / tau_eta
write(*,*)
write(*,'(A)') '--- ENERGY CASCADE ------------------------------------------'
write(*,'(A,ES14.6,A)') ' E_L ~ u_rms^2 = ', E_L, ' m2/s2'
write(*,'(A,ES14.6,A)') ' E_eta ~ sqrt(nu*eps) = ', E_eta, ' m2/s2'
write(*,'(A,F12.4)') ' Kolmogorov constant C_K = ', C_K
write(*,'(A)') ' E(kappa) = C_K * eps^(2/3) * kappa^(-5/3)'
write(*,*)
! Profile
write(*,'(A)') '--- PROFILE vs DISSIPATION RATE -----------------------------'
write(*,'(A)') ' epsilon eta tau_eta lambda L_int Re_lam (L/eta)'
write(*,'(A)') ' ---------------------------------------------------------------------------------'
n_pts = 40
do i = 1, n_pts
ec = eps / 100d0 * 10d0**(3d0 * dble(i-1) / dble(n_pts-1))
etc = (nu**3 / ec)**0.25d0
ttc = sqrt(nu / ec)
lmc = sqrt(10d0 * nu * k_tke / ec)
Lc = k_tke**1.5d0 / ec
rlc = u_rms * lmc / nu
nc = Lc / etc
write(*,'(ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,F10.2,2X,ES12.4)') &
ec, etc, ttc, lmc, Lc, rlc, nc
end do
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
write(*,'(A)') " eta = (nu^3/eps)^(1/4) (Kolmogorov length)"
write(*,'(A)') " tau_eta = (nu/eps)^(1/2) (Kolmogorov time)"
write(*,'(A)') " u_eta = (nu*eps)^(1/4) (Kolmogorov velocity)"
write(*,'(A)') " lambda = sqrt(10*nu*k/eps) (Taylor microscale)"
write(*,'(A)') " L = k^(3/2) / eps (integral scale)"
write(*,'(A)') " Re_lambda = u_rms*lambda/nu"
write(*,'(A)') " N_dns ~ (L/eta)^3 ~ Re_L^(9/4)"
write(*,'(A)') '============================================================'
end program kolmogorov_scales
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 kolmogorov_scales.f90 -o kolmogorov_scales
2. Execution with input.txt redirection:
kolmogorov_scales < input.txt
📄 Sample input.txt File Structure
Sample Data:
100 0 1.225 1.789e-5 2.5 30 0.5
Parameter Description:
eps_i\n[m/s]\n[kg/m]\n[Pas]\nk [m/s]\nU_i\nL_i