πŸ” 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.

Sedimentation & Terminal Velocity

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

sedimentation.f90 
! =========================================================================
! Source File: sedimentation.f90
! =========================================================================

program sedimentation
    implicit none
    integer :: i, iostat_val, regime_code
    double precision :: dp, rho_p, rho_f, mu, g, phi_shape
    double precision :: Vt, Re_p, CD, Ar, dp_star, Vt_star
    double precision :: F_drag, F_buoy, F_grav, W_net, tau_relax
    double precision :: dp_i, Vt_i, Re_i, CD_i, Ar_i
    double precision :: conc, Vt_hindered, n_RZ
    double precision, parameter :: PI = 3.141592653589793d0
    character(len=60) :: regime_name

    read(*,*,iostat=iostat_val) dp
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid particle diameter input.'
        stop
    end if
    read(*,*,iostat=iostat_val) rho_p
    read(*,*,iostat=iostat_val) rho_f
    read(*,*,iostat=iostat_val) mu
    read(*,*,iostat=iostat_val) g
    read(*,*,iostat=iostat_val) phi_shape
    read(*,*,iostat=iostat_val) conc
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all sedimentation inputs.'
        stop
    end if
    if (dp <= 0.0d0) then
        write(*,*) 'ERROR: Particle diameter must be positive.'
        stop
    end if
    if (rho_p <= 0.0d0 .or. rho_f <= 0.0d0 .or. mu <= 0.0d0) then
        write(*,*) 'ERROR: Densities and viscosity must be positive.'
        stop
    end if
    if (g <= 0.0d0) g = 9.81d0
    if (phi_shape <= 0.0d0 .or. phi_shape > 1.0d0) phi_shape = 1.0d0
    if (conc < 0.0d0 .or. conc >= 1.0d0) conc = 0.0d0

    ! Archimedes number
    Ar = rho_f * (rho_p - rho_f) * g * dp**3 / mu**2

    ! Dimensionless particle diameter
    dp_star = dp * (rho_f*(rho_p - rho_f)*g / mu**2)**(1.0d0/3.0d0)

    ! Terminal velocity by iterative drag balance
    call compute_terminal_velocity(dp, rho_p, rho_f, mu, g, phi_shape, &
                                   Vt, Re_p, CD, regime_code, regime_name)

    ! Dimensionless terminal velocity
    Vt_star = Vt * (rho_f**2 / (mu*(rho_p-rho_f)*g))**(1.0d0/3.0d0)

    ! Forces on particle
    F_grav = PI/6.0d0 * dp**3 * rho_p * g
    F_buoy = PI/6.0d0 * dp**3 * rho_f * g
    W_net = F_grav - F_buoy
    F_drag = 0.5d0 * CD * rho_f * Vt**2 * PI/4.0d0 * dp**2

    ! Particle relaxation time
    tau_relax = rho_p * dp**2 / (18.0d0 * mu)

    ! Hindered settling (Richardson-Zaki)
    if (conc > 0.0d0) then
        if (Re_p < 0.2d0) then
            n_RZ = 4.65d0
        else if (Re_p < 1.0d0) then
            n_RZ = 4.35d0 * Re_p**(-0.03d0)
        else if (Re_p < 500.0d0) then
            n_RZ = 4.45d0 * Re_p**(-0.1d0)
        else
            n_RZ = 2.39d0
        end if
        Vt_hindered = Vt * (1.0d0 - conc)**n_RZ
    else
        n_RZ = 0.0d0
        Vt_hindered = Vt
    end if

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   SEDIMENTATION & TERMINAL VELOCITY ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Particle Diameter dp      = ', dp, ' m'
    write(*,'(A,ES12.4,A)') '  Particle Density rho_p    = ', rho_p, ' kg/m3'
    write(*,'(A,ES12.4,A)') '  Fluid Density rho_f       = ', rho_f, ' kg/m3'
    write(*,'(A,ES12.4,A)') '  Fluid Viscosity mu        = ', mu, ' Pa.s'
    write(*,'(A,ES12.4,A)') '  Gravity                   = ', g, ' m/s2'
    write(*,'(A,ES12.4)')   '  Shape Factor phi          = ', phi_shape
    write(*,'(A,ES12.4)')   '  Volume Concentration      = ', conc
    write(*,*)
    write(*,'(A)') '--- DIMENSIONLESS PARAMETERS --------------------------------'
    write(*,'(A,ES12.4)')   '  Archimedes Number Ar      = ', Ar
    write(*,'(A,ES12.4)')   '  Dimensionless dp*         = ', dp_star
    write(*,'(A,ES12.4)')   '  Dimensionless Vt*         = ', Vt_star
    write(*,*)
    write(*,'(A)') '--- TERMINAL VELOCITY RESULTS -------------------------------'
    write(*,'(A,ES12.4,A)') '  Terminal Velocity Vt      = ', Vt, ' m/s'
    write(*,'(A,ES12.4)')   '  Particle Reynolds Re_p    = ', Re_p
    write(*,'(A,ES12.4)')   '  Drag Coefficient CD       = ', CD
    write(*,'(A,A)')        '  Settling Regime           = ', trim(regime_name)
    write(*,*)
    write(*,'(A)') '--- FORCE BALANCE -------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Gravity Force             = ', F_grav, ' N'
    write(*,'(A,ES12.4,A)') '  Buoyancy Force            = ', F_buoy, ' N'
    write(*,'(A,ES12.4,A)') '  Net Weight                = ', W_net, ' N'
    write(*,'(A,ES12.4,A)') '  Drag Force at Vt          = ', F_drag, ' N'
    write(*,'(A,ES12.4,A)') '  Relaxation Time           = ', tau_relax, ' s'
    write(*,*)
    write(*,'(A)') '--- HINDERED SETTLING ---------------------------------------'
    write(*,'(A,ES12.4)')   '  Richardson-Zaki n         = ', n_RZ
    write(*,'(A,ES12.4,A)') '  Hindered Velocity         = ', Vt_hindered, ' m/s'
    write(*,'(A,ES12.4)')   '  Ratio Vt_hind / Vt        = ', Vt_hindered/max(Vt,1.0d-30)
    write(*,*)

    ! Vt vs dp sweep
    write(*,'(A)') '--- VT VS PARTICLE SIZE SWEEP -------------------------------'
    write(*,'(A)') '  dp[m]         Vt[m/s]       Re_p          CD            regime'
    write(*,'(A)') '  -----------------------------------------------------------------------'
    do i = 1, 60
        dp_i = dp * 0.01d0 * (1000.0d0)**(dble(i-1)/59.0d0)
        call compute_terminal_velocity(dp_i, rho_p, rho_f, mu, g, phi_shape, &
                                       Vt_i, Re_i, CD_i, regime_code, regime_name)
        write(*,'(ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,A)') &
            dp_i, Vt_i, Re_i, CD_i, trim(regime_name)
    end do
    write(*,*)

    ! CD vs Re sweep
    write(*,'(A)') '--- CD VS RE SWEEP ------------------------------------------'
    write(*,'(A)') '  Re_p          CD_Stokes     CD_inter      CD_Newton'
    write(*,'(A)') '  -----------------------------------------------------------'
    do i = 1, 60
        Re_i = 0.01d0 * (1.0d6/0.01d0)**(dble(i-1)/59.0d0)
        write(*,'(ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') &
            Re_i, 24.0d0/Re_i, &
            24.0d0/Re_i*(1.0d0+0.15d0*Re_i**0.687d0), &
            0.44d0
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Stokes (Re<0.1): CD = 24/Re; Vt = (rho_p-rho_f)g dp^2/(18 mu).'
    write(*,'(A)') '  Intermediate: CD = 24/Re (1 + 0.15 Re^0.687) Schiller-Naumann.'
    write(*,'(A)') '  Newton (Re>1000): CD = 0.44.'
    write(*,'(A)') '  Shape correction: CD = CD_sphere / phi.'
    write(*,'(A)') '  Richardson-Zaki hindered: Vt_h = Vt (1-c)^n.'

