New Refrigerant Molecules from Structure Optimization

In the present work, various objective functions were formulated and optimized using the mixed integer nonlinear programming and the generalized reduced gradient nonlinear method from the solver tool of Microsoft® Excel 2016, respectively. The CH3FO2, C2H4F2O, CH2F2O2, CH2F2O, C3H4F2, and the C2H2F2O molecules were found to meet structural feasibility constraints and physical properties from refrigerant molecules and have not previously been reported in the literature. These new refrigerants present global warming potential values similar to that from the R-134a and Freon 12 refrigerants and null ozone depletion potential. Moreover, these molecules are normally flammable, as similar as to R-134a refrigerant. The CH3FO2, C2H4F2O, CH2F2O2, C2H2F2O, and CH2F2O show toxicity values similar to R-134a and Freon 12 refrigerants.


Introduction
The first comprehensive study of small molecule design for refrigerants was carried out by Thomas Midgely, Jr. in 1937, from which the refrigerants called Freons ® were synthesized [1]. By the action of ultraviolet light at the stratosphere, Freon 12 (R-12, CCl2F2) releases a chlorine atom which acts as a catalyst in the decomposition of ozone to produce ClO. This compound reacts with monatomic oxygen from the decomposition of ozone, or nitric oxide, thus producing chlorine atoms and enabling a cycle of chain reactions, with the corresponding depletion of the ozone layer [2]. Other chlorine containing refrigerants such as chlorofluorocarbons and hydrochlorofluorocarbons are also related to ozone depletion, whose emissions are regulated by the Montreal protocol [3]. These molecules were then replaced by environmentally friendly refrigerant molecules, such as R-134a [1]. It is worth noting that refrigerant molecules can cause global warming, which is controlled by the Kyoto protocol [4].
The United States and China are the major refrigerant manufacturers, accounting for most of the corresponding annual production [5]. According to the report issued by Markets and Markets [6], the refrigerant market size was USD 21.3 billion in 2020. It will reach USD 30.8 billion in 2025 due to the growing demand for different applications, such as domestic, commercial, and industrial refrigeration and mobile air conditioning. In Latin America, the refrigerant market led by Brazil, Chile, Argentina, Mexico, and Colombia is expected to grow at an annual rate of 4.1% and reach USD 2.74 billion in 2022 [7]. This growth is attributed to factors such as the increased use of refrigerants in cars, the food industry, refrigeration, and air conditioning for homes. In Ecuador, the refrigerant market is based on refrigerant import [8], for example, refrigerants with HFCs, perfluorocarbons, or hydrofluorocarbons. The dairy sector is fundamental and strategic for the country and represents 1.4% of Ecuador's agri-food [9]. Refrigerants are needed to transport dairy products, for example, those derived from milk. The province of Cotopaxi is one of the referents of the dairy sector of this country.
Hydrofluoroolefins (HFOs), hydrochlorofluoroolefins, hydrocarbons, carbon dioxide, ammonia, siloxanes, and hydrofluoroethers are being investigated [3]. HFOs, called the fourth-generation refrigerants, are composed of hydrogen, fluorine, and carbon atoms and contain at least one double bond in their structure. These kinds of compounds are eco-friendly and present cost-effectiveness and high energy efficiency [10].
Various authors have dealt with identifying structural formulas of working fluids and the synthesis of new compounds, which is carried out by applying computer-aided molecular design (CAMD) [11]. In general, a CAMD method finds molecular structures with desirable properties using group contribution methods [12].
Joback and Reid [13] implemented a set of group contribution estimation methods using molecular groups that make up a wide variety of organic compounds. Molecular design is performed by developing a mixed integer nonlinear programming (MINLP) model, whose solution provides a molecule with desirable refrigerant properties. It uses an expression involving the enthalpy of vaporization and the liquid heat capacity as the objective function [14][15][16]. By employing this methodology, 53 refrigerant molecules [15,16] were added to the literature, accounting for 328 molecules [17,18].
New refrigerant molecules can be found by increasing the regions of search space or changing the objective function. Therefore, this article presents non-reported refrigerant molecules by solving an optimization problem using an MINLP methodology and the generalized reduced gradient (GRG) nonlinear method from the Microsoft Excel solver tool.

