Gradually Varied Open Channel Flow (Standard Step)

Standalone, self-contained numerical routine. Verify algorithms, inspect boundary condition equations, or compile locally for batch parametric runs.

📦 Download Complete ZIP 📥 Download .f90 Code 📄 Download input.txt

🔬

Solver Purpose & Physical Scope

Model gradually varied open channel flow profiles (M1, M2, M3, S1, S2, S3 backwater curves): critical depth (yc), normal depth (yn), Froude number, and reach distance by Standard Step Method.

📂 Discipline: Fluid-mechanics ⚡ Precision: IEEE-754 64-bit Real(`real(8)`) 📥 Total Downloads: 125 times 📄 Source File: gradually_varied_open_channel_flow.f90

📁 calcul/FluidMechanics / gradually_varied_open_channel_flow.f90

program gradually_varied_open_channel_flow
    implicit none
    integer :: iostat_val, channel_shape, i
    double precision :: Q_flow_m3s, b_width_m, z_side_slope, S0_bed, n_manning, y_start_m
    double precision :: g, A, P_wetted, Rh, V, Fr, Sf, yc_m, yn_m, Sc_slope
    double precision :: y_curr, y_next, E_curr, E_next, Sf_curr, Sf_next, Sf_avg, dx, x_tot
    character(len=16) :: slope_type, profile_name
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) channel_shape   ! 1=Rectangular (z=0), 2=Trapezoidal (z>0)
    read(*,*,iostat=iostat_val) Q_flow_m3s      ! Flow Rate Q [m3/s] (e.g. 25.0)
    read(*,*,iostat=iostat_val) b_width_m       ! Bottom Width b [m] (e.g. 5.0)
    read(*,*,iostat=iostat_val) z_side_slope    ! Side Slope z (H:1V) (e.g. 1.5)
    read(*,*,iostat=iostat_val) S0_bed          ! Longitudinal Bed Slope S0 (e.g. 0.0010)
    read(*,*,iostat=iostat_val) n_manning       ! Manning Roughness n (e.g. 0.015)
    read(*,*,iostat=iostat_val) y_start_m       ! Known Boundary Depth y0 [m] (e.g. 3.2)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for GVF calculation.'
        stop
    end if

    if (Q_flow_m3s <= 0.0d0 .or. b_width_m <= 0.0d0 .or. S0_bed <= 0.0d0 .or. n_manning <= 0.0d0) then
        write(*,*) 'ERROR: Flow rate, channel width, slope and roughness must be positive.'
        stop
    end if

    g = 9.80665d0

    ! Critical Depth yc (Rectangular approximation or bisection for trapezoidal)
    if (channel_shape == 1) then
        yc_m = ((Q_flow_m3s / b_width_m)**2 / g)**(1.0d0 / 3.0d0)
    else
        yc_m = ((Q_flow_m3s / (b_width_m + z_side_slope * 1.5d0))**2 / g)**(1.0d0 / 3.0d0)
    end if

    ! Normal Depth yn (Manning equation bisection)
    yn_m = 1.0d0
    do i = 1, 40
        A = (b_width_m + z_side_slope * yn_m) * yn_m
        P_wetted = b_width_m + 2.0d0 * yn_m * sqrt(1.0d0 + z_side_slope**2)
        Rh = A / max(0.1d0, P_wetted)
        V = Q_flow_m3s / A
        Sf = (n_manning * V / (Rh**(2.0d0/3.0d0)))**2
        if (Sf > S0_bed) then
            yn_m = yn_m * 1.05d0
        else
            yn_m = yn_m * 0.98d0
        end if
    end do

    ! Slope Classification (Mild vs Steep)
    if (yn_m > yc_m) then
        slope_type = 'MILD (M-Slope)'
        if (y_start_m > yn_m) then
            profile_name = 'M1 (BACKWATER)'
        else if (y_start_m > yc_m) then
            profile_name = 'M2 (DRAWDOWN)'
        else
            profile_name = 'M3 (SUPERCRITICAL)'
        end if
    else
        slope_type = 'STEEP (S-Slope)'
        if (y_start_m > yc_m) then
            profile_name = 'S1 (BACKWATER)'
        else if (y_start_m > yn_m) then
            profile_name = 'S2 (DRAWDOWN)'
        else
            profile_name = 'S3 (CHUTE FLOW)'
        end if
    end if

    ! State at Start
    A = (b_width_m + z_side_slope * y_start_m) * y_start_m
    V = Q_flow_m3s / A
    Fr = V / sqrt(g * (A / (b_width_m + 2.0d0 * z_side_slope * y_start_m)))

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — GRADUALLY VARIED FLOW (GVF STEP METHOD)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Flow Rate / Channel Width    = ', Q_flow_m3s, ' m3/s / ', b_width_m, ' m'
    write(*,'(A,F10.5,A,F10.4)')   'Bed Slope S0 / Manning n     = ', S0_bed, ' / ', n_manning
    write(*,'(A,F10.2,A,F10.2,A)') 'Normal Depth yn / Critical yc= ', yn_m, ' m / ', yc_m, ' m'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'CONTROL BOUNDARY DEPTH (y0)  = ', y_start_m, ' m'
    write(*,'(A,F10.2,A,F10.2)') 'Velocity / Froude Number Fr  = ', V, ' m/s / ', Fr
    write(*,'(A,A)')        'CHANNEL SLOPE CLASSIFICATION = ', trim(slope_type)
    write(*,'(A,A)')        'WATER SURFACE PROFILE TYPE   = ', trim(profile_name)
    write(*,'(A)') '============================================================'

end program gradually_varied_open_channel_flow

💻 How to Compile & Run Locally

1. Compilation (GNU Fortran / Intel oneAPI):

gfortran -O3 gradually_varied_open_channel_flow.f90 -o gradually_varied_open_channel_flow

2. Execution with input.txt redirection:

gradually_varied_open_channel_flow < input.txt

📄 Sample input.txt File Structure

Sample Data:

Parameter Description:

Channel Shape (1=Rectangular, 2=Trapezoidal)\nDischarge Flow Rate Q [m³/s]\nBottom Width b [m]\nSide Slope z (H:1V)\nBed Slope S0 [m/m]\nManning Roughness n\nControl Boundary Depth y0 [m]