Mathematical simulation of heat and mass transfer during controlled depressurization of supercritical CO2 in extraction vessels

J. of Supercritical Fluids 122 (2017) 43–51

Contents lists available at ScienceDirect

The Journal of Supercritical Fluids journal homepage: www.elsevier.com/locate/supflu

Mathematical simulation of heat and mass transfer during controlled depressurization of supercritical CO2 in extraction vessels ˜ b M. Soledad Murias a , José M. del Valle a,∗ , Gonzalo A. Núnez a

Department of Chemical and Bioprocess Engineering, Pontificia Universidad Católica de Chile, Avda. Vicu˜ na Mackenna, 4860, Macul, Santiago, Chile Department of Chemical and Environmental Engineering, Universidad Técnica Federico Santa María, Avda. Vicu˜ na Mackenna, 3939, San Joaquín, Santiago, Chile b

a r t i c l e

i n f o

Article history: Received 27 April 2016 Received in revised form 13 November 2016 Accepted 14 November 2016 Available online 15 November 2016 Keywords: Carbon dioxide Depressurization Heat transfer Mathematical simulation Packed bed

a b s t r a c t Even after almost forty years of industrial application, companies are still reluctant to use supercritical (sc) CO2 as a solvent for extractions due to the perceived high production costs. Literature on the matter suggests that, because extraction at high pressure needs to be done in batches, using multiple extraction vessels with simulated-countercurrent flow could reduce operational costs. However, the more extraction vessels used, the less time there is to recondition them in order to have a semi-continuous operation; and if the reconditioning, and particularly the depressurization, is done too fast, the vessel could become brittle and permanently damaged. With the goal of optimizing the depressurization process in mind, numerical simulation of temperature and mass was carried considering a 1-dm3 vessel filled with a packed bed made with model materials using a mass flow function that depends on the chocked mass flux and a valve opening area. The correlation Nu = 0.0777Da−0.0373 Ra0.397 was obtained for convective heat transfer at the vessel wall, and temperature, pressure, and vented mass flow were simulated with about 20% improvement in predictions in comparison to our previous correlation. To explore the use of the model for practical purposes, it was used to simulate depressurization processes with volumes up to 1 m3 and with different initial conditions and vessel geometries so as to have a first approach on the effect of these parameters on the depressurization time. Simulated depressurization times reached a maximum value of 54.5 min for depressurizations of a 1-m3 extraction vessel starting at 60 ◦ C and 70 MPa, which are very plausible extraction conditions. This model can be used to determine optimal reconditioning time in industrial plants for cost minimization. © 2016 Published by Elsevier B.V.

1. Introduction Solvent extraction using supercritical (sc) CO2 is currently a well stablished operation in the food industry, in processes such as decaffeination of coffee, removal of bitter flavors from hops, and extraction of high-value compounds [1]. Even though scCO2 has many advantages over other organic solvents, and has been used as an extraction solvent for almost four decades, companies still hesitate when it comes to investing in this technology [1]. A reason for these doubts may be the perceived high costs associated with operating an extraction plant. However, cost optimization could prove to be the tool needed in order to refute such beliefs. ˜ and del Valle [2] optimized production costs in a scCO2 Núnez extraction plant. Some of the relevant variables for optimization

∗ Corresponding author. E-mail address: [email protected] (J.M. del Valle). http://dx.doi.org/10.1016/j.supflu.2016.11.009 0896-8446/© 2016 Published by Elsevier B.V.

were mass flow rate (inversely proportional to production costs; production costs decreased when increasing mass flow rate), aspect ratio of the extraction vessel (directly proportional to production costs for a constant superficial velocity of CO2 ), and number of extraction vessels with simulated-countercurrent flow (inversely proportional to production costs). One limiting aspect of increasing the number of extraction vessels is the reconditioning time. For semi-continuous operation, an exhausted batch must be replaced for a fresh one in less time than required another vessel to become exhausted [3]. This means that as the number of extraction vessels of the plant, the time available for reconditioning decreases if a semi-continuous operation is desired. The reconditioning process consists of four steps: (i) depressurization of the vessel, (ii) unloading of the exhausted material, (iii) loading of fresh material, and (iv) pressurization of the vessel. Of these stages, the least studied and the one with largest opportunities for optimization is the initial depressurization of the vessel. The focus of previous research has been mostly on safety aspects [4–6]

44

M.S. Murias et al. / J. of Supercritical Fluids 122 (2017) 43–51 Table 1 Operation conditions used for simulation.

Nomenclature

Temperature of heating water (◦ C) Initial temperature (◦ C) Initial pressure (MPa) Radium of extraction vessel (mm) Thickness of extraction vessel wall (mm) Height of extraction vessel (mm) Volume of extraction vessel (m3 )

