Abstract

Acetyl chloride hydrolysis is a highly sensitive exothermic reaction that has presented several industrial safety issues. In the present study, a multiparameter mathematical model, previously developed and applied to simulate the oscillatory thermal behavior of an experimental continuous stirred tank reactor, was used to determine the static/dynamic bifurcation behavior of this reactive system. The values predicted by the model showed good agreement with the experimental data reported in the literature. Full topological classification of its fixed points and iterative maps was obtained: unique solutions (stable and unstable), multiple solutions, cyclic envelope, and bifurcation objects of codimension 1 (e.g., fold and Hopf’s points) and codimension 2 (e.g., cusp and generalized Hopf and Bogdanov-Takens points) have been uncovered. The emphasis of the analysis is to determine safe operating conditions through understanding these topological features and manipulating the reactor design and operating parameters.

1. Introduction

Acetyl chloride, CH3COCl, is a highly volatile and corrosive organic acid chloride, derived from acetic acid. Due to these properties, its synthesis is usually carried out near to the place of its application. It is used as an acetylating agent for the production of esters and amides, in the pharmaceutical, dye and liquid crystal industries, and in the synthesis of acetophenone by the Friedel-Crafts acylation of benzene [1, 2]. In the presence of water, the highly exothermic and spontaneous reaction of acetyl chloride hydrolysis can occur ( cal/mol acetyl chloride). Consequently, a dangerous phase transition (from liquid to vapor) can take place. A clear example of the potential hazard of this reaction occurred in 1997 at Caird Environmental Ltd., located in Minworth, Birmingham (UK). During the cleaning operation of an empty 45-gallon drum, residual acetyl chloride was exposed to water, rapidly releasing a gas that caused an explosion. As a result, the drum flew like a rocket towards the surrounding buildings [3]. The total number of such tanks in the company was 80. So what would be the magnitude of the explosion if more acetyl chloride had been involved? Unfortunately, there is no precise information on the amount of acetyl chloride handled by the company or the record of related accidents and/or prevention methods.

