A Review on Process Modeling and Simulation of Cryogenic Carbon Capture for Post-Combustion Treatment

The cryogenic carbon capture (CCC) process is a promising post-combustion CO2 removal method. This method is very novel compared with conventional and well-developed methods. However, cryogenic carbon capture is not yet commercially available despite its techno-economic benefits. Thus, a model-based design approach for this process can provide valuable information. This paper will first introduce the cryogenic carbon capture process. Then, a comprehensive literature overview that focuses on different methods for modeling the process at the component level will be given. The modelling methods which are deemed most effective are presented more in depth for each of the key system components. These methods are compared with each other in terms of complexity and accuracy and the simplest methods with an acceptable level of precision for modelling a specific component in the CCC process are recommended. Furthermore, potential research areas in modeling and simulation of the CCC process are also highlighted.


Introduction
Due to population expansion, the world's growing economy, and increasing urbanization, the global energy demand has been increasing exponentially [1]. Although renewables are expanding at high rates, they still cannot satisfy the global energy demand, and because of that, 80% of the current energy demand is still supplied by fossil fuels [2]. The large-scale use of fossil fuels has resulted in air pollution and has deteriorated human health. Based on a recent analysis, global CO 2 emissions increased from 33.3 Gigatons/year in 2020 to 34.9 Gigatons/year in 2021 [3]. This prompt increase in global CO 2 emissions has brought about numerous intense problems including climate change, global warming, rising sea level, etc. Using renewable energy systems, increasing the efficiency of energy conversion systems, and utilizing CO 2 capture technologies are the most practical solutions to regulate CO 2 emissions [4]. The latter method is very advantageous, especially in cases such as cement production where the emissions are not avoidable. This method directly targets the emissions at the source and may also remove other pollutants such as NO x . Moreover, the integration of CO 2 capture technologies within large industrial plants is convenient and, in some cases, no fundamental change in the plant structure is required. For instance, the cryogenic carbon capture process is a bolt-on retrofit technology since the process just uses electricity, while amine-based processes just need heat and process integration. This technology can also provide an opportunity to supply and store CO 2 with high purities which can be used in fuel production or oil recovery processes. Pre-combustion, oxyfuel combustion, and post-combustion are three available CO 2 capture methods.
In the case of using the pre-combustion capture method, CO 2 is removed from a gas mixture prior to combustion [5]. In other words, hydrocarbon fuel is converted into carbon monoxide and syngas. Then, the generated carbon monoxide reacts with steam to produce hydrogen and CO 2 . Finally, the produced CO 2 is separated from hydrogen that can be combusted without generating any pollutant. This technology is beneficial since it is verified for industrial-scale oil refineries, and it can recover 90-95% of CO 2 . However, it requires high investment costs. Moreover, high NO x emissions need expensive scrubbers.
In the oxyfuel process, as the name indicates, pure oxygen is used rather than air in the combustion chamber, which results in generating flue gases that consist of almost pure CO 2 and vapor. As a result, the separation, and purification of CO 2 in this method is very simple [6]. Therefore, CO 2 is captured and stored just by the condensation of vapor content in the flue gas and low-temperature purification. When this method is applied, a plant should be equipped with an air separation unit (for oxygen production), a flue gas processing unit as well as a CO 2 processing unit. This method can be considered rewarding since it results in very low emissions of NO x , and it has high flexibility for using different types of fuels. Furthermore, achieving 100% combustion of biomass, better behavior of the combustion process, and using a smaller furnace are other attractive features of oxyfuel. However, replacing O 2 with air may bring about some challenges. In other words, the energy requirement for oxygen production is high, and keeping the flame stable is challenging [6,7]. Therefore, the net power production from the plant will be lower and flue gas should be recycled in large quantities to keep the temperature at reasonable levels.
Post-combustion CO 2 capture methods include chemical absorption, physisorption, membrane-based capture processes, and cryogenic carbon capture (CCC) processes. In the CCC process, gaseous CO 2 can be separated from other components in the flue gas due to their different condensation and desublimation temperatures. This method is gaining more popularity by using less energy [8][9][10][11][12][13][14][15], capturing CO 2 with higher rates and purities, and having faster responses to fluctuations in electricity demand [8,16,17]. Moreover, investigations have indicated that CO 2 capture costs can be decreased by 20-40% when the CCC process is used, and it distinctly indicates the economic advantages of the CCC process over other methods [18]. Furthermore, this technology can be easily added to the current industrial emission facilities without any concern regarding chemical solvents or physical sorbents [19]. This method can use natural gas as a refrigerant and because of that, it can store energy in the form of liquid natural gas (LNG). For capturing CO 2 using this method, only 12-18% of the power generated in the plant is consumed if the energy storage operates with efficiencies higher than 90% [8]. The energy supplied to the process is mostly used to drive compressors in refrigeration cycles mainly to liquefy natural gas and provide the required LNG for the process. The Cascade liquefaction cycle, mixed refrigerant cycle, and gas expander cycle are the available scenarios to liquefy natural gas and they are discussed in reference [20]. Despite all benefits and advantages of the CCC process, this technology is still in its early stages of development and deals with several challenges such as ice plugging when the flue gas is not water-free [15]. Moreover, this technology can be considered uneconomical when the CO 2 mole fraction in the gas mixture is very low. In this case, a large amount of gas other than CO 2 should be cooled to condense a very low amount of CO 2 .
It should be noted that capture conditions have a great impact on the total costs and the required energy used for the separation and compression of CO 2 when the CCC process is used [21]. These conditions are the CO 2 concentration in flue gases, pressure, temperature, the chemical composition of flue gases, and the required purity of captured CO 2 . However, among all aforementioned factors, the concentration of CO 2 in the flue gas affects the total costs and required energy more effectively.
Techno-economic evaluations have indicated that for large-scale applications, the CCC process costs are comparable with mature conventional monoethanolamine (MEA) absorption technology, while in smaller scales, the costs are significantly lower in the case of using the CCC process rather than MEA [22]. This process can generate extra LNG and store it when renewable electricity is available, or the plant power demand is low. The cooling energy stored in the form of LNG is then recovered during the peak hours when electricity available models provided in this study can be effectively used to simulate, evaluate, and optimize the process. For each component in the process, the modeling approach along with the accuracy, advantages, and disadvantages are highlighted.
In this study, first the flue gas pre-processing is discussed. Then, the operation of the CCC process integrated with the energy storage is explained and different components and elements in the process are introduced. Afterwards, the available models and equation of states which can be used to model and investigate these elements including heat exchanger, bubbler, LNG storage tank, and compressor are discussed thoroughly. In particular, the differential method is discussed, which can be used to model the heat exchangers in the CCC process while the particle velocity model and pseudo-homogeneous one-dimensional plug flow model are given to model the desublimation heat exchanger. Different equation of states and models such as PR EoS, the Rachford-Rice equation, and SRK EoS are discussed which can be used to calculate the solid-vapor equilibrium and liquid-vapor equilibrium of CO 2 mixtures in the bubbler. The minimum energy model is also presented here which can be applied to calculate the minimum energy required to separate CO 2 from the flue gas in the bubbler. Furthermore, the RGibbS reactor block which is available in ASPEN Plus process simulator software is also introduced and its accuracy to calculate solid-vapor equilibrium is evaluated. Moreover, the volume of fluid and BWRS methods which can be used to model the LNG tank as an energy storage in the CCC process and calculate the boil-off-gas are also discussed. The compression process model of scroll compressors and the linear compressor model are also given to model and study the scroll compressors and linear compressors, respectively. This section is followed by the accuracy analysis section where a comprehensive comparison of the accuracy of different models in the literature is provided. Finally, this study is concluded, and the research gap is presented.