Latin letters A Valve opening area (m2 ) Da Darcy number, dp2 ␧3 150−1 (1 − ␧)−2 L−2 (−) Particle diameter (m) dp F Mass flow rate (kg s−1 )h Heat transfer coefficient (W m-−2 K−1 ) Thermal conductivity (W m−1 K−1 ) k l Packed bed material height (m) Vessel height (m) L m System fluid mass (kg) Nusselt number, hLk−1 (−) Nu Pr Prandtl number, ˛−1 (−) r Radial position (m) Vessel internal radius(m)  R Ra Rayleigh number, gˇ Tw − T¯ L3 (˛)−1 (−) Re t T u x z

Although an improvement from previous work, the experimental design proposed by Richter et al. [9] had some deficiencies. We believe the number of thermocouples placed in the vessel to be insufficient to stablish the existence of axial or radial temperature profiles throughout the depressurization. Also, it is difficult to measure the temperature of the fluid near the vessel wall; mass flow could have caused vibrations or movement in the thermocouple placed in there that could have in turn affected measurements. Because of this, we believe Eq. (1) should be reevaluated through simulation and optimization. We believe it is possible to generate a mathematical model that can simulate the changes in temperature and pressure throughout the depressurization process. This model can be used to optimize the parameters of the correlation for convective heat transfer at the vessel wall, and to simulate the previous experimental results [9]. The model can be later used to obtain an optimal mass flow that will allow safe and fast depressurizations. Its ability to simulate depressurization processes at industrial scale will need to be assessed once experimental results are available at those conditions. By providing tools for optimizing the extraction process, our aim is to reduce uncertainties that prevent improving the productivity of industrial scCO2 extraction processes, which in turn precludes this technology from being a widespread choice for the food industry.

Reynold number, udp −1 (−) Time (s) Temperature (◦ C or K) Superficial velocity (m s−1 ) Vessel steel wall thickness (m) Axial position (m)

Greek letters ␣ Thermal diffusivity (m2 s−1 ) ␳ Density (kg m−3 ) Porosity ␧  Speed of sound (m s−1 ) Subscripts eff Effective f Fluid J Jacket Lid l s Solid t Time t w Wall

2. Model development

and not on the mechanics of the depressurization or how they affect industrial operation. The one practical limitation that this process has is that, when done at high speeds, it results in low temperatures that turn the vessel become brittle [7] and damage it permanently. Richter et al. [8] studied the controlled depressurization of an extraction vessel filled with pressed or pelletized organic material at pilot scale in order to identify the principles that rule this process. Their findings show that total packed bed porosity was fundamental in determining temperature and pressure changes within the vessel during depressurization, and that all the components of the extraction vessel (substrate, walls, and lids of the extraction vessel, CO2 ) were relevant for the heat transfer in the system. In a follow-up work, Richter et al. [9] studied this same problem using an extraction vessel with a well-defined geometry and model materials made from glass or sintered steel. The results of this work were the validation of previous results [8] and a new correlation to estimate the convective heat transfer coefficient at the vessel wall (Eq. (1)). This correlation was better than the one proposed before because it included the Darcy (Da) dimensionless number as a parameter, which represents the relative effect of permeability of the substrate. This allowed for one equation for different materials instead of one for each. Nu = 0.086Da−0.037 Ra0.357

60 60 26 38 18.5 230 0.001

(1)

This section is divided in three parts. The first one shows the assumptions and considerations made in order to simplify the model. The second part shows the equations used to model heat and mass transfer in the extraction vessel. Finally, we describe the numerical method used, how error and confidence intervals were calculated, and how we tested the effect of some operational parameters on depressurization processes. 2.1. Model assumptions In order to produce a model for controlled CO2 depressurization, several assumptions had to be made. First, differences in pressure within the extraction vessel were considered negligible. Physical properties of the packed bed (bulk density and porosity) and of the AISI 316 steel wall (thermal conductivity and heat capacity) were assumed to be constant during depressurization because of the small changes in them in the temperature interval between 250 and 350 K. The thermal diffusivity of steel was obtained from Bergman et al. [10]. Experimental results obtained by Richter et al. [9] were used to validate the results obtained through simulation. These experiments were carried out in a 1 dm3 extraction vessel with a heating jacket and with experimental conditions shown in Table 1. The vessel was loaded with a packed bed made from three possible model materials: sintered stainless steel cylinders, glass beads, and glass Raschig rings. Experimental depressurizations with packed beds were carried out in two stages depending on the valve opening. Initially, the valve opening was small (flow coefficient fc = 0.0005 provided by the needle valve manufacturer; flow coefficient is a measure of the efficiency of the valve in letting fluid flow) to avoid

M.S. Murias et al. / J. of Supercritical Fluids 122 (2017) 43–51

Energy balance for the vessel wall was calculated using the same equations used in Richter et al. [9]. Bulk porosity was calculated using the correlations proposed by Dixon [12] for spheres (Eq. (4a)) and cylinders (Eq. (4b)).

60

Temperature (ºC)

55

50

45

␧= 40

35

30 0

2

4

6

8

10

12

14

16

Time (min)

␧=

