Modeling and Multi-Criteria Optimization of a Process for H₂O₂ Electrosynthesis

Abstract

This article introduces a novel laboratory-scale process for the electrochemical synthesis of hydrogen peroxide (H${}_2$O${}_2$). The process aims at an energy-efficient, decentralized production, and a mathematical optimization of it is presented. A dynamic, zero-dimensional mathematical model of the reactor is set up in Aspen custom modeler^®. The proposed model constitutes a reasonable compromise between complexity and convergence. After thoroughly determining the reaction kinetics by adjustment to experimental data, the reactor unit is embedded in an Aspen Plus^® flowsheet in order to investigate its interaction with other unit operations. The downstream contains another custom module for membrane distillation. Electricity appears as a resource in the process, and optimization shows that it reaches product purities of up to 3 wt.-%. Both the process optimization and the adjustment of the reaction kinetics are treated as multi-criteria optimization (MCO) problems.

1. Introduction

Electrochemical applications have undergone a remarkable increase in importance during the last decades due to the need for efficient modern technologies, e.g., in syntheses of substances or mobility, but also due to the rapidly evolving availability of new experimental and numerical tools for precise study and characterization of electrochemical systems [1]. One example of such an electrochemical device is the fuel cell (FC) for generation of electrical energy from oxidation of hydrogen or short-chained hydrocarbons like ethanol or methanol [2,3]. The most promising candidate for the large-scale production of hydrogen for use in fuel cells is the electrolyzer for water electrolysis [4], which is itself based on electrochemistry. Moreover, batteries are currently facing strong interest from researchers, since they are an essential component, e.g., in modern microelectronic devices, in electromobility and in energy storage. Latz and Zausch gave thermodynamically consistent theoretical descriptions of transport phenomena [5], reaction kinetics [6] and thermal aspects [7] in these devices. Another application of electrolysis is the conversion of carbon dioxide (CO${}_2$) [8] to valuable products such as organic acids or short- and long-chained alcohols in order to reduce the global carbon footprint. A work realizing such a synthesis route experimentally to produce ethylene using copper(I) oxide catalysts can be found in [9].

The process introduced in this work comprises the oxygen reduction reaction (ORR) mechanism, with hydrogen peroxide (H${}_2$O${}_2$) as the desired product. This important commodity chemical has a wide range of applicability in decentralized or industrial processes. In the chemical industry, it is synthesized via the anthraquinone process, delivering product purities of 1–2 wt%, which can be enhanced up to 70 wt% by downstream processing [10]. Due to its high oxidation potential, it is commonly used for treatment of waste water or bleaching of paper pulp or textiles [11,12]. Furthermore, it is employed extensively in the production of fine chemicals [13]. In the encymatic halogenation of the phenolic monoterpenes thymol and carvacrol to different products for anti-inflammatory, anti-microbial and anti-cancer applications, H${}_2$O${}_2$ accompanied by the enzyme chloroperoxidase is used as a bio-catalyst [14]. Getrey et al. showed that the highest conversion rates are reached with electrochemical in-situ synthesis of H${}_2$O${}_2$ directly in the active zone for the following reaction [15]. Thus, there are no transport problems for the catalyst as in the case of external injection of hydrogen peroxide. Since oxygen is used here as a reactant to produce H${}_2$O${}_2$, maximum efficiencies are achieved using gas diffusion electrodes [16] to avoid the problem of limited solubility of O${}_2$ in the electrolyte.

Decentralized electrochemical production of H${}_2$O${}_2$ has received growing interest from researchers in recent years. Yang et al. perceive this as part of a general trend [17]. They predict that the increased availability and decreased cost of renewable electricity will transform the chemical industry towards on-site production. This could lead to logistic challenges, e.g., if larger amounts of product have to be gathered. On the other hand, the decentralized approach has the potential advantages of more efficient energy conversion, easier product storage and lower capital expenditure according to Sehrish et al. [18]. In this respect, Li et al. recently stated that gas diffusion electrodes, such as the one used in the current work (see Section 2.3), are crucial to exploit the advantages of an electrochemical production of H${}_2$O${}_2$ [19]. They quantitatively assessed their reactor with respect to energy consumption, production capacity, cost and long-term stability. They conclude that the technology is already feasible for on-site abatement of pharmaceutical residues in ground water. As for long-term stability, they report a life time of their electrode limited to about 46 days due to flooding of pores. In a real-world application, the electrodes need to be replaced before expiration to ensure a site-independent product quality.

Inspiration for modeling electrochemical reactors can be obtained from theoretical works on fuel cell research, the quantity of which is higher than that for electrolyzers. Since FCs convert substances into electric power, they can be considered as the inverse operation of electrosynthesis. Schultz et al. set up a rigorous, one-dimensional model of a direct methanol fuel cell using a Flory-Huggins activity model [20]. Later, the same group investigated the dynamic response of the system with respect to reaction kinetics [21]. The basic features to be covered by a model in any case are mass balances, material conversion, potential drops and membrane flows. In the case of spatial resolution (dim. > 0), local transport effects need to be described, too. They may be Knudsen diffusion or Maxwell–Stefan diffusion, the latter being governed by molecular friction in porous media [22], which is taken into account in [20]. Furthermore, Weinzierl and Krewer used numerical simulations to analyze the water management in direct methanol fuel cells [23,24]. For a more sophisticated approach, see e.g., [25], where a version of the 2D Navier–Stokes equations is solved for a mixture in order to investigate its coupling to chemical surface reactions via boundary conditions.

A mathematical model of an electrochemical cell for copper electrolysis, composed of 0D and 2D submodels, was given by Pohjoranta et al. [26], who validated their theoretical predictions by comparison to cyclic voltammetry measurements. However, this work is concerned with the overall process. Details about unit operations are secondary as long as their main dependencies are captured with reasonable accuracy. Hence, we stick to a 0D reactor model such as the one of Görgün et al. [27], who published a balance-based mathematical description of an electrolyzer with a proton exchange membrane (PEM), or the one of Argyropoulos et al. [28], who set up a semi-empirical 0D model of a direct methanol FC in a similar fashion. The aspired philosophy of the model equations is to aim for fast convergence while realistically reproducing measured data for the most important output quantities, i.e., the product weight fraction and the Faraday efficiency. This is additionally beneficial, because it does not imply a high number of transport coefficients as in the case of spatial resolution. Such coefficients (e.g., the binary Maxwell–Stefan diffusional coefficients) are often unknown in praxi due to the limited availability of experimental data for parameter estimation. However, it will be necessary to thoroughly determine the kinetics of the reaction mechanism under consideration.

A means of studying electrochemical reaction kinetics, which relies on experimental data, is the rotating ring-disk electrode (RRDE) system [29,30]. The method is based on the rotation of circular electrodes, inducing a forced flow to drag reactants to the electrode. Then, the detected electric currents at the disk and at the ring, with the latter being located at higher radii, contain information about the reaction kinetics. Yet, a micro-scale simulation is required in order to extract this information from experimental results. An application to determine electro-kinetic rate constants of the ORR mechanism at gold catalysts can be found in [31]. Recent reviews about catalysts for electrosynthesis of H${}_2$O${}_2$ by oxidation of water can be found in [32,33]. Recently, two-dimensional RRDE simulations based on a finite-difference discretization were successfully adapted to corresponding measurements in a detailed and quantitative electro-kinetic study [34]. Such a way of determining parameters of a unit operation by simulations on a smaller scale can be considered as a part of a bottom-up modeling approach, as sketched in Figure 1. Alternatively, the kinetic parameters can be calculated by direct adjustment of the unit model to experimental data, as in this work.

