Process Intensification Through Forced Periodic Operation: Nonlinear Frequency Response Analysis of an Isothermal CSTR for Methanol Synthesis

1. Introduction

Methanol is widely recognised as one of the most important platform chemicals in the modern chemical industry and one of the most promising liquid energy carriers for sustainable energy systems. Besides its conventional applications in the production of formaldehyde, acetic acid, methyl tert-butyl ether, dimethyl ether, and numerous petrochemical intermediates, methanol has recently gained substantial importance within the concept of carbon-neutral energy technologies and renewable energy storage systems [1,2]. Owing to its favourable physicochemical properties, including a high hydrogen-to-carbon ratio, liquid state under ambient conditions, relatively high volumetric energy density, and compatibility with existing transportation and storage infrastructure, methanol represents a highly attractive candidate for future low-carbon energy systems.

Particularly important is the possibility of producing so-called green methanol through catalytic hydrogenation of carbon dioxide using hydrogen generated by renewable electricity-driven water electrolysis. Such processes enable the simultaneous utilisation of captured carbon dioxide and the storage of intermittent renewable energy in chemical form [2]. Consequently, methanol synthesis technologies are increasingly viewed not only as conventional catalytic processes but also as integral components of integrated renewable energy systems and sustainable chemical manufacturing.

Despite substantial progress in catalyst development and reactor design, conventional industrial methanol synthesis is still predominantly operated under steady-state conditions. Steady-state reactor operation is generally preferred because of operational simplicity, easier process control, and stable product quality. However, for strongly nonlinear catalytic systems, steady-state operation does not necessarily correspond to optimal reactor performance. Numerous theoretical and experimental investigations performed during the last several decades demonstrated that periodically operated reactors may achieve significantly improved time-averaged conversion, selectivity, and productivity compared to corresponding steady-state systems [3,4,5,6,7,8].

Forced periodic operation (FPO) represents a process intensification strategy based on intentional time-dependent modulation of one or more reactor inputs. Typical periodically modulated variables include inlet concentration, total flow rate, reactor temperature, pressure, and catalyst activity. The underlying principle of FPO is the exploitation of nonlinear interactions between reaction kinetics, adsorption phenomena, catalyst surface dynamics, heat transfer, and mass transport. Depending on the nature and degree of the nonlinearities involved, periodic forcing may lead either to process intensification or process deterioration relative to corresponding steady-state conditions [5,6].

Theoretical foundations of periodically operated reactors were established during the 1960s and 1970s through pioneering investigations of Douglas, Bailey, Horn, and other researchers [3,9,10,11,12]. These studies demonstrated that periodically perturbed nonlinear reactors may exhibit improved time-average performance compared to steady-state systems. Subsequent investigations expanded these concepts toward generalised optimisation criteria, nonlinear dynamic analysis, waveform optimisation, and frequency-domain approaches for evaluation of periodically operated systems [7,13,14,15].

One of the most important analytical approaches developed for the evaluation of periodically operated nonlinear systems is the nonlinear frequency response (NFR) methodology. The NFR method is based on Volterra series representation of nonlinear systems and higher-order frequency response functions (FRFs), enabling analytical prediction of time-averaged process improvements without extensive transient simulations [16,17,18,19]. The methodology enables rapid identification of favourable forcing strategies, forcing frequencies, modulation amplitudes, waveform characteristics, and phase differences between simultaneously modulated inputs.

The NFR methodology has been successfully applied to chemical and electrochemical systems, including adsorption systems, membrane reactors, electrochemical cells, heterogeneous catalytic reactors, fixed-bed reactors, and continuous stirred tank reactor (CSTR) operating under both isothermal and non-isothermal conditions [16,17,18,19,20,21]. Methanol synthesis reactors were identified as particularly suitable candidates for periodic operation because of pronounced nonlinearities associated with catalytic kinetics, reversible reactions, adsorption equilibria, and catalyst surface transformations [20].

Previous investigations specifically analysed FPO of methanol synthesis reactors under both single input modulation and simultaneous modulation of two inputs [16,17]. These studies demonstrated that periodic modulation of inlet concentration and flow rate may substantially improve methanol productivity compared to steady-state operation. However, the majority of previously analysed cases focused primarily on sinusoidal forcing or generalised waveform analysis, while square-wave forcing strategies received considerably less attention despite their practical advantages and industrial applicability [22].

Square-wave forcing possesses several characteristics which make it highly attractive for industrial implementation. First, square-wave signals can be relatively easily generated using standard industrial control systems and switching equipment. Second, abrupt transitions between forcing states provide stronger nonlinear excitation compared to smooth sinusoidal forcing. Third, square-wave modulation often allows longer residence time near extreme operating conditions, potentially producing larger time-average process improvements. Finally, square-wave forcing is compatible with practical process operation strategies involving periodic switching between operating states.

Although the FPO of methanol synthesis reactors has already been investigated using the NFR methodology for both single input modulation and simultaneous modulation of two inputs under co-sinusoidal forcing conditions [16,17]. NFR analysis of FPO of an isothermal CSTR for methanol synthesis using square-wave modulation of inlet process variables has not yet been reported in the available literature. This research gap is particularly important because square-wave forcing possesses several practical and theoretical advantages compared to harmonic co-sinusoidal forcing strategies, as stated above. Consequently, investigation of square-wave periodic operation may provide new insights into nonlinear process intensification of methanol synthesis reactors.

Therefore, the present study focuses on detailed NFR analysis of an isothermal CSTR carrying out methanol synthesis under square-wave FPO. Both single-input modulation and simultaneous modulation of two inputs are systematically investigated. Furthermore, the obtained results are directly compared with previously published analyses performed under co-sinusoidal forcing conditions [16,17], allowing evaluation of similarities and differences between harmonic and square-wave modulation strategies, as well as identification of potentially superior operating conditions for methanol productivity improvement.

2. Nonlinear Frequency Response Methodology

