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. Pressure drop limit may be too restrictive.' 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