π» Fortran Source Code Library
Run calculations locally on your own machine. View code structure, read technical explanations, and download compilation packages including sample input files.
1D Isentropic Compressible Flow
Core Numerical Engine in Fortran 90 β’ 0 total downloads
compressible_flow.f90
! =========================================================================
! Source File: 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
Solver Description
Computes 1D isentropic compressible flow properties. Evaluates static-to-stagnation ratios for temperature ($T/T_0$), pressure ($P/P_0$), density ($\rho/\rho_0$), and checks sonic limit conditions (Mach number $M=1.0$).
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 -ffree-line-length-none compressible_flow.f90 -o comp_flow
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
comp_flow < input.txt
π₯ Downloads & Local Files
Preview of the required input file (input.txt):
! Mach Number ($M$)
2.0
! Specific Heat Ratio ($\gamma$)
1.4
! Stagnation Temperature ($T_0$) [K]
300.0
! Stagnation Pressure ($P_0$) [kPa]
100.0
2.0
! Specific Heat Ratio ($\gamma$)
1.4
! Stagnation Temperature ($T_0$) [K]
300.0
! Stagnation Pressure ($P_0$) [kPa]
100.0