Modelling of the Process of Extrusion of Dry Ice through a Single-Hole Die Using the Smoothed Particle Hydrodynamics (SPH) Method

This article presents the outcome of research on modelling the process of the extrusion of crystalline dry ice. The purpose of this process is to densify the material and obtain pellets of several millimeters in diameter. This reduces the sublimation rate in ambient conditions of the material whose temperature in a solid state is 195 K. A lower sublimation rate means a reduction of the loss of product in its final applications, which include refrigeration and reduction of atmospheric emissions of gaseous CO2. A ram-type extruder was considered in this analysis, in which dry ice was extruded through a single-hole die of varying geometry. The article presents the results of numerical analyses of the extrusion process, using a simulation method based on the Smoothed Particle Hydrodynamics (SPH) approach. The results from simulations were verified by the experimental data in terms of the maximum force required to complete the process, in order to assess the applicability of the proposed method in further research on dry ice compression.


Introduction
The issue of reuse of waste generated throughout the entire lifecycle of products is becoming a more important issue in the area of process engineering [1][2][3][4]. It is often practicable to recover waste materials generated during various production processes. Górecki observed that huge amounts of carbon dioxide are emitted into the atmosphere during production of ammonia compound, due to the high degree of purity of waste material from which it is recovered [5]. This material is condensed through compression and then stored at a pressure of 20 bar and a temperature of 216 K to 250 K [6,7]. Rapid expansion of atmospheric pressure results in a liquid-solid phase transition. In ambient conditions, crystallised carbon dioxide (CCD) has a temperature of 195 K and changes its state in the process of sublimation [8,9].
These peculiar properties are desired in a number of production processes, including refrigeration [10,11] and cleaning of surfaces [12,13]. These processes utilise the two specific characteristics of CCD, i.e., low temperature and sublimation. CCD obtained through expansion is a very fine grain material of ca. 100 µm size of particles [14,15]. The high specific surface characteristics of dry ice snow results in a high sublimation rate, affecting the process efficiency in refrigeration application [14]. For this reason, the most common commercial form of dry ice are pellets of 3-16 mm in diameter featuring a lower sublimation rate [16].
Dry ice snow is compressed in ram-based extrusion machines. The working system of such machines is shown schematically in Figure 1. The raw material is expanded inside the compression cavity (1). Then, crystalline carbon dioxide (3) is compressed by a ram pressing it against the residual, already compressed CCD (4) that fills up the die cavity (5).
The graph in Figure 2 represents the change of Fe as a function of ram displacement s, subdivided into three main phases of the change of force applied to extrude dry ice pellets.
In Phase 1, dry ice snow is compressed without any significant increase in the value of force Fe. This is because the displacement of ram reduces spaces between the particles. The exerted pressure does not induce elastic or plastic strains, due to ignorable contact between the particles. In Phase 2, the particles are packed closer together and the internal friction and elastic and plastic strain come into play as a result. The touching and relative displacement of particles increases the value of the exerted force Fe, due to elastic and plastic strain, and internal friction. In addition, the compressed material interacts with the compression cavity walls, due to elastic strain. This results in friction between the compressed material and the compression cavity walls, increasing the compression resistance. The dissipation of energy during compression, in cylindrical chambers with tapered ends, has been extensively described in the literature [18,19].
The value of force Fe increases until it equals the maximum resistance force of the process of compression FL. It has been demonstrated that the value of this force depends on the die parameters and the coefficient of friction between the extruded material and the side walls of the die [18]. The value of Fe decreases as the process of extrusion proceeds and the density of the material ρ no longer increases as a result [20]. Main parts of the ram-type extruder: 1-compression cavity, 2-ram, 3-die, 4-loose dry ice, 5-compressed dry ice [17].
The graph in Figure 2 represents the change of F e as a function of ram displacement s, subdivided into three main phases of the change of force applied to extrude dry ice pellets.
In Phase 1, dry ice snow is compressed without any significant increase in the value of force F e . This is because the displacement of ram reduces spaces between the particles. The exerted pressure does not induce elastic or plastic strains, due to ignorable contact between the particles. In Phase 2, the particles are packed closer together and the internal friction and elastic and plastic strain come into play as a result. The touching and relative displacement of particles increases the value of the exerted force F e , due to elastic and plastic strain, and internal friction. In addition, the compressed material interacts with the compression cavity walls, due to elastic strain. This results in friction between the compressed material and the compression cavity walls, increasing the compression resistance. The dissipation of energy during compression, in cylindrical chambers with tapered ends, has been extensively described in the literature [18,19].
The value of force F e increases until it equals the maximum resistance force of the process of compression F L . It has been demonstrated that the value of this force depends on the die parameters and the coefficient of friction between the extruded material and the side walls of the die [18]. The value of F e decreases as the process of extrusion proceeds and the density of the material ρ no longer increases as a result [20].
The graph in Figure 2 represents the change of Fe as a function of ram displacement s, subdivided into three main phases of the change of force applied to extrude dry ice pellets.
In Phase 1, dry ice snow is compressed without any significant increase in the value of force Fe. This is because the displacement of ram reduces spaces between the particles. The exerted pressure does not induce elastic or plastic strains, due to ignorable contact between the particles. In Phase 2, the particles are packed closer together and the internal friction and elastic and plastic strain come into play as a result. The touching and relative displacement of particles increases the value of the exerted force Fe, due to elastic and plastic strain, and internal friction. In addition, the compressed material interacts with the compression cavity walls, due to elastic strain. This results in friction between the compressed material and the compression cavity walls, increasing the compression resistance. The dissipation of energy during compression, in cylindrical chambers with tapered ends, has been extensively described in the literature [18,19].
The value of force Fe increases until it equals the maximum resistance force of the process of compression FL. It has been demonstrated that the value of this force depends on the die parameters and the coefficient of friction between the extruded material and the side walls of the die [18]. The value of Fe decreases as the process of extrusion proceeds and the density of the material ρ no longer increases as a result [20]. Compression of the material reduces the distances between the particles of the material observed as an increase in density ρ. As the value of ρ increases, so do the values of Young's modulus E, Poisson's ratio ν, and coefficient of friction µ [17,20,22]. The behaviour of the compressed material in Phase 2 may be accurately simulated with Drucker-Prager/Cap, Cam-Clay, and Mohr-Coulomb models [23][24][25].
The decrease of F e observed in Phase 3 marks the transition from compression to extrusion. The correlation coefficient defining the similarity of the F e vs. s curve, and a straight-line relationship is 0.9, which indicates a proportional decrease of force in relation to displacement. This is due to the linear change in the surface area of the side walls of the compression cavity (1 on Figure 1) that comes into contact with the material. This allows us to assume the constant values of E, ν, and µ, and use an elastic-plastic model to represent the process in numerical studies, as it has been done in previous studies reported in the literature [26][27][28].
According to Jankowiak [29] the Smoothed Particle Hydrodynamics (SPH) method can be used to accurately describe the geometry of materials featuring high plastic strain during compression. Originally, SPH was developed to simulate the flow of fluids. Since then, it was adapted to simulate processes that involve high strains, such as compression and extrusion of solids [16,29]. Sakai observed that the accuracy of predictions with the SPH method decreases in areas close to the boundary conditions, or free surfaces [30]. Even though the predictions obtained with SPH are close to the results obtained with the conventional mesh-based method, these are subject to the assumption of appropriate spacing between particles [16].
The SPH method is a meshless way of carrying out numerical simulations, in particular, fluid dynamics, based on Lagrangian description of the modeled liquid domain. Therefore, the SPH model is characterized by the fact that individual particles, which meet the roles corresponding to the nodes of Euler's FEM methods, move freely in terms of bonds imposed during the simulation. In this case, however, the mesh is not a restriction for their movement [28,30]. Single particles during simulated flow carry the necessary information, such as the value of speed or density. During the simulated flow, they interact with each other in accordance with the principles of fluid dynamics [31].
This calculation method is based on the theory of interpolation. Continuous distribution of such parameters as speed or liquid pressure are replaced by appropriate estimates with a certain interpolation kernel. Calculations are performed using a discreet set of a specific number of fluid particles. In the general approach, the estimation of any Fi in a certain position of the i particle is described as: where m j -the mass of the j particle, ρ j -density of the j particle, and W ij -kernel function, depend on the smoothing length h ij and distance between particles r ij . This method, using the general Equation (1) can be used to calculate typical tensor quantities, such as speed, pressure or stress. It should be noted, however, that the basis for calculations is information about density at specific points, represented by particles located in specific positions. Information about this parameter is used to estimate the physical quantity in request. Therefore, the authors decided that this method can be successfully used to estimate the force of extrusion of dry ice by the matrix, as shown by analyses of large deformations of material with elastic-plastic features.
This article presents the outputs of computer simulations of the process of dry ice extrusion using four different extrusion dies. The objective was to find the maximum value of F e for different shapes of the extrusion die. The predicted values were compared with the experimental data to assess the accuracy of representation of the process by numerical methods, primarily SPH based. The primary goal was to develop a relatively simple technique for estimating the maximum resistance during compression of solid dry ice that could be used to maximise the process efficiency when designing new dies. To this end, it was required to find a numerical solution, for example with the application of SPH, to model the process of extrusion of dry ice without needing to use sophisticated material models, which require extensive testing to determine the material properties. With this approach it will be possible to expand the research also on other materials.
In addition, the outcome of this research allowed the assessment of the influence of die shape parameters on the maximum force applied in the dry ice extrusion process.

