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