Impact of polygeneration topology on the technoeconomic performance of green hydrogen production utilizing integrated CPVT systems

System topology descriptions

The integrated concentrating solar polygeneration system consists of four interconnected subsystems: (1) a parabolic trough solar concentrating collector with a flat thermal receiver available in two technology variants (CPVT / C3JT); (2) a parabolic trough solar concentrating collector with a flat thermal receiver for thermal heating (PTC); (3) a proton exchange membrane electrolyzer for green hydrogen production; and (4) a multi-effect evaporation desalination unit (MEE). Figures 1 and 2 show the system schematic for polygeneration-I and polygeneration-II. The system is designed for distributed deployment at Egyptian Sea coast locations with all thermal energy co-products directed to the desalination unit.

For Polygeneration-I, the system operates as a dual output polygeneration plant. The parabolic trough collector delivers concentrated solar flux to the receiver, which simultaneously generates electrical power and extracts useful heat through the coolant loop. The full electrical output is routed to the PEM electrolyzer that receives external deionized water as feedstock and consumes the full available electrical output of the receiver to produce hydrogen. The MEE unit operates in parallel, independently consuming the receiver thermal output to produce freshwater. The thermal output is directed to the MEE desalination unit, where it serves as the motive for heat source driving successive evaporation effects. Both hydrogen and freshwater are produced continuously during daylight hours, with production rates governed by the instantaneous solar resource at each location.

For Polygeneration-II, the MEE subsystem operates identically to polygeneration-I except that the distillated produced by the MEE unit is fed directly into the PEM electrolyzer as the hydrogen production feedstock. At each hourly timestep, the electrical power required to electrolyze the MEE distillate flow rate is computed and the remaining electrical power is exported. The hydrogen production rate is therefore controlled by the MEE distillate output.

Two receiver technology variants are evaluated in parallel, defining the A/B subscript throughout the case naming:

CPVT receiver (Cases A): a linear parabolic trough concentrating collector focuses direct normal irradiance onto a flat receiver incorporating silicon-based multi-crystalline photovoltaic cells laminated to the absorber surface. A single-pass water coolant loop extracts waste heat from the receiver, maintaining cell temperatures at acceptable operating levels while simultaneously generating usable thermal output. The geometrical concentration ratio of 32x is moderate, consistent with established parabolic trough CPVT designs. Multi-crystalline silicon cells were selected instead of mono-crystalline cells for cost reasons. At a low concentration ratio of 32x, the system performance is mainly limited by optical and thermal losses rather than by cell efficiency, so the lower cost of multi-crystalline silicon makes it the more practical choice for a distributed-scale system.

C3JT receiver (Cases B): the same parabolic trough geometry is equipped with a high-efficiency GaInP/GaInAs/Ge triple-junction receiver operating at the same geometric concentration ratio. The triple-junction cell exploits current matching across three series-connected sub-cells, capturing a substantially broader portion of the solar spectrum than any single-junction device.

Schematic of the integrated concentrating solar polygeneration-I system.

Schematic of the integrated concentrating solar polygeneration-II system.

Parabolic trough collector optical model

The solar collection subsystem is a single-axis north-south tracking parabolic trough. The aperture area and incident solar power are determined by the collector geometry. The key geometric and optical parameters are given in Table 1 for both receiver variants. The instantaneous optical power delivered to the receiver is:

where DNI(t) is the direct normal irradiance (W m^− 2) and A_ap is the aperture area of the PTC.

The overall optical efficiency of the concentrating system, \(\:{\eta\:}_{\text{o}\text{p}\text{t}}\left(\theta\:\right)\), accounts for mirror reflectivity ρ, cover transmissivity τ, receiver absorptivity α, incidence angle modifier IAM(θ), and intercept factor γ and is given by:

The geometric concentration ratio \(\:{CR}_{g}\) relates the aperture width (W) to the receiver width (w) and is given by:

Concentrated photovoltaic model