Solid Carbon Dioxide
The compressed material dry ice snow was obtained by expansion of liquid carbon dioxide stored in a sealed container at a pressure of 20 bar and at −18 • C. Liquid to solid state transition occurs in an adiabatic process when liquid carbon dioxide undergoes a rapid expansion to atmospheric pressure. This results in crystallisation. The resulting dry ice snow has a density of 550 kg/m 3 [30].
Before the empirical tests were carried out, as part of this research, dry ice was stored in a foamed polystyrene container. The insulated container allowed the minimisation of the sublimation rate. These storage containers were also used to cool down the test assembly to a temperature close to the temperature of the tested material.

Single-Hole Dies
The tests and numerical simulations were carried out for four different single-hole dies shown in Figure 3, each having a different inlet section. The shape of the compression cavity of the tested dies was described with geometric parameters defining its lengthwise section. The following dimensions were the same for all the dies in consideration:  Spherical-cylindrical (WKWP) with the inlet section made up of two spherical frustums arranged in a series, convex followed by concave. This shape is defined by the sphere radii R1 and R2 and the heights of the spherical frustums h1 and h2 ( Figure  3d). The geometric parameters of the above-described types of compression cavities ( Figure 3) can be interrelated by mathematical relationships. These relationships can be written as follows: − for CS cavity: − and for cavity WKWP, with an additional parameter H calculated as the sum of h 1 and h 2 as follows: we can introduce a proportionality coefficient a: and using this coefficient to relate radii R 1 and R 2 to the other geometric parameters of the cavity: The values of the key geometric parameters, measured on the longitudinal section of the compression cavity are given in Tables 1-4 below.

