Robust Design and Operation of a Multistage Reactor for Methanol Synthesis from Renewable Resources

Abstract

Methanol is an important raw material widely used in the chemical industry. This article addresses the challenge of fluctuations in green hydrogen as a feed stream for methanol production from renewable feedstock. For a staged reactor design, robust operating conditions are generated through the simultaneous steady-state optimization of 50 process scenarios. The feed can be split and fed separately to the different reactor stages. However, neglected transient effects may render this design infeasible under dynamic conditions concerning carbon conversion and reactor temperature constraints. To overcome this, an additional dynamic optimization is conducted to ensure a feasible operation by an optimized feed-forward control of feed distribution and shell temperatures. In practice, this is possible because the disturbance, i.e., fluctuation, is measurable and predictable in a short time frame. The optimization yields optimal operating conditions, resulting in a reactor that is dynamically feasible for measurable fluctuating inlet conditions.

1. Introduction

Methanol is a widely used and important raw material for the chemical industry [1,2]. It can also be used for chemical energy storage [1,3]. Today, methanol is produced large-scale using fossil-based synthesis gas in heterogeneously catalyzed reactions. The worldwide annual production is 98 million tonnes [4]. Of that, only a small fraction, i.e., $0.2$ million tonnes, is produced from renewable resources [4]. In the light of climate change, it becomes clear that green methanol will play an important role [5,6]. For the production of green methanol, the fossil-based raw materials need to be substituted by renewables. Possible carbon sources are, among others, biogas [7,8], industrial fumes [9], and direct CO${}_2$ air capture [10]. A suitable hydrogen source is electrolysis using electricity produced by solar, wind, or water power [8,11]. Due to the transformation of the energy supply system, there are already a number of pilot-scale plants for the production of green methanol in operation [12]. However, as of today, the costs for renewable methanol are still much higher than for its fossil counterpart. This might change in the future due to legislation or more efficient processes [11,13,14].

One of the challenges in the production of green methanol is to cope with fluctuations in the feed streams [15,16,17], e.g., a varying availability of regenerative electricity and hence hydrogen. Conventional processes are operated under steady-state conditions with little to no fluctuations [2], therefore they can be operated very efficiently. For the design of green processes, the fluctuations need to be taken into account and hence the processes need to be more robust, i.e., the possible window of operation needs to be bigger. Robustness always comes with sub-optimality in a sense that the design needs to be more conservative [18]. This problem can be tackled by means of a multistage reactor design, as shown by Fischer et al. for methanation processes [19,20]. In the present paper, this concept is explored for methanol production.

This study aims at the robust optimization [18,21] of the operating conditions of a multistage reactor, consisting of three continuously stirred tank reactors, for the production of green methanol, where fluctuations in the hydrogen supply are taken into account. As stated above, a robust design always leads to conservative processes. Furthermore, neglecting the dynamic behavior of the process during robust optimization may lead to infeasible dynamic processes due to transient effects. A feed-forward control scheme consisting of a dynamic optimization in a short time horizon is proposed to overcome this limitation by manipulating the feed distribution and the shell temperatures.

The paper is structured as follows: The required model equations, assumptions, parameters, and the optimization problems are presented in Section 2. The obtained results are shown and discussed in Section 3. Finally, a short summary of the study, the results, and concluding remarks are given in Section 4.

2. Materials and Methods

The process investigated in this work is a heterogeneously catalyzed methanol synthesis using a conventional Cu/ZnO/Al₂O₃ catalyst. The process is widely used in industry, with operating conditions between [50, 60] bar and around 250 °C, typically in fixed-bed reactors [2,22,23].

2.1. Kinetic Model

The following three overall reactions take place:

$$\mathrm{CO}+2{\mathrm{H}}_2\stackrel{\rightharpoonup }{\leftharpoondown }{\mathrm{CH}}_3\mathrm{OH},$$

(1)

$${\mathrm{CO}}_2+3{\mathrm{H}}_2\stackrel{\rightharpoonup }{\leftharpoondown }{\mathrm{CH}}_3\mathrm{OH}+{\mathrm{H}}_2\mathrm{O},$$

(2)

$$\mathrm{CO}+{\mathrm{H}}_2\mathrm{O}\stackrel{\rightharpoonup }{\leftharpoondown }{\mathrm{CO}}_2+{\mathrm{H}}_2,$$

(3)

where Equation (1) is the CO hydrogenation, Equation (2) is the CO${}_2$ hydrogenation, and Equation (3) is the water gas shift reaction.

Established catalysts and kinetics are efficient for feeds consisting mainly of CO, which is also applied in this paper. This implies a two step process, where first CO₂ is converted to CO, which is then converted to methanol [24].

The reaction kinetics is taken from Seidel et al. [25,26], with parameters taken from [27] and given in

Table 1. Parameter values for reaction kinetics model [27].

Parameter	Value	Parameter	Value	Parameter	Value	Unit
${\alpha }_{1,1}$	$13.814$	${\alpha }_{1,2}$	$15.0921$	${\alpha }_{1,3}$	$1.2777$	-
${\alpha }_{2,1}$	3784.4	${\alpha }_{2,2}$	$1581.7$	${\alpha }_{2,3}$	$-2.167$	-
${\alpha }_{3,1}$	$-9.2833$	${\alpha }_{3,2}$	$-8.7639$	${\alpha }_{3,3}$	$0.5194$	-
${\alpha }_{4,1}$	$3.1475·{10}^{-3}$	${\alpha }_{4,2}$	$2.1105·{10}^{-3}$	${\alpha }_{4,3}$	$-1.037·{10}^{-3}$	-
${\alpha }_{5,1}$	$-4.2613·{10}^{-7}$	${\alpha }_{5,2}$	$-1.9303·{10}^{-7}$	${\alpha }_{5,3}$	$2.331·{10}^{-7}$	-
${\mathrm{A}}_{\mathrm{CO}}$	$-5.001$	${\mathrm{A}}_{{\mathrm{CO}}_2}$	$-3.145$	${\mathrm{A}}_{\mathrm{WGS}}$	$-4.4526$	$\mathrm{mol}{\mathrm{bar}}^{-3}{\mathrm{kg}}_{\mathrm{cat}}^{-1}{\mathrm{s}}^{-1}$
${\mathrm{B}}_{\mathrm{CO}}$	26.455	${\mathrm{B}}_{{\mathrm{CO}}_2}$	1.5308	${\mathrm{B}}_{\mathrm{WGS}}$	15.615	-
${\mathsf{\beta }}_{{\mathrm{H}}_2}^{\left(\mathrm{red}\right)}$	1.1064	-	-	-	-	bar${}^{-0.5}$
${\mathsf{\beta }}_{\mathrm{CO}}^{\left(\mathrm{oxi}\right)}$	0.14969	${\mathsf{\beta }}_{{\mathrm{CO}}_2}^{\left(\mathrm{het}\right)}$	0.062881	-	-	bar${}^{-1}$
$\Delta {\mathrm{G}}_1$	$335.7$	$\Delta {\mathrm{G}}_2$	$21841.5$	-	-	J mol${}^{-1}$
${\mathrm{k}}_1^{+}$	$79.174·{10}^{-4}$	${\mathrm{k}}_2^{+}$	$0.188·{10}^{-4}$	-	-	s${}^{-1}$
${\mathsf{\varphi }}_{\max }$	0.9	-	-	-	-	-

