πŸ” Search & Filter All 178+ Engineering Calculators β†’

πŸ’» Fortran Source Code Library

We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.

Turbulent Kinetic Energy (TKE)

Core Numerical Engine in Fortran 90 β€’ 22 total downloads

turbulent_ke.f90 
! =========================================================================
! Source File: 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

Solver Description

Estimate turbulent kinetic energy (k) and dissipation rate (epsilon) boundary conditions for solver setup.

Key Numerical Methods & Architecture

  • Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
  • Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
  • Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.

πŸ› οΈ Local Compilation

To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).

Compilation Command:

gfortran -O3 turbulent_ke.f90 -o turbulent_ke

Execution Command:

Execute the program by feeding the sample input file into the program using stdin redirection:

turbulent_ke < input.txt

πŸ“₯ Downloads & Local Files

Preview of the required input file (input.txt):

! up_init\nvp_init\nwp_init\nUmean_init\nδ (BL thickness) [m]\nρ [kg/m³]\nμ [Pa·s]\nΡ [m²/s³]\nΡ [m²/s³]
1.5
! Parameter 2
1.2
! Parameter 3
1.0
! Parameter 4
30
! Parameter 5
0.05
! Parameter 6
1.225
! Parameter 7
1.789e-5
! Parameter 8
0
! Parameter 9
0
πŸ“¦ Download Offline Bundle (ZIP) 8 dl πŸ“„ Download input.txt 6 dl πŸ“₯ Download turbulent_ke.f90 8 dl