The NFR method is an analytical and approximate theoretical framework used for evaluating the potential improvement of FPO processes [23]. The approach is founded on the analysis of the frequency response of weakly nonlinear systems subjected to periodic forcing around a corresponding steady-state [24]. By applying the NFR methodology, it is possible to determine: (i) whether periodic modulation of one or more process inputs can result in process enhancement; (ii) which operating conditions should be fulfilled to achieve such improvement, including the choice of manipulated input(s), waveform(s), forcing amplitudes, forcing frequencies, and phase difference in the case of simultaneous modulation of two inputs; and (iii) the magnitude of the expected improvement for a selected reactor performance criterion, thereby indicating whether the enhancement is practically significant or not.

For a stable nonlinear system, the response to periodic input modulation around the corresponding steady-state operating point is defined as a quasi-stationary response in the form of a complex periodic function (Figure 1) [25]. The NFR method is based on the concept of higher-order FRFs, which originate from the generalised Fourier transform and the Volterra series representation of nonlinear systems [16,17,18,19].

The response of a weakly nonlinear system subjected to periodic forcing can be represented as a sum of DC component, first, second and higher harmonics. For weakly nonlinear systems, the dominant contributions to the time-averaged reactor response are commonly determined by first-order and second-order FRFs, while higher-order contributions are often significantly smaller and may be neglected without substantial loss of accuracy [22].

For a stable nonlinear system subjected to periodic forcing around the corresponding steady-state operating point, the system response can be represented using the Volterra series expansion (Equation A1 in Supplement A) [26]. Then the output response of the nonlinear system can be presented as the sum of the steady-state component and an infinite number of higher-order contributions. Each term of the Volterra series corresponds to a particular order of nonlinearity of the system. The frequency-domain representation of these terms leads to generalised higher-order FRFs. The generalised FRFs represent the fundamental quantities of the NFR methodology and characterise the nonlinear dynamic response of the system in the frequency domain.

2.1. Periodic Forcing with Single and Simultaneous Input(s) Modulation(s)

If the input (x(t)) of the nonlinear system is periodically modulated in a general form around a previously established steady-state (x_s) (Equation (1)). The frequency response of the weakly nonlinear system output is a complex periodic function which can be represented as the sum of the output steady-state value (y_s), the basic harmonic (y_I), a non-periodic term (y_DC, so-called DC component) and an infinite number of higher harmonics (Equation (2)) [16,17].

$$x\left(t\right)=\underset{generalformofaperiodicfunction}{\underbrace{{\displaystyle \sum _{k=1}^n}{A}_k{e}^{j{\omega }_{k}t}}}$$

(1)

Where, A_k is the forcing amplitude, and ω is the forcing frequency.

$$\begin{gathered} y\left(t\right)=\underset{steady-stateterm}{\underbrace{y_s}}+\underset{non-periodicDCterm}{\underbrace{y_{DC}}}+\underset{Iharmonic}{\underbrace{B_I\mathit{cos}\left(\omega t+{\phi }_I\right)}}+\underset{IIharmonic}{\underbrace{B_{II}\mathit{cos}\left(2\omega t+{\phi }_{II}\right)}} \\ +\underset{IIIharmonic}{\underbrace{B_{III}\mathit{cos}\left(3\omega t+{\phi }_{III}\right)}+\dots } \end{gathered}$$

(2)

Where, B_i and φ_i are the amplitude and the phase shift of the ith output harmonic, respectively.

On the other hand, the nonlinear system, based on the Volterra series [27], generalised Fourier transform, and the concept of higher order frequency response functions can be replaced with a set of linear FRFs of different orders [5]. Those FRFs are directly related to the different harmonics and the DC component of the frequency response of the nonlinear system [15].

To evaluate the possible improvement of the reactor system subjected to the FPO, only the non-periodic term, i.e., the DC component, which describes the time-average behaviour, needs to be evaluated [14].

If the input of the reactor system is modulated in a cosine shape:

$$x\left(t\right)=x_{\mathrm{s}}+A\cos \left(\omega t\right)=x_s+{\displaystyle \frac{A}{2}}{e}^{j\omega t}+{\displaystyle \frac{A}{2}}{e}^{-j\omega t}$$

(3)

Based on the NFR method, the DC component of output y_DC can be approximately evaluated only based on a second-order approximation and thus, the corresponding asymmetrical second-order (ASO) FRF, $G_x^{\left(2\right)}\left(\omega ,-\omega \right)$, which corresponds to the output y to the modulated input x, as shown [5]:

$$y_{x,DC}\approx 2{\left({\displaystyle \frac{A}{2}}\right)}^2{G}_x^{\left(2\right)}\left(\omega ,-\omega \right)$$

(4)

For weakly nonlinear systems subjected to periodic modulation of multiple manipulated variables, the dynamic model of the nonlinear system in the frequency domain can be represented by several sets of FRFs [14]. In the present study, the analysis is given for a nonlinear system with two periodically modulated inputs and one output variable. Then, three series of FRFs are required for the complete evaluation of the influence of FPO on system performances [14]. Two of these series describe the individual influence of each input on the system output, while the third series consists of cross-FRFs that account for the cross-effect of simultaneous modulation of two inputs [5]. Unlike the individual FRF series, which contains first, second, etc. FRFs, the cross-series contains second- and higher-order functions [5].

The DC component of output y from the system, where two inputs (x and z) are periodically modulated around the previously established steady states (x_s and z_s), corresponds to the sum of the DC components of the individual inputs and the contribution corresponding to the cross-influence of both inputs [14]:

$$y_{DC}=y_{x,DC}+y_{z,DC}+y_{xz,DC}$$

(5)

The DC component of output z (y_z,DC) is defined in an analogous way as for the modulation of input x (Equation (4)), while the cross-DC component, y_xz;DC, for simultaneous modulation of two inputs can be given by the following expression [14]:

$$y_{xz,DC}\approx 2\left({\displaystyle \frac{A_x}{2}}\right)\left({\displaystyle \frac{A_z}{2}}\right)G_{x,z}^{\left(2*\right)}\left(\omega ,\phi \right)$$

(6)

The cross term ($G_{x,z}^{\left(2*\right)}\left(\omega ,\phi \right)$), which indicates the cross effect, is introduced, and it depends on the phase difference between two modulated inputs ($\phi$) as well as the cross asymmetrical second-order frequency response function (cross-ASO FRF, $G_{x,z}^{\left(2\right)}\left(\omega ,-\omega \right)$). This cross-ASO FRF is defined as follows [14]:

$$G_{x,z}^{\left(2*\right)}\left(\omega ,\phi \right)=\cos \phi Re\left(G_{x,z}^{\left(2\right)}\left(\omega ,-\omega \right)\right)+\sin \phi Im\left(G_{x,z}^{\left(2\right)}\left(\omega ,-\omega \right)\right)$$

(7)

We have investigated the simultaneous modulation of two inputs with the same forcing frequency, considering that the highest contribution of the cross term to the overall DC component can be achieved in such conditions, as concluded by Nikolić et al. [16]. Therefore, the overall expression of the non-periodic DC term of the output for simultaneous modulations of two inputs can be written as follows:

$$y_{DC}\approx 2{\left({\displaystyle \frac{A_x}{2}}\right)}^2{G}_x^{\left(2\right)}\left(\omega ,-\omega \right)+2{\left({\displaystyle \frac{A_z}{2}}\right)}^2{G}_z^{\left(2\right)}\left(\omega ,-\omega \right)+2\left({\displaystyle \frac{A_x}{2}}\right)\left({\displaystyle \frac{A_z}{2}}\right)G_{x,z}^{\left(2*\right)}\left(\omega ,\phi \right)$$

(8)

The derivation procedure of FRFs using the NFR method is standardised and can be found in our previous publications [16,17].

2.2. Application of the NFR Methodology for Square-Wave Modulation

In the present study, periodic operation with square-wave forcing is considered. Square-wave modulation is particularly attractive from the practical point of view because it can be experimentally implemented by simple switching between predefined operating levels. Also, previous studies have shown that the use of square-wave modulation results in greater improvements to the reactor system than the use of co-sinusoidal modulation [22,28].

Knowing that a square-wave periodic function can be approximated and presented in Fourier series form as a sum of harmonic functions of different frequencies [28], the periodic modulation of inputs x(t) and z(t) in a square-wave shape can be given as follows:

$$x\left(t\right)\approx {\displaystyle \frac{4A_x}{\pi }}{\sum }_{n=\mathrm{1,3},\dots }^N{\displaystyle \frac{1}{n}}\sin \left(n\omega t\right)={\displaystyle \frac{4A_x}{\pi }}{\sum }_{n=\mathrm{1,3},\dots }^N{\displaystyle \frac{1}{n}}\cos (n\omega t-{\displaystyle \frac{\pi }{2}})$$

(9)

$$z\left(t\right)\approx {\displaystyle \frac{4A_z}{\pi }}{\sum }_{n=\mathrm{1,3},\dots }^N{\displaystyle \frac{1}{n}}\sin \left(n(\omega t+\phi )\right)={\displaystyle \frac{4A_z}{\pi }}{\sum }_{n=\mathrm{1,3},\dots }^N{\displaystyle \frac{1}{n}}\cos (n(\omega t+\phi )-{\displaystyle \frac{\pi }{2}})$$

(10)

The Fourier series for the square-wave function contains only the odd harmonics, while even harmonics are equal to zero [28]. For practical application, approximated Fourier series of square-wave inputs are used by taking into account only the first N harmonics given in cosine form [28]. In accordance with this and previously introduced Equation (1), the square-wave inputs can be written in the complex form in the following way [28]:

$$x\left(t\right)\approx {\sum }_{-N}^N{\displaystyle \frac{1}{2}}\underset{A_{x\left(n\right)}}{\underbrace{{\displaystyle \frac{4A_x}{n\pi }}{e}^{-j{\displaystyle \frac{\pi }{2}}}}}{e}^{jn\omega t}, \mathrm{n}=1, 3,\dots , N$$

(11)

$$z\left(t\right)\approx {\sum }_{-N}^N{\displaystyle \frac{1}{2}}\underset{A_{z\left(n\right)}}{\underbrace{{\displaystyle \frac{4A_z}{n\pi }}{e}^{-j{\displaystyle \frac{\pi }{2}}}}}{e}^{j\left(n\left(\omega t+\phi \right)\right)}, \mathrm{n}=1, 3,\dots , N$$

(12)

Where, A_x(n) and A_z(n) are the complex amplitudes of the inputs x and z, respectively.

The overall expression of the non-periodic DC term of the output for simultaneous modulation of two inputs in a square-wave form, approximated with N harmonics by analogy with Equation (8), is given as follows:

$$\begin{gathered} y_{DC}\approx {\sum }_{n=\mathrm{1,3},\dots }^{N}2{\left({\displaystyle \frac{A_{x\left(n\right)}}{2}}\right)}^2\underset{ASOFRF}{\underbrace{G_x^{\left(2\right)}\left(n\omega ,-n\omega \right)}}+{\sum }_{n=1.3,\dots }^{N}2{\left({\displaystyle \frac{A_{z\left(n\right)}}{2}}\right)}^2\underset{ASOFRF}{\underbrace{G_z^{\left(2\right)}\left(n\omega ,-n\omega \right)}}+ \\ {\sum }_{n=\mathrm{1,3},\dots }^{N}2\left({\displaystyle \frac{A_{x\left(n\right)}}{2}}\right)\left({\displaystyle \frac{A_{z\left(n\right)}}{2}}\right)\underset{cross-ASOFRF}{\underbrace{\left[\cos \left(n\phi \right)Re\left(G_{x,z}^{\left(2\right)}\left(n\omega ,-n\omega \right)\right)+\sin \left(n\phi \right)Im\left(G_{x,z}^{\left(2\right)}\left(n\omega ,-n\omega \right)\right)\right]}} \end{gathered}$$

(13)

For the evaluation of the DC component for square-wave simultaneous modulation of two inputs (Equation (13)), the ASO-FRFs for each single input modulation, as well as the cross-ASO FRF for the Nth harmonic of the output, should be derived and evaluated. The resulting sets of algebraic equations are convenient to present and solve in matrix form. These equations have already been derived in a previous study for cosine modulation (first and second harmonic) [16,17], and now only the third, fifth, to ninth harmonics have been added.

2.3. Mathematical Model of the Methanol Synthesis in the Isothermal CSTR