The CPVT receiver photovoltaic output is computed using a double-diode equivalent circuit, which is the standard framework for silicon-based concentrating PV cells and has been validated in CPVT system dynamic simulations¹⁹. The current-voltage equation at cell temperature T_c and irradiance G is^7,23,24:

where, \(\:I\) is the cell current, \(\:{I}_{ph}\) is the photocurrent, \(\:{I}_{0,i}\) is the diode saturation current,, \(\:V\) is the cell voltage, \(\:{R}_{s}\) and \(\:{R}_{sh}\) are the series and shunt resistances, \(\:{n}_{i}\) is the diode ideality factor (\(\:{n}_{1}=1.5\), \(\:{n}_{2}=2)\), \(\:{V}_{t}\) is thermal voltage at cell temperature and equals to \(\:\frac{k\:{T}_{c}}{q}\), \(\:k\) is Boltzmann’s constant, \(\:{T}_{c}\) is the cell temperature, \(\:q\) is the electron charge, \(\:{I}_{sc,ref}\) is the reference short-circuit current at standard test conditions and \(\:{\alpha\:}_{Isc}\) is the temperature coefficient of short-circuit current.

The photogenerated current is corrected for irradiance deviation from standard test conditions (\(\:{DNI}_{ref}\:\)= 1,000 W m^− 2, \(\:{T}_{ref}\) = 25 °C) and for temperature²⁴:

where \(\:{I}_{sc,ref}\) is the reference short-circuit current at standard test conditions and \(\:{\alpha\:}_{Isc}\) is the temperature coefficient of short-circuit current.

The diode reverse saturation current is the most temperature-sensitive parameter in the model and is described as follows²⁴:

where, \(\:{I}_{0,\text{\:ref,}i}\:\)is the saturation current at reference temperature, \(\:{E}_{g}\left({T}_{c}\right)\) is the semiconductor bandgap energy at cell temperature (E_g,0 K = 1.12 eV) and \(\:{T}_{ref}\) is the reference temperature in kelvin.

Triple-junction photovoltaic model

The C3JT receiver employs a double-diode model for each of three series-connected sub-cells (j = 1,2 and 3, where top is 1; GaInP, middle is 2; GaInAs and bottom is 3; Ge), following the approach and validated against the experimental datasets of Helmers et al.¹⁶ and under concentration validation against Catelani et al.⁷. The double-diode formulation for sub-cell (j) is:

with,

where \(\:{n}_{1,\varvec{j}}\) and \(\:{n}_{2,\varvec{j}}\) is the diode ideality factor at each sub-cell. For top cell \(\:{n}_{\text{1,1}}=1.5\) and \(\:{n}_{\text{2,1}}=2\), mid cell \(\:{n}_{\text{1,2}}=1.97\) and \(\:{n}_{\text{2,2}}=2\) and bottom cell \(\:{n}_{\text{1,3}}=1\) and \(\:{n}_{\text{2,3}}=2\) that is consistent with the near-ideal behavior of the Ge bottom sub-cell confirmed experimentally by Roensch et al.¹⁸. These values are consistent with the parameterization of Catelani et al.⁷ and Anaty et al.⁶ adopted in the present model.

The series connection imposes current matching: the terminal current equals the minimum sub-cell photocurrent, and the terminal voltage is the algebraic sum of all sub-cell voltages:

The sub-cell open-circuit voltage dependence on temperature and geometric concentration ratio follows the analytical model of Helmers et al.¹⁶:

where \(\:{\beta\:}_{Voc,\varvec{j}}\) is the measured subcell-specific V_OC temperature coefficient (mV K^− 1): -2.8 for GaInP, -2.4 for GaInAs, and − 2.1 for Ge¹⁶. Bandgap temperature dependence for each sub-cell uses the Varshni equation²⁶, with material parameters for GaInP (E_g,0 K = 1.90 eV), GaInAs (1.49 eV), and Ge (0.75 eV)¹⁵.