Fig. 1. Comparison between simulated temperatures using Eq. (14) for heat transfer through the wall (solid line), simulated temperatures using our previous correlation [9] for heat transfer through the wall (dotted line), and experimental results [9] (symbols) at the center of the vessel for depressurizations starting at 60 ◦ C and 26 MPa packed with sintered steel cylinders. Results using the correlation proposed in this work are closer to experimental results in the first part of the process but later reach higher temperatures than those obtained through experimental procedures.

phase separation due to experimental restrictions. Valve opening was increased (fc = 0.0086) when mass flow was no longer detectable by the flow meter. The geometry and dimensions of the extraction vessel can be found in Table 1, and the characteristics of the packed material used in the experimental design can be found in Richter et al. [9]. 2.2. Mass and heat transfer in the system A single vessel filled with CO2 at high pressure and temperature and a packed bed made from model materials was considered. Also, the solid temperature was calculated separately from that of the fluid. The equations that represent the mass and energy balances along the vessel were adapted from the ones used by Slimi et al. [11]. The mass balance is given by Eq. (2). dm = −F (t) , dt

(2)

where F (t) = A␳␯, and F (t) is the product between the chocked mass flux at the inlet of the exit valve (calculated as the product between the fluid density ␳ and the speed of sound ) and the area of the opening of the valve (A), as suggested by Richter et al. [8]. The values of A were obtained experimentally. The components of the energy balances for the fluid (Eq. (3a)) and the solid (Eq. (3b)) are energy loss due to fluid motion, heat conduction within the solid phase, heat convection within the fluid phase expressed as effective conductivity, heat convection between solid and fluid phases, and energy loss due to mass being evacuated.

 ␣f 

u ∂Tf ∂Tf =− + ␧ ∂z ␧ ∂t



+ kfeff

 ∂2 Tf z

45

∂z 2



−

∂Ts ␣s kseff = 1−␧ ∂t

 1 ∂ eff

kf

r

r ∂r

 r

∂Tf ∂r

2

2

∂ Ts 1 ∂Ts ∂ Ts + + r ∂r ∂r 2 ∂z 2

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩



dp 0.4 + 0.05 + 0.412 2R



0.528 + 2.464

 1 − 0.667

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩



dp − 0.5 2R

dp 2R

3 

 0.677 − 9



 +

(3a)

1 h ␣ (T − Ts ) ks fs fs f



,

dp + 0.7 2R

1 − 0.763

dp ≤ 0.536, and 2R

dp 2R

, 0.536 ≤

2 ,

2

dp − 0.625 2R



dp ≤ 0.5, 2R

−0.5

dp 2R

, 0.6 ≤

(4a)

dp , 2R

dp ≤ 0.6, 2R dp ≤ 0.7, and 2R

2 , 0.7 ≤

(4b)

dp . 2R

The fluid-to-solid heat transfer (hfs ) coefficient was based on an empirical correlation proposed by Wakao et al. [13] (Eq. (5)) and the solid-fluid exchange area (␣fs ) was calculated using Eq. (6).





hfs =

kf dP

˛fs =

6 (1 − ε) for sphere, and dp

(6a)

˛fs =

4(d + l) (1 − ε) for cylinders. dl

(6b)

1/3

2 + 1.1 Prf

Ref0.6

(5)

,

where d and l are the diameter and height of the packed cylinder, respectively. Effective thermal conductivity is a dimensionless correction factor that can be applied to thermal conductivities to account for the irregularities of the system, such as the presence of a packed bed or multiple heat transfer phenomena occurring at the same time. Effective thermal conductivities for the fluid and solid phase were calculated using the correlations proposed by Dixon and Creswell [14] that are shown in Eq. (7).

⎛



kfeff

 r

=1+

ks kf

⎜ ⎝

(1 − ε) 1+

2 2R dp

16 k 3 s

1 hfs dp

⎛ 

kfeff

 z

=1+

ks kf

⎜ ⎝

(1 − ε)2 1+

16 k 3 s

1+ +

8kf 2Rhw 0.1 ks

4 2R dp

1 hfs dp

+

0.1 ks

⎞

 ⎟ ⎠,

(7a)

⎞ ⎟ 2 ⎠ , and

(7b)

(7c)

The initial and boundary conditions used to solve the mathematical model are shown in Eqs. (8) and (9), respectively. Initial conditions are not chosen arbitrarily but given by the extraction conditions; the depressurization starts just after the extraction ends, so the extraction temperature and pressure are the starting temperature and pressure of the depressurization. Tf |t=0 = Ts |t=0 = TJ , t = 0

(3b)

2

dp −1 R



0.36 + 0.1

dp 2R

, 0.5 ≤

kseff = (1 − ε) .



F (t) h hfs ␣fs (Tf − Ts ) − kf 2␲rRLkf



⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨

Tf |r=R =



hw hw + kw /xw

 Tf |r=0 ∀z, t

(8) (9a)

46

ks

M.S. Murias et al. / J. of Supercritical Fluids 122 (2017) 43–51

∂Ts ∂Tw = kw ∀ z, t , | | ∂r r=R ∂r r=R