Functional Groups and Structural Constraints
Functional groups selected for this study are expected in existing refrigerants, whose properties were reported in the literature [1,13]. The total number of considered groups was equal to 16, which are shown in Table 1. It is worth noting that alkene groups were included in the search to find HFOs.  (1) represents the octet rule, which provides a simple relationship to obtain molecules without free attachments. For achieving new molecules, each groups' availability must be greater or equal to free attachments from each group type, as indicated in Equations (2) and (3) [15].

Halogen Groups
where n j and v j are the number and valence, respectively, of groups of type j. Refrigerant molecules of low molecular weight will be searched. Therefore, an additional constraint is assigned to the number of groups of the same type to keep it within a range limited by a lower (n j l ) and an upper (n j u ) limit (in Equation (4)), which were equal to 2 and 7, respectively. The total number of groups and those containing halogen, oxygen, and nitrogen were maintained between 2 to 14 (in Equation (5)) and 0 to 10 (in Equation (6)), respectively, as reported by Duvedi and Achenie [15].
n Alcohols , and n Amines ) ≤ 10 (6) where n R indicates the number of halogen-containing groups in the molecule. Due to the strong tendency to explode, low boiling molecules with halogen and nitrogen atoms will not be considered [15,19], which can be identified by Equation (7).
Moreover, the search for double bond molecules is enabled by Equations (8) and (9).

Property Targets
Additional restrictions are related to properties of Freon 12 such as the evaporation (T e ) and condensing (T cd ) temperature (−1.1 • C and 43.3 • C, respectively), the enthalpy of vaporization (∆H v ) at the lowest temperature, and liquid heat capacity (c pl ) at the average temperature. ∆H v and c pl will be greater or equal to 18.4 kJ mol −1 and lesser or equal to 32.2 cal mol −1 K −1 , respectively [1,20]. By using these values, the amount of the refrigerant will be decreased. Moreover, the vapor pressure of a refrigerant at the evaporation temperature P s {T e } must be higher than the atmospheric pressure, which reduces the possibility of air and humidity infiltrating the refrigeration system [20]. The value assigned to search for a refrigerant molecule is equal to 1.4 bar. However, by decreasing the value of P s {T e }, the feasible region for the search of new refrigerant is expanded, increasing the possibility of finding out new molecules similar to the corresponding for Freon 12 [15]. Therefore, to search refrigerant molecules, a P s {T e } value equal to 1 bar was considered. Moreover, the value of the refrigerant's vapor pressure at the condensation temperature P s {T cd } must not exceed the value of 14 bar due to the equipment's cost [16].

Modeling of Physical Properties
Estimation methods reported in the literature [1,13] were used to calculate both the normal boiling point (T b ) and the critical temperature (T c ) in K and the critical pressure (P c ) in bar (Equations (10)-(12), respectively). Previously, to estimate the vapor pressure, the reduced pressure was calculated (Equations (13)-(16)) [1,15,16].
where n A is the total number of atoms in the molecule, and the r subscript means the corresponding reduced property. For calculating the liquid heat capacity (c pl ) in J mol −1 K −1 , the Chueh and Swanson method (Equation (17)) was used [21].
For estimating the latent heat of vaporization (∆H v b ) at the normal boiling point, Vetere's modification of the Kistiakowsky equation (Equations (18) and (19)) was employed. For calculating the heat of vaporization (in J mol −1 at any temperature T, the Watson relation (Equations (20) and (21)) was used [15].
S vb is the entropy of vaporization at the normal boiling point, and M is the molecular weight.
The performance was calculated by Equation (22), which represents the refrigeration cycle's estimated efficiency [16].

