Numerical Analysis of Dry Ice Blasting Convergent-Divergent Supersonic Nozzle

There are several well-known and widely used industrial cleaning methods in the market today. One of them is dry ice blasting. In this method, moisture-free air is compressed, mixed with solid CO2 particles, and blasted though a nozzle; in the process, the gas expands, propelling its velocity. The high-speed, two-phase flow cleans by supercooling and crushing particles on the surface, causing dry ice sublimation. As the nozzle is a crucial component of the system, the authors conducted a numerical analysis of the geometry of the proposed convergent-divergent nozzle. A mathematical model of the supersonic, two-phase flow was developed and implemented in commercial Computational Fluid Dynamics (CFD) code. Various operating parameters, such as inlet pressure and dry ice mass flow, were taken into consideration.

Industrial pollution is a huge challenge for plants operators. Impurities can cause a decrease in device efficiency, improper operation, or damage [1]. Various pollutant removal methods have been designed and successfully operated. One of them is dry ice blasting, which is characterized by moisture-less operating conditions [2], significantly lower amount of waste after cleaning (in comparison to i.e., sand blasting), and easy control of operational parameters by mass flow and pressure change. It is a relatively cheap and safe method [3]; however, the cooling mechanism of the cleaning surface is one of its limitations [4].
In dry ice blasting, the first step is air compression and drying. Dry ice in the form of pellets (3 mm in diameter and 5-15 mm in length) [5] is then loaded into a shredder. If pellets of smaller dimensions are required, they are run through a strainer. Dry ice is produced in a separate system [6], i.e., using the process described in [7] or on site [8] by conversion of liquid CO 2 to pellets with the use of a pelletizer. Once the pellets are loaded in the two-phase dry ice blasting machine, compressed air blasts them through a nozzle on to the cleaning surface [9]. The flow parameters are highly dependent on air pressure, dry ice mass flow rate, and nozzle-surface distance. The schematics of the system is shown in Figure 1.
The dry ice cleaning mechanism is based on three main phenomena: thermal effect (cooling), abrasion using flow kinetic energy, and sublimation [10]. The final force impacting the cleaning surface is a sum of three components: force of compressed air, force of solid CO 2 particles caused by their velocity, and sublimation force caused by sudden phase change supported by rapid velocity growth [11]. Dry ice blasting can be used for cleaning of electrical devices in the automobile, aircraft, and railway industries [12]. Although dry ice blasting was developed for industrial cleaning purposes, it can be successfully used as a surface pre-treatment method i.e., to increase the adhesive strength of aluminum bonding joints [13], in aluminum oxide coatings [14] or titanium coatings [15]. The dry ice cleaning mechanism is based on three main phenomena: thermal effect (cooling), abrasion using flow kinetic energy, and sublimation [10]. The final force impacting the cleaning surface is a sum of three components: force of compressed air, force of solid CO2 particles caused by their velocity, and sublimation force caused by sudden phase change supported by rapid velocity growth [11]. Dry ice blasting can be used for cleaning of electrical devices in the automobile, aircraft, and railway industries [12]. Although dry ice blasting was developed for industrial cleaning purposes, it can be successfully used as a surface pre-treatment method i.e., to increase the adhesive strength of aluminum bonding joints [13], in aluminum oxide coatings [14] or titanium coatings [15].
A key component of the dry ice blasting system is the nozzle. Its geometry has a significant influence on the two-phase flow outlet parameters and, consequently, on cleaning speed. In this study, a mathematical model of nozzle flow via a convergent-divergent channel was built and implemented in a numerical environment. The phenomena to be modeled was a supersonic twophase flow. Very few researchers have attempted dry ice blasting nozzle optimization from a cleaning efficiency point of view. There are papers describing the two-phase flow via a convergencedivergence nozzle. In [16], a model of air-water mixture via the Laval nozzle is presented. In [17] and [18], nozzles made specifically for dry ice blasting were examined; the authors aimed to optimize the nozzle shape to decrease operational noise. A model of a dry ice blasting nozzle is also presented in [19], where the researchers focused on a plasma spraying application.