. It is modeled as a Langmuir–Hinshelwood kinetics with three types of catalyst centers: heterolytic (het), reduced (red), and oxidized (oxi). The availability of those centers depends on the oxidization state of the catalyst, $\varphi$. As different centers are necessary for the three reactions, there is a direct influence of the catalyst state on the reaction rates $\mathit{r}$,

$$r_{\mathrm{CO}}=(1-\varphi )·k_1·\left(p_{\mathrm{CO}}·p_{{\mathrm{H}}_2}^2-{\displaystyle \frac{p_{{\mathrm{CH}}_3\mathrm{OH}}}{K_1}}\right)·{\Theta }^{\left(\mathrm{oxi}\right)}·{{\Theta }^{\left(\mathrm{het}\right)}}^4,$$

(4)

$$r_{{\mathrm{CO}}_2}={\varphi }^2·k_2·\left(p_{{\mathrm{CO}}_2}·p_{{\mathrm{H}}_2}-{\displaystyle \frac{p_{{\mathrm{CH}}_3\mathrm{OH}}·p_{{\mathrm{H}}_2\mathrm{O}}}{K_2·p_{{\mathrm{H}}_2}}}\right)·{{\Theta }^{\left(\mathrm{red}\right)}}^2·{{\Theta }^{\left(\mathrm{het}\right)}}^4,$$

(5)

$$r_{\mathrm{WGS}}={\displaystyle \frac{\varphi }{1-\varphi }}·k_3·\left(p_{{\mathrm{CO}}_2}-{\displaystyle \frac{p_{\mathrm{CO}}·p_{{\mathrm{H}}_2\mathrm{O}}}{K_3·p_{{\mathrm{H}}_2}}}\right)·{\Theta }^{\left(\mathrm{red}\right)}·{\Theta }^{\left(\mathrm{oxi}\right)},$$

(6)

where $\mathit{p}$ are the partial pressures of components $i\in [{\mathrm{CH}}_3\mathrm{OH},{\mathrm{CO}}_2,\mathrm{CO},{\mathrm{H}}_2,{\mathrm{H}}_2\mathrm{O}]$, $\mathit{K}$ are temperature dependent equilibrium constants,

$$\begin{array}{cc}\hfill log\left(\mathit{K}\right)& ={\mathbf{\alpha }}_{\mathbf{1}}+{\mathbf{\alpha }}_{\mathbf{2}}/T+{\mathbf{\alpha }}_{\mathbf{3}}·log\left(T\right)+{\mathbf{\alpha }}_{\mathbf{4}}·T+{\mathbf{\alpha }}_{\mathbf{5}}·T^2,\hfill \end{array}$$

(7)

and $k_j$ are temperature dependent reaction rate constants for reactions $j\in [\mathrm{CO},{\mathrm{CO}}_2,\mathrm{WGS}]$,

$$\begin{array}{c}\hfill k_j=exp\left({\mathrm{A}}_j-{\mathrm{B}}_j·\left({\displaystyle \frac{{\mathrm{T}}_{\mathrm{ref}}}{T}}-1\right)\right),\end{array}$$

(8)

with ${\mathrm{T}}_{\mathrm{ref}}=523.15$ K. The relative number of available catalyst centers can be calculated with

$${\Theta }^{\left(\mathrm{oxi}\right)}={\left(1+{\mathsf{\beta }}_{\mathrm{CO}}^{\left(\mathrm{oxi}\right)}·p_{\mathrm{CO}}\right)}^{-1},$$

(9)

$${\Theta }^{\left(\mathrm{red}\right)}={\left(1+{\mathsf{\beta }}_{{\mathrm{H}}^2}^{\left(\mathrm{red}\right)}·p_{{\mathrm{H}}^2}^{0.5}\right)}^{-1},$$

(10)

$${\Theta }^{\left(\mathrm{het}\right)}={\left(1+{\mathsf{\beta }}_{{\mathrm{CO}}^2}^{\left(\mathrm{het}\right)}·p_{{\mathrm{CO}}^2}\right)}^{-1}.$$

(11)

The parameters $\mathbf{\beta }$ were fit to experimental data by Seidel et al. [27]. Note that terms with parameters equal to 0 in [27] are neglected here. For a refit of the kinetics, the full version would be required.

Finally, the differential equation to describe the oxidization state of the catalyst reads

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\varphi }{\mathrm{d}t}& ={\mathrm{k}}_1^{+}·\left(y_{\mathrm{CO}}·\Delta \varphi -{\gamma }_1^{-1}·y_{{\mathrm{CO}}_2}·\varphi \right)+{\mathrm{k}}_2^{+}·\left(y_{{\mathrm{H}}_2}·\Delta \varphi -{\gamma }_2^{-1}·y_{{\mathrm{H}}_2\mathrm{O}}·\varphi \right),\hfill \end{array}$$

(12)

where $y_i$ are the mole fractions of component i, $\Delta \varphi ={\mathsf{\varphi }}_{\max }-\varphi$, and

$$\begin{array}{c}\mathbf{\gamma }=exp\left({\displaystyle \frac{-\Delta \mathrm{G}}{\mathrm{R}·T}}\right).\end{array}$$

(13)

The kinetics were fitted to steady-state and dynamical experiments conducted in a micro-Berty reactor. This reactor type resembles a good approximation of an ideally mixed continuous stirred tank reactor (CSTR).

In contrast to traditional kinetic models for methanol synthesis [28,29,30], dynamic changes of the catalyst are taken into account (Equation (12)), therefore the model is also able to reproduce the transient behavior [25].

2.2. Reactor Model

We consider a non-isothermal CSTR under constant pressure. In the following description of the reactor model, the species are denoted by subscripts i and k, the reactions are denoted by subscript j. The gas and solid phases are denoted by superscript G and S, respectively. The derivation of the material balances follows that of Seidel et al. [27] closely. However, note that in the present work the temperature also changes with time, hence the equations have to be adopted accordingly.

We start with the component material balance,

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\left(n_i^{\mathrm{G}}+n_i^{\mathrm{S}}\right)={\dot{n}}_{\mathrm{in}}·y_{i,\mathrm{in}}-{\dot{n}}_{\mathrm{out}}·y_i+m_{\mathrm{cat}}·\sum _j{\nu }_{i,j}·r_j,\end{array}$$

(14)

with,

$$\begin{array}{c}{n}_i^{\mathrm{G}}=n^{\mathrm{G}}·y_i,\hfill \end{array}$$

(15)

$$\begin{array}{c}{n}^{\mathrm{G}}={\displaystyle \frac{\mathrm{p}·\mathrm{V}}{\mathrm{R}·T}},{\dot{n}}_{\mathrm{out}}={\displaystyle \frac{\mathrm{p}·{\dot{V}}_{\mathrm{out}}}{\mathrm{R}·T}},{\dot{n}}_{\mathrm{in}}={\displaystyle \frac{\mathrm{p}·{\dot{V}}_{\mathrm{in}}}{\mathrm{R}·T_{\mathrm{in}}}},\hfill \end{array}$$

(16)