(9b)

∂Tf ∂Ts = = 0 ∀ z, t , | | ∂r r=0 ∂r r=0 Tf |r=0 =

hw hw + kw /xl



Tf |

r=0

(9c)

∀ t andz = {0, L}, and

(9d)

Ts |r=0 =

hw hw + kw /xl



Ts |

r=0

∀ t andz = {0, L},

(9e)

where hw is the heat transfer coefficient between the steel wall and the fluid, calculated using a correlation like the one proposed by Richter et al. [9], shown in Eq. (10). These conditions come from balancing heat conduction from the vessel wall or lids and heat convection through the packed bed. Border condition (9a) represents heat transfer through the vessel wall, and border conditions (9d) and (9e) represent heat transfer through the upper and lower lids of the vessel, for the fluid and for the solid, respectively. The parameters of this correlation were obtained by minimizing the least square error between the experimental and simulated data.

 kf  aDab Rac L

(10)

The dimensionless Rayleigh (Ra), and Darcy (Da) numbers were calculated as proposed by Richter et al. [9] with the height of the vessel L as the characteristic length, using the following equations:



Ra = Da =



gˇ Tw − T¯ L3 ␯␣ ␬ L2

, and

(11a) (11b)

where ␣, ␤, and  are the thermal diffusivity, volumetric thermal expansion coefficient, and kinematic viscosity of the CO2 , respectively, and ␬ is the permeability of the packed bed, which we calculated using the equation of Carman-Kozeny [15], Eq. (12): =

dp2 ␧3e 150(1 − εe )2

(13)

i=1

z = h/2

hw =

1  |yi − f (xi )| n yi n

MAPE =

z = h/2

For optimization, only the results of sintered steel cylinders and glass beads were used. This was so, because a Raschig ring is not a commonly used shape of packed bed materials in the industry, unlike cylinders or spheres. In order to measure the difference between the simulated and experimental results, a Mean Absolute Percentage Error (MAPE) was calculated using Eq. (13).

,

(12)

where dp is the particle size or equivalent diameter and ␧e is the interparticle porosity. 2.3. Numerical method, parameter optimization, and sensibility Both solid and fluid phases were divided in 10 nodes in the radial direction and 10 nodes in the axial direction. This number was chosen because it was the smallest at which adding one more node made no difference in the results. The resulting equations were 201 in total: (i) one hundred differential equations to account for the heat transfer in the fluid; (ii) one hundred differential equation to represent the heat transfer in the solid phase of the packed bed; (iii) one equation to represent the mass transfer between the system and the surroundings. Central finite differences where used for the spatial derivatives as discretization method. The resulting time-dependent system was resolved using the ode23s solver in MATLAB R2012a (MathWorks, Natick, MA). Parameter optimization was carried out in the same software using the nlinfit function and confidence intervals were obtained with the nlparci function. This last function calculates the 95% confidence interval of optimized parameters using their value, residual value, and Jacobian matrix, all of which are given by the nlinfit function.

System pressure was calculated in each iteration using mean temperature and density (calculated dividing the total CO2 mass in the vessel by the effective volume) into the state equation proposed by Huang et al. [16]. Physical properties, such as ␳ and cP , of pure CO2 were obtained using the NIST Database [17] with local temperature and system pressure as inputs. All of these properties were updated in every iteration of the ode solver by introducing the updated local temperature and system pressure into the NIST function for MATLAB, refpropm. In the event of phase separation, the NIST database needs the user to specify if the properties given are those of a gas or a liquid. Because of this, physical properties when phase separation occurs where calculated as the weighted mean between their values in the liquid and gas state and the weight was the quality of the mixture. To explore the use of the model for practical purposes, we simulated different depressurization processes with varying parameters, like vessel volume and geometry, and initial pressure and temperature. We tested five different vessel volumes filled with CO2 and a packed bed made of 1-cm glass beads following a progression from laboratory scale, to pilot scale, and to industrial scale. Tested volumes were 0.005 m3 , 0.08 m3 , 0.25 m3 , 0.5 m3 and 1 m3 with an aspect ratio equal to 5 (i.e. the height-to-diameter ratio meaning that, for an aspect ratio of 5, if the vessel has a diameter of 1 m then its height will be 5 m) chosen by the authors based ˜ on a similar assumption made by Núnez et al. [18]. Valve opening for each case was chosen so that the minimum temperature in the system did not decrease below 0 ◦ C as a precaution so that the vessel material does not become brittle [7]. Four different aspect ratios (4, 4.5, 5, and 6) were also tested in a 1-m3 vessel using the same valve opening found in the previous step. The same was done with different initial temperature and pressure, where tested values were 60, 70, and 80 ◦ C at 30 MPa, and 30, 50, and 70 MPa at 60 ◦ C, respectively. These values were chosen because they are plausible temperatures of the thermal fluid in the heating jacket and extraction pressures, respectively. All of these simulations were run until the pressure inside the vessel reached near-atmospheric values. The thickness of the steel wall and lids in these tests were calculated using standard formulae in literature for pressurized vessel design [19]. The number of nodes was not increased for these tests due to the lack of experimental results to compare the simulations with, so we chose to keep the time needed to run each test to a minimum in order to be able to test the effect of additional parameters.