General Equations
The numerical model for the continuous phase was based on the following equations (i = 1,2,3): Energy: Momentum: where viscous stress tensor can be written as: Continuity: State: A key component of the dry ice blasting system is the nozzle. Its geometry has a significant influence on the two-phase flow outlet parameters and, consequently, on cleaning speed. In this study, a mathematical model of nozzle flow via a convergent-divergent channel was built and implemented in a numerical environment. The phenomena to be modeled was a supersonic two-phase flow. Very few researchers have attempted dry ice blasting nozzle optimization from a cleaning efficiency point of view. There are papers describing the two-phase flow via a convergence-divergence nozzle. In [16], a model of air-water mixture via the Laval nozzle is presented. In [17] and [18], nozzles made specifically for dry ice blasting were examined; the authors aimed to optimize the nozzle shape to decrease operational noise. A model of a dry ice blasting nozzle is also presented in [19], where the researchers focused on a plasma spraying application.

General Equations
The numerical model for the continuous phase was based on the following equations (i = 1,2,3): Energy: Momentum: where viscous stress tensor can be written as: Continuity: State: Effective fluid viscosity can be expressed as a sum of molecular viscosity and turbulent eddy viscosity: The total energy heat transfer model was created by identifying high-speed flow value (Mach number > 0.3) and compressibility effects of the system.

Turbulence Model
A variety of turbulence models have been developed and tested for a wide range of flow applications. Among the most popular for high-speed flow modeling are the k-ε and k-ω (both with some modifications) models [20,21]. The standard k-ε model can be easily applied to describe most turbulent flows. Its advantages include reasonable accuracy for a wide variety of flows and robustness. A main drawback, however, is its poor performance in non-equilibrium boundary layers [22]. k-ω is based on a similar approach for turbulence modeling like k-ε, but turbulent energy dissipation rate (ω) is used instead of dissipation equation. In contrast to k-ε, the k-ω model can be successfully applied for near-wall region flows calculations, even when high pressure gradients occur. Its main disadvantage is strong sensitivity to values of ω in the freestream far from the boundary layer [23]. The Shear Stress Transport (SST) version of the k-ω model helps overcome this drawback by application of the k-ε model in free stream regions [24]. That is why the k-ω SST model was applied in these simulations.

Discrete Phase Model
The dry ice particles are assumed to be spherical. Their particle diameters are set as fulfilling the Rosin-Rammler distribution with the given Rosin-Rammler power and mean diameter. The drag function is defined as [25]: where the drag coefficient is calculated using the Schiller-Naumann equation [26]:

Geometry
Two geometries were analyzed in this paper. Geometry A is the nozzle used in the dry ice blasting application. As nozzle length is significant-can cause operational difficulties in some applications and its shape seems to be complicated-a new solution was proposed (geometry B). The analyzed geometries consist of the convergent-divergent nozzle canal, directly connected to a box where flow can be observed. As the domain has symmetrical planes, a quarter of the whole geometry was considered. In both cases, the nozzle inlet cross-section area (A in ) equals 2 cm 2

Computational Fluid Dynamics (CFD) Grid
As numerous simulations were planned, a grid test was done to obtain proper results at reasonable computational costs. The test results are presented for Geometry B. Four grids were taken into account, three of them generated using an automatic method, with face sizing on the nozzle wall. In the last one, an additional layering on the nozzle wall was implemented. The grid's details are presented in Table 1. Comparing the Y+ values on the nozzle's wall (Figure 4), grid III was selected for further simulations. An analogous process was conducted for Geometry A. One mesh from each geometry was selected for further processing ( Figures 5 and 6). Main mesh parameters are presented in Table 2.