CCC Process Description
In this section, the pre-processing of the flue gas before entering the CCC process and the performance of the CCC process are concisely discussed. As Figure 1 indicates, the flue gas which leaves the plant should be cooled down to 40 • C (state 1) using a vapor compression refrigeration system or other methods (this temperature is considered as flue gas reference temperature for CO 2 capture processes according to NETL). Afterward, the flue gas temperature is further reduced by water cooling which results in reducing the flue gas temperature to 25 • C (state 2). The pressure of the flue gas is then increased to 1.4 bar using a blower (state 3). The flue gas is then cooled to 2 • C (state 4) and its hazardous components including mercury and sulfur dioxide will be washed away (state 11) as it comes to direct contact with chilled water at 1 • C (state 10) and it is directed to a dryer to remove its vapor content. Herein, the pre-processed flue gas (state 5) enters the CCC process to separate its CO 2 content. As a result, the separated CO 2 is stored in the form of liquid (state 6), while the cold nitrogen-rich gas mixture (state 7) might be used to generate chilled water before being released into the atmosphere (state 8). The CCC process is shown in Figure 2. The process flow diagram shown by Figure 2 is similar to the process studied and evaluated by Jensen [8] and Fazlollahi et al. [23]. The dry flue gas after desulfurization enters the CCC process and it is cooled to 175 K by a nitrogen-rich gas mixture which comes from the desublimation heat exchanger and liquid The CCC process is shown in Figure 2. The process flow diagram shown by Figure 2 is similar to the process studied and evaluated by Jensen [8] and Fazlollahi et al. [23]. The dry flue gas after desulfurization enters the CCC process and it is cooled to 175 K by a nitrogen-rich gas mixture which comes from the desublimation heat exchanger and liquid CO 2 stream. Then, the flue gas enters a desublimation heat exchanger system, and its temperature is reduced to 154 K by contact liquid. As a result, the CO 2 in the flue gas undergoes a phase change and solid particles of CO 2 are formed. The clean flue gas leaves the desublimation heat exchanger system from the top and cools the incoming flue gas, while the CO 2 /contact liquid slurry leaves the heat exchanger from the bottom and it is pressurized by a pump. The solid CO 2 particles are separated from the contact liquid by filtration. As Figure 2 indicates, the process requires two refrigeration loops termed the internal and external refrigeration cycles which have the highest energy consumption in the process. The internal refrigeration cycle which may use mixed refrigerant is responsible for providing a part of cooling energy that is used to cool down the contact liquid. As can be seen, the compressed mixed refrigerant is condensed by melting the solid CO 2 and increasing its temperature to 233 K. The liquid CO 2 stream flashes to generate pure liquid CO 2 . The liquid-rich CO 2 stream cools the incoming flue gas before leaving the CCC process. The external refrigeration loop is also used to produce liquified natural gas to operate the process and provide the rest of the cooling effect for reducing the temperature of contact liquid. The C3MR natural gas liquefaction process is discussed in detail by Wang et al. [35]. The liquefied natural gas is used as an intermediate refrigerant to cool the contact liquid which is responsible for the desublimation of CO 2 in a heat exchanger. As shown in Figure 2, the LNG leaves the natural gas liquefaction process at a temperature of 179 K. When renewable energy is available and the energy price is low, extra LNG can be generated and stored in a tank. The stored LNG can be used during peak hours when the energy price is high. The LNG which comes from the tank/natural gas liquefaction process is expanded, and its temperature is reduced to 153.6 K; then, it enters a multi-stream heat exchanger to cool the contact liquid. The contact liquid is generally a hydrocarbon mixture, but it can be any other heat transfer fluid that has low viscosity and vapor pressure at cryogenic temperatures [8]. As shown in Figure 2, this liquid is used to separate CO 2 in a direct heat exchanger where it is mixed with the flue gases. Using a direct heat exchanger is operationally simple and can desublimate CO 2 continuously, while as Figure 2 indicates, a solid separation unit is required to separate solid particles of CO 2 from the liquid. In the case of using an indirect heat exchanger, the solid separation unit is not required. However, using this type of heat exchanger adds complexity to the process by causing challenges in cleaning the surface of heat exchangers, requiring maintenance cycles, and avoiding a continuous operation. Energies 2023, 16, x FOR PEER REVIEW 6 of 35

Modeling Approaches for the CCC Process
In this part, the most important models which are available in the literature and can be used to simulate the CCC process are discussed, and their advantage and accuracy are also evaluated. Figure 3 indicates different components and elements in the CCC process. The available methods which can be used to model and investigate each component are shown by this figure. In the figure, only the models that are deemed more effective for components modelling in the CCC processing are presented. In this section, a differential method for modeling the heat exchangers in the CCC process is discussed. This model can be used to simulate counter flow heat exchangers as well as multi-stream heat exchangers in the CCC process. Then, the particle velocity model and pseudo-homogeneous one-dimensional plug flow model, which can be used to study and model the desublimation heat exchanger in the CCC process, are presented. Furthermore, different equations which can be used to calculate the solid-vapor equilibrium and liquid-vapor equilibrium of CO2 mixtures in the bubbler are presented in this section. These equations include PR EoS, the Rachford-Rice equation, and SRK EoS. The minimum energy model is also discussed in this section which can be used to calculate the minimum energy required for the separation of CO2 from the flue gas using the CCC process under the ideal operating conditions. The volume of fluid and BWRS methods are also presented to model the LNG tank and calculate the boil-off-gas. Finally, the compression process model of scroll compressors and the linear compressor model are discussed to model the scroll compressors and the linear compressors, respectively.

Modeling Approaches for the CCC Process
In this part, the most important models which are available in the literature and can be used to simulate the CCC process are discussed, and their advantage and accuracy are also evaluated. Figure 3 indicates different components and elements in the CCC process. The available methods which can be used to model and investigate each component are shown by this figure. In the figure, only the models that are deemed more effective for components modelling in the CCC processing are presented. In this section, a differential method for modeling the heat exchangers in the CCC process is discussed. This model can be used to simulate counter flow heat exchangers as well as multi-stream heat exchangers in the CCC process. Then, the particle velocity model and pseudo-homogeneous onedimensional plug flow model, which can be used to study and model the desublimation heat exchanger in the CCC process, are presented. Furthermore, different equations which can be used to calculate the solid-vapor equilibrium and liquid-vapor equilibrium of CO 2 mixtures in the bubbler are presented in this section. These equations include PR EoS, the Rachford-Rice equation, and SRK EoS. The minimum energy model is also discussed in this section which can be used to calculate the minimum energy required for the separation of CO 2 from the flue gas using the CCC process under the ideal operating conditions. The volume of fluid and BWRS methods are also presented to model the LNG tank and calculate the boil-off-gas. Finally, the compression process model of scroll compressors and the linear compressor model are discussed to model the scroll compressors and the linear compressors, respectively.

Heat Exchangers
As Figure 2 indicates, the CCC process deals with heat exchangers of different types. In this process and other cryogenic processes in general, plate-fin heat exchangers are widely used due to their exceptional advantages including compactness, low weight, and high efficiency. It is also worth noting that plate-fin heat exchangers are the only type of heat exchanger that can carry multiple streams [36]. Therefore, using an efficient, accurate, and fast method to model plate-fin heat exchangers are of paramount importance for modeling and investigating the CCC process. Different methods are available and can be used to simulate the heat exchangers. For instance, You et al. [37] developed a CFD method to model the shell and tube heat exchanger. They used porosity and permeability concepts in their model to reduce the computational costs. They modeled the heat exchanger over a wide range of Reynolds numbers, and they indicated that the model can predict the velocity and temperature fields with a maximum of 15% relative deviation. He et al. [38] also developed a modified porous medium model to evaluate flow characteristics and heat transfer on the shell side of shell and tube heat exchangers. They simulated three types of shell and tube heat exchangers with vertical baffles, helical baffles, and finned tube banks. Comparing the model results with experimental data showed that the model can predict a pressure drop with a maximum of 25.1% relative deviation. Haider et al. [39] also used Open-FOAM v6 to simulate a multi-stream plate-fin heat exchanger. In this study, they directly modeled and solved parting sheets and sidebars. However, to reduce computational costs and avoid the fine computational grid, the fins were not resolved directly. Therefore, they also used a porous medium approach and considered fins as porous solids. Korzeń et al. [40] also developed a mathematical model which can simulate the transient operation of a plate and fin heat exchanger to cool water with cold air. It can