contains

    subroutine compute_terminal_velocity(d, rhop, rhof, muf, grav, phi, &
                                         Vout, Reout, CDout, rcode, rname)
        implicit none
        double precision, intent(in) :: d, rhop, rhof, muf, grav, phi
        double precision, intent(out) :: Vout, Reout, CDout
        integer, intent(out) :: rcode
        character(len=60), intent(out) :: rname
        double precision :: Vt_old, Vt_new, Re_loc, CD_loc, Ap
        integer :: iter

        Ap = PI/4.0d0 * d**2

        ! Initial guess: Stokes
        Vt_new = (rhop - rhof) * grav * d**2 / (18.0d0 * muf)
        if (Vt_new < 0.0d0) Vt_new = abs(Vt_new)

        do iter = 1, 500
            Vt_old = Vt_new
            Re_loc = rhof * Vt_old * d / muf
            if (Re_loc < 1.0d-10) Re_loc = 1.0d-10

            if (Re_loc < 0.1d0) then
                CD_loc = 24.0d0 / Re_loc
            else if (Re_loc < 1000.0d0) then
                CD_loc = 24.0d0/Re_loc * (1.0d0 + 0.15d0*Re_loc**0.687d0)
            else if (Re_loc < 2.0d5) then
                CD_loc = 0.44d0
            else
                CD_loc = 0.1d0   ! drag crisis for spheres
            end if

            ! Shape correction
            CD_loc = CD_loc / phi

            ! New velocity from force balance: W = Drag
            ! (rho_p-rho_f)*g*Vol = 0.5*CD*rho_f*V^2*Ap
            Vt_new = sqrt(4.0d0*d*(rhop-rhof)*grav / (3.0d0*CD_loc*rhof))

            if (abs(Vt_new - Vt_old) < 1.0d-12*max(Vt_new,1.0d-30)) exit
        end do

        Vout = Vt_new
        Reout = rhof * Vt_new * d / muf
        CDout = CD_loc

        if (Reout < 0.1d0) then
            rcode = 1; rname = 'Stokes (creeping flow)'
        else if (Reout < 1000.0d0) then
            rcode = 2; rname = 'Intermediate (Schiller-Naumann)'
        else if (Reout < 2.0d5) then
            rcode = 3; rname = 'Newton (turbulent wake)'
        else
            rcode = 4; rname = 'Supercritical (drag crisis)'
        end if
    end subroutine compute_terminal_velocity

end program sedimentation

Solver Description

Compute terminal settling velocity for single or hindered particles through Stokes, intermediate, and Newton drag regimes. Includes shape factors.

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

Execution Command:

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

sedimentation < input.txt

πŸ“₯ Downloads & Local Files

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

! Particle diameter dp\nrp_init\nrf_init\nViscosity μ [Pa·s]\nGravity g [m/s²]\nShape factor Ï‒ (1 = sphere)\nVolume concentration c (0Γ’β‚¬β€œ0.65)
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
πŸ“¦ Download Offline Bundle (ZIP) 10 dl πŸ“„ Download input.txt 10 dl πŸ“₯ Download sedimentation.f90 7 dl