$$\begin{array}{c}{n}_i^{\mathrm{S}}={\mathrm{m}}_{\mathrm{cat}}·{\mathrm{q}}_{\mathrm{sat}}·{\Theta }_i,\hfill \end{array}$$

(17)

$$\begin{array}{c}{\Theta }_i=p_i·\sum _l{\mathsf{\beta }}_i^l·{\Theta }^l,l\in [\left(\mathrm{oxi}\right),\left(\mathrm{red}\right),\left(\mathrm{het}\right)],\hfill \end{array}$$

(18)

where $\mathit{y}$ are the molar fractions, $\mathbf{\nu }$ are the stoichiometric coefficients, the reactor’s pressure, volume, and temperature are given by $\mathrm{p}$, $\mathrm{V}$, and T, respectively, ${\dot{V}}_{\mathrm{out}}$ is the volumetric outlet flow-rate, ${\dot{V}}_{\mathrm{in}}$ is the feed flow rate, $T_{\mathrm{in}}$ is the feed temperature, ${\mathrm{m}}_{\mathrm{cat}}$ is the catalyst mass, ${\mathrm{q}}_{\mathrm{sat}}$ is the adsorption capacity of the catalyst, and $\Theta$ are competitive Langmuir adsorption isotherms, which depend on the partial pressures of the species. Details are given in [31]. The reactor temperature T is equal to the outlet temperature, hence the subscript “out” will be neglected for better readability. Note that all ${\mathbf{\beta }}_{\mathit{i}}$ in Equation (18) not used in Equations (9)–(11) are equal to 0.

The time derivatives for $n_i^{\mathrm{G}}$, and $n_i^{\mathrm{S}}$ are given by

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{n}_i^{\mathrm{S}}}{\mathrm{d}t}& ={\mathrm{m}}_{\mathrm{cat}}·{\mathrm{q}}_{\mathrm{sat}}·\sum _k\frac{\partial {\Theta }_i}{\partial p_k}·\frac{\mathrm{d}{p}_k}{\mathrm{d}t},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{n}_i^{\mathrm{G}}}{\mathrm{d}t}& =\frac{\partial n^{\mathrm{G}}}{\partial T}·\frac{\mathrm{d}T}{\mathrm{d}t}·y_i+n^{\mathrm{G}}·\frac{\mathrm{d}{y}_i}{\mathrm{d}t};\hfill \end{array}$$

(20)

an explicit calculation of the partial derivatives of the adsorption isotherms, $\frac{\partial {\Theta }_i}{\partial p_k}$, can be found in the Appendix of Nikolić et al. [31].

Introducing Equations (19) and (20) into Equation (14),

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{\partial n^{\mathrm{G}}}{\partial T}·\frac{\mathrm{d}T}{\mathrm{d}t}·y_i+n^{\mathrm{G}}·\frac{\mathrm{d}{y}_i}{\mathrm{d}t}+{\mathrm{m}}_{\mathrm{cat}}·q_{\mathrm{sat}}·\mathrm{p}·\sum _k\frac{\partial {\Theta }_i}{\partial p_k}·\frac{\mathrm{d}{y}_k}{\mathrm{d}t}=\hfill \\ \hfill & {\dot{n}}_{\mathrm{in}}·y_{i,\mathrm{in}}-{\dot{n}}_{\mathrm{out}}·y_i+{\mathrm{m}}_{\mathrm{cat}}·\sum _j{\nu }_{i,j}·r_j,\hfill \end{array}\end{array}$$

(21)

and summation over the single species yields the total material balance,

$$\begin{array}{c}\hfill \frac{\partial n^{\mathrm{G}}}{\partial T}·\frac{\mathrm{d}T}{\mathrm{d}t}+{\mathrm{m}}_{\mathrm{cat}}·{\mathrm{q}}_{\mathrm{sat}}·\mathrm{P}·\sum _i\sum _k\frac{\partial {\Theta }_i}{\partial p_k}·\frac{\mathrm{d}{y}_k}{\mathrm{d}t}={\dot{n}}_{\mathrm{in}}-{\dot{n}}_{\mathrm{out}}+{\mathrm{m}}_{\mathrm{cat}}·\sum _i\sum _j{\nu }_{i,j}·r_j.\end{array}$$

(22)

Solving Equation (22) for ${\dot{n}}_{\mathrm{out}}$ to eliminate it in Equation (21) yields the component material balance in its final form,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & n^{\mathrm{G}}·\frac{\mathrm{d}{y}_i}{\mathrm{d}t}+{\mathrm{m}}_{\mathrm{cat}}·{\mathrm{q}}_{\mathrm{sat}}·\mathrm{p}·\left(\sum _k\frac{\partial {\Theta }_i}{\partial p_k}·\frac{\mathrm{d}{y}_k}{\mathrm{d}t}-\sum _i\sum _k\left(\frac{\partial {\Theta }_i}{\partial p_k}·\frac{\mathrm{d}{y}_k}{\mathrm{d}t}\right)·y_i\right)=\hfill \\ \hfill & {\dot{n}}_{\mathrm{in}}·\left(y_{i,\mathrm{in}}-y_i\right)+{\mathrm{m}}_{\mathrm{cat}}·\left(\sum _j{\nu }_{i,j}·r_j-\sum _i\sum _j\left({\nu }_{i,j}·r_j\right)·y_i\right).\hfill \end{array}\end{array}$$

(23)

Note that Equation (23), in contrast to Equations (21) and (22), is identical to that obtained assuming a constant reactor temperature [27], i.e., the component material balance is not temperature-dependent. However, the volumetric outlet flow rate (which can be derived from Equation (22)) differs, as it is temperature dependent,

$$\begin{array}{c}\hfill {\dot{V}}_{\mathrm{out}}={\displaystyle \frac{\mathrm{V}}{T}}·\frac{\mathrm{d}T}{\mathrm{d}t}+{\displaystyle \frac{T}{T_{\mathrm{in}}}}{\dot{V}}_{\mathrm{in}}+{\displaystyle \frac{{\mathrm{m}}_{\mathrm{cat}}·\mathrm{R}·T}{\mathrm{p}}}·\left(\sum _i\sum _j{\nu }_{i,j}·r_j-{\mathrm{q}}_{\mathrm{sat}}·\mathrm{p}·\sum _i\sum _k\frac{\partial {\Theta }_i}{\partial p_k}·\frac{\mathrm{d}{y}_k}{\mathrm{d}t}\right).\end{array}$$

(24)

2.3. Energy Balance

To describe the change in temperature, an energy balance is required. We start with the total energy balance,

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\left(\sum _i{n}_i^{\mathrm{G}}·h_i^{\mathrm{G}}+\sum _i{n}_i^{\mathrm{S}}·h_i^{\mathrm{S}}+{\mathrm{m}}_{\mathrm{cat}}·h_{\mathrm{cat}}^{*}\right)={\dot{n}}_{\mathrm{in}}·\sum _i{h}_{i,\mathrm{in}}^{\mathrm{G}}·y_{i,\mathrm{in}}-{\dot{n}}_{\mathrm{out}}·\sum _i{h}_i^{\mathrm{G}}·y_i-{\dot{Q}}_{\mathrm{c}},\end{array}$$

(25)