The process model proposed in this article is set up in Aspen Plus V10^®. However, this commercial flowsheet simulator does not provide a model of an electrochemical reactor as a unit operation and hence, an electrolysis cell model was implemented in Aspen custom modeler V10^® (ACM) here. Note that this has been done before by Redissi and Bouallou [35], who performed an economic evaluation of a high-temperature co-electrolysis process for conversion of CO${}_2$/H${}_2$O to synthesis gas, based on simulations in Aspen Plus^®. Apart from that, the process comprises a membrane distillation module [36] for separation of acidic components from the product solution in the downstream processing, which is neither available in Aspen Plus^®. Implementations of such unit operations in Aspen Plus^® can be found in [37,38]. A critical review on the modeling of membrane distillation technologies was given by Hitsov et al. [39]. A species-selective flow through a membrane is forced by the thermal gradient between a heated feed side and a cooler permeate side. The pores are hydrophobic and therefore the liquid phases from each side only partially penetrate the membrane. Material transport through the gas-filled segments of the pores is initiated by the difference in vapor pressures resulting from the different temperatures on both sides.

In most practical applications, electrochemical modules are part of higher-level processes with potential up- and downstream processing. However, the availability of studies considering electrochemistry on a process scale, i.e., in a coupled network of interacting modules rather than on the level of a unit operation, is limited. One of the few works on this topic was published by Karst et al. [40], who employed a multi-scale approach to design a fuel cell process with methane reforming in the upstream part for automotive application. Furthermore, Sanchéz et al. modeled an electrolysis process for hydrogen production in Aspen Plus^® [41]. The current article contributes to this point by applying multi-criteria optimization (MCO) [42,43] with an adaptive scalarization scheme to a process containing an electrochemical reactor. The resulting optimal operating points apply for the entire process rather than a single device. For the MCO, Aspen Plus V10^® is controlled via a programming interface. For this purpose, it offers an ActiveX automation server to be used within a VBA (visual basic for applications) code or a VB (visual basic) code. Previous work exploiting this interface for the optimization of separation processes can be found in [44,45]. The VBA code from [42] is used in the current work, since it already comes with an automated graphical user interface in MS Excel. Only the flowsheet simulator is exchanged here.

The model equations for the electrochemical reactor and the membrane distillation module are discussed in detail in Section 2.1 and Section 2.2 of this article. Subsequently, in Section 3.1, particular attention is paid to the determination of electro-kinetic rate constants for the reactor model. This is realized here by direct adjustment of the model to experimental data, using MS Excel as a central interface coupling the simulation with an external optimization solver. Since two measured variables are available, the adjustment procedure is treated as an MCO problem, as suggested in [46]. A validation of the parametrized model equations is presented, too. Finally, according to the second modeling step in Figure 1, the adjusted model is embedded in a flowsheet for the entire process in Section 3.2. Pareto-optimal working points are then adaptively calculated [43,47], with particular respect to the electrical power required per product amount and the heat input.

2. Materials and Methods

2.1. Electrochemical Reactor Model

The electrochemical cell consists of a proton-exchange membrane (PEM), separating the two half cells, and two electrodes, each of which is separated from the PEM by a thin gap filled with electrolyte. Since the process incorporates a weakly soluble gas (oxygen in this case) as one of the educts, the reactor is designed with a gas diffusion electrode as cathode adjacent to a gas compartment. The model equations given in this section describe this as a zero-dimensional set-up. Yet, they contain several geometry parameters. The model captures time-dependent behavior, so that it can be applied in dynamic and stationary flowsheet simulations. The employed species are water, hydrogen peroxide, oxygen and hydrogen. In addition to that, sulfuric acid (H${}_2$SO${}_4$) as solute in the water is considered, to become dissociated with the ions H${}_3$O${}^{+}$, HSO${}_4^{-}$ and SO${}_4^{2-}$. The models are implemented using Aspen custom modeler V10^® with ENRTL-RK as the employed set of thermodynamic methods for the electrolyte.

2.1.1. Molar Balances in the Cathode Half Cell

The amounts of substances $n_{c,i}$ of species i in the cathode half cell are determined by the first-order differential equation

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}{n}_{c,i}}{\mathrm{d}t}& =& F_{cin}{x}_{cin,i}-F_{cout}{x}_{c,i}+F_{gin}{x}_{gin,i}\hfill \\ & -& F_{gout}{x}_{g,i}+N_{m,i}+\sum _j{\nu }_{c,ij}{v}_{c,j}(i\ne H_{2}O),\hfill \end{array}$$

(1)

where the last term corresponds to reactive conversions of substances. For an explanation of each quantity, refer to the list of symbols in Appendix A. The mole fractions in the cathode half cell and the gas compartment are related to the amounts of substances by

$$\begin{array}{ccc}\hfill x_{c,i}& =& n_{c,i}/\sum _{i=1}^N{n}_{c,i},x_{g,i}=0(i\in liq),\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill x_{g,i}& =& n_{g,i}/\sum _{i=1}^N{n}_{g,i},x_{c,i}=0(i\in gas).\hfill \end{array}$$

(3)

From Equation (1), it can be seen that the balance is drawn around the cathode half cell and the gas compartment together. We neglect the solution of oxygen from the gas compartment in the bulk of the electrolyte in the cathode compartment, because the solubility of O${}_2$ in water is only about 9 mg/L. Thus, at a moderate flow rate of 10 mL/min, there would be only 0.09 mg/min of dissolved oxygen leaving the cathode compartment. This can be neglected in the component balance because the system will be fed with about 5.3 g/min O${}_2$ ($\gg 0.09$ mg/min). Evaporation of liquids into the gas compartment is neglected for similar reasons. Both effects would require a sub-model for mass transport through the gas diffusion electrode, which would be very difficult to parametrize realistically. So, the cost-benefit relation is insufficient for the type of overall process model that is set up here.

The concentrations in the bulk of the cathode half cell $c_{c,i}$ and at the cathode surface are given by the equations

$$\begin{array}{ccc}& & c_{c,i}=\frac{n_{c,i}}{V_c},\hfill \end{array}$$

(4)

$$\begin{array}{ccc}& & \frac{D_i}{{\delta }_{c,i}}A\left(c_{c,i}^s-c_{c,i}\right)=\sum _j{\nu }_{c,ij}{v}_{c,j}(i\notin gas),\hfill \end{array}$$

(5)

$$\begin{array}{ccc}& & c_{c,i}^s=H_{g,i}{p}_{g,i}(i\in gas),\hfill \end{array}$$

(6)