Computational Fluid Dynamics (CFD) Grid
As numerous simulations were planned, a grid test was done to obtain proper results at reasonable computational costs. The test results are presented for Geometry B. Four grids were taken into account, three of them generated using an automatic method, with face sizing on the nozzle wall. In the last one, an additional layering on the nozzle wall was implemented. The grid's details are presented in Table 1.
Comparing the Y+ values on the nozzle's wall (Figure 4), grid III was selected for further simulations. An analogous process was conducted for Geometry A. One mesh from each geometry was selected for further processing (Figures 5 and 6). Main mesh parameters are presented in Table 2.

Computational Fluid Dynamics (CFD) Grid
As numerous simulations were planned, a grid test was done to obtain proper results at reasonable computational costs. The test results are presented for Geometry B. Four grids were taken into account, three of them generated using an automatic method, with face sizing on the nozzle wall. In the last one, an additional layering on the nozzle wall was implemented. The grid's details are presented in Table 1. Comparing the Y+ values on the nozzle's wall (Figure 4), grid III was selected for further simulations. An analogous process was conducted for Geometry A. One mesh from each geometry was selected for further processing (Figures 5 and 6). Main mesh parameters are presented in Table 2.

Boundary Conditions and General Settings
The calculations were conducted using the Ansys CFX commercial software. A flow regime at the inlet was set as mixed, normal speed (63,053 m/s), and relative static pressure was introduced to the model. High turbulence (10% of intensity) was chosen. Static temperature at the level of 281.17 was applied. Reference pressure for the domain was set as 0 Pa, which is the recommended setting for supersonic flows. Opening boundary condition was selected at the outlet, with the relative pressure =1 bar. Medium intensity turbulence and static temperature (281.17 K) were also selected for this surface. Dry ice was modeled as a solid with molar mass equal to 44.01 kg/kmol and density of 1.55 g/cm 3 [11]. Rosin-Rammler particle size at the level of 10 −5 m and Rosin-Rammler power of 3.5 were introduced to the model. Analysis was done in the steady state mode. Convergence criteria was set at the level of 10 −5 for the RMS residual type. For an advection scheme and turbulence numerics, a high-resolution option was selected, which means that the blending factor between the first and second order upwind schemes, used to calculate the advection terms, varies in the domain, depending on the local solution field used to enforce a boundedness criterion. For low variable gradients, the blend factor takes the value of 1.0 (second order) to provide high accuracy. In the area of rapid parameter change, the blend factor decreases to 0.0 (first order) to avoid overshoots and undershoots and maintain robustness. Auto timescale was chosen for fluid timescale control. The conservative length scale option was set. Boundary conditions implemented in the geometry are shown in Figure 7 and presented in Table 3. The calculations were done on two processors (Intel ® Xeon ® CPU E5-2600 2.20 GHz) with 64 GB RAM memory using eight computing nodes.

Boundary Conditions and General Settings
The calculations were conducted using the Ansys CFX commercial software. A flow regime at the inlet was set as mixed, normal speed (63,053 m/s), and relative static pressure was introduced to the model. High turbulence (10% of intensity) was chosen. Static temperature at the level of 281.17 was applied. Reference pressure for the domain was set as 0 Pa, which is the recommended setting for supersonic flows. Opening boundary condition was selected at the outlet, with the relative pressure =1 bar. Medium intensity turbulence and static temperature (281.17 K) were also selected for this surface. Dry ice was modeled as a solid with molar mass equal to 44.01 kg/kmol and density of 1.55 g/cm 3 [11]. Rosin-Rammler particle size at the level of 10 −5 m and Rosin-Rammler power of 3.5 were introduced to the model. Analysis was done in the steady state mode. Convergence criteria was set at the level of 10 −5 for the RMS residual type. For an advection scheme and turbulence numerics, a high-resolution option was selected, which means that the blending factor between the first and second order upwind schemes, used to calculate the advection terms, varies in the domain, depending on the local solution field used to enforce a boundedness criterion. For low variable gradients, the blend factor takes the value of 1.0 (second order) to provide high accuracy. In the area of rapid parameter change, the blend factor decreases to 0.0 (first order) to avoid overshoots and undershoots and maintain robustness. Auto timescale was chosen for fluid timescale control. The conservative length scale option was set. Boundary conditions implemented in the geometry are shown in Figure 7 and presented in Table 3. The calculations were done on two processors (Intel ® Xeon ® CPU E5-2600 2.20 GHz) with 64 GB RAM memory using eight computing nodes.