Risk assessment has already been shown to be of paramount importance in managing reactive systems. In fact, safety has become an essential factor in the design and operation of reactive process. Furthermore, learning to predict and prevent the hazards of chemical process is expected to be an essential part of the education of chemical engineers [46]. In this sense, different methods of analysis have been developed. The parametric sensitivity test is a practical tool for thermal hazard analysis. However, only some critical points can be identified with this technique (e.g., maximum temperature points or hot spots) [79]. However, their criteria are very conservative and are not sufficient to fully describe the phenomenon of thermal instability of a reactive system [5]. Another method is dynamic (bifurcation) analysis, which allows rigorous tracking of different variables over a wide range of operating conditions. This implies the mathematical application of the singularity theory and/or bifurcation theory. The first defines the possible topological shapes and features that a dynamic model can take (i.e., using bifurcation diagrams). Furthermore, it provides a useful framework to classify bifurcation phenomena (evidencing the different types of multiplicity) and their mathematical singularities. The research groups of Uppal et al. [10, 11] and Balakotaiah et al. [7, 1216] applied this theory for the study of the steady-state multiplicity and dynamic behavior of hypothetical reactions in a CSTR. They used parameter-parameter diagrams and phase diagrams to qualitatively explain and characterize the singular behavior of each of the formed regions. Later, Golubitsky and Schaeffer [17] developed a special variant of this theory; they called it the singularity theory with a distinguished parameter. It defines an “organizing center” as a reference point at which the highest number of singularities of the system occurs. From this point, unfolding and parametric projections of space can be made, which makes it possible to predict all kinds of bifurcation diagrams. Gray and Roberts [18] presented an excellent summary of this theory and show the most representative singular points, their respective nondegeneracy conditions, unfolding, and bifurcation diagrams. Subsequently, they applied it to study the nonlinear behavior of carbon monoxide oxidation [19]. Chang et al. [20] published a review of this theory and its application to the study of the stability of hypothetical reactive systems carried out in a recirculating CSTR. Stable and unstable equilibria and limit cycles were identified for single reactions, while quasiperiodic, periodic, and chaotic attractors were recognized for two consecutive reactions. Cho et al. [21] used this theory to study the stability of a hypothetical, autocatalytic, and nonisothermal CSTR and outlined its behavior through qualitative bifurcation diagrams. More recently, Ball and Gray [22], Elnashaie and Elshishini [23], and Ajbar and Alhumaizi [24] used this theory to study the stability of 2.3-epoxy-1-propanol hydrolysis, reactions involved in the petrochemical industry (e.g., benzene oxidation and ethylene hydrogenation) and enzymatic reactions, respectively, predicting different behaviors around the arisen bifurcations and topologically discriminating the generated regions. On the other hand, the bifurcation theory characterizes the topological structural changes of a system, which is represented by a set of equations (e.g., ordinary differential equations). The locus where such changes arise is called a bifurcation and can be classified according to its mathematical properties (e.g., the nature of Jacobian matrix eigenvalues and the normal form). Their behaviors can be verified, that is, by constructing bifurcation diagrams (state variables versus bifurcation parameter) and phase diagrams. Guckenheimer and Holmes [25] proposed a relevant theoretical framework on the principles of bifurcations and placed special emphasis on the study of nonlinear, quasiperiodic, and chaotic oscillations. Kubíček and Marek [26] proposed very useful mathematical tools and numerical techniques to determine the different types of bifurcations and to construct the solution curves for equilibrium points and limit cycles (i.e., using continuation methods). Some of their algorithms have been implemented in software such as Fortran and C language. Kim et al. [27] used this theory to demonstrate the occurrence of period doubling bifurcations, irregular oscillations, and chaos in a system of two chemical reactions in series carried out in a tubular reactor. Troger and Steindl [28] contextualized this theory from an engineering point of view and applied it to solve mechanical problems. More recently, Strogatz [29] proposed a systematic development of this theory. It covers essential topics such as first-order differential equations and their bifurcations, analysis in phase planes, limit cycles and their bifurcations, chaos, and strange attractors. Elnashaie and Elshishini [30] explicitly used this theory in real cases of heterogeneous chemical reactors at the petrochemical industry. Govaerts [31] detailed numerical methods for solving problems with bifurcations, including improvements to equilibrium continuation algorithms (e.g., arc length continuation). Kim et al. [32] constructed bifurcation and phase diagrams to explain the highly nonlinear behavior of a bioreactor with immobilized enzymes, considering diffusion and substrate inhibition phenomena. Kim and Chang [33] tested the chaotic behavior of two consecutive first-order irreversible reactions carried out in a nonisothermal CSTR, using bifurcation diagrams and phase portraits. Furthermore, its dynamical model was validated by an artificial neural network method. Wiggins [34], Kuznetsov [35], and Seydel [36] have shown the evolution of bifurcation theory in recent decades, based on excellent graphical patterns of system behavior around different types of bifurcations. Kuznetsov’s work stands out, since in a very peculiar and elegant way the mathematical principles of this theory were synthesized, introducing important concepts such as local and global bifurcations, codimension of a bifurcation, bifurcations in continuous and discrete times, and continuation methods of equilibria and limit cycles. In addition, he developed the MATCONT’s graphical interface, available for MATLAB®, which incorporates numerical continuation methods for the construction of bifurcation and continuation diagrams. Izhikevich [37] adopted the notions of this theory and applied them to neural problems. He also proposed adequate and very assertive diagrams about the appearance of different types of bifurcation. Ajbar and Alhumaizi [24] used this theory to study the dynamic behavior of chemostats (enzymatic reactions) with stability problems. It was possible to identify different bifurcations and multiplicity phenomenon characteristic of these reactions. Finally, in the works of Ball and Gray [38], Ojeda Toro et al. [5, 6], and Gómez García et al. [39, 40], both theories were synergistically combined. In these, the study of thermal stability and dynamic behavior of specific reactive systems was systematically addressed, such as hydrolysis of methyl isocyanate, acid-catalyzed hydrolysis of glycidol, acetic anhydride hydrolysis, and hydrogen peroxide decomposition. All possible bifurcation diagrams and two parameter continuation diagrams were constructed. Therefore, the regions of stability (safe operating conditions) and instability could be topologically discriminated for each case, and the dynamic behavior was verified by computed time series (sensitivity analysis). It should be noted that the bifurcation and dynamical analyses and the aforementioned mathematical tools can be considered as key strategies for the design and control of a nonisothermal reactor with a priori stability problems for its large-scale commissioning and operating start-up.

