Non-Newtonian Yield-Pseudoplastic Hydraulics (Herschel-Bulkley)

program non_newtonian_herschel_bulkley
    implicit none
    integer :: iostat_val, fluid_type, iter
    double precision :: tau_yield_Pa, K_consist, n_index, rho_kgm3, D_pipe_mm, L_pipe_m, Q_flow_m3h
    double precision :: D_m, R_m, A_pipe, V_mean, tau_w, tau_w_low, tau_w_high, Q_calc, f_error, df_dtau
    double precision :: a_ratio, rp_mm, dp_pipe_kPa, Re_gen, mu_app_wall, wall_shear_rate
    character(len=32) :: flow_regime
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) fluid_type      ! 1=Drilling Mud, 2=Mineral Slurry, 3=Toothpaste/Gel, 4=Tomato Paste / Food
    read(*,*,iostat=iostat_val) tau_yield_Pa    ! Yield Stress tau_y [Pa] (e.g. 12.0)
    read(*,*,iostat=iostat_val) K_consist       ! Flow Consistency K [Pa.s^n] (e.g. 1.80)
    read(*,*,iostat=iostat_val) n_index         ! Flow Behavior Index n (e.g. 0.65)
    read(*,*,iostat=iostat_val) rho_kgm3        ! Fluid Density [kg/m3] (e.g. 1250.0)
    read(*,*,iostat=iostat_val) D_pipe_mm       ! Pipe Inner Diameter [mm] (e.g. 80.0)
    read(*,*,iostat=iostat_val) L_pipe_m        ! Pipe Length [m] (e.g. 100.0)
    read(*,*,iostat=iostat_val) Q_flow_m3h      ! Target Flow Rate [m3/h] (e.g. 24.0)

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

    if (D_pipe_mm <= 0.0d0 .or. Q_flow_m3h <= 0.0d0 .or. K_consist <= 0.0d0 .or. n_index <= 0.0d0) then
        write(*,*) 'ERROR: Diameter, flow rate, consistency and n must be positive.'
        stop
    end if

    D_m = D_pipe_mm / 1000.0d0
    R_m = D_m / 2.0d0
    A_pipe = PI * (R_m**2)
    V_mean = (Q_flow_m3h / 3600.0d0) / A_pipe

    ! Solve for Wall Shear Stress tau_w using Newton-Raphson / Bisection
    tau_w_low = tau_yield_Pa * 1.0001d0
    tau_w_high = max(100.0d0, tau_yield_Pa * 50.0d0)

    ! Bracket upper bound
    do iter = 1, 100
        a_ratio = tau_yield_Pa / tau_w_high
        Q_calc = PI * (R_m**3) * ((tau_w_high / K_consist)**(1.0d0 / n_index)) * &
                 ((1.0d0 - a_ratio)**(1.0d0 + 1.0d0/n_index)) * &
                 ( ((1.0d0 - a_ratio)**2) / (3.0d0 + 1.0d0/n_index) + &
                   (2.0d0 * a_ratio * (1.0d0 - a_ratio)) / (2.0d0 + 1.0d0/n_index) + &
                   (a_ratio**2) / (1.0d0 + 1.0d0/n_index) )
        if (Q_calc >= (Q_flow_m3h / 3600.0d0)) exit
        tau_w_high = tau_w_high * 2.0d0
    end do

    ! Bisection for tau_w
    do iter = 1, 60
        tau_w = 0.5d0 * (tau_w_low + tau_w_high)
        a_ratio = tau_yield_Pa / tau_w
        Q_calc = PI * (R_m**3) * ((tau_w / K_consist)**(1.0d0 / n_index)) * &
                 ((1.0d0 - a_ratio)**(1.0d0 + 1.0d0/n_index)) * &
                 ( ((1.0d0 - a_ratio)**2) / (3.0d0 + 1.0d0/n_index) + &
                   (2.0d0 * a_ratio * (1.0d0 - a_ratio)) / (2.0d0 + 1.0d0/n_index) + &
                   (a_ratio**2) / (1.0d0 + 1.0d0/n_index) )
        if (Q_calc < (Q_flow_m3h / 3600.0d0)) then
            tau_w_low = tau_w
        else
            tau_w_high = tau_w
        end if
    end do

    ! Plug radius (rp = R * a_ratio)
    rp_mm = (R_m * a_ratio) * 1000.0d0

    ! Pressure drop Delta P [kPa]
    dp_pipe_kPa = (4.0d0 * tau_w * L_pipe_m / D_m) / 1000.0d0

    ! Wall shear rate & Apparent viscosity
    wall_shear_rate = ((max(0.0d0, tau_w - tau_yield_Pa)) / K_consist)**(1.0d0 / n_index)
    mu_app_wall = tau_w / max(1.0d-3, wall_shear_rate)

    ! Generalized Reynolds number (Metzner-Reed style)
    Re_gen = (rho_kgm3 * (V_mean**(2.0d0 - n_index)) * (D_m**n_index)) / &
             (K_consist * ((3.0d0*n_index + 1.0d0) / (4.0d0*n_index))**n_index * 8.0d0**(n_index - 1.0d0))

    if (Re_gen < 2100.0d0) then
        flow_regime = 'LAMINAR (PLUG CORE PRESENT)'
    else if (Re_gen < 3500.0d0) then
        flow_regime = 'TRANSITIONAL YIELD FLOW'
    else
        flow_regime = 'TURBULENT NON-NEWTONIAN'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — NON-NEWTONIAN HERSCHEL-BULKLEY HYDRAULICS'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Yield Stress / Consistency K = ', tau_yield_Pa, ' Pa / ', K_consist, ' Pa.s^n'
    write(*,'(A,F10.3,A,F10.2,A)') 'Behavior Index n / Flow Rate = ', n_index, ' / ', Q_flow_m3h, ' m3/h'
    write(*,'(A,F10.1,A,F10.1,A)') 'Pipe ID / Pipe Length        = ', D_pipe_mm, ' mm / ', L_pipe_m, ' m'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'TOTAL PIPE PRESSURE DROP (ΔP)= ', dp_pipe_kPa, ' kPa'
    write(*,'(A,F10.2,A)')  'Wall Shear Stress (tau_w)    = ', tau_w, ' Pa'
    write(*,'(A,F10.2,A)')  'Unyielded Plug Radius (rp)   = ', rp_mm, ' mm'
    write(*,'(A,F10.2,A)')  'Plug Area Fraction           = ', (rp_mm / (D_pipe_mm/2.0d0))**2 * 100.0d0, ' %'
    write(*,'(A,F10.1)')    'Generalized Reynolds (Re_gen)= ', Re_gen
    write(*,'(A,A)')        'PREDICTED FLOW REGIME        = ', trim(flow_regime)
    write(*,'(A)') '============================================================'

end program non_newtonian_herschel_bulkley