Fortran Source Code & Offline Testing

Surge Tank & Surge Analysis

Core Numerical Engine in Fortran 90 • 31 total downloads

surge_tank.f90

! =========================================================================
! Source File: surge_tank.f90
! =========================================================================

program surge_tank
    implicit none
    integer :: i, n_steps, iostat_val
    double precision :: L_tunnel, A_tunnel, A_tank, f_darcy, D_tunnel
    double precision :: Q0, H0, g, t_end, dt
    double precision :: z, Vt, Q_t, z_max, z_min, t_max, t_min, period_est
    double precision :: z_prev, Vt_prev, dz_dt, dVt_dt, Sf
    double precision :: z_k1, Vt_k1, z_k2, Vt_k2, z_k3, Vt_k3, z_k4, Vt_k4
    double precision :: thoma_crit, thoma_actual, z_first_peak
    double precision :: t
    integer :: peak_count, first_peak_found
    character(len=40) :: stability

    read(*,*,iostat=iostat_val) L_tunnel
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid tunnel length input.'
        stop
    end if
    read(*,*,iostat=iostat_val) A_tunnel
    read(*,*,iostat=iostat_val) D_tunnel
    read(*,*,iostat=iostat_val) A_tank
    read(*,*,iostat=iostat_val) f_darcy
    read(*,*,iostat=iostat_val) Q0
    read(*,*,iostat=iostat_val) H0
    read(*,*,iostat=iostat_val) g
    read(*,*,iostat=iostat_val) t_end
    read(*,*,iostat=iostat_val) n_steps
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all surge tank inputs.'
        stop
    end if
    if (L_tunnel<=0.0d0.or.A_tunnel<=0.0d0.or.A_tank<=0.0d0.or.D_tunnel<=0.0d0) then
        write(*,*) 'ERROR: Tunnel and tank dimensions must be positive.'
        stop
    end if
    if (Q0<0.0d0.or.f_darcy<0.0d0) then
        write(*,*) 'ERROR: Q0 and friction factor must be non-negative.'
        stop
    end if
    if (g<=0.0d0) g = 9.81d0
    if (t_end<=0.0d0) t_end = 600.0d0
    if (n_steps<100) n_steps = 2000

    dt = t_end / dble(n_steps)

    ! Thoma stability criterion
    ! A_tank_critical = f L V0^2 / (2 g H0)
    ! where V0 = Q0/A_tunnel
    if (H0 > 0.0d0) then
        thoma_crit = f_darcy * L_tunnel * (Q0/A_tunnel)**2 / (2.0d0*g*H0)
        thoma_actual = A_tank
    else
        thoma_crit = 0.0d0
        thoma_actual = A_tank
    end if
    if (A_tank > thoma_crit) then
        stability = 'STABLE (A_tank > A_Thoma)'
    else
        stability = 'UNSTABLE (A_tank < A_Thoma)'
    end if

    ! Natural oscillation period estimate
    ! T = 2*pi*sqrt(L*A_tank/(g*A_tunnel))
    period_est = 2.0d0 * 3.141592653589793d0 * &
                 sqrt(L_tunnel*A_tank/(g*A_tunnel))

    ! Initial conditions: sudden full load rejection (Q drops to 0)
    ! z = water level deviation in surge tank from steady state
    ! Vt = tunnel velocity
    z = 0.0d0
    Vt = Q0 / A_tunnel
    z_max = 0.0d0
    z_min = 0.0d0
    t_max = 0.0d0
    t_min = 0.0d0
    peak_count = 0
    first_peak_found = 0
    z_first_peak = 0.0d0

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   SURGE TANK & SURGE ANALYSIS ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Tunnel Length             = ', L_tunnel, ' m'
    write(*,'(A,ES12.4,A)') '  Tunnel Area               = ', A_tunnel, ' m2'
    write(*,'(A,ES12.4,A)') '  Tunnel Diameter           = ', D_tunnel, ' m'
    write(*,'(A,ES12.4,A)') '  Tank Area                 = ', A_tank, ' m2'
    write(*,'(A,ES12.4)')   '  Darcy Friction Factor     = ', f_darcy
    write(*,'(A,ES12.4,A)') '  Steady-State Flow Q0      = ', Q0, ' m3/s'
    write(*,'(A,ES12.4,A)') '  Net Head H0               = ', H0, ' m'
    write(*,'(A,ES12.4,A)') '  Simulation Time           = ', t_end, ' s'
    write(*,*)
    write(*,'(A)') '--- STABILITY & PERIOD --------------------------------------'
    write(*,'(A,ES12.4,A)') '  Thoma Critical Area       = ', thoma_crit, ' m2'
    write(*,'(A,ES12.4,A)') '  Actual Tank Area          = ', A_tank, ' m2'
    write(*,'(A,A)')        '  Stability Assessment      = ', trim(stability)
    write(*,'(A,ES12.4,A)') '  Natural Period Estimate   = ', period_est, ' s'
    write(*,*)

    write(*,'(A)') '--- SURGE OSCILLATION PROFILE -------------------------------'
    write(*,'(A)') '  t[s]          z[m]          Vt[m/s]       Q_tunnel[m3/s]'
    write(*,'(A)') '  ----------------------------------------------------------------'

    do i = 0, n_steps
        t = dble(i) * dt
        Q_t = Vt * A_tunnel

        ! Output every n_steps/200 steps or first/last
        if (mod(i, max(n_steps/200,1)) == 0 .or. i == n_steps) then
            write(*,'(F12.3,2X,F12.6,2X,F12.6,2X,ES12.4)') t, z, Vt, Q_t
        end if

        ! Track extremes
        if (z > z_max) then
            z_max = z; t_max = t
        end if
        if (z < z_min) then
            z_min = z; t_min = t
        end if

        ! Detect first peak (dz changes sign from + to -)
        if (i > 0 .and. first_peak_found == 0) then
            if (z < z_prev .and. z_prev > 0.01d0) then
                first_peak_found = 1
                z_first_peak = z_prev
            end if
        end if

        if (i == n_steps) exit

        z_prev = z
        Vt_prev = Vt

        ! RK4 integration of surge tank equations
        ! dz/dt = (A_tunnel/A_tank) * Vt
        ! dVt/dt = (g/L) * (H0 - z - f*L*Vt*|Vt|/(2*g*D))
        ! After load rejection: driving head = -z (simplified)
        call surge_rhs(z, Vt, z_k1, Vt_k1)
        call surge_rhs(z+0.5d0*dt*z_k1, Vt+0.5d0*dt*Vt_k1, z_k2, Vt_k2)
        call surge_rhs(z+0.5d0*dt*z_k2, Vt+0.5d0*dt*Vt_k2, z_k3, Vt_k3)
        call surge_rhs(z+dt*z_k3, Vt+dt*Vt_k3, z_k4, Vt_k4)

        z  = z  + dt/6.0d0*(z_k1  + 2.0d0*z_k2  + 2.0d0*z_k3  + z_k4)
        Vt = Vt + dt/6.0d0*(Vt_k1 + 2.0d0*Vt_k2 + 2.0d0*Vt_k3 + Vt_k4)
    end do

    write(*,*)
    write(*,'(A)') '--- SURGE RESULTS -------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Maximum Surge zmax        = ', z_max, ' m'
    write(*,'(A,ES12.4,A)') '  Time of zmax              = ', t_max, ' s'
    write(*,'(A,ES12.4,A)') '  Minimum Surge zmin        = ', z_min, ' m'
    write(*,'(A,ES12.4,A)') '  Time of zmin              = ', t_min, ' s'
    write(*,'(A,ES12.4,A)') '  Oscillation Period        = ', period_est, ' s'
    write(*,'(A,A)')        '  Stability                 = ', trim(stability)
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  dz/dt = (A_tunnel/A_tank) Vt.'
    write(*,'(A)') '  dVt/dt = (g/L)(−z − f L Vt|Vt|/(2gD)).'
    write(*,'(A)') '  Thoma criterion: A_crit = f L V0^2/(2gH0).'
    write(*,'(A)') '  Period: T = 2 pi sqrt(L A_tank/(g A_tunnel)).'
    write(*,'(A)') '  Integration: 4th-order Runge-Kutta.'