In this study, a systematic methodology was used to study the dynamic behavior of acetyl chloride hydrolysis. It was based on the combination of the principles of the singularity and bifurcation theories and was supported by continuation methods incorporated in the MATCONT (MATLAB®) package. The experimental results reported by Baccaro et al. [41] were initially reviewed. Therefore, both reaction kinetics of the reaction and the characteristics of the reactor and experimental conditions were adopted. They demonstrated that the hydrolysis of acetyl chloride, in the presence of acetone as solvent, presents thermal oscillations in a continuous stirred tank reactor (CSTR). Here, the characteristics of the CSTR were considered to derive an appropriate dimensionless mathematical model. To examine in detail the topological features of the acetyl chloride hydrolysis, different bifurcation and continuation diagrams were mapped and analyzed. As a result, stable (safe) and unstable operating conditions were defined.

2. Kinetics of Acetyl Chloride Hydrolysis

The acetyl chloride hydrolysis can be expressed as follows:

The background on its reaction rate law and its experimental conditions were reported elsewhere [4143]. Briefly, a 2-liter polyester CSTR was used to carry out the hydrolysis reaction. The CSTR was continuously fed by two liquid currents: the first, acetyl chloride with a volumetric flow rate of 27.5 cm3/min, a concentration of 4.21 M and a temperature of 284 K, and the second, a mixture of 8.33 to 10 wt.% of water (5.06 M) in acetone with a controlled flow rate of 64 cm3/min and a temperature of 284 K. In addition, the reactor included cooling coils that also functioned as baffles. A cold mixture of 10% methanol, 25% ethylene glycol, and 65% water by volume was pumped through the cooling coils with a flow rate of 200 cm3/min and a temperature of 270 K. The temperature change of the refrigerant was negligible due to its high volumetric flow rate. Subsequently, combining the mass and energy balances and using nonlinear regression, kinetic parameters were linear fitted to the following reaction rate law: where the subscripts A and B correspond to acetyl chloride and water, respectively. The experimental conditions as well as the values of the fitted kinetic parameters ( and ) are presented in Table 1.

3. Dynamic Model for the Thermal Stability Analysis

For the thermal stability analysis, a homogeneous and un-steady-state CSTR model was considered. The start-up of the reactor corresponds with a CSTR of fixed volume, with constant overflow () and constant density. Furthermore, all thermal and composition gradients were neglected. Therefore, the material and energy balances can be represented as follows: where .

Differential equations (3) and (4) are restricted to the following initial conditions:

Equations (3) and (4) were reparametrized, correlating the selected parameters with those of the original data reported by Baccaro et al. [41]. Thus, the following dimensionless variables and parameters were introduced:

Therefore, the dimensionless dynamic model can be expressed as follows: with the following initial conditions:

Thermal stability analysis for the acetyl chloride hydrolysis can be performed by simultaneously solving equations (7) and (8) and varying the parameters (, , , , , and ). Experimental conditions, from Baccaro et al. [41] (Table 1), were used as the starting point of the numerical routine. They are listed in dimensionless form in Table 2.

