{'fraction_non_alphanumeric': 0.032505772638588074, 'fraction_numerical': 0.014230225356650616, 'mean_word_length': 4.799045048785551, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Design of a Pilot SOFC System for the Combined Production of Hydrogen and Electricity under Refueling Station Requirements†

Abstract The objective of the current work is to support the design of a pilot hydrogen and electricity producing plant that uses natural gas (or biomethane) as raw material, as a transition option towards a 100% renewable transportation system. The plant, with a solid oxide fuel cell (SOFC) as principal technology, is intended to be the main unit of an electric vehicle station. The refueling station has to work at different operation periods characterized by the hydrogen demand and the electricity needed for supply and self‐consumption. The same set of heat exchangers has to satisfy the heating and cooling needs of the different operation periods. In order to optimize the operating variables of the pilot plant and to provide the best heat exchanger network, the applied methodology follows a systematic procedure for multi‐objective, i.e. maximum plant efficiency and minimum number of heat exchanger matches, and multi‐period optimization. The solving strategy combines process flow modeling in steady state, superstructure‐based mathematical programming and the use of an evolutionary‐based algorithm for optimization. The results show that the plant can reach a daily weighted efficiency exceeding 60%, up to 80% when considering heat utilization.


Introduction
Up to 2016, greenhouse gas (GHG) emissions have decreased more than 22% since 1990 in the European Union (EU)-28, due to a change of the energy mix and to better energy efficiencies, technological changes and innovation [1]. This change is driven by all GHG source sectors, except transportation, which is the only fuel combustion category which increased its GHG emissions from 1990 to 2016. This reflects that fuel efficiency has not improved significantly to counteract the overall volume increase [1]. Final energy consumption in 2015 by transportation sector (excluding international maritime) represents 33% of the total share, surpassed by households and services (42%) and followed by the industry sector (25%) [2]. Road transportation contributes with 73% of the emissions, equivalent to 862 MtCO 2e per year (2015) [2]. The need to reduce transportation emissions is revealed through actions like the European Strategy for low-emission mobility [3], where there are funding allocated and measures to (i) increase transportation efficiency, (ii) speed-up the deployment of low-emission alternative energy for transportation (like H 2 and electricity), and (iii) move towards low-and zeroemission vehicles. The recent Regulation on Binding Annual GHG Emission Targets by Member States for 2021-2030 for the Sectors Not Regulated under the EU Emissions Trading System [4] (therefore, including transportation) establishes that these sectors must reduce their emissions by 30% by 2030 compared to 2005 as their contribution to the overall target.
Electro-mobility fleets comprise all type of electric vehicles, i.e., fuel cell electric vehicles (FCEVs), battery electric vehicles (BEVs) and plug-in hybrid electric vehicles (PHEVs), crucial enablers for carbon emissions reduction [5]. There exists a correlation between electric vehicles uptake and refueling infrastructure development: together with more exposure to new powertrains, the user will purchase an electric vehicle also if she/he perceives the easiness and convenience when refueling the vehicle [5]. In this context, a transition towards a renewable and more efficient transportation sector requires the installation of refueling stations of the next generation.

Hydrogen Refueling Stations and SOFC Technology
The roll-out of hydrogen refuelling station (HRS) networks has been mainly initiated in Japan, California and Germany. Currently, H 2 is generally delivered to a fuelling station in the same way as it is distributed to industry: in gaseous state, through pressurized tanks on trucks or tube trailers and pipelines, and it is usually produced from fossil fuels. Alternatively, H 2 can be produced on-site, commonly by a steam methane reformer (SMR) or an electrolyzer using preferably renewable electricity, and stored until further consumption on-site [6]. Hydrogen production (on-site or off-site) is a usual characterization criterion for HRSs, together with the state of the H 2 that is supplied to the station in an off-site configuration (gas or liquid, in which case it is stored in an onsite cryogenic tank, refilled by a liquid H 2 tanker) [7]. Every HRS responds to particular sizing and operating needs [8], as it can be dimensioned to fulfill different H 2 demands and requirements, depending on its locations and builder's criteria: for instance, seasonal, daily, hourly-following HRS demand, different selected levels of storage pressure (and thus, a number of compression steps) or the use of a cascade FCEV charging system [7,9]. Multiple HRS configurations are therefore possible [10,11]. It is recognized that HRS modeling and optimization are convenient tools to size the components of the HRS, due to different refueling profiles and possible states of charge of the FCEV [12].
Within this framework, the current work aims at optimally designing one crucial section of the HRS, the in-situ H 2 generation plant. The higher the efficiency of it, the lower the H 2 price in order to be competitive in the marketplace. As a difference from previous optimization works about in-situ H 2 generation, mainly based on SMR or electrolysis [8,[13][14][15], the current work uses a solid oxide fuel cell (SOFC)-based system not only for the production of H 2 , but also for heat and electricity. Polygeneration and modulation of SOFC-based systems have been previously studied, demonstrating feasibility, flexibility and high efficiency of multi-purpose SOFC-systems [16] and potentially of their associated supply chain, if compared with centralized and distributed SMR systems [17]. The current approach takes advantage of the SOFC high electricity efficiency, modularity (to enable a staged development), work at high temperature, capability for flexible operation at different fuel utilization (FU) values and for internal reforming.
The current work proposes a novel and transitional HRS concept, within the framework of the EU H2020 CH2P project (Cogeneration of Hydrogen, Heat and Power using solid oxide based system fed by methane rich gas). The SOFC-based CH2P system provides (i) hydrogen (H 2 ) for FCEVs, (ii) hot water that can be used in the retail station (for instance, in car wash facilities) and (iii) electricity that is required by the retail station and the HRS itself, with an excess that can be injected into the grid, stored and/or supplied to nearby stationary users or to BEVs and PHEVs. Figure 1 is a block flow diagram (BFD) that summarizes the CH2P plant layout. It consists of natural gas (NG) and water pre-processing steps, fuel prereforming, (internal reforming) SOFC stack, H 2 separation (through a combination of water-gas shift -WGS reactors and pressure swing adsorption -PSA) and combustion. The refueling station is assumed to work at discrete periods of operation characterized by different H 2 and electricity needs. The major challenge is the identification of a single heat exchanger network (HEN) that is able to operate at all the specified periods, not only without compromising the efficiency of the plant, but increasing it thanks to an optimal heat integration. A multi-period HEN provides the heating and cooling requirements for the different periods, through the same set of heat exchangers (HEXs). To reach this objective and to elucidate the best operating conditions for the CH2P plant, the novelty of the current work lies in the application of a multiobjective (MOO) and multi-period optimization approach for the conceptual design of the SOFC-based plant. Process simulation and an iterative and systematic optimization procedure [18] are used in the current paper to estimate the performance of the SOFC-based plant. The current work is the starting point of the detailed design of the plant prototype.
The following sections are structured as follows: the methodology section (Section 2) divides and explains the optimization procedure, up to the plant layout suggestion. Sections 3 and 4 focus on the application of the methodology to the CH2P plant optimization and on the description and discussion of results.

Methodology
The present work adapts the systematic and combinatorial approach for multi-period optimization and the SYNHEAT superstructure developed in [19,20], with the queuing multiobjective optimization algorithm (QMOO) from [21]. With this strategy, the aim is to optimize process design specifications as well as the HEN structure (i.e. connection among hot and cold streams), compatible to all the working periods. A similar approach has been proposed using a sequential resolution in [22].
The optimal design of the CH2P plant aims at determining the value of design variables (specifically in the current case, inlet/outlet temperatures, flow rates, FU) and the structure of the process HEN that is suitable to all the predefined operation periods. Integer (binary) variables are needed to select the final matches among hot and cold streams and utilities, among all possible connections; and a linear problem is posed taking into account the SYNHEAT superstructure to represent the possible HEX's arrangements; i.e., heat transfer flow, energy balances and particular matches' conditions.
The methodology follows the systematic approach proposed in the OSMOSE platform developed by the group of Industrial Process and Energy Systems Engineering from the É cole Polytechnique Fédérale de Lausanne (EPFL) [23], comprising the resolution of mass and energy balances, system energy integration and system performance evaluation. OSMOSE is a computer platform that bridges models and data flows from different software and tools to implement a systematic decision-making task. Figure 2 summarizes the main steps of the MOO methodology developed for the purpose of the current work. Aspen Plus software is used for the CH2P process flow modeling in steady state. The SYNHEAT superstructure stands for the type of streams connections and layout, originally proposed in [24]. The SYNHEAT algorithm from [20], is modified here to take into account the multi-period optimization strategy from [19]. The core of the MOO methodology is the QMOO algorithm developed at EPFL for energy system optimizations. It is a steady state algorithm, elitist to improve convergence and it uses clustering to preserve the diversity of the candidate population [21]. It forms the master or upper optimization level of the iterative solution process, with the subsequent steps executed by MATLAB routines (for data pre-and post-processing -including calculation of objective functions, filtering and communication among software's) as the slave or lower optimization level. MOO algorithm is used to identify the most appropriate input values of the defined decision variables, by manipulations driven by the results of the objective functions. The developed workflow comprises calls to the mathematical programming language AMPL and the mixed-integer linear programming (MILP) solver CPLEX. The results of the optimization were the starting step of the CH2P plant HEN design; from the proposed set of stream's connections, the next step was to determine a suitable HEN with common HEX areas. Successive steps derived from further CH2P project consortium discussions and decisions, and the subsequent modeling work, as practical decisions from conceptual design towards the real implementation of the plant.

Steady State Modeling of the Components of the Plant and Analysis of PFD Options
The steady state process models of the CH2P system units (see the BFD in Figure 1) are implemented in Aspen Plus software. The first process simulations are performed for three specific process flow diagrams (PFDs), so as to have a preliminary understanding of the plant of study. The heating and cooling needs are represented in this first step by simple heaters and coolers, without further integration. The purpose of this preliminary analysis is also to select different process architectures that will be optimized by the systematic optimization approach.
The variables that are imported from/exported to Aspen Plus and used for optimization, are introduced in Aspen Plus Calculator Blocks. The utility streams to be optimized, which are modeled in Aspen Plus, are set to zero. The property set HXDESIGN is selected to extract the information needed for the calculation of the film transfer coefficient of each individual stream (h, W m -2°C-1 ). The evaporator is divided into three bodies: economizer, evaporator and superheater, so as to take into account the different values of h (even if these functions will be finally provided by a single piece of equipment in the real plant design). The type of HEX considered for the calculation of h is a plate-type (planar, plate-and-frame) HEX, counter-current, with an average equivalent hydraulic diameter (L) of 0.006 m, defined as two times the gap between plates [25], and a fluid velocity (v ¥ ) through the channel of 6 m s -1 [26]. The expressions to calculate the h factor, considering turbulent flow (turbulent flow is taken into account, as it is beneficial for plate exchangers, from a convection heat transfer standpoint) can be found in [25]. The h values for water streams are estimated based on bibliographic data [27].

Multi-Period Process Design with Sequential Synthesis of the HEN
Given a plant layout, mainly three factors can determine the system performance and operational flexibility among the different periods: (i) operating variables, (ii) size of the productive components, and (iii) the multi-period HEN. As depicted in Figure 2, (1) the MOO optimization algorithm controls a sequence that consists in the resolution of the process system model in (2) Aspen Plus, followed by (3) the SYNHEAT algorithm (in AMPL), and (4) the computation of the weighted efficiency (over all the periods, in MATLAB) and the minimum number of heat exchangers for each multi-period configuration generated by the evolutionary optimization algorithm. The design variables of the candidate solution from MOO are first employed to simulate the process without a HEN in Aspen Plus to obtain the stream information for each period, which is employed in the lower-level mathematical programming. The multi-period SYNHEAT based HEN synthesis proposes HEN configurations and scheduling of HEXs at each period. Then, the objective functions are evaluated as the fitness of the candidate solutions. The evaluated solutions will be compared to be preserved or discarded. To keep the problem linear, (i) the approach considers the optimization of the design specifications at the QMOO upper level, while the optimal matching of streams is optimized at lower level using an AMPL-CPLEX [28,29] simultaneous solving approach, and (ii) HEXs areas are calculated in an subsequent step (see Section 2.3 for further explanation).
The ad hoc SYNHEAT superstructure for the optimization of the proposed CH2P plant has the following characteristics [30,24]: (i) constant h values for each stream, along the different periods (calculated from the steady state simulations); (ii) the minimum temperature (ΔT min ) is given as an input, thus, not considered for optimization; (iii) each HEX is assigned to the same pair of hot and cold streams at all the periods; (iv) countercurrent HEX configuration; (v) a fix number of stages established by the user (below the maximum number of either the hot and cold streams), as defined in the SYNHEAT superstructure [24]; (vi) for each stage, streams are split to cover all hot and cold streams matches; (vii) the outlet temperatures of each stage, for each stream, are treated as variables; (viii) in order to simplify the model formulation, utilities are placed at the outlet of the superstructure and streams are mixed isothermally; (ix) the approach from [31] is used to model condensers, evaporators and isothermal reactors; (x) consideration of one hot and one cold utility; (xi) the hot utility load is treated as a decision variable of the master (upper) level optimization, and therefore, it is treated as a process stream; (xii) the cold utility is considered a cold end, reducing the size of the combinatorial (superstructure) problem; (xiii) the HEN design can restrict stream splitting, can forbid matches and restrict/oblige other matches (e.g., avoid the contact of a combustible stream with an O 2 -rich stream); parallel and in-series configurations are allowed; (xiv) the dynamics when changing flow rates or temperatures between periods, e.g., residence times and thermal inertia, are neglected; (xv) the SYNHEAT approach does not deal with streams' temperatures change, but these are modified by the MOO algorithm.
Summarizing, the multi-period sequential synthesis of the HEN for this paper can be stated as follows: Given a set of hot and cold process streams with given period dependent (fixed or variable) mass flow rates and period dependent inlet and outlet temperatures (fixed or variable); and a set of utility options, with corresponding reference hot and cold thermal streams inlet and outlet temperatures; determine the optimal HEX configuration structure allowing to satisfy the thermal energy needs of the system for each period of operation; and the energy exchanged among streams and number of HEXs; with constraints on process flow intrinsic characteristics, stream matches, minimum temperature difference between hot and cold stream, heat transferred and stream splitting, to optimize the selected objective functions that best describe the user needs.

Problem Decomposition and Solution Methodology
The detailed sequential approach is presented as well in Figure 2. The connection among the different steps, i.e., the communication among software's and tools, calculation of objective functions and filtering (underlined in Figure 2) are all coded in an ad hoc set of MATLAB routines. In an evolutionary algorithm, the constraints (or conditions that must be fulfilled) are not enforced, but verified by filtering (see the last step before the loop). The solution procedure follows three big steps, in the light of decreasing the number of decision variables, and thus, of decreasing the computational time.
First step, individual period optimization. It considers the total amount of selected variables per iteration. The best performing variable's values are selected for each one of the periods by filtering the highest efficiencies. The selected constraints and gross power and H 2 needs (as defined by each period needs) are monitored.
Second step, bi-period optimization. The HEN proposed in the bi-period optimization only considers the two most relevant periods, according to the operating hours in the weighted efficiency. Only one variable is optimized in the current run: the hot utility flow rate and load, by mode. Therefore, a total of one variable multiplied by two period's variables are considered in each iteration. Once a NG flow rate is determined, the multi-period SYNHEAT algorithm calculates the amount of cooling water and the minimum number of stream connections to satisfy the heating/cooling needs. Analogously to the previous case, the best performing variable's values are selected for each one of the periods by filtering the highest efficiencies. Constraints, gross power and H 2 needs are monitored.
Third step, multi-period optimization. In this case, the proposed HEN meets all period's requirements. The followed approach is the same as that in the bi-period optimization, and thus, one variable multiplied by the number of periods equals the number of variables considered in each iteration. Again, constraints, gross power and H 2 needs are monitored.
By using this approach, the identification of a global optimum is not ensured. However, it allows to select, systematically and considering practical considerations that cannot be implemented in the model, a range of operating conditions and a HEN. It allows to quickly explore many process and HEN configurations with engineering and economic constraints provided by the CH2P consortium. Note that, for practical designs, it is not always the goal to find global optimality but improved (near-optimum) solutions.
The optimization procedure results in a HEN, based on a common set of connections for the selected periods, together with the corresponding set of operating conditions. The optimization result does not provide with common HEX's areas, to maintain the linearity of the mathematical problem. In the following step after optimization, the practical target goes towards the heuristic minimisation of the number of bypasses and parallel exchangers, to propose a common HEN with common UA factors for all the periods.
The calculation of the overall heat transfer coefficients (U, W m -2°C-1 ) accounts for the film transfer coefficient of each individual cold or hot stream (h, W m -2°C-1 ) evaluated at mean properties values. Film transfer coefficients are evaluated following the expressions and hypotheses from [25][26][27]32]. The use of UA factors, which stands for the overall heat transfer coefficient (U) multiplied by the exchange area (A ex , in m 2 ), instead of uniquely areas of exchange, aims at mitigating the uncertainty related to the estimation of heat transfer coefficients and thus, alleviating imprecision when using the calculated data for pilot plant HEN design. Note that the HeatX models in Aspen Plus consider user's A ex and U values.
The calculation of U proceeds, as shown in Eqs. (1) and (2): where, in the calculation of the corrected U factor, 0.95 corresponds to the temperature correction factor [25] and 1.1 accounts for a 10% of area increase due to fouling [32]. The wall resistance term d l in Eq. (1), plate thickness/thermal conductivity, is considered negligible [26].
The calculation of the logarithmic mean temperature difference (LMTD) [25] follows Eq. (3), to be used in the calculation of the area of exchange (Eq. (4)): where temperatures of the hot (T) and cold streams (t) are characterized at their inlet (in) and outlet (out) values (to/from the HEX). After the calculation of the HEX's areas, the HEN is modeled in Aspen Plus, and calculations are performed for every period. The final pilot plant configuration responds to subsequent layout adaptations, to conclude realistic and economic plant conditions.

Plant Description, Steady State Modeling and Selected Layout
The CH2P project aims at building up a prototype plant producing 20 kg of H 2 per day and 25 kW of gross power, at full capacity. This SOLIDpower ''large stack module concept'' consists of four stacks. The final step of the project is the conception of a full size plant, producing 400 kg of H 2 per day and 500 kW of gross power, based on, among others, the lessons learnt from the design and operation of the pilot plant. The plant has two well differentiated subsystems: (i) the hot one, formed by the reformer, SOFC, WGS reactors and water knockout unit, and (ii) the cold one, composed by the sulfur removal unit, water treatment unit (mainly a reverse osmosis process plus an electrodeionization module) and PSA, connected by the HEN. Three system architectures or layout configurations are preliminarily assessed, without HEN but with heaters and coolers that represented the needed temperature changes in the system, to determine the most suitable layout options to be considered for optimization.

Periods of Operation
The prototype is designed to follow the same typical requirements as a full plant in a HRS. The electricity produced by the system should be sufficient to meet the internal requirements, including self-electricity needs and HRS demand. Hydrogen production should mainly meet the necessities of the refueling station in terms of demand and storage capacity. Accordingly, the relevant operation points or periods for the prototype were defined. Herewith, the SOFC stack is intended to work at an extended range of FU, different NG flow rates (to the burner and to the SOFC), different SOFC stack air flow rates and air fed to the burner. This requires a wide dynamic range of operation for the CH2P plant, from full load (100% production of electricity and H 2 ) down to 15% of H 2 (minimum flow that the PSA unit can manage) and 17% of electricity (corresponding to an average value of the scaled down needs of the retail station that is present in the HRS). Table 1 and Figure 3 summarize the operation periods of the system. The nominal operation point, period 3, corresponds to 100% electricity production (25 kW) and 100% of H 2 production (20 kg day -1 ). If H 2 is the most important product of the plant, period 2 will be used. Otherwise, if electricity production is the priority, period 5 will be the selected operation point.

System Modeling
The feedstock preparation block consists of a sulfur removal unit, a water purification unit, a steam generator and a pre-reformer or steam reformer. Sulfur compounds, which are usually present in the NG distributing pipelines, can be harmful for the pre-reformer catalyst and the Ni-YSZ cermet anode in the SOFC. Hence, a sulfur removal unit is required to remove sulfur additives from the NG stream. The allowable concentration limit for the sulfur is less than 1 ppm [33]. The water purification unit demineralizes and cleans the water stream from biological impurities. The water/steam flow rate required for the system is determined by the steam to carbon ratio (S/C) which is required after mixing steam with the NG stream. The S/C ratio is a critical parameter to prevent carbon deposition in the SOFC system and stack module. The value of S/C ratio depends on the composition of the NG mixture supplied to the system. Typical values for the S/C ratio are in the range of 1.5 to 2.5 [34,35]. For the current study an S/C ratio of 2.2 is used. The degree of superheating of the steam is determined so as to prevent condensation when steam is mixed with the cold NG stream. The mixture of steam and cleaned NG is fed to the steam reformer. An isothermal and non-adiabatic reforming unit is envisaged. The conversion ratio of the methane and higher hydrocarbons present in the NG stream is determined by the operation temperature of the reformer. The operation temperature of the reformer is a variable that is controlled by the thermal balance of the system during the MOO. For the system under consideration, a conversion ratio in the range of 40 to 90% is expected.
The reformate mixture is then heated up to the SOFC stack inlet temperature and supplied to it. Air to the SOFC stack enters the system from ambient conditions, to be then heated up as well to the SOFC stack inlet temperature and introduced into the stack. The fuel and air inlet temperatures do not need to be the same, and can vary depending on the operation points and thermal energy balance of each period. Additional NG feed is provided to the burner for thermal management purposes depending on the operation points.
The FU and mass flow rate of fuel are determined by the H 2 and power requirements for each operation period. Air mass flow rate is determined by two factors; the temperature difference between the SOFC inlet and outlet and/or the oxygen utilization ratio in the SOFC. The anode and cathode outlet streams are then fed to the downstream (SOFC) process. Downstream the SOFC, the anode off-gas is further purified to recover the H 2 fraction, in a system formed by two WGS reac-  tors (at high and low temperature -HT, LT) and a PSA. The WGS reactors convert CO into CO 2 , aiming at increasing the H 2 content and reducing the CO concentration, which is a critical impurity that affects the performance of the PSA downstream unit; the lower the CO concentration, the higher the yield of the PSA unit (i.e., more H 2 is recovered). Cathode offgas can be cooled down and discharged, or further used in the system burner. The system architecture varies based on the integration of upstream and downstream processes; i.e. the utilization of cathode and anode outlet streams. A preliminary analysis, before optimization, of three PFD options is performed without a HEN and without optimization of the split ratios for the cathode split (in Options 1 and 2) or reformate split (in Options 1 and 3), in Section 2.3. Find below the model description of each one of the process components. Specific units' heat losses follow SOLIDpower experience, and are assumed to be independent of temperature.
The Aspen Plus plant model uses the Peng-Robinson Boston-Mathias property method, which applies the Peng Robinson cubic equation of state with the Boston-Mathias alpha function for all thermodynamic properties. Apart from the S/C ratio, the operating temperature and the area specific resistance (ASR) of the SOFC stack are also fixed, at 740°C and 0.47 O cm 2 , respectively. The value of the ASR, which is indicative of the SOFC performance, is obtained from experimental analysis of the SOLIDpower SOFC stack (see further explanation below). The system works at atmospheric pressure, except the PSA unit and the final H 2 conditioning (pressurization) step. Pressure losses are neglected. Dutch NG composition (that is 0.814 of CH 4 , 0.144 of N 2 , 0.03 of C 2 H 6 , 0.01 of CO 2 and 0.002 of C 3 H 8 , in mole fraction) is considered as reference fuel composition for the CH2P plant.

SOFC Model
A lumped and non-isothermal SOFC model is utilized. The model is based on a commercially available, anode-supported cell design SOFC stack. The following assumptions are made for model implementation: (i) the electrochemical and chemical reactions occur at an average reactor temperature; this assumption, though not entirely accurate, is acceptable for a 0-D SOFC stack model in process system simulations [36][37][38] (ii) the reactions reach chemical equilibrium at that temperature, and the outlet gas composition is at equilibrium; and (iii) the temperature of the outlet anode and cathode streams are equal.
The global mass and energy balances are solved to obtain the outlet temperature of the anode and cathode streams. The thermodynamic ideal voltage is calculated following the methodology presented in [39]. The electrochemical model for the SOFC indicating the stack performance is lumped into an ASR parameter. The ASR value encompasses all the loss mechanisms such as the ohmic losses due to charge transport, activation losses and losses related to the mass diffusion through the electrode. For the current study, the ASR parameter is obtained from the performance data provided by SOLIDpower, in the framework of CH2P project, from the experimental investigation of their anode supported SOFC stacks. The performance data were obtained for the SOFC stack operating at its nominal average stack temperature of 740°C, and for varying fuel flow rates and utilization. It is assumed that heat loss of the SOFC stack is 2.5 kW.

Pre-reformer
A steam reforming process is considered. The inlet amount of NG depends mainly on the amount of electricity and H 2 to produce or the selected FU. The steam fed to the pre-reformer is determined by the fixed S/C ratio. The NG, depending on the quality and source, contains low concentrations of the higher hydrocarbons (such as butane or pentane) apart from methane. The generalized steam reforming reaction is given by equation Eq. (9): The pre-reformer is an isothermal reactor, modeled under the assumption that it reaches equilibrium compositions. The thermal energy required for the reforming reaction and operation temperature of the reactor (and thus, conversion ratio of methane) is determined by the thermal management of the SOFC stack. The pre-reformer is implemented in Aspen Plus using the RGibbs reactor model available. A constant heat loss of 0.25 kW is assumed.

WGS Reactors
The HT and LT WGS reactors are modeled as a non-isothermal reactors, and only reaction from Eq. (7) takes place. The outlet composition and temperature are obtained based on the assumption that equilibrium is reached at the outlet. The reactor is simulated in Aspen Plus through an REquil reactor block available in the standard library. A constant heat loss of 0.1 kW is assumed for the model.

Pressure Swing Absorption (PSA) Unit
The PSA unit separates the H 2 from the anode off-gas mixture which also contains CO 2 , N 2 , unconverted CH 4 and a very low fraction (less than 1 mol%) of CO, to produce an H 2 stream that complies with ISO 14687-2 (fuel cell quality; CO concentration below 200 ppb). The PSA consists of 4 columns that go through cycles of adsorption, desorption, pressurization and depressurization. A detailed explanation about the PSA functioning can be found in [40,41]. The PSA operates at 7 bar(a). The gas stream after the WGS process and water knockout (i.e. water separation and recycling), is compressed up to 7 bar(a) in a 4-stage compressor with interstage cooling, as the up to 30 bar(a) compressor. The water vapor in the gas stream that condenses during compression is removed. The dry compressed gas is fed into the PSA unit. The PSA performance can be characterized in terms of the H 2 separation ratio, defined as the ratio of moles of H 2 separated from the feed stream to the total moles of H 2 in the feed stream (between 0.8-0.9 in the current work). The PSA is modeled using a separator block (Sep) from Aspen Plus library. The separation ratio is provided as an input to the model.

Burner
The burner combusts the off-gas from the PSA, a make-up NG stream, the anode off-gas split and H 2 (provided the PFD option). Cathode outlet from the SOFC may be fed to the burner as the oxidizer (depending on the PFD option), and fresh air is also supplied based on the combustion requirements and temperature limits. The combustion chamber is modeled in Aspen Plus using an RStoic reactor model, which considers all the possible combustion reactions. A constant heat loss value of 0.6 kW is assumed for this unit.

Water Knockout Unit
After the WGS reactors, the gas stream is further cooled down to a temperature below 100°C to allow water condensation. In this unit, the condensed water is removed and recirculated to the inlet. It is simulated as a flash separator (Flash2 unit) from Aspen Plus library.

Sulfur Removal Unit, Pumps, Compressors and Blowers
The sulfur removal unit is modeled as a separator block (Sep), which splits all sulfur compounds. The auxiliary components such as pumps and blowers are implemented using the standard models available in the Aspen Plus library. The mechanical and polytropic efficiencies for the pumps, compressors and blowers are provided as user inputs, with values between 0.7 and 0.9 [25].

Preliminary Analysis of Three PFD Options
This preliminary analysis has as purpose to select the process layout or layouts which have a larger efficiency and performance, to be further optimized in the next Section. The result is based on the energy efficiency.
In PFD Option 1 (Figure 4), the dry anode off-gas stream is split into two streams, one stream is fed to the PSA to produce H 2 , while the other stream is fed to the burner for thermal management. The ratio of this split is determined by the operation mode and thermal management needs. The cathode outlet stream is also split into two streams, one of which is sent to the burner and the other is released to the atmosphere after it is cooled down in the heat recovery units.
In PFD Option 2 ( Figure 5), the dry anode off-gas stream is entirely sent to the PSA unit. The H 2 produced from the PSA is split into a product stream and a minor stream, which enters into the burner. The split ratio is determined by the H 2 production volume and thermal management needs. The cathode outlet stream downstream process is similar to PFD Option 1.
In PFD Option 3 (Figure 6), the anode outlet stream follows the same path as in Option 1, whereas the complete cathode outlet stream is released after cooling in the heat recovery unit. This configuration provides an active pressure control on the cathode side to regulate the differential pressure between the anode and cathode chambers in the stack.
A summary of achievable efficiencies for the three PFD options are presented in Figure 7, for the different operation periods (as summarized in Table 1 and Figure 3). Note that these results do not come from a process optimization, but from discrete variables adjustment performed to fulfill the outlet needs. The highest performance is achieved by Option 3, followed by Options 1 and 2, in all the periods. The efficiency values for the three options are pretty close to each other; however, Option 2 has the lowest value due to the conversion of the anode off-gas into a fraction of H 2 that is burnt downstream.
This preliminary analysis is concluded in the PFD represented in Figure 8, which represents two design options, i.e., two cathode off-gas possibilities, to be further studied in Section 4. In Option A, the complete cathode outlet stream is released after thermal energy recovery. In Option B, the cathode outlet stream is sent to the burner to combust the off-gas stream from the PSA unit together with an additional NG stream supplied to the burner. Option A has an additional heat source, Q9, at a temperature above 700°C. Streams Q5 and Q6 deal with high temperatures as well, above 700°C, being placed at the outlet of the SOFC stack.

Multi-Period and Multi-Objective HEN Optimization
The described steady state model is adapted to (i) reflect the requirements of the calculation of the film transfer coefficients and of the SYNHEAT superstructure algorithm, (ii) automatize the iterations and (iii) accomplish the requirements of the Aspen Plus handler in the MATLAB routines described in Section 2. In the SYNHEAT superstructure, the reformer (C102, as named in Figure 8) and evaporator (Q2EV) streams are considered as isothermal. The NG stream fed into the burner (as the hot utility of the system) is set at zero, as it is calculated during the optimization.
The calculations performed using the described approach in Section 2, are summarized in the following Sections 4.1 and

Input Data for HEN Optimization
The following input conditions are implemented in the SYNHEAT superstructure algorithm (unit's nomenclature refers to Figure 8): (i) All five operation periods are consid-ered; (ii) Q13 and the pre-reformer (C102) are physically the same unit, and thus they must be heated up by the same stream; (iii) the water knock-out process occurs at 35°C; i.e., T out of Q8 is set at 35°C; (iv) T out of Q5 and Q9 are fixed at 160°C (i.e. temperature of vented gases); (v) T out of Q6 (HT-WGS) is 325°C and of Q7 (LT-WGS) is 200°C; (vi) there are five hot streams: Q5, Q6, Q7, Q8, and Q9, and seven cold streams: Q1, Q4, Q13, Q2SH, Q2E, Q2EV (note that the evaporator is divided three bodies, economizer, evaporator and superheater), C102; (vii) four temperature stages (as defined by the SYNHEAT superstructure); (viii) Q2EV and C102 are streams at constant temperature; (ix) forbidden matches among the following streams: Q6-Q1, Q7-Q1, Q8-Q1, Q9-Q13, Q9-C102, to avoid contact between a combustible gas and oxygen rich streams; (x) no splitting of high temperature streams (above 500°C), so no parallel heat exchange is allowed in the streams: Q5, Q6, Q9, Q4, and C102, to practically avoid the use of high temperature valves; (xi) constant heat losses in SOFC stack, reformer, burner and WGS reactors are considered as reported in Section 3; moreover, the air recuperator, Q1, has a constant heat loss of 0.25 kW; (xii) the stack inlet cathodic line suffers from a temperature loss of 10°C; (xiii) the cold utility circuit for the cooling of the compression stages is not modeled (cooling is included though within the MCompr model from Aspen Plus); (xiv) ΔT min of 50°C for the hot and cold streams exchange, ΔT min of 100°C for the hot and cold streams exchange at temperatures above 500°C (Q4, C102, Q5, Q6, Q9)  and ΔT min of 5°C for the hot stream-cold utility exchange; (xv) for the calculation of a weighted efficiency, the time of operation of periods 2 and 3 is 9 h day -1 for each one of them; of periods 1, 4 and 5 it is 2 h day -1 . For bi-period optimization, it is assumed that periods 2 and 3 (the two most used periods) operate each during 12 h day -1 .
Following the methodology explained in Section 2, the value of the film transfer coefficients h for all the streams of the CH2P plant are summarized in Table 2. This estimation, as cal-culated before the optimization procedure, takes into account the steady state conditions derived from the modeling work in Section 2, for the five operation periods. Average values are summarized in Table 2 (variation among periods of +/-7%).
The h values are kept constant along the HEN design process. Note that an exhaustive calculation of temperature dependent h is not within the purposes of the current work. The h values for water streams, Q2E, Q2EV, Q8 and cooling water (cold utility), are estimated based on bibliographic data [27], based on shell and tube geometry: 4,000, 6,000, 5,000, and 4,000 W m -2 C -1 , respectively.
The composite curve for the CH2P plant, period 3 and Option B, before heat integration and with a unique ΔT min value for all the streams (15°C), is shown in Figure 9 as example. It shows that the case study is in essence a threshold problem requiring only a cold utility, as the hot needs are provided by the NG sent to the burner inside the plant layout. The cold composite curve shows the two isothermal processes: reforming and evaporation. The hot composite curve is characterized by the high temperature of the flue gas, which is the main stream that provides heat to all the system. The cold utility is needed at low temperature, more specifically, for water condensation. One cold utility was therefore selected: water that can be used in the retail station as hot water (with a temperature of 65°C; heated up from 20°C), being part of the cogeneration solution that the CH2P plant introduces in a HRS.

Decision Variables and Constraints
The decision variables used in the current MOO problem are of two types: (i) binary variables, as for the two process layouts considered (Option A and B); this consideration results in two sets of optimization results; and (ii) continuous variables. Eight independent variables are used, and 6 constraints (solutions filtering) are taken into account, so as to respect most relevant process conditions. The list of independent variables and constraints are provided in Table 3 and  Table 4, respectively. The different values ranges have been set based on the operating experience, needs and recommendations from HyGear and SOLIDpower. The ΔT in the stack (as the difference between the inlet and the outlet temperature of the cathodic air) is also calculated and monitored, aiming for minimum values.
The optimization problem has in total 8 decision variables. A first optimization attempt, aimed at optimizing the complete set of 8 decision variables, for the 5 periods at the same time, yielding a 40 variables problem (see Figure 2). The solution procedure runtime was too large; that is the reason why the solution procedure follows the three steps described in Section 2.3, i.e., the sequence of individual optimization, bi-period optimization and multi-period optimization.

Objective Functions
The performance evaluation for HEN synthesis counts with two objectives: (i) Maximization of the weighted efficiency. It takes into account all the single periods' efficiencies, as in Eq. (10), multiplied by the average number of working hours per period per day (as described in Section 4.1), and divided by 24 h. Hydrogen and electricity (net SOFC power, i.e., gross power minus plant electricity needs and H 2 compression up to 30 bar(a)) are the considered products in the numerator.
where _ m is the mass flow rates in kg s -1 , LHV is the lower heating value in kJ kg -1 . Net electricity is expressed in kW. In the following Results section, the weighted efficiency values reported after optimization also take into account heat recovery as product in the numerator. (ii) Minimization of the number of heat exchanger matches. This is a discrete objective. As such, only integer numbers can be considered, and the range of variability is small, as the same heat-cooling needs have to be fulfilled.

Performance and Selection of System Layout -Common Set of Connections
The process system layouts for the Options A and B, with the optimized HEN based on a comment set of HEX connections for MOO and multi-period optimization, are depicted in Figures 10 and 11, respectively. The ''QX'' nomenclature corresponds to the HEX nomenclature from Figure 8. Note that one cold or hot stream can exchange heat with more than one hot or cold stream; for instance, the hot stream represented by Q5 exchange heat with streams Q4, Q1, C102, Q13 and cooling water in Figure 10. Connections (stream splits and mixes) that are active only during specific operation points are represented by dashed lines in the figures. The resulting HEN for Option A has 13 HEXs, while that for Option B has 12 HEXs.   The numeric difference corresponds to the use of a single or multi-body Q1. Moreover, steam generation happens differently in both schemes (Q2 layout), and specific periods operating needs are also different (only Q2 in Option A; Q2, Q5 and Q6 in Option B). The resulting values for the independent process variables for Options A and B at different operation periods are provided in Tables 5 and 6, respectively. The reported values for variables from 1 to 7 corresponds to the single period optimization step. The reported value for variable 8 correlates with the NG flow rate obtained in the last multi-period optimization step. Table 7 summarizes the efficiency values for both options at every step of the optimization. It points out that a common HEN induces higher use of NG as hot utility, exchanging heat at the hottest stage, triggering lower efficiency values, if the values for individual optimization are compared with those of multi-period optimization. Note that the weighted efficiency can be penalized by the low efficiency of periods that are less used, resulting in larger values for the bi-period optimizations.
The multi-period weighted efficiencies of the two layout options for the CH2P prototype system are summarized in Figure 12 for all the operation periods. It can be observed that the overall system efficiency for Option B is higher than for Option A. The weighted system efficiency for Option B is 9 percentage points higher than that for Option A. Option B achieved a weighted system efficiency of 59.3% compared to 50.9% for Option A. The higher performance of Option B compared to Option A is mainly due to the effect of cathode offgas supplied to the burner. In Option A, air at ambient temperature is supplied to the burner, lowering the specific enthalpy of the flue gas from the burner and the flue gas temperature. Lower flow rate of NG to the burner is required in Option B than in Option A in order to meet the thermal management  The selected layout, Option B, is further studied towards the overall efficiency. The sensitivity analysis from Figure 13, formed by the different points of study considered in the optimization, shows different operating zones for different H 2 production targets. The results from Figure 13 represent the efficiency against the ratio of H 2 calorific value to the gross electrical power produced. The colored scheme describes the amount of H 2 produced. As compiled in Table 8, each point represented in the graph, corresponds to a set of independent variables for the cogeneration plant, with different outputs of H 2 and gross power. Points 1, 2, 3, and 4 correspond to the Pareto front of sets of working conditions that maximizes efficiency and H 2 produced. Table 8 points out that higher FU's (thus, larger production of power -Points 1 and 2), perform the highest system's efficiencies with a common set of heat exchange connections. Figure 13 shows that for the integrated system, higher efficiencies (above 0.6) are obtained for a ratio

Common HEN with Common UA Factors
The selected multi-period Option B is further adapted with stream's splitters and mixers to conveniently accommodate the different heating/cooling needs. The steps following the determination of a common set of HEN connections and driven to find a common set of HEX for all 5 periods, are: (i) the selection of the largest UA factors per HEX among the results obtained for each period after MOO and multi-period optimization (see the selected UA factors in Table 9); (ii) the elimination of the HEXs that are bypassed in any of the periods; (iii) the change of intermediate temperatures (stage temperatures), while keeping set point temperatures, and the adjustment of the variables burner outlet temperature and NG flow rate that goes to the burner (and the consequent stream of air going to the burner), i.e., independent variable number 6 and 8 in Table 3.
The considered constraints are summarized in Table 10, as optimization is not performed in this step (as mentioned before, the different values ranges have been set based on the operating experience, needs and recommendations from HyGear and SOLIDpower). Note that the non-systematic adaptation of the proposed HEN connections requires the relaxation of some temperature values that were considered as fixed during the previous step. These values are listed in Table 10 as constraints and correspond to the inlet temperatures of the PSA, water treatment plant and WGS reactors, as well as the outlet temperatures of flue gases and superheated steam (constraints 11 to 16). The non-systematic adaptation of the proposed HEN connections requires as well the conversion of some independent variables (Table 3) into dependent variables; these are the SOFC inlet fuel and air temperatures and the temperature of pre-reforming (variables 1, 2, and 4). Overall, the practical target goes towards the heuristic minimisation of the number of bypasses and parallel heat exchangers.
The operating conditions and efficiencies for each period are listed in Table 11. Figure 14 shows the proposed solution, with 8 HEX. It can be observed  Sensitivity analysis (search space used in the optimization) based on independent variables ranges (except for FU, which has been varied between 0.01 and 0.85, for illustrative purposes), including S/C. H 2 power and SOFC power are not subjected to specific operating modes. Color scheme refers to the amount of H 2 produced / kg s -1 . Table 8 Set of decision variables (including S/C) corresponding to the points highlighted in Figure 13.  Table 11 are higher than the efficiencies reported in Table 7 (for Option B) mainly due to the relaxation of different temperature values that were considered fixed before, and due to the layout changes. Moreover, higher efficiencies are obtained in modes where more H 2 is produced. The CH2P plant layout weighted efficiency is of 63%, while the efficiency considering hot water generation (thus, the tri-generation of H 2 , electricity and heat) goes up to 75%. Period 3 is the one with the highest efficiencies (as it was the one used as reference along the whole design process); i.e., 67% and 81%. Efficiencies taking into account tri-generation, range between 61-81%.
In the light of the presented results, the system design is further adapted, taking into account engineering and economic constraints. Further modifications derived into the description of the operating envelope of the system, understood as the operating limits of the plant, so as to specifically know the capabilities of the system, beyond the 5 operation periods (see Figure 15). In order to know this working area, several operation periods are calculated (periods from 6 to 12). It is considered that the minimum power that the plant has to provide must cover the electricity needed to fulfil (i) H 2 compression from 1 up to 200 bar(a) (assuming hypothetical storage needs at that pressure), (ii) retail station, and (iii) HRS needs (items (i), (ii), and (iii) correspond to the three distinguished sections of Frontier 1 or F1 in Figure 15). The limit of this F1 gray zone represents the operation points that produce zero net power. The minimum power that the CH2P plant is able to produce depends on the amount of NG that is sent to the SOFC and the FU used. This value is around 8 kW. The Table 9 Calculated UA factors for each on the hot-cold matches resulting in the final HEN proposed in Figure 14.   minimum H 2 produced corresponds to the 15% of PSA maximum load. It is also limited by the maximum FU that the SOFC stack can handle (two distinguished sections in Frontier 2 or F2 in Figure 15). No plant operation is expected in the gray zones from Figure 15. The zone between gray areas sum-marizes all the operation points where the CH2P prototype could work, with the system configuration depicted in Figure 14, and operating conditions summarized in Table 11. Trigeneration efficiencies around and above 70% are expected.

Conclusions
A methodology for the design of a HEN has been implemented within the framework of the EU H2020 CH2P project. The novelty of the current work lies in the application of a multi-objective and multi-period optimization approach for the conceptual design of a SOFC-based electric vehicle station, to allow the feasible and efficient operation of the system under multiple demand situations. It required extending the capability of the software tools available at the start of the project to include in the workflow methods recently developed for multi-period optimization. The interface functionalities required specifically for the CH2P project were coded in a separate and dedicated standalone set of MATLAB routines. The present study is an illustration of the capability of the developed tool that can be adapted and used for the HEN conceptual design of other systems, for which the principles of the SYNHEAT algorithm are relevant. The sequence in the identification of viable solutions started with HEN topologies pro-  posed by the MOO. Several layouts and unit's assessments were then made to achieve a single set of HEX's, including practical considerations that cannot be implemented in a MILP model for MOO. The preliminary outcome for the CH2P prototype is a configuration comprising 8 planar HEX (characterized by their UA factor), yielding period efficiencies (considering H 2 compression up to 30 bar(a)) that range from 53% (period 1) to 67% (periods 3 and 4). The cold utility considered is cooling water heated up to 65°C, which can be used in the HRS as a cogeneration option. Tri-generation efficiencies range between 61% (period 1) to 81% (period 3). The current paper concludes the conceptual design of the CH2P plant; the next steps are directed towards practical decisions towards real implementation of the plant design. for the HEN drawing. CH2P project has received funding from the Fuel Cells and Hydrogen 2 Joint Undertaking under grant agreement number 735692. This Joint Undertaking receives support from the European Union's Horizon 2020 research and innovation program, from Hydrogen Europe and N.ERGHY. Swiss partners has received funding from the State Secretariat of Education, Research and Innovation (SEFRI), under contract number 16.0223.