๐Ÿ” Search & Filter All 178+ Engineering Calculators โ†’

๐Ÿ’ป 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.

Non-Newtonian Pipe Flow

Core Numerical Engine in Fortran 90 โ€ข 22 total downloads

non_newtonian_pipe.f90 
! =========================================================================
! Source File: non_newtonian_pipe.f90
! =========================================================================

program non_newtonian_pipe
    implicit none
    integer :: model, i, n_profile, iostat_val
    double precision :: D, L_pipe, Q, rho, n_pl, K_pl, tau_y
    double precision :: R, A_pipe, V_bulk, mu_app, Re_gen, f, dP, tau_wall
    double precision :: V_max, r_plug, r_i, V_i, gamma_wall, Re_MR
    double precision :: Q_i, dP_i, V_i2, Re_i, f_i, mu_i
    double precision, parameter :: PI = 3.141592653589793d0

    read(*,*,iostat=iostat_val) model
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid model input.'
        stop
    end if
    read(*,*,iostat=iostat_val) D
    read(*,*,iostat=iostat_val) L_pipe
    read(*,*,iostat=iostat_val) Q
    read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) n_pl
    read(*,*,iostat=iostat_val) K_pl
    read(*,*,iostat=iostat_val) tau_y
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all non-Newtonian pipe flow inputs.'
        stop
    end if
    if (D <= 0.0d0 .or. L_pipe <= 0.0d0 .or. Q <= 0.0d0) then
        write(*,*) 'ERROR: Diameter, length and flow rate must be positive.'
        stop
    end if
    if (rho <= 0.0d0 .or. K_pl <= 0.0d0 .or. n_pl <= 0.0d0) then
        write(*,*) 'ERROR: Density, K and n must be positive.'
        stop
    end if
    if (tau_y < 0.0d0) tau_y = 0.0d0

    R = D / 2.0d0
    A_pipe = PI * R**2
    V_bulk = Q / A_pipe

    select case(model)
    case(1)
        ! Power-law model: tau = K * gamma_dot^n
        call power_law_pipe(R, V_bulk, rho, n_pl, K_pl, &
            mu_app, Re_gen, Re_MR, f, dP, tau_wall, gamma_wall, V_max)
        r_plug = 0.0d0
    case(2)
        ! Bingham plastic: tau = tau_y + K * gamma_dot
        call bingham_pipe(R, V_bulk, rho, K_pl, tau_y, &
            mu_app, Re_gen, Re_MR, f, dP, tau_wall, gamma_wall, V_max, r_plug)
    case(3)
        ! Herschel-Bulkley: tau = tau_y + K * gamma_dot^n
        call herschel_bulkley_pipe(R, V_bulk, rho, n_pl, K_pl, tau_y, &
            mu_app, Re_gen, Re_MR, f, dP, tau_wall, gamma_wall, V_max, r_plug)
    case default
        write(*,*) 'ERROR: Model must be 1 power-law, 2 Bingham, or 3 Herschel-Bulkley.'
        stop
    end select

    dP = f * (L_pipe/D) * 0.5d0 * rho * V_bulk**2

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   NON-NEWTONIAN PIPE FLOW ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,I8)')       '  Rheological Model        = ', model
    write(*,'(A,ES12.4,A)') '  Pipe Diameter            = ', D, ' m'
    write(*,'(A,ES12.4,A)') '  Pipe Length              = ', L_pipe, ' m'
    write(*,'(A,ES12.4,A)') '  Volume Flow Rate         = ', Q, ' m3/s'
    write(*,'(A,ES12.4,A)') '  Fluid Density            = ', rho, ' kg/m3'
    write(*,'(A,ES12.4)')   '  Flow Index n             = ', n_pl
    write(*,'(A,ES12.4,A)') '  Consistency K            = ', K_pl, ' Pa.s^n'
    write(*,'(A,ES12.4,A)') '  Yield Stress tau_y       = ', tau_y, ' Pa'
    write(*,*)
    write(*,'(A)') '--- FLOW RESULTS --------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Bulk Velocity            = ', V_bulk, ' m/s'
    write(*,'(A,ES12.4,A)') '  Wall Shear Stress        = ', tau_wall, ' Pa'
    write(*,'(A,ES12.4,A)') '  Wall Shear Rate          = ', gamma_wall, ' 1/s'
    write(*,'(A,ES12.4,A)') '  Apparent Viscosity       = ', mu_app, ' Pa.s'
    write(*,'(A,ES12.4)')   '  Generalized Re           = ', Re_gen
    write(*,'(A,ES12.4)')   '  Metzner-Reed Re          = ', Re_MR
    write(*,'(A,ES12.4)')   '  Fanning Friction Factor  = ', f/4.0d0
    write(*,'(A,ES12.4)')   '  Darcy Friction Factor    = ', f
    write(*,'(A,ES12.4,A)') '  Pressure Drop            = ', dP, ' Pa'
    write(*,'(A,ES12.4,A)') '  Maximum Velocity         = ', V_max, ' m/s'
    if (r_plug > 0.0d0) then
        write(*,'(A,ES12.4,A)') '  Plug Radius              = ', r_plug, ' m'
        write(*,'(A,F10.3,A)')  '  Plug Fraction r_plug/R   = ', r_plug/R*100.0d0, ' percent'
    end if
    write(*,*)

    ! Velocity profile
    n_profile = 50
    write(*,'(A)') '--- VELOCITY PROFILE ----------------------------------------'
    write(*,'(A)') '  r/R           V(r)/V_max    r[m]          V[m/s]'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 0, n_profile
        r_i = R * dble(i) / dble(n_profile)
        call velocity_at_r(model, r_i, R, V_bulk, V_max, n_pl, &
                           tau_y, tau_wall, r_plug, V_i)
        write(*,'(F10.5,2X,F12.6,2X,ES12.4,2X,ES12.4)') &
            r_i/R, V_i/max(V_max,1.0d-30), r_i, V_i
    end do
    write(*,*)

    ! Q vs dP parametric sweep
    write(*,'(A)') '--- Q VS DELTAP SWEEP ---------------------------------------'
    write(*,'(A)') '  Q[m3/s]       V[m/s]        Re_gen        dP[Pa]        mu_app[Pa.s]'
    write(*,'(A)') '  -----------------------------------------------------------------------'
    do i = 1, 40
        Q_i = Q * 0.1d0 * dble(i)
        V_i2 = Q_i / A_pipe
        select case(model)
        case(1)
            call power_law_pipe(R, V_i2, rho, n_pl, K_pl, &
                mu_i, Re_i, Re_MR, f_i, dP_i, tau_wall, gamma_wall, V_max)
        case(2)
            call bingham_pipe(R, V_i2, rho, K_pl, tau_y, &
                mu_i, Re_i, Re_MR, f_i, dP_i, tau_wall, gamma_wall, V_max, r_plug)
        case(3)
            call herschel_bulkley_pipe(R, V_i2, rho, n_pl, K_pl, tau_y, &
                mu_i, Re_i, Re_MR, f_i, dP_i, tau_wall, gamma_wall, V_max, r_plug)
        end select
        dP_i = f_i * (L_pipe/D) * 0.5d0 * rho * V_i2**2
        write(*,'(ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') &
            Q_i, V_i2, Re_i, dP_i, mu_i
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Power-law: tau = K gamma_dot^n.'
    write(*,'(A)') '  Bingham: tau = tau_y + K gamma_dot (n=1).'
    write(*,'(A)') '  Herschel-Bulkley: tau = tau_y + K gamma_dot^n.'
    write(*,'(A)') '  Metzner-Reed Re: Re_MR = rho V^(2-n) D^n / (K 8^(n-1) ((3n+1)/(4n))^n).'
    write(*,'(A)') '  Laminar f = 64/Re_MR (generalized).'