respectively. The concentrations $c_{c,i}^s$ of the non-gaseous components at the surface are calculated from the linear relation (5), which is a local material conservation condition at the surface with a simplified version of Fick’s law in its differential form describing the molar flow rate through the boundary layer. It includes the assumption of spatially constant concentration gradients for the species diffusing towards the surface or away from it. Furthermore, it implies that the concentration profile adapts instantaneously to local concentration changes. The distinction between bulk and surface concentration in 0D is used as a feature to make the model more versatile by applying different values for the thickness of the boundary layer ${\delta }_{c,i}$. For water as the solvent, the limit ${\delta }_{c,H_{2}O}\to 0$ is reasonably used (i.e., $c_{c,H_{2}O}=c_{c,H_{2}O}^s$). The concentrations of the gaseous components in the reactive zone (only O${}_2$ in the process) follow from Henry’s law (Equation (6)) for solution of gases in liquids, since those species are fed to the reaction zone from the gas compartment side.

In order to match the number of equations and the number of calculated variables, another equation for the determination of $F_{cout}$ is required. For a realistic description of the reactor, it is necessary to consider the Equation

$$V_c=\sum _{i\in liq}{V}_{cm,i}{n}_{c,i},$$

(7)

requiring the gap in the half cell to be steadily brimfull with liquid. However, with this equation used to calculate $F_{cout}$, the model converged poorly in our tests. Instead, the equation is used to determine $n_{c,H_{2}O}$, and hence, the dynamic Equation (1) is not used for water. To obtain an equation for $F_{cout}$, Equation (7) is differentiated with respect to time, and the resulting expression for ${\dot{n}}_{c,H_{2}O}$ is set equal to its analogue arising from the unused molar balance (1) for water. This yields the incompressibility condition

$$\begin{array}{ccc}\hfill 0& =& F_{cin}{x}_{cin,H_{2}O}-F_{cout}{x}_{c,H_{2}O}+N_{m,H_{2}O}\hfill \\ & +& \sum _j{\nu }_{c,H_{2}O,j}{v}_{c,j}+\frac{1}{V_{cm,H_{2}O}}\sum _{\stackrel{\mathrm{i}\in \mathrm{liq}}{i\ne H_{2}O}}{V}_{cm,i}{\dot{n}}_{c,i}\hfill \end{array}$$

(8)

(with ${\dot{n}}_{c,i}$ in the last sum from Equation (1)) for calculation of $F_{cout}$. Thus, the molar balance for water is contained implicitly in the equations. $F_{gout}$ follows from the phenomenological relation

$$F_{gout}=\frac{Kv}{V_{gm}}\sqrt{\frac{{\rho }_0}{{\rho }_g}\left(\frac{p_g}{p_0}-1.0\right)}$$

(9)

with a slight pressure drop being assumed across the outlet valve of the gas compartment.

2.1.2. Molar Balances in the Anode Half Cell

The balance equations on the anode side are analogue to the cathode balances, except for the gas compartment, which does not exist for the anode half cell. The molar component balance for the anode reads:

$$\frac{\mathrm{d}{n}_{a,i}}{\mathrm{d}t}=F_{ain}{x}_{ain,i}-F_{aout}{x}_{a,i}-N_{m,i}+\sum _j{\nu }_{a,ij}{v}_{a,j}(i\ne H_{2}O),$$

(10)

with the mole fractions given by

$$x_{a,i}=n_{a,i}/\sum _{i=1}^N{n}_{a,i}.$$

(11)

Again, the bulk and surface concentrations are obtained from

$$\begin{array}{ccc}& & c_{a,i}=\frac{n_{a,i}}{V_a},\hfill \end{array}$$

(12)

$$\begin{array}{ccc}& & \frac{D_i}{{\delta }_{a,i}}A\left(c_{a,i}^s-c_{a,i}\right)=\sum _j{\nu }_{a,ij}{v}_{a,j},\hfill \end{array}$$

(13)

with $c_{a,H_{2}O}=c_{a,H_{2}O}^s$. In contrast to the cathode side, dissolved oxygen emerging from oxidation of water occurs on the anode side in the process. The anode half cell is supposed to be entirely filled with liquid, too:

$$V_a=\sum _{i\in liq}{V}_{am,i}{n}_{a,i},$$

(14)

meaning that oxygen bubble formation is not considered here. Such bubbles increase the overpotential by blocking active surface sites. Nouri-Khorasani et al. report changes in the overpotential of 28…43 mV in this respect [48]. Neglecting this effect is reasonable because the cell in this article works at about 1.6 V ($\gg 43$ mV) under moderate operating conditions. The molar balance (10) for water is used only implicitly during the derivation of the incompressibility condition

$$\begin{array}{ccc}\hfill 0& =& F_{ain}{x}_{ain,H_{2}O}-F_{aout}{x}_{a,H_{2}O}-N_{m,H_{2}O}\hfill \\ & +& \sum _j{\nu }_{a,H_{2}O,j}{v}_{a,j}+\frac{1}{V_{am,H_{2}O}}\sum _{\stackrel{\mathrm{i}\in \mathrm{liq}}{i\ne H_{2}O}}{V}_{am,i}{\dot{n}}_{a,i}\hfill \end{array}$$

(15)

for calculation of $F_{aout}$.

There are two contributions to $N_m$. The oxonium ions H${}_3$O${}^{+}$ are the only charged particles traversing the PEM. Moreover, water is the only non-charged component being conducted through the membrane by electro-osmotic drag. Then, we have

$$\begin{array}{ccc}\hfill N_{m,H_3{O}^{+}}& =& \frac{I}{F},\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill N_{m,H_{2}O}& =& \frac{I}{F}\lambda ,\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill N_{m,i}& =& 0(i\ne H_3{O}^{+},H_{2}O)\hfill \end{array}$$

(18)

with the drag coefficient $\lambda$.

2.1.3. Reaction Kinetics and Voltage

The electric current $I_j$ across the electrode surface resulting from reaction j is calculated using the Butler–Volmer equation [6,49]:

$$I_j=AF{l}_j{k}_j\left(c_{ox,j}^s{\mathrm{e}}^{-\frac{l_j\alpha F}{RT}(E-E_j^0)}-c_{red,j}^s{\mathrm{e}}^{\frac{l_j(1-\alpha )F}{RT}(E-E_j^0)}\right).$$

(19)

The sum of this equation over all electrochemical reactions for one electrode yields the condition

$$I=\sum _j{I}_j,$$

(20)

which determines the electrode potential E for the case of galvanostatic operation. Equation (20) only holds if all occurring electrochemical reactions are known. Equation (19) enables us to give the Faraday efficiencies

$${\eta }_j=\frac{I_j}{I}$$

(21)

explicitly. Note that these Faraday efficiencies are not identical to the $\eta$ usually resulting from experiments, where the molar flow rate of the product species is measured and divided by its theoretical maximum according to the electric current. These experimental values correspond to a species rather than a reaction. They are averaged over all reactions involving the species, and they become equal to ${\eta }_j$ in case of only one reaction, or more precisely, in the case where the respective component occurs only in reaction j. Apart from that, the electric reaction currents $I_j$ are used to calculate the reaction rates from Faraday’s law of electrolysis:

$$v_j=\frac{I_j}{l_{j}F}.$$

(22)