Methanol is still predominantly industrially synthesised from synthesis gas, a gaseous mixture primarily composed of CO, CO₂, and H₂, and represents one of the most important bulk chemical production processes. It is important to note that these raw materials can also be from renewable sources [2]. The synthesis proceeds through the hydrogenation of both carbon monoxide and carbon dioxide according to the following reactions:

$$CO+2H_2\leftrightarrows C{H}_{3}OH$$

(14)

$$C{O}_2+3H_2\leftrightarrows C{H}_{3}OH+H_{2}O$$

(15)

Both methanol synthesis reactions are exothermic and involve a reduction in the total number of moles, i.e., a decrease in system volume. Consequently, thermodynamic equilibrium favours methanol formation at elevated pressures and relatively low temperatures. Simultaneously with methanol synthesis, the reverse water-gas shift (RWGS) reaction also occurs:

$${CO}_2+H_2\leftrightarrows CO+H_{2}O$$

(16)

Industrial methanol production is still predominantly based on low-pressure catalytic synthesis from syngas [29]. In these processes, Cu/ZnO/Al₂O₃ catalysts are most commonly employed due to their high activity and selectivity toward methanol formation.

Application of the NFR methodology requires, as an initial step, the formulation of an appropriate mathematical model of the investigated system. A detailed derivation and description of the model were presented in our previous studies [16,17], while only the essential aspects are briefly summarised here. The developed model is based on a reactor model incorporating a detailed kinetic description of methanol synthesis previously developed and validated against both steady-state and dynamic experimental data obtained in the Micro-Berty reactor system [20]. The kinetic formulation follows the Langmuir–Hinshelwood reaction mechanism, assuming that the overall process proceeds through adsorption of reactants on the catalyst surface, surface reactions between adsorbed species, and subsequent desorption of the reaction products [20].

The catalytic surface is assumed to contain three different categories of active sites: oxidised surface centres, reduced surface centres, and active sites responsible for heterolytic hydrogen dissociation. The fraction of reduced active centres is represented by the variable ϕ, whose temporal evolution is governed by catalyst oxidation–reduction dynamics with the corresponding dynamic equation [17]:

$${\displaystyle \frac{d\varphi }{dt}}=k_1^{+}\left(y_{CO}\left({\varphi }_{max}-\varphi \right)-{\displaystyle \frac{1}{K_1}}{y}_{{CO}_2}\varphi \right)+k_2^{+}\left(y_{H_2}\left({\varphi }_{max}-\varphi \right)-{\displaystyle \frac{1}{K_2}}{y}_{H_{2}O}\varphi \right)$$

(17)

Where, K₁ and K₂ are the equilibrium constants, and k₁⁺, k₂⁺ are the dynamic rate constants.

The kinetic model further assumes ideal gas-phase behaviour, negligible catalyst deactivation, and absence of side reactions other than the methanol synthesis and RWGS reactions introduced previously. Under these assumptions, the resulting kinetic expressions for CO hydrogenation, CO₂ hydrogenation, and the RWGS reaction are summarised in Supplement B together with the corresponding adsorption relationships, equilibrium expressions, kinetic, and adsorption parameters.

The mathematical model of the methanol synthesis reactor was formulated using several assumptions: methanol synthesis is carried out in an isothermal and isobaric CSTR, the gas mixture behaves ideally within the investigated operating range, adsorption equilibrium between the fluid and catalyst phases is established, adsorption phenomena follow the Langmuir–Hinshelwood mechanism, catalyst deactivation effects are neglected, and only the reactions defined by Equations (14)–(16) are considered.

Under these assumptions, the dynamic behaviour of the reactor system is described by a set of eight ordinary differential equations. These equations include the equation describing the catalyst dynamics (Equation (17)), the material balances for all six species present in the system (CH₃OH, CO₂, CO, H₂, H₂O, and inert N₂) (Equation (18)), and the overall material balance used for the determination of the outlet volumetric flow rate (Equation (19)).

$$V_G{\displaystyle \frac{d{p}_i}{dt}}+m_{cat}{q}_{sat}RT{\displaystyle \sum _{l=1}^6}{\displaystyle \frac{\partial {\theta }_i}{\partial p_l}}{\displaystyle \frac{d{p}_l}{dt}}=\dot{V_0}{p}_{i,0}-\dot{V}{p}_i+m_{cat}RT{\displaystyle \sum _{j=1}^3}{\nu }_{ij}{r}_j,i=1,\dots ,6$$

(18)

$$m_{cat}{q}_{sat}RT{\displaystyle \sum _{i=1}^6}{\displaystyle \sum _{l=1}^6}{\displaystyle \frac{\partial {\theta }_i}{\partial p_l}}{\displaystyle \frac{d{p}_l}{dt}}=\dot{V_0}{p}_{tot}-\dot{V}{p}_{tot}+m_{cat}RT{\displaystyle \sum _{i=1}^6}{\displaystyle \sum _{j=1}^3}{\nu }_{ij}{r}_j$$

(19)

Since the NFR analysis is more conveniently performed using a dimensionless formulation of the model, the process variables are transformed into dimensionless form by defining them as relative deviations from the corresponding steady-state values (Supplement C). These dimensionless variables are subsequently introduced into the governing model equations.

In the next step, all nonlinear terms appearing in the mathematical model are expanded into Taylor series around the selected steady-state operating point [5]. After introducing the dimensionless variables and incorporating the Taylor-expanded nonlinear terms into the governing equations, the complete dimensionless reactor model is obtained (Supplement C). The resulting dimensionless formulation serves as the basis for the derivation of higher-order frequency response functions and subsequent NFR analysis of periodically operated methanol synthesis reactors.

2.4. Process Improvement Evaluation

The outlet molar flow rate of methanol is chosen to be an indicator of improvement for square-wave single and simultaneous modulation of inlet fraction of CO (x=CO) and flow rate (z=F). The dimensionless DC component of the outlet molar flow rate should be evaluated to calculate the improvement. Based on (Equation (13)), for simultaneous modulation, it can be given as follows:

$$\begin{gathered} N_{DC}\approx {\sum }_{n=\mathrm{1,3},\dots }^{N}2{\left({\displaystyle \frac{A_{CO\left(n\right)}}{2}}\right)}^2{H}_{CO}^{\left(2\right)}\left(n\omega ,-n\omega \right)+{\sum }_{n=1.3,\dots }^{N}2{\left({\displaystyle \frac{A_{F\left(n\right)}}{2}}\right)}^2{H}_F^{\left(2\right)}\left(n\omega ,-n\omega \right) \\ +{\sum }_{n=\mathrm{1,3},\dots }^{N}2\left({\displaystyle \frac{A_{CO\left(n\right)}}{2}}\right)\left({\displaystyle \frac{A_{F\left(n\right)}}{2}}\right)(\mathit{cos}n\phi Re\left(H_{CO,F}^{\left(2\right)}\left(n\omega ,-n\omega \right)\right) \\ +\mathit{sin}n\phi Im\left(H_{CO,F}^{\left(2\right)}\left(n\omega ,-n\omega \right)\right)) \end{gathered}$$

(20)

For the single input modulation case of molar fraction of CO or total inlet volumetric flow rate, dimensionless DC components are reduced to:

$$N_{DC,CO}\approx {\sum }_{n=\mathrm{1,3},\dots }^{N}2{\left({\displaystyle \frac{A_{CO\left(n\right)}}{2}}\right)}^2{H}_{CO}^{\left(2\right)}\left(n\omega ,-n\omega \right)$$

(21)

$$N_{DC,F}\approx {\sum }_{n=1.3,\dots }^{N}2{\left({\displaystyle \frac{A_{F\left(n\right)}}{2}}\right)}^2{H}_F^{\left(2\right)}\left(n\omega ,-n\omega \right)$$

(22)

The DC components for all analysed cases are approximated with N first harmonics. H-ASO FRFs, which are directly related to the previously derived G-FRFs (given in [17]), correlate the outlet molar flow rate of methanol with the modulation of the inlet fraction of CO ($H_{CO}^{\left(2\right)}\left(n\omega ,-n\omega \right)$), flow rate ($H_F^{\left(2\right)}\left(n\omega ,-n\omega \right)$) and the cross-effect for simultaneous modulation of those inputs $H_{CO,F}^{\left(2*\right)}\left(n\omega ,-n\omega \right)$.

The mean, i.e., time-average, outlet molar flow rate of methanol can then be evaluated from derived DC components (Equation (20)-22) and the steady-state molar flow rate of methanol:

$${\dot{n}}_{mean}\left(C{H}_{3}OH\right)\approx {\dot{n}}_s\left(C{H}_{3}OH\right)\left(1+N_{DC}\right)$$

(23)

In the case of square-wave simultaneous modulation of the mole fraction of CO and inlet volumetric flow rate, the value of the inlet molar flow rate of CO is, in principle, different from the steady-state value, and it can be calculated using the following expression:

$${\dot{n}}_{CO,mean}\left(i\right)={\dot{n}}_{CO,s}\left(1+{\displaystyle \sum _{n=\mathrm{1,3},\dots }^N}\left({\displaystyle \frac{A_{CO\left(n\right)}}{2}}\right)\left({\displaystyle \frac{A_{F\left(n\right)}}{2}}\right)\cos \left(n\phi \right)\right)$$

(24)

Based on the time-average value of the molar flow rate of methanol (Equation (23)), the following additional performance indicators can be defined [16]:

• Yield of methanol based on total carbon

$$Y_{totC}\left({\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}\right)={\displaystyle \frac{{\dot{n}}_{mean}\left(C{H}_{3}OH\right)}{{\dot{n}}_{0,ss}\left({CO}_2\right)+{\dot{n}}_{0,mean}\left(CO\right)}}$$

(25)

3.

Yield of methanol based on hydrogen

$$Y_{H_2}\left({\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}\right)=2{\displaystyle \frac{{\dot{n}}_{mean}\left(C{H}_{3}OH\right)}{{\dot{n}}_{0,ss}\left(H_2\right)}}$$

(26)

4.

Normalised methanol production per unit of catalyst

$${\dot{n}}_{norm}\left({\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}\right)={\displaystyle \frac{{\dot{n}}_{mean}\left(C{H}_{3}OH\right)}{m_{cat}}}$$

(27)

For simultaneous modulation of two inputs, the phase difference has a decisive role, and it is possible to obtain improvement even in cases when single-input modulations would lead to deterioration of reactor performances [14]. The forcing amplitudes and phase difference should be optimised in a way that the maximal value of outlet molar flow rate is achieved throughout maximisation of the corresponding DC component (Equation (20)). Considering that for a square-wave, the DC component of the first harmonic is dominant [28], in optimisation of forcing parameters, only the first harmonic and corresponding ASO FRFs are included.

The optimal phase difference was analytically determined as a function of forcing frequency, as follows [16]:

$${\phi }_{opt,CO,F}=arctang\left({\displaystyle \frac{Im\left(H_{CO,F}^{\left(2*\right)}\left(\omega ,-\omega \right)\right)}{Re\left(H_{CO,F}^{\left(2*\right)}\left(\omega ,-\omega \right)\right)}}\right)$$

(28)

Then, the optimal forcing amplitudes were determined using the fminimax solver in MATLAB based on derived cross-ASO H-FRFs, a defined DC component with Equation (20), and values of optimal phased difference (Equation (28)).

3. Results and Discussion

The present study focuses on square-wave FPO of the isothermal CSTR for methanol synthesis previously investigated under co-sinusoidal forcing conditions by Nikolić et al. [16,17]. In the preceding study, simultaneous modulation of the inlet CO concentration and inlet volumetric flow rate was identified as the most favourable forcing scenario among all investigated combinations of manipulated variables. Therefore, the current analysis is restricted to this previously established optimal pair of inputs, while the influence of the forcing waveform on the attainable reactor performance is investigated in detail.

The same reactor model, kinetic scheme, catalyst dynamics, and steady-state operating point as those reported by Nikolić et al. in Part I and Part II are used throughout this study [16,17]. Consequently, all differences observed between the previously reported results and those presented here originate exclusively from the replacement of harmonic forcing by square-wave forcing. Such an approach enables a direct assessment of the influence of the forcing waveform on the nonlinear dynamic response of the methanol synthesis reactor.

