π» 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-Glauert Compressibility Correction
Core Numerical Engine in Fortran 90 β’ 19 total downloads
prandtl_glauert.f90
! =========================================================================
! Source File: prandtl_glauert.f90
! =========================================================================
program prandtl_glauert
implicit none
double precision, parameter :: PI = 3.141592653589793d0
double precision, parameter :: DEG2RAD = PI / 180.0d0
double precision :: M_inf, gamma, Cp_inc, CL_inc, alpha_deg
double precision :: chord, p_inf, T_inf, R_gas
double precision :: gm1, gp1, alpha_rad
double precision :: rho_inf, a_inf, u_inf, q_inf
double precision :: beta, beta_sq, inv_beta
double precision :: Cp_PG, Cp_KT, Cp_Laitone, denom_kt, denom_la
double precision :: CL_PG, dCLda_PG
double precision :: beta_sup, CL_Ack, Cd_wave, Cp_upper, Cp_lower
double precision :: L_span, D_wave_span, Cp_star
logical :: is_sub
integer :: i, n_points, iostat_val
double precision :: dM, M_cur, b_cur, cp_pg_val, cp_kt_val, cp_la_val, cl_cur, dk, dl
read(*,*,iostat=iostat_val) M_inf
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid Mach.'; stop; end if
read(*,*,iostat=iostat_val) gamma
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid gamma.'; stop; end if
read(*,*,iostat=iostat_val) Cp_inc
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid Cp_inc.'; stop; end if
read(*,*,iostat=iostat_val) CL_inc
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid CL_inc.'; stop; end if
read(*,*,iostat=iostat_val) alpha_deg
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid alpha.'; stop; end if
read(*,*,iostat=iostat_val) chord
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid chord.'; stop; end if
read(*,*,iostat=iostat_val) p_inf
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid pressure.'; stop; end if
read(*,*,iostat=iostat_val) T_inf
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid temperature.'; stop; end if
read(*,*,iostat=iostat_val) R_gas
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid R_gas.'; stop; end if
if (M_inf <= 0.0d0) then; write(*,*) 'ERROR: Mach must be > 0.'; stop; end if
if (gamma <= 1.0d0) then; write(*,*) 'ERROR: gamma must be > 1.'; stop; end if
if (p_inf <= 0.0d0) then; write(*,*) 'ERROR: Pressure must be > 0.'; stop; end if
if (T_inf <= 0.0d0) then; write(*,*) 'ERROR: Temperature must be > 0.'; stop; end if
if (R_gas <= 0.0d0) then; write(*,*) 'ERROR: R_gas must be > 0.'; stop; end if
if (chord <= 0.0d0) then; write(*,*) 'ERROR: Chord must be > 0.'; stop; end if
gm1 = gamma - 1.0d0
gp1 = gamma + 1.0d0
alpha_rad = alpha_deg * DEG2RAD
is_sub = (M_inf < 1.0d0)
rho_inf = p_inf / (R_gas * T_inf)
a_inf = sqrt(gamma * R_gas * T_inf)
u_inf = M_inf * a_inf
q_inf = 0.5d0 * rho_inf * u_inf**2
if (is_sub) then
beta_sq = 1.0d0 - M_inf**2
beta = sqrt(beta_sq)
inv_beta = 1.0d0 / beta
Cp_PG = Cp_inc / beta
CL_PG = CL_inc / beta
dCLda_PG = 2.0d0 * PI / beta
denom_kt = sqrt(beta_sq + M_inf**4 * ((Cp_inc + 1.0d0)**2) / 4.0d0)
if (abs(denom_kt) > 1.0d-15) then
Cp_KT = Cp_inc / denom_kt
else
Cp_KT = Cp_PG
end if
denom_la = beta + (M_inf**2 * (Cp_inc / 2.0d0)) / (1.0d0 + beta)
if (abs(denom_la) > 1.0d-15) then
Cp_Laitone = Cp_inc / denom_la
else
Cp_Laitone = Cp_PG
end if
L_span = q_inf * chord * CL_PG
D_wave_span = 0.0d0
Cp_star = (2.0d0 / (gamma * M_inf**2)) * &
( ((2.0d0/gp1)*(1.0d0+0.5d0*gm1*M_inf**2))**(gamma/gm1) - 1.0d0 )
Cp_upper = 0.0d0; Cp_lower = 0.0d0; Cd_wave = 0.0d0; CL_Ack = 0.0d0; beta_sup = 0.0d0
else
beta_sup = sqrt(M_inf**2 - 1.0d0)
beta = beta_sup; inv_beta = 1.0d0/beta
Cp_upper = -2.0d0 * alpha_rad / beta_sup
Cp_lower = 2.0d0 * alpha_rad / beta_sup
Cp_PG = Cp_upper; Cp_KT = Cp_upper; Cp_Laitone = Cp_upper
CL_Ack = 4.0d0 * alpha_rad / beta_sup
Cd_wave = 4.0d0 * alpha_rad**2 / beta_sup
CL_PG = CL_Ack
dCLda_PG = 4.0d0 / beta_sup
L_span = q_inf * chord * CL_Ack
D_wave_span = q_inf * chord * Cd_wave
Cp_star = 0.0d0
end if
write(*,'(A)') '============================================================'
write(*,'(A)') ' PRANDTL-GLAUERT COMPRESSIBILITY CORRECTION CALCULATOR'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- FREESTREAM CONDITIONS -----------------------------------'
write(*,'(A,F12.6)') ' Freestream Mach (M_inf) = ', M_inf
write(*,'(A,F12.6)') ' Specific Heat Ratio (g) = ', gamma
write(*,'(A,ES14.6,A)') ' Static Pressure (p_inf) = ', p_inf, ' Pa'
write(*,'(A,F12.2,A)') ' Static Temperature = ', T_inf, ' K'
write(*,'(A,F12.4,A)') ' Density (rho_inf) = ', rho_inf, ' kg/m3'
write(*,'(A,F12.4,A)') ' Speed of Sound (a_inf) = ', a_inf, ' m/s'
write(*,'(A,F12.4,A)') ' Freestream Velocity = ', u_inf, ' m/s'
write(*,'(A,ES14.6,A)') ' Dynamic Pressure (q) = ', q_inf, ' Pa'
write(*,*)
write(*,'(A)') '--- INPUT PARAMETERS ----------------------------------------'
write(*,'(A,F12.6)') ' Cp (incompressible) = ', Cp_inc
write(*,'(A,F12.6)') ' CL (incompressible) = ', CL_inc
write(*,'(A,F12.4,A)') ' Angle of Attack (alpha) = ', alpha_deg, ' deg'
write(*,'(A,F12.4,A)') ' Chord Length = ', chord, ' m'
write(*,*)
if (is_sub) then
write(*,'(A)') '--- COMPRESSIBILITY FACTOR (SUBSONIC) -----------------------'
write(*,'(A,F12.8)') ' beta = sqrt(1-M^2) = ', beta
write(*,'(A,F12.8)') ' 1/beta (amplification) = ', inv_beta
write(*,*)
write(*,'(A)') '--- PRESSURE COEFFICIENT CORRECTIONS ------------------------'
write(*,'(A,F14.8)') ' Cp Prandtl-Glauert = ', Cp_PG
write(*,'(A,F14.8)') ' Cp Karman-Tsien = ', Cp_KT
write(*,'(A,F14.8)') ' Cp Laitone = ', Cp_Laitone
write(*,'(A,F14.8)') ' Cp* (sonic reference) = ', Cp_star
write(*,*)
write(*,'(A)') '--- LIFT COEFFICIENT ----------------------------------------'
write(*,'(A,F12.8)') ' CL Prandtl-Glauert = ', CL_PG
write(*,'(A,F12.8)') ' dCL/dalpha (per rad) = ', dCLda_PG
write(*,'(A,F12.8)') ' dCL/dalpha (per deg) = ', dCLda_PG * DEG2RAD
write(*,*)
else
write(*,'(A)') '--- COMPRESSIBILITY FACTOR (SUPERSONIC) ---------------------'
write(*,'(A,F12.8)') ' beta = sqrt(M^2-1) = ', beta_sup
write(*,*)
write(*,'(A)') '--- ACKERET LINEARIZED THEORY -------------------------------'
write(*,'(A,F14.8)') ' Cp upper surface = ', Cp_upper
write(*,'(A,F14.8)') ' Cp lower surface = ', Cp_lower
write(*,'(A,F12.8)') ' CL (Ackeret) = ', CL_Ack
write(*,'(A,F12.8)') ' Cd_wave (Ackeret) = ', Cd_wave
write(*,'(A,F12.8)') ' dCL/dalpha (per rad) = ', dCLda_PG
write(*,'(A,F12.8)') ' dCL/dalpha (per deg) = ', dCLda_PG * DEG2RAD
write(*,*)
end if
write(*,'(A)') '--- FORCE ESTIMATES (per unit span) -------------------------'
write(*,'(A,ES14.6,A)') ' Lift / span = ', L_span, ' N/m'
if (.not. is_sub) then
write(*,'(A,ES14.6,A)') ' Wave Drag / span = ', D_wave_span, ' N/m'
end if
write(*,*)
write(*,'(A)') '--- CORRECTION PROFILE vs MACH NUMBER -----------------------'
if (is_sub) then
write(*,'(A)') ' M_inf beta Cp_PG Cp_KT Cp_Laitone CL_PG'
else
write(*,'(A)') ' M_inf beta_s Cp_upper Cp_lower CL_Ack Cd_wave'
end if
write(*,'(A)') ' ------------------------------------------------------------------'
n_points = 40
if (is_sub) then
dM = 0.90d0 / dble(n_points)
do i = 0, n_points
M_cur = 0.05d0 + dble(i) * dM
b_cur = sqrt(1.0d0 - M_cur**2)
cp_pg_val = Cp_inc / b_cur
dk = sqrt(1.0d0 - M_cur**2 + M_cur**4*((Cp_inc+1.0d0)**2)/4.0d0)
if (abs(dk)>1.0d-15) then; cp_kt_val = Cp_inc/dk; else; cp_kt_val = cp_pg_val; end if
dl = b_cur + (M_cur**2*(Cp_inc/2.0d0))/(1.0d0+b_cur)
if (abs(dl)>1.0d-15) then; cp_la_val = Cp_inc/dl; else; cp_la_val = cp_pg_val; end if
cl_cur = CL_inc / b_cur
write(*,'(F8.4,2X,F10.6,2X,F11.6,2X,F11.6,2X,F11.6,2X,F11.6)') &
M_cur, b_cur, cp_pg_val, cp_kt_val, cp_la_val, cl_cur
end do
else
dM = (min(M_inf*1.5d0,5.0d0) - 1.05d0) / dble(n_points)
do i = 0, n_points
M_cur = 1.05d0 + dble(i) * dM
b_cur = sqrt(M_cur**2 - 1.0d0)
write(*,'(F8.4,2X,F10.6,2X,F11.6,2X,F11.6,2X,F11.6,2X,F11.8)') &
M_cur, b_cur, &
-2.0d0*alpha_rad/b_cur, 2.0d0*alpha_rad/b_cur, &
4.0d0*alpha_rad/b_cur, 4.0d0*alpha_rad**2/b_cur
end do
end if
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
write(*,'(A)') ' Prandtl-Glauert: Cp(M) = Cp_0 / sqrt(1-M^2)'
write(*,'(A)') ' Karman-Tsien: Cp(M) = Cp_0 / sqrt(1-M^2+M^4*(Cp_0+1)^2/4)'
write(*,'(A)') ' Laitone: Cp(M) = Cp_0 / [beta + M^2*(Cp_0/2)/(1+beta)]'
write(*,'(A)') ' Ackeret (super): Cp = +/-2*alpha/sqrt(M^2-1)'
write(*,'(A)') ' CL_sub = CL_0/beta, CL_sup = 4*alpha/sqrt(M^2-1)'
write(*,'(A)') '============================================================'
end program prandtl_glauert
Solver Description
Apply Prandtl-Glauert rule to correct subsonic pressure coefficients for compressibility effects.
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_glauert.f90 -o prandtl_glauert
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
prandtl_glauert < input.txt
π₯ Downloads & Local Files
Preview of the required input file (input.txt):
! Mach ($M_\infty$)\nΓΒ³\ncpinc_init\nclinc_init\nΓΒ± [deg]\nChord [m]\npΓ’ΛΕΎ [Pa]\nTΓ’ΛΕΎ [K]\nGas
0.6
! Parameter 2
1.4
! Parameter 3
-1.0
! Parameter 4
0.5
! Parameter 5
5
! Parameter 6
1.5
! Parameter 7
101325
! Parameter 8
300
! Parameter 9
287.058
0.6
! Parameter 2
1.4
! Parameter 3
-1.0
! Parameter 4
0.5
! Parameter 5
5
! Parameter 6
1.5
! Parameter 7
101325
! Parameter 8
300
! Parameter 9
287.058