Numerical Analysis
The numerical model was prepared in Abaqus 2020 software by Dassault Systèmes, with the simulations performed in Abaqus/Explicit. Figure 4 shows the model of CS die with a FEM mesh superimposed on it.

Numerical Analysis
The numerical model was prepared in Abaqus 2020 software by Dassault Systèmes, with the simulations performed in Abaqus/Explicit. Figure 4 shows the model of CS die with a FEM mesh superimposed on it.  The model is composed of three parts represented by finite elements (Figure 4). The main element is the single-hole die (1), integrated with the feed barrel to simplify the model, with a good effect on the uniformity of the finite element mesh. Inside the cylindrical part there is a compression ram (2), which moves by sliding with zero friction. The extruded material (3), i.e., dry ice snow, was modelled as an axisymmetric solid, corresponding in shape to the integrated single-hole die. Two reference points were defined: − Reference point of the single-hole die (4), fixed on the die (1), and located in the bottom plane on the axis of symmetry; − Reference point of the ram (5), fixed on the ram (2), and located on its top surface and on the axis of symmetry.
These reference points were used in the calculations to set the initial-boundary conditions relative to the ram and to the die.
The analysis was performed with the following parameters adopted for the respective aspects of the numerical model: The same concept was applied in the model of the ram that pressed the material through the die. Its characteristic parameters were:

−
Outside diameter of 27.5 mm and 2 mm thickness. This allowed the avoidance of friction between the ruled surface and the sides of the feed barrel. Furthermore, it was now possible to represent the process of extrusion with the ram moving inside a The same concept was applied in the model of the ram that pressed the material through the die. Its characteristic parameters were: − Outside diameter of 27.5 mm and 2 mm thickness. This allowed the avoidance of friction between the ruled surface and the sides of the feed barrel. Furthermore, it was now possible to represent the process of extrusion with the ram moving inside a cylindrical chamber, on the specially built test bench. − FEM mesh of R3D4 4-node, 3-D, quadrilateral, and infinitely rigid elements of 1 mm averaged edge length, uniformly sized over the whole extrusion ram surface.
A slightly different approach was applied in modelling the domain associated with the compressed dry ice ( Figure 6): − The compressed dry ice was modelled as a deformable elasto-plastic body; − The shape and size of this domain was assumed to be commensurate to the internal cavity of the integrated die ( Figure 5). This shape was assumed to change in the process of extrusion, as is typical of a deformable body; − FEM mesh of C3D8R linear, 6-node, 3-D, reduced-integration elements. The FEM elements had 2 mm averaged edge length and were uniformly sized over the whole domain; − SPH based approach was applied to define the properties of the compressed dry ice to allow the use of the SPH method. This ensures undisturbed flow of the compressed material through the die, leaving out nodal interactions. Thus, the individual finite elements are modelled as non-interacting particles (Figure 7).

Experimental Verification
The experimental tests were done using a test frame mounted on a specially designed The following contact conditions were assigned between all contacting surfaces: − Tangentially: µ = 0.1 friction, in the direction normal to the surface-no penetration; − With infinitely rigid walls of the integrated single-hole die (1) and ram (2) and "no penetration" boundary condition in place, the compressed dry ice (3) devoid of such conditions flew out through the opening in the bottom of the integrated die.

Experimental Verification
The experimental tests were done using a test frame mounted on a specially designed test bench (Figure 9). This test set-up included the MTS Insight 50 kN (MTS Systems Corporation, Eden Prairie, MN, USA) test frame (1), and the ram mounting plate, (4) supported by linear guides (5), was attached to the load cell (2) and grips (3) assembly. These guides are connected to the base (8), resting on the bottom support (9). The ram (6), which interacts with the compression barrel assembly (7), is fixed to the ram mounting plate (4). This assembly is composed of the feed barrel integrated with the compression cavity (1 in Figure 1), including the die (3 in Figure 1). In the experimental tests, loose dry ice snow, of a known weight, was forced through the assembly by the movement of the MTS crosshead displaced with a uniform speed. The observed parameter was the force vs. displacement.
The test procedure included the following conditions: In the experimental tests, loose dry ice snow, of a known weight, was forced through the assembly by the movement of the MTS crosshead displaced with a uniform speed. The observed parameter was the force vs. displacement.
The test procedure included the following conditions: − After every three test measurements, the compression barrel and ram were cooled in a container with dry ice at 195 K in order to reduce sublimation, during the process of extrusion, which could increase due to the ambient temperature of the laboratory (ca. 293 K), which was higher than the temperature of dry ice (195 K); − Coaxial alignment of the ram, feed barrel, and compression cavity to ensure no contact between the ruled surface of the ram and the feed barrel opening. The ram diameter was 2% smaller than the feed barrel diameter to avoid the risk of the metallic surfaces coming in contact due to thermal expansion; − The die was filled up with dry ice pellets before the first extrusion cycle, and after each cooling of the compression assembly, in order to stabilise the distribution of stress (the value measured in the first cycle was always left out); − Constant speed of extrusion of v e = 5 mm/s.; − 31 ± 1 g of dry ice snow was fed into the feed barrel each time; − Measurement of the applied force to an accuracy of 0.5% of the maximum measurement range of the test frame (0.5 accuracy class of the load cell fitted on the MTS test frame).
The cohesion of the produced pellets was assessed through sensory analysis, with the criterion being the lack of cracking, and glassy structure over the entire pellet ( Figure 10). The extrusion producing pellets affected by cracking or visible disfigurement were rejected. The test was repeated 16 times for each shape of the compression cavity, with each case producing acceptable pellets. The cohesion of the produced pellets was assessed through sensory analysis, with the criterion being the lack of cracking, and glassy structure over the entire pellet ( Figure  10). The extrusion producing pellets affected by cracking or visible disfigurement were rejected. The test was repeated 16 times for each shape of the compression cavity, with each case producing acceptable pellets.