3. Results The new optimized equation for convective heat transfer at the vessel wall is shown in Eq. (14). Confidence intervals for all nine parameters can be found in Table 2. Adjustment between simulated and experimental temperatures improved by 18% when using these new parameters; the mean absolute percentage error decreased from 8.1% to 6.6%, as seen in Fig. 1. Fig. 2 shows simulated temperatures (Fig. 2A), pressures (Fig. 2B), and mass venting rate (Fig. 2C)

M.S. Murias et al. / J. of Supercritical Fluids 122 (2017) 43–51

47

Table 2 Parameters and confidence intervals for the regression Nu = aDab Rac between Nusselt dimensionless number (Nu), and Rayleigh (Ra) and Darcy (Da) dimensionless numbers. Parameter

Lower value

Estimated value

Upper value

a b c

0.0772 −0.0374 0.3940

0.0777 −0.0373 0.3967

0.0783 −0.0372 0.3990

using Eq. (14). Adjustment is also very satisfactory for pressure drop and for mass venting rate. Nu = 0.0777Da−0.0373 Ra0.397

(14)

Fig. 3 shows the simulated temperatures for all points of the extraction vessel at the start of the process (1 min), after increasing the valve opening, and at the middle of the depressurization (12 min). Although temperature changes throughout the process, the distribution of the temperature along the vessel does not drastically change. This is not the case with the temperature distribution in the solid (Fig. 3A–C). Simulations with different vessel volume and geometry, and initial pressure and temperature were run until the vessel pressure neared atmospheric values. However, in practical applications it is preferable to stop before this point as the time needed to decrease 1 MPa rises considerably after saturation conditions are met. Because of this, Table 3 shows the time needed for the vessel to reach 1 MPa for all tests. Note that this time is an overestimation of the actual depressurization time in an industrial plant, because it is not practical to work with a constant valve opening area. In practical applications, workers will further open the valve if dry ice is not observed. As referential values, Fiori [3] suggests that the total time for reconditioning should be around one hour for optimal operation of a 0.81-m3 extraction vessel operating at 60 ◦ C and 55 MPa. Quirin & Gerard [20] estimate this time to be around 45 to 90 min in an industrial plant depending on the size of the extractor. Fig. 4 shows results of tests with different volumes. Simulation of a 0.5-m3 extraction vessel are not shown because they were virtually the same as simulations for a 1-m3 vessel. Total depressurization times went from 8 to 19 min. The valve opening used in each case was chosen so that temperature reached 0 ◦ C at the lowest. These openings are shown as a function on vessel volume in Fig. 5. The correlation log(A) = 1.04 log(V ) − 4.81 was found between the valve opening chosen and vessel volume with a R2 coefficient of 0.97. Fig. 6 shows the results of tests for a 1-m3 extraction vessel with different aspect ratios. Even though the largest ratio resulted in the shortest depressurization time for the same opening, the temperature recovery was best when the aspect ratio was smallest. This last thing was not expected because of the shorter distance between the vessel wall and its center at larger aspect ratios. Pressure drop seems to be relatively unaffected by aspect ratio through the depressurization, except in the last increase. Fig. 7 shows the results of tests for a 1-m3 extraction vessel with different initial pressures. The pressure at which the system reached saturation conditions decreases with higher initial pressures, as observed also by Gebbeken & Eggers [5]. At higher initial pressure, total depressurization time with the same valve opening was higher, as expected. Temperature recovery was better at lower pressures but the minimum temperature reached was not greatly affected by the change in initial pressure. Results of depressurizations for a 1-m3 extraction vessel with different initial temperatures are shown in Fig. 8 and did not turn out as expected. While depressurizations at 70 ◦ C and 80 ◦ C follow similar patterns, the one starting at 60 ◦ C has a less steep temper-

Fig. 2. Comparison between simulated temperature (A), pressure (B) and mass venting rate (C), and experimental results [9] at the center of the vessel for depressurizations starting at 60 ◦ C and 26 MPa packed with sintered steel cylinders or glass beads. Temperature, pressure, and mass flow are simulated in an accurate way in the first part of the depressurization process but later deviate from the changes occurring in the experimental data.

ature drop and even reaches higher temperatures than the other ones. 4. Discussion The developed model was used to simulate depressurizations of extraction vessels as large as 1 m3 , with different vessel geometries,

48

M.S. Murias et al. / J. of Supercritical Fluids 122 (2017) 43–51

A

D 60 55 50 45 40 35 30

B

E 60

50 45 40 35

Temperature (ºC)

55

30

F

C

60 55 50 45 40 35 30 1.0

0.8

Dime

0.6

0.4

nsion

0.2

les s a

0.0

0.0

xial p

ositio

n (–)

Dim

ens