Heat Exchangers
As Figure 2 indicates, the CCC process deals with heat exchangers of different types. In this process and other cryogenic processes in general, plate-fin heat exchangers are widely used due to their exceptional advantages including compactness, low weight, and high efficiency. It is also worth noting that plate-fin heat exchangers are the only type of heat exchanger that can carry multiple streams [36]. Therefore, using an efficient, accurate, and fast method to model plate-fin heat exchangers are of paramount importance for modeling and investigating the CCC process. Different methods are available and can be used to simulate the heat exchangers. For instance, You et al. [37] developed a CFD method to model the shell and tube heat exchanger. They used porosity and permeability concepts in their model to reduce the computational costs. They modeled the heat exchanger over a wide range of Reynolds numbers, and they indicated that the model can predict the velocity and temperature fields with a maximum of 15% relative deviation. He et al. [38] also developed a modified porous medium model to evaluate flow characteristics and heat transfer on the shell side of shell and tube heat exchangers. They simulated three types of shell and tube heat exchangers with vertical baffles, helical baffles, and finned tube banks. Comparing the model results with experimental data showed that the model can predict a pressure drop with a maximum of 25.1% relative deviation. Haider et al. [39] also used Open-FOAM v6 to simulate a multi-stream plate-fin heat exchanger. In this study, they directly modeled and solved parting sheets and sidebars. However, to reduce computational costs and avoid the fine computational grid, the fins were not resolved directly. Therefore, they also used a porous medium approach and considered fins as porous solids. Korzeń et al. [40] also developed a mathematical model which can simulate the transient operation of a plate and fin heat exchanger to cool water with cold air. It can be used to model the transient response of the heat exchanger when the cold stream mass flow rate decreases suddenly and at the same time, the hot stream mass flow rate increases. Comparing the numerical results with experimental data showed that the model can predict the temperature of the outlet streams accurately.
In this part, a differential method for modeling heat exchangers in the CCC process is described and evaluated. This model can be used to simulate counter flow heat exchangers with acceptable accuracy and can also be extended to simulate multi-stream heat exchangers. Therefore, this method is advantageous, and because of that, it is discussed here.

Differential Method
The differential method is an accurate technique for modeling heat exchangers since, unlike the integral method which considers fixed values for thermo-hydraulic properties of hot and cold streams, this method considers the variation in thermo-hydraulic properties of process streams and integrates the heat transfer and pressure loss functions in a stepwise approach from one end to another end of heat exchangers. It has been shown that this method is very effective and efficient when the structure of the studied heat exchanger is very complex [41].
In the differential method, as shown in Figure 4, the length of the heat exchanger is divided into N sections where the thermo-hydraulic properties of streams are assumed to be invariant in each section. The temperature variation in the first section of the counter-flow heat exchanger for hot and cold streams can be calculated using the following equations [36]: where (UA) is the overall heat transfer rate, while C c and C h represent the capacity rate of cold and hot streams, respectively, and ∆T 1 is the temperature difference between inlet hot and cold streams, which is always known. This value can be used to calculate the temperature difference at the starting point of the second partition, as follows: Energies 2023, 16, x FOR PEER REVIEW 8 of 35 be used to model the transient response of the heat exchanger when the cold stream mass flow rate decreases suddenly and at the same time, the hot stream mass flow rate increases. Comparing the numerical results with experimental data showed that the model can predict the temperature of the outlet streams accurately. In this part, a differential method for modeling heat exchangers in the CCC process is described and evaluated. This model can be used to simulate counter flow heat exchangers with acceptable accuracy and can also be extended to simulate multi-stream heat exchangers. Therefore, this method is advantageous, and because of that, it is discussed here.

Differential Method
The differential method is an accurate technique for modeling heat exchangers since, unlike the integral method which considers fixed values for thermo-hydraulic properties of hot and cold streams, this method considers the variation in thermo-hydraulic properties of process streams and integrates the heat transfer and pressure loss functions in a stepwise approach from one end to another end of heat exchangers. It has been shown that this method is very effective and efficient when the structure of the studied heat exchanger is very complex [41].
In the differential method, as shown in Figure 4, the length of the heat exchanger is divided into N sections where the thermo-hydraulic properties of streams are assumed to be invariant in each section. The temperature variation in the first section of the counterflow heat exchanger for hot and cold streams can be calculated using the following equations [36]: where (UA) is the overall heat transfer rate, while and represent the capacity rate of cold and hot streams, respectively, and Δ is the temperature difference between inlet hot and cold streams, which is always known. This value can be used to calculate the temperature difference at the starting point of the second partition, as follows: (3) Figure 4. Discretization of a counter-flow heat exchanger into N sections for temperature calculation using the differential method.
Consequently, the temperature variation in the second section of the heat exchanger for hot and cold streams can be calculated using the following equations: Therefore, temperature variations for hot and cold streams in the third partition of the heat exchanger can be found by multiplying Equations (4) and (5) by the second term on the right side of Equation (3). Under this situation, temperature variation in the last partition for hot and cold streams can be calculated using the following equations: The temperature difference between hot and cold streams at the cold end of the heat exchanger can also be obtained using the following equation: where ∆T ce is the temperature difference between hot and cold streams at the cold end of the counter flow heat exchanger which is the right wall of section N in Figure 4. At this point, the temperature distribution of cold and hot streams along a counter-flow heat exchanger can be obtained using Equations (1)-(8), having the inlet temperature of hot and cold streams. Considering Equations (1)-(8), temperature variation sets along the heat exchanger forms a geometric series. Therefore, summation of the series and letting N tend to infinity generates the following equations: In the case of using multi-stream heat exchangers, Equations (9)-(11) can be extended in the following forms to calculate the temperature difference between hot and cold streams in both hot and cold ends: where nc and nh denote the number of cold and hot streams, respectively, while h i and A i are the heat transfer coefficient and net effective surface area of the stream i, respectively. The surface temperature of a multi-stream heat exchanger can also be calculated using the following equation: As Equation (14) indicates, it is assumed that the surface temperature of the heat exchanger in all cross sections has a constant value. This assumption brings about significant simplicity in calculations by limiting the unknown parameters to stream temperatures and surface temperature. It is also worth noting that when a satisfactory stacking pattern is used, especially in multi-stream heat exchangers, the model which uses this assumption can predict temperature with very low errors [39,40].
The differential method is discussed in more detail in [36,42] by providing and designing a computer program for evaluating the performance of a counter flow heat exchanger.
Wang et al. [38] developed a novel segmented differential method to model and optimize a multi-stream compact heat exchanger. In this method, a multi-stream compact heat exchanger is fragmented into several sections along the direction of the mainstream. Each section is considered as a sub-heat exchanger which is made of N layers and N − 1 separating plates along the height of the heat exchanger. This model is developed to determine the temperature and pressure of streams along with the temperature of separating plates at each computational node. Therefore, heat transfer equations along with energy conservation equations in the entire computational domain of a sub-heat exchanger should be solved to calculate the unknown parameters. The calculated parameters are then used as inputs for the next sub-heat exchanger. The energy conservation for streams is presented by the following equation [41]: On the left side of Equation (15), the energy gradient at i th node is indicated where D represents the sign symbol of flow direction while T and F are the stream temperature and heat capacity flow rate, respectively. On the right side of this equation, however, convective heat transfer on the separating plates and fins is provided where ω, δ, b, and h stand for fin number per unit width, fin thickness, fin height, and convective heat transfer coefficient, respectively. Subscriptions p, f , and i also indicate separating plate, fin, and the number of nodes, respectively.
Energy conservation for separating plates can also be presented by the following equation: where λ represents thermal conductivity. Finally, energy conservation for fins is given by making a balance between the convective heat transfer rate (the left side of equation) and conduction heat transfer rate (the right side of the equation), as follows: The differential method is an accurate model which can precisely determine the variation in the temperature and pressure field inside the heat exchanger [42]. However, this model has a higher computational cost when it is compared with the integral method, especially when it is used for dynamic models. Therefore, the integral method would be more applicable for dynamic applications when the internal pressure and temperature fields are not required.

Bubbler (CO 2 Separator)
In this part, the available models and equations of states which can be used to model the desublimation heat exchangers and packed beds in the CCC process and calculate the solid-vapor equilibrium and liquid-vapor equilibrium are given. In the analyzed studies, the models and equations of states which have mainly been considered for their accuracy and computational effectiveness can be listed as:

•
Particle velocity model; The last model has been used to calculate the solid-vapor equilibrium temperature by Guido et al. [43]. However, in this study, it is extended to calculate the liquid-vapor equilibrium. The model's accuracy compared to the experimental value will be later provided in Table A1.
In this part, the particle velocity model and pseudo homogeneous one-dimensional plug flow model are discussed first, and reference to the relevant studies is provided. These two models can be used to study, investigate, and simulate the desublimation heat exchanger system and packed beds in the CCC process. Then, the equations used to perform the phase equilibrium analysis of CO 2 mixtures at cryogenic temperatures are given.

