💻 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.
Rankine-Hugoniot Shock Relations
Core Numerical Engine in Fortran 90 • 30 total downloads
rankine_hugoniot.f90
! =========================================================================
! Source File: rankine_hugoniot.f90
! =========================================================================
program rankine_hugoniot
implicit none
! ---------------------------------------------------------------
! Normal Shock Relations (Rankine-Hugoniot)
! Reads: M1, gamma, p1 [Pa], T1 [K], R_gas [J/kg.K]
! ---------------------------------------------------------------
! Inputs
double precision :: M1, gamma, p1, T1, R_gas
! Derived upstream quantities
double precision :: rho1, a1, u1
! Shock ratios
double precision :: p_ratio, rho_ratio, T_ratio
double precision :: M2_sq, M2
double precision :: p0_ratio, ds_over_R, u_ratio
! Absolute downstream quantities
double precision :: p2, T2, rho2, a2, u2
! Stagnation quantities
double precision :: p01, p02, T01, T02
! Profile sweep variables
integer :: i, n_points
double precision :: dM, M_cur
double precision :: pr_cur, rr_cur, Tr_cur, M2s_cur, M2_cur, p0r_cur, ds_cur
! Helpers
double precision :: gp1, gm1, gp1h, gm1h ! gamma+1, gamma-1, halves
integer :: iostat_val
! ------------------------------------------------------------------
! Read inputs from stdin
! ------------------------------------------------------------------
read(*,*,iostat=iostat_val) M1
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid upstream Mach number input.'
stop
end if
read(*,*,iostat=iostat_val) gamma
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid specific heat ratio input.'
stop
end if
read(*,*,iostat=iostat_val) p1
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid upstream pressure input.'
stop
end if
read(*,*,iostat=iostat_val) T1
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid upstream temperature input.'
stop
end if
read(*,*,iostat=iostat_val) R_gas
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid gas constant input.'
stop
end if
! ------------------------------------------------------------------
! Validate inputs
! ------------------------------------------------------------------
if (M1 < 1.0d0) then
write(*,*) 'ERROR: Upstream Mach number must be >= 1.0 for a normal shock.'
stop
end if
if (gamma <= 1.0d0) then
write(*,*) 'ERROR: Specific heat ratio gamma must be > 1.0.'
stop
end if
if (p1 <= 0.0d0) then
write(*,*) 'ERROR: Upstream pressure must be positive.'
stop
end if
if (T1 <= 0.0d0) then
write(*,*) 'ERROR: Upstream temperature must be positive.'
stop
end if
if (R_gas <= 0.0d0) then
write(*,*) 'ERROR: Gas constant must be positive.'
stop
end if
! ------------------------------------------------------------------
! Convenience variables
! ------------------------------------------------------------------
gp1 = gamma + 1.0d0
gm1 = gamma - 1.0d0
gp1h = gp1 / 2.0d0
gm1h = gm1 / 2.0d0
! ------------------------------------------------------------------
! Upstream state
! ------------------------------------------------------------------
rho1 = p1 / (R_gas * T1)
a1 = sqrt(gamma * R_gas * T1)
u1 = M1 * a1
! ------------------------------------------------------------------
! Normal shock ratios
! ------------------------------------------------------------------
! Static pressure ratio p2/p1
p_ratio = 1.0d0 + (2.0d0 * gamma / gp1) * (M1**2 - 1.0d0)
! Density ratio rho2/rho1 (= u1/u2 from continuity)
rho_ratio = (gp1 * M1**2) / (gm1 * M1**2 + 2.0d0)
! Temperature ratio T2/T1
T_ratio = p_ratio / rho_ratio
! Downstream Mach number squared
M2_sq = (gm1 * M1**2 + 2.0d0) / (2.0d0 * gamma * M1**2 - gm1)
M2 = sqrt(M2_sq)
! Velocity ratio u2/u1 (inverse of density ratio from continuity)
u_ratio = 1.0d0 / rho_ratio
! Total pressure ratio p02/p01
p0_ratio = (rho_ratio ** (gamma / gm1)) &
* (p_ratio ** (-1.0d0 / gm1))
! Entropy change Delta_s / R
ds_over_R = -log(p0_ratio)
! ------------------------------------------------------------------
! Absolute downstream state
! ------------------------------------------------------------------
p2 = p_ratio * p1
T2 = T_ratio * T1
rho2 = rho_ratio * rho1
a2 = sqrt(gamma * R_gas * T2)
u2 = M2 * a2
! ------------------------------------------------------------------
! Stagnation (total) conditions
! ------------------------------------------------------------------
p01 = p1 * (1.0d0 + gm1h * M1**2) ** (gamma / gm1)
T01 = T1 * (1.0d0 + gm1h * M1**2)
p02 = p0_ratio * p01
T02 = T01 ! Total temperature is constant across normal shock
! ------------------------------------------------------------------
! Print results
! ------------------------------------------------------------------
write(*,'(A)') '============================================================'
write(*,'(A)') ' NORMAL SHOCK RELATIONS (Rankine-Hugoniot)'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- UPSTREAM CONDITIONS -------------------------------------'
write(*,'(A,F12.4)') ' Upstream Mach (M1) = ', M1
write(*,'(A,F12.6)') ' Specific Heat Ratio (g) = ', gamma
write(*,'(A,F12.2,A)') ' Gas Constant (R) = ', R_gas, ' J/kg.K'
write(*,'(A,ES14.6,A)') ' Upstream Pressure (p1) = ', p1, ' Pa'
write(*,'(A,F12.2,A)') ' Upstream Temperature(T1)= ', T1, ' K'
write(*,'(A,F12.4,A)') ' Upstream Density (rho1) = ', rho1, ' kg/m3'
write(*,'(A,F12.4,A)') ' Speed of Sound (a1) = ', a1, ' m/s'
write(*,'(A,F12.4,A)') ' Flow Velocity (u1) = ', u1, ' m/s'
write(*,*)
write(*,'(A)') '--- STAGNATION (TOTAL) CONDITIONS ---------------------------'
write(*,'(A,ES14.6,A)') ' Total Pressure p01 = ', p01, ' Pa'
write(*,'(A,F12.2,A)') ' Total Temperature T01 = ', T01, ' K'
write(*,*)
write(*,'(A)') '--- NORMAL SHOCK RATIOS -------------------------------------'
write(*,'(A,F12.6)') ' p2/p1 (pressure) = ', p_ratio
write(*,'(A,F12.6)') ' rho2/rho1 (density) = ', rho_ratio
write(*,'(A,F12.6)') ' T2/T1 (temperature) = ', T_ratio
write(*,'(A,F12.6)') ' u2/u1 (velocity) = ', u_ratio
write(*,'(A,F12.6)') ' M2 (downstream Mach)= ', M2
write(*,'(A,F12.8)') ' p02/p01 (total pressure)= ', p0_ratio
write(*,'(A,F12.8)') ' Ds/R (entropy rise) = ', ds_over_R
write(*,*)
write(*,'(A)') '--- DOWNSTREAM CONDITIONS -----------------------------------'
write(*,'(A,ES14.6,A)') ' Downstream Pressure p2 = ', p2, ' Pa'
write(*,'(A,F12.2,A)') ' Downstream Temp T2 = ', T2, ' K'
write(*,'(A,F12.4,A)') ' Downstream Density rho2 = ', rho2, ' kg/m3'
write(*,'(A,F12.4,A)') ' Speed of Sound a2 = ', a2, ' m/s'
write(*,'(A,F12.4,A)') ' Flow Velocity u2 = ', u2, ' m/s'
write(*,'(A,F12.6)') ' Downstream Mach M2 = ', M2
write(*,'(A,ES14.6,A)') ' Total Pressure p02 = ', p02, ' Pa'
write(*,'(A,F12.2,A)') ' Total Temperature T02 = ', T02, ' K'
write(*,*)
! ------------------------------------------------------------------
! Profile sweep: M1 from 1.0 to user M1 in 40 steps
! ------------------------------------------------------------------
write(*,'(A)') '--- SHOCK PROPERTY PROFILE vs MACH NUMBER -------------------'
write(*,'(A)') ' M1 p2/p1 rho2/rho1 T2/T1 M2' // &
' p02/p01 Ds/R'
write(*,'(A)') ' ----------------------------------------------------------' // &
'----------------------------'
n_points = 40
if (M1 <= 1.0d0) then
dM = 0.0d0
else
dM = (M1 - 1.0d0) / dble(n_points)
end if
do i = 0, n_points
M_cur = 1.0d0 + dble(i) * dM
if (M_cur < 1.0d0) M_cur = 1.0d0
pr_cur = 1.0d0 + (2.0d0 * gamma / gp1) * (M_cur**2 - 1.0d0)
rr_cur = (gp1 * M_cur**2) / (gm1 * M_cur**2 + 2.0d0)
Tr_cur = pr_cur / rr_cur
M2s_cur = (gm1 * M_cur**2 + 2.0d0) / (2.0d0 * gamma * M_cur**2 - gm1)
M2_cur = sqrt(M2s_cur)
p0r_cur = (rr_cur ** (gamma / gm1)) * (pr_cur ** (-1.0d0 / gm1))
ds_cur = -log(p0r_cur)
write(*,'(F8.4,2X,F11.6,2X,F11.6,2X,F11.6,2X,F11.6,2X,F11.8,2X,F11.8)') &
M_cur, pr_cur, rr_cur, Tr_cur, M2_cur, p0r_cur, ds_cur
end do
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
write(*,'(A)') ' p2/p1 = 1 + 2*gamma/(gamma+1) * (M1^2 - 1)'
write(*,'(A)') ' rho2/rho1 = (gamma+1)*M1^2 / ((gamma-1)*M1^2 + 2)'
write(*,'(A)') ' T2/T1 = (p2/p1) / (rho2/rho1)'
write(*,'(A)') ' M2^2 = ((gamma-1)*M1^2 + 2) / (2*gamma*M1^2 - (gamma-1))'
write(*,'(A)') ' p02/p01 = (rho2/rho1)^(g/(g-1)) * (p2/p1)^(-1/(g-1))'
write(*,'(A)') ' Ds/R = -ln(p02/p01)'
write(*,'(A)') '============================================================'
end program rankine_hugoniot
Solver Description
Compute normal shock wave property jumps (pressure, density, temperature) using Rankine-Hugoniot relations.
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 rankine_hugoniot.f90 -o rankine_hugoniot
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
rankine_hugoniot < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Upstream Mach Number ($M_1$)\nSpecific Heat Ratio ($\gamma$)\nUpstream Pressure ($p_1$) [Pa]\nUpstream Temperature ($T_1$) [K]\nGas Preset
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0