Square-wave forcing introduces a significantly richer frequency spectrum compared to co-sinusoidal forcing. While harmonic forcing contains only the fundamental frequency, square-wave signals consist of an infinite series of odd harmonics whose amplitudes decrease proportionally to the inverse harmonic order. As demonstrated previously for generic nonlinear systems and adiabatic CSTRs, these additional harmonics may contribute substantially to the non-periodic component of the reactor response and consequently modify the attainable process improvement [22,28]. The present section, therefore, combines NFR analysis with productivity and yield calculations in order to evaluate the potential of square-wave forcing for intensification of methanol synthesis.

3.1. Reference Productivity Results Under Cosine Forcing

Before discussing the reactor performance obtained under square-wave forcing, it is useful to briefly recall the productivity trends previously reported for harmonic forcing of the same methanol synthesis reactor. The corresponding productivity trends are summarised in Figure D1, reproduced from the earlier studies and included in Supplement D for completeness [16,17]. The performance indicators of the chosen optimal steady-state are as follows [17]:

• Methanol productivity of 336.91 mmol min^-1 ${kg}_{cat}^{-1}$;
• Yield of methanol based on total carbon of 61.05%;
• Yield of methanol based on total hydrogen of 39.09%.

As shown in Figure D1, periodic modulation of the inlet CO alone produces only a minor deviation from the steady-state productivity after which methanol productivity asymptotically approaches the steady-state value. A more pronounced effect was reported for periodic modulation of the inlet volumetric flow rate. In this case, the productivity decreases significantly below the steady-state value over a broad range of forcing frequencies. The strongest deterioration occurs in the intermediate-frequency region, where the negative contribution of the corresponding ASO FRFs reaches its maximum magnitude [17].

In contrast to the individual forcing scenarios, simultaneous harmonic modulation of the inlet CO and inlet volumetric flow rate produces a substantial increase in methanol productivity. As demonstrated by Nikolić et al. [16], the positive nonlinear interaction between the two manipulated variables overcomes the unfavourable individual contributions and generates a significant positive DC component of the reactor response. As a result, methanol productivity increases considerably above the steady-state value and reaches a maximum enhancement of approximately 33% under the optimal forcing conditions [16].

Figure D1 (Supplement D), therefore, serves two important purposes in the context of the present study. First, it establishes the previously reported productivity benchmark obtained under harmonic forcing. Second, it provides the rationale for restricting the current analysis exclusively to simultaneous modulation of the inlet CO concentration and inlet volumetric flow rate. Since this forcing pair was already identified as the optimal simultaneous modulation scenario [16], the present work focuses on determining whether the productivity enhancement previously achieved under harmonic forcing can be further increased by replacing the forcing waveform with a square-wave signal.

3.1. Nonlinear Frequency Response Analysis of Square-Wave CO and Flow Rate Modulation

The ASO H-FRFs, which correlate the outlet molar flow rate to inlet modulation of CO fraction, corresponding to square-wave modulation for up to the 9th harmonics, are presented in Figure 2. It can be observed that all investigated odd harmonics produce negative contributions over the entire frequency range. The first harmonic exhibits the largest absolute magnitude and therefore dominates the overall response, while the contributions of higher harmonics decrease progressively with increasing harmonic order. Such behaviour is fully consistent with previous theoretical analyses of square-wave forcing, where it was shown that the contribution of individual harmonics generally decreases with harmonic number [28].

The ASO H-FRFs, which correlate the outlet molar flow rate of methanol to modulated inlet volumetric flow rate, for square-wave modulation for up to the 9th harmonics, are presented in Figure 3. The overall behaviour is qualitatively similar to that observed for CO modulation. Again, all harmonics exhibit negative values over the entire frequency domain, indicating that single input modulation of the inlet flow rate would also lead to deterioration of the time-average reactor performance. The first harmonic remains dominant and reaches significantly larger absolute values than higher harmonics, meaning that its influence is highest.

Since ASO H-FRFs for single CO and flow rate modulations are negative throughout the investigated frequency range, they would not generate an improvement relative to steady-state operation. This observation is consistent with the conclusions previously obtained for harmonic forcing [16,17]. In those studies, both CO and flow rate single modulations would lead to a worsening of reactor performance.

The real and imaginary parts of the cross-ASO H-FRFs for simultaneous modulation of those two inputs in a square-wave shape, for up to 9th harmonics, are presented in Figure 4. According to the NFR theory for systems with two periodically modulated inputs, the cross contribution depends simultaneously on the values of cross-ASO H-FRFs as well as the phase difference between the input signals, as explained above (Equations (7) and (28)) [5]. The real and imaginary parts of the cross-ASO H-FRF for the first, dominant harmonic, provide the information necessary for the determination of the optimal phase difference (Equation (28)) [5,22].

The optimal phase difference and the corresponding optimal forcing amplitudes obtained from the NFR analysis are shown in Figure 5. The optimal phase shift exhibits an s-profile curve decreasing from 0 rad to -2π rad as the frequency increases. At higher frequencies, the phase difference asymptotically approaches an approximately constant value (2π rad).

At lower frequencies (ω<0.62), the optimal amplitude of the CO and Flow rate modulations is equal to 0, meaning that there is no sense to subject reactor system to FPO using those frequencies, i.e., FPOs would not lead to improvement. As the forcing frequency increases, the amplitudes increase sharply, reaching maximum values (A_CO of 0.81 and A_f of 1.0). This trend is fully consistent with the methanol productivity enhancement observed later in Figure 6.

3.2. Process Intensification Under Square-Wave Forcing

The predicted methanol productivity under square-wave forcing is analysed in this section. The methanol productivity calculated for single inputs and simultaneous modulations of CO and flow rate is presented in Figure 6. The figure summarises the cumulative effect of all harmonic contributions and directly reflects the overall DC component calculated by the NFR methodology (Equations (20)–(22)).

The productivity obtained under periodic single modulation of the inlet CO remains close to the steady-state value over the entire investigated frequency range. A slight reduction is observed at intermediate frequencies (ω≈0.62), corresponding to the frequency region where the first harmonic ASO H-FRF reaches its minimum value (Figure 2). The overall magnitude of the productivity change remains relatively small, confirming that modulation of the inlet CO concentration alone does not provide a practically significant route toward reactor intensification nor deterioration. This result is fully consistent with the conclusions previously obtained for co-sinusoidal forcing of the same reactor system, where concentration modulation alone also failed to generate substantial productivity enhancement [17].