4. Characteristics of the Dynamic Behavior

The mathematical model was solved by an iterative numeric method using MATCONT (MATLAB®) package. At first, the steady state was determined over a long integration time. It represents the starting point of the numerical continuation for the subsequent construction of the bifurcation and continuation diagrams in different conditions. Singularity and bifurcations theories were combined to locate and interpret the geometric points in a multiparameter space. The eigenvalues of the continuation points have a direct correspondence with the stability of the system (the appendix shows the relationship between the dynamic reactor model and the eigenvalues of the characteristic or Jacobian matrix), as shown below [23, 35]: (1)Stable Points. These can be classified as stable nodes or stable foci. The former is present when the eigenvalues of a solution are negative real numbers. Meanwhile, the second ones appear when the eigenvalues of a solution are a pair of complex conjugates with negative real part. When the reactive system is operated in a stable node and a disturbance is applied to it, the system variables (conversion and temperature) quickly tend to return to the conditions of the stable node. On the other hand, when operating on a stable focus and a disturbance is applied to it, the system variables present oscillatory behavior that decreases in amplitude with time until reaching the conditions of the stable focus(2)Unstable Points. These can be classified as unstable nodes, unstable foci, and unstable saddles. The first arises when the eigenvalues are positive real numbers. For their part, the seconds appear when the eigenvalues are a pair of complex conjugates with a positive real part. And the third parties arise when a pair of eigenvalues are real numbers with opposite sign. Contrary to stable points, when the reactive system operates in an unstable node, the conversion and temperature tend to deviate rapidly from these conditions towards a more stable region. This trend persists when the reactive system operates in an unstable focus, with the variables deviating from this condition in an oscillatory manner and increasing their amplitude over time. However, if the system operates near a saddle point (generating an imaginary separatrix of stable regions), the conversion and temperature tend rapidly towards the nearest stable region(3)Fold Bifurcation (FB). In many references, it can also be found as static limit point, turning point, or saddle node. Its occurrence must satisfy the following condition: a real eigenvalue must change sign (e.g., a positive real eigenvalue becomes negative or vice versa). This is the point responsible for the appearance of the hysteresis phenomenon. If a perturbation is applied when the reactive system operates near this bifurcation, the conversion and temperature tend to jump across the separatrix generated by a saddle point until they again reach a stable point (stable node or stable focus). The occurrence of this phenomenon can be highly risky in the operation of highly reactive systems(4)Hopf Bifurcation (HB). In several references it can be also found as the Andronov-Hopf bifurcation. In this case, the following condition must be met: the real part of a complex conjugate eigenvalue pair changes sign (e.g., negative real part of a complex conjugate pair changes to positive real part). When the reactive system is operated in the vicinity of this bifurcation, it is exposed to both conversion and temperature exhibiting continuous oscillatory behavior (with fixed amplitude). This region is also considered out of control and highly dangerous for highly exothermic reactive systems(5)First Periodic Orbits or Limit Cycle Oscillation. These arise just after a Hopf bifurcation appears, that is, when the eigenvalues are pure imaginary numbers. This condition is associated with the first cycle or oscillatory behavior that forms after the appearance of the Hopf bifurcation. In other words, it corresponds to the initial cycle with lower emerging amplitude(6)Generalized Hopf (GH) Bifurcation. It can be also found as Bautin bifurcation. It arises when a Hopf bifurcation appears (i.e., it has a pair of pure imaginary eigenvalues) and the first Lyapunov coefficient has a value of zero. This type of bifurcation occurs when two Hopf bifurcations (one supercritical and one subcritical) and a fold bifurcation of limit cycles collide. When the reactive system operates near this bifurcation, both the conversion and temperature are exposed to a hysteresis phenomenon, where the separatrix is an unstable limit cycle and the stable solutions can be a stable focus and a stable limit cycle (described above). It should be noted that under these conditions, the operation of exothermic reactive system is highly unsafe(7)Bogdanov-Takens (BT) Bifurcation. It can also be found as double zero bifurcation. It arises when a pair of eigenvalues has zero values. This type of bifurcation occurs when twofold bifurcations, a supercritical Hopf bifurcation, and a homoclinic saddle orbit collide simultaneously. When the reactive system operates in the vicinity of this point, the variables are exposed to undergo the hysteresis phenomenon, where the separatrix is a saddle point, and the solutions can be two stable foci or an oscillatory envelope (stable limit cycles). Again, to ensure a stable and safe operation of exothermic reactive systems, it is important to avoid these types of conditions(8)Cusp Bifurcation (CB). This arises when the following conditions are met: an eigenvalue has zero value, and the second derivative of a state function with respect to a state variable has zero value (this restriction is also known as the center manifold). This bifurcation arises when twofold bifurcations collide. Consequently, the reactive system in these conditions can be exposed to multiple steady-state regimes (hysteresis phenomenon), which makes it a highly unsafe operation in the vicinity of this point

