Fortran Source Code & Offline Testing

Turbulent Prandtl Number Solver

Core Numerical Engine in Fortran 90 • 29 total downloads

turbulent_prandtl.f90

! =========================================================================
! Source File: turbulent_prandtl.f90
! =========================================================================

program turbulent_prandtl
    implicit none
    double precision :: nu_t, alpha_t_in, Pr_t_in
    double precision :: rho, mu, cp_val, k_cond, U_mean, dTdy, dudy
    double precision :: nu, alpha_mol, Pr_mol
    double precision :: Pr_t, alpha_t, nut_ratio
    double precision :: Pr_t_KC, Pr_t_JR, Pr_t_K94, Pr_t_RA, Pr_t_Ceb, Pr_t_const
    double precision :: q_lam, q_turb, q_total, k_eff, keff_ratio
    double precision :: ra_factor
    double precision :: x, arg
    integer :: i, n_pts, ios

    double precision :: nr, p1, p2, p3, p4, p5, p6, ker, ar

    read(*,*,iostat=ios) nu_t;       if(ios/=0)then;write(*,*)'ERROR: Invalid nu_t.';stop;end if
    read(*,*,iostat=ios) alpha_t_in;  if(ios/=0)then;write(*,*)'ERROR: Invalid alpha_t.';stop;end if
    read(*,*,iostat=ios) Pr_t_in;     if(ios/=0)then;write(*,*)'ERROR: Invalid Pr_t.';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) cp_val;      if(ios/=0)then;write(*,*)'ERROR: Invalid cp.';stop;end if
    read(*,*,iostat=ios) k_cond;      if(ios/=0)then;write(*,*)'ERROR: Invalid k_cond.';stop;end if
    read(*,*,iostat=ios) U_mean;      if(ios/=0)then;write(*,*)'ERROR: Invalid U_mean.';stop;end if
    read(*,*,iostat=ios) dTdy;        if(ios/=0)then;write(*,*)'ERROR: Invalid dTdy.';stop;end if
    read(*,*,iostat=ios) dudy;        if(ios/=0)then;write(*,*)'ERROR: Invalid dudy.';stop;end if

    if(nu_t<=0d0)then;write(*,*)'ERROR: nu_t 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(cp_val<=0d0)then;write(*,*)'ERROR: cp must be > 0.';stop;end if
    if(k_cond<=0d0)then;write(*,*)'ERROR: k_cond must be > 0.';stop;end if

    nu = mu / rho
    alpha_mol = k_cond / (rho * cp_val)
    Pr_mol = nu / alpha_mol

    nut_ratio = nu_t / nu

    ! Determine Pr_t and alpha_t
    if (Pr_t_in > 0d0) then
        Pr_t = Pr_t_in
        alpha_t = nu_t / Pr_t
    else if (alpha_t_in > 0d0) then
        alpha_t = alpha_t_in
        Pr_t = nu_t / alpha_t
    else
        Pr_t = 0.85d0
        alpha_t = nu_t / Pr_t
    end if

    ! 6 correlations for Pr_t
    ! 1. Kays & Crawford (simplified)
    x = nut_ratio
    if (x > 0.1d0) then
        arg = 5.165d0 / x
        if (arg > 500d0) arg = 500d0
        Pr_t_KC = 1d0 / (0.5882d0 + 0.228d0*x - 0.0441d0*x*x*(1d0-exp(-arg)))
        if (Pr_t_KC < 0.5d0) Pr_t_KC = 0.85d0
        if (Pr_t_KC > 3.0d0) Pr_t_KC = 0.85d0
    else
        Pr_t_KC = 0.85d0
    end if

    ! 2. Jischa & Rieke (1979): Pr_t = 0.85 + 0.015/Pr
    Pr_t_JR = 0.85d0 + 0.015d0 / Pr_mol

    ! 3. Kays (1994): Pr_t = 0.85 + 0.7/(Pr*(nu_t/nu))
    if (Pr_mol * nut_ratio > 1d-10) then
        Pr_t_K94 = 0.85d0 + 0.7d0 / (Pr_mol * nut_ratio)
    else
        Pr_t_K94 = 0.85d0
    end if

    ! 4. Reynolds analogy: Pr_t = 1.0
    Pr_t_RA = 1.0d0

    ! 5. Cebeci: near wall ~1.0, outer ~0.9 (average)
    Pr_t_Ceb = 0.9d0

    ! 6. Constant = 0.85 (standard CFD default)
    Pr_t_const = 0.85d0

    ! Heat transfer
    q_lam  = -k_cond * dTdy
    q_turb = -rho * cp_val * alpha_t * dTdy
    q_total = q_lam + q_turb
    k_eff = k_cond + rho * cp_val * alpha_t
    keff_ratio = k_eff / k_cond

    ! Reynolds analogy factor
    if (Pr_t > 0d0) then
        ra_factor = 1d0 / Pr_t
    else
        ra_factor = 1d0
    end if

    ! Output
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   TURBULENT PRANDTL NUMBER (Pr_t) CALCULATOR'
    write(*,'(A)') '============================================================'
    write(*,*)

    write(*,'(A)') '--- INPUT CONDITIONS ----------------------------------------'
    write(*,'(A,ES14.6,A)') '  nu_t (eddy viscosity)   = ', nu_t, ' m2/s'
    write(*,'(A,ES14.6,A)') '  Density (rho)           = ', rho, ' kg/m3'
    write(*,'(A,ES14.6,A)') '  Dyn. Viscosity (mu)     = ', mu, ' Pa.s'
    write(*,'(A,F12.4,A)')  '  Spec. Heat (cp)         = ', cp_val, ' J/kg.K'
    write(*,'(A,F12.6,A)')  '  Conductivity (k)        = ', k_cond, ' W/m.K'
    write(*,'(A,F12.4,A)')  '  U_mean                  = ', U_mean, ' m/s'
    write(*,'(A,F12.4,A)')  '  dT/dy                   = ', dTdy, ' K/m'
    write(*,'(A,F12.4,A)')  '  du/dy                   = ', dudy, ' 1/s'
    write(*,*)

    write(*,'(A)') '--- MOLECULAR PROPERTIES ------------------------------------'
    write(*,'(A,ES14.6,A)') '  Kinematic Visc (nu)     = ', nu, ' m2/s'
    write(*,'(A,ES14.6,A)') '  Thermal Diffus (alpha)  = ', alpha_mol, ' m2/s'
    write(*,'(A,F12.4)')    '  Prandtl Number (Pr)     = ', Pr_mol
    write(*,'(A,F14.4)')    '  nu_t / nu               = ', nut_ratio
    write(*,*)

    write(*,'(A)') '--- TURBULENT PRANDTL NUMBER --------------------------------'
    write(*,'(A,F12.6)')    '  Pr_t (used)             = ', Pr_t
    write(*,'(A,ES14.6,A)') '  alpha_t = nu_t/Pr_t     = ', alpha_t, ' m2/s'
    write(*,'(A,F12.4)')    '  alpha_t / alpha          = ', alpha_t / alpha_mol
    write(*,*)

    write(*,'(A)') '--- Pr_t CORRELATIONS (6 models) ----------------------------'
    write(*,'(A,F12.6)')    '  Kays-Crawford           = ', Pr_t_KC
    write(*,'(A,F12.6)')    '  Jischa-Rieke            = ', Pr_t_JR
    write(*,'(A,F12.6)')    '  Kays (1994)             = ', Pr_t_K94
    write(*,'(A,F12.6)')    '  Reynolds analogy         = ', Pr_t_RA
    write(*,'(A,F12.6)')    '  Cebeci                  = ', Pr_t_Ceb
    write(*,'(A,F12.6)')    '  Constant (default)      = ', Pr_t_const
    write(*,*)

    write(*,'(A)') '--- HEAT TRANSFER -------------------------------------------'
    write(*,'(A,ES14.6,A)') '  q_laminar               = ', q_lam, ' W/m2'
    write(*,'(A,ES14.6,A)') '  q_turbulent             = ', q_turb, ' W/m2'
    write(*,'(A,ES14.6,A)') '  q_total                 = ', q_total, ' W/m2'
    write(*,'(A,ES14.6,A)') '  k_eff (effective cond)  = ', k_eff, ' W/m.K'
    write(*,'(A,F14.4)')    '  k_eff / k               = ', keff_ratio
    write(*,*)

    write(*,'(A)') '--- REYNOLDS ANALOGY ----------------------------------------'
    write(*,'(A,F12.6)')    '  2*St/Cf = 1/Pr_t        = ', ra_factor
    if (abs(ra_factor - 1d0) < 0.05d0) then
        write(*,'(A)')      '  Status: FULL REYNOLDS ANALOGY (Pr_t ~ 1)'
    else
        write(*,'(A)')      '  Status: MODIFIED REYNOLDS ANALOGY'
    end if
    write(*,*)

    ! Profile sweep
    write(*,'(A)') '--- PROFILE vs nu_t/nu --------------------------------------'
    write(*,'(A)') '  nut/nu    Pr_t(KC)    Pr_t(JR)    Pr_t(K94)   keff/k      at/a'
    write(*,'(A)') '  ------------------------------------------------------------------'

    n_pts = 40
    do i = 1, n_pts
        nr = 1d0 + (999d0) * dble(i-1) / dble(n_pts-1)

        ! KC
        if (nr > 0.1d0) then
            arg = 5.165d0 / nr
            if (arg > 500d0) arg = 500d0
            p1 = 1d0 / (0.5882d0 + 0.228d0*nr - 0.0441d0*nr*nr*(1d0-exp(-arg)))
            if (p1 < 0.5d0 .or. p1 > 3d0) p1 = 0.85d0
        else
            p1 = 0.85d0
        end if

        ! JR
        p2 = 0.85d0 + 0.015d0 / Pr_mol

        ! K94
        if (Pr_mol * nr > 1d-10) then
            p3 = 0.85d0 + 0.7d0 / (Pr_mol * nr)
        else
            p3 = 0.85d0
        end if

        ! k_eff/k
        ker = 1d0 + Pr_mol * nr / p1

        ! alpha_t / alpha
        ar = nr / p1

        write(*,'(F10.2,2X,F10.4,2X,F10.4,2X,F10.4,2X,F10.2,2X,F10.2)') &
            nr, p1, p2, p3, ker, ar
    end do

    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
    write(*,'(A)') '  Pr_t = nu_t / alpha_t = epsilon_m / epsilon_H'
    write(*,'(A)') '  Jischa-Rieke: Pr_t = 0.85 + 0.015/Pr'
    write(*,'(A)') '  Kays (1994): Pr_t = 0.85 + 0.7/(Pr*nu_t/nu)'
    write(*,'(A)') '  k_eff = k * (1 + Pr*nu_t/(Pr_t*nu))'
    write(*,'(A)') '  q_turb = -rho*cp*alpha_t*dT/dy'
    write(*,'(A)') '  2*St/Cf = 1/Pr_t  (Reynolds analogy)'
    write(*,'(A)') '============================================================'

end program turbulent_prandtl

Solver Description

Calculate turbulent Prandtl number (Prt) profiles and eddy diffusivity ratios for heat transfer closure.

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_prandtl.f90 -o turbulent_prandtl

Execution Command:

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

turbulent_prandtl < input.txt

📥 Downloads & Local Files

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

! nut_i\nat_i\nprt_i\n[kg/m]\n[Pas]\ncp [J/kgK]\nk_i\nU_i\ndT_i\ndu_i
5e-3
! Parameter 2
0
! Parameter 3
0
! Parameter 4
1.225
! Parameter 5
1.789e-5
! Parameter 6
1006
! Parameter 7
0.0262
! Parameter 8
30
! Parameter 9
100
! Parameter 10
500

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

💻 Fortran Source Code Library

Turbulent Prandtl Number Solver

Solver Description

Key Numerical Methods & Architecture

🛠️ Local Compilation

📥 Downloads & Local Files