ss r ionle

0.2

adia

l

0.4

0.6

itio pos

0.8

n (–

)

Fig. 3. Simulated solid (A–C) and fluid (D–F) temperatures in the vessel for depressurizations with sintered steel cylinders starting at 60 ◦ C and 30 MPa after 1 min (A, D), after increasing the valve opening (B, E), and after 12 min of starting the process (C, F). The dimensionless radial position goes from zero at the vessel axis to one at its lateral wall, whereas the dimensionless axial position goes from zero at the bottom of the vessel to one at its top.

and with different starting temperatures and pressures. Simulated depressurization times went from 4 to 55 min, which is in accordance to the reference values reported previously [3,20]. However,

some of these times were too short (fast depressurizations), reaching undesirable lower temperatures (less than −10 ◦ C). Simulated times presented in this work are overestimated respect of the times

M.S. Murias et al. / J. of Supercritical Fluids 122 (2017) 43–51

49

Table 3 Summary table for simulated time required for an extraction vessel to decrease its pressure to 1 MPa at different vessel geometries, starting conditions, and constant valve opening area. Vessel volume (L)

Aspect ratio (L/D,–)

Initial temperature (◦ C)

Initial pressure (MPa)

Valve opening area (mm2 )

Time to reach 1 MPa (min)

5 80 250 500 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000

5.0 5.0 5.0 5.0 5.0 4.0 4.5 5.0 6.0 5.0 5.0 5.0 5.0 5.0 5.0

60 60 60 60 60 60 60 60 60 60 60 60 60 70 80

30 30 30 30 30 30 30 30 30 30 50 70 30 30 30

0.08 0.60 5.00 11.00 16.00 16.00 16.00 16.00 16.00 16.00 16.00 16.00 16.00 16.00 16.00

5.00 11.17 11.67 13.50 12.83 13.83 13.16 12.83 10.83 12.83 26.83 54.50 12.83 3.83 7.33

A

-5

Temperature (ºC)

2

Valve opening area (m )

10

80 L

5L 1000 L

B

30

-7

10

-2

10

25

Pressure (MPa)

-6

10

-1

10

0

10

Extraction vessel volume (m3)

20

Fig. 5. Valve opening area required for minimum temperature not to decrease below 0 ◦ C (symbols) in depressurizations starting at 60 ◦ C and 30 MPa, and correlation of that opening as a function of vessel volume (solid line).

15

10

250 L

5

0 0

2

4

6

8

10

12

14

16

18

Time (min) Fig. 4. Temperature and pressure changes for depressurizations starting at 60 ◦ C and 30 MPa with different volume of vessels packed with glass beads.

required to depressurize an extraction vessel in an industrial plant at conditions studied where the valve opening is not kept at a constant value, but opened and closed according to what is needed in the operation of the plant. So, our ultimate goal is applying the model to obtain an optimal mass flow of CO2 leaving the extraction vessel by varying the valve opening in a control loop that incorporates this information. Temperature and pressure behavior presented the following five stages during the depressurization of extraction vessels of different volumes: (i) an abrupt decrease of the variable; then, (ii) a plateau in which temperature and pressure became stable; (iii) another decrease; (iv) a rise and fall of temperature (and pressure in larger

vessels); and finally, (v) a steady rise in temperature, and fall in pressure. This differs from the experimental observations of Richter et al. [8] in a small 1-gallon vessel whose experiments only showed stages (i), (ii), and (v) because their attempts to avoid a gas-liquid transition. This implies that stages (iii) and (iv) are due to gas-liquid transitions phenomena within the extraction vessel during faster depressurization. The plateau reached at stage (ii) could be related to CO2 evaporation occurring within the vessel, which would be in accordance with the results of Gebbeken & Eggers [5], who identified phase separation after the initial drop of pressure was over. Local density plays a large role in these changes. In smaller vessels (0.005 and 0.080 m3 ), stage (ii) coincides with density increases in the vessel core from 600 to 900 kg m−3 , which might indicate liquid formation in there. Stage (iii) is initiated by the decrease in density from values around 900 to nearly 100 kg m−3 , which could mean that all the previously mentioned liquid has been evacuated. In larger vessels, density rapidly decreases from 800 to 50 kg m−3 , which explains the longer plateau in temperature (stage ii), indicating the onset of evaporation. The increase in pressure in that stage is also explained, as gas takes up more space inside the vessel. Depressurizations with different initial conditions showed interesting results. In the case of depressurizations of extraction vessels at different starting temperatures, results did not go as

50

M.S. Murias et al. / J. of Supercritical Fluids 122 (2017) 43–51

A

60

Temperature (ºC)

50

40

30

30 MPa

20

50 MPa

10

70 MPa

0

-10

B

14

Pressure (MPa)

12 10 8 6 4 2 0 0

Fig. 6. Temperature at the center of the vessel and pressure changes for depressurizations starting at 60 ◦ C and 30 MPa in 1-m3 vessels of different aspect ratio (L/D) packed with glass beads. Higher aspect ratio appears to decrease depressurization time but the vessel reaches lower temperatures, which could potentially lead to the vessel becoming brittle.