contains

    subroutine surge_rhs(z_in, Vt_in, dz_out, dVt_out)
        implicit none
        double precision, intent(in) :: z_in, Vt_in
        double precision, intent(out) :: dz_out, dVt_out
        double precision :: friction_head
        dz_out = (A_tunnel / A_tank) * Vt_in
        friction_head = f_darcy * L_tunnel * Vt_in * abs(Vt_in) / &
                        (2.0d0 * g * D_tunnel)
        dVt_out = (g / L_tunnel) * (-z_in - friction_head)
    end subroutine surge_rhs

end program surge_tank

Solver Description

Simulate water level oscillations in surge tanks after sudden load rejection using RK4 integration. Evaluates Thoma stability.

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

Execution Command:

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

surge_tank < input.txt

📥 Downloads & Local Files

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

! Tunnel length L [m]\nTunnel area At [mÂ²]\nTunnel diameter Dt [m]\nSurge tank area As [mÂ²]\nDarcy friction factor f\nSteady-state flow Qâ‚€ [mÂ³/s]\nNet head Hâ‚€ [m]\nGravity g [m/sÂ²]\nNumber of steps\nns_init
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0
! Parameter 6
0.0
! Parameter 7
0.0
! Parameter 8
0.0
! Parameter 9
0.0
! Parameter 10
0.0

📦 Download Offline Bundle (ZIP) 13 dl 📄 Download input.txt 11 dl 📥 Download surge_tank.f90 7 dl

💻 Fortran Source Code Library

Surge Tank & Surge Analysis

Solver Description

Key Numerical Methods & Architecture

🛠️ Local Compilation

📥 Downloads & Local Files