For each dimensionless parameter, dimensionless temperature () and conversion () diagrams can be obtained as a function of dimensionless time (), where critical points and stable/unstable regions are adequately discriminated. Likewise, continuation diagrams can be mapped by simultaneously varying two parameters.

5. Results and Discussion

As a first stage in the numerical analysis, the thermal profile for the hydrolysis of acetyl chloride was calculated using the reaction rate law experimentally verified by Baccaro et al. [41]. Figure 1(a) reveals a hot spot at and and a sustained cyclical behavior before it approaches steady state. From another perspective, a temperature-conversion phase plane was plotted (Figure 1(b)). Notice that it exhibits spiral trajectories before reaching steady state.

Starting from steady-state conditions (i.e., and ), dimensionless temperature and conversion bifurcation diagrams were calculated as function of , , and parameters. For each case, the envelop of periodic solutions, Hopf bifurcation points, and the stable and unstable conditions are shown in Figure 2. Each diagram can be divided in three regions as follows: (i)The first one corresponds to the ranges of ( cm3/min), ( mol/cm3), and ( K): in this region, global stability is achieved since there is only one stable solution. Conversion and dimensionless temperature profiles decrease with increasing values of values, from (adiabatic case) till (Figures 2(a) and 2(b)). The opposite effect was observed with the heat released and cooling temperature (e.g., conversion and dimensionless temperature increase with an increase of and , Figures 2(c)2(f))(ii)The second region corresponds to the following ranges: ( cm3/min1), ( mol/cm3), and (): this region begins with a supercritical Hopf bifurcation (at , , and ). From this point on, the state variables exhibit high sensitivity to fluctuations in operational conditions. Indeed, the state variables present a high sensitivity to fluctuations in operating conditions (i.e., they become unstable and a branch of periodic solutions emerges. Only the maximum and minimum values of cyclic orbits are shown in Figure 2). The maximum values of conversion and dimensionless temperature reach the following values: and ( K) for ; and ( K) for , and and ( K) for . In the first two cases (i.e., or ), temperatures are higher than the boiling points of single reactive species: acetyl chloride ( K or ), water ( K or ), acetic acid ( K or ), and hydrochloric acid ( K or ). Therefore, the reactive system can undergo a phase transition (liquid-vapor) within an oscillating thermal envelope from and parameters. Consequently, the pressure in the reactor can increase to the point of causing an explosion. Finally, in this region, the oscillatory dynamic behavior prevails up to a second supercritical Hopf bifurcation (at , , and )(iii)In the third region, which corresponds to the values of ( cm3⋅min-1), ( mol⋅cm-3), and ( K), steady-state stability is restored. Note that as the values of and increase, dimensionless temperature also increases proportionally, reaching quite high values

The dynamic behaviors generated from the variation in the water feed ratio () and the feed temperature () are similar to those explained in Figure 2. They were included as supplementary material (Figure S1a-d).