expected. The simulated temperature in the depressurization starting at 60 ◦ C differed from the other two. Even more, the final temperature was highest in the depressurization starting at 60 ◦ C, and lowest in the depressurization starting at 70 ◦ C. However, this can be explained by the properties of CO2 during these depressurizations. In the depressurizations starting at 70 and 80 ◦ C, CO2 went from a supercritical state (density around 800 kg m−3 ), to a subcooled liquid state (density around 1100 kg m−3 ), and then to a superheated gas state (density around 50 kg m−3 ) during the process. In the case of depressurizations starting at 60 ◦ C, CO2 transitions almost instantly from a supercritical state to a superheated gas state (density rapidly goes from 800 to 40 kg m−3 ), which allows it to reach stage (ii) (the plateau of temperature) faster than in the other cases. In other words, lower initial temperatures reduce the magnitude of the initial drop by allowing early CO2 evaporation. In depressurizations with different starting pressures, the ones with higher starting pressure and the same valve opening area were expected to last longer. However, that seemed to be the only visible effect, as the variations of temperature followed the same patterns in each case. Pressures at the end of stage (i) were lower in depressurizations at higher starting pressures, which was also seen by Gebbeken & Eggers [5] and attributed to higher initial specific entropy. At higher initial pressures, the pressure needed for CO2 to transition from a superheated gas state to a subcooled liquid state lowers. This means that it will take longer to reach that pressure inside the vessel and, because of that, temperature reached lower values as initial pressure increased.

6

12

18

24

30

36

42

48

54

60

Time (min) Fig. 7. Temperature at the center of the vessel and pressure changes for depressurizations in a 1-m3 vessel packed with glass beads starting at 60 ◦ C and with different initial pressures. Lower initial pressure seems to shorten depressurization time and improve temperature recovery, but the minimum temperature reached is not greatly affected.

In depressurization of extraction vessels with different aspect ratios, it can be seen that, initially, both temperature and pressure seem independent of this factor. This means that even reducing the distance between the wall and the packed bed core does not help compensate the energy loss due to mass flow. However, the temperature at the end of the depressurization is higher when the aspect ratio is smallest. Even though the largest aspect ratio had the shortest depressurization time, the lower final temperature makes further analysis necessary in order to make a recommendation as to which value is best for optimization. If temperature is still close to the minimum temperature allowed, it would not be possible to further open the valve. This means that a smaller aspect ratio could shorten depressurization times by allowing CO2 to leave the vessel faster. The simulations described previously were possible after we determined a correlation between Nu as a function of Ra and Da to estimate the convective heat transfer at the vessel wall. The coefficients of this correlation are similar to those proposed by Richter et al. [9] but predictions were improved by almost 20% when using the parameters obtained through simulation. Confidence intervals support the relevancy of each parameter in the model, even though parameter b in Eq. (14) (the exponent of Da) is in the order of magnitude of 10−2 . This means that Da is appropriate to account for the effects of the geometry of the vessel and the properties of the packed bed on convective heat transfer at the vessel wall.Temperature, pressure, and mass changes during the depressurization experiments carried out by Richter et al. [9] were simulated.

M.S. Murias et al. / J. of Supercritical Fluids 122 (2017) 43–51

A

Temperature (ºC)

80

60

40

60 ºC

20

51

to give first insight of depressurization times at different volumes, geometries, and initial conditions. Phase separation was roughly considered in this model but further investigation is needed in order to properly simulate it, and industrial results are necessary in order to validate simulations. After these modifications, this model could be used to create a process control capable of adjusting valve openings at each instant so that depressurizations can be optimized. Acknowledgments

80 ºC

0

This research was funded by FONDECYT projects1120827 and 1150623. M. Soledad Murias would like to thank Chilean agency Conicyt for a scholarship to carry out MSc. studies in Engineering Science in UC (CONICYT-PCHA/Magíster Nacional/2015 – folio 22151040). We would also like to thank Simula UC, for lending us a computer powerful enough to run our simulations, and Prof. Claudio Gelmi, for his insight on Matlab and heat transfer.

70 ºC

-20

B

30

25

Pressure (MPa)

References 20

15

10

5

0 0

2

4

6

8

10

12

14

16

Time (min) Fig. 8. Temperature at the center of the vessel and pressure changes for depressurizations in a 1-m3 vessel packed with glass beads starting at 30 MPa and with different initial temperatures. Lower initial temperature seems to lengthen the depressurization process but the lower temperature reached is higher than in the other cases, which could lead to further opening of the valve and shortening the process time.