Particle Velocity Model
The desublimation process in the bubbler shown in Figure 2 is simulated by James et al. [44] using the particle velocity method. This model was applied to simulate the phase change processes including the desublimation of CO 2 . Therefore, this model solves heat and mass transfer equations simultaneously along with the momentum transfer equation to track the phase changes.
A one-dimensional numerical model of the desublimation heat exchanger system is used to predict temperature distribution, velocity, pressure, composition, mass, and heat transfer rates inside the heat exchanger. It was assumed that the local solid-vapor equilibrium exists in the particle interface and these particles are distributed uniformly inside the heat exchanger and are displaced equally. The CO 2 -rich flue gas stream enters the bottom of the heat exchanger. In this case, as discussed before, a cool exchange medium enters the top of the heat exchanger and cools the flue gas from −98 • C to very low temperatures (−138 • C) [44]. Therefore, CO 2 is separated from the flue gas in the form of solid particles. As a result, the CO 2 /exchange liquid slurry leaves the bottom of the heat exchanger, while CO 2 lean gas leaves the top of the heat exchanger.
In this model, the following equation is used to describe convection heat transfer from or to droplets.
where A, N, h, and θ represent the droplet surface area, the number of solid particles, convective heat transfer coefficient, and blowing factor, respectively. It is assumed that droplets have a constant surface temperature and composition. The convective heat transfer coefficient is given as follows: where k, d p , and Nu represent the conductive heat transfer coefficient, particle diameter, and Nusselt number, respectively. The latter parameter can be calculated using the following equation: The rate of evaporation and desublimation can also be calculated using the following equation: where x A∞ and W B0,i represent the mole fraction of CO 2 in upstream flow and the rate of evaporation, respectively. P vap and k xm are also vapor pressure and mass transfer coefficient, respectively, which are given as follows [44,45]: where D AB is binary diffusivity and Sh is the Sherwood number, which can be calculated as follows: where Sc denotes the Schmidt number, which is defined using the following equation: The particle velocity model represented here can be effectively used to study and evaluate the motion of solid CO 2 particles inside the heat exchanger and to determine their velocity. Therefore, the model considers the main forces including buoyancy and drag which affect the particles and can be calculated as follows: where F b , F d , C D , and υ rel are buoyancy force, drag force, drag coefficient, and velocity of particle relative to the cool exchange medium, respectively, while g and p subscriptions represent the gas phase and solid phase for particles, respectively. Moreover, g stands for gravity and possesses a negative value in Equation (26), and the direction of drag force is also in the opposite direction of particle relative velocity. Therefore, considering Equations (26) and (27), the overall force applied to particle i can be calculated as follows: Considering Equation (28) and Newton's third law, the velocity of particle i is given as follows: where m p indicates the mass of particle. It should be noted that the velocity of particles calculated using Equation (29) cannot exceed terminal velocity which is defined as the highest velocity limit which a particle can approach. Therefore, when a particle reaches this velocity, it neither accelerates nor decelerates. The value of terminal velocity can be calculated by assigning a zero value for overall force in Equation (28) and solving the equation for the velocity. This model which considers different heat and mass transfer relationships can be used to investigate, optimize, and simulate the desublimation heat exchanger system which is made of vertical countercurrent heat exchangers. This heat exchanger system can be used to cool down the rich CO 2 flue gas from −98 to −138 • C where the pressure may vary from ambient pressure to several bars. To model this heat exchanger system, momentum, species continuity, and energy equations are solved at the same time for both solid and gas phases in a one-dimensional heat exchanger which operates under steady-state conditions.

Pseudo-Homogeneous One-Dimensional Plug Flow Model
This model is applied to study the dynamic behavior of packed beds used in the CCC process [46,47]. In this case, flue gas, which is composed of N 2 , CO 2 , and H 2 O at a comparatively high temperature is fed to an initially refrigerated packed bed which results in the separation of the components since they have different dew and sublimation points. The separation process consists of three steps including the capture cycle, CO 2 recovery cycle, and cooling cycle. These three cycles are discussed in detail in reference [47]. The process for capturing CO 2 can be simulated using a Pseudo homogeneous one-dimensional plug flow. To simulate the process, component mass balances for the gaseous phase are given as follows [47]: where g , ω i,g , υ g , D e f f , a s , z, and . m i represent bed void fraction, mass fraction of i th gaseous component, superficial velocity, effective diffusion coefficient, specific solid surface area per unit bed volume, the axial coordinate of the bed, and mass deposition rate per unit surface area for component i, respectively.
The following equation also provides component mass balance for the solid phase: where m i denotes the mass deposition of component i per unit volume of the bed. The total continuity equation for the gaseous phase is given as follows: In this model, the energy balance for gas and solid phases is given as follows: . m i a s ∆h i (33) where C p,s and C p,g represent the specific heat for solid and gas phases, respectively, while ∆h i and λ e f f are enthalpy change due to phase change of component i and effective conductivity, respectively. The latter parameter can be calculated as follows: λ e f f = λ bed, 0 + RePrλ g Pe ax + Re 2 Pr 2 λ g where λ bed, 0 is the effective bed conductivity at no flow condition, while Pe ax is the Peclet number for axial heat dispersion and can be calculated as follows: where p is defined as follows: Nusselt number can be calculated to obtain gas-to-particle heat transfer coefficient as follows: Axial mass dispersion can also be calculated as follows: Mass deposition/condensation and evaporation parameters used in Equations (30)-(33) are given in the following equation: where y i,s and g represent the molar fraction of component ith in solid phase and mass deposition rate constant, respectively, while σ superscript indicates the equilibrium condition. To calculate the axial temperature and mass deposition profiles, the partial differential equations given by Equations (30)- (33) can be converted to ordinary differential equations using the method of lines (MOL). The utilization and implementation of this method are discussed in reference [48]. The resulting system of ordinary differential equations can be solved by ODE solvers such as ODE 15 s, which are available in MATLAB ® . Furthermore, the above-mentioned one-dimensional convection-dominated PDEs can also be solved using a numerical algorithm based on WENO principles by considering a higher-order discretization for the convection terms and automatic local grid adaptation. More details about the numerical algorithm can be found in Smit et al. [49].

Peng Robinson
Peng Robinson Equation of State (PR EoS) is one of the most popular and accurate applied models for thermodynamic and volumetric calculations since 1976 and can be used to obtain thermo-physical and phase equilibrium properties of CO 2 . The most convenient form of PR EoS is given by the following equation: where T, υ, and p represent temperature, specific volume, and pressure, respectively, while a and b are empirical parameters. These two parameters are used to consider the effect of attraction forces between the molecules and the molecular volume. For pure material, they can be calculated as follows [50]: where T c , P c , and ω represent critical temperature, critical pressure, and acentric factor, respectively. This model has been applied, investigated, and modified by numerous researchers and more than 220 modifications have been reported for this model which can be found in a paper published by Lopez-Echeverry et al. [51]. Although in many cases the modifications applied to PR EoS follow the same ideas, differences in details are substantial, and because of that, each version generates a definite improvement compared to other versions of the model [52]. PR EoS can be used to calculate the fugacity coefficient of components in a mixture that is involved in the equilibrium condition to predict the solidvapor or liquid-vapor equilibrium temperature. Therefore, to apply the PR EoS represented by Equations (40)- (43), it is necessary to use mixing rules to calculate parameters a and b. There are several mixing rules available including Wong-Sandler (WS) and Modified Huron-Vidal second-order model (MHV2), and they are discussed by Yang et al. [53]. Here, the Vander Waals (vdW) mixing rule is discussed and used to calculate a and b parameters in a mixture, as follows [53]: where NC and x stand for the number of components in the mixture and mole fraction, respectively, and the value of a i and b i can be calculated using Equations (41) and (42). k ij also represents the binary interaction parameter between i and j species. This parameter can be calculated as a function of temperature by minimizing the difference between experimental and calculated freezing temperature as discussed by reference [54]. It is also worth noting that ASPEN Plus ® V9.0 (AspenTech [55]) can also be used to obtain this parameter. Once the empirical parameters are calculated using Equations (44) and (45), the following equation can be applied to obtain the fugacity coefficient of a specific component in the mixture [53]: where z in Equation (46) represents the compressibility factor. Equation (49) represents the solid-vapor equilibrium and can be used to calculate the freezing point of CO 2 [54]: where y CO 2 , ϕ V CO 2 , P, P Sat CO 2 Solid , ϕ Sat CO 2 , and v CO 2 Solid are the mole fraction of CO 2 in the vapor phase, the vapor phase partial fugacity coefficient for CO 2 , the system pressure, solid CO 2 vapor pressure, the fugacity of pure CO 2 vapor at P Sat CO 2 Solid , and the molar volume of solid CO 2 , respectively. However, in the case of liquefying CO 2 and separating it from the gas mixture, the following VLE relationship can be used to calculate the dew point of CO 2 [56]: where y CO 2 and x CO 2 represent the mole fraction of CO 2 in vapor and liquid phases, respectively. It has been reported that the PR EoS is more accurate than SRK (Soave-Redlich-Kwong) and RK (Redlich-Kwong) equations of states and the calculated results fit the experimental data with higher accuracy [53]. However, the accuracy of this model decreases as the pressure increases. As discussed above, PR EoS or other equations of state such as SRK EoS along with a mixing rule are mostly used to calculate the solid-vapor equilibrium and liquid-vapor equilibrium of CO 2 mixtures. In other words, the equation of state and the selected mixing rule are used together in order to perform the phase equilibrium analysis of CO 2 in a gas mixture and predict the CO 2 dew point and freezing point.