MINLP Solution
To solve the MINLP problem, the generalized reduced gradient nonlinear (GRG nonlinear) method from the solver tool of Microsoft ® Excel 2016, which uses the branch and bound (BB) algorithm, was employed. Its configuration was the default settings of the program. For each group, an initial guess equal to 1 was used. If a variable has a non-integer solution, the BB method chooses it and assigns integer values according to the corresponding restriction, thus formulating subproblems and solving them. Therefore, subproblems are created from the problem with non-integer solutions, whose relation is similar to the tree branches. This relation is shown in Figure 1, whose intersection points of the BB tree are MINLP problems solved by the GRG code [22]. If a feasible solution for each subproblem is found out, called a candidate incumbent, this is stored and changed by a better optimal solution during the subproblems resolution. Moreover, the branches will not be generated from this subproblem. If there is no feasible solution or the algorithm tolerance is not achieved, or all of the subproblems were solved, the resolution is stopped [22,24]. The integer values for the variables corresponding to the updated candidate incumbent will be the solution.
It is worth mentioning that the MINLP problem is non-convex [15,16]. The BB algorithm implemented in the GRG nonlinear method from the solver tool of Microsoft ® Excel 2016 can find an optimal solution even when the MINLP does not satisfy the convexity conditions [22]. Therefore, the molecules reported in this work represent local optimal solutions.

Environmental and Safety Evaluations
The safety and the environmental effects were investigated in terms of the ozone depletion potential (ODP), the global warming potential (GWP), the flammability, and the toxicity.
The ODP and the GWP were evaluated using the number of chlorine and fluorine atoms, respectively, that the refrigerant molecule contains [25].
The flammability (F) applied to organic compounds was determined using the F index [26], which can be calculated using Equation (23 where p1 to p17 are constants, C1 is equal to one for the case of monocarbonic compounds. Otherwise, it is equal to zero. ROE, RCO, RCOO, RNH correspond to the number of ether, carbonyl, ester, and imine groups, RRNG and RARM are equal to aliphatic and aromatic rings, respectively, and RUS represents the number of double bonds in the carbon skeleton, including aliphatic and aromatic rings. RF, RCl, RBr, ROH, RNO2, RNH2, RCN, and RCOOH are equal to F, Cl, and Br atoms, and OH, NO2, NH2, CN, and COOH groups, respectively. The R factors from p3 to p9 are divided by the total number of skeletal carbons, while the corresponding factors from p10 to p17 are divided by the number of hydrogen atoms in the respective pure hydrocarbon molecule. The toxicity (X) was determined by the group contribution methods [27,28], as indicated in Equation (24). If a feasible solution for each subproblem is found out, called a candidate incumbent, this is stored and changed by a better optimal solution during the subproblems resolution. Moreover, the branches will not be generated from this subproblem. If there is no feasible solution or the algorithm tolerance is not achieved, or all of the subproblems were solved, the resolution is stopped [22,24]. The integer values for the variables corresponding to the updated candidate incumbent will be the solution.
It is worth mentioning that the MINLP problem is non-convex [15,16]. The BB algorithm implemented in the GRG nonlinear method from the solver tool of Microsoft ® Excel 2016 can find an optimal solution even when the MINLP does not satisfy the convexity conditions [22]. Therefore, the molecules reported in this work represent local optimal solutions.

Environmental and Safety Evaluations
The safety and the environmental effects were investigated in terms of the ozone depletion potential (ODP), the global warming potential (GWP), the flammability, and the toxicity.
The ODP and the GWP were evaluated using the number of chlorine and fluorine atoms, respectively, that the refrigerant molecule contains [25].
The flammability (F) applied to organic compounds was determined using the F index [26], which can be calculated using Equation (23).
where p 1 to p 17 are constants, C 1 is equal to one for the case of monocarbonic compounds. Otherwise, it is equal to zero. R OE , R CO , R COO , R NH correspond to the number of ether, carbonyl, ester, and imine groups, R RNG and R ARM are equal to aliphatic and aromatic rings, respectively, and R US represents the number of double bonds in the carbon skeleton, including aliphatic and aromatic rings. R F , R Cl , R Br , R OH , R NO2 , R NH2 , R CN , and R COOH are equal to F, Cl, and Br atoms, and OH, NO 2 , NH 2 , CN, and COOH groups, respectively. The R factors from p 3 to p 9 are divided by the total number of skeletal carbons, while the corresponding factors from p 10 to p 17 are divided by the number of hydrogen atoms in the respective pure hydrocarbon molecule.
The toxicity (X) was determined by the group contribution methods [27,28], as indicated in Equation (24).
where n j and α j are the number and the toxicity, respectively, of groups of type j. Table 2 lists the optimized objective functions together with the molecular structures obtained by solving the MINLP formulation. The size of the problem was 18 continuous variables. Thirteen continuous variables are the boiling temperature (T b ), critical temperature (T c ), critical pressure (P c ), reduced boiling temperature (T br ), reduced evaporating temperature (T er ), reduced condensing temperature (T cdr ), reduced vapor pressure at the evaporating temperature (P sr {T e }), reduced vapor pressure at the condensing temperature (P sr {T cd }), the vapor pressure at the evaporating temperature (P s {T e }), the vapor pressure at the condensing temperature (P s {T cd }), liquid heat capacity (c pl ), the entropy of vaporization at the normal boiling point (S vb ), and enthalpy of vaporization (∆H v ). Five additional continuous variables are G, h, k, ∆H v b , and n, defined by Equations (14)- (16), (18) and (21), respectively.  Table 2), the alkene groups shown in Table 1 were not included in the search. Moreover, the SHF molecule groups were considered new guesses to find out the CH 3 FO 2 molecule. The use of the number of groups from undesired molecules as a new guess for refrigerant searching was reported by Duvedi and Achenie [15]. The SHF refrigerant has been reported in the literature [16].