where $\mathit{h}$ are enthalpies, and ${\dot{Q}}_{\mathrm{c}}$ is the cooling. Applying the product rule, introducing $\mathrm{d}{h}_{\mathrm{cat}}^{\ast }={\mathrm{c}}_{\mathrm{p}}^{\mathrm{cat}}\mathrm{d}T$, and assuming constant ${\mathrm{m}}_{\mathrm{cat}}$, and ${\mathrm{c}}_{\mathrm{p}}^{\mathrm{cat}}$ yields

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \sum _i\frac{\mathrm{d}{n}_i^{\mathrm{G}}}{\mathrm{d}t}·h_i^{\mathrm{G}}+\sum _i{n}_i^{\mathrm{G}}·\frac{\mathrm{d}{h}_i^{\mathrm{G}}}{\mathrm{d}t}+\sum _i\frac{\mathrm{d}{n}_i^{\mathrm{S}}}{\mathrm{d}t}·h_i^{\mathrm{S}}+\sum _i{n}_i^{\mathrm{S}}·\frac{\mathrm{d}{h}_i^{\mathrm{S}}}{\mathrm{d}t}+{\mathrm{m}}_{\mathrm{cat}}·{\mathrm{c}}_{\mathrm{p}}^{\mathrm{cat}}\frac{\mathrm{d}T}{\mathrm{d}t}=& \\ \hfill {\dot{n}}_{\mathrm{in}}·\sum _i{h}_{i,\mathrm{in}}·y_{i,\mathrm{in}}-{\dot{n}}_{\mathrm{out}}·\sum _i{h}_i^{\mathrm{G}}·y_i-{\dot{Q}}_{\mathrm{c}}.\end{array}\end{array}$$

(26)

To simplify the equation, the component material balance, Equation (14), is multiplied with the component gas enthalpies, $h_i^{\mathrm{G}}$, summed over i, and solved for

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \sum _i\frac{\mathrm{d}{n}_i^{\mathrm{G}}}{\mathrm{d}t}·h_i^{\mathrm{G}}=& -\sum _i\frac{\mathrm{d}{n}_i^{\mathrm{S}}}{\mathrm{d}t}·h_i^{\mathrm{G}}+{\dot{n}}_{\mathrm{in}}·\sum _i{h}_i^{\mathrm{G}}·y_{i,\mathrm{in}}\hfill \\ \hfill & -{\dot{n}}_{\mathrm{out}}·\sum _i{h}_i^{\mathrm{G}}·y_i+{\mathrm{m}}_{\mathrm{cat}}\sum _j\Delta h_{\mathrm{R},j}·r_j,\hfill \end{array}\end{array}$$

(27)

where $\Delta h_{\mathrm{R},j}$ is the heat of reaction j. Putting this into Equation (26) yields

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \sum _i{n}_i^{\mathrm{G}}·\frac{\mathrm{d}{h}_i^{\mathrm{G}}}{\mathrm{d}t}+\sum _i\frac{\mathrm{d}{n}_i^{\mathrm{S}}}{\mathrm{d}t}·\left(h_i^{\mathrm{S}}-h_i^{\mathrm{G}}\right)+\sum _i{n}_i^{\mathrm{S}}·\frac{\mathrm{d}{h}_i^{\mathrm{S}}}{\mathrm{d}t}+{\mathrm{m}}_{\mathrm{cat}}·c_{\mathrm{p}}^{\mathrm{cat}}\frac{\mathrm{d}T}{\mathrm{d}t}=\hfill \\ \hfill & {\dot{n}}_{\mathrm{in}}·\sum _i\left(h_{i,\mathrm{in}}^{\mathrm{G}}-h_i^{\mathrm{G}}\right)·y_{i,\mathrm{in}}-{\mathrm{m}}_{\mathrm{cat}}\sum _j\Delta h_{\mathrm{R},j}·r_j-{\dot{Q}}_{\mathrm{c}}.\hfill \end{array}\end{array}$$

(28)

The enthalpies can be described by

$$\begin{array}{c}\hfill h_i={\int }_{{\mathrm{T}}_{\mathrm{ref}}}^T{\mathit{c}}_{\mathrm{p},i}\left(T\right)\mathrm{d}T+\Delta {\mathrm{h}}_{\mathrm{f},\mathrm{ref},i},\end{array}$$

(29)

where $\Delta {\mathrm{h}}_{\mathrm{f},\mathrm{ref},i}$ is the heat of formation and ${\mathit{c}}_{\mathrm{p}}$ is the isobaric heat capacity. We assume that ${\mathit{c}}_{\mathrm{p}}^{\mathrm{S}}$ does not depend on the temperature, a correlation for ${\mathit{c}}_{\mathrm{p}}^{\mathrm{G}}$ in J mol${}^{-1}$ K${}^{-1}$,

$$\begin{array}{cc}\hfill c_{\mathrm{p}},{\mathrm{i}}^{\mathrm{G}}& ={\Gamma }_{1,i}+{\Gamma }_{2,i}·T+{\Gamma }_{3,i}·T^2+{\Gamma }_{4,i}·T^3,\hfill \end{array}$$

(30)

as well as values for $\Delta {\mathbf{h}}_{\mathrm{f},\mathrm{ref}}$ and $\mathbf{\Gamma }$, given in

Table 2. Thermodynamic data, taken from Appendix A of Reid et al. [32]; $\Delta {\mathrm{h}}_{\mathrm{f},\mathrm{ref},i}$ in J mol${}^{-1}$; reference temperature for $\Delta {\mathrm{h}}_{\mathrm{f},\mathrm{ref},i}$: ${\mathrm{T}}_{\mathrm{ref}}=298.2$ K.

Species	${\mathbf{\Gamma }}_{1,\mathit{i}}$	${\mathbf{\Gamma }}_{2,\mathit{i}}$	${\mathbf{\Gamma }}_{3,\mathit{i}}$	${\mathbf{\Gamma }}_{4,\mathit{i}}$	$\Delta {\mathit{h}}_{\mathbf{f},\mathbf{ref},\mathit{i}}$
CH3OH	$2.115·{10}^1$	$7.092·{10}^{-2}$	$2.587·{10}^{-5}$	$-2.852·{10}^{-8}$	$-2.013·{10}^5$
CO2	$1.980·{10}^1$	$7.344·{10}^{-2}$	$-5.602·{10}^{-5}$	$1.715·{10}^{-8}$	$-3.938·{10}^5$
CO	$3.087·{10}^1$	$-1.285·{10}^{-2}$	$2.789·{10}^{-5}$	$-1.272·{10}^{-8}$	$-1.106·{10}^5$
H2	$2.714·{10}^1$	$9.274·{10}^{-3}$	$-1.381·{10}^{-5}$	$7.645·{10}^{-9}$	0
H2O	$3.224·{10}^1$	$1.924·{10}^{-3}$	$1.055·{10}^{-5}$	$-3.596·{10}^{-9}$	$-2.420·{10}^5$

, can be found in [32].

