💻 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.
Critical Pressure Coefficient
Core Numerical Engine in Fortran 90 • 35 total downloads
critical_cp.f90
! =========================================================================
! Source File: critical_cp.f90
! =========================================================================
program critical_cp
implicit none
double precision, parameter :: PI = 3.141592653589793d0
double precision :: M_inf, gamma, Cp_min_incomp
double precision :: gp1, gm1
double precision :: Cp_star, Cp_vacuum, Cp_stag
double precision :: p_star_ratio, T_star_ratio, rho_star_ratio, a_star_ratio
double precision :: factor_inf, factor_sonic, p_ratio_sonic
double precision :: M_cr, Cp_star_cr, Cp_min_cr
logical :: has_Mcr
double precision :: M_lo, M_hi, M_mid, f_lo, f_mid, f_hi
integer :: iter, i, n_points
integer, parameter :: MAX_ITER = 500
double precision, parameter :: TOL = 1.0d-12
double precision :: dM, M_cur, cps_cur, cpv_cur, cpst_cur, fac_cur, fac_s
integer :: iostat_val
read(*,*,iostat=iostat_val) M_inf
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid Mach number.'; 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_min_incomp
if (iostat_val /= 0) then; write(*,*) 'ERROR: Invalid Cp_min.'; 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
gp1 = gamma + 1.0d0
gm1 = gamma - 1.0d0
factor_inf = 1.0d0 + 0.5d0 * gm1 * M_inf**2
factor_sonic = 2.0d0 / gp1 * factor_inf
p_ratio_sonic = factor_sonic ** (gamma / gm1)
Cp_star = (2.0d0 / (gamma * M_inf**2)) * (p_ratio_sonic - 1.0d0)
Cp_vacuum = -2.0d0 / (gamma * M_inf**2)
Cp_stag = (2.0d0 / (gamma * M_inf**2)) * (factor_inf ** (gamma / gm1) - 1.0d0)
T_star_ratio = factor_sonic
p_star_ratio = p_ratio_sonic
rho_star_ratio = p_star_ratio / T_star_ratio
a_star_ratio = sqrt(T_star_ratio)
has_Mcr = .false.
M_cr = 0.0d0
Cp_star_cr = 0.0d0
Cp_min_cr = 0.0d0
if (Cp_min_incomp < -1.0d-10) then
M_lo = 0.01d0
M_hi = 0.999d0
f_lo = cp_star_func(M_lo, gamma) - Cp_min_incomp / sqrt(1.0d0 - M_lo**2)
f_hi = cp_star_func(M_hi, gamma) - Cp_min_incomp / sqrt(1.0d0 - M_hi**2)
if (f_lo * f_hi < 0.0d0) then
do iter = 1, MAX_ITER
M_mid = 0.5d0 * (M_lo + M_hi)
f_mid = cp_star_func(M_mid, gamma) - Cp_min_incomp / sqrt(1.0d0 - M_mid**2)
if (abs(f_mid) < TOL .or. abs(M_hi - M_lo) < TOL) exit
f_lo = cp_star_func(M_lo, gamma) - Cp_min_incomp / sqrt(1.0d0 - M_lo**2)
if (f_lo * f_mid < 0.0d0) then
M_hi = M_mid
else
M_lo = M_mid
end if
end do
M_cr = M_mid
Cp_star_cr = cp_star_func(M_cr, gamma)
Cp_min_cr = Cp_min_incomp / sqrt(1.0d0 - M_cr**2)
has_Mcr = .true.
end if
end if
write(*,'(A)') '============================================================'
write(*,'(A)') ' CRITICAL PRESSURE COEFFICIENT CALCULATOR (Cp*)'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- FREESTREAM CONDITIONS -----------------------------------'
write(*,'(A,F12.4)') ' Freestream Mach (M_inf) = ', M_inf
write(*,'(A,F12.6)') ' Specific Heat Ratio (g) = ', gamma
write(*,'(A,F12.6)') ' Cp_min (incompressible) = ', Cp_min_incomp
write(*,*)
write(*,'(A)') '--- CRITICAL PRESSURE COEFFICIENT ---------------------------'
write(*,'(A,F14.8)') ' Cp* (sonic condition) = ', Cp_star
write(*,'(A,F14.8)') ' Cp_vacuum (p_local=0) = ', Cp_vacuum
write(*,'(A,F14.8)') ' Cp_stagnation (M_loc=0) = ', Cp_stag
write(*,*)
write(*,'(A)') '--- SONIC CONDITION RATIOS ----------------------------------'
write(*,'(A,F12.8)') ' p*/p_inf = ', p_star_ratio
write(*,'(A,F12.8)') ' T*/T_inf = ', T_star_ratio
write(*,'(A,F12.8)') ' rho*/rho_inf = ', rho_star_ratio
write(*,'(A,F12.8)') ' a*/a_inf = ', a_star_ratio
write(*,*)
if (has_Mcr) then
write(*,'(A)') '--- CRITICAL MACH NUMBER ------------------------------------'
write(*,'(A,F12.8)') ' Critical Mach M_cr = ', M_cr
write(*,'(A,F14.8)') ' Cp* at M_cr = ', Cp_star_cr
write(*,'(A,F14.8)') ' Cp_min at M_cr (P-G) = ', Cp_min_cr
write(*,*)
else
write(*,'(A)') '--- CRITICAL MACH NUMBER ------------------------------------'
write(*,'(A)') ' Critical Mach M_cr = N/A (no intersection found)'
write(*,*)
end if
write(*,'(A)') '--- Cp PROFILE vs FREESTREAM MACH ---------------------------'
write(*,'(A)') ' M_inf Cp* Cp_vacuum Cp_stag'
write(*,'(A)') ' ----------------------------------------------------------'
n_points = 40
if (M_inf > 1.0d0) then
dM = (min(M_inf * 1.5d0, 5.0d0) - 0.1d0) / dble(n_points)
else
dM = 0.95d0 / dble(n_points)
end if
do i = 0, n_points
M_cur = 0.1d0 + dble(i) * dM
if (M_cur < 0.01d0) M_cur = 0.01d0
fac_cur = 1.0d0 + 0.5d0 * gm1 * M_cur**2
fac_s = (2.0d0 / gp1 * fac_cur) ** (gamma / gm1)
cps_cur = (2.0d0 / (gamma * M_cur**2)) * (fac_s - 1.0d0)
cpv_cur = -2.0d0 / (gamma * M_cur**2)
cpst_cur = (2.0d0 / (gamma * M_cur**2)) * (fac_cur ** (gamma / gm1) - 1.0d0)
write(*,'(F8.4,2X,F12.8,2X,F12.8,2X,F12.8)') M_cur, cps_cur, cpv_cur, cpst_cur
end do
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
write(*,'(A)') ' Cp* = (2/(g*M^2)) * [ (2/(g+1) * (1+(g-1)/2*M^2))^(g/(g-1)) - 1 ]'
write(*,'(A)') ' Cp_vacuum = -2/(g*M^2)'
write(*,'(A)') ' Cp_stag = (2/(g*M^2)) * [ (1+(g-1)/2*M^2)^(g/(g-1)) - 1 ]'
write(*,'(A)') ' M_cr: solve Cp*(M) = Cp_min_inc / sqrt(1 - M^2) (Prandtl-Glauert)'
write(*,'(A)') '============================================================'
contains
function cp_star_func(M, gam) result(cps)
double precision, intent(in) :: M, gam
double precision :: cps, fac
fac = 1.0d0 + 0.5d0*(gam-1.0d0)*M**2
cps = (2.0d0/(gam*M**2)) * ((2.0d0/(gam+1.0d0)*fac)**(gam/(gam-1.0d0)) - 1.0d0)
end function cp_star_func
end program critical_cp
Solver Description
Calculate critical pressure coefficient (Cp*) for airfoils in compressible flow systems.
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 critical_cp.f90 -o critical_cp
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
critical_cp < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Min Cp (incompressible) ($C_{p,\min}$)\nSpecific Heat Ratio ($\gamma$)\ncpmin_init
0.7
! Parameter 2
1.4
! Parameter 3
-1.2
0.7
! Parameter 2
1.4
! Parameter 3
-1.2