🔬
Solver Purpose & Physical Scope
For steady, 1D isentropic flow of an ideal gas with constant specific heats, local properties are related to stagnation properties via Mach number ($M$):
📂 Discipline: Fluid
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 502 times
📄 Source File:
compressible_flow.f90
📁 calcul/FluidMechanics /
compressible_flow.f90
program compressible_flow
implicit none
double precision :: M, gamma, T0, P0, rho0
double precision :: T_ratio, P_ratio, rho_ratio, A_ratio
double precision :: a, V, T_s, P_s, rho_s, a_star
double precision :: pi
integer :: mode ! 1=isentropic, 2=normal shock
character(len=30) :: flow_type
! Normal shock variables
double precision :: M2, P2_P1, T2_T1, rho2_rho1, P02_P01
double precision :: dS_R
integer :: i
double precision :: M_i
pi = 3.14159265358979d0
! Read inputs
read(*,*) mode ! 1=isentropic, 2=normal shock
read(*,*) M ! Mach number
read(*,*) gamma ! Specific heat ratio (1.4 for air)
read(*,*) T0 ! Stagnation temperature [K]
read(*,*) P0 ! Stagnation pressure [kPa]
! Stagnation density
rho0 = P0 * 1000.0d0 / (287.0d0 * T0) ! Assuming air (R=287 J/kg.K)
! ============================================
! ISENTROPIC FLOW RELATIONS
! ============================================
T_ratio = 1.0d0 + (gamma - 1.0d0) / 2.0d0 * M**2
T_s = T0 / T_ratio
P_ratio = T_ratio**(gamma / (gamma - 1.0d0))
P_s = P0 / P_ratio
rho_ratio = T_ratio**(1.0d0 / (gamma - 1.0d0))
rho_s = rho0 / rho_ratio
! Speed of sound
a = sqrt(gamma * 287.0d0 * T_s)
V = M * a
a_star = sqrt(2.0d0 * gamma * 287.0d0 * T0 / (gamma + 1.0d0))
! Area ratio A/A*
A_ratio = (1.0d0/M) * ((2.0d0/(gamma+1.0d0)) * (1.0d0 + (gamma-1.0d0)/2.0d0 * M**2)) &
**((gamma+1.0d0)/(2.0d0*(gamma-1.0d0)))
if (M < 1.0d0) then
flow_type = 'Subsonic'
else if (abs(M - 1.0d0) < 1.0d-6) then
flow_type = 'Sonic'
else
flow_type = 'Supersonic'
end if
! ============================================
! OUTPUT — ISENTROPIC
! ============================================
write(*,'(A)') '============================================================'
write(*,'(A)') ' COMPRESSIBLE FLOW CALCULATOR'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- ISENTROPIC FLOW RELATIONS ---'
write(*,'(A,F12.4)') ' Mach Number (M) = ', M
write(*,'(A,F12.4)') ' Specific Heat Ratio (gamma) = ', gamma
write(*,'(A,A)') ' Flow Type = ', trim(flow_type)
write(*,*)
write(*,'(A)') ' Stagnation Conditions:'
write(*,'(A,F12.2,A)') ' T₀ = ', T0, ' K'
write(*,'(A,F12.4,A)') ' P₀ = ', P0, ' kPa'
write(*,'(A,F12.4,A)') ' rho₀ = ', rho0, ' kg/m3'
write(*,*)
write(*,'(A)') ' Static Conditions:'
write(*,'(A,F12.2,A)') ' T = ', T_s, ' K'
write(*,'(A,F12.4,A)') ' P = ', P_s, ' kPa'
write(*,'(A,F12.4,A)') ' rho = ', rho_s, ' kg/m3'
write(*,*)
write(*,'(A)') ' Ratios:'
write(*,'(A,F12.6)') ' T/T₀ = ', 1.0d0/T_ratio
write(*,'(A,F12.6)') ' P/P₀ = ', 1.0d0/P_ratio
write(*,'(A,F12.6)') ' rho/rho₀ = ', 1.0d0/rho_ratio
write(*,'(A,F12.6)') ' A/A* = ', A_ratio
write(*,*)
write(*,'(A)') ' Velocity & Speed of Sound:'
write(*,'(A,F12.2,A)') ' Speed of Sound (a) = ', a, ' m/s'
write(*,'(A,F12.2,A)') ' Flow Velocity (V) = ', V, ' m/s'
write(*,'(A,F12.2,A)') ' Critical a* (sonic) = ', a_star, ' m/s'
write(*,*)
! ============================================
! NORMAL SHOCK (if supersonic)
! ============================================
if (mode == 2 .and. M > 1.0d0) then
! Normal shock relations
M2 = sqrt(((gamma-1.0d0)*M**2 + 2.0d0) / (2.0d0*gamma*M**2 - (gamma-1.0d0)))
P2_P1 = 1.0d0 + 2.0d0*gamma/(gamma+1.0d0) * (M**2 - 1.0d0)
T2_T1 = (1.0d0 + 2.0d0*gamma/(gamma+1.0d0)*(M**2-1.0d0)) * &
((2.0d0+(gamma-1.0d0)*M**2)/((gamma+1.0d0)*M**2))
rho2_rho1 = ((gamma+1.0d0)*M**2) / ((gamma-1.0d0)*M**2 + 2.0d0)
P02_P01 = (((gamma+1.0d0)*M**2)/((gamma-1.0d0)*M**2+2.0d0))**(gamma/(gamma-1.0d0)) &
* (1.0d0/(1.0d0+2.0d0*gamma/(gamma+1.0d0)*(M**2-1.0d0)))**(1.0d0/(gamma-1.0d0))
! Entropy change
dS_R = -log(P02_P01)
write(*,'(A)') '--- NORMAL SHOCK RELATIONS ---'
write(*,'(A,F12.4)') ' Upstream Mach (M1) = ', M
write(*,'(A,F12.6)') ' Downstream Mach (M2) = ', M2
write(*,*)
write(*,'(A)') ' Across Shock Ratios:'
write(*,'(A,F12.6)') ' P2/P1 = ', P2_P1
write(*,'(A,F12.6)') ' T2/T1 = ', T2_T1
write(*,'(A,F12.6)') ' rho2/rho1 = ', rho2_rho1
write(*,'(A,F12.6)') ' P₀2/P₀1 = ', P02_P01
write(*,*)
write(*,'(A)') ' Downstream Conditions:'
write(*,'(A,F12.2,A)') ' P2 = ', P_s * P2_P1, ' kPa'
write(*,'(A,F12.2,A)') ' T2 = ', T_s * T2_T1, ' K'
write(*,'(A,F12.4)') ' Entropy Change deltaS/R = ', dS_R
write(*,'(A,F12.2,A)') ' Stag. Pressure Loss = ', (1.0d0-P02_P01)*100, ' %'
write(*,*)
end if
! ============================================
! ISENTROPIC TABLE
! ============================================
write(*,'(A)') '--- ISENTROPIC FLOW TABLE ---'
write(*,'(A)') ' M T/T₀ P/P₀ rho/rho₀ A/A*'
write(*,'(A)') ' ---'
do i = 1, 30
M_i = 0.05d0 + dble(i-1) * 2.95d0 / 29.0d0
T_ratio = 1.0d0 + (gamma-1.0d0)/2.0d0 * M_i**2
write(*,'(F6.2,4X,F10.6,2X,F10.6,2X,F10.6,2X,F10.6)') M_i, &
1.0d0/T_ratio, &
1.0d0/T_ratio**(gamma/(gamma-1.0d0)), &
1.0d0/T_ratio**(1.0d0/(gamma-1.0d0)), &
(1.0d0/M_i) * ((2.0d0/(gamma+1.0d0)) * T_ratio) &
**((gamma+1.0d0)/(2.0d0*(gamma-1.0d0)))
end do
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ---'
write(*,'(A)') ' T₀/T = 1 + (gamma-1)/2 x m2'
write(*,'(A)') ' P₀/P = (T₀/T)^(gamma/(gamma-1))'
write(*,'(A)') ' rho₀/rho = (T₀/T)^(1/(gamma-1))'
write(*,'(A)') ' A/A* = (1/M)[(2/(gamma+1))(1+(gamma-1)/2 m2)]^((gamma+1)/(2(gamma-1)))'
end program compressible_flow
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 compressible_flow.f90 -o comp_flow
2. Execution with input.txt redirection:
comp_flow < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 1.5 1.4 300.0 101.325
Parameter Description:
Solver mode (1=Isentropic Relations, 2=Normal Shock Relations) Inlet Mach Number Specific Heat Ratio (gamma) Stagnation Temperature [K] Stagnation Pressure [kPa]