Modeling of Isocyanate Synthesis by the Thermal Decomposition of Carbamates

: The presented work is devoted to isocyanate synthesis by the thermal decomposition of carbamates model. The work describes the existing isocyanate-obtaining processes and the main problems in the study of isocyanate synthesis by the thermal decomposition of carbamates, which can be solved using mathematical and computer models. Experiments with carbamates of various structures were carried out. After processing the experimental data, the activation energy and the pre-exponential factor for isocyanate synthesis by the thermal decomposition of carbamates were determined. Then, a mathematical model of the reactor for the thermal decomposition of carbamates using the COMSOL Multiphysics software was developed. For this model, computational experiments under different conditions were carried out. It was shown that the calculation results correspond to the experimental ones, so the suggested model can be used in the design of the equipment for isocyanate synthesis by the thermal decomposition of carbamates.


Introduction
Isocyanates are currently one of the most in demand products in the chemical industry, since they are used as raw materials in polyurethane production [1]. Polyurethanes are used in many industries, including construction; in the manufacture of automotive parts; and as insulation materials, paints, varnishes, adhesives, fillers for upholstered furniture, etc. [1]. Isocyanates are also valuable intermediates in the fine organic synthesis of pesticides and other biologically active substances. However, isocyanate production is associated with a number of problems because intermediates are highly toxic and explosive substances, and, therefore, ensuring industrial safety requires great investment.

Isocyanates: An Overview of Preparation Methods
The phosgenation of amines is the most popular method for isocyanate synthesis [2]. The toxicity of phosgene and the risk of its leakage, however, create a number of disadvantages in this method. Moreover, phosgenation generates a large volume of hydrochloric acid, which causes the corrosion of the equipment, thereby increasing the maintenance costs. Notably, the inconvenience of transporting phosgene due to its toxicity also limits the applicability of the method, especially in small-scale production. The global trend of green chemistry calls for alternative phosgene-free methods for producing isocyanates [3].
The thermal decomposition of carbamates is a promising method for isocyanate production as an alternative to the phosgene method. Methods for the production of industrially important aliphatic monoisocyanates with low boiling points are prioritized. These isocyanates are consumed in large volumes and have a very wide range of application.
The currently existing methods for isocyanate production by the thermolysis of carbamates can be divided into gas and liquid phases.