The formal potential $E_j^0$ is calculated from the Nernst equation, and the equilibrium potential $E_j^{eq}$, i.e., the electrode potential in case of only reaction j and zero electric current, is derived from $E_j^0$ accordingly:

$$\begin{array}{ccc}\hfill E_j^0& =& E_j^o-\frac{RT}{l_{j}F}\sum _{i}ln({\gamma }_{ij}^s{)}^{{\nu }_{ij}},\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill E_j^{eq}& =& E_j^0-\frac{RT}{l_{j}F}\sum _{i}ln{\left(\frac{c_{ij}^s}{c_0}\right)}^{{\nu }_{ij}}.\hfill \end{array}$$

(24)

Values for potentials are given with reference to the standard hydrogen electrode (SHE) in this work. Equations (19)–(24) are evaluated separately for each electrode. However, they are only written down here once, because they are the same in both cases.

To determine the overall voltage drop over the cell, we require the Ohmic voltage drop

$$U_{Ohm}=I(R_{liq}+R_m)$$

(25)

with the resistance of the electrolyte in the gaps and the resistance of the membrane:

$$\begin{array}{ccc}\hfill R_{liq}& =& \frac{d_a+d_c}{{\sigma }_{liq}A},\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill R_m& =& \frac{d_m}{{\sigma }_{m}A}.\hfill \end{array}$$

(27)

Furthermore, the overpotential associated with reaction j is defined as

$$E_j^{over}=E-E_j^{eq},$$

(28)

which is, just like the Faraday efficiencies from (21), a measure of efficiency for the cell.

In the literature, the overall cell voltage U is usually composed of the open circuit voltage $U_0$, the overpotentials and the Ohmic losses. However, in the case of multiple active reactions, it is not straightforward to calculate $U_0$ and the total overpotentials at each electrode. Using Equations (23) and (26) from [20] instead, U can be written in terms of $E_a$ and $E_c$:

$$U=U_{Ohm}+E_c-E_a,$$

(29)

where the potential difference $E_c-E_a$ consists of the open circuit voltage and the total overpotentials. This form can be evaluated because $E_a$ and $E_c$, which are not reaction-specific, are calculated from (19) together with (20).

Finally, the electric power consumption of the cell is

$$P_{el}=UI.$$

(30)

In this convention, the electric current I has a positive algebraic sign. The applied electrical power enables the necessary increase in the Gibbs free energy and thus, the conversion of chemicals. However, it partly gets lost due to Joule heating and overpotentials, which are caused by polarization effects at the electrodes [28,49].

2.1.4. Reactions

The model equations given so far are generic in the sense that they do not correspond to a particular electrochemical reaction mechanism. In order to refer to the process of H${}_2$O${}_2$ production, the following reactions were implemented together with the model.

Cathode:

$$\begin{array}{ccc}\hfill {\mathrm{O}}_2+2{\mathrm{H}}_3{\mathrm{O}}^{+}+2{\mathrm{e}}^{-}& \stackrel{k_1}{⟶}& 2{\mathrm{H}}_2\mathrm{O}+{\mathrm{H}}_2{\mathrm{O}}_2\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill {\mathrm{O}}_2+4{\mathrm{H}}_3{\mathrm{O}}^{+}+4{\mathrm{e}}^{-}& \stackrel{k_2}{⟶}& 6{\mathrm{H}}_2\mathrm{O}\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill {\mathrm{H}}_2{\mathrm{O}}_2+2{\mathrm{H}}_3{\mathrm{O}}^{+}+2{\mathrm{e}}^{-}& \stackrel{k_3}{⟶}& 4{\mathrm{H}}_2\mathrm{O}\hfill \end{array}$$

(33)

Anode:

$$\begin{array}{ccc}\hfill 6{\mathrm{H}}_2\mathrm{O}& \stackrel{k_a}{⟶}& 4{\mathrm{H}}_3{\mathrm{O}}^{+}+4{\mathrm{e}}^{-}+{\mathrm{O}}_2\hfill \end{array}$$

(34)

The desired reaction (31) is the oxygen 2-electron reduction to hydrogen peroxide, which is often accompanied by the oxygen 4-electron reduction to water (Equation (32)) and a subsequent reduction of H${}_2$O${}_2$ (Equation (33)). Together, (31)–(33) constitute the oxygen reduction reaction (ORR) mechanism [30,31,50]. At the counter electrode, the anodic reaction of water electrolysis (Equation (34)) is used for the production of H${}_3$O${}^{+}$. These are the four active reaction channels for the process.

2.2. Membrane Distillation Model

The process includes downstream processing employing a module for membrane distillation to remove the added acid from the product solution. This separation process is principally based on differences in the chemical potential of species i on both sides of the membrane. Assuming the trans-membrane flux $N_{m,i}$ to be proportional to this difference, it can be expressed as [51]

$$N_{m,i}=B_i\left(p_{i,f}^0{\gamma }_{i,f}{x}_{i,f}-p_{i,p}^0{\gamma }_{i,p}{x}_{i,p}\right),$$

(35)

neglecting any non-ideal behavior of the gas phases on each side. The transport coefficients $B_i$ depend on the membrane thickness and on how efficiently species i is transferred through the membrane. The flux (35) is forced by a thermal gradient across the membrane. The feed stream to the feed side of the module is heated to a temperature $T_f$ well above the feed temperature $T_p$ on the permeate side. Consequently, the saturation vapor pressure $p_{i,f}^0$ of species i in the feed side of the pores is higher than the saturation vapor pressure $p_{i,p}^0$ of species i in the permeate side of the pores. The solutions on each side penetrate the membrane only to a certain extent but never reach the liquid coming from the other side due to the hydrophobic design of the pores. The difference in vapor pressures on each side causes a material flux through the gas-filled section of the pores. Again, the vapor pressures of the liquid components and the activity coefficients ${\gamma }_{i,f}$, ${\gamma }_{i,p}$ at compositions $x_{i,f}$, $x_{i,p}$ are obtained using the corresponding methods from Aspen’s ENRTL-RK set. The pressures of gaseous species are calculated using Henry’s law.

The outlet compositions $x_{out,i,f}$, $x_{out,i,p}$ are determined by the component balances

$$\begin{array}{ccc}\hfill F_{out,f}{x}_{out,i,f}& =& F_{in,f}{x}_{i,f}-N_{m,i},\hfill \end{array}$$

(36)

$$\begin{array}{ccc}\hfill F_{out,p}{x}_{out,i,p}& =& F_{in,p}{x}_{i,p}+N_{m,i}\hfill \end{array}$$

(37)

on the feed and the permeate side, respectively, and the outlet flow rates $F_{out,f}$, $F_{out,p}$ follow from the sum of Equation (36) and (37) over all components, exploiting the summation conditions:

$$\begin{array}{ccc}\hfill F_{out,f}& =& F_{in,f}-\sum _i{N}_{m,i},\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill F_{out,p}& =& F_{in,p}+\sum _i{N}_{m,i}.\hfill \end{array}$$

(39)

2.3. Experimental Details