Statistical Analysis of the Test Data
One-way analysis of variance was performed using the post hoc Tukey's test. Two populations of test data were compared within each type of die, i.e., CS, WK, WP, WKWP. The analysis was carried out using Statistica (version 13.3, TIBCO Software Inc., Palo Alto, CA, USA). All the comparisons were one-way ANOVAs, with statistical significance defined by the value of p < 0.05.

Discussion
Examples of the resistance force Fe vs. displacement curves obtained with the abovementioned test set-up are presented in Figure 11. Considering a relatively small scatter of data, these curves, based on their shape, can be sub-divided into phases, as was the case in previous studies (Figure 2), which attests to the accuracy of representation of the

Statistical Analysis of the Test Data
One-way analysis of variance was performed using the post hoc Tukey's test. Two populations of test data were compared within each type of die, i.e., CS, WK, WP, WKWP. The analysis was carried out using Statistica (version 13.3, TIBCO Software Inc., Palo Alto, CA, USA). All the comparisons were one-way ANOVAs, with statistical significance defined by the value of p < 0.05.

Discussion
Examples of the resistance force F e vs. displacement curves obtained with the abovementioned test set-up are presented in Figure 11. Considering a relatively small scatter of data, these curves, based on their shape, can be sub-divided into phases, as was the case in previous studies (Figure 2), which attests to the accuracy of representation of the already used test method [21]. ANOVAs were carried out for the observed data for each die type, in turn. As the first step, normality of the respective populations was checked using the Shapiro-Wilk test, taking into account only the maximum extrusion resistance force FL. The minimum value of p was 0.069 or more.
Next, homogeneity of variances in comparable populations was checked for the respective die shapes. To this end, the value of p was determined using the Brown-Forsythe test. The hypothesis of the homogeneity of variance was confirmed as the obtained values exceeded the 0.05 criterion.
In the final step of ANOVA, the value of p was calculated for the cases which obtained a value of index below 0.05. This allowed the rejection of the hypothesis that compared dies had the same value of the parameter under analysis, i.e., the extrusion resistance force FL. This supports the hypothesis with statistically significant differences between the compared populations.
Tukey's post hoc test was used to determine the statistical significance of the differences between the respective populations, the results of which are presented in Tables 5-8. Based on the data presented in Table 8 we see that in the WKWP group statistical significance is observed only between the WKWP70-50 and WKWP70-25 populations. The other pairwise comparisons had a probability of error higher than 5%, and thus should not be made.  ANOVAs were carried out for the observed data for each die type, in turn. As the first step, normality of the respective populations was checked using the Shapiro-Wilk test, taking into account only the maximum extrusion resistance force F L . The minimum value of p was 0.069 or more.
Next, homogeneity of variances in comparable populations was checked for the respective die shapes. To this end, the value of p was determined using the Brown-Forsythe test. The hypothesis of the homogeneity of variance was confirmed as the obtained values exceeded the 0.05 criterion.
In the final step of ANOVA, the value of p was calculated for the cases which obtained a value of index below 0.05. This allowed the rejection of the hypothesis that compared dies had the same value of the parameter under analysis, i.e., the extrusion resistance force F L . This supports the hypothesis with statistically significant differences between the compared populations.
Tukey's post hoc test was used to determine the statistical significance of the differences between the respective populations, the results of which are presented in Tables 5-8. Based on the data presented in Table 8 we see that in the WKWP group statistical significance is observed only between the WKWP70-50 and WKWP70-25 populations. The other pairwise comparisons had a probability of error higher than 5%, and thus should not be made.  The key statistical data obtained for the respective groups of dies are compiled in Table 8 and illustrated on the box plots in Figure 12, relating the median of the maximum observed extrusion force F L_EMP with the associated key statistical data.
The medians, which slightly differed from the relevant average values, were used in the analysis of differences between the observed and predicted data. Figure 13 shows an example of a numerical simulation of the process of extrusion of dry ice snow as a F e vs. s e relationship. As it can be seen on the graph, two phases of the simulated process can be distinguished, which can be assigned to the relevant experimental data: − Phase 2, in which the value of F e decreases due to elastic strain of the material and its frictional interaction with the internal surfaces of the compression cavity. This increase continues up to the maximum compression force F L, which marks the commencement of extrusion, − Phase 3, in which the value F e decreases due to the plastic strain of the extruded material, whose value depends on the shape of the compression cavity. The medians, which slightly differed from the relevant average values, were used in the analysis of differences between the observed and predicted data. Figure 13 shows an example of a numerical simulation of the process of extrusion of dry ice snow as a Fe vs. se relationship. As it can be seen on the graph, two phases of the simulated process can be distinguished, which can be assigned to the relevant experimental data: − Phase 2, in which the value of Fe decreases due to elastic strain of the material and its frictional interaction with the internal surfaces of the compression cavity. This increase continues up to the maximum compression force FL, which marks the commencement of extrusion, − Phase 3, in which the value Fe decreases due to the plastic strain of the extruded material, whose value depends on the shape of the compression cavity. Note that this simulation with elasto-plastic model lacks a distinguishable Phase 1, identified during actual extrusion of the material. This is due to the simplification characteristics of the adopted elasto-plastic model, where the material is modelled as compressed right from the start, thus leaving out the effect of compression. That said, the Note that this simulation with elasto-plastic model lacks a distinguishable Phase 1, identified during actual extrusion of the material. This is due to the simplification characteristics of the adopted elasto-plastic model, where the material is modelled as compressed right from the start, thus leaving out the effect of compression. That said, the maximum force applied to extrude the material through the die F L is the most interesting parameter at this stage of research.
An interesting phenomenon is the oscillation of the value of the F e material force by the matrix, observed at the end of stage 2 and throughout Stage 3, while after exceeding the displacement s e about 30 mm, this oscillation significantly reduces its amplitude. This phenomenon is due to the stick-slip effect, which is associated with the elasto-plastic representation of the extruded material. In the area of pure occurrence of this phenomenon (approx. 8 mm < s e < approx. 30 mm), the extruded material does not extend from the matrix at total speed, equal to the speed of the piston movement v e . In this respect, s e displacement, due to the elastic deformation of its volume, the speed of the material flowing from the matrix is significantly lower. Considering the coincidence of this phenomenon, with friction of the material on the walls of the matrix and the internal friction of the deposit, are good conditions for the resulting effect of the stick-slip effect-i.e., periodic fluctuations in movement speed, as a result of the mutual impact of friction force and force responsible for pressing. The reduction of its intensity (decrease in the amplitude of the oscillation of force) is present after exceeding a certain displacement (s e = approx. 30 mm), when a portion of the compressed material (which is deformed in an elastic way), which is initially placed in front of the coincidence, begins to fill its volume. At the same time, the material originally present, leaves the matrix. Therefore, part of the material previously deformed in an elastic way (already in the range of plastic deformation) is less susceptible to this effect. Therefore, this parameter was analysed using the data predicted by the numerical simulation F L_FEM and the median of the experimental data F L_EMP . The results are compiled in Table 9 below. Table 9. Maximum predicted extrusion force F L_FEM and its average value calculated from the experimental data F L_EMP . The results are also represented graphically in Figure 14 as relationships of the maximum extrusion resistance force F L , both F L_FEM (predicted) and F L_EMP (observed), and the relevant geometric parameters of the respective die shapes (Tables 1-3). A relatively high similarity of the obtained curves was observed in terms of monotonicity. Unfortunately, the values noticeably differ between the test points in all cases. This can be explained by a fixed value of yield stress σ pl used in all the numerical simulations, which was determined experimentally on compressed material. This implies inadequacy of the applied procedure as in the analysed system yield stress apparently depends on the type of die used in the process. However, this assumption appears correct as long as the above-mentioned value is taken as a process-specific parameter rather than an intrinsic property of the material. of die used in the process. However, this assumption appears correct as long as the abovementioned value is taken as a process-specific parameter rather than an intrinsic property of the material. This hypothesis needed to be verified. To this end, at all the test points (i.e., the values of the geometric parameters of the respective dies and their corresponding values of FL_FEM and FL_EMP) mathematical functions were derived, describing the percent difference between the obtained values in relation to the observed extrusion resistance:

Die Designation Predicted Maximum Extrusion Force
and were related to the geometric parameters of the dies under analysis ( Figure 14). The relationship obtained in this way showed a difference reaching as much as 40%, which excludes the data predicted with this numerical model from further analysis and, more importantly, from further work on a simplified model of the dry ice extrusion process with the application of SPH. However, this difference served to derive the corrective function for the yield stress σpl used in the numerical elasto-plastic model from the following equation: In this way, corrected yield stresses σpl_COR were obtained for the test points set on all the analysed dies ( Figure 15). Relationships describing the variation of this parameter as This hypothesis needed to be verified. To this end, at all the test points (i.e., the values of the geometric parameters of the respective dies and their corresponding values of F L_FEM and F L_EMP ) mathematical functions were derived, describing the percent difference between the obtained values in relation to the observed extrusion resistance: and were related to the geometric parameters of the dies under analysis ( Figure 14). The relationship obtained in this way showed a difference reaching as much as 40%, which excludes the data predicted with this numerical model from further analysis and, more importantly, from further work on a simplified model of the dry ice extrusion process with the application of SPH. However, this difference served to derive the corrective function for the yield stress σ pl used in the numerical elasto-plastic model from the following equation: In this way, corrected yield stresses σ pl_COR were obtained for the test points set on all the analysed dies ( Figure 15). Relationships describing the variation of this parameter as a function of the key geometric parameters of the analysed dies, were derived in addition, considering the future application of this approach in optimisation of the geometric parameters of extrusion dies over the entire variation domain.  The obtained corrected yield stress values σpl_COR were used in the repeated numerical calculations performed with the use of the SPH method. Other settings remained unchanged. The results of the comparison between the corrected extrusion resistance force FL_FEM_COR and the observed values of this force FL_EMP are shown in Figure 16, together with the percent difference between these two ΔFL_COR. As it can be seen, this treatment improved the accuracy of representation of the process by numerical analysis, with regards to the obtaining of the maximum extrusion resistance force. After correction of the initial yield stress the percent difference did not exceed 10%, and the curves of FL_FEM_COR and FL_EMP became more similar in terms of monotonicity (as compared to the curves obtained for the constant value of σpl). σ pl_COR = -0.0096·α 2 + 0.2527·α + 3.0632 The obtained corrected yield stress values σ pl_COR were used in the repeated numerical calculations performed with the use of the SPH method. Other settings remained unchanged. The results of the comparison between the corrected extrusion resistance force F L_FEM_COR and the observed values of this force F L_EMP are shown in Figure 16, together with the percent difference between these two ∆F L_COR . As it can be seen, this treatment improved the accuracy of representation of the process by numerical analysis, with regards to the obtaining of the maximum extrusion resistance force. After correction of the initial yield stress the percent difference did not exceed 10%, and the curves of F L_FEM_COR and F L_EMP became more similar in terms of monotonicity (as compared to the curves obtained for the constant value of σ pl ).