Flash Model Based on the Rachford-Rice Equation
The Rachford-Rice equation, which is mostly used for the flash calculation, can be used to model SVE conditions and the separation of CO 2 from a gas mixture in the solid form. In this case, the feed gas with a molar flow rate of F and composition of z which may contain N 2 , O 2 , and CO 2 enters a flash tank. As a result, solidified CO 2 with a molar flow rate of S leaves the tank from the bottom, while uncondensed gases with a molar flow rate of V and composition of x V leave the tank from the top. The CO 2 recovery for the separation process can be defined as follows when it is assumed that pure CO 2 is formed in the solid phase [43]: The following equations can also be used to discuss the mass balance in the separation process [43]: The Rachford-Rice equation, which is used for the SVE calculation, is given as follows: where z i and K i represent the mole fraction and equilibrium constant of i-th component in the feed gas. The latter parameter can be calculated as follows: As mentioned before, only CO 2 will be formed in the solid phase and because of that the equilibrium constant for other components rather than CO 2 would reach infinity. Under this situation, considering Raoult's law, the equilibrium constant for CO 2 can be calculated by the following equation: where P and P sub represent the system pressure and sublimation pressure, respectively. The variation in sublimation pressure against temperature can be calculated as follows where the temperature is in K and the sublimation pressure is in Pa [57]: Taking the limit K i → ∞ for i = CO 2 , the Rachford-Rice equation takes the following form which can be used to calculate the SVE temperature [43]: Comparing Equation (59) with Equations (44)- (49) demonstrates the simplicity of the model. Figure 5 also compares the calculated SVE temperature using the Rachford-Rice equation and those obtained by Baxter et al. [18] when the CO 2 recovery is 90%. As can be seen, the model is accurate, especially at low pressures (<15 bar). It is worth noting that the CCC process desublimates CO 2 at low pressures and normally it processes the flue gas at ambient pressure [58]. Therefore, the Rachford Rice model can be used to calculate the SVE temperature with sufficient accuracy for cryogenic carbon capture applications.   The Rachford-Rice model, which is developed by Guido et al. [43] for desublimation and the separation of CO2 in solid form, is extended here to calculate the liquid-vapor equilibrium temperature and separation of CO2 from the flue gas in liquid form. In this case, it is also assumed that pure CO2 is formed in liquid form. Therefore, the equilibrium constant for CO2 can be calculated using Equation (60), while the equilibrium constant for other components tends to be infinite.
where represents the saturation pressure of CO2 and is given as follows: Therefore, the LVE temperature can be calculated by using Equations (60) and (61) in Equation (59). Figure 7 indicates the LVE temperature, which is calculated for different CO2 recovery values using the Rachford-Rice equation. To verify the model, the feed gas  The Rachford-Rice model, which is developed by Guido et al. [43] for desublimation and the separation of CO2 in solid form, is extended here to calculate the liquid-vapor equilibrium temperature and separation of CO2 from the flue gas in liquid form. In this case, it is also assumed that pure CO2 is formed in liquid form. Therefore, the equilibrium constant for CO2 can be calculated using Equation (60), while the equilibrium constant for other components tends to be infinite.
where represents the saturation pressure of CO2 and is given as follows: Therefore, the LVE temperature can be calculated by using Equations (60) and (61) in Equation (59). Figure 7 indicates the LVE temperature, which is calculated for different CO2 recovery values using the Rachford-Rice equation. To verify the model, the feed gas The Rachford-Rice model, which is developed by Guido et al. [43] for desublimation and the separation of CO 2 in solid form, is extended here to calculate the liquid-vapor equilibrium temperature and separation of CO 2 from the flue gas in liquid form. In this case, it is also assumed that pure CO 2 is formed in liquid form. Therefore, the equilibrium constant for CO 2 can be calculated using Equation (60), while the equilibrium constant for other components tends to be infinite.
where P sat represents the saturation pressure of CO 2 and is given as follows: Therefore, the LVE temperature can be calculated by using Equations (60) and (61) in Equation (59). Figure 7 indicates the LVE temperature, which is calculated for different CO 2 recovery values using the Rachford-Rice equation. To verify the model, the feed gas composition was selected according to Yang et al. [53], and for 20-bar pressure, the calculated LVE temperatures were compared with the reference. As can be seen, the model predicts the LVE temperature with acceptable accuracy. However, for higher pressures (>30 bar), the model is not very accurate. composition was selected according to Yang et al. [53], and for 20-bar pressure, the calculated LVE temperatures were compared with the reference. As can be seen, the model predicts the LVE temperature with acceptable accuracy. However, for higher pressures (>30 bar), the model is not very accurate. In fact, at 2 MPa pressure when the concentration of CO2 in the gas mixture is less than 21.7%, CO2 will not form a liquid phase at any temperature [53]. In this case, the dew point of CO2 is lower than the desublimation temperature. Therefore, CO2 forms a pure solid phase rather than a liquid solution.
Like the previous equation of states, the Rachford-Rice equation can also be used to calculate solid-vapor equilibrium temperature and liquid-vapor equilibrium temperature for an ideal system. However, this model is not recommended for high pressures (>15 bar) since the assumption of the ideal system is not valid anymore and Equation (59) is not enough to describe the physics of the problem.

Minimum Energy Model
This model can be used to calculate the minimum energy required for the separation of CO2 from the flue gases in the cryogenic process. Therefore, this model describes an ideal cryogenic carbon capture process that does not have any energy loss in the form of work or heat [59]. In this model, as Equation (62) indicates, the minimum energy required to separate a mixed gas stream into constituent pure components considering the ideal gas assumption equals to the difference in Gibbs Energy between the mixed and separated gas streams.
where , , T, , and R represents the minimum energy required for separating the mixed stream gas into its pure components, the Gibbs energy of mixing, temperature, the entropy of mixing, and universal gas constant, respectively, while subscription denotes the mole fraction of ith species. Therefore, for the separation of CO2 from the flue gases with a specified capture percentage, the minimum energy required for capturing CO2 can be presented as follows: In fact, at 2 MPa pressure when the concentration of CO 2 in the gas mixture is less than 21.7%, CO 2 will not form a liquid phase at any temperature [53]. In this case, the dew point of CO 2 is lower than the desublimation temperature. Therefore, CO 2 forms a pure solid phase rather than a liquid solution.
Like the previous equation of states, the Rachford-Rice equation can also be used to calculate solid-vapor equilibrium temperature and liquid-vapor equilibrium temperature for an ideal system. However, this model is not recommended for high pressures (>15 bar) since the assumption of the ideal system is not valid anymore and Equation (59) is not enough to describe the physics of the problem.

