🔬
Solver Purpose & Physical Scope
Evaluate the transonic similarity parameter (K) for thin airfoils at near-sonic flow speeds.
📂 Discipline: Cfd
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 509 times
📄 Source File:
transonic_similarity.f90
📁 calcul/CFD /
transonic_similarity.f90
program transonic_similarity
implicit none
double precision, parameter :: PI = 3.141592653589793d0
double precision :: M_inf, gamma, tau, Cp_inc, CL_inc
double precision :: chord, p_inf, T_inf, R_gas
double precision :: gm1, gp1
double precision :: rho_inf, a_inf, u_inf, q_inf
double precision :: K_param, tau_23, beta_val, inv_beta
double precision :: Cp_PG, Cp_bar, Cp_star, CL_PG
double precision :: M_dd_conv, M_dd_super
character(len=30) :: regime
logical :: is_sonic
integer :: i, n_points, iostat_val
double precision :: dM, M_cur, K_cur, b_cur, ib_cur, cp_cur
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) tau
if (iostat_val/=0) then; write(*,*) 'ERROR: Invalid tau.'; 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) 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 (tau<=0.0d0) then; write(*,*) 'ERROR: tau must be > 0.'; stop; end if
if (chord<=0.0d0) then; write(*,*) 'ERROR: chord must be > 0.'; 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
gm1 = gamma - 1.0d0; gp1 = gamma + 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
tau_23 = tau**(2.0d0/3.0d0)
K_param = (1.0d0 - M_inf**2)/tau_23
beta_val = sqrt(abs(1.0d0 - M_inf**2))
if (beta_val < 1.0d-12) beta_val = 1.0d-12
inv_beta = 1.0d0/beta_val
if (abs(1.0d0 - M_inf**2) > 1.0d-6) then
Cp_PG = Cp_inc/sqrt(abs(1.0d0 - M_inf**2))
else
Cp_PG = -9999.0d0
end if
Cp_bar = Cp_PG/tau_23
if (M_inf < 1.0d0 .and. beta_val > 1.0d-6) then
CL_PG = CL_inc/beta_val
else
CL_PG = CL_inc
end if
if (M_inf > 1.0d-6) then
Cp_star = (2.0d0/(gamma*M_inf**2))* &
(((2.0d0/gp1)*(1.0d0+0.5d0*gm1*M_inf**2))**(gamma/gm1) - 1.0d0)
else
Cp_star = -999.0d0
end if
is_sonic = .false.
if (M_inf > 0.01d0 .and. M_inf < 1.0d0) then
if (Cp_PG < Cp_star) is_sonic = .true.
end if
if (K_param > 2.0d0) then
regime = 'SUBSONIC-LIKE'
else if (K_param < -2.0d0) then
regime = 'SUPERSONIC-LIKE'
else
regime = 'TRANSONIC'
end if
M_dd_conv = 0.87d0 - CL_inc/10.0d0 - tau
M_dd_super = 0.95d0 - CL_inc/10.0d0 - tau
write(*,'(A)') '============================================================'
write(*,'(A)') ' TRANSONIC SIMILARITY PARAMETER 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)') ' Velocity (u_inf) = ', u_inf, ' m/s'
write(*,'(A,ES14.6,A)') ' Dynamic Pressure (q) = ', q_inf, ' Pa'
write(*,*)
write(*,'(A)') '--- AIRFOIL PARAMETERS --------------------------------------'
write(*,'(A,F12.6)') ' Thickness Ratio (tau) = ', tau
write(*,'(A,F12.8)') ' tau^(2/3) = ', tau_23
write(*,'(A,F12.6)') ' Cp (incompressible) = ', Cp_inc
write(*,'(A,F12.6)') ' CL (incompressible) = ', CL_inc
write(*,'(A,F12.4,A)') ' Chord Length = ', chord, ' m'
write(*,*)
write(*,'(A)') '--- TRANSONIC SIMILARITY PARAMETER --------------------------'
write(*,'(A,F14.8)') ' K = (1-M^2)/tau^(2/3) = ', K_param
write(*,'(A,A)') ' Flow Regime = ', trim(regime)
write(*,*)
write(*,'(A)') '--- COMPRESSIBILITY CORRECTIONS -----------------------------'
write(*,'(A,F12.8)') ' beta = sqrt(|1-M^2|) = ', beta_val
write(*,'(A,F12.6)') ' 1/beta (amplification) = ', inv_beta
write(*,'(A,F14.8)') ' Cp Prandtl-Glauert = ', Cp_PG
write(*,'(A,F14.8)') ' Cp_bar = Cp_PG/tau^2/3 = ', Cp_bar
write(*,'(A,F12.8)') ' CL Prandtl-Glauert = ', CL_PG
write(*,*)
write(*,'(A)') '--- CRITICAL CONDITION CHECK --------------------------------'
write(*,'(A,F14.8)') ' Cp* (sonic condition) = ', Cp_star
write(*,'(A,F14.8)') ' Cp_PG at this Mach = ', Cp_PG
if (is_sonic) then
write(*,'(A)') ' Status: LOCALLY SUPERSONIC (Cp_PG < Cp*)'
else
write(*,'(A)') ' Status: SUBSONIC EVERYWHERE'
end if
write(*,*)
write(*,'(A)') '--- DRAG DIVERGENCE MACH (KORN EQUATION) --------------------'
write(*,'(A,F12.6)') ' M_dd (conventional) = ', M_dd_conv
write(*,'(A,F12.6)') ' M_dd (supercritical) = ', M_dd_super
write(*,'(A)') ' Korn: M_dd + CL/10 + t/c = kappa'
write(*,'(A)') ' kappa = 0.87 (conventional), 0.95 (supercritical)'
if (M_inf > M_dd_conv) then
write(*,'(A)') ' *** DRAG DIVERGENCE: M_inf > M_dd(conv) ***'
end if
write(*,*)
write(*,'(A)') '--- FORCE ESTIMATES (per unit span) -------------------------'
write(*,'(A,ES14.6,A)') ' Lift / span = ', q_inf*chord*CL_PG, ' N/m'
write(*,*)
write(*,'(A)') '--- PROFILE vs MACH NUMBER ----------------------------------'
write(*,'(A)') ' M_inf K beta 1/beta Cp_PG'
write(*,'(A)') ' ----------------------------------------------------------'
n_points = 40
dM = 1.0d0/dble(n_points)
do i = 0, n_points
M_cur = 0.5d0 + dble(i)*dM
K_cur = (1.0d0 - M_cur**2)/tau_23
b_cur = sqrt(abs(1.0d0 - M_cur**2))
if (b_cur < 1.0d-8) b_cur = 1.0d-8
ib_cur = 1.0d0/b_cur
if (abs(1.0d0 - M_cur**2) > 1.0d-6) then
cp_cur = Cp_inc/sqrt(abs(1.0d0 - M_cur**2))
else
cp_cur = -9999.0d0
end if
write(*,'(F8.4,2X,F12.6,2X,F10.6,2X,F10.4,2X,F12.6)') &
M_cur, K_cur, b_cur, ib_cur, cp_cur
end do
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
write(*,'(A)') ' K = (1 - M^2) / tau^(2/3) (transonic similarity parameter)'
write(*,'(A)') ' Cp_PG = Cp_0 / sqrt(|1 - M^2|) (Prandtl-Glauert)'
write(*,'(A)') ' Cp_bar = Cp / tau^(2/3) (scaled pressure coefficient)'
write(*,'(A)') ' Korn: M_dd + CL/10 + t/c = kappa'
write(*,'(A)') ' TSD: (1-M^2)phi_xx + phi_yy = 0 (linear, breaks at M~1)'
write(*,'(A)') '============================================================'
end program transonic_similarity
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 transonic_similarity.f90 -o transonic_similarity
2. Execution with input.txt redirection:
transonic_similarity < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.8 1.4 0.12 -1.0 0.5 1.5 101325 300 287.058
Parameter Description:
Mach ($M_\infty$)\nγ\nThickness Ratio (τ = t/c)\ncpinc_init\nclinc_init\nChord [m]\np∞ [Pa]\nT∞ [K]\nGas