Experimental data were generated in order to adjust and validate the model. The measurements were conducted with the real electrochemical cell, which had a cross-sectional area of $A=100$ cm${}^2$. The thicknesses of the anode and cathode compartment are $d_a=7$ mm and $d_c=1.2$ mm, respectively, and the corresponding compartment volumes are calculated accordingly as $V_a=d_{a}A=70$ cm${}^3$ and $V_c=d_{c}A=12$ cm${}^3$. The half cells are separated by a fumasep^® F-10120-PK proton exchange membrane (PEM). The cathode was made of a hydrophobized commercial Freudenberg H23C8 gas diffusion layer, which was sputtered with platinum as catalyst in order to obtain a gas diffusion electrode. The gas compartment is fed with 175 Nml/min synthetic air. For the anode half cell, an electrolyte solution of 2 mol/L H${}_2$SO${}_4$ in water is used, which is fed with a rate of 70 mL/min and which is recycled after it passed the cell. The experiments were conducted at room temperature (about 22.5 ${}^{\circ }$C in the laboratory), and no temperature control or external cooling was applied to the cell. Furthermore, the measurements were executed using galvanostatic operation with an electric current of $I=2.37$ A. Due to the low values of the product purity x (below 1 %), the corresponding relative measurement error is rather large (up to 50 %). Thus, values of model parameters derived from the measured data are considered as estimates.

3. Results and Discussion

3.1. Model Adjustment

In this subsection, the reactor model is adjusted to the experimental data. Sensitivity analyses of the model performed in Aspen Plus^® showed that the electro-kinetic parameters from the Butler–Volmer Equation (19) have the highest impact on the output quantities of the model and hence, they are chosen as free variables in an optimization problem designed to adjust the model to the measurements. The applicability of the determined reaction kinetics after a potential upscaling of the process then depends on how realistically the other model parameters are chosen. Some of the remaining input parameters can be determined from literature, e.g., the standard potentials under standard conditions $E_j^o$.

Table 1. Parameter values used in the model adjustment calculations of the electrochemical reactor for determination of the electro-kinetic rate constants $k_1$, $k_2$ and $k_3$.

Parameter	Value	Unit	Type/Source
A	100	cm${}^2$	geometry
$d_a$	7	mm
$d_c$	1.2	mm
$V_a$	70	cm${}^3$
$V_c$	12	cm${}^3$
$c_{a,H_{2}S{O}_4}$	2	mol/L	experimental design
$c_{c,H_{2}S{O}_4}$	2	mol/L
$F_{ain}$	70	ml/min
$F_{gin}$	35	ml/min
I	2.37	A
p	1.01325	bar
T	22.5	${}^{\circ }$C
$E_1^o$	0.685	V vs. SHE	[31]
$E_2^o$	1.229	V vs. SHE	[30,50]
$E_3^o$	1.76	V vs. SHE	[50]
$E_a^o$	0	V vs. SHE
$H_{g,O_2}$	1.3	mol/(bar·m${}^3$)	[52] (H${}_2$O as solvent)
$d_m$	127.5	$\mathsf{\mu }$m	from data sheet
${\sigma }_m$	5.44	S/m	for 0.5 M H${}_2$SO${}_4$
${\alpha }_1$	0.2		estimate
${\alpha }_2$	0.01
${\alpha }_3$	0.01
${\alpha }_a$	0.5
$D_i/{\delta }_{a,i}$ ($i\ne$H${}_2$O)	0.05	cm/s
$D_i/{\delta }_{c,i}$ ($i\ne$H${}_2$O, O${}_2$)	0.05	cm/s
$k_a$	0.001	cm/s
$Kv$	0.1	m${}^3$/h
$\lambda$	6
${\sigma }_{liq}$	1.5	S/m
NLPQLP tolerance	0.001		convergence
gradient accuracy	0.0001
simulation accuracy	${10}^{-8}$
sandwiching quality	1.05

summarizes the parameter values used in the model adjustment of the electrochemical reactor. At the cathode, reactions (31)–(33) are active and they are enumerated consecutively, i.e., the reaction index j (used e.g., for ${\alpha }_j$ and $E_j^o$) is 1 for the reaction from (31), 2 for (32) and 3 for (33). At the anode, the water electrolysis reaction (34) is active. The parameters $d_m$, ${\sigma }_m$ concerning the membrane are taken from the data sheet provided by the manufacturer (see Table 1). The reference quantities are set to $p_0=1.01325$ bar, $c_0=1$ mol/L and ${\rho }_0=1$ g/cm${}^3$. The true component approach is used to specify the acid concentration, i.e., only the mole fraction of H${}_2$SO${}_4$ is entered and the employed property methods provided by Aspen calculate the corresponding ionic concentations. All feed temperatures and pressures to the reactor are set to 22.5 ${}^{\circ }$C and 1.01325 bar, except for the feed pressure to the gas compartment, which is set to 1.02 bar. The gas compartment is fed with pure oxygen in the simulation, and the associated feed rate is set to 20 % of the experimental value, which is the percentage of oxygen in synthetic air. Still, there are some model parameters to be estimated, e.g., the electro-osmotic drag coefficient is chosen as $\lambda =6$ from experimental experience. In previous work for the ORR mechanism at PdAu${}_3$ surfaces, the charge transfer coefficients were calculated [34]. Here, they must be in the same order of magnitude since the $E_j^o$ do not differ.

The model parameters obtained from model adjustments here result from the solution of a least squares problem. For the electrochemical reactor, the product purity x and the Faraday efficiency $\eta$ were measured as functions of the feed rate to the cathode half cell. The objective function f to be minimized is the weighted sum over all normalized and squared differences between the measured values and the values predicted by the model at given cathode feed rate $F_{cin}$:

$$f=w_x{S}_x+w_{\eta }{S}_{\eta },$$

(40)

with $S_x$, $S_{\eta }$ given by

$$\begin{array}{ccc}\hfill S_x& =& \sum _{l=1}^{N_x}{\left(\frac{x_{sim,l}-x_{exp,l}}{{\overline{x}}_{exp}}\right)}^2,\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill S_{\eta }& =& \sum _{l=1}^{N_{\eta }}{\left(\frac{{\eta }_{sim,l}-{\eta }_{exp,l}}{{\overline{\eta }}_{exp}}\right)}^2.\hfill \end{array}$$

(42)

The terms are normalized to their mean experimental value. Here, the simulated Faraday efficiency is calculated as in the experiments, i.e., the product purity x is recorded, multiplied by the set point for the flow rate $F_{cin}$ and finally divided by the theoretical maximum $I/\left(2F\right)$ for the case without side reactions:

$${\eta }_{sim}=\frac{x{F}_{cin}}{I/\left(2F\right)}.$$

(43)

The solver NLPQLP from the Schittkowski suite [53] is used to minimize f by varying the optimization variables $k_1$, $k_2$ and $k_3$.

Since electro-kinetic rate constants may vary by many orders of magnitude [49], they are replaced using the nonlinear transformation

$$k_j={10}^{-b_j}$$

(44)

and the $b_j$ are used as optimization variables instead. This way, the optimization solver calculates steps in a logarithmically scaled design space and it does not need to distinguish between steps of different orders of magnitude. Without this transformation, no meaningful results could be found.