Minimum Energy Model
This model can be used to calculate the minimum energy required for the separation of CO 2 from the flue gases in the cryogenic process. Therefore, this model describes an ideal cryogenic carbon capture process that does not have any energy loss in the form of work or heat [59]. In this model, as Equation (62) indicates, the minimum energy required to separate a mixed gas stream into constituent pure components considering the ideal gas assumption equals to the difference in Gibbs Energy between the mixed and separated gas streams.
where E min , ∆G mix , T, ∆S mix , and R represents the minimum energy required for separating the mixed stream gas into its pure components, the Gibbs energy of mixing, temperature, the entropy of mixing, and universal gas constant, respectively, while subscription x i denotes the mole fraction of ith species. Therefore, for the separation of CO 2 from the flue gases with a specified capture percentage, the minimum energy required for capturing CO 2 can be presented as follows: where η stands for CO 2 recovery. Figure 8 indicates the minimum energy which is needed to separate CO 2 from a gas mixture for different temperatures and CO 2 recovery values when the mole fraction of CO 2 in the flue gas is 13%. As can be seen, reducing the temperature results in decreasing the minimum separation energy, while increasing the CO 2 recovery results in increasing the minimum separation energy, especially when it increases from 0.9 to 0.99. Figure 8 justifies very low energy consumption in the CCC process as it separates CO 2 at low temperatures.
where stands for CO2 recovery. Figure 8 indicates the minimum energy which is needed to separate CO2 from a gas mixture for different temperatures and CO2 recovery values when the mole fraction of CO2 in the flue gas is 13%. As can be seen, reducing the temperature results in decreasing the minimum separation energy, while increasing the CO2 recovery results in increasing the minimum separation energy, especially when it increases from 0.9 to 0.99. Figure 8 justifies very low energy consumption in the CCC process as it separates CO2 at low temperatures. The minimum amount of work required for compressing the separated CO2 in the form of solid or liquid can also be calculated as follows: where , is the minimum compression work, while P represents pressure and denotes the volume of the compressed substance from state 1 to state 2. It is also worth noting that in the case of near-isochoric solids and liquids, the amount of compression work would be very low and with a good approximation, it is negligible as it is less than 1 kJ/kg for CO2.

Soave-Redlich-Kwong (SRK) Model
SRK is one of the most popular equations of states in the hydrocarbon industry [60]. It can be used to predict a solid-vapor and liquid-vapor equilibrium with adequate accuracy. However, SRK is not a perfect method for predicting liquid compressibility. Although this model has been modified and improved for different conditions and applications since it was proposed by Soave in 1972 [61][62][63], the standard form of SRK EoS is given as follows [50]: The minimum amount of work required for compressing the separated CO 2 in the form of solid or liquid can also be calculated as follows: where W com,min is the minimum compression work, while P represents pressure and υ denotes the volume of the compressed substance from state 1 to state 2. It is also worth noting that in the case of near-isochoric solids and liquids, the amount of compression work would be very low and with a good approximation, it is negligible as it is less than 1 kJ/kg for CO 2 .

Soave-Redlich-Kwong (SRK) Model
SRK is one of the most popular equations of states in the hydrocarbon industry [60]. It can be used to predict a solid-vapor and liquid-vapor equilibrium with adequate accuracy. However, SRK is not a perfect method for predicting liquid compressibility. Although this model has been modified and improved for different conditions and applications since it was proposed by Soave in 1972 [61][62][63], the standard form of SRK EoS is given as follows [50]: where V m denotes molar volume. Like PR EoS, when SRK EoS is used to calculate the equilibrium conditions in a mixture, mixing rules should be used to calculate a and b. Therefore, these two parameters can be calculated as follows when the van der Waals mixing rule is applied: where the value of a i and b i can be calculated for each component in the mixture as follows: The parameter m in Equation (68) can be calculated as a function of an acentric factor for each component as follows: The binary interaction parameter used in Equation (66) is given by Shen et al. [50] for several pairs including CH 4 , CO 2 , N 2 , and C 2 H 6 which are derived from HYSIS. Once a and b are determined using Equations (66)- (70), SRK EoS can be used to calculate the fugacity coefficient for each component as follows [53]: where A and B can be calculated by Equations (47) and (48). Once the fugacity coefficient is calculated for each component in the mixture, Equations (49) and (50) can be used for liquid-vapor equilibrium and solid-vapor equilibrium calculations.

ASPEN Plus Process Simulation Software
ASPEN Plus is the leading chemical process simulator and has been widely used to build process models. It can be used to improve the existing processes or design new processes. This software has also been used to model the CCC process and predict solidvapor and liquid-vapor equilibria using the RGibbs reactor block [8,23,64,65]. In this case, the RGibbs reactor can predict the equilibrium conditions by minimizing the Gibbs energy [65,66]. This software can also perform exergy analysis and calculate the exergy loss in the process. Moreover, the energy sensitivity and genetic algorithm available in ASPEN Plus can also be used for system optimization. Table 1 compares the accuracy of RGibbs block in the ASPEN Plus simulator for the calculation of CO 2 frost point temperature with the experimental data from the different literature. As can be seen, this software predicts the equilibrium conditions accurately and can be used for modeling the CCC process.

LNG Storage Tank
As mentioned before, the CCC process can store cooling energy in the form of LNG during the off-peak hours when extra energy is available in the power plant. The stored energy can be recovered with a very fast response during the peak hours to supply more electric power and meet the residential demands. Although high-quality insulation used in the storage tank can limit the ambient heat input into the LNG tank significantly, trivial heat leak from the surrounding is always unavoidable due to the temperature gradient between LNG and the ambient [71]. It is worth noting that the value of heat leakage is dependent on the structure of the tank. For instance, smaller tanks with lower storage capacities lose more cooling energy to the environment since they have the large surface area to volume ratio [72]. The heat leakage results in minor evaporation of the liquid and causes boil-off gas (BOG) generation in the tank. Moreover, heat leakage increases the temperature of the vapor phase faster than the liquid phase since the vapor phase possesses lower heat capacity than the liquid phase [73]. This temperature difference between the two phases inside the tank brings about heat conduction at the interface and stratification phenomenon [74]. The minor heat penetration from the ambient to the LNG tank also contributes to increasing the temperature of the liquid phase which is in contact with the tank walls, making it warmer than part of the liquid located in the center of the tank. This temperature gradient between the liquid phase in the center of the tank and the liquid phase close to the tank walls causes a local convection heat transfer inside the tank. As a result of local convection heat transfer, heated liquid rises to the surface and results in increasing the surface temperature and intensifying the BOG phenomenon by elevating the evaporation rate, thermal stratification, and vapor pressure.
The generation of BOG, which is discussed above, can increase the pressure inside the storage tank, and in some cases, it may cause safety problems if the BOG is not removed from the tank for a long time and vapor pressure inside the tank passes the safety limits. It has been shown that the amount of BOG formation is highly affected by design and operating conditions of the LNG tanks [71]. Hence, the modeling and simulation of storage tanks is very important since the available experimental data are very limited [73].

