🔬
Solver Purpose & Physical Scope
Model aqua-ammonia (NH3-H2O) and lithium bromide (LiBr-H2O) absorption refrigeration systems. Evaluates cycle thermodynamic state properties, circulation ratio, heat duty at each component (generator, condenser, evaporator, absorber), pump work, cooling Coefficient of Performance (COP), Carnot COP limit, and second-law exergy efficiency.
📂 Discipline: Thermo
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 494 times
📄 Source File:
absorption_refrig.f90
📁 calcul/Thermodynamics /
absorption_refrig.f90
program absorption_refrig
implicit none
integer :: system_type, i, iostat_val
double precision :: T_gen, T_abs, T_cond, T_evap
double precision :: Q_gen, Q_abs, Q_cond, Q_evap, W_pump
double precision :: COP, COP_carnot, eta_exergy
double precision :: x_strong, x_weak, f_circ, m_ref, m_sol
double precision :: h_ref_vap, h_ref_liq, h_sol_strong, h_sol_weak
double precision :: P_high, P_low, dT_overlap
double precision :: T_i, COP_i, f_i
character(len=40) :: system_name
read(*,*,iostat=iostat_val) system_type
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid system type input.'
stop
end if
read(*,*,iostat=iostat_val) T_gen
read(*,*,iostat=iostat_val) T_abs
read(*,*,iostat=iostat_val) T_cond
read(*,*,iostat=iostat_val) T_evap
read(*,*,iostat=iostat_val) m_ref
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all absorption cycle inputs.'
stop
end if
if (T_gen<=0.0d0.or.T_abs<=0.0d0.or.T_cond<=0.0d0.or.T_evap<=0.0d0) then
write(*,*) 'ERROR: All temperatures must be positive (K).'
stop
end if
if (T_gen <= T_cond) then
write(*,*) 'ERROR: Generator temp must exceed condenser temp.'
stop
end if
if (T_evap >= T_abs) then
write(*,*) 'ERROR: Evaporator temp must be below absorber temp.'
stop
end if
if (m_ref <= 0.0d0) m_ref = 1.0d0
select case(system_type)
case(1)
system_name = 'Aqua-Ammonia (NH3-H2O)'
call aqua_ammonia_cycle(T_gen, T_abs, T_cond, T_evap, m_ref, &
x_strong, x_weak, f_circ, Q_gen, Q_abs, Q_cond, Q_evap, &
W_pump, COP, P_high, P_low, &
h_ref_vap, h_ref_liq, h_sol_strong, h_sol_weak)
case(2)
system_name = 'LiBr-H2O (water is refrigerant)'
call libr_h2o_cycle(T_gen, T_abs, T_cond, T_evap, m_ref, &
x_strong, x_weak, f_circ, Q_gen, Q_abs, Q_cond, Q_evap, &
W_pump, COP, P_high, P_low, &
h_ref_vap, h_ref_liq, h_sol_strong, h_sol_weak)
case default
write(*,*) 'ERROR: System type must be 1 (NH3-H2O) or 2 (LiBr-H2O).'
stop
end select
! Carnot COP for absorption
! COP_Carnot = (T_evap/(T_cond-T_evap)) * ((T_gen-T_abs)/T_gen)
COP_carnot = (T_evap/(T_cond - T_evap)) * ((T_gen - T_abs)/T_gen)
if (COP_carnot < 0.0d0) COP_carnot = 0.0d0
eta_exergy = COP / max(COP_carnot, 1.0d-10)
dT_overlap = T_abs - T_evap ! temperature lift
write(*,'(A)') '============================================================'
write(*,'(A)') ' ABSORPTION REFRIGERATION CYCLE ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUTS --------------------------------------------------'
write(*,'(A,A)') ' System Type = ', trim(system_name)
write(*,'(A,F12.2,A)') ' Generator Temperature = ', T_gen, ' K'
write(*,'(A,F12.2,A)') ' Absorber Temperature = ', T_abs, ' K'
write(*,'(A,F12.2,A)') ' Condenser Temperature = ', T_cond, ' K'
write(*,'(A,F12.2,A)') ' Evaporator Temperature = ', T_evap, ' K'
write(*,'(A,ES12.4,A)') ' Refrigerant Flow = ', m_ref, ' kg/s'
write(*,*)
write(*,'(A)') '--- SOLUTION PROPERTIES -------------------------------------'
write(*,'(A,F10.4)') ' Strong Solution x_s = ', x_strong
write(*,'(A,F10.4)') ' Weak Solution x_w = ', x_weak
write(*,'(A,F10.4)') ' Circulation Ratio f = ', f_circ
write(*,'(A,ES12.4,A)') ' Solution Flow (strong) = ', m_ref*f_circ, ' kg/s'
write(*,'(A,ES12.4,A)') ' High Pressure = ', P_high, ' Pa'
write(*,'(A,ES12.4,A)') ' Low Pressure = ', P_low, ' Pa'
write(*,*)
write(*,'(A)') '--- HEAT DUTIES ---------------------------------------------'
write(*,'(A,ES12.4,A)') ' Generator Heat Qgen = ', Q_gen, ' W'
write(*,'(A,ES12.4,A)') ' Absorber Heat Qabs = ', Q_abs, ' W'
write(*,'(A,ES12.4,A)') ' Condenser Heat Qcond = ', Q_cond, ' W'
write(*,'(A,ES12.4,A)') ' Evaporator Heat Qevap = ', Q_evap, ' W'
write(*,'(A,ES12.4,A)') ' Pump Work Wpump = ', W_pump, ' W'
write(*,*)
write(*,'(A)') '--- PERFORMANCE ---------------------------------------------'
write(*,'(A,F10.4)') ' COP (cooling) = ', COP
write(*,'(A,F10.4)') ' COP Carnot (absorption) = ', COP_carnot
write(*,'(A,F10.4)') ' Exergetic Efficiency = ', eta_exergy
write(*,'(A,F10.2,A)') ' Temperature Lift = ', dT_overlap, ' K'
write(*,*)
write(*,'(A)') ' Energy Balance Check:'
write(*,'(A,ES12.4,A)') ' Qgen + Qevap + Wpump = ', Q_gen+Q_evap+W_pump, ' W'
write(*,'(A,ES12.4,A)') ' Qcond + Qabs = ', Q_cond+Q_abs, ' W'
write(*,*)
! COP vs Generator Temperature sweep
write(*,'(A)') '--- COP VS GENERATOR TEMPERATURE SWEEP ----------------------'
write(*,'(A)') ' T_gen[K] COP f_circ COP_Carnot'
write(*,'(A)') ' -----------------------------------------------------------'
do i = 1, 50
T_i = T_cond + 5.0d0 + (T_gen + 50.0d0 - T_cond - 5.0d0)*dble(i-1)/49.0d0
select case(system_type)
case(1)
call aqua_ammonia_cop(T_i, T_abs, T_cond, T_evap, COP_i, f_i)
case(2)
call libr_h2o_cop(T_i, T_abs, T_cond, T_evap, COP_i, f_i)
end select
write(*,'(F12.2,2X,F10.4,2X,F10.4,2X,F10.4)') T_i, COP_i, f_i, &
(T_evap/(T_cond-T_evap))*((T_i-T_abs)/T_i)
end do
write(*,*)
! COP vs Evaporator Temperature sweep
write(*,'(A)') '--- COP VS EVAPORATOR TEMPERATURE SWEEP ---------------------'
write(*,'(A)') ' T_evap[K] COP COP_Carnot'
write(*,'(A)') ' -------------------------------------------'
do i = 1, 40
T_i = 240.0d0 + (T_abs - 5.0d0 - 240.0d0)*dble(i-1)/39.0d0
select case(system_type)
case(1)
call aqua_ammonia_cop(T_gen, T_abs, T_cond, T_i, COP_i, f_i)
case(2)
call libr_h2o_cop(T_gen, T_abs, T_cond, T_i, COP_i, f_i)
end select
write(*,'(F12.2,2X,F10.4,2X,F10.4)') T_i, COP_i, &
(T_i/max(T_cond-T_i,1.0d0))*((T_gen-T_abs)/T_gen)
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' COP = Qevap / (Qgen + Wpump).'
write(*,'(A)') ' COP_Carnot = (Te/(Tc-Te))*(Tg-Ta)/Tg.'
write(*,'(A)') ' f = x_strong / (x_strong - x_weak).'
write(*,'(A)') ' Simplified property correlations for educational use.'
contains
subroutine aqua_ammonia_cycle(Tg, Ta, Tc, Te, mref, xs, xw, f, &
Qg, Qa, Qc, Qe, Wp, cop, Ph, Pl, hrv, hrl, hss, hsw)
implicit none
double precision, intent(in) :: Tg, Ta, Tc, Te, mref
double precision, intent(out) :: xs, xw, f, Qg, Qa, Qc, Qe, Wp, cop
double precision, intent(out) :: Ph, Pl, hrv, hrl, hss, hsw
double precision :: h_evap_nh3, Tg_C, Ta_C, Tc_C, Te_C
Tg_C = Tg - 273.15d0
Ta_C = Ta - 273.15d0
Tc_C = Tc - 273.15d0
Te_C = Te - 273.15d0
! Simplified NH3 saturation pressures
Ph = 1.0d6 * exp(11.0d0 - 3100.0d0/Tc) ! condenser
Pl = 1.0d6 * exp(11.0d0 - 3100.0d0/Te) ! evaporator
if (Ph < Pl) Ph = Pl * 2.0d0
! Simplified concentration correlations
xs = 0.60d0 - 0.002d0*(Ta_C - 30.0d0) ! strong solution (rich in NH3)
xw = 0.30d0 + 0.002d0*(Tg_C - 100.0d0) ! weak solution (lean)
if (xs > 0.70d0) xs = 0.70d0
if (xs < 0.30d0) xs = 0.30d0
if (xw > xs - 0.02d0) xw = xs - 0.02d0
if (xw < 0.05d0) xw = 0.05d0
f = xs / max(xs - xw, 1.0d-10)
! Simplified enthalpies (kJ/kg → J/kg)
h_evap_nh3 = 1370.0d3 ! latent heat NH3
hrv = 1450.0d3 + 2100.0d0*(Te_C - (-33.0d0)) ! vapor from evaporator
hrl = 300.0d3 + 4500.0d0*(Tc_C - 25.0d0) ! liquid from condenser
hss = -100.0d3 + 3800.0d0*(Ta_C - 25.0d0) ! strong solution
hsw = 100.0d3 + 4200.0d0*(Tg_C - 100.0d0) ! weak solution
Qe = mref * h_evap_nh3
Qc = mref * (hrv - hrl)
Qg = mref * (hrv - hrl) + mref * f * (hsw - hss)
if (Qg < Qe) Qg = Qe * 1.5d0
Qa = Qg + Qe - Qc
Wp = mref * f * 0.001d0 * (Ph - Pl) / 0.80d0 ! pump work
cop = Qe / max(Qg + Wp, 1.0d-10)
end subroutine aqua_ammonia_cycle
subroutine libr_h2o_cycle(Tg, Ta, Tc, Te, mref, xs, xw, f, &
Qg, Qa, Qc, Qe, Wp, cop, Ph, Pl, hrv, hrl, hss, hsw)
implicit none
double precision, intent(in) :: Tg, Ta, Tc, Te, mref
double precision, intent(out) :: xs, xw, f, Qg, Qa, Qc, Qe, Wp, cop
double precision, intent(out) :: Ph, Pl, hrv, hrl, hss, hsw
double precision :: Tg_C, Ta_C, Tc_C, Te_C
Tg_C = Tg - 273.15d0
Ta_C = Ta - 273.15d0
Tc_C = Tc - 273.15d0
Te_C = Te - 273.15d0
! Water saturation pressures (simplified Antoine)
Ph = 610.78d0 * exp(17.27d0*Tc_C/(Tc_C+237.3d0))
Pl = 610.78d0 * exp(17.27d0*Te_C/(Te_C+237.3d0))
if (Ph < Pl) Ph = Pl * 2.0d0
! LiBr concentrations (mass fraction LiBr)
! Strong = high LiBr (leaves generator), Weak = low LiBr (leaves absorber)
xw = 0.55d0 + 0.001d0*(Ta_C - 30.0d0)
xs = 0.62d0 + 0.001d0*(Tg_C - 80.0d0)
if (xs > 0.70d0) xs = 0.70d0
if (xs < xw + 0.02d0) xs = xw + 0.02d0
if (xw < 0.40d0) xw = 0.40d0
! f = xs/(xs-xw) for LiBr system (based on LiBr mass balance)
f = xs / max(xs - xw, 1.0d-10)
! Simplified water enthalpies
hrv = 2500.0d3 + 1900.0d0*(Te_C - 5.0d0)
hrl = 105.0d3 + 4186.0d0*(Tc_C - 25.0d0)
hss = 150.0d3 + 2000.0d0*(Tg_C - 80.0d0)
hsw = 80.0d3 + 2200.0d0*(Ta_C - 30.0d0)
Qe = mref * (hrv - hrl)
Qc = mref * (hrv - hrl)
Qg = mref * hrv + mref*(f-1.0d0)*hss - mref*f*hsw
if (Qg < 0.0d0) Qg = Qe * 1.4d0
Qa = Qg + Qe - Qc
Wp = mref * f * 0.001d0 * (Ph - Pl) / 0.75d0
cop = Qe / max(Qg + Wp, 1.0d-10)
end subroutine libr_h2o_cycle
subroutine aqua_ammonia_cop(Tg, Ta, Tc, Te, cop_out, f_out)
implicit none
double precision, intent(in) :: Tg, Ta, Tc, Te
double precision, intent(out) :: cop_out, f_out
double precision :: xs, xw, Qg, Qa, Qc, Qe, Wp, Ph, Pl
double precision :: hrv, hrl, hss, hsw
call aqua_ammonia_cycle(Tg, Ta, Tc, Te, 1.0d0, xs, xw, f_out, &
Qg, Qa, Qc, Qe, Wp, cop_out, Ph, Pl, hrv, hrl, hss, hsw)
end subroutine aqua_ammonia_cop
subroutine libr_h2o_cop(Tg, Ta, Tc, Te, cop_out, f_out)
implicit none
double precision, intent(in) :: Tg, Ta, Tc, Te
double precision, intent(out) :: cop_out, f_out
double precision :: xs, xw, Qg, Qa, Qc, Qe, Wp, Ph, Pl
double precision :: hrv, hrl, hss, hsw
call libr_h2o_cycle(Tg, Ta, Tc, Te, 1.0d0, xs, xw, f_out, &
Qg, Qa, Qc, Qe, Wp, cop_out, Ph, Pl, hrv, hrl, hss, hsw)
end subroutine libr_h2o_cop
end program absorption_refrig
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 absorption_refrig.f90 -o absorption_refrig
2. Execution with input.txt redirection:
absorption_refrig < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 403.15 308.15 308.15 263.15 1.0
Parameter Description:
Working pair (1=NH3-H2O, 2=LiBr-H2O) Generator temperature T_gen [K] Absorber temperature T_abs [K] Condenser temperature T_cond [K] Evaporator temperature T_evap [K] Refrigerant mass flow rate [kg/s]