Model Validation
Indirect model validation was conducted due to difficulties in the direct measurement of flow parameters. In the experiment, the nozzle in geometry A was used ( Figure 8). As described in [27], outlet flow velocity has a crucial influence on dry ice blasting cleaning time. It can be assumed that cleaning speed depends on the kinetic energy of the stream; therefore, ~u 2 . During the experiment, the given surface was cleaned from two distances: 15 cm and 30 cm, and cleaning time was measured. The first test involved a steel sheet covered in synthetic paint ( Figure 9): inlet pressure of 5 bar, dry ice mass flow rate of 48 kg/h (test I). In the second test, a ceramic surface was covered in flour with wallpaper glue (Figure 10): inlet pressure of 3 bar, dry ice mass flow rate of 40 kg/h (test II). The validation results are presented in Table 4.

Model Validation
Indirect model validation was conducted due to difficulties in the direct measurement of flow parameters. In the experiment, the nozzle in geometry A was used ( Figure 8). As described in [27], outlet flow velocity has a crucial influence on dry ice blasting cleaning time. It can be assumed that cleaning speed depends on the kinetic energy of the stream; therefore,~u 2 . During the experiment, the given surface was cleaned from two distances: 15 cm and 30 cm, and cleaning time was measured. The first test involved a steel sheet covered in synthetic paint ( Figure 9): inlet pressure of 5 bar, dry ice mass flow rate of 48 kg/h (test I). In the second test, a ceramic surface was covered in flour with wallpaper glue (Figure 10): inlet pressure of 3 bar, dry ice mass flow rate of 40 kg/h (test II). The validation results are presented in Table 4.

Model Validation
Indirect model validation was conducted due to difficulties in the direct measurement of flow parameters. In the experiment, the nozzle in geometry A was used ( Figure 8). As described in [27], outlet flow velocity has a crucial influence on dry ice blasting cleaning time. It can be assumed that cleaning speed depends on the kinetic energy of the stream; therefore, ~u 2 . During the experiment, the given surface was cleaned from two distances: 15 cm and 30 cm, and cleaning time was measured. The first test involved a steel sheet covered in synthetic paint ( Figure 9): inlet pressure of 5 bar, dry ice mass flow rate of 48 kg/h (test I). In the second test, a ceramic surface was covered in flour with wallpaper glue (Figure 10): inlet pressure of 3 bar, dry ice mass flow rate of 40 kg/h (test II). The validation results are presented in Table 4.

Surface
Boundary Condition Inlet Normal speed and relative static pressure, injection of surface particles Nozzle wall Adiabatic wall Sym1, Sym2 Symmetry Outlet Opening

Model Validation
Indirect model validation was conducted due to difficulties in the direct measurement of flow parameters. In the experiment, the nozzle in geometry A was used ( Figure 8). As described in [27], outlet flow velocity has a crucial influence on dry ice blasting cleaning time. It can be assumed that cleaning speed depends on the kinetic energy of the stream; therefore, ~u 2 . During the experiment, the given surface was cleaned from two distances: 15 cm and 30 cm, and cleaning time was measured. The first test involved a steel sheet covered in synthetic paint ( Figure 9): inlet pressure of 5 bar, dry ice mass flow rate of 48 kg/h (test I). In the second test, a ceramic surface was covered in flour with wallpaper glue (Figure 10): inlet pressure of 3 bar, dry ice mass flow rate of 40 kg/h (test II). The validation results are presented in Table 4.