If ideally all operating parameters could be fixed (strict control) except for one (one degree of freedom) and the objective of the operation was to generate the greatest amount of acetic acid and hydrogen chloride in a stable manner, then the operating parameters should be the following: ( cm3/min1), ( mol/cm3), and ( K). On the other hand, if the idea is to prevent the chemical reaction from occurring or to control it safely, the operating parameters should be the following: ( cm3/min1), ( mol/cm3), and ( K). However, any real process could be exposed to the simultaneous disturbance or change of more than one operating parameter. For this purpose, the analysis of the system dynamics is performed using a continuation diagrams, as demonstrated below.

Bifurcation diagrams, as a function of two parameters, also called two-parameter continuation diagram (TPCD), is another way of presenting the dynamic behavior of reactive systems [5, 6, 23, 39, 40]. TPCD is especially useful if two codimension bifurcations arise (this type of bifurcation can be observed if two parameters are perturbed simultaneously) [35]. Figures 3 and 4 show the possible TPCD combinations for acetyl chloride hydrolysis (additional figures for - and - are included as supplementary material, Figure S2 a and b).

The region enclosed by the Hopf bifurcation curve represents the oscillatory behavior of the state variables. Each pair of the Hopf bifurcation loci is connected by a periodic orbit with different amplitudes. In other words, in this region, both conversion and dimensionless temperature experience a sustained cyclical behavior (similar to those shown in Figure 2), with the amplitude of the oscillation changing as the parameters described above are modified. Beyond the limits of the Hopf curve, a transition regimen (stable to cyclical envelope) can be expected. Likewise, the limit point curve defines the region of multiplicity of stationary states (multiple solutions). As mentioned above in Section 4, in this region, the system variables are prone to operate in multiple steady-state regimes, causing abrupt jumps in them even when the reactive system experiences small perturbations in its parameters.

Additional topological features of this reactive system are evident in the TPCDs. They include the presence of generalized Hopf (GH), cusp bifurcation (CB), and Bogdanov-Takens (BT) bifurcation points. The occurrence of a GH bifurcation (the first Lyapunov coefficient disappears of the Hopf bifurcation) greatly affects the qualitative topology of the continuation diagram and therefore the dynamic behavior of the state variables of the system. From its coordinates (point marked as 0 in Figure 5), two separate branches emerge (subcritical (H+) and supercritical (H-) Hopf bifurcations corresponding to the Hopf bifurcation points with positive and negative Lyapunov coefficient, respectively). They generate a separatrix (which divides the vs. plane in two regions marked as 1 and 2 in Figure 5).

To disaggregate the bifurcation diagram, the quantitative and qualitative behavior of the reactive system in the vicinity of the GH bifurcation point within the multiparametric space was examined, for instance, a locus in the region 1 where the system presents a stable focus in the absence of limit cycles (Figure 5, plot 1, H-). Crossing the separatrix from region 1 to region 2 infers the appearance of a unique stable limit cycle (Figure 5, plot 2, H+). Note that in the conditions of the GH point, an extra unstable cycle is evident (Figure 5, plot 0). This implies that, for parameter values close to the GH point, the system presents two limit cycles, which collide and disappear through a saddle-node bifurcation of periodic orbits. However, for this particular case, the saddle-node bifurcation (fold or turning point) of the cycles overlaps the HB curve with great precision, making the missing region to be imperceptible.

On the other hand, two branches of the fold bifurcation curve (T1 and T2) meet tangentially at point CB (marked as 0 in Figure 6) dividing the bifurcation diagram into two regions (marked as 1 and 2 in Figure 6). Make a round trip from region 2 where limit cycles are not possible (i.e., stable node can be reached when the reactor operates in this region, Figure 6, plot 2). Then, the node moves to another coordinate on the T1 component of the fold curve (Figure 6, plot T1). Entering region 1 produces other stable steady-state nodes that divide the phase plane into two attraction domains (Figure 6, plot 1). These roots or equilibrium values merged and disappeared at point CB. If we continue the clockwise journey, the T2 component of the fold curve reverts to stable single-node behavior (Figure 6, plot T2). To complete our round trip, notice that a stable cycle limit is present at the coordinates of point CB (plot 0 in Figure 6). In summary, the incidence of CB point implies the presence of hysteresis or multiplicity phenomena in steady state.