contains

    subroutine power_law_pipe(Rin, Vb, rhoin, nin, Kin, &
                              mu_out, Re_out, Re_MR_out, f_out, dP_out, &
                              tw_out, gw_out, Vm_out)
        implicit none
        double precision, intent(in) :: Rin, Vb, rhoin, nin, Kin
        double precision, intent(out) :: mu_out, Re_out, Re_MR_out, f_out
        double precision, intent(out) :: dP_out, tw_out, gw_out, Vm_out
        double precision :: gamma_nom, nprime, Kprime, Din

        Din = 2.0d0*Rin
        gamma_nom = 8.0d0*Vb/Din
        nprime = nin
        Kprime = Kin * ((3.0d0*nin+1.0d0)/(4.0d0*nin))**nin

        gw_out = gamma_nom * (3.0d0*nin+1.0d0)/(4.0d0*nin)
        tw_out = Kin * gw_out**nin
        mu_out = tw_out / max(gw_out, 1.0d-30)

        Re_MR_out = rhoin * Vb**(2.0d0-nin) * Din**nin / &
                    (Kprime * 8.0d0**(nin-1.0d0))
        Re_out = rhoin * Vb * Din / mu_out

        if (Re_MR_out < 2100.0d0) then
            f_out = 64.0d0 / max(Re_MR_out, 1.0d-30)
        else
            f_out = 0.316d0 / Re_MR_out**0.25d0
        end if

        dP_out = f_out * (1.0d0/Din) * 0.5d0*rhoin*Vb**2
        Vm_out = Vb * (3.0d0*nin+1.0d0) / (nin+1.0d0)
    end subroutine power_law_pipe

    subroutine bingham_pipe(Rin, Vb, rhoin, Kin, tyin, &
                            mu_out, Re_out, Re_MR_out, f_out, dP_out, &
                            tw_out, gw_out, Vm_out, rp_out)
        implicit none
        double precision, intent(in) :: Rin, Vb, rhoin, Kin, tyin
        double precision, intent(out) :: mu_out, Re_out, Re_MR_out, f_out
        double precision, intent(out) :: dP_out, tw_out, gw_out, Vm_out, rp_out
        double precision :: Din, gamma_nom, He, xi

        Din = 2.0d0*Rin
        gamma_nom = 8.0d0*Vb/Din
        gw_out = gamma_nom
        tw_out = tyin + Kin*gw_out
        mu_out = tw_out / max(gw_out, 1.0d-30)
        Re_out = rhoin * Vb * Din / mu_out
        Re_MR_out = Re_out
        He = tyin * rhoin * Din**2 / (Kin**2)

        if (Re_out < 2100.0d0) then
            f_out = 64.0d0 / max(Re_out, 1.0d-30)
        else
            f_out = 0.316d0 / Re_out**0.25d0
        end if

        ! Buckingham-Reiner correction for laminar
        xi = tyin / max(tw_out, 1.0d-30)
        rp_out = Rin * xi
        Vm_out = Vb * 2.0d0 / (1.0d0 - xi**2) * &
                 (1.0d0 - 4.0d0/3.0d0*xi + 1.0d0/3.0d0*xi**4) / &
                 max(1.0d0 - 4.0d0/3.0d0*xi + 1.0d0/3.0d0*xi**4, 1.0d-30)
        if (Vm_out < Vb) Vm_out = Vb * 2.0d0
        dP_out = f_out * (1.0d0/Din) * 0.5d0*rhoin*Vb**2
    end subroutine bingham_pipe

    subroutine herschel_bulkley_pipe(Rin, Vb, rhoin, nin, Kin, tyin, &
                                     mu_out, Re_out, Re_MR_out, f_out, dP_out, &
                                     tw_out, gw_out, Vm_out, rp_out)
        implicit none
        double precision, intent(in) :: Rin, Vb, rhoin, nin, Kin, tyin
        double precision, intent(out) :: mu_out, Re_out, Re_MR_out, f_out
        double precision, intent(out) :: dP_out, tw_out, gw_out, Vm_out, rp_out
        double precision :: Din, gamma_nom, Kprime, xi

        Din = 2.0d0*Rin
        gamma_nom = 8.0d0*Vb/Din
        gw_out = gamma_nom * (3.0d0*nin+1.0d0)/(4.0d0*nin)
        tw_out = tyin + Kin*gw_out**nin
        mu_out = tw_out / max(gw_out, 1.0d-30)

        Kprime = Kin * ((3.0d0*nin+1.0d0)/(4.0d0*nin))**nin
        Re_MR_out = rhoin * Vb**(2.0d0-nin) * Din**nin / &
                    (Kprime * 8.0d0**(nin-1.0d0))
        Re_out = rhoin * Vb * Din / mu_out

        if (Re_MR_out < 2100.0d0) then
            f_out = 64.0d0 / max(Re_MR_out, 1.0d-30)
        else
            f_out = 0.316d0 / Re_MR_out**0.25d0
        end if

        xi = tyin / max(tw_out, 1.0d-30)
        rp_out = Rin * xi
        Vm_out = Vb * (3.0d0*nin+1.0d0) / (nin+1.0d0)
        dP_out = f_out * (1.0d0/Din) * 0.5d0*rhoin*Vb**2
    end subroutine herschel_bulkley_pipe

    subroutine velocity_at_r(mod, r_pos, Rin, Vb, Vm, nin, tyin, tw, rp, Vout)
        implicit none
        integer, intent(in) :: mod
        double precision, intent(in) :: r_pos, Rin, Vb, Vm, nin, tyin, tw, rp
        double precision, intent(out) :: Vout
        double precision :: eta

        eta = r_pos / Rin
        if (mod == 1) then
            ! Power-law velocity profile
            Vout = Vm * (1.0d0 - eta**((nin+1.0d0)/nin))
        else
            ! Bingham / HB: plug core + sheared annulus
            if (r_pos <= rp) then
                Vout = Vm
            else
                Vout = Vm * ((1.0d0 - eta) / max(1.0d0 - rp/Rin, 1.0d-30)) &
                       **((nin+1.0d0)/nin)
            end if
        end if
        if (Vout < 0.0d0) Vout = 0.0d0
    end subroutine velocity_at_r

end program non_newtonian_pipe

Solver Description

Analyze pipe flow of power-law, Bingham plastic, and Herschel-Bulkley fluids. Compute Metzner-Reed Re, pressure drop, and velocity profile.

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

Execution Command:

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

non_newtonian_pipe < input.txt

๐Ÿ“ฅ Downloads & Local Files

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

! Pipe diameter D\nPipe diameter D\nPipe length L\nVolume flow Q [mร‚ยณ/s]\nDensity รย [kg/mร‚ยณ]\nFlow index n\nConsistency K [Paร‚ยทsรขยยฟ]\nty_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
๐Ÿ“ฆ Download Offline Bundle (ZIP) 6 dl ๐Ÿ“„ Download input.txt 10 dl ๐Ÿ“ฅ Download non_newtonian_pipe.f90 6 dl