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