Conclusions
The main purpose of the on-going research by the authors is to develop a relatively simple technique of estimating the maximum resistance during compression of solid dry ice that could be used in the design of extrusion dies, helping to improve the efficiency of the extrusion process, i.e., to minimise the input of energy while ensuring high quality product. The above-mentioned overall objective motivated the authors to look for an appropriate numerical modelling method that could be used in this connection. Therefore, it was assumed that it would be possible to represent dry ice with an elasto-plastic material assigned with the properties of dry ice pellets, thus avoiding complex models used in simulations of the process of compression. This approach also minimises the number of material properties that must be determined experimentally as an input for the calculations.
The FEM model prepared as part of this research represents the extrusion of dry ice snow through the die, with omission of the compression phase, designated here as Phase 1. Still, owing to the correct representation of phases 2 and 3, this approach can be successfully used in research dealing with the determination of the maximum compression resistance. Therefore, we can conclude that for this aspect (determination of the maximum extrusion resistance), the process of compression may be simulated by extrusion of compressed material through the die.
Yet another issue to be dealt with is the adoption of an appropriate parameter defining the plastic strain of the domain extruded through the die. The proposed SPH model   Figure 16. Comparison of the curves representing the relationship of the maximum force during dry ice extrusion and the adopted variable geometric parameters of the analysed single-hole dies for: corrected maximum predicted force (F L_FEM_COR ), calculated with the corrected yield stress σ pl_COR , maximum observed force (F L_EMP ), and their percent difference ∆F L_COR for the die groups CS (a), WK (b), WP (c), and WKWP (d).

