💻 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.
TEMA Shell-and-Tube HX Sizing
Core Numerical Engine in Fortran 90 • 50 total downloads
shell_tube_tema.f90
! =========================================================================
! Source File: shell_tube_tema.f90
! =========================================================================
program shell_tube_tema
implicit none
! Inputs
integer :: N_shell, N_p, layout, hot_side
double precision :: T_hot_in, T_hot_out, t_cold_in, t_cold_out
double precision :: m_s, m_t
double precision :: rho_s, Cp_s, mu_s, k_s
double precision :: rho_t, Cp_t, mu_t, k_t
double precision :: do_val, tw, L, Pt, B
double precision :: Rfs, Rft, kw
! Calculated properties
double precision :: di, K1, n1, Ds, C, As, Gs, Res, Prs, Nus, hs, fs
double precision :: Aat, vt, Ret, Prt, Nut, ht, ft
double precision :: Rw, U_clean, U_fouled, Q_s, Q_t, Q
double precision :: dT1, dT2, LMTD_cf, F, P, R, K, P_1
double precision :: denom_val, val_num, val_den, num_val
double precision :: A_req, A_prov, dP_shell, dP_tube, De
integer :: Nt, Nb, iostat_val, i
logical :: has_cross_error
! Shell and tube temperatures mapping
double precision :: T_hi, T_ho, T_ci, T_co
double precision :: T_s_in, T_s_out, T_t_in, T_t_out
double precision, parameter :: pi = 3.141592653589793d0
double precision :: U_clean_inv, U_fouled_inv
! Read all input variables
read(*,*,iostat=iostat_val) N_shell
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid N_shell.'
stop
end if
read(*,*,iostat=iostat_val) N_p
read(*,*,iostat=iostat_val) layout
read(*,*,iostat=iostat_val) T_hot_in
read(*,*,iostat=iostat_val) T_hot_out
read(*,*,iostat=iostat_val) t_cold_in
read(*,*,iostat=iostat_val) t_cold_out
read(*,*,iostat=iostat_val) m_s
read(*,*,iostat=iostat_val) m_t
read(*,*,iostat=iostat_val) rho_s
read(*,*,iostat=iostat_val) Cp_s
read(*,*,iostat=iostat_val) mu_s
read(*,*,iostat=iostat_val) k_s
read(*,*,iostat=iostat_val) rho_t
read(*,*,iostat=iostat_val) Cp_t
read(*,*,iostat=iostat_val) mu_t
read(*,*,iostat=iostat_val) k_t
read(*,*,iostat=iostat_val) do_val
read(*,*,iostat=iostat_val) tw
read(*,*,iostat=iostat_val) L
read(*,*,iostat=iostat_val) Pt
read(*,*,iostat=iostat_val) B
read(*,*,iostat=iostat_val) Rfs
read(*,*,iostat=iostat_val) Rft
read(*,*,iostat=iostat_val) kw
read(*,*,iostat=iostat_val) hot_side
! Basic sanity checks
if (N_shell < 1) N_shell = 1
if (N_p < 1) N_p = 2
if (layout /= 1 .and. layout /= 2) layout = 1
if (hot_side /= 1 .and. hot_side /= 2) hot_side = 1
! Map process temperatures to hot/cold streams
T_hi = T_hot_in
T_ho = T_hot_out
T_ci = t_cold_in
T_co = t_cold_out
! Map to shell and tube sides
if (hot_side == 1) then
T_s_in = T_hi
T_s_out = T_ho
T_t_in = T_ci
T_t_out = T_co
else
T_s_in = T_ci
T_s_out = T_co
T_t_in = T_hi
T_t_out = T_ho
end if
! Compute heat duty (average of shell and tube sides)
Q_s = m_s * Cp_s * abs(T_s_in - T_s_out)
Q_t = m_t * Cp_t * abs(T_t_in - T_t_out)
Q = (Q_s + Q_t) / 2.0d0
! Counter-flow LMTD
dT1 = T_hi - T_co
dT2 = T_ho - T_ci
has_cross_error = .false.
if (dT1 <= 0.0d0 .or. dT2 <= 0.0d0) then
LMTD_cf = 0.0d0
has_cross_error = .true.
else if (abs(dT1 - dT2) < 1.0d-6) then
LMTD_cf = dT1
else
LMTD_cf = (dT1 - dT2) / log(dT1 / dT2)
end if
! P and R calculations
if (T_hi - T_ci <= 0.0d0) then
has_cross_error = .true.
else
P = (T_co - T_ci) / (T_hi - T_ci)
if (abs(T_co - T_ci) < 1.0d-8) then
R = 1.0d0
else
R = (T_hi - T_ho) / (T_co - T_ci)
end if
end if
F = -1.0d0
if (.not. has_cross_error) then
if (P < 0.0d0 .or. P >= 1.0d0 .or. R < 0.0d0) then
has_cross_error = .true.
else
if (abs(R - 1.0d0) < 1.0d-5) then
R = 1.0001d0
end if
if (1.0d0 - P * R <= 0.0d0 .or. 1.0d0 - P <= 0.0d0) then
has_cross_error = .true.
else
K = ((1.0d0 - P * R) / (1.0d0 - P))**(1.0d0 / dble(N_shell))
if (abs(K - R) < 1.0d-8) then
P_1 = P / dble(N_shell)
else
P_1 = (1.0d0 - K) / (R - K)
end if
denom_val = 2.0d0 - P_1 * (R + 1.0d0 + sqrt(R**2 + 1.0d0))
num_val = 2.0d0 - P_1 * (R + 1.0d0 - sqrt(R**2 + 1.0d0))
if (denom_val <= 0.0d0 .or. num_val <= 0.0d0 .or. (1.0d0 - P_1) <= 0.0d0 .or. (1.0d0 - R*P_1) <= 0.0d0) then
has_cross_error = .true.
else
val_num = sqrt(R**2 + 1.0d0) * log((1.0d0 - P_1) / (1.0d0 - R * P_1))
val_den = (R - 1.0d0) * log(num_val / denom_val)
if (val_den /= 0.0d0) then
F = val_num / val_den
if (F <= 0.0d0 .or. F > 1.0d0) then
has_cross_error = .true.
end if
else
has_cross_error = .true.
end if
end if
end if
end if
end if
! Iterative design logic
if (.not. has_cross_error) then
! Lookup K1 and n1
call get_bundle_constants(layout, N_p, K1, n1)
di = do_val - 2.0d0 * tw
! De: Equivalent diameter
if (layout == 1) then
De = (2.0d0 * sqrt(3.0d0) * Pt**2 - pi * do_val**2) / (pi * do_val)
else
De = 4.0d0 * (Pt**2 - pi * do_val**2 / 4.0d0) / (pi * do_val)
end if
! Start sizing loop
Nt = N_p
do while (Nt <= 100000)
! Shell diameter
Ds = do_val * (dble(Nt) / K1)**(1.0d0 / n1)
! Shell cross-flow area
C = Pt - do_val
As = (Ds * C * B) / Pt
! Shell-side velocity, Re, Pr
Gs = m_s / As
Res = (Gs * De) / mu_s
Prs = (mu_s * Cp_s) / k_s
! Shell film coefficient
Nus = 0.36d0 * (Res**0.55d0) * (Prs**(1.0d0/3.0d0))
hs = (Nus * k_s) / De
! Tube-side area, velocity, Re, Pr
Aat = (dble(Nt) * pi * di**2) / (4.0d0 * dble(N_p))
vt = m_t / (rho_t * Aat)
Ret = (rho_t * vt * di) / mu_t
Prt = (mu_t * Cp_t) / k_t
! Tube film coefficient
if (Ret >= 2100.0d0) then
Nut = 0.027d0 * (Ret**0.8d0) * (Prt**(1.0d0/3.0d0))
else
Nut = 1.86d0 * (Ret * Prt * di / L)**(1.0d0/3.0d0)
if (Nut < 3.66d0) Nut = 3.66d0
end if
ht = (Nut * k_t) / di
! Wall conduction resistance
Rw = (do_val * log(do_val / di)) / (2.0d0 * kw)
! Overall heat transfer coefficients (based on Ao)
U_clean_inv = 1.0d0/hs + Rw + do_val / (di * ht)
U_clean = 1.0d0 / U_clean_inv
U_fouled_inv = 1.0d0/hs + Rfs + Rw + (do_val/di)*Rft + do_val / (di * ht)
U_fouled = 1.0d0 / U_fouled_inv
! Required Area based on fouled U
A_req = Q / (U_fouled * F * LMTD_cf)
! Provided Area
A_prov = dble(Nt) * pi * do_val * L
if (A_prov >= A_req) then
exit
end if
Nt = Nt + N_p
end do
! Final outputs computations
if (Ret >= 2100.0d0) then
ft = 0.184d0 * (Ret**(-0.2d0))
else
ft = 64.0d0 / Ret
end if
dP_tube = (ft * (L * dble(N_p) / di) + 4.0d0 * dble(N_p)) * (rho_t * vt**2) / 2.0d0
fs = 1.44d0 * (Res**(-0.3d0))
Nb = idint(L / B) - 1
if (Nb < 0) Nb = 0
dP_shell = (fs * Gs**2 * Ds * dble(Nb + 1)) / (2.0d0 * rho_s * De)
else
Nt = 0
Ds = 0.0d0
De = 0.0d0
hs = 0.0d0
ht = 0.0d0
U_clean = 0.0d0
U_fouled = 0.0d0
dP_shell = 0.0d0
dP_tube = 0.0d0
A_req = 0.0d0
A_prov = 0.0d0
Nb = 0
Res = 0.0d0
Ret = 0.0d0
Gs = 0.0d0
vt = 0.0d0
end if
! Print results to stdout in KEY=value format for easy parsing
write(*,'(A,I2)') 'N_SHELL=', N_shell
write(*,'(A,I2)') 'N_PASSES=', N_p
write(*,'(A,I2)') 'LAYOUT=', layout
write(*,'(A,F14.4)') 'T_HI=', T_hi
write(*,'(A,F14.4)') 'T_HO=', T_ho
write(*,'(A,F14.4)') 'T_CI=', T_ci
write(*,'(A,F14.4)') 'T_CO=', T_co
write(*,'(A,F14.4)') 'M_S=', m_s
write(*,'(A,F14.4)') 'M_T=', m_t
write(*,'(A,F14.2)') 'Q=', Q
write(*,'(A,I6)') 'N_T=', Nt
write(*,'(A,F14.6)') 'D_S=', Ds
write(*,'(A,F14.4)') 'H_S=', hs
write(*,'(A,F14.4)') 'H_T=', ht
write(*,'(A,F14.4)') 'U_CLEAN=', U_clean
write(*,'(A,F14.4)') 'U_FOULED=', U_fouled
write(*,'(A,F14.4)') 'DP_SHELL=', dP_shell
write(*,'(A,F14.4)') 'DP_TUBE=', dP_tube
write(*,'(A,F14.4)') 'A_REQ=', A_req
write(*,'(A,F14.4)') 'A_PROV=', A_prov
write(*,'(A,F14.4)') 'LMTD_CF=', LMTD_cf
write(*,'(A,F14.6)') 'F_FACTOR=', F
write(*,'(A,F14.4)') 'RE_S=', Res
write(*,'(A,F14.4)') 'RE_T=', Ret
write(*,'(A,F14.4)') 'V_S=', Gs / rho_s
write(*,'(A,F14.4)') 'V_T=', vt
write(*,'(A,I4)') 'N_BAFFLES=', Nb
if (has_cross_error) then
write(*,'(A,I1)') 'HAS_CROSS_ERROR=', 1
else
write(*,'(A,I1)') 'HAS_CROSS_ERROR=', 0
end if
! Print temperature profiles for visualization
write(*,'(A)') '--- TEMPERATURE PROFILE ALONG EXCHANGER --------------------'
write(*,'(A)') ' Position % T_hot [C] T_cold [C]'
write(*,'(A)') ' --------------------------------------------------'
do i = 0, 20
write(*,'(F10.1,4X,F10.2,4X,F10.2)') dble(i)/20.0d0*100.0d0, &
T_hi - (T_hi - T_ho) * dble(i)/20.0d0, &
T_co - (T_co - T_ci) * dble(i)/20.0d0
end do
contains
subroutine get_bundle_constants(lay, Np, K1, n1)
integer, intent(in) :: lay
integer, intent(in) :: Np
double precision, intent(out) :: K1, n1
if (lay == 1) then ! Triangular 30
select case (Np)
case (1)
K1 = 0.319d0; n1 = 2.142d0
case (2)
K1 = 0.249d0; n1 = 2.207d0
case (4)
K1 = 0.175d0; n1 = 2.285d0
case (6)
K1 = 0.0743d0; n1 = 2.499d0
case default ! 8 or more
K1 = 0.0365d0; n1 = 2.675d0
end select
else ! Square 90
select case (Np)
case (1)
K1 = 0.215d0; n1 = 2.207d0
case (2)
K1 = 0.156d0; n1 = 2.291d0
case (4)
K1 = 0.158d0; n1 = 2.311d0
case (6)
K1 = 0.081d0; n1 = 2.473d0
case default ! 8 or more
K1 = 0.0321d0; n1 = 2.721d0
end select
end if
end subroutine get_bundle_constants
end program shell_tube_tema
Solver Description
Perform full TEMA-level shell-and-tube sizing. This iteratively solves the required tube count $N_t$, shell inside diameter $D_s$, local heat transfer coefficients on both shell and tube sides, and pressure drops.
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 shell_tube_tema.f90 -o shell_tube_tema
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
shell_tube_tema < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Shell Type (1=E, 2=F, etc.)
1
! Tube Layout (1=Triangular, 2=Square)
2
! Tube Pass Count (1, 2, 4, 6)
1
! Hot fluid inlet temp Thi [°C]
90.0
! Hot fluid outlet temp Tho [°C]
40.0
! Cold fluid inlet temp Tci [°C]
20.0
! Cold fluid outlet temp Tco [°C]
35.0
! Hot fluid mass flow rate mh [kg/s]
2.0
! Cold fluid mass flow rate mc [kg/s]
3.0
! Hot fluid density [kg/m3]
990.0
! Hot fluid specific heat [J/kg-K]
4180.0
! Hot fluid dynamic viscosity [Pa-s]
0.0008
! Hot fluid thermal conductivity [W/m-K]
0.6
! Cold fluid density [kg/m3]
998.0
! Cold fluid specific heat [J/kg-K]
4180.0
! Cold fluid dynamic viscosity [Pa-s]
0.001
! Cold fluid thermal conductivity [W/m-K]
0.6
! Tube outer diameter [m]
0.01905
! Tube wall thickness [m]
0.00165
! Active tube length [m]
3.0
! Tube pitch [m]
0.02381
! Baffle spacing [m]
0.15
! Hot fluid fouling factor [m2-K/W]
0.0002
! Cold fluid fouling factor [m2-K/W]
0.0002
! Tube wall thermal conductivity [W/m-K]
50.0
! Shell-side fluid selection (1=Hot fluid in shell, 2=Cold fluid in shell)
1
1
! Tube Layout (1=Triangular, 2=Square)
2
! Tube Pass Count (1, 2, 4, 6)
1
! Hot fluid inlet temp Thi [°C]
90.0
! Hot fluid outlet temp Tho [°C]
40.0
! Cold fluid inlet temp Tci [°C]
20.0
! Cold fluid outlet temp Tco [°C]
35.0
! Hot fluid mass flow rate mh [kg/s]
2.0
! Cold fluid mass flow rate mc [kg/s]
3.0
! Hot fluid density [kg/m3]
990.0
! Hot fluid specific heat [J/kg-K]
4180.0
! Hot fluid dynamic viscosity [Pa-s]
0.0008
! Hot fluid thermal conductivity [W/m-K]
0.6
! Cold fluid density [kg/m3]
998.0
! Cold fluid specific heat [J/kg-K]
4180.0
! Cold fluid dynamic viscosity [Pa-s]
0.001
! Cold fluid thermal conductivity [W/m-K]
0.6
! Tube outer diameter [m]
0.01905
! Tube wall thickness [m]
0.00165
! Active tube length [m]
3.0
! Tube pitch [m]
0.02381
! Baffle spacing [m]
0.15
! Hot fluid fouling factor [m2-K/W]
0.0002
! Cold fluid fouling factor [m2-K/W]
0.0002
! Tube wall thermal conductivity [W/m-K]
50.0
! Shell-side fluid selection (1=Hot fluid in shell, 2=Cold fluid in shell)
1