πŸ’» 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.

Transonic Similarity Parameter

Core Numerical Engine in Fortran 90 β€’ 24 total downloads

transonic_similarity.f90
! =========================================================================
! Source File: 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


Solver Description

Evaluate the transonic similarity parameter (K) for thin airfoils at near-sonic flow speeds.

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 transonic_similarity.f90 -o transonic_similarity

Execution Command:

Execute the program by feeding the sample input file into the program using stdin redirection:

transonic_similarity < input.txt

πŸ“₯ Downloads & Local Files

Preview of the required input file (input.txt):

! Mach ($M_\infty$)\nγ\nThickness Ratio (Γβ€ž = t/c)\ncpinc_init\nclinc_init\nChord [m]\npΓ’Λ†ΕΎ [Pa]\nTΓ’Λ†ΕΎ [K]\nGas
0.8
! Parameter 2
1.4
! Parameter 3
0.12
! Parameter 4
-1.0
! Parameter 5
0.5
! Parameter 6
1.5
! Parameter 7
101325
! Parameter 8
300
! Parameter 9
287.058