Gas-Phase Methods for the Isocyanate Synthesis by Thermal Decomposition of Carbamates
The process of the thermolysis of carbamates (urethanes), esters of carbamic acid, in the gas phase is widely described in patent literature [4][5][6][7][8][9][10]. Urethane degradation can be carried out in a wide range of operating pressures and temperatures, with [8,[11][12][13] or without catalysts [14,15], according to the patent data. Almost any carbamate synthesized from thermally stable alcohols and amines can be degraded to produce isocyanates. The thermolysis of carbamates in the gas phase is carried out at temperatures of some 400 °C. It has a number of disadvantages, including a high amount of byproducts, which leads to contaminated devices, high equipment requirements, and increased energy consumption.

Liquid-Phase Methods for the Isocyanate Synthesis by Thermal Decomposition of Carbamates
In order to lower the process temperature and reduce undesirable side reactions, liquid-phase methods for the isocyanate synthesis by the thermal decomposition of carbamates were developed.
A review of patent sources shows that there are several directions for the research of the optimal process in the liquid phase. For instance, the thermolysis of urethanes can be carried out in high-boiling inert solvents [16][17][18][19][20][21]. This research direction is divided into several categories according to the type of catalyst used for thermolysis. Thermolysis in those solvents that exhibit catalytic activity during carbamate thermolysis looks promising [18]. The development using ionic liquids as a solvent, which is described in [22,23], is also of interest.
We studied patented methods that used high-boiling solvents and came to the conclusion that, in most cases, aliphatic or aromatic hydrocarbons are used as solvents, and transition metals and their compounds are used as thermolysis catalysts [24][25][26]. For example, [24] proposes a method for isocyanate production by the thermal decomposition of urethanes in a high-boiling solvent with a wide range of metals or metal compounds, such as manganese, molybdenum, tungsten, vanadium, iron, cobalt, chromium, copper, and nickel, as a catalyst.

Mathematical Modeling: An Overview of Software Products
Based on the review of the literature discussed above, gas-phase non-catalytic thermolysis is the most industrially applicable phosgene-free method of synthesizing isocyanates from carbamates. In order to shift the equilibrium of the main reaction and synthesize a purer isocyanate, the process should be carried out with an inert carrier gas or under vacuum. It is easier to control an inert gas than a vacuum along the entire reactor length, and an inert gas does not depend on the hydraulic resistance of the reactor. Carrier gas processes, apart from being easily scaled up, can reduce the partial pressure of carbamates in the gas phase, which helps to lower their boiling points. This ensures an easy and complete evaporation of carbamates upon entry into the thermolysis reactor, which decreases the amount of by-products.
The development of new efficient processes and optimization of existing technological facilities provides the development of chemical technology. Therefore, computer modelling realized by running specialized programs that provide a highly accurate mathematical description of the processes under study is of interest at the present stage of chemical technology development. Modern approaches to mathematical and computer modeling are widely used both in the study of a specific process and production development. Mathematical modeling makes it possible to study the properties of objects on mathematical models, predict optimal schemes and operating modes of commercial plants, and develop automated control systems for technological processes.
Mathematical modeling also allows the scaling up of the process and optimizing critical parameters at the design stage of a semi-commercial plant under non-isothermal conditions in the reactor [27].
Such computer models are based on: − heat and material balances; − component composition; − the nature of the interacting substances.
When studying any chemical technological process, material balances are introduced in the form of differential equations that represent the law of conservation of matter.
Heat balance includes the thermodynamic properties of the components, such as the density and temperature at the boiling point or under normal conditions, as well as the heat of combustion and formation, viscosity, etc.
This simulation system usually consists of: Those methods, widely described in [28,29], allow solving a significant part of the issues with chemical technological processes.
Modelling helps to study the following processes [30]: This paper studies methods for the synthesis of three aliphatic isocyanates by the thermal decomposition of carbamates, researches the kinetics of the thermolysis, and proposes a mathematical model of the process developed with the COMSOL Multiphysics software. In addition, this paper assesses the efficacy of the model and establishes its limits.

Statement of the Problem
The main problem in the study of isocyanate production by the thermal decomposition of both aliphatic and aromatic carbamates is the absence of the described kinetic regularities. However, there are a large number of patents on this topic, including descriptions of pilot plants and their operating modes. The latter suggests that this area is promising, because the thermal decomposition of carbamates allows us to abandon the use of highly toxic phosgene and to avoid problems associated with the resulting hydrogen chloride.
The main goal of this work was to develop a universal method for obtaining isocyanates from carbamates using modern modeling methods, with the possibility of the further integration of this technology into enterprises of the low-tonnage production of biologically active substances from groups of pesticides, fungicides, insecticides, and medicines. Equipment using this technology can be installed in close proximity to the production facilities discussed above.
To reach this goal, we created a laboratory plant for the thermal decomposition of carbamates, in which the displacement reactor was the most important unit. We carried out a large number of experiments with carbamates of various structures and studied the dependence of the initial substance conversion on the directly controlled parameters-the temperature and flow rate of the carrier gas, which determine the residence time in the reactor. After processing the experimental data, the reaction rate constant and the pre-exponential multiplier were determined. Then, we created a mathematical model of the carbamate thermal decomposition reactor in accordance with the characteristics, determined the process conditions, and calculated the conversion at the specified values of variables. At the next stage, the adequacy of the model was checked, which consisted of checking the experimental and calculated values using the Fisher statistical method, the positive result of which indicates the sufficient reliability of the model. It is important to highlight that the proven model can be used in the design of the equipment.
An approach for determining the reaction rate constant and the pre-exponential multiplier was proposed. A model of a carbamate thermal decomposition reactor was developed using the computational fluid dynamics (CFD )method realized in the COMSOL Multiphysics software. The adequacy of the developed model was verified by comparing the experimental and calculated values of the transformation degree. It is important to highlight that the proven model can be used in the design of the equipment.

Experimental Studying
The research objects are O-methyl-N-(n-butyl) carbamate, O-methyl-N-cyclohexyl carbamate, and O-methyl-N-benzyl carbamate, which are the raw materials for the synthesis of industrially significant aliphatic isocyanates n-butyl isocyanate, cyclohexyl isocyanate, and benzyl isocyanate.

Methods of Conducting Experiments
The chemical reaction network for a kinetic study is shown in Figure 1. The reaction involves the interaction of a primary amine A and an electrophilic agent B (dimethyl carbonate) with the production of the required carbamate C, which is then thermally decomposed in a displacement reactor at a temperature of between 250 and 600 °C. Then. the reaction mass containing the target isocyanate D and the unreacted carbamate C as the main products is absorbed into a sorption solution containing N-methyl-N-benzylamine E and 3,5-dibromopyridine as an internal standard. As a result of the reaction between isocyanate D and N-methyl-N-benzylamine E, asymmetric urea F is produced. In order to improve the separation of peaks in HPLC analysis, a molar excess of amine E with respect to carbamate C binds with acetic anhydride G, forming the corresponding acetamide H during sample preparation. The resulting sample can be sent for analysis in order to obtain the molar ratio of the reaction products. All the conditions of analytical measurements were determined for a qualitative and quantitative analysis of the resulting mixtures of products after thermolysis. Moreover, the design strategy for a thermal decomposition laboratory plant was chosen.

Methods of Analytical Research
The samples were analyzed using the Shimadzu LC20 Promince high-performance liquid chromatography system (LC-20AD pump; MZ PerfectSil Target C18 250 mm × 4.6 mm, 5 microns column; RID-20A refractometric detector and SPD-20A UV detector installed in series).
Conditions for the chromatographic measurements:

Description of the Laboratory Plant
A laboratory facility was constructed, since the gas-phase thermolysis of carbamates was chosen. The facility was based on a reactor, with the hydrodynamic regime closest to that of a plug-flow reactor. First, the main units of the experimental facility were determined, as shown in 1. Feed unit for initial reagents; 2. Preheating system; 3. Reactor unit; 4. Sorption unit; 5. Facility control system. A laboratory facility was developed to carry out thermolysis; the diagram is shown in Figure 3. The facility consists of a feed unit for initial reagents, a preheating system, a thermal decomposition reactor unit, a sorption unit for reaction products, and a facility control system. The characteristics of the most important components of the facility are shown in Table 1. A thermostat (TS) is the basis of the preheating and thermolysis units. The thermostat is a heat-insulated box of 250 × 300 × 330 mm made of stainless steel AISI304 (08X18H10). This section is divided into 2 thermal zones according to its function. From outside, the feed units for the initial reagents are connected to the thermostat on one side and the elements of the sorption unit on the other.

The Results of the Experiments of the Carbamates Thermolysis
A series of experiments were carried out for the decomposition of N-alkyl-O-methyl carbamates in a displacement reactor in the gas phase to obtain experimental data in a wide range of temperatures and residence times for various carbamates-O-methyl-N-benzyl carbamate, O-methyl-N-(n-butyl) carbamate, and O-methyl-N-cyclohexyl carbamate-which allow us to calculate the kinetic parameters of the process. Experiments with a carrier gas flow rate of 0.05 mL/min over the entire temperature range were mandatory for all three carbamates. The remaining experiments were conducted within the selected zone in Table 2 for different carbamates. Depending on the degree of transformation obtained, the experiments could vary or be duplicated. Carrier gas flow, L/min The conditions of the experiments and the fractional conversions for all the thermal decomposition kinetics experiments are shown in Tables 3-5. The results show the rise in the degree of conversion of the original starting material with the increasing maximum temperature in the reactor and/or the increasing residence time of the reaction mixture in the reactor, while the formation of such by-products as cyanurates and carbodiimides and other possible products was not observed. This fact indicates a high selectivity of obtaining the corresponding isocyanates from various carbamates.

Processing and Analysis of Experimental Data
A greater number of experiments at high temperatures with high fractional conversions were carried out for N-benzyl-O-methyl carbamate, which led to a greater number of experiments in general compared to other carbamates. The data obtained were used later to assess the efficacy of the developed model. In addition, the obtained experimental data were used to calculate the values of the activation energies Eа and the pre-exponential factors k0. To obtain the calculated values for the degree of conversion, the function was set in a software product MathCAD, allowing us to find the value of the activation energy EA and the pre-exponential factor k0 based on the experimental data (temperature profile of the reactor, the carrier gas flow, and the degree of conversion). The function was compiled based on the equation of unsteady mass transfer, with the assumption of the absence of mass transfer by diffusion and back-mixing because of the insignificant contribution. The function also includes non-isothermic character of the temperature profile and the dimensions of the reactor [30].
where is the degree of transformation in the i-th calculation cell; is the the cross-sectional area of the reactor; is the pre-exponential factor; is the activation energy, J/mol; is the gas constant, 8.31 J/(mol*K); ( ) is the the temperature in the reactor, which is described by the function from Section 2.2.3, K; is the total pressure in the reactor, equal to 1 atm; is the the molar flow rate of the gas mixture, mol; is the the degree of transformation determined experimentally; ℎ is the estimated length of the cell, ℎ = ; is the reactor length, m; is the number of splits along the length of the reactor. To obtain adequate results using the analytical model of the reactor, n = 1000 was taken as the required number of splits along the length of the reactor. This number of splits does not require much computing power and allows one to obtain a value with an error of less than 1%.
The results of the function calculations for the three carbamates are in Table 6. Table 6. Activation energies and pre-exponential factors of the studied carbamates. The obtained values were used later as the parameters of the developed model.

Simulation
COMSOL Multiphysics was chosen to study the decomposition process of O-methyl-N-benzyl carbamate, O-methyl-N-(n-butyl) carbamate, and O-methyl-N-cyclohexyl carbamate due to its ease of use and the fact that it is the preferred choice for complex calculations.
The process of the gas-phase decomposition of carbamate highly depends on the temperature regimes used and the hydrodynamic situation in the reactor. Modern computer simulation is the only means by which this conditionality can be ascertained with a high accuracy.
The following assumptions are made for the development of the model:

−
The system considered is a homogeneous, multi-component, and viscous compressible fluid; − This system undergoes a chemical reaction of decomposition, the order of which is equal to one; − The thermal effect of the decomposition reaction is not taken into account due to the low reagent content; − The density of the medium is determined by the density of the carrier gas, argon, due to the low reagent content; − The temperature of the reactor wall is taken as a constant. − The carbamate decomposition process was calculated and studied with the CFD method using the COMSOL Multiphysics software.

Mathematical Model Development
The mathematical model of the process is a complex of differential equations: the conservation of mass over two components (the source substance and the product of the chemical reaction) and over the whole mixture, a momentum conservation equation, and an energy conservation equation. The hydraulic resistance of the column is calculated using the Brinkman equation (the third and fourth components of the right-hand part of the momentum conservation equation). It is assumed that all the equations take into account the column porosity.  The following equations are used to calculate the chemical reaction rate: r = −M kC , r = −r .
The Arrhenius equation is used to calculate the chemical reaction rate constant: The Peng-Robinson state equation was used in the calculation: The reactor model boundary conditions: T(x , y , z , t) = T , T(x , y , z , t) = T , v ⃗(x , y , z , t) = v ⃗ , v ⃗(x , y , z , t) = 0, Y (x , y , z , t) = 0, where ρ is the density, kg/m3; v ⃗ is the rate vector, m/s; T is the temperature, K; p is the pressure, Pa; Y1 is the source substance mass fraction, kg/kgcm; Y2 is the product mass fraction, kg/kg; τ kl is the viscous stress tensor kg/m•s 2 ; C is the heat capacity, J/K; µ is the dynamic viscosity, Pa • s; v is the velocity, m/s; c is the concentration of i product (i = 1 or 2), mol/L; r is the reaction rate, mol/s, Dj is the diffusivity, m 2 /s; k is the reaction rate constant; k0 is the pre-exponential factor; T0 is the inlet flow temperature, K; Tf is the function for temperature on the wall, K; R is the gas constant 8.31 J/mol•K; K is the permeability factor; vm is the molar volume, m 3 /kmol; a, b are the empirical factors conditional on the nature of the substance, and in multi-component systems on the composition; indexes: in-input, w-wall. The "Brinkman Equations", "Heat Transfer in Porous Media", and "Transport of Diluted Species" interfaces were used when setting up the mathematical model.

Virtual Geometry Design and Mesh Generation
This stage of the simulation requires a virtual geometry of the apparatus used for the thermal decomposition of the carbamates. The virtual geometry is shown in Figure 4. The resulting geometry is a hollow cylindrical tube which is bounded from the outside by the surface of the reactor wall and from the inside by the outer surface of the thermocouple sleeve. The outer diameter is 0.01 m and the inner is 0.007 m. The flow input is from the side where L = 0 m and the flow output is from the side where L = 0.14 m.
The geometrical model design is followed by mesh generation. Using COMSOL Multiphysics allowed us to generate combined calculation mesh, consisting of tetrahedral-and prism-shaped integral elements. Moreover, prism-shaped elements are concentrated in the wall region, where the mesh inflation was added for a more accurate modelling of the wall effects. To evaluate the mesh quality, the option "Skewness" is selected. It details the skew of the mesh elements compared to the regular shape (the regular tetrahedron and the regular prism). For the specified parameter, a value of 1 corresponds to the regular shape, and a value of 0 indicates that a 3D element has become a flat 2D shape. The generated mesh is shown in Figure 5. The proposed mesh consists of 156,879 elements. The resulting mesh has an average element skewness of 0.7014. The skewness is not lower than 0.132 relative to all mesh elements. Thus, the quality of the mesh is sufficient for computer modeling.

Temperature Profile Determination on Reactor Wall
Some features that are associated with the design of the experimental plant ( Figure 3) characterize the reactor heating. The reactor zone is directly heated by a separate heating element to high temperatures of up to 600 K. All the pipelines and fittings that are connected to the reactor, the preheating zone, the mixing zone, and the inlet and outlet of the reactor are located in a specially air-heated zone. The temperature at these zones is maintained at 200 °C. The entire setup has a high specific quantity of metal per structure; thus, a significant part of the heat supplied directly to the reactor by its heater is quickly dissipated. Thus, the temperature profile on the reactor wall is determined not only by its heating element, but also by the entire setup as a whole. Predicting such a profile requires an accurate and complex calculation of the heating element, taking into account the entire geometry of the setup, which significantly complicates the task. Therefore, it was decided to determine the temperature profile along the length of the reactor experimentally at different gas flow rates. As a result, a temperature function with the determination coefficient R2 = 99.8 was obtained: where G is the gas flow rate, l/min; l is the reactor length, sm; Tmax is the maximum of temperature, °C.
The coefficients of the obtained relationship are presented in Table 7. The resulting function is used as a boundary condition on the reactor wall when calculating the process using the proposed model. Thus, the simulated temperature profile within the whole reactor volume will be close enough to the experimental one.

Results and Discussion
The computer simulation based on the calculated model aims to predict the course of carbamate thermal decomposition in the previously shown tubular reactor conditional on the process parameters. In the future, the proposed model can be used to increase the process efficiency, to optimize the geometry of laboratory and industrial equipment, and for scaling the process to an industrial level.
The result of each calculation is the fields of the gas flow rate, pressure, temperature, and composition of the system under study within the indicated virtual geometry. A series of calculations was carried out, yielding parameters that correspond to the previously noted experimental data. Based on the results obtained, the efficacy of the proposed model was determined by comparing the calculated and experimental fractional conversions at the reactor outlet.
The input parameters for the computational experiments are presented in Table 8. Carrier gas flow, L/min 0.05-3.0

Temperature Profile Investigation
The obtained function for temperature on the reactor wall (Section 2.2.3) was used for process simulation. The results of comparing the experimental data with the calculated ones are shown in Figure 6. The temperature profiles along the longitudinal section of the reactor, which allow one to see the profile of temperature changes over the entire volume of equipment, with carrier gas flow rates of 0.05, 1.5, and 3.0 L/min at 600 °C are shown in Figure 7.  The COMSOL tools allowed the model-based calculation of data that yielded the calculated temperature profiles and fractional conversions along the reactor length and the profiles of rate and fractional conversions of the thermal decomposition reaction. Figures 8-21 show some of the data received.              The results show an effect that indicates a high fractional conversion in the near-wall region of the reactor. This is due to a low carrier gas flow rate and, as a result, an increase in the residence time of the initial carbamate at those points. Cross-mixing enables the leveling of the concentration of components, and the fractional conversion is leveled at the reactor outlet. In addition, the diagrams represent the regions with the highest rate of the chemical reaction. At low carrier gas flow rates and high temperatures, the same fractional conversion can be observed at the reactor outlet, while the pattern inside the plug-flow reactor can be completely different. This shall be taken into consideration when designing larger reactors, where the uniformity of temperature across the reactor cross-section is crucial. On the curves, it can be seen that the area with a high chemical reaction rate is near the inlet of the reactor and represents a relatively narrow interval in the reactor at low flow rates of carrier gas. However, while the flow rate increases, this zone becomes the zone of maximum temperatures and there is an expansion of the zone, resulting in a decrease in the value of the maximum reaction rate.

Comparative Analysis of Experimental Values of Conversion with the Calculated Data
Based on the data of the simulation experiments with integrated COMSOL Multiphysics tools, the fractional conversions of the carbamate thermal decomposition product were calculated and are shown in Tables 9-11. The values were determined by averaging the fractional conversion at the reactor cross-section at the output at L = 0.14.

Conclusions
A mathematical model for carbamate thermolysis developed with the help of the COMSOL Multiphysics modeling environment was proposed. This model represents the process of the thermal decomposition of O-methyl-N-alkyl carbamates to alkyl isocyanates in a displacement reactor under non-isothermal conditions and can be used to design industrial reactor equipment.
A study was made in a non-isothermal displacement reactor with varying residence times and temperature profiles along the reactor length in a constructed laboratory facility for the thermal decomposition of O-methyl-N-alkyl carbamates. The following kinetic parameters were obtained: the activation energy and the pre-exponential factor of the thermal decomposition of O-methyl-N-butyl carbamate, O-methyl-N-cyclohexyl carbamate, and O-methyl-N-benzyl carbamate, which were used as input parameters of the model during the computational experiments.
The experimental and calculated temperature profiles along the reactor length at different carrier gas rates were compared.
A comparative analysis of the fractional conversion obtained through the experiments and the fractional conversion calculated using COMSOL Multiphysics was carried out. The deviation of the calculated values from the experimental ones did not exceed 7%. The adequacy of the developed model was assessed using Fisher's test. Thus, the proposed mathematical model in the COMSOL Multiphysics modeling environment makes it possible to obtain the values of the process parameters, such as the temperature profiles, flow rates, reaction rates, and component concentrations in the reactor. This can contribute to the development of a new type of equipment for the thermal decomposition of carbamates, since it will partially replace full-scale experiments with computational ones and requires fewer resources during development.
The obtained kinetic parameters and the created mathematical model can be used in the future for modeling and designing pilot plants for the production of biologically active substances from the class of carbamates (Iodocarb, Propamocarb), urea (Benomyl), sulfonylurea (Glibenclamide, Glipizide), and other substances through the phosgene-free method.