Results and Analysis
Five cases for each geometry were considered to determine the impact of operational parameters on the outlet particles' velocity, which seems to be a crucial variable from a cleaning efficiency point of view. Details of the simulations are presented in Table 5. Infinite static pressure is the hypothetic pressure measured in the infinite tank connected to the nozzle inlet. Inlet relative static pressure is the pressure in the nozzle inlet, recalculated using formulas describing the flow via a convergent-divergent nozzle [28]. Examples of particle tracks and velocities are presented in Figures 11 and 12. The example of vector plot is presented in the Figure 13.

Results and Analysis
Five cases for each geometry were considered to determine the impact of operational parameters on the outlet particles' velocity, which seems to be a crucial variable from a cleaning efficiency point of view. Details of the simulations are presented in Table 5. Infinite static pressure is the hypothetic pressure measured in the infinite tank connected to the nozzle inlet. Inlet relative static pressure is the pressure in the nozzle inlet, recalculated using formulas describing the flow via a convergentdivergent nozzle [28]. Examples of particle tracks and velocities are presented in Figures 11 and 12. The example of vector plot is presented in the Figure 13.

Inlet Pressure Impact
The impact of the first inlet pressure was taken into consideration. In each analyzed case, a supersonic flow occurs. As Geometry A consists of two subsequent convergent-divergent parts, axial parameters distribution differs from values obtained for the simple convergent-divergent nozzle (Geometry B) (Figures 14-17). In the geometry A the pressure along the axis initially decrease, then reach the nozzle throat which is connected with rapid pressure drop, but this phenomena is predictable for the convergent-divergent nozzles. Further pressure growth and second decrease results from the nozzle shape. This geometry consists of two throats, what explains the functions shape (Pressure and Mach number along the axis). The flow via canal A reaches the maximum Mach number in the axis (equals to 1.5090) for the highest analyzed pressure (5 bar) ( Figure 14). For Geometry B, maximum Mach number in the nozzle axis varies from 1.207 for 3 bar, up to 1.679 for 5 bar ( Figure 16). As shown in Figure 15 and 17, in each case, a pressure drop occurs in the critical cross-section of the nozzle. Beyond this point, no significant discrepancies in total pressure is noticed in cases having the same geometry.
The lowest particle velocity was observed for inlet pressure of 3 bar. Particles speed up with the increase of initial pressure. A correlation between maximum particle velocity and initial total pressure is presented in Figures 18 and 19.

Inlet Pressure Impact
The impact of the first inlet pressure was taken into consideration. In each analyzed case, a supersonic flow occurs. As Geometry A consists of two subsequent convergent-divergent parts, axial parameters distribution differs from values obtained for the simple convergent-divergent nozzle (Geometry B) (Figures 14-17). In the geometry A the pressure along the axis initially decrease, then reach the nozzle throat which is connected with rapid pressure drop, but this phenomena is predictable for the convergent-divergent nozzles. Further pressure growth and second decrease results from the nozzle shape. This geometry consists of two throats, what explains the functions shape (Pressure and Mach number along the axis). The flow via canal A reaches the maximum Mach number in the axis (equals to 1.5090) for the highest analyzed pressure (5 bar) ( Figure 14). For Geometry B, maximum Mach number in the nozzle axis varies from 1.207 for 3 bar, up to 1.679 for 5 bar ( Figure 16). As shown in Figures 15 and 17, in each case, a pressure drop occurs in the critical cross-section of the nozzle. Beyond this point, no significant discrepancies in total pressure is noticed in cases having the same geometry.            Figure 17. Total pressure in the nozzle axis, Geometry B.
The lowest particle velocity was observed for inlet pressure of 3 bar. Particles speed up with the increase of initial pressure. A correlation between maximum particle velocity and initial total pressure is presented in Figures 18 and 19.

Dry Ice Mass Flow Rate Impact
In this section, the influence of dry ice mass flow rate on the flow parameters will be analyzed.

Dry Ice Mass Flow Rate Impact
In this section, the influence of dry ice mass flow rate on the flow parameters will be analyzed. Three cases were considered for each geometry: 30 kg/h, 40 kg/h, and 50 kg/h. They were calculated