Receiver thermal model

The energy balance on the receiver aperture couples optical input, electrical generation, thermal extraction, and thermal loss:

Where \(\:{S}_{absorbed}\) is the absorbed solar power, \(\:{P}_{electric}\) is electric power generated, \(\:{Q}_{th,useful}\) is the recoverable thermal energy and \(\:\:{Q}_{loss}\) is thermal losses energy.

The recoverable thermal power extracted from the coolant fluid is:

Receiver thermal loss to the environment combines convective and radiative contributions and is modelled as:

where U_L is the overall heat loss coefficient (W m^− 2 K^− 1), A_rec is the receiver outer surface area, \(\:{T}_{c}\) is the mean cell or tube temperature, and \(\:{T}_{amb}\) is the ambient temperature.

The coolant outlet temperature \(\:{T}_{f,out}\)_t is determined iteratively from the coupled heat transfer and PV electrical model at each simulation timestep. The overall heat loss coefficient U_L is computed dynamically at each hourly timestep within the SIMSCAPE thermal domain, as a function of the instantaneous wind speed, ambient temperature, and receiver surface temperature. This dynamic treatment ensures that seasonal variations in wind speed and ambient conditions are fully captured in the thermal loss calculation throughout the year.

PEM electrolyzer model

The PEM electrolyzer model follows the comprehensive electrochemical-thermal framework developed and validated against experimental data by Bessarabov et al.²⁷ and Aouali et al.⁸ and implemented in the MATLAB/SIMSCAPE environment. The electrolyzer cell voltage is expressed as the sum of the Nernst reversible equilibrium potential and three overpotential contributions:

Where \(\:{V}_{cell}\) is the electrolyzer cell voltage, \(\:{E}_{rev}\) is the reversible potential, \(\:{\eta\:}_{act}\) is the activation overpotential, \(\:{\eta\:}_{ohm}\) is the ohmic overpotential and \(\:{\eta\:}_{conc}\) is the concentration overpotential.

The temperature reversible Nernst potential accounts for operating temperature \(\:{T}_{el}\) and product gas partial pressures²⁷:

where, \(\:{\text{E}}_{th}^{0}\) is standard reversible potential (1.229 V at 25 °C), \(\:{R}_{u}\) is universal gas constant (8.314 J/(mol⋅K)), \(\:{T}_{el}\) is cell operating temperature (kelvin), \(\:F\) is Faraday’s constant (96485 C/mol), \(\:{p}_{{\text{H}}_{2}},{p}_{{\text{O}}_{2}}\) are partial pressures of hydrogen and oxygen (bar), \(\:{a}_{{\text{H}}_{2}\text{O}}\) is activity of liquid water (≈ 1 for pure water).

Activation overpotential at the anode (a) and cathode (c) is modelled by:

where, the anode or cathode activation overpotential is calculated as²⁸:

where, \(\:x\) stands for anode or cathode, \(\:{\alpha\:}_{x}\) is the anode (\(\:{\alpha\:}_{a}=0.51\)) or cathode (\(\:{\alpha\:}_{c}=0.627\)) charge transfer coefficient, \(\:{z}_{x}\) is the stoichiometric coefficient of electrons for the anode (\(\:{z}_{a}=3\)) or cathode (\(\:{z}_{c}=-1\)) reaction, \(\:i\) is the operating current density (A/cm²), \(\:{i}_{0,x}\) is the exchange current density at the anode or cathode (A/cm²).

The anode or cathode exchange current density (\(\:{i}_{0,x}\)) is a critical parameter that is highly dependent on temperature and catalyst properties. It is modeled using an Arrhenius-type expression²⁷:

