🔬
Solver Purpose & Physical Scope
Chevron Plate Heat Exchanger (PHE) rating and sizing. Evaluates enlargement factor $\phi$, hydraulic diameter $D_h$, convection coefficients using Shah & Martin correlations, and port/channel pressure drops.
📂 Discipline: Exchangers
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 540 times
📄 Source File:
hx_plate.f90
📁 calcul/HeatExchangers /
hx_plate.f90
program hx_plate
implicit none
! Inputs
integer :: mode, Nplates
double precision :: T_hi, T_ho, T_ci, T_co
double precision :: m_h, m_c
double precision :: rho_h, Cp_h, mu_h, k_h, mu_wh
double precision :: rho_c, Cp_c, mu_c, k_c, mu_wc
double precision :: W, Lp, t_plate, b_depth, pc_pitch, beta_deg, d_port_mm
double precision :: dP_max_kpa
double precision :: Rfh, Rfc, kw
! Calculated properties
double precision :: X, phi, Dh, Q_h, Q_c, Q, dT1, dT2, LMTD_cf
double precision :: hs_h, hs_c, dP_h, dP_c, dP_port_h, dP_port_c, U_f, U_c, A_prov, A_req
double precision :: Re_h, Re_c, vch_h, vch_c
integer :: Nplates_curr, iostat_val
logical :: has_cross_error
double precision, parameter :: pi = 3.141592653589793d0
double precision :: dP_max
! Read all inputs from stdin
read(*,*,iostat=iostat_val) mode
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid inputs'
stop
end if
read(*,*,iostat=iostat_val) Nplates
read(*,*,iostat=iostat_val) T_hi
read(*,*,iostat=iostat_val) T_ho
read(*,*,iostat=iostat_val) T_ci
read(*,*,iostat=iostat_val) T_co
read(*,*,iostat=iostat_val) m_h
read(*,*,iostat=iostat_val) m_c
read(*,*,iostat=iostat_val) rho_h
read(*,*,iostat=iostat_val) Cp_h
read(*,*,iostat=iostat_val) mu_h
read(*,*,iostat=iostat_val) k_h
read(*,*,iostat=iostat_val) mu_wh
read(*,*,iostat=iostat_val) rho_c
read(*,*,iostat=iostat_val) Cp_c
read(*,*,iostat=iostat_val) mu_c
read(*,*,iostat=iostat_val) k_c
read(*,*,iostat=iostat_val) mu_wc
read(*,*,iostat=iostat_val) W
read(*,*,iostat=iostat_val) Lp
read(*,*,iostat=iostat_val) t_plate
read(*,*,iostat=iostat_val) b_depth
read(*,*,iostat=iostat_val) pc_pitch
read(*,*,iostat=iostat_val) beta_deg
read(*,*,iostat=iostat_val) d_port_mm
read(*,*,iostat=iostat_val) dP_max_kpa
read(*,*,iostat=iostat_val) Rfh
read(*,*,iostat=iostat_val) Rfc
read(*,*,iostat=iostat_val) kw
! Defaults handling
if (mu_wh <= 0.0d0) mu_wh = mu_h
if (mu_wc <= 0.0d0) mu_wc = mu_c
! Basic sanity checks
if (W <= 0.0d0 .or. Lp <= 0.0d0 .or. b_depth <= 0.0d0 .or. pc_pitch <= 0.0d0 .or. beta_deg <= 0.0d0) then
write(*,*) 'ERROR: Geometrical plate values must be positive and non-zero.'
stop
end if
! Compute heat duty (average of hot and cold sides)
Q_h = m_h * Cp_h * abs(T_hi - T_ho)
Q_c = m_c * Cp_c * abs(T_co - T_ci)
Q = (Q_h + Q_c) / 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
if (has_cross_error) then
write(*,*) 'ERROR: Temperature cross error detected. LMTD is undefined.'
stop
end if
! Convert dP_max from kPa to Pa
dP_max = dP_max_kpa * 1000.0d0
! Initial geometrical constants
X = (pi * b_depth) / pc_pitch
phi = (1.0d0 / 6.0d0) * (1.0d0 + sqrt(1.0d0 + X**2) + 4.0d0 * sqrt(1.0d0 + X**2 / 2.0d0))
Dh = 2.0d0 * b_depth / phi
if (mode == 1) then
! Verification mode: evaluate at input plate count
call calculate_phe(Nplates, hs_h, hs_c, dP_h, dP_c, &
dP_port_h, dP_port_c, U_f, U_c, A_prov, A_req, &
Re_h, Re_c, vch_h, vch_c)
else
! Sizing mode: iterate plate count to meet area and pressure limits
Nplates_curr = Nplates
if (mod(Nplates_curr, 2) == 0) Nplates_curr = Nplates_curr + 1
if (Nplates_curr < 5) Nplates_curr = 5
do while (Nplates_curr <= 50000)
call calculate_phe(Nplates_curr, hs_h, hs_c, dP_h, dP_c, &
dP_port_h, dP_port_c, U_f, U_c, A_prov, A_req, &
Re_h, Re_c, vch_h, vch_c)
if (A_prov >= A_req .and. dP_h <= dP_max .and. dP_c <= dP_max) then
Nplates = Nplates_curr
exit
end if
Nplates_curr = Nplates_curr + 2
end do
if (Nplates_curr > 50000) then
write(*,*) 'ERROR: Plate sizing loop did not converge within 50,000 plates.'
stop
end if
end if
! Print outputs to stdout
write(*,'(A,I5)') 'N_PLATES=', Nplates
write(*,'(A,F18.4)') 'ENLARGEMENT_PHI=', phi
write(*,'(A,F18.6)') 'DH=', Dh
write(*,'(A,F18.4)') 'Q=', Q
write(*,'(A,F18.4)') 'LMTD=', LMTD_cf
write(*,'(A,F18.4)') 'A_PROV=', A_prov
write(*,'(A,F18.4)') 'A_REQ=', A_req
write(*,'(A,F18.4)') 'U_CLEAN=', U_c
write(*,'(A,F18.4)') 'U_FOULED=', U_f
write(*,'(A,F18.4)') 'H_HOT=', hs_h
write(*,'(A,F18.4)') 'H_COLD=', hs_c
write(*,'(A,F18.4)') 'RE_HOT=', Re_h
write(*,'(A,F18.4)') 'RE_COLD=', Re_c
write(*,'(A,F18.4)') 'V_HOT=', vch_h
write(*,'(A,F18.4)') 'V_COLD=', vch_c
write(*,'(A,F18.4)') 'DP_HOT_TOTAL=', dP_h
write(*,'(A,F18.4)') 'DP_COLD_TOTAL=', dP_c
write(*,'(A,F18.4)') 'DP_HOT_PORT=', dP_port_h
write(*,'(A,F18.4)') 'DP_COLD_PORT=', dP_port_c
write(*,'(A,F18.4)') 'DP_HOT_CHAN=', dP_h - dP_port_h
write(*,'(A,F18.4)') 'DP_COLD_CHAN=', dP_c - dP_port_c
contains
subroutine calculate_phe(Np_val, hs_h, hs_c, dP_h, dP_c, &
dP_port_h, dP_port_c, U_f, U_c, A_prov_val, A_req_val, &
Re_h_val, Re_c_val, vch_h_val, vch_c_val)
integer, intent(in) :: Np_val
double precision, intent(out) :: hs_h, hs_c, dP_h, dP_c, dP_port_h, dP_port_c, U_f, U_c, A_prov_val, A_req_val
double precision, intent(out) :: Re_h_val, Re_c_val, vch_h_val, vch_c_val
integer :: Nch
double precision :: Ach, Pr_h, Pr_c, Kf_h, p_h, f_turb_h, f_lam_h, f_h
double precision :: Kf_c, p_c, f_turb_c, f_lam_c, f_c
double precision :: beta_rad, Nu_h, Nu_c, Rw
double precision :: U_clean_inv, U_fouled_inv, d_port, v_port_h, v_port_c
double precision :: dP_channel_h, dP_channel_c
Nch = (Np_val - 1) / 2
if (Nch < 1) Nch = 1
Ach = b_depth * W
vch_h_val = m_h / (rho_h * Ach * dble(Nch))
vch_c_val = m_c / (rho_c * Ach * dble(Nch))
Re_h_val = rho_h * vch_h_val * Dh / mu_h
Re_c_val = rho_c * vch_c_val * Dh / mu_c
Pr_h = mu_h * Cp_h / k_h
Pr_c = mu_c * Cp_c / k_c
! Chevron friction factors smooth continuous fit
Kf_h = 0.15d0 + 0.08d0 * exp(0.06d0 * beta_deg)
p_h = 0.30d0 - 0.002d0 * beta_deg
f_turb_h = Kf_h * Re_h_val**(-p_h)
f_lam_h = (10.0d0 + 0.7d0 * beta_deg) / Re_h_val
f_h = f_turb_h
if (f_lam_h > f_h) f_h = f_lam_h
Kf_c = 0.15d0 + 0.08d0 * exp(0.06d0 * beta_deg)
p_c = 0.30d0 - 0.002d0 * beta_deg
f_turb_c = Kf_c * Re_c_val**(-p_c)
f_lam_c = (10.0d0 + 0.7d0 * beta_deg) / Re_c_val
f_c = f_turb_c
if (f_lam_c > f_c) f_c = f_lam_c
! Nusselt numbers
beta_rad = beta_deg * pi / 180.0d0
Nu_h = 0.205d0 * Pr_h**(1.0d0/3.0d0) * (mu_h / mu_wh)**0.14d0 * &
(f_h * Re_h_val**2 * sin(2.0d0 * beta_rad))**(1.0d0/3.0d0)
Nu_c = 0.205d0 * Pr_c**(1.0d0/3.0d0) * (mu_c / mu_wc)**0.14d0 * &
(f_c * Re_c_val**2 * sin(2.0d0 * beta_rad))**(1.0d0/3.0d0)
hs_h = Nu_h * k_h / Dh
hs_c = Nu_c * k_c / Dh
Rw = t_plate / kw
U_clean_inv = 1.0d0 / hs_h + Rw + 1.0d0 / hs_c
U_c = 1.0d0 / U_clean_inv
U_fouled_inv = 1.0d0 / hs_h + Rfh + Rw + Rfc + 1.0d0 / hs_c
U_f = 1.0d0 / U_fouled_inv
A_req_val = Q / (U_f * LMTD_cf)
A_prov_val = dble(Np_val - 2) * W * Lp * phi
! Port velocity
d_port = d_port_mm / 1000.0d0
v_port_h = m_h / (rho_h * pi * d_port**2 / 4.0d0)
v_port_c = m_c / (rho_c * pi * d_port**2 / 4.0d0)
! Pressure drops
dP_channel_h = f_h * (Lp / Dh) * (rho_h * vch_h_val**2 / 2.0d0) * (mu_h / mu_wh)**(-0.14d0)
dP_channel_c = f_c * (Lp / Dh) * (rho_c * vch_c_val**2 / 2.0d0) * (mu_c / mu_wc)**(-0.14d0)
dP_port_h = 1.4d0 * 1.0d0 * (rho_h * v_port_h**2 / 2.0d0)
dP_port_c = 1.4d0 * 1.0d0 * (rho_c * v_port_c**2 / 2.0d0)
dP_h = dP_channel_h + dP_port_h
dP_c = dP_channel_c + dP_port_c
end subroutine calculate_phe
end program hx_plate
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 hx_plate.f90 -o hx_plate
2. Execution with input.txt redirection:
hx_plate < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 29 90.0 40.0 20.0 35.0 2.0 3.0 990.0 4180.0 0.0008 0.6 0.0008 998.0 4180.0 0.001 0.6 0.001 0.35 0.85 0.0005 0.003 0.008 45.0 50.0 150.0 0.0001 0.0001 50.0
Parameter Description:
Mode (1=Verify / Rating, 2=Sizing) Initial plate count Nplates Hot fluid inlet temp [°C] Hot fluid outlet temp [°C] Cold fluid inlet temp [°C] Cold fluid outlet temp [°C] Hot fluid mass flow rate [kg/s] Cold fluid mass flow rate [kg/s] Hot fluid density [kg/m3] Hot fluid specific heat [J/kg-K] Hot fluid viscosity [Pa-s] Hot fluid conductivity [W/m-K] Hot wall viscosity [Pa-s] Cold fluid density [kg/m3] Cold fluid specific heat [J/kg-K] Cold fluid viscosity [Pa-s] Cold fluid conductivity [W/m-K] Cold wall viscosity [Pa-s] Plate width W [m] Thermal height Lp [m] Plate thickness t_plate [m] Pressing depth b [m] Corrugation wavelength pc [m] Chevron angle beta [deg] Port diameter [mm] Max allowable pressure drop [kPa] Hot side fouling factor [m2-K/W] Cold side fouling factor [m2-K/W] Plate wall thermal conductivity [W/m-K]