Table 2. Starting values and boundaries for the optimization variables $b_j$ estimated by adapting the electrochemical reactor model to measured data.

Variable	Starting Value	Lower Boundary	Upper Boundary
$b_1$	7.0	3.0	8.0
$b_2$	7.0	6.7	8.0
$b_3$	7.0	5.0	8.0

shows the employed starting values and boundaries for the optimization variables.

Particular attention needs to be paid to the different convergence parameters when working in the framework of flowsheet optimization. The error tolerance of the optimization solver and the accuracy of gradient calculations must be less strict than the accuracy of the simulation. Otherwise, the information the solver is looking for in the function evaluations will simply disappear in numerical noise. Here, the tolerance of NLPQLP is set to ${10}^{-3}$, the relative accuracy of the gradients is ${10}^{-4}$ and the flowsheet tolerance is ${10}^{-8}$. The latter is used as accuracy requirement for mass balances, flash calculations and for the system of equations itself, which is solved in equation-oriented mode.

A common choice for $w_x$, $w_{\eta }$ is the reciprocal variance of the corresponding experimental data [54], ensuring that more precisely measured quantities are granted more impact on the objective function. However, the variances are often unknown, as is the case here, too. Hence, the weighting factors can be chosen out of a continuum between the extreme compromises $w_x=0$, $w_{\eta }>0$ and $w_x>0$, $w_{\eta }=0$, where the optimization variables are adapted to only one quantity, respectively. We use the sandwiching algorithm [43] to adaptively determine distinct sets of weighting factors, for each of which the optimization problem described above is solved. The result is a finite number of solutions, each corresponding to a different estimate for the $k_j$. These points represent a Pareto set [47] with well-defined sandwiching quality in the objective space spanned by the criteria $S_x$ and $S_{\eta }$. A similar procedure has also been presented in [46] for solving reconciliation problems.

Aspen Plus V10^® is coupled to the NLPQLP solver and the sandwiching subroutines via MS Excel. A VBA code is used to manage the data transfer, and a self-designed (Excel-based) graphical interface is used for problem specification and automated formatting of results. The framework is similar to the one published in [42], except for CHEMCAD as a flowsheet simulator being replaced by Aspen here. The applied HappLS class from the Aspen Plus GUI 36.0 Type Library exposes all required properties and methods to access flowsheet data from within VBA.

Figure 2 shows the Pareto frontier of the MCO problem for adjustment of the electrochemical reactor model to measured data displayed in the objective space spanned by the sums $S_x$ and $S_{\eta }$. Within the given approximation quality (see bottom of Table 1), this frontier is represented by the calculated Pareto points and the straight lines connecting them. It is the set separating the feasible and infeasible parts in the objective space. In the plot, the extreme compromises are marked corresponding to $w_x=0$ (bottom right) and $w_{\eta }=0$ (top left). For all six calculated points, the results including the respective weighting factors $w_x$, $w_{\eta }$ determined by the sandwiching algorithm are listed in

Table 3. Results including the respective weighting factors $w_x$, $w_{\eta }$ for the Pareto points from the multi-criteria optimization problem for adjustment of the electrochemical reactor model to measured data. The points are plotted in the objective space in Figure 2.

Solution no.	${\mathit{w}}_{\mathit{x}}$	${\mathit{w}}_{\mathit{\eta }}$	${\mathit{S}}_{\mathit{x}}$	${\mathit{S}}_{\mathit{\eta }}$	${\mathit{b}}_1$	${\mathit{b}}_2$	${\mathit{b}}_3$
1	0.2253	0.0000	0.1079	0.2894	7.274	7.546	6.012
2	0.0000	3.8629	0.2671	0.0198	7.270	7.130	6.323
3	0.2252	0.1329	0.1595	0.0376	7.326	7.223	6.172
4	0.2253	0.0461	0.1334	0.0909	7.341	7.314	6.110
5	0.2125	1.2819	0.1958	0.0225	7.330	7.163	6.252
6	0.2253	0.0289	0.1244	0.1299	7.237	7.351	6.081

.

A comparison of the experimental data to the predictions of the adapted models is exhibited in Figure 3. The measured data represented by black dots are displayed as a function of the feed rate to the cathode half cell $F_{cin}$. The extreme compromises $w_x=0$, $w_{\eta }=0$ corresponding to solutions 1 and 2 in Table 3 were chosen for the plot. I.e., the $b_j$ for the two solutions were inserted into Equation (44) to obtain the respective electro-kinetic rate constants $k_j$ for utilization in the cell model. Subsequently, the models were solved for 40 equidistant values of $F_{cin}$ in the interval 0.3…7 mL/min, in order to reach a quasi-continuous representation for the lines in Figure 3. For both parametrizations, the model from Section 2.1 achieves a reasonable approximation of the experimental data. The simulation results for $w_x=0$ (solid line) deviate more strongly from the measurements for x, especially at low feed rates, and the simulation results for $w_{\eta }=0$ agree less well with the measurements for $\eta$, particularly at larger feed rates. This is the expected result because e.g., for $w_x=0$, the deviations accumulated in $S_x$ from (41) are not taken into account in the objective function f from (40). Hence, the corresponding measurement points cannot force the curve to pass them by more closely in the optimization procedure. For $F_{cin}$ between 0.3…2 mL/min, the model reproduces the experimental values for $\eta$ almost exactly for both parametrizations. The results of solution 1 ($w_{\eta }=0$) are chosen for the model used in the subsequent section about process optimization, because there, a more accurate prediction of the product purity is desired, since the considered criteria partly depend on x, but not on $\eta$. Using $j_0=Fk{\left(c_O^b\right)}^{1-\alpha }{\left(c_R^b\right)}^{\alpha }$ [49], an exchange current density of $j_0\approx 5$ $\mathsf{\mu }$A/cm${}^2$ can be found for the corresponding kinetic parameters. This is comparable to the lowest values found in [55] for the Pt/Nafion interface.

3.2. Process Optimization

3.2.1. Flowsheet Parameter Settings