Periodic modulation of the inlet volumetric flow rate produces a considerably larger effect. In agreement with the negative values of the corresponding ASO H-FRFs shown in Figure 3, a pronounced productivity decrease occurs over a broad frequency range (maximum decrease of 13,8% relative to the steady-state occurs in a very narrow frequency range). The minimum productivity is observed in the intermediate-frequency region (ω∈[0.62, 8.5]) where the negative contribution of the first harmonic is largest. Although the productivity asymptotically returns toward the steady-state value at higher forcing frequencies, no frequency region exists where flow rate modulation alone becomes beneficial. Consequently, the results clearly demonstrate that periodic modulation of the inlet flow rate cannot be considered an attractive forcing strategy when applied independently.

A fundamentally different behaviour is obtained when both inputs are modulated simultaneously. In this case, at the frequencies lower than 0.62, the system behaves quasi-stationary, after which the reactor productivity increases significantly above the steady-state value over a wide range of forcing frequencies. The improvement becomes particularly pronounced in the intermediate- and high-frequency regions, where the positive cross contribution dominates the overall DC component.

The maximum productivity predicted under simultaneous square-wave forcing reaches approximately 553.6 mmol min⁻¹ ${kg}_{cat}^{-1}$ for the frequencies higher than 10, compared to approximately 336.9 mmol min⁻¹ ${kg}_{cat}^{-1}$ under steady-state operation. This corresponds to a productivity increase of approximately 64.3%. Such an enhancement is remarkably large considering that no changes in catalytic reaction, reactor design, operating pressure, or operating temperature are introduced. The productivity increase is achieved exclusively through the dynamic operation of the reactor around the same optimal steady-state operating point. It is important to note that this enhancement is achieved around exactly the same steady-state operating point employed previously for harmonic forcing studies [17]. Therefore, the observed improvement cannot be attributed to a more favourable selection of operating conditions but must originate exclusively from the different nonlinear reactor response generated by the square-wave forcing signal.

Another important feature visible in Figure 6 is the relatively broad frequency range over which significant productivity enhancement is obtained. The productivity maximum is not restricted to a narrow frequency interval but extends over a wider region of forcing frequencies. From a practical perspective, such behaviour is advantageous because it reduces the sensitivity of reactor performance to small deviations from the optimal operating conditions.

The results presented in Figure 6 demonstrate that square-wave forcing preserves the favourable interaction mechanism previously identified under harmonic cosine operation while simultaneously amplifying its effect through the contribution of higher harmonics. As a consequence, substantially larger productivity improvements become attainable under otherwise identical operating conditions. To evaluate the practical significance of this improvement, a direct comparison between square-wave and co-sinusoidal forcing is presented in the following section.

3.3. Comparison with Cosine Forcing

The results presented in the previous sections demonstrate that square-wave forcing can substantially enhance methanol productivity through the combined action of favourable nonlinear interactions and the cumulative contribution of higher harmonics. However, the practical significance of these findings can only be fully assessed through direct comparison with the previously investigated co-sinusoidal forcing strategy. Since the same reactor model, kinetic scheme, catalyst dynamics, operating conditions, manipulated variables, and steady-state operating point were employed in both studies [16,17], the comparison presented in this section allows the influence of the forcing waveform itself to be isolated and evaluated.

3.3.1. Methanol Productivity

The comparison between the maximum methanol productivities obtained under steady-state operation, co-sinusoidal forcing, and square-wave forcing is presented in Figure 7. The figure clearly demonstrates that the forcing waveform exerts a substantial influence on the attainable process improvement.

As previously reported by Nikolić et al.l. [16], simultaneous co-sinusoidal modulation of the inlet CO and inlet volumetric flow rate produces a significant increase in methanol productivity compared to steady-state operation. The productivity increases from approximately 336.9 mmol min⁻¹ ${kg}_{cat}^{-1}$ under optimal steady-state conditions to approximately 450 mmol min⁻¹ ${kg}_{cat}^{-1}$ under optimal harmonic forcing. This corresponds to a productivity improvement of approximately 33%, confirming the strong potential of periodic operation for methanol synthesis intensification. When the harmonic forcing signal is replaced by a square-wave signal, the attainable productivity increases even further. The maximum productivity reaches an improvement of approximately 64.3% relative to steady-state operation. Therefore, the productivity enhancement obtained under square-wave forcing is almost twice as large as that previously achieved under co-sinusoidal forcing.

The comparison presented in Figure 7 suggests that the additional productivity enhancement obtained under square-wave forcing originates primarily from amplification of the positive nonlinear interaction terms (Figure 4). Since the same manipulated variables, operating conditions, and optimisation procedure were employed in both studies, the observed difference must be associated with the higher harmonic content of the forcing signal. Consequently, the additional productivity enhancement observed under square-wave operation must originate exclusively from the different spectral composition of the forcing waveform. The present results, therefore, provide strong support for earlier theoretical predictions indicating that generalised periodic waveforms may generate larger non-periodic shifts than purely harmonic forcing functions [22,28].

Under cosine forcing, the reactor response is governed primarily by the contribution of the first and second harmonics. In contrast, square-wave forcing introduces a sequence of odd harmonics, each contributing to the overall DC component through the corresponding H-ASO- and cross-ASO H-FRFs. Although the magnitude of these individual contributions decreases with harmonic order, their cumulative effect eventually becomes significant. The higher harmonics, therefore, provide additional positive contributions to the nonlinear interaction term and increase the overall productivity enhancement. The observed behaviour is also consistent with the conclusions of previous theoretical studies with simple ideal reaction conditions [22,28]. The present work demonstrates that the same conclusion remains valid for a considerably more complex catalytic reactor system describing industrially relevant methanol synthesis.

From a process intensification perspective, the increase from approximately 33% to approximately 64% is particularly significant. More importantly, the square-wave forcing is much easier to implement in practice.

3.3.2. Methanol Yield Based on Total Carbon

Although productivity is generally the primary performance indicator for industrial methanol synthesis, evaluation of the methanol yield provides additional insight into the utilisation efficiency of reactants. Therefore, the influence of the forcing waveform on the carbon-based methanol yield was also investigated.