Results and Discussion
It is worth mentioning that the C 2 H 4 F 2 O compound (row 2, Table 2) can be assumed to be an isomer of the R-E152a refrigerant [29], and the C 3 H 4 F 2 [30] and CH 2 F 2 O [18] formulas were used in the literature to report compounds with a molecular structure different to the corresponding reported in this work.
Related to the molecules' environmental characteristics found in this work, all of the molecules have null ODP due to the absence of chlorine atoms [25]. Since the GWP values depend on the number of fluorine atoms in the molecular structure [25], all of the molecules have a GWP value similar to that from the R-134a and R-12, which have four and two fluorine atoms, respectively.
The F values range from 0.4 to 0.6 for all new refrigerant molecules, which is considered normally flammable [26]. It is worth mentioning that the F value for the case of R-134a (C 2 H 2 F 4 , Table 3) and R-12 (CCl 2 F 2 , Table 3) was equal to 0.4478 and 0.1517, respectively. By minimizing the cosine of the ratio between the c pl and ∆H v and using a value equal to 1.4 bar for the P s {T e }, the F 2 O 2 molecule can be find out. However, this molecule decomposes above 195 K to produce O 2 and F 2 [31]. Fluorine gas is corrosive to most common materials and is toxic [32]. Therefore, the F 2 O 2 molecule was not considered as a new refrigerant. Moreover, the toxicity evaluation is determinant for searching molecules with physical properties similar to R-12 molecule.

Conclusions
New molecules instead of Freon 12 can be found using new objective functions, which account for Freon 12 properties. Using the solver tool of Microsoft ® Excel 2016 through the GRG nonlinear method, which uses the branch and bound algorithm, the formulated objective functions were optimized. The CH 3 FO 2 , C 2 H 4 F 2 O, CH 2 F 2 O 2 , CH 2 F 2 O, C 3 H 4 F 2 , and C 2 H 2 F 2 O molecules meet structural feasibility and physical constraints found for Freon 12. These compounds show global warming potential values similar to those corresponding to R-134a and Freon 12 refrigerants and present null ozone depletion potential values. Moreover, these molecules are normally flammable. The CH 3 FO 2 , C 2 H 4 F 2 O, CH 2 F 2 O 2 , C 2 H 2 F 2 O, and CH 2 F 2 O molecules show toxicity values similar to R-134a and Freon 12 refrigerants.