🔬
Solver Purpose & Physical Scope
Calculate critical pressure coefficient (Cp*) for airfoils in compressible flow systems.
📂 Discipline: Cfd
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 305 times
📄 Source File:
critical_cp.f90
📁 calcul/CFD /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 critical_cp.f90 -o critical_cp
2. Execution with input.txt redirection:
critical_cp < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.7 1.4 -1.2
Parameter Description:
Min Cp (incompressible) ($C_{p,\min}$)\nSpecific Heat Ratio ($\gamma$)\ncpmin_init