With that, and by introducing the heat of cooling, the following equation can be obtained for the change in temperature:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \left(\underset{\text{gaseous species}}{\underbrace{n^{\mathrm{G}}·{\overline{c}}_{\mathrm{p}}^{\mathrm{G}}}}+{\underset{\text{adsorbed species}}{\underbrace{\overline{){\mathrm{m}}_{\mathrm{cat}}·q_{\mathrm{sat}}·{\overline{c}}_{\mathrm{p}}^{\mathrm{S}}}}}}^0+\underset{“\mathrm{dead}” \mathrm{catalyst}}{\underbrace{{\mathrm{m}}_{\mathrm{cat}}·c_{\mathrm{p}}^{\mathrm{cat}}}}\right)·\frac{\mathrm{d}T}{\mathrm{d}t}=\hfill \\ \hfill & -\underset{\text{heat of transport}}{\underbrace{{\dot{n}}_{\mathrm{in}}·{\sum }_i{\int }_{T_{\mathrm{in}}}^T{c}_{\mathrm{p},i}^{\mathrm{G}}\mathrm{d}T·y_{i,\mathrm{in}}}}-\underset{\text{heat of reaction}}{\underbrace{{\mathrm{m}}_{\mathrm{cat}}{\sum }_j\Delta h_{\mathrm{R},j}·r_j}}-\underset{\text{heat of cooling}}{\underbrace{{\mathrm{K}}_{\mathrm{W}}·{\mathrm{A}}_{\mathrm{W}}·\left(T-T_{\mathrm{c}}\right)}}-\underset{\text{heat of adsorption}}{\underbrace{{\sum }_i\frac{\mathrm{d}{n}_i^{\mathrm{S}}}{\mathrm{d}t}·\Delta {\mathrm{h}}_{\mathrm{ads},i}}},\hfill \end{array}\end{array}$$

(31)

with

$$\begin{array}{c}\hfill {\overline{c}}_{\mathrm{p}}^{\mathrm{G}}=\sum _i{c}_{\mathrm{p},i}^{\mathrm{G}}·y_i,\Delta {\mathrm{h}}_{\mathrm{ads},i}=h_i^{\mathrm{S}}-h_i^{\mathrm{G}}.\end{array}$$

(32)

The adsorption enthalpy is assumed to be constant and identical for each component, with $\Delta {\mathrm{h}}_{\mathrm{ads}}=1·{10}^3$ J mol${}^{-1}$. Note that, since the heat of adsorption is relatively small compared to the heat of reaction, only a rough estimate is used.

2.4. Multistage Reactor

The considered process is depicted in Figure 1. The multistage reactor consists of three CSTR with side-feeds for each reactor, which can be seen as an approximation of a staged multibed reactor. For CO₂ methanation, it is demonstrated in [19,20] that such a feed split increases robustness against fluctuations. The idea is transferred here to methanol production from CO. A different concept of a staged fixed bed reactor for methanol production from CO₂ was studied by Torcida et al. [33]. There, intermediate cooling without feed split between the different stages is used to increase reactor efficiency.

We consider 3 reactors within the setup; they will be denoted by the subscript $r\in [1,2,3]$ in the following.

The feed is split by a split-ratio, $\mathit{\kappa }$, with

$$\begin{array}{c}{\dot{n}}_{\mathrm{in},1}={\kappa }_1·{\dot{n}}_{\mathrm{in},\mathrm{tot}},\hfill \end{array}$$

(33)

$$\begin{array}{c}{\dot{n}}_{\mathrm{in},\overline{r}}={\kappa }_{\overline{r}}·{\dot{n}}_{\mathrm{in},\mathrm{tot}}+{\dot{n}}_{\mathrm{out},\overline{r}-1},\overline{r}\in [2,3],\hfill \end{array}$$

(34)

$$\begin{array}{c}\sum _r{\kappa }_r=1,\hfill \end{array}$$

(35)

where ${\dot{n}}_{\mathrm{in},\mathrm{tot}}$ is the total volumetric inlet flow into the system. It consists of the fluctuating hydrogen stream, $\xi ·{\dot{n}}_{{\mathrm{H}}_2,in,tot}$, and a carbon stream, ${\dot{n}}_{\mathrm{C},in,tot}$, its composition and amount is calculated with

$$\begin{array}{c}{\dot{n}}_{\mathrm{in},\mathrm{tot}}=\xi ·{\dot{n}}_{{\mathrm{H}}_2,\mathrm{in},\mathrm{tot}}+{\dot{n}}_{\mathrm{C},\mathrm{in},\mathrm{tot}},\hfill \end{array}$$

(36)

$$\begin{array}{c}{y}_{{\mathrm{H}}_2,\mathrm{in},\mathrm{tot}}={\displaystyle \frac{\xi ·{\dot{n}}_{{\mathrm{H}}_2,\mathrm{in},\mathrm{tot}}}{{\dot{n}}_{\mathrm{in},\mathrm{tot}}}},\hfill \end{array}$$

(37)

$$\begin{array}{c}{y}_{{\mathrm{CO}}_2,\mathrm{in},\mathrm{tot}}={\displaystyle \frac{{\dot{n}}_{\mathrm{C},\mathrm{in},\mathrm{tot}}·{\psi }_{\mathrm{C}}}{{\dot{n}}_{\mathrm{in},\mathrm{tot}}}},\hfill \end{array}$$

(38)

$$\begin{array}{c}{y}_{\mathrm{CO},\mathrm{in},\mathrm{tot}}={\displaystyle \frac{{\dot{n}}_{\mathrm{C},\mathrm{in},\mathrm{tot}}·\left(1-{\psi }_{\mathrm{C}}\right)}{{\dot{n}}_{\mathrm{in},\mathrm{tot}}}},\hfill \end{array}$$

(39)

The inlet compositions and temperatures for reactors 2 and 3 are calculated accordingly with

$$\begin{array}{c}{\mathit{y}}_{\mathrm{in},\overline{r}}={\displaystyle \frac{{\mathit{y}}_{\overline{r}-1}·{\dot{n}}_{\mathrm{out},\overline{r}-1}+{\kappa }_{\overline{r}}·{\dot{n}}_{\mathrm{in},\mathrm{tot}}·{\mathit{y}}_{\mathrm{in},\mathrm{tot}}}{{\dot{n}}_{\mathrm{out},\overline{r}-1}+{\kappa }_{\overline{r}}·{\dot{n}}_{\mathrm{in},\mathrm{tot}}}},\hfill \end{array}$$

(40)

$$\begin{array}{c}{T}_{\mathrm{in},\overline{r}}={\displaystyle \frac{T_{\overline{r}-1}·{\dot{n}}_{\mathrm{out},\overline{r}-1}+{\kappa }_{\overline{r}}·{\dot{n}}_{\mathrm{in},\mathrm{tot}}·T_{\mathrm{in},\mathrm{tot}}}{{\dot{n}}_{\mathrm{out},\overline{r}-1}+{\kappa }_{\overline{r}}·{\dot{n}}_{\mathrm{in},\mathrm{tot}}}},\hfill \end{array}$$

(41)

where $T_{\mathrm{in}},\mathrm{tot}$ is the temperature of the inlet into the system. Note that the inlet temperature and inlet composition for stage $\overline{r}$ are mean values based on the mole amounts of the entering streams, hence those values are estimates.

The reactors are parameterized such that the overall system resembles an industrial scale reactor. The heat transfer coefficient, ${\mathrm{K}}_{\mathrm{W}}$, was approximated according to [34,35]. The parameters are given in

Table 3. Reactor parameters.

Parameter	Value	Unit
${\mathrm{V}}_{\mathrm{tot}}$	$2.83$	m${}^3$
P	70	bar
${\mathrm{K}}_{\mathrm{W}}$	250	W m${}^{-2}$ K ${}^{-1}$
${\mathrm{A}}_{\mathrm{W}}$	18.85	m${}^2$
${\mathrm{q}}_{\mathrm{sat}}$	0.98	mol kg${}^{-1}$
${\mathrm{c}}_{\mathrm{p}}^{\mathrm{cat}}$	1063	J kg${}^{-1}$ K${}^{-1}$
${\mathrm{m}}_{\mathrm{cat},\mathrm{tot}}$	2002	kg
${\mathsf{\rho }}_{\mathrm{cat}}$	1770	kg m${}^{-3}$

. The overall reactor volume and catalyst mass are split equally,

$$\begin{array}{cc}\hfill {\mathrm{V}}_r& ={\displaystyle \frac{{\mathrm{V}}_{\mathrm{tot}}}{3}},\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill m_{\mathrm{cat},r}& ={\displaystyle \frac{{\mathrm{m}}_{\mathrm{cat},\mathrm{tot}}}{3}}.\hfill \end{array}$$

(43)

2.5. Optimization Approach

We consider three optimization problems. The first one is the nominal steady-state optimization problem,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{\mathit{x}}{max}& J=STY{\left(\mathit{x}\right)}^2,\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \text{Equations (4)–(43) in steady-state},\hfill \\ \hfill & {\mathit{x}}^{\mathrm{lo}}\le \mathit{x}\le {\mathit{x}}^{\mathrm{up}},\hfill \\ \hfill & X_{\mathrm{C}}\left(\mathit{x}\right)\ge 0.6,\hfill \\ \hfill & {\left(\mathit{T}\left(\mathit{x}\right)-{\mathit{T}}_c\right)}^2\le {30}^2,\hfill \\ \hfill & {\displaystyle \sum _i}{y}_i\left(\mathit{x}\right)=1,\hfill \\ \hfill & {\displaystyle \sum _i}{y}_{i,\mathrm{in},\mathrm{tot}}=1,\hfill \\ \hfill & y_{\mathrm{CO},\mathrm{in},\mathrm{tot}}+y_{{\mathrm{CO}}_2,\mathrm{in},\mathrm{tot}}\ge 0.01,\hfill \end{array}\end{array}$$

(44)

where $\mathit{x}\in G\subseteq {\mathbb{R}}^d$ are variables with upper (${\mathit{x}}^{\mathrm{up}}$) and lower bounds (${\mathit{x}}^{\mathrm{lo}}$). The space–time yield is defined as

$$\begin{array}{c}STY=\left({\displaystyle \frac{{\dot{n}}_{\mathrm{out},4}·y_{{\mathrm{CH}}_3\mathrm{OH},4}}{{\mathrm{V}}_{\mathrm{tot}}}}\right),\hfill \end{array}$$

(45)

and the carbon conversion is defined as

$$X_{\mathrm{C}}={\displaystyle \frac{y_{{\mathrm{CH}}_3\mathrm{OH},4}·{\dot{n}}_{\mathrm{out},4}}{{\dot{n}}_{\mathrm{C},\mathrm{in},\mathrm{tot}}}},$$

(46)

upper and lower bounds for the variables are given in

Table 4. Variables and respective bounds. Variables below the mid rule are degrees of freedom. Degrees of freedom marked with “*” may change during the dynamic optimization.

Variable	Domain	Unit	Description
$\mathit{T}$	$[480,573]$	K	Reactor temperature
${\dot{\mathit{V}}}_{\mathrm{out}}$	$[1·{10}^{-4},5·{10}^{-1}]$	m${}^3$ s${}^{-1}$	Volumetric outlet flowrate
$y_{\mathrm{CO},\mathrm{in}}$	$[0,0.5]$	-	CO inlet mole fraction
$y_{{\mathrm{CO}}_2,\mathrm{in}}$	$[0,0.5]$	-	CO2 inlet mole fraction
$y_{{\mathrm{H}}_2,\mathrm{in}}$	$[0.5,0.85]$	-	H2 inlet mole fraction
${\mathit{T}}_{\mathrm{c}}^{\ast }$	$[490,550]$	K	Cooling temperature
${\mathit{T}}_{\mathrm{c},\mathrm{s}}^{\ast }$	$[490,550]$	K	Cooling temperature setpoint
$T_{\mathrm{in},\mathrm{tot}}$	$[490,520]$	K	System feed temperature
${\dot{n}}_{\mathrm{C},\mathrm{in},\mathrm{tot}}$	$[1,300]$	mol s${}^{-1}$	System inlet carbon molar flow
${\dot{n}}_{{\mathrm{H}}_2,\mathrm{in},\mathrm{tot}}$	$[1,300]$	mol s${}^{-1}$	System inlet hydrogen molar flow
${\psi }_{\mathrm{C}}$	$[0,1]$	-	CO, CO2 split ratio
${\mathit{\kappa }}^{\ast }$	$[0,1]$	-	Feed split ratio
${\mathit{\kappa }}_s^{\ast }$	$[0,1]$	-	Feed split ratio setpoint
$s_1^{\ast }$	$[0,0.3]$	-	Carbon conversion slack variable
$s_2^{\ast }$	$[0,10]$	-	Cooling temperature slack variable

. We enforce a specific carbon conversion rate, while carbon conversion and the objective of a maximized space–time yield are concurrent. This also leads to higher reactor temperatures, as the reaction is exothermic. To prevent hotspot formation, there is a restriction on the differences between the reactor temperatures and the shell temperatures in Problem 1 (44) [36].

The second problem is the robust steady-state optimization,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{\mathit{x}}{max}& J={\mathbb{E}}_{\xi }\left(STY{(\mathit{x},\xi )}^2\right),\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \text{Equations (4)–(43) in steady-state},\hfill \\ \hfill & {\mathit{x}}^{\mathrm{lo}}\le \mathit{x}\le {\mathit{x}}^{\mathrm{up}},\hfill \\ \hfill & X_{\mathrm{C}}(\mathit{x},\xi )\ge 0.6,\hfill \\ \hfill & {\left(\mathit{T}(\mathit{x},\xi )-{\mathit{T}}_c\right)}^2\le {30}^2,\hfill \\ \hfill & {\displaystyle \sum _i}{y}_i(\mathit{x},\xi )=1,\hfill \\ \hfill & {\displaystyle \sum _i}{y}_{i,\mathrm{in},\mathrm{tot}}=1,\hfill \\ \hfill & y_{\mathrm{CO},\mathrm{in},\mathrm{tot}}+y_{{\mathrm{CO}}_2,\mathrm{in},\mathrm{tot}}\ge 0.01,\hfill \end{array}\end{array}$$

(47)