To complete the analysis, the dynamic behavior around the BT bifurcation point (marked as 0 in Figure 7) was evaluated. It separates two branches T- and T+ from the fold curve. Notice that the Hopf bifurcation curve (H in Figure 7) is tangential to the fold curve. When we move clockwise around the coordinates of the bifurcation point BT in the parameter plane, four regions can be distinguished. Regions 1 and 2 are enclosed by the loci of limit points. Two equilibria (saddle point and a stable node) are predicted in these regions, colliding and vanishing via fold bifurcation curve (plots 1 and 2 in Figure 7). The stable node becomes a stable focus and loses stability after crossing the H boundary of the Hopf bifurcation. A stable limit cycle is present in regions 3 and 4 (plot 3, 4, T+ in Figure 7). Going back to point 0 produces a node transition (plot 0, T-, H in Figure 7).

In summary, for a safe operation of the reactor (e.g., to avoid oscillations of the state variables and the hysteresis phenomena), the following operating ranges must be avoided: (i) ( cm3/min1) and ( K) (Figure 3(a))(ii) ( cm3/min1) and ( mol/cm3) (Figure 3(b))(iii) ( cm3/min1) and ( K) (Figure 3(c))(iv) ( mol/cm3) and ( cm3/min1) (Figure 3(d))(v) ( mol/cm3) and ( mol/cm3) (Figure 4(a))(vi) ( mol/cm3) and ( K) (Figure 4(b))(vii) ( K) and ( mol/cm3) (Figure 4(c))

Note that the analysis of each parameter allows for a comprehensive risk assessment of the system. This defines the safe and unsafe regions of reactor operation and characterizes the atypical or complex behavior of the reactive system.

At this point in the discussion, an important question for practical applications arises: what is the location of operating points where the reactor could normally be operated? Considering that the hydrolysis of acetyl chloride is carried out at high temperatures, the operating conditions must be set other than the practical limit of stability (i.e., they must be chosen based on low conversion). However, choosing a safe operating point is not easy because very small differences between parameter values can result in large differences in dynamic behavior and ultimately in steady-state quality. Experimental data and simulations of dynamic behavior have shown that the oscillations shift the reactor to a lower steady state. Additional data on reactor performance can be found in bifurcation diagrams involving dimensionless reactor temperature and acetyl chloride conversion (Figure 2). Operation within these multiple steady-state regions bounded by turning points and oscillation regime regions can take place only with a high flow rate of acetyl chloride, which is not common. In the case of reduced feed flow (e.g., rapid or steady decrease in acetyl chloride flow rate), conversion and temperature will fall following the branch of steady states. No regime of multiplicity or oscillation is located there. Conversely, an increase in the flow rate of the cooling medium will increase heat removal by changing, through branch shifting, the operation point from a region of sustained oscillations to a lower steady state characterized by low acetyl chloride conversion. The operator must be aware of these situations in the event of a change in acetyl chloride flow rate and/or refrigerant flow rate.

6. Conclusions

This work has made it possible to mathematically verify the behavior of the hydrolysis of acetyl chloride. Based on the bifurcation theory, applied to the multiparameter dynamic system, safe operating conditions were identified. Bifurcation and multiparameter continuation diagrams were constructed. All dynamic and critical states were located using MATLAB® software. The oscillatory behavior of the variables and the multiplicity of steady states were demonstrated, in a wide range of parameters. Furthermore, the quantitative dynamic behavior around two codimension bifurcations was predicted satisfactorily. Therefore, the following conditions are recommended for safe operation: low values of (<0.0330), (<0.0452), (<0.04869), and (<1.2439) and high values of (>3.1940), that is, in dimensional way,  mol/cm3,  K,  K, and  cm3/min1.

