Fortran Source Code & Offline Testing

Eddy Viscosity & Reynolds Stress

Core Numerical Engine in Fortran 90 • 21 total downloads

eddy_viscosity.f90

! =========================================================================
! Source File: eddy_viscosity.f90
! =========================================================================

program eddy_viscosity
    implicit none
    double precision, parameter :: C_mu = 0.09d0, Cs_smag = 0.1d0
    double precision, parameter :: A_plus = 26.0d0, Pr_t = 0.85d0

    double precision :: k_in, eps_in, omega_in, l_in, dudy
    double precision :: rho, mu, delta, y_pos, U_inf
    double precision :: nu, Re_delta, Cf, tau_w, u_tau, y_plus, f_mu
    double precision :: k_val, eps_val, omega_val, l_val
    double precision :: nt_ke, nt_kw, nt_ml, nt_1eq, nt_bl, nt_smag
    double precision :: mu_eff, nut_ratio, mueff_ratio
    double precision :: tau_turb, tau_lam, tau_total
    double precision :: nt_damped, Re_t, alpha_t
    double precision :: Delta_les
    integer :: i, n_pts, ios
    double precision :: kc, ec, oc, ntc1, ntc2, nrc, rtc, mrc

    read(*,*,iostat=ios) k_in;    if(ios/=0)then;write(*,*)'ERROR: Invalid k.';stop;end if
    read(*,*,iostat=ios) eps_in;  if(ios/=0)then;write(*,*)'ERROR: Invalid epsilon.';stop;end if
    read(*,*,iostat=ios) omega_in;if(ios/=0)then;write(*,*)'ERROR: Invalid omega.';stop;end if
    read(*,*,iostat=ios) l_in;    if(ios/=0)then;write(*,*)'ERROR: Invalid l_mix.';stop;end if
    read(*,*,iostat=ios) dudy;    if(ios/=0)then;write(*,*)'ERROR: Invalid dudy.';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) delta;   if(ios/=0)then;write(*,*)'ERROR: Invalid delta.';stop;end if
    read(*,*,iostat=ios) y_pos;   if(ios/=0)then;write(*,*)'ERROR: Invalid y.';stop;end if
    read(*,*,iostat=ios) U_inf;   if(ios/=0)then;write(*,*)'ERROR: Invalid U_inf.';stop;end if

    if(k_in<=0d0)then;write(*,*)'ERROR: k 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
    if(U_inf<=0d0)then;write(*,*)'ERROR: U_inf must be > 0.';stop;end if

    nu = mu / rho
    k_val = k_in

    if (eps_in > 0d0) then
        eps_val = eps_in
    else
        if (l_in > 0d0) then
            eps_val = C_mu**0.75d0 * k_val**1.5d0 / l_in
        else
            eps_val = C_mu**0.75d0 * k_val**1.5d0 / (0.09d0 * delta)
        end if
    end if

    if (omega_in > 0d0) then
        omega_val = omega_in
    else
        omega_val = eps_val / (C_mu * k_val)
    end if

    if (l_in > 0d0) then
        l_val = l_in
    else
        l_val = C_mu**0.75d0 * k_val**1.5d0 / eps_val
    end if

    nt_ke   = C_mu * k_val**2 / eps_val
    nt_kw   = k_val / omega_val
    nt_ml   = l_val**2 * abs(dudy)
    nt_1eq  = 0.55d0 * sqrt(k_val) * l_val
    nt_bl   = l_val**2 * abs(dudy)
    Delta_les = delta / 20.0d0
    nt_smag = (Cs_smag * Delta_les)**2 * abs(dudy)

    mu_eff = mu + rho * nt_ke
    nut_ratio = nt_ke / nu
    mueff_ratio = mu_eff / mu

    Re_delta = rho * U_inf * delta / mu
    if (Re_delta > 1d0) then
        Cf = 0.027d0 * Re_delta**(-1d0/7d0)
    else
        Cf = 0.01d0
    end if
    tau_w = 0.5d0 * Cf * rho * U_inf**2
    u_tau = sqrt(tau_w / rho)
    if (y_pos > 0d0) then
        y_plus = y_pos * u_tau / nu
    else
        y_plus = 0d0
    end if
    f_mu = 1d0 - exp(-y_plus / A_plus)
    nt_damped = nt_ke * f_mu**2

    tau_lam   = mu * dudy
    tau_turb  = rho * nt_ke * dudy
    tau_total = tau_lam + tau_turb

    Re_t = k_val**2 / (nu * eps_val)
    alpha_t = nt_ke / Pr_t

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   EDDY VISCOSITY (TURBULENT VISCOSITY) CALCULATOR'
    write(*,'(A)') '============================================================'
    write(*,*)

    write(*,'(A)') '--- INPUT CONDITIONS ----------------------------------------'
    write(*,'(A,ES14.6,A)') '  k (TKE)                 = ', k_val, ' m2/s2'
    write(*,'(A,ES14.6,A)') '  epsilon (dissipation)   = ', eps_val, ' m2/s3'
    write(*,'(A,ES14.6,A)') '  omega (specific dissip)  = ', omega_val, ' 1/s'
    write(*,'(A,ES14.6,A)') '  l_mix (mixing length)   = ', l_val, ' m'
    write(*,'(A,ES14.6,A)') '  du/dy (velocity grad)   = ', dudy, ' 1/s'
    write(*,'(A,ES14.6,A)') '  Density (rho)           = ', rho, ' kg/m3'
    write(*,'(A,ES14.6,A)') '  Dyn. 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(*,'(A,ES14.6,A)') '  y (wall distance)       = ', y_pos, ' m'
    write(*,'(A,F12.4,A)')  '  U_inf                   = ', U_inf, ' m/s'
    write(*,*)

    write(*,'(A)') '--- EDDY VISCOSITY MODELS -----------------------------------'
    write(*,'(A,ES14.6,A)') '  nu_t (k-epsilon)        = ', nt_ke, ' m2/s'
    write(*,'(A,ES14.6,A)') '  nu_t (k-omega)          = ', nt_kw, ' m2/s'
    write(*,'(A,ES14.6,A)') '  nu_t (mixing length)    = ', nt_ml, ' m2/s'
    write(*,'(A,ES14.6,A)') '  nu_t (Prandtl 1-eq)     = ', nt_1eq, ' m2/s'
    write(*,'(A,ES14.6,A)') '  nu_t (Baldwin-Lomax)    = ', nt_bl, ' m2/s'
    write(*,'(A,ES14.6,A)') '  nu_t (Smagorinsky LES)  = ', nt_smag, ' m2/s'
    write(*,*)

    write(*,'(A)') '--- EFFECTIVE VISCOSITY -------------------------------------'
    write(*,'(A,ES14.6,A)') '  mu_eff = mu + rho*nt    = ', mu_eff, ' Pa.s'
    write(*,'(A,F14.4)')    '  nu_t / nu               = ', nut_ratio
    write(*,'(A,F14.4)')    '  mu_eff / mu             = ', mueff_ratio
    write(*,*)

    write(*,'(A)') '--- WALL DAMPING (Van Driest) -------------------------------'
    write(*,'(A,F14.4)')    '  y+                      = ', y_plus
    write(*,'(A,F12.6)')    '  f_mu = 1-exp(-y+/A+)   = ', f_mu
    write(*,'(A,ES14.6,A)') '  nu_t_damped = nt*fmu^2  = ', nt_damped, ' m2/s'
    write(*,'(A,F14.4)')    '  nu_t_damped / nu        = ', nt_damped / nu
    write(*,*)

    write(*,'(A)') '--- SHEAR STRESS -------------------------------------------'
    write(*,'(A,ES14.6,A)') '  tau_laminar             = ', tau_lam, ' Pa'
    write(*,'(A,ES14.6,A)') '  tau_turbulent           = ', tau_turb, ' Pa'
    write(*,'(A,ES14.6,A)') '  tau_total               = ', tau_total, ' Pa'
    if (abs(tau_total) > 1d-30) then
        write(*,'(A,F12.6)')    '  tau_turb / tau_total    = ', tau_turb / tau_total
    end if
    write(*,*)

    write(*,'(A)') '--- DERIVED QUANTITIES --------------------------------------'
    write(*,'(A,ES14.6)')   '  Re_t = k^2/(nu*eps)     = ', Re_t
    write(*,'(A,F12.4)')    '  Pr_t (assumed)          = ', Pr_t
    write(*,'(A,ES14.6,A)') '  alpha_t = nu_t/Pr_t     = ', alpha_t, ' m2/s'
    write(*,*)

    write(*,'(A)') '--- PROFILE vs k -------------------------------------------'
    write(*,'(A)') '  k           nt_ke        nt_kw        nt/nu       Re_t         mueff/mu'
    write(*,'(A)') '  --------------------------------------------------------------------------'

    n_pts = 40
    do i = 1, n_pts
        kc = 0.1d0 + (max(k_val*3d0, 50d0) - 0.1d0) * dble(i-1) / dble(n_pts-1)
        ec = C_mu**0.75d0 * kc**1.5d0 / l_val
        oc = ec / (C_mu * kc)
        ntc1 = C_mu * kc**2 / ec
        ntc2 = kc / oc
        nrc = ntc1 / nu
        rtc = kc**2 / (nu * ec)
        mrc = (mu + rho * ntc1) / mu

        write(*,'(ES12.4,2X,ES12.4,2X,ES12.4,2X,F10.2,2X,ES12.4,2X,F10.4)') &
            kc, ntc1, ntc2, nrc, rtc, mrc
    end do

    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
    write(*,'(A)') '  k-eps:   nu_t = C_mu * k^2 / epsilon  (C_mu=0.09)'
    write(*,'(A)') '  k-omega: nu_t = k / omega'
    write(*,'(A)') '  Mixing:  nu_t = l^2 * |du/dy|'
    write(*,'(A)') '  1-eq:    nu_t = 0.55 * sqrt(k) * l'
    write(*,'(A)') '  Smag:    nu_t = (Cs*Delta)^2 * |S|  (Cs=0.1)'
    write(*,'(A)') '  Damping: f_mu = 1 - exp(-y+/26)'
    write(*,'(A)') '============================================================'

end program eddy_viscosity

Solver Description

Calculate turbulent eddy viscosity (mu_t) and Reynolds stress tensors for RANS turbulence models.

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 eddy_viscosity.f90 -o eddy_viscosity

Execution Command:

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

eddy_viscosity < input.txt

📥 Downloads & Local Files

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

! k [m/s]\ne_i\no_i\nl_i\nd_i\n[kg/m]\n[Pas]\ndel_i\ny_i\nU [m/s]
2.5
! Parameter 2
100
! Parameter 3
0
! Parameter 4
0
! Parameter 5
500
! Parameter 6
1.225
! Parameter 7
1.789e-5
! Parameter 8
0.05
! Parameter 9
0.005
! Parameter 10
30

📦 Download Offline Bundle (ZIP) 6 dl 📄 Download input.txt 8 dl 📥 Download eddy_viscosity.f90 7 dl

💻 Fortran Source Code Library

Eddy Viscosity & Reynolds Stress

Solver Description

Key Numerical Methods & Architecture

🛠️ Local Compilation

📥 Downloads & Local Files