Figure 4 shows the flowsheet of the process for the electrosynthesis of H${}_2$O${}_2$ introduced in this work. The reactor model described in Section 2.1 and parametrized in Section 3.1 is integrated as the core unit EC (1). The design specifications ACICONTA and ACICONTC control the acid feeds ACIFEEDA and ACIFEEDC, such that the respective direct feeds FIRI-16 and FIRI-18 to the anode and cathode half cells have H${}_2$SO${}_4$ concentrations of 0.5 mol/L. Both acid feed streams have 10 wt.-% H${}_2$SO${}_4$. Accordingly, H2OCONTA and H2OCONTC control the water feed streams H2OFEEDA and H2OFEEDC. For the cathode, this flow rate is variable, and for the anode it is set to 0.1 kmol/h. Due to the anode recycle loop, the streams ACIFEEDA and H2OFEEDA are reduced to very small flow rates by the solvers, because most material is recycled. H${}_3$O${}^{+}$ and H${}_2$O can leave this loop through the PEM in EC, and the gas separation O2SEP (3) was inserted to remove all oxygen emerging from reaction (34) at the anode, which would otherwise become highly accumulated. HEAT1 (4) establishes a variable temperature difference between both sides of the membrane distillation module MD1 (2). PFCONT varies the stream PERMFEED, such that stream AUX05 has a total flow rate of 0.01 kmol/h. The flash H2O2CONC (4) for concentration of the product solution is specified with a temperature of $T_c=110{}^{\circ }$C and a vapor fraction of $y_c=0.5$ mol/mol. H${}_2$O${}_2$ enters the downstream recycle loop through the membrane in MD1 and leaves it as the high-boiling component through the liquid outlet TOFPROC from H2O2CONC. Both cooling units diminish the temperature to 22 ${}^{\circ }$C and all pumps raise the pressure by 5 hPa except for P04, which is specified to a 10 hPa pressure increase. The oxygen feed O2FEEDC has a flow rate of 0.01 kmol/h at 22 ${}^{\circ }$C and 1.025 bar, which is decreased to 1.02 bar by the valve V3. All other feeds have a pressure of 1.01325 bar, also at 22 ${}^{\circ }$C. The $B_i$ for MD1 are taken as 1 kmol/(bar$·$h), except for water and H${}_2$O${}_2$, for which they were estimated to $B_{H_{2}O}=0.01$ kmol/(bar$·$h) and $B_{H_2{O}_2}=7$ kmol/(bar$·$h) from experimental data. The geometry, electro-kinetic and membrane parameters, as well as the estimated parameters for EC, are taken as in Table 1, and the electro-kinetic rate constants for the Butler–Volmer equation are derived from Table 3 (row with $w_{\eta }=0$) using Equation (44). The flowsheet simulation takes into account the downstream section, recycle streams and the impact of the interactions between unit operations. With all of its structure, specification and parametrization, it is supposed to be the digital twin of the laboratory process, capable of replacing experiments by numerical simulations for finding optimal points of operation.

3.2.2. Flowsheet Convergence

The flowsheet is solved in sequential-modular mode. Two steps turned out to be crucial in order to obtain convergence. First, O2SEP had to be inserted for the above-mentioned reason. Second, the design specifications were nested inside the loops converging the tear streams FBA1 and AUX04. Otherwise, the tear streams cannot converge, because, e.g., the anode recycle loop would have to be closed with the initial flow rate of H2OFEEDA, which is far too large for the amount of substance in the loop to remain finite. Tear streams as well as design specifications are converged with Broyden solvers. The relative tolerance for AUX04 is 10${}^{-4}$, while it could be reduced to 10${}^{-8}$ for FBA1. The accuracies of the design specifications in the upstream processing are 10${}^{-5}$, and for PFCONT it is 10${}^{-3}$. The error tolerances for flash convergence, entropy balance and for fugacity calculations are set to 10${}^{-8}$, and the accuracy of mass balances is 10${}^{-4}$.

3.2.3. Multi-Criteria Optimization Problem

For the optimization problem, the solver accuracy was set to 10${}^{-3}$, numerical derivatives were evaluated using a relative perturbation of 3% and the quality of the sandwiching approximation was 1.05.

Table 4. Configuration of the MCO problem for optimization of working points for the process producing H${}_2$O${}_2$ represented by the flowsheet displayed in Figure 4.

Variable	Type	Unit	Init. Value
I	var.	A	2.37
$F_{cin}$	var.	kmol/h	0.025
$T_{MD}$	var.	10 ${}^{\circ }$C	70
$x_{out}$	max.	wt.-%	0.7164
$F_{cout}$	max.	g/min	7.9769
$P_{el,tot}$/$N_{out}$	min.	kWh/kg	4.618
Heat consumption	min.	kW	0.2099
O${}_2$ consumption	min.	mg/min	6.955

provides an overview over the MCO problem designed to analyze the process and to calculate Pareto-optimal working points for it. The optimization variables are I, $F_{cin}$ and the temperature $T_{MD}$ of the stream to the feed side of the membrane distillation unit MD1. The latter quantity is a specification for HEAT1. I is varied in the interval 0.5…10 A. At $I=0.5$ A there is almost only 0.05 wt.-% H${}_2$O${}_2$ left in the main outlet TOFPROC, and we do not wish to go below that. Apart from that, it should be noted that for I as high as 10 A, a notable production of H${}_2$, which is not considered here, sets in in the real system. Furthermore, at such high electric currents, external cooling of the reactor is required, which is not considered here, either. The boundaries for $F_{cin}$ are 0.01…1 kmol/h according to a reasonable interval of $\approx 3\dots 300$ mL/min. $T_{MD}$ is kept between 23…98 ${}^{\circ }$C, so that the feed temperature for MD1 is always higher than the permeate temperature but without evaporating the solution. The optimization solver works with $T_{MD}$ measured in units of 10 ${}^{\circ }$C because then the actual optimization variable is $\mathcal{O}\left(1\right)$ within its boundaries, which turned out to be more stable. Among the five objectives, three quantities have not been introduced yet: the overall electric power consumption $P_{el,tot}$ is taken as the sum of the cell power $P_{el}$ and the powers of all pumps, and the ratio of this quantity to the component flow rate $N_{out}$ of H${}_2$O${}_2$ in TOFPROC is a measure for the electrical energy dissipated per kg peroxide output. $P_{el}$ takes its lowest value of 0.4 W for $I=0.5$ A, $F_{cin}=0.01$ kmol/h. $x_{out}$ is the overall product purity, i.e., the mass fraction of H${}_2$O${}_2$ in the main outlet. The heat consumption is the sum of the net duties of HEAT1 and COOL1 for separation of acidic components and the net duties of the flash H2O2CONC and COOL2 in the recycle loop for product concentration. The O${}_2$ consumption is taken as the mass flow rates of O2FEEDC minus O2OUTA and O2OUTC, since the two latter streams are considered as recycled.

3.2.4. Multi-Criteria Optimization Results and Discussion

The highest H${}_2$O${}_2$ purities in TOFPROC were found at about 3 wt.-%. The outcome of the Pareto approximation shall be discussed here within the context of two-dimensional projections in the objective space. Figure 5a shows such projections for the heat consumption and the O${}_2$ consumption of the entire process with $P_{el,tot}/N_{out}$ on the horizontal axis. The data sets for Figure 5 are given in

Table 5. Objectives and optimization variables selected from the Pareto set resulting from MCO for the process. Only non-dominated points, plotted in Figure 5, are given here.