Conclusions
The main purpose of the on-going research by the authors is to develop a relatively simple technique of estimating the maximum resistance during compression of solid dry ice that could be used in the design of extrusion dies, helping to improve the efficiency of the extrusion process, i.e., to minimise the input of energy while ensuring high quality product. The above-mentioned overall objective motivated the authors to look for an appropriate numerical modelling method that could be used in this connection. Therefore, it was assumed that it would be possible to represent dry ice with an elasto-plastic material assigned with the properties of dry ice pellets, thus avoiding complex models used in simulations of the process of compression. This approach also minimises the number of material properties that must be determined experimentally as an input for the calculations.
The FEM model prepared as part of this research represents the extrusion of dry ice snow through the die, with omission of the compression phase, designated here as Phase 1. Still, owing to the correct representation of phases 2 and 3, this approach can be successfully used in research dealing with the determination of the maximum compression resistance. Therefore, we can conclude that for this aspect (determination of the maximum extrusion resistance), the process of compression may be simulated by extrusion of compressed material through the die.
Yet another issue to be dealt with is the adoption of an appropriate parameter defining the plastic strain of the domain extruded through the die. The proposed SPH model with the adopted constant yield stress σ pl (determined by testing the compressed material) failed to accurately predict the maximum extrusion resistance force for the different die parameters. To cope with this problem, yield stress was treated as a process specific variable, depending on the geometric parameters of the die rather than a constant property of the extruded material. The corrected yield stress σ pl_COR variable in the domain of geometric parameters of the die was used, and the corrected FEM model reduced the difference between the predicted and observed values to less than 10%, thus solving the problem. Being a gross simplification of the actual process of extrusion, with the above-mentioned level of error, the accuracy of the proposed elasto-plastic model with the use of SPH should be considered satisfactory. This is especially true when the focus is on the maximum extrusion force rather than on the simulation of the entire process, as it is the case in this research. In addition, the tests performed on dry ice snow do not yield stable results, mainly due to technical constraints, with sublimation of this material being the main challenge. This is evidenced by the problem with reducing the probability of error to an acceptable level when comparing the subsequent data populations. This problem observed in laboratory experiments also applies to dry ice production on an industrial scale, i.e., under generally less controlled conditions.
Considering the above constraints, it is reasonable to conclude that dry ice may be correctly represented by elasto-plastic material in numerical simulations of the process of dry ice extrusion for the purposes of studies on the efficiency of this process. Furthermore, subject to appropriate input assumptions (value of the yield stress parameter), the above simulation can yield satisfactory results confirmed by experimental data. The authors plan to apply this approach to modelling the process of extrusion in future research, analysing the influence of the shape of the die cavity on the energy requirements of the process (by analysing the resistance force), over continuously varying domains, covering the geometric parameters of the die and the associated shape optimisation efforts. However, an additional parameter should be introduced for this purpose, that would allow verification of the quality of the produced pellets. Funding: This research is a part of the project: "Developing an innovative method using the evolutionary technique to design a shaping die used in the extrusion process of crystallized CO 2 to reduce consumption of electricity and raw material", number: "LIDER/3/0006/L-11/19/NCBR/2020" financed by the National Centre for Research and Development in Poland, https://www.gov.pl/web/ncbr (accessed on 15 December 2021).
Institutional Review Board Statement: Not applicable.