Volume of Fluid (VOF) Method
Modeling the storage of LNG in a full-scale tank requires multiphase flow models, and because of that, the VOF method is selected and discussed here. The VOF model is the Eulerian method, of which is a very popular and powerful technique for tracking the vapor-liquid interface and can be applied to a fixed Eulerian mesh [75]. This method can be used for two or more immiscible phases to determine the intersection position between the fluids. The formulation of this method is based on the fact that two or more fluids do not interpenetrate each other. This model solves a single set of momentum equations in each control volume and estimates the value of volume fraction for each phase or fluid in the cells. The value of the volume fraction in each cell for a specific phase can be 1 or 0, or a value in the range of 0-1 when the control volume is completely occupied by the phase, it is vacant of the phase, or it is partially filled with the phase, respectively. Therefore, according to Equation (72), the summation of volume fractions of all phases in each cell should be unity.
where α is the volume fraction, while i and n represent the ith phase and the total number of phases, respectively. The interface only appears if a control volume is partially filled with pth phase. When this method is used to track the interface, the variation in density is just considered since the temperature varies slightly and because of that other fluid properties are considered invariable [76]. Herein, the variation in density against the temperature is estimated using the Boussinesq approximation in conservation equations for mass, momentum, and energy, given as follows: ∂ρ ∂t where ρ, → V, p, → g , and k represent density, velocity vector, pressure, gravity vector, and conductivity, respectively, while T and T 0 denote local temperature and reference temperature, respectively. → F is also body force generated by surface tension at the interface and S h also denotes the energy which causes the phase change. It should be noted that in Equation (75), the effect of the gravitational field is included in the pressure term. Furthermore, in this equation, E represents the energy and can be calculated by mass-average using the following equation [77]: where subscription q denotes the qth phase, while n represents the number of phases. The value of energy for each phase, i.e., E q , can be calculated using the temperature and specific heat of each phase. Moreover, the average method used for calculating the energy can be also used for calculating the temperature using the following equation: In Equation (74), β also stands for thermal expansion coefficient and can be calculated using the following equation: To track the interface between two phases, the continuity equation should be solved for volume fraction of the second phase as follows: where . m is the phase change flow rate resulting from evaporation or condensation at the interface and it gets positive and negative values for evaporation and condensation processes, respectively. Subscription v also denotes vapor phase. In this case, the phase change flow rate can be presented using the Lee phase change model, as follows [78]: .
where T sat is saturation temperature which changes by pressure inside the tank according to Clausius-Clapeyron equation, while r is the mass transfer intensity factor. The variation in the former parameter against the pressure can be measured for liquid methane, as follows: Considering an appropriate value for r has a great impact on converging the solution by keeping the interfacial temperature close to the saturation temperature [78]. According to the previous CFD studies, 0.1 s −1 could be the best value for r, which may guarantee a stable solution [79][80][81]. However, there are several studies which considered different values for r. For instance, 100 s −1 was also considered for r by [82,83]. Once the value of the phase change flow rate is obtained using Equations (80) and (81), energy source in Equation (75) can be calculated as follows: where L H denotes the latent heat of evaporation.
In this method, material properties that appear in transport equations are calculated based on the volume fraction of the component phases in each control volume. In twophase liquid-vapor flow systems, these properties can be determined as follows if the volume fraction of the liquid phase is tracked [77]: The finite volume method can be used to discretize the partial differential equations into a system of algebraic equations. These algebraic equations can be numerically solved to provide the solution field. ANSYS FLUENT can be used as a commercial solver software for the simulation of LNG tank using VOF method discussed above [71,84]. This commercial software utilizes the finite volume method for the discretization of partial differential equations. In this case, the SIMPLE algorithm can be used for pressure-velocity coupling, while QUICK schemes are used for solving the convective terms [84]. Moreover, time discretization can be performed using a First-Order Implicit scheme and the volume fraction equation can be solved using the Geo Reconstruct scheme. The value of α l in Equations (84)- (87) can be set if the filling level is known.
To set the required boundary conditions for the simulation, a convective heat transfer with air which is surrounding the tank can be imposed on the exterior boundary [84]: where T ext represents external temperature. The value considered for the heat transfer coefficient in Equation (88) is also 7 W m 2 K , while k also represents the thermal conductivity of evacuated perlite used as thermal insulation and can be calculated using the following equation for different temperatures in Kelvin: On the inner wall of the tank shown which represents the interface between solid and liquid subdomains, the nonslip condition (v = 0) can be imposed. On this boundary, the continuity of heat flux and temperature should also be satisfied. To reduce the computational costs, axis boundary conditions may be imposed on the symmetry axis. To set initial conditions for the fluid subdomain, initial temperature and vapor pressure can be set to 113 K and 1.13 bar, respectively, while the initial velocity for both phases should be set to 0.
The numerical solution can be obtained in a two-step process. First, a steady state problem is solved to obtain a temperature profile tank wall by setting the tank internal wall temperature to 113 K. The calculated temperature profile is then used as the initial temperature of the wall subdomain to solve the transient problem. In this case, a smaller time step is recommended since it contributes to solution stability and convergence. Therefore, the first 20 time steps can be solved in fluid and wall subdomains with ∆t = 0.001 s as the time step using the temperature profile calculated in the steady-state problem for the wall subdomain. Afterward, the value of the time step can be increased to ∆t = 0.01 s to reach the end time of simulation.
The results of previous simulations using the VOF model in ANSYS FLUENT commercial solver have distinctly proven the precision of the model by predicting the BOG generation rate and temperature distribution inside the tank with high accuracy [85]. In other words, it has been shown that the simulation results of the 3D model discussed above are almost the same as the experimental results, and using a 2D planar model has also an acceptable accuracy with almost 5% error [86].

BWRS Model
Benedict-Webb-Rubin Equation of State (BWR EoS) can be used to model and investigate the LNG tank in the CCC process. In other words, it can be used to calculate density and enthalpy, which are two critical parameters in the calculation of BOG [72]. This model is very convenient for computer programming. In this case, the equation of state is given as follows: where ρ, T, P are molar density, temperature, and pressure, respectively, while B 0 , A 0 , C 0 , D 0 , R, γ, b, a, and α are coefficients that can be found in [87]. This equation has been modified by Kenneth Starling and has resulted in the following equation called the BWRS model [72]: A comparison between Equations (90) and (91) indicates that the BWRS model possesses one more coefficient, i.e., d, and because of that, the modified model has higher accuracy than the BWR model. Equation (91) has also been revised by Nishiumi et al. [88], which resulted in the following equation of state: Although Equation (92) is more complicated than Equation (91), previous studies have indicated that errors in predicting thermodynamic properties can be reduced from 10 to 0.35% or lower when Equation (92) is used at reduced temperatures [88].
For convenience, the following substitutions can be made to solve Equation (92) for density: As a result, Equation (92) will be changed to the following equation: Therefore, Equation (93) can be solved for density using the bisection method or the Newton-Raphson method.
Once density is calculated for LNG, enthalpy can also be obtained, as follows [88]: where H 0 is enthalpy in the ideal gas state. After the calculation of density and enthalpy using Equations (93) and (94), the value of BOG can be calculated using the following equation: where V, m, and t represent the volume of the tank, rate of LNG evaporation, and time span for measuring the BOG. The following equation can be used to calculate the rate of LNG evaporation: where ∆h represents enthalpy change, while Φ is the tank heat leakage. The results of the BWRS EoS model indicate that increasing the percentage of methane in LNG intensifies the BOG generation [72]. Moreover, for different methane percentages in the LNG composition, smaller tanks have higher BOG generation since they possess a higher surface area to volume ratio. Although the method discussed here for evaluating the LNG storage tank and BOG calculation assumes a constant temperature and density for LNG and it might not be as precise as the VOF method discussed before, it is very simple and has acceptable accuracy. As result, when details such as stream function counters, fluid velocity inside the tank, and streamlines are not needed, they can be used to evaluate the LNG storage tank in the CCC process and calculate the BOG.

Compressor
The compressor is one of the most critical components in the CCC process, consuming over 80% of the required energy in the process for compressing CO 2 , natural gas, and mixed refrigerant [8]. Therefore, investigation and modeling of compressors are of paramount importance. Single-shaft turbo-compressors have been modeled and optimized in reference [89] and can be used for LNG production and cryogenic applications. However, in this part, the available models for evaluating scroll compressors and linear compressors are discussed since they are more efficient than conventional compressors. These two types of compressors, especially the novel oil-free linear compressors, are widely used to compress refrigerants in cryogenic applications.

Compression Process Model of Scroll Compressors
In this model, variations in temperature, mass, and pressure in different chambers can be calculated as functions of an orbiting angle. Therefore, the variation in refrigerant temperature against the orbiting angle can be specified using the following equation [90]: where m, C V , θ, υ, ω, h, V and .
Q represent the mass of refrigerant, specific heat at constant volume, scroll orbiting angle, specific volume, angular speed of compressor shaft, specific enthalpy of refrigerant, volume of the control volume, and heat flow rate, respectively. In this equation, . m in denotes the mass flow rate which enters the control volume. Based on the refrigerant mass conservation, the mass balance in the compressor is provided by the following equation [90]: Equations (97) and (98) are both valid for real gases [90]. In fact, the first-order differential equation given by Equation (97) has two autonomous variables, i.e., temperature and mass. This equation needs to be integrated numerically to reach a solution. In this case, expressions used to take into account the variation in pressure with the temperature at constant specific volume and enthalpies can be obtained from an equation of the state. The required expressions in this equation including the variation in pressure with respect to temperature at a constant specific volume and enthalpies can be derived from the equation of state. However, Equation (97) requires models which can predict the mass flow rates which enter/exit each chamber and the heat transfer rate. Considering Equation (97) along with Equation (98) and an equation of state, thermophysical properties can be calculated for different values of orbiting angles in whole compression chambers.
To model the compression process for the scroll compressor, the first-order differential equations given by Equations (97) and (98) should be solved in a step-by-step procedure. Therefore, the Euler method discussed by Conte et al. [91] can be used to solve equations at any desired orbiting angle θ j , as follows: where ∆θ represents the step width. To solve the above-mentioned first-order equations and obtain the pressure and temperature as a function of orbiting angle, it should be noted that several differential equations should be solved at the same time since there are several compressor chambers, and each compressor chamber possesses its own differential equations.
The following equation can be used to calculate the real amount of compressional work which needs to be carried out on the refrigerant [90]: where h dis , h suc , . Q avg,alum , and . Q avg,steel represent the average enthalpy of discharge gas, enthalpy of suction gas at suction pressure and suction temperature, and average heat transfer rates from steel scroll and the aluminum scroll to the refrigerant, respectively. In this case, the average heat transfer rates are given as follows: . Q alum,j denote the heat transfer rate at step j which can be calculated as discussed by Chen et al. [91]. The average enthalpy of discharge gas which is needed in Equation (101) can also be calculated as follows: where h dis,j and . m dis,j are the enthalpy and mass flow rate (kg/radiant) of the discharged refrigerant at the orbiting angle of θ j , respectively.
The isentropic compression work can also be calculated as follows [92]: where h isentropic represents the isentropic enthalpy of discharge refrigerant. This parameter can be calculated based on the real gas properties and considering the inlet entropy and discharge pressure [92]. The isentropic compression work and the real compression work given by Equations (101) and (105), respectively, can be used to obtain the isentropic efficiency of the compressor, as follows: The overall compressor efficiency can also be calculated using the following equation: where P is the power input to the compressor and can be calculated as follows: W compression η motor−mechanical (108) η motor−mechanical used in Equation (108) represents motor-mechanical efficiency and takes into account both mechanical losses and motor losses. This parameter can be calculated by regression carried out on the performance testing results of the compressor and using the following equation: where a 1 , a 2 , a 3 , and a 4 are constants in Equation (109).