Dry Ice Mass Flow Rate Impact
In this section, the influence of dry ice mass flow rate on the flow parameters will be analyzed. Three cases were considered for each geometry: 30 kg/h, 40 kg/h, and 50 kg/h. They were calculated assuming 4 bar for static inlet pressure. As can be observed from Figures 20 and 21 no significant discrepancies in Mach number occur between cases for the selected geometry. Dry ice mass flow rate has a negligible influence on the pressure in the axis (Figures 22 and 23). Maximum particle velocity decreases with growth in dry ice mass flow rate (Figures 24 and 25). The flow is little dependent on the dry-ice mass flow rate as its volumetric concentration in the air is very low but it should be mentioned that even such a small amount of the dispersed phase can influence on the sound velocity in the medium, what means that Mach number in each flow corresponds with a little bit different flow velocities expressed in m/s. This is the explanation for the maximum particles velocity change with the dry ice mass flow rate growth (Figures 24 and 25).

Dry Ice Mass Flow Rate Impact
In this section, the influence of dry ice mass flow rate on the flow parameters will be analyzed. Three cases were considered for each geometry: 30 kg/h, 40 kg/h, and 50 kg/h. They were calculated assuming 4 bar for static inlet pressure. As can be observed from Figures 20 and 21 no significant discrepancies in Mach number occur between cases for the selected geometry. Dry ice mass flow rate has a negligible influence on the pressure in the axis (Figures 22 and 23). Maximum particle velocity decreases with growth in dry ice mass flow rate (Figures 24 and 25). The flow is little dependent on the dry-ice mass flow rate as its volumetric concentration in the air is very low but it should be mentioned that even such a small amount of the dispersed phase can influence on the sound velocity in the medium, what means that Mach number in each flow corresponds with a little bit different flow velocities expressed in m/s. This is the explanation for the maximum particles velocity change with the dry ice mass flow rate growth (Figures 24-25).

Outcomes
Dry ice blasting is a promising cleaning method, even for sophisticated applications. The most important part of the system is the nozzle. In this study, two geometries of convergent-divergent nozzles were examined. Geometry A is the shape of the nozzle used during dry ice blasting. Geometry B was developed to simplify the canal shape and shorten it to allow operation in narrow spaces, which often occur in real applications. A mathematical model of the two-phase supersonic flow was built and implemented in a numerical environment. The simulations show the impact of the inlet pressure and dry ice mass flow rate on the particles and flow behavior in the nozzle. The results showed a high dependency of inlet pressure on flow and particle velocities; it was found to be lower for particle concentration. This outcome agrees with the experiment in [27], which proved that the crucial parameter, from a cleaning speed point of view, is inlet air pressure.
As shown in [29], sound velocity decreases with the growing number of particles in the mixture.

Outcomes
Dry ice blasting is a promising cleaning method, even for sophisticated applications. The most important part of the system is the nozzle. In this study, two geometries of convergent-divergent nozzles were examined. Geometry A is the shape of the nozzle used during dry ice blasting. Geometry B was developed to simplify the canal shape and shorten it to allow operation in narrow spaces, which often occur in real applications. A mathematical model of the two-phase supersonic flow was built and implemented in a numerical environment. The simulations show the impact of the inlet pressure and dry ice mass flow rate on the particles and flow behavior in the nozzle. The results showed a high dependency of inlet pressure on flow and particle velocities; it was found to be lower for particle concentration. This outcome agrees with the experiment in [27], which proved that the crucial parameter, from a cleaning speed point of view, is inlet air pressure.
As shown in [29], sound velocity decreases with the growing number of particles in the mixture. This implies that flow characterized by the same Mach number will move more slowly when containing more particles. This conclusion agrees with the simulation results.
The dry ice sublimation phenomenon was not included. Further development of the model will include dry ice phase change modeling. In this case, the ideal-gas equation should be replaced by real gas assumption.