The results are presented in Figure 8. Under steady-state conditions, the reactor exhibits the highest carbon-based methanol yield among the three operating modes considered. Both periodic operating strategies lead to a slight reduction of this performance indicator relative to steady-state operation. The decrease in carbon-based yield is not unexpected. Increasing reactor productivity generally requires operation under conditions that favour larger overall reaction rates, which do not necessarily coincide with conditions maximising conversion efficiency of the carbon-containing reactants. Consequently, a trade-off between productivity and carbon utilisation efficiency appears.

Nevertheless, several important observations can be made. First, the reduction in carbon-based yield remains relatively moderate compared to the corresponding increase in productivity. While productivity increases by more than 64% under square-wave forcing, the associated decrease in carbon-based yield is considerably smaller (19%). Second, the differences between harmonic and square-wave forcing remain relatively limited. Both forcing strategies exhibit similar trends, indicating that the fundamental reactor behaviour considering this indicator remains unchanged despite the different waveforms. From a practical perspective, this observation is encouraging because it indicates that the productivity improvement is achieved without severely compromising overall reactor selectivity toward methanol formation.

3.3.3. Methanol Yield Based on Hydrogen Utilisation

The influence of the forcing waveform on hydrogen utilisation is presented in Figure 9. In contrast to the carbon-based yield, the hydrogen-based methanol yield exhibits a markedly different behaviour. Under steady-state operation, the hydrogen-based yield attains the lowest value among the three investigated operating modes. The introduction of periodic forcing significantly improves hydrogen utilisation, and the magnitude of this improvement depends strongly on the forcing waveform.

As previously reported for co-sinusoidal forcing [16], simultaneous modulation of the inlet CO and flow rate increases the hydrogen-based yield substantially above the steady-state value. However, square-wave forcing produces an even larger improvement. The maximum hydrogen-based yield obtained under square-wave operation reaches approximately 64.2% (ω>8), whereas cosine forcing achieves approximately 52.2%.

The improvement of hydrogen utilisation is particularly important in the context of green methanol synthesis. In emerging Power-to-Methanol processes, hydrogen generated by water electrolysis typically represents the most energy-intensive and economically valuable reactant [2]. Consequently, efficient utilisation of hydrogen becomes one of the most important performance criteria of the overall process assessment.

Although a detailed economic analysis is beyond the scope of the present work, the observed increase in hydrogen-based yield clearly indicates that square-wave forcing improves not only methanol productivity but also the effectiveness of hydrogen conversion into the desired product. Therefore, the benefits of square-wave operation extend beyond simple productivity enhancement and include improved utilisation of one of the most valuable reactants involved in sustainable methanol production. The simultaneous increase of productivity and hydrogen-based yield is particularly noteworthy because such improvements often exhibit conflicting trends in catalytic reactor systems. In the present case, however, both indicators improve simultaneously under optimal square-wave forcing conditions.

Taken together, the results presented in Figure 7, Figure 8 and Figure 9 demonstrate that replacing co-sinusoidal forcing with square-wave forcing significantly improves the overall performance of the methanol synthesis reactor (CSTR). However, the richer harmonic structure of the square-wave signal produces larger positive nonlinear interactions and consequently larger non-periodic shifts of the reactor response. As a result, substantially greater productivity enhancements and improved hydrogen utilisation become attainable under otherwise identical operating conditions.

4. Conclusions

In the present study, the NFR methodology was applied to investigate square-wave FPO of an isothermal CSTR for methanol synthesis. The analysis was focused exclusively on simultaneous modulation of the inlet CO concentration and inlet volumetric flow rate, corresponding to the most favourable pair of manipulated variables previously identified for harmonic forcing [16]. The objective was to determine whether replacement of the forcing waveform could further intensify reactor performance while maintaining the same reactor model, operating conditions, catalyst characteristics, and steady-state operating point.

The frequency-domain analysis showed that the overall reactor response under square-wave forcing is determined not only by the fundamental harmonic but also by the cumulative contribution of higher odd harmonics. Although the magnitude of individual harmonic contributions decreases with increasing harmonic order, their combined effect produces an additional positive shift of the time-average reactor response. Consequently, the attainable process improvement under square-wave forcing exceeds that obtained under conventional harmonic forcing.

The predicted methanol productivity under optimal square-wave forcing reached approximately 553.6 mmol min⁻¹ ${kg}_{cat}^{-1}$, compared with approximately 336.9 mmol min⁻¹ ${kg}_{cat}^{-1}$ under optimal steady-state operation. This corresponds to a productivity enhancement of more than 64%. A direct comparison with previously published cosine forcing showed that replacing harmonic forcing with square-wave forcing significantly increases productivity by as much as two times [16]. In addition to the productivity enhancement, square-wave forcing also improved hydrogen utilisation efficiency. The hydrogen-based methanol yield increased substantially compared with both steady-state operation and co-sinusoidal forcing (64.3% increase compared to steady-state and 23.1% increase compared to cosine modulation [16]). Considering the growing importance of hydrogen utilisation in sustainable and Power-to-Methanol technologies, this observation further emphasises the potential practical significance of the proposed operating strategy.

Overall, the results demonstrate that square-wave FPO represents a highly effective intensification strategy for methanol synthesis in an isothermal CSTR. From a broader perspective, the present results demonstrate that optimisation of the forcing waveform can be as important as optimisation of the manipulated variables themselves. While previous studies identified the most favourable forcing variables for methanol synthesis in an isothermal CSTR, the present study shows that substantial additional improvements can be achieved through appropriate waveform design. This finding further extends the applicability of NFR methodology as a systematic tool for process intensification and dynamic reactor optimisation and highlights the potential of waveform engineering as a complementary strategy for improving the performance of catalytic reactor systems. Especially considering that the application of square-wave modulation is in reactor engineering, practically easier and more affordable than, for example, cosine modulation.

Future work should focus on rigorous numeric time-domain validation of the predicted optimal operating conditions, assessment of the influence of forcing waveform imperfections, and extension of the proposed methodology to non-isothermal reactor configurations and industrially relevant Power-to-Methanol process schemes.