where, \(\:{i}_{0,\text{\:ref\:},x}\) is reference exchange current density for the anode (\(\:{i}_{0,\text{\:ref\:},a}=1{\times\:10}^{-7}\)) or cathode (\(\:{i}_{0,\text{\:ref\:},c}=1{\times\:10}^{-3}\)), \(\:{C}_{{H}_{2}O,\text{\:ratio\:}}\) is the ratio of water concentration to reference concentration anode (\(\:{C}_{{H}_{2}O,\text{\:ratio\:}}=0.22\)), \(\:{E}_{excit,x}\) is the excitation energy for the anode (\(\:{E}_{excit,a}=76{\times\:10}^{3}\)) or cathode (\(\:{E}_{excit,c}=4.3{\times\:10}^{3}\)) reaction (J/mol), \(\:{T}_{\text{ref\:}}\) is the reference temperature for the excitation parameters and equals to 353 K (K).

And the catalyst surface area factor (\(\:{\gamma\:}_{m,x}\)) for the platinum catalyst is calculated as^27,29:

where, \(\:{\phi\:}_{x}\) is the fraction of active metal catalyst for anode or cathode anode (\(\:{\phi\:}_{x}=0.75\)), \(\:{m}_{cat,x}\) is the catalyst loading at the anode (\(\:{m}_{cat,a}=1{\times\:10}^{-3}\)) or cathode (\(\:{m}_{cat,c}=0.3{\times\:10}^{-3}\)) (g/cm²), \(\:{\rho\:}_{x}\) is the density of the anode which is platinum catalyst or of the cathode which is iridium catalyst (g/cm³), \(\:{d}_{cat,x}\) is the catalyst crystallite diameter for anode (\(\:{d}_{cat,a}=2.9{\times\:10}^{-7}\)) or cathode (\(\:{d}_{cat,a}=2.7{\times\:10}^{-7}\)) (cm).

The ohmic overpotential accounts for the voltage losses due to resistance to the flow of charge within the electrolyzer cell

The ionic overpotential \(\:{\eta\:}_{Ionic}\) arises from the resistance of the membrane to proton transport. The electronic overpotential, which accounts for the resistance in the solid components of the cell, is typically much smaller than the ionic overpotential and it is assumed to be a small, fixed fraction of the ionic overpotential equals to 1% of ionic overpotential. Ionic overpotential is calculated as:

where, \(\:i\) is the operating current density (A/cm²), \(\:{\text{t}}_{memb}^{total}\) is the total thickness of the membrane and gas diffusion layers (cm), \(\:{\sigma\:}_{memb}\) is the ionic conductivity of the membrane (S/cm).

The membrane conductivity (\(\:{\sigma\:}_{memb}\)) is a sensitive parameter that depends on the membrane’s hydration level and operating temperature²⁹

where, \(\:{\delta\:}_{\text{memb\:}}\) is a dimensionless fitting constant for the conductivity model (\(\:{\delta\:}_{\text{memb\:}}=1.65\)), \(\:{E}_{\text{viscosity\:},{H}_{2}O}\) is the activation energy related to water viscosity (\(\:{E}_{\text{viscosity\:},{H}_{2}O}=14000\)) (J/mol), \(\:{\epsilon}_{\text{mem\:}}\), is the membrane porosity which is a function of the number of sorbed water molecules per \(\:{\lambda\:}_{\text{sorbed\:}}\) acid site (\(\:{\lambda\:}_{\text{sorbed\:}}=20\:,\:{\lambda\:}_{0,\text{sorbed\:}}=1.8\)) and the ratio of membrane molar volume to water (r = 30)^27,30:

\(\:\beta\:\) is the degree of acid site dissociation, given by:

And \(\:{K}_{\text{acid\:}}\) is the dissociation equilibrium constant^27,30:

Concentration total overpotential at high current densities is calculated as:

The concentration overpotential at the anode or cathode is calculated as²⁷:

where, \(\:{i}_{limit}\) is the limiting current density (A/cm²), a key parameter representing the maximum current that can be sustained before mass transport limitations become critical (\(\:{i}_{limit}=4\:\text{A}/\text{c}\text{m}^{2}\)).