Linear Compressor Model
The linear compressor can be successfully used in cryogenic applications. In a linear compressor, the piston moves along a linear track, and because of that, the friction and energy loss are very low. The variation in internal energy of the linear compressor is given as follows [93]: where . Q step , . m in , . m out , h in , and h out denote the heat transfer between the cylinder wall and refrigerant, inlet mass flow rate, outlet mass flow rate, and specific enthalpy of the refrigerant at the compressor inlet and outlet, respectively.
. W sha f t also represents the shaft power and can be calculated as follows: .
where P c and V c are the pressure in the cylinder and volume, respectively. The variation in in-cylinder temperature can be calculated using the following equation: where m c represents the mass of refrigerant in the compressor, while T c denotes the temperature of the refrigerant in the cylinder. The motor efficiency of linear compressors is usually assumed to be 90% for numerical studies [93], while the volumetric efficiency for linear compressors can be calculated as follows [94]: where . m exp , ρ suction , x P , and A P represent the experimental mass flow rate of refrigerant, the density of suction gas, displacement of the piston, and area of the piston, respectively. The overall isentropic efficiency can also be defined using the following equation [95]: where h 2,s and h 1,s denote the enthalpy of refrigerant in discharge and suction stages, respectively, while P is the input power to the system.

Accuracy Comparison among Models of Different Component in the CCC Process
In this investigation, the available knowledge for modeling and simulating the main components in the CCC process is presented. Table A1 summarizes the explored models and methods for different components in the CCC process. The accuracy of each model is given in terms of the average absolute deviation (AAD%) in the reviewed references. This parameter can be obtained for the desired variable when the calculated and experimental value at the ith point is represented by variable cal,i and variable exp,i . Although all presented models are accurate, they possess different AAD values due to the assumptions made in the models and the utilized method for modeling. This parameter is provided by some references and calculated for others in the present investigation.
As Appendix A indicates, the differential method is mostly used to design and model the heat exchangers in the CCC process. According to this table, ASPEN Plus is a very popular commercial software to model the CCC process integrated with power plants or other industrial processes due to its exceptional capability to model thermal and chemical processes. The RGibbs reactor available in Aspen Plus can be used to accurately calculate equilibrium conditions as well as the dew point and frost point temperatures of CO 2 , solid CO 2 solubility in hydrocarbon mixture, etc. Furthermore, ANSYS FLUENT as a CFD tool has also been widely used by researchers to model and evaluate the LNG storage tank. They used the VOF method in 2D and 3D dynamic models to evaluate the BOG phenomenon with acceptable accuracy. Furthermore, it was shown that the flash model based on the Rachford-Rice equation can be used to predict the SVE conditions with acceptable accuracy, especially at low pressures. This model was also extended to predict LVE conditions and calculate the dew point of CO 2 . It was shown that in both cases, this simple model can be applied to predict the equilibrium conditions at low pressures without involving the complexity of mixing rules.

Conclusions
Although the CCC process is a new technology for CO 2 emission mitigation, different aspects of this process have been modeled and simulated using numerical methods, but several fields need to be investigated more, as follows:

•
A numerical model which can predict the pressure drop in the process accurately may contribute to reaching an optimum size for tubes and heat exchangers to reduce the pressure drop without the generation of a maldistribution phenomenon. • A dynamic model of the power plant and the CCC process, even in the case of using the integral method for modeling the heat exchangers, would provide helpful information about the transient behavior of the process. Funding: This research received funding from Energy Cluster Denmark for COLDCO2 project. More information about the funding and project can be found in https://www.energycluster.dk/en/ projects/acomar-2/ (accessed on 14 November 2022).

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article. was not involved in study design, in the collection, analysis and interpretation of data. We thank Larry Baxter (Brigham Young University and Aalborg University) for his great support and helpful information about the modelling of the CCC process.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Modeling the overall energy consumption in the CCC process based on Stirling cooler [17] A steady state model using ASPEN Plus 6.6% Simulation of the CCC process using energy storage system [8] A steady state model using ASPEN Plus 1% [22] Modeling the energy consumption in the process and sizing the required heat exchangers [22] 1-D pseudo homogeneous axially dispersed plug flow model 2.2% Dynamic simulation of packed beds used in the CCC process. [46,47] ASPEN HYSYS commercial process simulator -Modeling four different natural gas liquefaction processes to supply the required LNG in the CCC process [23]  Modeling a small-scale natural gas liquefaction process [24] Differential method -3D simulation of multi-stream heat exchangers [36] A 3D model based on segmented differential method -Modeling multi-stream heat exchangers [41] Differential method -Modeling and sizing multi-stream heat exchangers [42] A 2D model using VOF method and ANSYS FLUENT commercial process simulator 6.5% [71] Modeling large-scale LNG storage and carbon capture systems [71] BWRS EoS model -Modeling BOG generation in LNG tanks [72] A 2D model using VOF method and ANSYS FLUENT commercial process simulator 4% [78] Modeling the evaporation and condensation phase-change processes [78] A 3D model using Eulerian multiphase flow method and ANSYS FLUENT commercial process simulator 6.2 -Modeling refrigerant flow boiling in horizontal serpentine tube [79] A 2D model using VOF method and ANSYS FLUENT commercial software, version 19 -Simulation of natural convection and BOG generation in small, pressurized LNG storage tank [84] A 3D model using VOF method and ANSYS FLUENT commercial process simulator 10% [85] Simulation of BOG generation rate for LNG tank [85] A 2D model using VOF method and ANSYS FLUENT commercial process simulator 5% [86] Simulation of BOG generation rate for LNG tank [86] Compression Process Model -Modeling scroll compressors which can be used in the CCC process [90] An improved BWRS EoS model 0.35% [88] Predicting thermodynamic properties at extreme low reduced temperature [88] Modified SRK EoS model 1% [63] Predicting the vapor-liquid equilibrium in binary asymmetric mixtures [63] Minimum energy model -Minimum energy required for separating CO 2 from the flue gas [59] Classical approach based on PR EoS model and ASPEN Plus commercial process simulator V9 Predicting multi-phase equilibrium conditions of CO 2 mixtures with hydrocarbon and non-hydrocarbon components and presenting stability analysis [66]  Simulation and economic analysis of a cryogenic carbon capture process assisted by a solar absorption refrigeration system [64] Rachford-Rice flash model 1.27% in the pressure range of 1-15 bar SVE calculation and predicting the frost point temperature of CO 2 [43] Linear compressor 2.47% and 8.49% for prediction of mass and power [93] Modeling two-stage linear compressor [93] ASPEN Plus -SVE calculations for CO 2 capture application [96] ASPEN HYSIS -Modeling and evaluation of a combined system of steam methanol reforming/steam methane reforming/high temperature PEM fuel cells/CO 2 capture/liquefaction [97] ASPEN HYSIS -Process design and thermoeconomic evaluation of a CO 2 liquefaction process (Energy, exergy, and exergoeconomic analyses) [98] Transient model using ASPEN HYSIS and PR EoS -Modeling and optimizing natural gas liquefaction process which provides coolant for CCC process [99] Aspen plus model using SRK -Modeling hybrid membrane-cryogenic capture process [