Projection of Heat Consumption on${\mathit{P}}_{\mathit{el},\mathit{tot}}/{\mathit{N}}_{\mathit{out}}$(Figure 5a, Left Vertical axis):
${\mathit{P}}_{\mathit{el},\mathit{tot}}/{\mathit{N}}_{\mathit{out}}$	Heat Consumption	I	${\mathit{F}}_{\mathit{cin}}$	${\mathit{T}}_{\mathit{MD}}$
[kWh/kg]	[kW]	[A]	[kmol/h]	[${}^{\circ }$C]
3.778	0.444	0.5	0.09	97.943
5.028	0.200	0.5	0.019	70.080
5.420	0.194	0.5	0.015	70.000
17.168	0.158	10.0	0.012	55.224
65.408	0.137	0.5	0.01	23
Projection of O${}_2$ Consumption on ${\mathit{P}}_{\mathit{el},\mathit{tot}}/{\mathit{N}}_{\mathit{out}}$ (Figure 5a, Right Vertical axis):
${\mathit{P}}_{\mathit{el},\mathit{tot}}/{\mathit{N}}_{\mathit{out}}$	O${}_{\mathbf{2}}$ Consumption	$\mathit{I}$	${\mathit{F}}_{\mathit{cin}}$	${\mathit{T}}_{\mathit{MD}}$
[kWh/kg]	[mg/min]	[A]	[kmol/h]	[${}^{\circ }$C]
3.778	0.750	0.5	0.09	97.943
5.028	0.600	0.5	0.019	70.080
5.420	0.558	0.5	0.015	70.000
6.210	0.505	0.5	0.011	70.727
65.408	0.469	0.5	0.01	23
Projection of ${\mathit{P}}_{\mathit{el},\mathit{tot}}/{\mathit{N}}_{\mathit{out}}$ on ${\mathit{F}}_{\mathit{cout}}$ (Figure 5b, Left Vertical axis):
${\mathit{F}}_{\mathit{cout}}$	${\mathit{P}}_{\mathit{el},\mathit{tot}}/{\mathit{N}}_{\mathit{out}}$	$\mathit{I}$	${\mathit{F}}_{\mathit{cin}}$	${\mathit{T}}_{\mathit{MD}}$
[g/min]	[MWh/kg]	[A]	[kmol/h]	[${}^{\circ }$C]
28.123	0.004	0.5	0.09	97.943
311.043	0.138	2.381	1	70
311.043	3.232	2.38	1	23
Projection of O${}_{\mathbf{2}}$ Consumption on ${\mathit{F}}_{\mathit{cout}}$ (Figure 5b, Right Vertical axis):
${\mathit{F}}_{\mathit{cout}}$	O${}_{\mathbf{2}}$ Consumption	$\mathit{I}$	${\mathit{F}}_{\mathit{cin}}$	${\mathit{T}}_{\mathit{MD}}$
[g/min]	[mg/min]	[A]	[kmol/h]	[${}^{\circ }$C]
3.151	0.469	0.5	0.01	23
310.876	0.903	0.5	1	40.452
311.043	9.766	2.38	1	23

. Within the given approximation quality, the plotted lines represent the boundary between feasible and infeasible regions. Resulting points appearing as dominated in the projections (i.e., points where an improvement in both objectives is possible without becoming infeasible) in the corresponding projection have been removed from the plots. In Figure 5a, the transition from high to low heat and O${}_2$ consumptions while simultaneously worsening the ratio $P_{el,tot}/N_{out}$ along the Pareto sets in this diagram is achieved by reducing the feed rate to the cathode half cell from 0.09 to 0.01 kmol/h and by lowering $T_{MD}$. Of course, this will lead to smaller heat consumptions, but it also results in less H${}_2$O${}_2$ crossing the membrane in MD1 and hence in lower values of $N_{out}$.

Very sharp bends in the Pareto frontiers are found in Figure 5b. For the O${}_2$ consumption, the bending point in the bottom right can be reached from the point in the top right, mainly by decreasing the electric current from 2.38 to 0.5 A. This drastically reduces the amount of oxygen converted in electrochemical reactions. In fact, both objectives on the vertical axes can be decreased by more than an order of magnitude in this range, with almost no reduction in $F_{cout}$. Such interesting regions are sought after by and can only be found using MCO, while the sandwiching algorithm prevents an unnecessary amount of points being calculated in these intervals. Mind that only two criteria are regarded in this reflection. e.g., $x_{out}$ drops significantly with decreasing I. From the solution in the bottom left of Figure 5b, the bending point is reached by increasing $F_{cin}$. This notably improves $F_{cout}$, while the other two objectives in the projection worsen only marginally.

In order to assess the relation of electrical and thermal energy consumption, the reactor and the membrane distillation module are compared. At moderate operating conditions ($I=2.37$ A, $F_{cin}=0.025$ kmol/h, $T_{MD}=70{}^{\circ }$C), the cell power $P_{el}$ is 3.84 W. It mainly depends on the electrical current and is increased to 56.87 W by increasing I to its upper boundary at 10 A. The heat for membrane distillation is injected in HEAT1. Under the moderate conditions just mentioned, this heater operates at a thermal input of 25.59 W. This heat input depends on $T_{MD}$, and it is increased to 40.54 W by increasing $T_{MD}$ to its uppper boundary at 98 ${}^{\circ }$C. Hence, depending on the operating conditions, $P_{el}$ can be smaller or larger than the energy required for membrane distillation, and its range comprises that of the heat input in HEAT1. The overall heat consumption (see Table 5) exceeds the electrical power, but such a contribution also exists in the conventional anthraquinone process.

In Table 5, there are several non-identical values of $T_{MD}$ very close to 70 ${}^{\circ }$C. This was investigated using a sensitivity analysis in the range $T_{MD}=69.9\dots 70.1{}^{\circ }$C with $I=0.5$ A and $F_{cin}=0.01$ kmol/h. $F_{cout}$ and the O${}_2$ consumption do not depend on $T_{MD}$, the latter because there is no oxygen in the downstream. $P_{el,tot}/N_{out}$ does not vary notably in this range. However, $x_{out}$ as well as the heat consumption both vary by about 0.2% in this range. On the one hand, this variation is large enough so that, together with the given solver tolerance of 0.001, our algorithm can distinguish between the respective values in the table and they do not result from numerical noise. On the other hand, the variation is small enough for a robust implementation of the Pareto-optimal operating points in the real process. The challenge will then be to accurately establish the optimal temperature all over the feed side of the membrane in MD1.

4. Conclusions

This work introduces a process for electrosynthesis of H${}_2$O${}_2$ and focuses on its modeling and multi-criteria optimization. The proposed model equations for the required custom units pursue the goal of rapid convergence and a reasonable mathematical mapping of the main output quantities, particularly the product purity but also the Faraday efficiency. For calculation of the outlet flow rates of the reactor, a numerically stable incompressibility condition is derived. The Butler–Volmer kinetics are determined from adjustment to experimental data, which is realized by solving an MCO problem, since it is not clear in general how to weight different measured quantities against each other. The model is implemented in Aspen custom modeler^®, and its validated form is embedded in an Aspen Plus^® flowsheet. For the application of MCO to both, model adjustment as well as electrosynthesis on a process level, the availability of studies is limited so far.

The Pareto set obtained from MCO can be exploited by process operators for a quick and simple selection of efficient working points without the need for further simulations, and especially without time-consuming experiments. The MCO results show that the product purities of the process introduced here are higher than those of the anthraquinone process (before downstream). Hence, the proposed process is a promising approach to sustainable production of hydrogen peroxide. However, the product flow rates in the range of several mL/min cannot compete with the industrial process. However, they need not, because the process considered here is designed for decentralized applications. Nevertheless, the product amount should be increased in future work. In this respect, upscaling cannot be the (only) solution because the process is supposed to remain decentralized. However, stacking multiple electrochemical cells can be one way to achieve higher product flow rates (at the cost of more energy input of course).