Appendix

The dynamic model of the CSTR (equations (7) and (8)) can also be expressed as follows:

They can be linearized considering only the first-order terms of the Taylor series as follows: where the subscript SS denotes at steady-state condition.

From equations (A.1) and (A.2), it can be inferred that

When equations (A.5) and (A.6) are replaced in (A.3) and (A.4) and ordered, two coupled homogeneous linear differential equations are obtained where and are defined as deviation variables.

In matrix form, they can be expressed as where

is the vector of deviation variables and is the Jacobian or characteristic matrix of the system.

The differential equation system (A.8) has the flowing analytic solution:

Using the Sylvester theorem, it is possible to express (A.10) as a function of the eigenvalues () of the Jacobian matrix [44], as shown below: where is the identity matrix.

Therefore, the solution vector equals to

It is clear that the behavior of depends on the nature of the eigenvalues. These are obtained from the following equation [45, 46]:

The determinant solution is a λ polynomial (characteristic equation):

Note that the degree of this polynomial will depend on the dimensional size of the system; it is grade two for the present case.

By solving the quadratic equation, the eigenvalues are equal to where

Thus, the eigenvalues depend on the trace values (tr) and the determinant (det) of Jacobian matrix . In turn, these will depend on the operational parameter values of the original system, according to equations (7) and (8).

Notation

Latin Symbols
:Frequency factor (cm3⋅mol-1⋅min-1)
:Cooling coil area (cm2)
:Reacting concentration (mol⋅cm-3)
:Specific heat of the reaction mixture (cal⋅mol-1⋅K-1)
:Activation energy (cal⋅mol-1)
:Reaction rate (mol⋅cm-3⋅min-1)
:Ideal gas constant (cal⋅mol-1⋅K-1)
:Time (min) or temperature (K)
:Residence time (min)
:Overall heat transfer coefficient of the coil (cal⋅cm-2⋅min-1⋅K-1)
:Reactor volume (cm3)
:Conversion.
Greek Symbols
:Dimensionless frequency factor
:Dimensionless heat transfer term
:Reaction enthalpy (cal⋅mol-1)
:Dimensionless temperature
:Density (g⋅cm-3)
:Dimensionless time
:Volumetric flow (cm3⋅min-1)
:Dimensionless heat released term.
Subscripts
:Initial condition
:Acetyl chloride
B:Water
:Cooling condition
:Feed condition.

Data Availability

The data associated with this research are available from the corresponding author upon request.

Disclosure

Juan Carlos Ojeda Toro actual address is Grupo de Investigación de Aprovechamiento Tecnológico de Materiales y Energía–GIATME, Unidad de Ciencias Básicas, Facultad de Ingeniería, Universidad ECCI, Medellín, Colombia.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this article.

Acknowledgments

Juan Carlos Ojeda Toro was a beneficiary of a COLCIENCIAS grant (Programa Beca Doctorados Nacionales, Convocatoria 647 de 2014). The authors gratefully acknowledge the financial support of the Universidad Nacional de Colombia through the projects HERMES-51167, HERMES-51225, and HERMES 55259 (Convocatoria para el Fortalecimiento de la Investigación, Creación, e Innovación Articulado con la Formación en la Universidad Nacional de Colombia 2020-2021).

Supplementary Materials

Fig. S1: conversion and temperature dynamic behavior of the acetyl chloride hydrolysis. (a and b) and (c and d) . The other parameter values are maintained constant as in Table 1. Fig. S2: multiparameter continuation diagrams showing the critical points. (a) and (b) . The other parameter values are maintained constant as in Table 1. (Supplementary Materials)