Temperature distributions in Fig. 3 help illustrate the role played by the packed bed on heat transfer. At times where fluid temperature drop is higher, solid temperature drops in order to act as a buffer for the fluid. On the contrary, when temperature is stabilizing, solid temperature is almost homogeneous. These results confirm previous belief [8,9] that the solid phase acts as a reservoir of energy for the fluid phase and that accuracy in the calculations of physical properties of the solid phase is very relevant for depressurizations. Initial temperature drop when opening the valve, which was done twice during experiments, is more abrupt in simulated than in experimental results. This could be because the valve opening area is instantly changed in the model but, in reality, opening area takes multiple values in the time it takes for it to change from its initial to its final setting. Also, thermal inertia of the vessel could have been overestimated, which would mean that the heat retained by the vessel walls and by the packed material is not being released into the system as fast as in the model. This difference seems to have a greater effect on pressure drop after the second stage of valve opening than in the first stage. One reason for this could be that pressure was calculated using the mean temperature of the vessel, which has a small deviation at the beginning of the depressurization but not after the second stage. In conclusion, the model here proposed simulated depressurizations reported in previous work reasonably well, and was able

[1] J.M. del Valle, Extraction of natural compounds using supercritical CO2 : going from the laboratory to the industrial application, J. Supercrit. Fluids 96 (2015) 180–199. ˜ [2] G.A. Núnez, J.M. del Valle, Supercritical CO2 oilseed extraction in multi-vessel plants. 2. Effect of number and geometry of extractors on production cost, J. Supercrit. Fluids 92 (2014) 324–334. [3] L. Fiori, Supercritical extraction of grape seed oil at industrial-scale: plant and process design modeling, economic feasibility, Chem. Eng. Process. 49 (2010) 866–872. [4] R. Eggers, V. Green, Pressure discharge from a pressure vessel filled with CO2 , J. Loss Prev. Process Ind. 3 (1990) 59–63. [5] B. Gebbeken, R. Eggers, Blowdown of carbon dioxide from initially supercritical conditions, J. Loss Prev. Process Ind. 9 (1996) 285–293. [6] J. Zhang, D. Zhu, W. Tian, S. Qiu, G. Su, D. Zhang, Depressurization study of supercritical fluid blowdown from simple vessel, Ann. Nucl. Energy 66 (2014) 94–103. [7] R.E. Smallman, A.H.W. Ngan, Physical Metallurgy and Advanced Materials, 7th edition, Butterworth-Heinemann, Oxford, 2007. ˜ [8] E.A. Richter, J.M. del Valle, G.A. Núnez, Thermodynamic properties of CO2 during controlled decompression of supercritical extraction vessels, J. Supercrit. Fluids 98 (2015) 102–110. [9] E.A. Richter, M.S. Murias, J.M. del Valle, Heat transfer and venting rate during controlled decompression of supercritical extraction vessels, J. Supercrit. Fluids 120 (2017) 275–284. [10] T.L. Bergman, A.S. Lavine, F.P. Incropera, D.P. Dewitt, Fundamentals of Heat and Mass Transfer, 7th edition, John Wiley & Sons, Inc., Danvers, MA, 2011. [11] K. Slimi, S. Ben Nasrallah, J.P. Fohr, Transient natural convection in a vertical cylinder opened at the extremities and filled with a fluid saturated porous medium: validity of Darcy flow model and thermal boundary layer approximations, Int. J. Heat Mass Transf. 41 (1997) 1113–1125. [12] A.G. Dixon, Correlations for wall and particle shape effects on fixed bed, Can. J. Chem. Eng. 66 (1988) 705–708. [13] N. Wakao, S. Kaguei, T. Funazkri, Effect of fluid dispersion coefficients on particle-to-fluid heat transfer coefficients in packed beds, Chem. Eng. Sci. 34 (1978) 325–336. [14] A.G. Dixon, D.L. Cresswell, Theoretical prediction of effective heat transfer parameters in packed beds, AIChE J. 25 (1979) 663–676. [15] E.U. Schlünder, E. Tsotsas, Wärmeübertragung in Festbetten, durchmischten Schüttgütern und Wirbelschichten, Thieme Verlag, Stuttgart-New York, 1988. [16] F. Huang, M. Li, L.L. Lee, K.E. Starling, F.T.H. Chung, An accurate equation of state for carbon dioxide, J. Chem. Eng. Japan 18 (1984) 490–496. [17] E.W. Lemmon, M.L. Huber, M.O. McLinden, NIST Standard Reference Database 23: Mini-reference Fluid Thermodynamic and Transport Properties (REFPROP), Version 9.0, National Institute of Standards and Technology, Standard Reference Data Program, Gaithersburg, ML, 2007. ˜ [18] G.A. Núnez, C.A. Gelmi, J.M. del Valle, Simulation of a supercritical carbon dioxide extraction plant with three extraction vessels, Comput. Chem. Eng. 35 (2011) 2687–2695. [19] H.H. Bednat, Pressure Vessel Design Handbook, 2nd edition, Krieger Publ. Co., Malabar, FL, 1996. [20] K.W. Quirin, D. Gerard, Supercritical fluid extraction (SFE), in: Flavourings. Production, Composition, Applications, Regulations, 2nd. edition, John Wiley & Sons, Weinheim, Germany, 2007, pp. 44–65.