💻 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.
Bernoulli & Extended Energy Equation
Core Numerical Engine in Fortran 90 • 31 total downloads
bernoulli_energy.f90
! =========================================================================
! Source File: bernoulli_energy.f90
! =========================================================================
! ============================================================================
! ThermoFluidCalc — Bernoulli & Extended Energy Equation Solver
! ============================================================================
program bernoulli_energy
implicit none
! Input variables
integer :: mode ! 1 = Simple Bernoulli, 2 = Extended Energy
integer :: solve_for ! 1=P2, 2=V2, 3=z2, 4=Q, 5=h_p, 6=h_L
double precision :: rho ! Fluid density [kg/m³]
double precision :: P1, V1, z1 ! Point 1: Pressure [Pa], Velocity [m/s], Elevation [m]
double precision :: P2, V2, z2 ! Point 2: Pressure [Pa], Velocity [m/s], Elevation [m]
double precision :: D1, D2 ! Pipe diameters [m]
double precision :: Q ! Volumetric flow rate [m³/s]
double precision :: h_p, h_t, h_L ! Pump head [m], Turbine head [m], Head loss [m]
! Constants
double precision, parameter :: g = 9.81d0
double precision, parameter :: pi = 3.141592653589793d0
! Intermediate and output variables
double precision :: hp1, hv1, hz1, H1
double precision :: hp2, hv2, hz2, H2
double precision :: C1, C2, RHS, Q2, V2_sq
double precision :: W_p, W_t ! Power [W]
character(len=50) :: solved_var_name
! Read input from stdin
read(*,*) mode
read(*,*) solve_for
read(*,*) rho
read(*,*) P1
read(*,*) V1
read(*,*) z1
read(*,*) P2
read(*,*) V2
read(*,*) z2
read(*,*) D1
read(*,*) D2
read(*,*) Q
read(*,*) h_p
read(*,*) h_t
read(*,*) h_L
! Override simple Bernoulli head parameters to zero if mode=1
if (mode == 1) then
h_p = 0.0d0
h_t = 0.0d0
h_L = 0.0d0
end if
! Calculate velocity from flow rate and diameters if available
! Point 1
if (Q > 0.0d0 .and. D1 > 0.0d0) then
V1 = 4.0d0 * Q / (pi * D1**2)
else if (V1 > 0.0d0 .and. D1 > 0.0d0 .and. Q <= 0.0d0) then
Q = V1 * pi * D1**2 / 4.0d0
end if
! Point 2 (only if we are not solving for V2 or Q)
if (solve_for /= 2 .and. solve_for /= 4) then
if (Q > 0.0d0 .and. D2 > 0.0d0) then
V2 = 4.0d0 * Q / (pi * D2**2)
else if (V2 > 0.0d0 .and. D2 > 0.0d0 .and. Q <= 0.0d0) then
Q = V2 * pi * D2**2 / 4.0d0
end if
end if
! Perform solving based on target variable
select case (solve_for)
case (1) ! Solve for P2
solved_var_name = "Downstream Pressure (P2)"
P2 = P1 + rho * g * ((V1**2 - V2**2) / (2.0d0 * g) + (z1 - z2) + h_p - h_t - h_L)
case (2) ! Solve for V2
solved_var_name = "Downstream Velocity (V2)"
V2_sq = V1**2 + 2.0d0 * g * ((P1 - P2) / (rho * g) + (z1 - z2) + h_p - h_t - h_L)
if (V2_sq < 0.0d0) then
write(*,*) "ERROR: No physical solution for V2 (energy state is impossible)"
stop
end if
V2 = sqrt(V2_sq)
if (D2 > 0.0d0) then
Q = V2 * pi * D2**2 / 4.0d0
end if
case (3) ! Solve for z2
solved_var_name = "Downstream Elevation (z2)"
z2 = z1 + (P1 - P2) / (rho * g) + (V1**2 - V2**2) / (2.0d0 * g) + h_p - h_t - h_L
case (4) ! Solve for Q
solved_var_name = "System Flow Rate (Q)"
! Coefficients for Q^2 relation
C1 = 0.0d0
C2 = 0.0d0
if (D1 > 0.0d0) C1 = 8.0d0 / (g * pi**2 * D1**4)
if (D2 > 0.0d0) C2 = 8.0d0 / (g * pi**2 * D2**4)
if (D1 > 0.0d0 .and. D2 > 0.0d0) then
if (abs(D1 - D2) < 1.0d-7) then
write(*,*) "ERROR: Cannot solve for Q when diameters are equal (V1 = V2)"
stop
end if
RHS = (P2 - P1) / (rho * g) + (z2 - z1) + h_t + h_L - h_p
Q2 = RHS / (C1 - C2)
if (Q2 < 0.0d0) then
write(*,*) "ERROR: No physical solution for Q (check energy terms direction)"
stop
end if
Q = sqrt(Q2)
V1 = 4.0d0 * Q / (pi * D1**2)
V2 = 4.0d0 * Q / (pi * D2**2)
else if (D1 > 0.0d0 .and. D2 <= 0.0d0) then
! V2 is fixed, solve for Q
RHS = (P2 - P1) / (rho * g) + V2**2 / (2.0d0 * g) + (z2 - z1) + h_t + h_L - h_p
Q2 = RHS / C1
if (Q2 < 0.0d0) then
write(*,*) "ERROR: No physical solution for Q"
stop
end if
Q = sqrt(Q2)
V1 = 4.0d0 * Q / (pi * D1**2)
else if (D1 <= 0.0d0 .and. D2 > 0.0d0) then
! V1 is fixed, solve for Q
RHS = (P1 - P2) / (rho * g) + V1**2 / (2.0d0 * g) + (z1 - z2) + h_p - h_t - h_L
Q2 = RHS / C2
if (Q2 < 0.0d0) then
write(*,*) "ERROR: No physical solution for Q"
stop
end if
Q = sqrt(Q2)
V2 = 4.0d0 * Q / (pi * D2**2)
else
write(*,*) "ERROR: Cannot solve for Q without at least one pipe diameter"
stop
end if
case (5) ! Solve for h_p
solved_var_name = "Required Pump Head (h_p)"
h_p = (P2 - P1) / (rho * g) + (V2**2 - V1**2) / (2.0d0 * g) + (z2 - z1) + h_t + h_L
case (6) ! Solve for h_L
solved_var_name = "Head Loss (h_L)"
h_L = (P1 - P2) / (rho * g) + (V1**2 - V2**2) / (2.0d0 * g) + (z1 - z2) + h_p - h_t
if (h_L < 0.0d0) then
write(*,*) "ERROR: Calculated Head Loss is negative (violates thermodynamics, check inputs)"
stop
end if
case default
write(*,*) "ERROR: Invalid solve target"
stop
end select
! Calculate final component heads
hp1 = P1 / (rho * g)
hv1 = V1**2 / (2.0d0 * g)
hz1 = z1
H1 = hp1 + hv1 + hz1
hp2 = P2 / (rho * g)
hv2 = V2**2 / (2.0d0 * g)
hz2 = z2
H2 = hp2 + hv2 + hz2
! Calculate power values
W_p = 0.0d0
W_t = 0.0d0
if (Q > 0.0d0) then
if (h_p > 0.0d0) W_p = rho * g * Q * h_p
if (h_t > 0.0d0) W_t = rho * g * Q * h_t
end if
! Output details
write(*, '(A, A)') "Solved Variable = ", trim(solved_var_name)
write(*, '(A, F14.4)') "Fluid Density (rho) = ", rho
write(*, '(A, F14.4)') "Upstream Pressure (P1) = ", P1 / 1000.0d0 ! Pa to kPa
write(*, '(A, F14.4)') "Upstream Velocity (V1) = ", V1
write(*, '(A, F14.4)') "Upstream Elevation (z1) = ", z1
write(*, '(A, F14.4)') "Downstream Pressure (P2) = ", P2 / 1000.0d0 ! Pa to kPa
write(*, '(A, F14.4)') "Downstream Velocity (V2) = ", V2
write(*, '(A, F14.4)') "Downstream Elevation (z2) = ", z2
write(*, '(A, F14.4)') "Flow Rate (Q_lmin) = ", Q * 1000.0d0 * 60.0d0 ! m3/s to L/min
write(*, '(A, F14.4)') "Flow Rate (Q_m3h) = ", Q * 3600.0d0 ! m3/s to m3/h
write(*, '(A, F14.4)') "Pipe 1 Diameter (D1) = ", D1 * 1000.0d0 ! m to mm
write(*, '(A, F14.4)') "Pipe 2 Diameter (D2) = ", D2 * 1000.0d0 ! m to mm
write(*, '(A, F14.4)') "Pump Head Added (h_p) = ", h_p
write(*, '(A, F14.4)') "Turbine Head Extracted (h_t) = ", h_t
write(*, '(A, F14.4)') "Head Loss (h_L) = ", h_L
write(*, '(A, F14.4)') "Pressure Head 1 (hp1) = ", hp1
write(*, '(A, F14.4)') "Velocity Head 1 (hv1) = ", hv1
write(*, '(A, F14.4)') "Elevation Head 1 (hz1) = ", hz1
write(*, '(A, F14.4)') "Total Head 1 (H1) = ", H1
write(*, '(A, F14.4)') "Pressure Head 2 (hp2) = ", hp2
write(*, '(A, F14.4)') "Velocity Head 2 (hv2) = ", hv2
write(*, '(A, F14.4)') "Elevation Head 2 (hz2) = ", hz2
write(*, '(A, F14.4)') "Total Head 2 (H2) = ", H2
write(*, '(A, F14.4)') "Pump Power (W_p) = ", W_p / 1000.0d0 ! W to kW
write(*, '(A, F14.4)') "Pump Power (HP) = ", W_p / 745.7d0 ! W to HP
write(*, '(A, F14.4)') "Turbine Power (W_t) = ", W_t / 1000.0d0 ! W to kW
write(*, '(A, F14.4)') "Turbine Power (HP) = ", W_t / 745.7d0 ! W to HP
end program bernoulli_energy
Solver Description
Calculate pressure, velocity, elevation, flow rate, pump head, or losses using Bernoulli & Extended Energy Equations. Interactive EGL/HGL diagram plotting.
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 bernoulli_energy.f90 -o bernoulli_energy_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
bernoulli_energy_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Solver Mode (1=Simple Bernoulli, 2=Extended Energy)
2
! Solve For Target (1=P2, 2=V2, 3=z2, 4=Q, 5=hp, 6=hL)
1
! Fluid Density [kg/m3]
1000.0
! Upstream Pressure P1 [Pa]
200000.0
! Upstream Velocity V1 [m/s]
2.0
! Upstream Elevation z1 [m]
10.0
! Downstream Pressure P2 [Pa]
100000.0
! Downstream Velocity V2 [m/s]
0.0
! Downstream Elevation z2 [m]
5.0
! Pipe 1 Diameter D1 [m]
0.1
! Pipe 2 Diameter D2 [m]
0.08
! Flow Rate Q [m3/s]
0.0
! Pump Head added hp [m]
15.0
! Turbine Head extracted ht [m]
0.0
! Head Loss hL [m]
2.0
2
! Solve For Target (1=P2, 2=V2, 3=z2, 4=Q, 5=hp, 6=hL)
1
! Fluid Density [kg/m3]
1000.0
! Upstream Pressure P1 [Pa]
200000.0
! Upstream Velocity V1 [m/s]
2.0
! Upstream Elevation z1 [m]
10.0
! Downstream Pressure P2 [Pa]
100000.0
! Downstream Velocity V2 [m/s]
0.0
! Downstream Elevation z2 [m]
5.0
! Pipe 1 Diameter D1 [m]
0.1
! Pipe 2 Diameter D2 [m]
0.08
! Flow Rate Q [m3/s]
0.0
! Pump Head added hp [m]
15.0
! Turbine Head extracted ht [m]
0.0
! Head Loss hL [m]
2.0