Multi-effect evaporation desalination model

The MEE desalination unit is modelled as a forward-feed system with N_eff = 6 serially connected effects operating at successively decreasing temperatures driven by the receiver thermal output. The MEE desalination unit was implemented as a physics-based dynamic model within the MATLAB/SIMSCAPE environment, constructed using standard blocks from the SIMSCAPE Fluids Two-Phase Fluid library. This approach provides a first-principles representation of the core thermodynamic phenomena of evaporation and condensation which govern MEE plant performance, rather than relying on simplified steady-state performance correlations. The primary output of the model at each simulation timestep is the instantaneous freshwater production rate, computed dynamically as a function of the thermal energy supplied by the solar field. The heat transfer within each effect is computed using the SIMSCAPE Two-Phase Fluid library’s internal two-phase convection model, which accounts for boiling and condensation phenomena beyond the single-phase Nusselt correlation of the pipe block. The complete MEE SIMSCAPE model was validated against the performance data of El-Dessouky and Ettouney³¹ before its integration into the polygeneration simulation platform, confirming the reliability of the mass and energy balance calculations across all six effects.

The model was built with a modular architecture in which each effect of the MEE train was developed as an independent SIMSCAPE subsystem. The six subsystems were then connected in series to form the complete multi-effect train, mirroring the physical construction of a real forward-feed MEE plant and allowing straightforward reconfiguration for a different number of effects. Within each effect, four physical processes are represented by dedicated SIMSCAPE sub-blocks: (i) simultaneous evaporation and condensation in the inter-effect heat exchanger, modelled using a pipe block with two-phase fluid domain; (ii) phase separation, implemented via a saturation properties sensor coupled with a vapor quality sensor to determine the thermodynamic state and mass fraction of the produced vapor; (iii) inter-stage thermal coupling, in which the brine and distillate thermodynamic properties computed at the outlet of each effect are used to calculate the motive heat input available to the subsequent effect; and (iv) distillate collection, in which the condensate produced in the heat exchanger of each effect is accumulated and summed to yield the total plant freshwater output as seen in Fig. 3.

Full model structure for MEE in SIMSCAPE.

Technoeconomic analysis

The techno-economic analysis translates the annual simulation outputs like hydrogen mass, freshwater, and exported electricity into standardized economic metrics that allow direct comparison of different system configurations. The methodology is based on a discounted cash flow analysis over the full plant operational lifetime n, accounting for all capital expenditure, operational costs, replacement costs, and co-product revenues.

The total capital expenditure (CAPEX) represents the upfront investment required to construct the integrated system, encompassing the solar field, concentrating receiver, PEM electrolyzer stack, MEE desalination unit, and balance-of-plant costs including installation works. Annual operational expenditure (OPEX) comprises a fixed component as scheduled maintenance, insurance, and labor, assumed as a fixed percentage of CAPEX and a variable component covering consumables and periodic component replacements such as in electrolyzer stack.

The present value of any uniform annual quantity A over the plant lifetime is calculated as:

where \(\:i\) is the real discount rate and \(\:n\) the plant lifetime in years. This formulation is applied consistently to all recurring annual cost and revenue streams.

The primary economic output is the Levelized Cost of Hydrogen (LCOH), representing the average cost of producing one kilogram of hydrogen over the plant lifetime after crediting all co-product revenues:

where: NPV (Total Plant Costs) = Total CAPEX + NPV (All Replacement Costs) + NPV(Total Annual OPEX) ; NPV (Product Revenue) = NPV(Annual Co-Product Production × Market Price of Co-Product); NPV (Hydrogen Production) = NPV(Annual Hydrogen Production in kg).

In polygeneration-I configurations, the co-product revenue term includes freshwater sales. In polygeneration-II configurations, it includes exported electricity revenues. The key financial and cost parameters adopted in this analysis are summarized in Table 3.

Finally, the polygeneration system efficiency on an HHV basis is defined as: