Optimization of Methanol Synthesis under Forced Periodic Operation

: Traditionally, methanol is produced in large amounts from synthesis gas with heterogeneous Cu/ZnO/Al 2 O 3 catalysts under steady state conditions. In this paper, the potential of alternative forced periodic operation modes is studied using numerical optimization. The focus is a well-mixed isothermal reactor with two periodic inputs, namely, CO concentration in the feed and total feed ﬂow rate. Exploiting a detailed kinetic model which also describes the dynamics of the catalyst, a sequential NLP optimization approach is applied to compare optimal steady state solutions with optimal periodic regimes. Periodic solutions are calculated using dynamic optimization with a periodicity constraint. The NLP optimization is embedded in a multi-objective optimization frame-work to optimize the process with respect to two objective functions and generate the corresponding Pareto fronts. The ﬁrst objective is the methanol outlet ﬂow rate. The second objective is the methanol yield based on the total carbon in the feed. Additional constraints arising from the complex methanol reaction and the practical limitations are introduced step by step. The results show that signiﬁcant improvements for both objective functions are possible through periodic forcing of the two inputs considered here.


Introduction
Methanol is one of the most important raw materials of the chemical industry. It is used to produce paraffins, olefins, other organic chemicals, and fuels [1]. Commonly, it is produced in large amounts from synthesis gas using heterogeneous catalysts like Cu/ZnO/Al 2 O 3 under steady state conditions. More recently, there has also been a growing interest in dynamic methanol reactor operation in the context of energy storage and power to methanol processes [2]. In these types of processes, green hydrogen is produced from excess wind or solar energy via electrolysis. It reacts with CO/CO 2 from biogas or waste streams to methanol [3][4][5]. This will result in a more flexible use of electrical energy from renewable resources and simultaneously in a reduced emission of greenhouse gases to the environment. However, due to the fluctuating availability of the reactants, the reactor will face dynamic varying feeds. The focus of this paper is on improving reactor performance compared to conventional steady state operation by forced periodic operation, which is a specific type of dynamic reactor operation. The feed can be provided by renewable resources, conventional sources, or a mixture of both. Forced periodic operation can be beneficial if the time averaged reactor performance under forced periodic operation is higher than the corresponding steady state value, which is only possible for nonlinear systems (see, e.g., [6]). In general, the idea is not new and has been discussed since the 1960s (see, e.g., [7][8][9][10]). For a recent overview, we refer to [6,11]. For a rigorous evaluation of new forced periodic operation modes, a comparison with optimal steady state conditions is essential. First, work on forced periodic operation of methanol synthesis conducted on some Cu/ZnO and Cu/ZnO/Al 2 O 3 catalysts at temperatures between 225 and 270 • C and relatively moderate pressures between 1.97 and 2.93 MPa was done by [12,13]. The results indicate a potential of the methanol synthesis for improvement through periodic operation. However, a rigorous optimization of forced periodic operation and corresponding steady states was not performed due to the experimental focus of this work. Hence, a rigorous evaluation of forced periodic operation of methanol synthesis is lacking. This gap is closed in the present paper. Numerical NLP optimization is applied to compare optimal steady state solutions with optimal periodic regimes. Multi-objective optimization is applied to trace out the Pareto fronts to simultaneously maximize the methanol flow rate and the methanol yield with respect to the total carbon feed. Additional constraints arising in practice are integrated. In the first step, a well-mixed isothermal reactor with multiple periodic inputs is considered, namely, CO partial pressure in the feed and total feed flow rate. The underlying mathematical model is based on the lumped kinetic model presented in [14,15]. It accounts for dynamic changes of the catalyst and shows good agreement with steady state and dynamic experimental data from Vollbrecht [16].

Kinetic Model
Methanol synthesis from CO/CO 2 , H 2 over a commercial Cu/ZnO/Al 2 O 3 catalyst considered in this paper comprises the following three reactions The first and the second reaction represent methanol formation by hydrogenation of CO and CO 2 . The third reaction is the reverse water-gas shift reaction. For the computations in this paper, we use the simplified reaction kinetics presented in [14,15]. It assumes a Langmuir-Hinshelwood mechanism with three different active surface centers. The resulting expressions for the three reaction rates are: with the reformulated Arrhenius equation: The reference temperature used is T ref = 523.15 K and equilibrium constants are according to [16,17] : The surface coverages needed in Equations (4)-(6) are: These three surface centers have certain properties and are important for different reactions: i : for oxidized surface centers, also assumed as active center for CO hydrogenation; ii : * for reduced surface centers, also assumed as active center for CO 2 hydrogenation; iii : ⊗ as active surface centers for heterolytic decomposition of hydrogen.
It is further assumed that the oxidized surface centers and the reduced surface centers can be reversibly transformed into each other while keeping the total amount of active centers constant [14,15]. To describe this transformation under dynamic conditions, we introduce the following differential equation: A maximum amount of 90 % of reduced centers, i.e., φ ≤ 0.9, was assumed, corresponding to φ max = 0.9 in the above equation [15]. The equilibrium constants K 1 and K 2 (Equation (14)) can be expressed as a function of the free energies according to Parameters used in this paper are given in Table 1. They have been refitted to the experimental data of Vollbrecht [16] using φ ≤ 0.9 in the constraint set of the nonlinear least squares problem instead of shifting φ by 0.1 as was done in [14,15] and also expressed in Equation (24) by [18]. Hence, the parameter values in this paper are slightly different from those reported in [14].

Reactor Model
In this paper, a well-mixed, isothermal CSTR with constant pressure and ideal gas phase is considered, corresponding to the lab-scale micro-Berty reactor described in [16]. In the following section, the subscripts i and k describe the species (i.e., CH 3 OH, CO 2 , CO, H 2 , H 2 O, N 2 , total number N k = 6) and the subscript j describes the reaction (total number N r = 3). The superscript G denotes the gas phase and S the solid phase. Under the abovementioned assumptions, the model equations follow from the overall material balances of each component i according to with Therein, Θ i is the total coverage of component i at the different surface centers of the solid phase. It is assumed that the Θ i s are in equilibrium with the gas phase and therefore depend on the partial pressures of the gas phase according to Θ i = Θ i (p 1 , p 2 , ...p N k ), as described in detail in Appendix A. Further, due to the assumption of an ideal gas phase and a constant pressure, n G is constant. Therefore, the time derivatives of n G i and n S i are The derivatives ∂Θ i /∂p k are also explained in Appendix A. Substituting Equations (20) and (21) into Equation (17) yields with p k = py k The total material balance is obtained by summation of the component material balances over all species: Equation (23) is used for the calculation of the flow rateṅ at the outlet, which is, in general, different from the flow rate at the inlet, because the number of moles and also the volumes are not preserved by chemical reaction network Equations (1)-(3) considered here. For the numerical calculations, the flow rate at the outlet in the component material balances (22) was eliminated using the total material balance (23). Finally, conversion of the material balances to volumetric quantities can be done using the ideal gas law with equal pressures and temperatures in the in-and outlet p 0 = p, T 0 = T for the isothermal and isobaric mode of operation considered here.

Steady State Optimization
For a rigorous evaluation of improvements by forced periodic operation, the calculation of optimal steady state conditions is required as a reference. Steady state optimization is formulated as a nonlinear program (NLP) according to where J represents a suitable performance criterion (production rate or yield), y the state variables and x additional optimization variables. In the present paper, pressure, temperature and the steady state flow rate at the inlet are fixed and the inlet composition of the reaction mixture is optimized. Details and bounds on the variables are given in the Results section. Equality constraints are represented by the model equations. Additional equality and inequality constraints are also discussed in the Results section. The NLP problem was solved in MATLAB 2020b using fmincon with a multi-start algorithm with 250 starting points to avoid suboptimal local minima, which are a well-known problem for the present type of nonlinear and hence nonconvex optimization problem (see, e.g., Horst and Tuy [19].

Optimization of Forced Periodic Operation
Following the ideas in [20][21][22][23], simultaneous periodic forcing of two input variables u periodic,1 and u periodic,2 is considered, which is often more promising than periodic forcing of a single input variable. In particular, the focus of this paper is two harmonic forcing functions according to u periodic,2 (t) = u SS,2 (1 + A 2 cos(ωt + ∆φ)).
The forcing parameters to be optimized are the amplitudes A 1 , A 2 , the forcing frequency ω and the phase shift ∆φ between the two input variables. Additionally, the op-erating point and the corresponding inputs u SS have to be optimized to find the most promising operating point for the forced periodic operation. In general, this operating point can differ from the optimal steady state operation. The corresponding optimization problem is formulated as where the first condition implies periodicity of the solution. The length of the period τ is determined by the forcing frequency ω according to τ = 2π/ω. The calculation of periodic solutions requires a solution of the dynamic model equations, as indicated by the dynamic constraints. Further inequality constraints will be discussed in more detail in the Results section. The dynamic optimization problem was solved in MATLAB 2020b with fmincon and a multi-start optimization approach with a sequential approach (see, e.g., Cervantes and Biegler [24]).

Multi-Objective Optimization
In the present case, multiple, rather than single, objectives are of interest for evaluating the reactor performance. The two objective functions to be considered in this paper are The first objective function is the time averaged molar outlet flow rate of methanol. The second objective function is the time averaged yield of methanol, based on the carbon species entering the reactor. Average values were calculated using a single period with 1000 equidistant time points.
The system for methanol synthesis is complex and nonlinear. A straightforward approach with a weighted sum of both objectives may give poor results for steep or even nonconvex Pareto fronts [25,26]. Therefore, the -constraint method is applied. This is a scalarization technique which reduces the multi-objective optimization problem to a single objective optimization with the second objective function as an additional constraint [27][28][29]. The quantity is used as a parameter to change this constraint between a given upper and lower bound of the second objective. For the optimization of periodic regimes, the multi-objective problem formulation reads, e.g.,  = h(x, y) g(x, y) < 0 lb ≤ x ≤ ub An analogous formulation applies to the multi-objective steady state optimization. Both multi-objective optimization problems are solved again in MATLAB 2020b with fmincon, and a multi-start heuristic with 250 starts per point on the Pareto front.

Steady State Multi-Objective Optimization
For steady state operation pressure, temperature and volumetric flow rate at the inlet are fixed and the composition of the reaction mixture at the inlet is optimized. Parameter values correspond to a lab-scale micro-Berty reactor according to Vollbrecht [16]. They are summarized together with the fixed operational parameters in Table 2. For steady state optimization, the following additional constraints are used: The first constraint represents the summation condition of all species in the feed. The second constraint guarantees at least 1% of carbon in the feed to ensure methanol production. Furthermore, lower and upper bounds for the feed species are implemented. Hydrogen is assumed to always be available in excess. A multi-objective optimization regarding Equations (27) and (28) was done and the resulting Pareto front is shown in Figure 1. Two reference points OP1 and OP4 are marked on the Pareto front of steady state operation represented by the crosses in Figure 1. At the Pareto optimal steady state operating point OP1, the methanol flow rate is 336 mmol/min/kg cat with a yield of 61%. At the Pareto optimal steady state operating point OP4, the methanol flow rate is 413 mmol/min/kg cat with a yield of 52%. The corresponding optimal steady state feed concentrations are shown in the first three diagrams of Figure 2. It can be seen from Figure 2 that the optimal feed contains more CO than CO 2 . This is a well-known fact, which follows from the equilibrium limitations of the reaction network and in particular the water inhibition for methanol production (see, e.g., Vollbrecht [16] and references therein). The optimal CO and CO 2 concentrations in the feed are continuously decreasing from the left to the right along the Pareto front.

Multi-Objective Optimization of Forced Periodic Operation
Using the nonlinear frequency response method [30,31], it was found that simultaneous periodic modulation of the CO feed concentration and the volumetric flow rate at the inlet according to y CO,0 (t) = y CO,0,SS (1 + A CO cos(ωt)), are most promising for forced periodic operation.
To satisfy the summation condition for all y i,0 at any time t, the N 2 inert gas feed concentration is also varied periodically in an opposite way in order to compensate the periodic change of the CO feed concentration y N 2 ,0 (t) = y N 2 ,0,SS (1 − A N 2 cos(ωt)).
(34) ForV 0.SS , a fixed value of 1.14 × 10 −7 m 3 /s, and for y N 2 ,0,SS , a fixed value of 0.15 was used (see also Table 2). The other parameter values in Table 2 were also used for the optimization of forced periodic operation. With the CO modulation, the carbon content is varied. The flow rate modulation affects the residence time.
Additional constraints used for the optimization of forced periodic operation are y CO,0,SS + y CO 2 ,0,SS + y H 2 ,0,SS + y N 2 ,0,SS − 1 = 0, Therein, Equations (32) and (34), together with constraints (35) and (36), guarantee that the summation conditions at the inlet are satisfied at each time t. Equation (38) represents the periodicity constraint, with an appropriate summation condition for the state variables at the beginning of the period (39).
Lower and upper bounds for the optimization variables are summarized in Equation (40). The amplitudes are bounded between 0 and 1. The phase shift is bounded between −π and π. It was found (data not shown) that the improvement through periodic forcing increases with increasing frequency. However, in view of practical implementation, the period time τ is limited here between 18 s and 3600 s.
An additional constraint arises when taking a look at the reaction Equations (1) -(3). The reactions of methanol synthesis are mole number and volume reducing. Therefore, an additional path constraint for the outlet flow rate is required to ensure a positive outlet flow at constant pressure at any time point according tȯ Results for the multi-objective optimization of forced periodic optimization are also shown in Figure 1 and compared to optimal steady state operation. Percentages of improvement for the individual objective functions (27) and (28) relative to the optimal steady state operating points OP1 and OP4 are also given in this figure. The improvement in the methanol flow rate is between 15% at OP5 and 24% at OP3 and increases from left to right along the Pareto front. The improvement for the yield is between 17% at OP6 and 5% at OP3 and decreases from left to right along the Pareto front. The corresponding optimized parameter values for forced periodic operation are shown in Figure 2. The trends of the mean inlet concentrations y i,0,SS follow the trends of the corresponding steady state values in the first three diagrams. The amplitude of CO forcing is increasing from left to right, whereas the amplitude of N 2 is decreasing correspondingly as explained above. The amplitude of the flow rate forcing is almost constant. The predicted optimal phase shift between CO concentration and flow rate is low over the whole range. The period τ is close to the lower boundary of 18 s most of the time. The outlier between a yield of 0.55 and 0.6 most likely represents a close by local optimum. The inlet concentrations and forcing parameters for the presented operating points in Figure 1 are summarized in Table  3. Finally, the values of the variables y i ,V and φ and the relative amount of reduced surface centers, at the beginning of the period at t = 0, are shown in Figure 3 to illustrate the evolution of the periodic solutions along the Pareto front.  In summary, we find that the Pareto front of optimal forced periodic operation always lies above the Pareto front of optimal steady state operation, with somewhat substantial improvements for the conditions considered here.

Conclusions
In this contribution, we have shown theoretically that the performance of methanol synthesis from synthesis gas with a commercial Cu/ZnO/Al 2 O 3 catalyst can be improved significantly by periodic forcing of the CO feed concentration and the phase-shifted feed flow rate compared to optimal steady state operation. Improvements were measured in terms of methanol flow rate and methanol yield relative to the total carbon in the feed. The focus was on harmonic forcing functions, and a well-mixed isothermal CSTR. Due to the complexity of the methanol synthesis using a Cu/ZnO/Al 2 O 3 catalyst, the isothermal CSTR should be considered as a first reasonable step to provide some fundamental insight into the effect of forced periodic operation on methanol synthesis.
In industry, methanol synthesis is usually performed in a cooled fixed bed reactor. Perfect isothermicity cannot be achieved in such a reactor and additional temperature effects play a role. They may put additional constraints on practical reactor operation but also offer additional degrees of freedom for periodic forcing, which we want to investigate in our future work. However, since the modeling of a fixed bed reactor leads to a set of partial differential equations, more advanced numerical methods are required for optimization purposes.
The results in this paper show that excess hydrogen in the feed is beneficial for the methanol flow rate and the yield in terms of total carbon for steady state and forced periodic operation. In practice, the cost of hydrogen can be a major factor, particularly for renewable methanol. For reduced hydrogen feed, the production rate of methanol and the yield will be reduced for steady state as well as forced periodic operation compared to the over-stoichiometric case considered in this paper. However, it is expected that forced periodic operation will be still superior to the steady state operation due to the additional degrees of freedom. To quantify the economic benefit, further calculations are required using an economic objective function accounting for the specific hydrogen supply cost. In view of practical applications, this will be done for a fixed bed reactor in our future work. Further, we are aiming at an experimental validation of the theoretical findings presented in this paper. Experiments, however, are not trivial due to the relatively high pressure and the properties of the reactants considered here and are therefore clearly beyond the scope of the present paper and a subject for future research.