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 = /(2k) - d_ij/3 (anisotropy)' write(*,'(A)') '============================================================' end program turbulent_ke