where the space–time yield is maximized with the same constraints as in Problem 1 (44), but with taking fluctuations into account (see the left plot of Figure 2): $\xi \in [0.75,1.25]$ is the disturbance of H₂ at the inlet, the problem is discretized along $\xi$ using 50 equally distant distributed support points. It is assumed that the carbon source is able to deliver CO and CO₂ as needed, but the hydrogen source is prone to fluctuations due to its renewable character, i.e., hydrogen is produced via electrolysis using wind energy. Each support point of $\xi$ refers to a different scenario with a given disturbance rate. The problem is solved simultaneously, i.e., the optimizer identifies an optimal solution that holds the constraints for each support point of $\xi$.

Note that we only consider steady-state scenarios of the process here. Transient behavior that may occur when changing from one disturbance to another is neglected in this step.

The third optimization problem is a dynamic optimization, where the design parameters of the reactor setup are fixed. Only operating conditions, i.e., the feed split and the shell temperatures, may be changed dynamically. The problem is defined as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{\mathit{x}}{max}& J=\int STY{\left(\mathit{x},\xi \right)}^2/t_{\mathrm{end}}^2-\mathrm{U}·(s_1+s_2)\mathrm{d}t\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathrm{Equations}\left(4\right)\text{–}\left(43\right),\hfill \\ \hfill & {\mathit{x}}^{\mathrm{lo}}\le \mathit{x}\le {\mathit{x}}^{\mathrm{up}},\hfill \\ \hfill & X_{\mathrm{C}}(\mathit{x},\xi )+s_1\ge 0.6,\hfill \\ \hfill & {\left(\mathit{T}(\mathit{x},\xi )-{\mathit{T}}_c\right)}^2\le {(30+s_2)}^2,\hfill \\ \hfill & {\displaystyle \sum _i}{y}_i(\mathit{x},\xi )=1,\hfill \\ \hfill & {\displaystyle \sum _i}{y}_{i,\mathrm{in}}=1,\hfill \\ \hfill & {\dot{\mathit{T}}}_c={\displaystyle \frac{-{\mathit{T}}_c+{\mathit{T}}_{c,s}}{50}},\hfill \\ \hfill & \dot{\mathit{\kappa }}={\displaystyle \frac{-\mathit{\kappa }+{\mathit{\kappa }}_s}{15}},\hfill \end{array}\end{array}$$

(48)

where $\mathit{s}$ are slack variables, and $\mathrm{U}=100$ is a heuristically set penalty constant. The slack variables allow for violation of the constraints at the expense of a highly reduced cost function, hence the optimizer tries to avoid this as much as possible.

Note that additional constraints are present in Problem 3 (48) to determine the dynamic behavior of the feed split and shell temperature. First order lag equations are used to obtain a more realistic result, with time constants of 15 and 50, respectively. Without those lags, the optimizer would be able to change the operating conditions instantly from one extreme to another in the feasible range, giving rise to so-called bang-bang solutions. Although bang-bang solutions may be optimal, they are not achievable in a real application. Variables ${\mathit{T}}_{\mathrm{c},\mathrm{s}}$, and ${\mathit{\kappa }}_{\mathrm{c},\mathrm{s}}$ are the setpoints of the lag elements, which are degrees of freedom for the optimizer.

All optimization problems are implemented in Julia [37], using the JuMP [38] based package InfiniteOpt [39]. For Problem 3 (48) an orthogonal collocation of third order [40] with 200 equally distant support points is used for the discretization of the infinite dimensional problem. All three problems are solved using the NLP solver KNITRO [41]; for Problems 1 (44) and 2 (47), the multi-start option is set to 50.

2.6. Test Function

A drawback of Problem 2 (47) is that only steady-state realizations are investigated, neglecting the transient behavior when fluctuations occur. The design obtained from Problem 1 (44) does not take fluctuations into account. Therefore, an additional test is required to check whether the imposed constraints are met under dynamically changing disturbances. The employed test function is depicted in the right subplot of Figure 2. It can be described by the following set of equations,

$$\begin{array}{cc}\hfill \Delta t& =t+226,t\in [0,700],\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill f& =\Re (exp(0.01\mathrm{i}·\Delta t)+exp(0.001\mathrm{i}·\Delta t)+exp(0.008\mathrm{i}·\Delta t\left)\right),\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill f_{\mathrm{norm}}& ={\displaystyle \frac{f+|min(f\left)\right|}{max(f+|min\left(f\right)\left)\right|}},\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill f_{\xi }& =\left|0.9+f_{\mathrm{norm}}·\left(0.3\right)\right|,\hfill \end{array}$$

(52)

where $\mathrm{i}$ is the imaginary number, and $\Re \left(f\right)$ is the real part of f. The time is shifted such that $f_{\xi }(t=0)=1$. The test function can be interpreted as a complex Fourier series, which is truncated after the third entry and normalized such that the resulting function lies within the desired interval.

Note that the employed test function for dynamic transient behavior does not resemble the worst-case scenario, which would be a series of rectangle steps. Instead, it is modeled such that it resembles real wind fluctuations [42,43] better. However, the presented disturbance fluctuates more rapidly than it would be observed in the real world. Hence, it is a model profile not related to an existing green electricity profile. A short time horizon of this modeled disturbance was arbitrarily chosen as the test function, marked with a black rectangle in the left plot of Figure 2.

The test function is used in a dynamic forward simulation of the obtained reactor setups from Problems 1 (44) and 2 (47), starting at the steady-state of the system.

A linear interpolation is used for the simulation of the design obtained from Problem 2 (47). The feed split ratios in between the support points are interpolated.

The simulation is implemented in Julia using the DifferentialEquations [44] package.

3. Results and Discussion

Two case studies are discussed in this section, one with the nominal reactor obtained from Problem 1 (44) and one with the robust reactor obtained from Problem 2 (47). Both reactor realizations are tested in a simulation framework and are compared to each other. In a last step, the dynamic optimization in Problem 3 (48) is solved for both designs.

3.1. Nominal Steady-State Design

The nominal design is the solution of Problem 1 (44). The obtained optimal values of the degrees of freedom are given in

Table 5. Obtained values for degrees of freedom for nominal design, solution of Problem 1 (44).

${\dot{\mathit{n}}}_{\mathbf{C},\mathbf{in},\mathbf{tot}}$	${\dot{\mathit{n}}}_{{\mathbf{H}}_2,\mathbf{in},\mathbf{tot}}$	${\mathit{\psi }}_{\mathbf{C}}$	${\mathit{T}}_{\mathbf{in},\mathbf{tot}}$	${\mathit{T}}_{\mathbf{c}}$	$\mathit{\kappa }$	$\mathit{STY}$
$3.67$	$20.81$	$0.0125$	490	$[496,537,528]$	$[0.195,0.805,0]$	$0.78$

.

The lowest possible inlet temperature is chosen. As the reaction is exothermic, the inlet stream acts as a sort of cooling mechanism.

The inlet stream consists mostly of hydrogen, with $y_{{\mathrm{H}}_2,\mathrm{in},\mathrm{tot}}=0.85$. Note that this value lies on the boundary. If no restriction on $y_{{\mathrm{H}}_2,\mathrm{in},\mathrm{tot}}$ is present, the value is chosen to be higher, which would be favorable for the temperature management of the reactor setup. However, most of the increased hydrogen stream would leave the reactor unreacted, increasing the overall process cost.

