π» 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.
Prandtl Mixing Length Solver
Core Numerical Engine in Fortran 90 β’ 27 total downloads
prandtl_mixing_length.f90
! =========================================================================
! Source File: prandtl_mixing_length.f90
! =========================================================================
program prandtl_mixing_length
implicit none
double precision, parameter :: PI = 3.141592653589793d0
double precision :: y, delta, U_inf, rho, mu, kappa, A_plus, dp_dx
double precision :: nu, Re_delta, Cf, tau_w, u_tau, delta_nu
double precision :: y_plus, u_plus, u_local, du_dy
double precision :: l_pr, l_vd, l_out, l_comb
double precision :: nu_t, tau_turb, nut_ratio
character(len=30) :: layer
integer :: i, n_pts, iostat_val
double precision :: eta, yc, ypc, upc, ldc, nrc, bc, duy
! ---- read inputs --------------------------------------------------
read(*,*,iostat=iostat_val) y; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid y.';stop;end if
read(*,*,iostat=iostat_val) delta; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid delta.';stop;end if
read(*,*,iostat=iostat_val) U_inf; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid U_inf.';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) kappa; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid kappa.';stop;end if
read(*,*,iostat=iostat_val) A_plus; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid A_plus.';stop;end if
read(*,*,iostat=iostat_val) dp_dx; if(iostat_val/=0)then;write(*,*)'ERROR: Invalid dp_dx.';stop;end if
if(y<=0d0)then;write(*,*)'ERROR: y 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
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(kappa<=0d0)then;write(*,*)'ERROR: kappa must be > 0.';stop;end if
if(A_plus<=0d0)then;write(*,*)'ERROR: A_plus must be > 0.';stop;end if
nu = mu / rho
Re_delta = rho * U_inf * delta / mu
! Skin friction (power-law estimate for turbulent BL)
if (Re_delta > 1.0d0) then
Cf = 0.027d0 * Re_delta**(-1.0d0/7.0d0)
else
Cf = 0.01d0
end if
tau_w = 0.5d0 * Cf * rho * U_inf**2
u_tau = sqrt(tau_w / rho)
delta_nu = nu / u_tau
! Wall-units at y
y_plus = y * u_tau / nu
! Velocity from law of the wall
if (y_plus <= 5.0d0) then
u_plus = y_plus
layer = 'VISCOUS SUBLAYER'
else if (y_plus <= 30.0d0) then
u_plus = 5.0d0 * log(y_plus) - 3.05d0
layer = 'BUFFER LAYER'
else
u_plus = (1.0d0/kappa)*log(y_plus) + 5.0d0
layer = 'LOG-LAW REGION'
end if
if (y > delta) then
u_plus = U_inf / u_tau
layer = 'OUTER / WAKE'
end if
u_local = u_plus * u_tau
! Mixing lengths
l_pr = kappa * y ! Prandtl original
l_vd = kappa * y * (1.0d0 - exp(-y_plus / A_plus)) ! Van Driest
l_out = 0.09d0 * delta ! Clauser outer
l_comb = min(l_vd, l_out) ! combined
! du/dy from law of the wall
if (y_plus <= 5.0d0) then
du_dy = u_tau**2 / nu ! du+/dy+ = 1 => du/dy = u_tau / delta_nu
else if (y_plus <= 30.0d0) then
du_dy = 5.0d0 * u_tau / y ! du+/dy+ = 5/y+
else
du_dy = u_tau / (kappa * y) ! du+/dy+ = 1/(kappa*y+)
end if
! Turbulent quantities
nu_t = l_comb**2 * abs(du_dy)
tau_turb = rho * l_comb**2 * du_dy**2
nut_ratio = nu_t / nu
! ---- output -------------------------------------------------------
write(*,'(A)') '============================================================'
write(*,'(A)') ' PRANDTL MIXING LENGTH CALCULATOR'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- FLOW CONDITIONS -----------------------------------------'
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,F12.4,A)') ' Freestream Velocity = ', U_inf, ' m/s'
write(*,'(A,ES14.6,A)') ' BL Thickness (delta) = ', delta, ' m'
write(*,'(A,ES14.6)') ' Re_delta = ', Re_delta
write(*,'(A,F12.6)') ' Von Karman kappa = ', kappa
write(*,'(A,F12.4)') ' Van Driest A+ = ', A_plus
write(*,'(A,ES14.6,A)') ' Pressure Gradient = ', dp_dx, ' Pa/m'
write(*,*)
write(*,'(A)') '--- WALL PARAMETERS -----------------------------------------'
write(*,'(A,ES14.6)') ' Skin Friction Cf = ', Cf
write(*,'(A,ES14.6,A)') ' Wall Shear Stress tau_w = ', tau_w, ' Pa'
write(*,'(A,F12.6,A)') ' Friction Velocity u_tau = ', u_tau, ' m/s'
write(*,'(A,ES14.6,A)') ' Viscous Length delta_nu = ', delta_nu, ' m'
write(*,*)
write(*,'(A)') '--- WALL UNITS AT y -----------------------------------------'
write(*,'(A,ES14.6,A)') ' y (position) = ', y, ' m'
write(*,'(A,F12.4)') ' y/delta = ', y/delta
write(*,'(A,F14.4)') ' y+ = ', y_plus
write(*,'(A,F14.6)') ' u+ = ', u_plus
write(*,'(A,F12.4,A)') ' u (local velocity) = ', u_local, ' m/s'
write(*,'(A,A)') ' Layer Classification = ', trim(layer)
write(*,*)
write(*,'(A)') '--- MIXING LENGTH -------------------------------------------'
write(*,'(A,ES14.6,A)') ' l (Prandtl kappa*y) = ', l_pr, ' m'
write(*,'(A,ES14.6,A)') ' l (Van Driest damped) = ', l_vd, ' m'
write(*,'(A,ES14.6,A)') ' l (outer Clauser) = ', l_out, ' m'
write(*,'(A,ES14.6,A)') ' l (combined min) = ', l_comb, ' m'
write(*,'(A,F12.6)') ' l_combined / delta = ', l_comb / delta
write(*,'(A,F12.6)') ' l_combined / y = ', l_comb / y
write(*,*)
write(*,'(A)') '--- TURBULENT QUANTITIES ------------------------------------'
write(*,'(A,ES14.6,A)') ' du/dy (velocity grad) = ', du_dy, ' 1/s'
write(*,'(A,ES14.6,A)') ' nu_t (eddy viscosity) = ', nu_t, ' m2/s'
write(*,'(A,F14.4)') ' nu_t / nu = ', nut_ratio
write(*,'(A,ES14.6,A)') ' tau_turb (Reynolds str) = ', tau_turb, ' Pa'
write(*,'(A,F12.6)') ' tau_turb / tau_w = ', tau_turb / tau_w
write(*,*)
! ---- profile sweep ------------------------------------------------
write(*,'(A)') '--- PROFILE vs y/delta --------------------------------------'
write(*,'(A)') ' y/delta y+ u+ l/delta nu_t/nu'
write(*,'(A)') ' ----------------------------------------------------------'
n_pts = 40
do i = 1, n_pts
eta = dble(i) / dble(n_pts) ! y/delta = 0.025 .. 1.0
yc = eta * delta
ypc = yc * u_tau / nu
! velocity
if (ypc <= 5.0d0) then
upc = ypc
else if (ypc <= 30.0d0) then
upc = 5.0d0 * log(ypc) - 3.05d0
else
upc = (1.0d0/kappa)*log(ypc) + 5.0d0
end if
! mixing length
bc = kappa * yc * (1.0d0 - exp(-ypc / A_plus))
ldc = min(bc, 0.09d0*delta) / delta
! du/dy
if (ypc <= 5.0d0) then
duy = u_tau**2 / nu
else if (ypc <= 30.0d0) then
duy = 5.0d0 * u_tau / yc
else
duy = u_tau / (kappa * yc)
end if
nrc = (min(bc, 0.09d0*delta))**2 * abs(duy) / nu
write(*,'(F8.4,2X,F12.2,2X,F12.4,2X,F12.6,2X,F12.2)') &
eta, ypc, upc, ldc, nrc
end do
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
write(*,'(A)') ' Prandtl: l = kappa * y'
write(*,'(A)') ' Van Driest: l = kappa*y*[1-exp(-y+/A+)]'
write(*,'(A)') ' Outer: l = 0.09 * delta (Clauser)'
write(*,'(A)') ' Combined: l = min(l_VD, l_outer)'
write(*,'(A)') ' tau_turb = rho * l^2 * (du/dy)^2'
write(*,'(A)') ' nu_t = l^2 * |du/dy|'
write(*,'(A)') ' Law of wall: u+ = y+ (visc), u+ = (1/k)ln(y+)+5 (log)'
write(*,'(A)') '============================================================'
end program prandtl_mixing_length
Solver Description
Solve wall boundary layer eddy viscosity using Prandtl's mixing length hypothesis and van Driest damping.
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 prandtl_mixing_length.f90 -o prandtl_mixing_length
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
prandtl_mixing_length < input.txt
π₯ Downloads & Local Files
Preview of the required input file (input.txt):
! y (distance from wall) [m]\nΓΒ΄ (BL thickness) [m]\nUinf_init\nΓΒ [kg/mΓΒ³]\nΓΒΌ [PaΓΒ·s]\nΓΒΊ (von KΓΒ‘rmΓΒ‘n)\nAplus_init\ndpdx_init
0.005
! Parameter 2
0.05
! Parameter 3
30
! Parameter 4
1.225
! Parameter 5
1.789e-5
! Parameter 6
0.41
! Parameter 7
26
! Parameter 8
0
0.005
! Parameter 2
0.05
! Parameter 3
30
! Parameter 4
1.225
! Parameter 5
1.789e-5
! Parameter 6
0.41
! Parameter 7
26
! Parameter 8
0