Regarding the feed split, approximately one fifth of the feed is fed into the first reactor and four fifths are fed into the second. The first reactor is relatively cold, which leads to a fast conversion and a high temperature increase; therefore, only a small fraction of the reactants can enter. With increasing temperature, the reaction and heat production slow down. Therefore, most of the feed is fed into the much hotter second reactor, which is not only cooled by the shell, but also by the entering streams. The feed to the third reactor is not used. Feeding the third reactor would reduce the residence time, as the first two reactors are bypassed. This would lead to the production of less methanol, which is unfavorable for the used objective function.

The dynamic simulation results for the nominal reactor design are shown in Figure 3. The figure consists of four subplots: the upper left corner shows the mole fraction of methanol over time, the upper right corner the reactor temperature over time, the lower left corner the molar outlet streams over time, and the lower right corner the carbon conversion rate over time. The fluctuations can be well observed in all of those subplots. It is clear from the lower right subplot that the nominal reactor is not feasible for the present disturbance, as it lies well below the desired carbon conversion of 60%. The obtained space–time yield is $STY=0.8047$.

3.2. Robust Design

The robust design is the solution of Problem 2 (47). The obtained values of the constant degrees of freedom are given in

Table 6. Obtained values for degrees of freedom for robust design, solution of Problem 2 (47).

${\dot{\mathit{n}}}_{\mathbf{C},\mathbf{in},\mathbf{tot}}$	${\dot{\mathit{n}}}_{{\mathbf{H}}_2,\mathbf{in},\mathbf{tot}}$	${\mathit{\psi }}_{\mathbf{C}}$	${\mathit{T}}_{\mathbf{in},\mathbf{tot}}$	${\mathit{T}}_{\mathbf{c}}$	${\mathbb{E}}_{\mathit{\xi }}\left(\mathit{STY}\right)$
$3.01$	$13.66$	$0.014$	490	$[543,526,519]$	$0.67$

, Figure 4 shows the feed split and the resulting carbon conversion rates. Note that the x-axis of Figure 4 is the disturbance.

The inlet composition is slightly different than for the nominal case. There is proportionally less H₂ fed into the reactor. However, the obtained value for ${\dot{n}}_{{\mathrm{H}}_2,\mathrm{in},\mathrm{tot}}$ leads to a mole fraction of 0.85 for the highest hydrogen inlet increase. There is also slightly more CO₂ present in the inlet stream.

The inlet temperature is, again, chosen to be as small as possible. There is a notable difference between the shell temperatures chosen for the nominal and robust designs. While the shell temperature of the first reactor in the nominal design is very low, it is the highest in the robust design. This can be explained by the feed split. In contrast to the nominal design, most of the feed is fed into the first reactor now.

The feed split changes for different disturbance values, although its main split remains constant: most of the feed goes into reactor 1. Interestingly, some of the feed is fed directly into reactor 3. Note that reactor 3 is the coldest and hence has the fastest reaction rate.

Generally, with increasing hydrogen flows, the steady-state carbon conversion increases and vice versa, as the temperature is reduced by a higher inlet flow, which increases the reaction rate due to the temperature dependence of the equilibrium of the water-gas-shift reaction.

Compared to the nominal design, the mean space–time yield of the robust design is reduced by 14%. This is to be expected, as the robust design has to be more conservative and hence inefficient.

The dynamic simulation results for the robust reactor design are shown in Figure 5 (the figure has the same layout as Figure 3, described above). Although much better than in the nominal case, the robust design is still infeasible with regard to the carbon conversion, as it lies below the desired 60%. This can be explained by the neglected transient behavior and non-linear dynamics of the model, as Problem 2 (47) only considers a steady-state optimization. The obtained space–time yield is $STY=0.7037$, which is 12.6% less than in the nominal case and the price we have to pay for the increased robustness.

3.3. Dynamic Optimization

The dynamic optimization is the solution of Problem 3 (48).

The result for the robust design is shown in Figure 6, where the collocation support points are plotted.

The feed split ratios are chosen differently than for the robust steady-state design. The feed bypass to the third reactor is completely shut down, while the second reactor gets more of the feed stream. For the first reactor, the value of the split is mostly close to that of the steady-state design, except for some drastic dynamic changes between 300 and 400 s, and 636 and 700 s. The first interval is located in the region where the simulation of the steady-state design became infeasible. Therefore, the optimizer has to choose a different control to achieve feasibility. The second interval comes due to the absence of terminal constraints, the objective function is increased by maximizing the output. In a practical application, the optimization would be repeated before this interval is reached.

The achieved space–time yield is $STY=0.6987$, which is slightly below that of the robust case. It is, however, feasible.

The dynamic optimization for the nominal design converges to an infeasible point. The optimizer is not able to control the process such that the reactor temperatures remain within their bounds.

4. Conclusions

Methanol synthesis is a crucial topic in light of the upcoming changes in energy supply. As traditional steady-state processes will have to cope with energy fluctuations, they will also have to deal with fluctuations in their operating conditions. In this paper, we have investigated the influence of fluctuations on the inlet flow of a staged methanol synthesis reactor.

To address this challenge, both a nominal and a robust steady-state optimization of the reactor setup’s operating conditions were conducted. As expected, the obtained optimal nominal reactor has a better performance than the robust reactor, which has to be more conservative in its design. However, for both designs, the dynamic behavior of the process was neglected, i.e., only steady-state optimizations were conducted. To take the transient behavior into account, both designs were tested with a predefined disturbance function. While the robust reactor is able to hold the imposed constraint on carbon conversion for a longer period, it ultimately also operates in an infeasible area.

In our case study, the disturbance is caused by fluctuating wind conditions. These wind conditions can be measured and predicted within a short time horizon. This provides an opportunity for a control scheme via dynamic optimization. To ensure feasible operation, an additional dynamic optimization was performed using fluctuation-specific split ratios and shell temperatures as degrees of freedom. While the nominal design converges to an infeasible point, the optimizer is able to control the operating conditions such that the robust reactor becomes feasible. Of course, the results depend on the imposed fluctuation function: its respective upper and lower bounds, as well as its dynamics. In a real application, additional buffer tanks would be used to make sure that the fluctuations exhibit desirable characteristics. Those tanks could also be included in the optimization problem to account for the buffer costs, which is, however, not the scope of this paper.

The approach presented in this paper can be interpreted as a feed-forward control scheme. In future works, we plan to enhance the dynamic optimization step used for calculating the control action with a machine learning-based approach [45]. For disturbances beyond the considered scenarios, an additional buffer tank for H₂ would be required. Simultaneous optimization of reactor and buffer designs under dynamic conditions is another interesting topic for future work. Furthermore, we intend to augment the overall control scheme with a feedback loop or, alternatively, conduct repetitive online optimizations to account for additional disturbances and plant model mismatch. The number of stages used in the present work was chosen exemplarily. Further work will be concerned with finding the optimal amount of stages for different scenarios. Additionally, the CSTRs will be replaced by industrially more relevant fixed bed reactor segments

In conclusion, this paper presents a novel approach to addressing the challenge of fluctuations in feed streams for methanol synthesis from renewable feedstock. By combining robust optimization with dynamic optimization, we have demonstrated that it is possible to generate a reactor design that is dynamically feasible for measurable fluctuating inlet conditions in a limited range.