Numerical Simulations of Polymer Devolatilization in a Steam Contactor: Correlating Particle Distribution With Volatile Content Removal

In polymer manufacturing processes, residual monomers, solvents, and other undesirable substances can deteriorate the quality of the final products and pose health and environmental risks. Therefore, the polymer industry strictly regulates the concentration of permitted compounds, especially in products like food packaging, surgical prostheses, and healthcare items, to ensure safety and compliance with health standards. Removal of volatiles from polymers, or devolatilization, is typically performed upstream of shaping operations such as extrusion, pultrusion, and injection molding, to ensure the quality of the final product. Polymers include highly chained hydrocarbons like polyethylene [1], polystyrene [2], rubber [3], and polydimethylsiloxane [4], whereas volatiles released in these processes are mostly monomers such as octene [1], cyclohexane [5,6], and trifluoroethane [4]. Common devolatilization methods include vacuum degassing, thermal degassing, and vapor–liquid stripping [2].

The current study investigates the polymer devolatilization of an industrial-grade thermoplastic polymer via the vapor–liquid stripping method, with cyclohexane as the volatile and steam as the stripping agent. A steam contactor is a device used to achieve devolatilization by mixing steam and cement, the latter, a term referred to as a mixture of a polymer and cyclohexane. The cement experiences shearing and fragmentation within the steam contactor due to the force exerted by the steam. Concurrently, the high pressures and temperatures within the contactor lead to the rapid evaporation of the cyclohexane solvent [7]. This rapid evaporation, known as flashing, leads to the formation of polymer particles, which are then carried by the steam toward the outlet of the contactor. In addition, the stability of the steam enables the effective removal of the polymer from the solvent without significant clustering [8].

The primary aim of a devolatilization device such as an extruder or a mixing chamber is to maximize the evaporation of volatile components from the polymer. The extruders are usually used for general polymers and can be of single screw [9–12] or twin screw type extruders [1,13–16]. Mixing chambers are observed to be useful for studying devolatilization in rubber-based compounds [3]. Many studies have primarily focused on the convective mass transfer aspects of the devolatilization process [3,4,17]. For instance, the mass transfer rate during devolatilization was characterized using a liquid phase mass transfer coefficient, a driving force, and the effective surface area involved in devolatilization [3]. Dierkes and Noordermeer [3] studied rubber mixing with silica and silane, where high shearing aids dispersion, but risks scorch, while silanization benefits from high temperatures. They developed a devolatilization model based on penetration theory that accounted for factors such as temperature, mixer volume, fill factor, rotor speed, reaction time, and the partial pressure of volatile components. Woyuan et al. [17] developed a mass transfer model for the devolatilization process of highly viscous syrup in a rotating packed bed based on penetration theory and mass conservation with acetone as the volatile organic compound. It was found that bubbling played an important role in the devolatilization. Hirschfeld and Wünsch [4] systematically studied bubble-free polymer devolatilization using experimental degassing tests and surface renewal analysis with a simplified extruder model. Their work highlighted the importance of mixing and gas-side mass transfer resistance, validated by numerical simulations with the open-source computational fluid dynamics (CFD) software, openfoam. More recently, a new devolatilization process in a twin screw extruder with additives was proposed and it was observed that the devolatilization efficiency increased with the vacuum degree and devolatilization temperature [1].

Given the complexities inherent in multi-physics processes involving fluid flow, heat transfer, and multiphase interactions within such manufacturing systems and the high temperatures and pressures involved, conducting laboratory experiments to investigate these issues is very challenging. In contrast, CFD simulations can be conducted under various operating conditions and parameters with relative ease [5]. Although computational studies on devolatilization are limited, there are related studies on similar principles such as spray drying and gasification, where heat is transferred into the polymer by passing it through a hot flow of steam or inert gases. The study of Lou et al. [18] reviewed the state-of-the-art progress in detailed numerical simulations of atomization and evaporation with the level set method. They reviewed improvements in methodology, such as better mass conservation and coupling with the volume of fluid method, and discussed challenges like droplet resolution and physical mechanism understanding. Xu et al. [19] compared the Lagrangian and Eulerian approaches under identical operating conditions, focusing on their performance in predicting particle deposition in turbulent flows. They found that the Eulerian approach for particles required less computational time and cost compared to the Lagrangian approach while providing higher accuracy. Al-Qadasi and Ozkan [20] conducted a parametric study on the drying stage of biomass steam gasification. They examined the effects of temperature and steam-to-biomass ratio on the components of the produced gas and carbon conversion efficiency. Although their study did not involve the evaporation of a polymer solution, it did address chemical decomposition at high temperatures. The model was developed based on a comprehensive pyrolysis reaction model by adopting the CFD-based Eulerian–Eulerian approach.

Other studies [21–24] have focused on multiphase simulations of shock-driven droplet breakup and evaporation. In Ref. [21], the effects of multiphase parameters on a shock-driven particle-laden hydrodynamic instability were considered through simulations in a shock tube environment as the computational domain. Multiphase physics included momentum and energy coupling, with the study discussing the addition of mass coupling through evaporation. Dahal and McFarland [22] developed a numerical method for predicting the interaction of active, phase-changing particles in a shock-driven flow. The particle-in-cell technique was used to couple particles in a Lagrangian coordinate system with a fluid in an Eulerian coordinate system. The effect of the particle size distribution (PSD) and particle evaporation was examined further for this case. They showed that larger particles reduced the vorticity deposition, while particle evaporation increased it. Finally, studies of Duke-Walker et al. [23,24] used experimental measurements and simulation results to highlight the importance of breakup parameters on the evaporation rate and large-scale mixing in shock-driven multiphase instability. Overall, it was shown that the large-scale hydrodynamics instability enhanced the evaporation and the shock-driven multiphase instability. It was also shown that the breakup of the droplets significantly increased the strength of the instability, and rate of droplet evaporation.

About the specific application in this study, there have been very few studies [5,6,25,26] involving polymer devolatilization in steam contactors. For instance, Gabor et al. [5] examined the impact of geometry modifications such as alterations in gap length, width, and mixing zone depth on various parameters like the extent of cyclohexane phase evaporation and PSD. It was observed that a narrower gap and a wider mixing zone resulted in improved performance. A parametric study [25] was carried out to determine the benefits of altering the operating steam pressure and several geometric parameters of the steam contactor. Later, Shindle and Chandy [6,26] through two different sequential calculations conducted detailed parametric studies to determine the effects of different contactor geometries. The first step involved the simulations of the primary breakup of the polymer sheet and the second step then used the resulting diameter distribution to model the multiphase heat and mass transfer of the polymer mixture including evaporation.

In the current study, 3D CFD calculations are carried out using an Eulerian–Lagrangian approach, wherein superheated steam is represented as the continuous phase tracked in an Eulerian frame, and polymer–solvent mixture droplets are handled via Lagrangian tracking, as was the case in the previous studies [5,6,25]. However, while these previous studies have explored how the inlet parameters impact the final drying of the polymer mixture, they have not—and in fact, no other study in the literature has—examined internal particle dynamics and its effect on devolatilization. For instance, the distribution of particles within such a system can be critical in achieving the required devolatilization. Achieving maximum cyclohexane evaporation requires each particle to receive sufficient heat from the steam. In an ideal scenario where all the particles are equally sized and have identical weights, if the particles are clustered together as lumps, then the uneven heat can result in insufficient evaporation of the volatiles, which then produces more wet polymer at the exit of the contactor. On the other hand, if the particles are evenly spaced within the contactor, the heat would be equally distributed. This argument can be further extended to particles with varying sizes and masses, as is the scenario in the contactor presented in the current study. It is expected that particles with larger masses would require a larger amount of steam for better evaporation, and here too, the extent of evaporation would depend on the spacing between the particles or in other words, the clustering of the particles. These features are further analyzed in detail through particle statistics presented here.

So the objective of the current study is to conduct a series of Eulerian–Lagrangian 3D computations at varying geometric and physical conditions to analyze the particle distribution and its effect on the devolatilization process. Two metrics related to the particle clustering within the contactor are developed as part of this study. One of them is the cluster distribution index (CDI) employed in previous polymer mixing-related simulations for assessment of distributive characteristics [27–32], and the second is a mass-based CDI calculated from the radial distribution of masses developed and presented here for the first time. The aim is to determine the factors and conditions that allow for the optimization of steam usage in terms of its cost and sustainability. By examining and studying these phenomena, valuable insights can be obtained to enhance the efficiency and effectiveness of steam utilization in the contactor system. The governing equations are explained in the next section followed by a description of the problem and details of the numerical method employed. Results from the computations consisting of grid-independence analysis, validation with sieve data, Eulerian predictions, and particle distribution are described next. The paper concludes with a summary of the findings.

The dimensionless parameter Reynolds number, $Re=ρVL/μ$⁠, is employed to characterize the flow, where $ρ$ is the density, $V$ is the fluid velocity, $L$ is the characteristic length, and $μ$ is the dynamic viscosity of the fluid. The Re calculated in this study ranges from $3×105$ to $1.3×106$⁠, based on the steam flow conditions and a characteristic length equal to the inlet diameter of the contactor. Hence, the flow is considered to be turbulent. In the current study, the Favre-averaged Navier–Stokes (FANS) approach [33], which models all the turbulent scales, is employed. Additionally, due to the high steam velocity in the steam contactor, the flow is compressible, with the steam inlet Mach number varying from $Ma=V/c=0.14$–$0.84$⁠, implying a subsonic inlet. The maximum steam Mach number in the domain is around $3.2$⁠, indicating supersonic flow.

The cyclohexane undergoes evaporation into a gas phase due to the high temperature and pressure from the steam inside of the contactor. This is incorporated through pressure-based boiling and multicomponent heat mass solutions enabled in ansys fluent [36]. Multicomponent particles are treated as a mixture of species (polymer and cyclohexane) within the cement particles. The vaporization or boiling model follows a convection/diffusion-controlled approach, and the equilibrium between vapor and particles is determined using Raoult’s law. So the correlation between the vapor concentration of a species $Ci,s$ over the surface and its mole fraction in the condensed phase $XiL$ is $Ci,s=pi/RT$$=XiLp/RT$⁠, where $p$ is the pressure, $R$ is the gas constant, and $T$ is the temperature. See Ref. [36] for more details.

The current investigation focuses on simulating a polymer devolatilization process within a 1.34-m-long steam contactor, generating distinct solid polymer particles. The steam contactor’s configuration, as depicted in Fig. 1, involves the introduction of superheated steam through axial inlets and the simultaneous injection of a cement mixture consisting of polymer and cyclohexane through a circumferential slit, as illustrated. Upon their interaction, the high momentum of the superheated steam leads to shearing, along with mass and heat transfer between the steam and the polymer. This triggers a flashing of the polymer mixture, resulting in solvent evaporation and the formation of polymer particles. These particles continue to mix with the steam throughout the length of the contactor, until they reach the outlet, eventually entering the polymer isolation coagulator (not included in the current analysis). The assumptions involved in the current computational study are the following:

1. Flow is assumed to be steady.

2. Steam is treated as an ideal gas and any real gas effects are neglected.

3. The polymer is treated as a Newtonian fluid.

4. Primary breakup of the polymer liquid is neglected and the liquid enters contactor as discrete parcels.

This section describes the computational details of the simulations presented in this paper. The study involves the simultaneous solution of the FANS and species transport equations coupled with the Lagrangian-based DPM equations (as described earlier) in ansys fluent. The simulations here also employ a species transport model in ansys fluent to solve for both the polymer and cyclohexane in cement.

A pressure-based solver is employed for all simulations, utilizing a semi-implicit method for pressure-linked equations (SIMPLE) pressure–velocity coupling scheme and second-order upwind schemes for the discretization of various governing equations including continuity, momentum, energy, species transport, and turbulence model equations. The DPM incorporates an accuracy control method, wherein the integration step error is calculated at each instance, and if deemed excessive, the integration step is adjusted accordingly. This adaptive integration time-step approach ensures adherence to a specified tolerance level [36]. Given the focus on analyzing the composition of particles exiting the system, an escape boundary condition is applied at the outlet. This allows for the assessment and comparison of statistics associated with escaped particles across different cases.

A DPM injection time-step of $10−3$ s was determined to be suitable for this application. This time-step, utilized during the computation of particle trajectories within the DPM framework, serves solely for numerical purposes, facilitating the advancement of discrete phase particles in the simulations. The trajectories are computed within a Lagrangian framework, necessitating unsteady integration of the particles [36]. Each particle is allotted a maximum of 45,000 time-steps, with a step length factor of 5, before its trajectory is terminated. This precaution ensures that particles are not indefinitely trapped in the contactor due to re-circulation zones or other obstructions. These settings guarantee that all injected particles are tracked throughout the entire domain, with minimal computational overhead.

The boundary conditions have been outlined and presented in Table 1 for all the surfaces depicted in Fig. 1. Mass flow inlet boundary conditions with values based on specific cases are applied to the steam and cement inlets. The walls of the contactor are modeled as non-penetrating and no-slip surfaces for the continuous phase. The exit of the contactor remains open to the ambient atmosphere, thereby allowing a pressure outlet boundary condition to be set to a gauge pressure of $0$ Pa. To accommodate the system’s symmetry, a symmetry boundary condition is applied to the plane of symmetry, as depicted in Fig. 1. In the context of particle dynamics, the cement inlet surface serves as the injection inlet for particles and the steam inlet allows for the particles to escape. The injection inlet employs a specific uniform mass flowrate and temperature (based on the case), but with a variable particle diameter spanning across a range established by the sieve analysis (see Sec. 4.4). The contactor walls are configured as reflective surfaces for the particles, i.e., the particles bounce back upon contacting the walls. Eventually, the particles exit via the outlet surface. A hybrid initialization method is employed in the simulation to determine the velocity and pressure fields by solving Laplace’s equation. Other variables such as temperature, turbulent kinetic energy and dissipation, and species mass fractions are automatically initialized by patching them with domain-averaged values.

In this study, the particles are introduced as liquid droplets, with diameters spanning from 0 to 3000 $μ$m and are segmented into suitable discrete intervals, each represented by a mean diameter for which trajectory calculations are performed. To characterize the distribution of droplet sizes, a Rosin-Rammler expression [36] is used. Rather than employing the entire range from 0 to 3000 $μ$m as a single injection, this study opts for multiple injections to enhance control over the particle injection process. The selected diameter tiers align with those utilized in the sieve analysis (presented later), i.e., 0–600 $μ$m, 600–725 $μ$m, 725–850 $μ$m, 850–1180 $μ$m, 1180–1700 $μ$m, and 1700–2360 $μ$m, with respective weighted values of 0.5, 0.456, 0.0122, 0.0068, 0.012, and 0.013 assigned to each tier. These weighted values, obtained through iterations to minimize the net error below 10%, are used for all the cases presented in this study.

FANS-based calculations of steam coupled with the DPM-based calculations of cement (as described in Sec. 2) are carried out for contactors of different sizes and a variety of flow inlet conditions using ansys fluent. A total of 69 cases constituting different combinations of steam inlet diameter, cement inlet gap size, steam inlet flowrate, steam inlet temperature, cement flowrate, inlet polymer content, cement inlet temperature, and cement injection velocity are carried out in this study. Note that with the cement being modeled as a discrete phase, both the flowrate and injection velocity need to be specified. This section begins with a grid-dependent study that assesses the validity of the mesh employed in the calculations. Thereafter, the computational model is validated by comparing results from 12 cases (with different flow conditions) to industrial-grade sieve analysis data. This is then followed by a presentation of the metrics, namely the CDI and mass-CDI, that can quantify the particle distribution to assess the efficiency of the contactor. The section ends with a demonstration of the effectiveness of CDI and mass-CDI in assessing devolatilization for the different cases.

Figures 2–5 present the geometry and flow conditions for all 69 cases in this paper. Figure 2 shows the contactor steam inlet diameter and cement inlet gap sizes on the left and right sides of the $y$-axis, respectively, against the case number on the $x$-axis. Diamond symbols represent the diameter and circle symbols represent the gap size. A similar presentation is done for steam inlet flowrate and steam inlet temperature in Fig. 3, cement inlet flowrate and inlet polymer content in Fig. 4, and cement inlet temperature and cement injection velocity in Fig. 5.

A grid-dependent study is necessary to validate the entire numerical approach, the employed mesh resolution, solution techniques, and the type of boundary conditions. To achieve this, a grid-convergence study is conducted, involving simulations with three different meshes but identical operating conditions, ensuring that the mesh is adequate to resolve all flow and temperature gradients. Each simulation is carried out using the methodology described before. A typical mesh used for the simulations in the current study is presented in Fig. 6. This part of the study compares results from three different grids: grid 1 with 34,500 cells, grid 2 with 81,260 cells, and grid 3 with 257,300 cells. The corresponding CPU hours consumed for these grids are 40, 180, and 260 h, respectively. The case considered here for grid independence is case 59 (see Figs. 2–5). Figure 7 compares the axial profiles of velocity magnitude, temperature, pressure, and density in the continuum phase between the simulations using the three different grid densities. It can be seen that the difference between the quantities keeps decreasing from grid 1 to grid 2 to grid 3. The differences between grid 2 and grid 3 are minimal, indicating that either grid 2 or grid 3 are appropriate choices for the mesh resolution in this study.

With the grid now chosen for the CFD calculations, the objective in this section is to validate the accuracy of our simulations by comparing results to sieve data more comprehensively. The focus is specifically on matching the PSD obtained from industrial-grade sieve analysis with the results obtained from the CFD simulations for several different cases. This can help in assessing the reliability of the computational model. For validation, two distinct contactors with fixed diameters and gap sizes, and flow inlet conditions as shown in Figs. 2–5 for cases 1–24, are employed. These contactors have diameters of 4 in. (cases 1–12) and 5 in. (cases 13–24), with cement inlet gap sizes of 0.19 in. (cases 1–12) and 0.35 in. (cases 13–24), respectively. To validate the model, the particle diameter size distribution is compared with the sieve data, utilizing the aforementioned two contactor sizes.

Sieve data determine the PSD of granular/particle materials. The process involves passing a representative sample obtained at a specific set of operating conditions, through a sequence of stacked sieves, each with progressively smaller openings as shown in Fig. 8(b). These sieves are arranged in descending order of mesh size, with the largest at the top and the smallest at the bottom. Following sieving, the material retained on each sieve is weighed, enabling the determination of cumulative and individual weight percentages for each particle size fraction. Since the production of the crumbs is continuous, the samples are taken in several minutes and $%$ of particle size collected in every sieve is measured. These data are then utilized to construct a histogram curve representing the PSD for a specific case.

Here, $Si$ and $Ei$ represent the simulated and experimental values for the PSD, respectively, while $E¯$ signifies the mean value of the experimental PSD values. The RRMSE reflects the average relative magnitude of discrepancies between predicted and observed values within a dataset. The RRMSE and $R2$ values for the validation data depicted in Fig. 9 are $13%$ and 0.97, respectively, which indicate that the model’s predictions closely align with the true values and a strong linear correlation between the predicted and measured values is obtained. These features further affirm the accuracy and reliability of the model.

Figures 10(a) and 10(b) show contours of the velocity and temperature of steam (continuum phase), respectively, whereas Figs. 10(c) and 10(d) show the cement particle distribution (discrete phase) colored by contours of the particle diameter and mass fraction of cyclohexane content, respectively, for case 2. The cement particles are introduced at a velocity significantly lower than that of the steam. The interaction of the axial flow of steam and the radial injection of the cement causes a significant increase in the velocity magnitude of the steam near the cement inlet. The reason for this behavior is as follows: The injection of the particles into the flow domain results in a reduction in the effective passage area of the supersonic flow of steam, thereby inducing a shock and drastic changes in properties. This is supported by Fig. 10(b), where a sudden temperature change is evident (the temperature changes from $150∘C$ to $60∘C$ drastically around the cement inlet region). Further evidence is provided by Fig. 7 presented earlier, where it can be seen that the shock at $X=0.43$ m resulting from the cement injection caused the axial profiles of velocity and density to increase, and temperature and pressure to decrease. Downstream of the cement injection region, there is a notable transfer of momentum and heat, causing particles to follow the path of the shock and leading to flashing and subsequent evaporation of the volatile cyclohexane within the particles. This is affirmed by the fact that the particles do not occupy all the spaces and the spread of particles is not uniform as shown in Figs. 10(c) and 10(d). Again, axial temperature profiles shown in Fig. 7(d) show a decrease at $X=1.1$ m, further corroborating this behavior as well.

Figure 10(c) shows the distribution of the particles colored by contours of particle diameter within the steam contactor. It can be observed that particles initially tend to move toward the core of the contactor after injection and then spread more across the cross section. Downstream, particles with less mass move closer to the walls, while larger particles tend to reside in the core. This behavior can be attributed to inertia, where larger particles move slower and are more resistant to sudden changes in direction compared to lighter particles. Additionally, as particles progress downstream, they undergo evaporation, resulting in a size reduction. Figure 10(d) reveals that evaporation is non-uniform and inconsistent. Comparing Figs. 10(c) and 10(d), it can be deduced that smaller particles undergo more rapid evaporation than larger particles. This phenomenon is attributed to the fact that smaller particles possess proportionally a larger surface area-to-volume ratio in comparison to larger particles. This facilitates higher evaporation by providing more sites for molecules to escape from the particle into the surrounding steam. Additionally, larger particles tend to have a significant portion of their volume shielded from direct contact with the steam, further accentuating the evaporation disparity between smaller and larger particles.

With the model successfully validated with sieve data using cases 1–24, an additional 45 cases (cases 25–69) with varying geometrical and flow configurations as seen in Figs. 2–5 are studied. To further illustrate the crucial role played by the distribution of particles in the extent of evaporation of volatile elements for various boundary and geometric conditions, Fig. 11 presents the variation of the increase in polymer content with the temperature gradient between the incoming steam and cement. This figure shows that as the temperature gradient increases, the polymer content also increases, indicating a higher rate of evaporation of the volatile component. However, there are instances where evaporation is lower despite a high-temperature difference, indicating that a more thorough examination of particle distribution is probably needed to analyze methods to control volatile content evaporation.

So a set of parameters that can enable a comparison between the real or actual and ideal or optimum distributions of particles, with the latter representing the maximum efficiency, is developed here. The ideal distribution refers to the scenario where particles of equal mass are uniformly spread throughout the domain, with the main focus being the exit region. This ideal distribution can be developed in two distinct ways: one involving mass-less particles with only positional significance and another involving particles with mass but less emphasis on their positions. Two quantities of interest can be formulated here based on these factors: CDI or cluster distribution index and mass-CDI or cluster distribution index based on mass. CDI depends on the position of the particles, namely the $x$ and $y$ coordinates, while mass-CDI is based on mass and radius measured from the centerline of the contactor (not individual $x$ and $y$⁠). Figures 2–5 illustrate the geometrical and boundary conditions for cases 25–36, displaying the 12 cases used in the analysis. These include steam contactors with diameters of 3, 4, and 5 in. Each diameter consists of four variations in gap sizes, as mentioned.

Figure 12 illustrates the positions of the final 10% of particles exiting the contactor across all 12 cases, projected onto the exit plane. The sub-figures are organized as follows: from left to right, i.e., from $a$–$b$–$c$⁠, $d$–$e$–$f$⁠, $g$–$h$–$i$⁠, and $j$–$k$–$l$⁠, the contactor diameter increases, and from top to bottom, i.e., $a$–$d$–$g$–$j$⁠, $b$–$e$–$h$–$k$⁠, and $c$–$f$–$i$–$l$⁠, the gap sizes increase. Analysis of Fig. 12 reveals that as the gap size increases, the particles tend to align themselves in an orderly manner and arrange in concentric circles near the exit region. As the particles are injected at the gap in a direction that is perpendicular to the flow, larger gap sizes result in reduced injection velocities, allowing more time for the particles to align themselves with the axial flow. Consequently, higher gap sizes exhibit a more clustering of particles. So for the same inlet diameter, higher gap sizes lead to low injection velocities of steam, causing much less mixing of particles with steam and resulting in particle clustering. This trend remains consistent across different contactor diameters (from top to bottom, i.e., $a$–$d$–$g$–$j$⁠, $b$–$e$–$h$–$k$⁠, and $c$–$f$–$i$–$l$⁠), specifically cases 25, 28, 31, 34, cases 26, 29, 32, 35, and cases 27, 30, 33, 36. The particles are injected from the wall of the contactor toward the center. Particles are expected to disperse away from the center as they align with the flow direction for the maximum spread of the particles in the domain. Cases with smaller inlet diameters exhibit particles closest to the center, as the high axial velocities prevent them from spreading throughout the exit (from left to right), i.e., $a$–$b$–$c$⁠, $d$–$e$–$f$⁠, $g$–$h$–$i$⁠, and $j$–$k$–$l$⁠. With an increase in contactor diameter (from top to bottom), i.e., $a$–$d$–$g$–$j$⁠, $b$–$e$–$h$–$k$⁠, and $c$–$f$–$i$–$l$⁠, the steam velocity decreases, facilitating easier penetration of particles into the steam flow.

To quantify the spatial distribution of particles as depicted in Fig. 12 in terms of CDI, the distribution of particle pair distances is utilized. The particle pair distance represents the pairwise distance between two individual particles. Figure 13 displays the histogram illustrating the number of particle pair combinations against their respective particle pair distances for each case, i.e., the real case scenario. Additionally, the histogram corresponding to the ideal distribution is presented on the same graph. It is important to note that the ideal distribution histogram remains identical across all contactor cases with the same diameter. Cases 27, 30, and 33 exhibit the highest overlap with the ideal particle pair distance distribution, as visually evident in Fig. 12. The CDI quantifies the disparity between the particle pair distances in the real and ideal particle distributions. Furthermore, cases 25, 28, 31, and 34 show the highest disparity between the ideal and real distribution which can be seen in Figs. 12 and 13.

Figure 14 illustrates the positions of the last 10% of particles exiting the contactor in all 12 cases, with their masses taken into consideration. Here, the higher the particle size, the higher the mass of the particles. In all the cases, it can be seen that the particles with the highest masses tend to stay near the centerline of the contactor and are surrounded by smaller masses. This increases the cumulative masses near the core but reduces as we move toward the wall of the contactor. As gap size increases (top to bottom), i.e., $a$–$d$–$g$–$j$⁠, $b$–$e$–$h$–$k$⁠, and $c$–$f$–$i$–$l$⁠, more number of larger particles tend to align themselves near the core. This can be again explained as the effect of the injection speed of the particles, where higher mass shows larger inertia and if injected at a low velocity, they tend to stay in a low shear region and closer to the core of the contactor.

Figure 15 shows the distribution of mass of particles with radius used to calculate the mass-CDI which is analogous to the particle pair distance distribution in the case of CDI. The hemispherical domain is divided into strips radially and the cumulative mass in individual radial strips is calculated. Again, the ideal distribution is identical in all the cases as the exit diameter is the same. In Fig. 15, it can be seen that the ideal distribution requires minimum mass at the core and maximum mass close to the walls. It can be observed that for small diameters (case 25, case 28, case 31, and case 34), the particles are more clustered to the minimum radius which corresponds to the center of the exit plane. As the diameter increases (left to right, i.e., $a$–$b$–$c$⁠, $d$–$e$–$f$⁠, $g$–$h$–$i$⁠, and $j$–$k$–$l$⁠), the mass tends to move toward the wall (maximum radius on the $x$-axis) of the contactor indicating the spread of the total mass in the domain. Increased contactor diameter results in lower axial steam velocity upon entering, leading to extended residence time for particles within the domain. This prolonged duration allows particles to disperse more readily throughout the domain in instances where steam inlet velocity is lower, as observed in cases with higher inlet diameters (case 27, case 30, case 33, and case 36). Conversely, instances characterized by high axial velocity exhibit less dispersion due to the shorter residence time of particles, particularly notable in scenarios involving contactors with smaller diameters (case 25, case 28, case 31, and case 34).

It can also be noted that when the gap size increases (from top to bottom), i.e., $a$–$d$–$g$–$j$⁠, $b$–$e$–$h$–$k$⁠, and $c$–$f$–$i$–$l$⁠, the injection velocity decreases and more mass tends to move closer to the core of the contactor. An increase in the cement gap size results in a reduction of the inlet cement particle velocity in the radial direction. When these particles mix with the steam, characterized by low radial particle velocity, but higher axial steam velocity, their movement in the radial direction becomes impeded. In case 27, despite having a high inlet diameter resulting in low axial velocity, the presence of a low cement gap size results in high radial velocity after the injection of particles. Consequently, particles experience extended residence time, facilitating their dispersion throughout the domain. Conversely, in case 34, characterized by a low contactor diameter, but high cement inlet gap size, the high steam velocity is counteracted by low particle velocity. This results in insufficient residence time for particles within the domain, reducing their spread.

To summarize all the findings, CDI and mass-CDI are calculated for all 69 cases and are presented in Fig. 16. Cases such as 24, 22, 13, and 14 belong to the low CDI category, whereas cases such as 24, 40, and 48 have low mass-CDI values. Examples of high CDI values include cases 66, 46, and 58, whereas those for high mass-CDI values include cases 66, 58, and 62. It is useful to classify the cases into four categories, namely (a) low CDI and mass-CDI, (b) low CDI and high mass-CDI, (c) high CDI and low mass-CDI, and (d) high CDI and high mass-CDI cases. One such example is shown in Fig. 17. The spatial distribution of particles, with and without considering mass, is depicted in Fig. 17. Case 24 (Fig. 17(a)) is an example of low CDI and low mass-CDI. Here, particles are well spread in the domain, resulting in a low CDI value and cumulative mass is low near the core and increases radially, resulting in a low mass-CDI value. The well-spread particles and uniform increase in cumulative mass would potentially lead to higher evaporation efficiency due to increased heat transfer between the particles and the continuum phase. In case 54 (Fig. 17(b)), the low CDI is caused by the larger spread of particles, and high mass-CDI is evident as the mass is clumped near the core. Comparing case 24 with case 54 shows that, even though the space in the domain is occupied by particles evenly to some extent, resulting in low CDI, the bulk of the larger size particles resides near the core and not closer to the walls resulting in higher mass-CDI for case 54. Case 69 (Fig. 17(c)), on the other hand, demonstrates that the larger radius is occupied by fewer particles, resulting in high CDI. However, since the space close to the wall is occupied by a larger mass, the mass-CDI is still low. Case 58 (Fig. 17(d)) exemplifies a scenario of high CDI and high mass-CDI, where particles are not well spread in the domain, and mass is concentrated near the core. This results in a significant amount of unoccupied space or a large space occupied by small particles in the domain, leading to very low devolatilization efficiency.

Finally, the increase in average (of the final 10% of the contactor) polymer content is examined again similar to Fig. 11, except this time, in correlation with the calculated CDI and mass-CDI. Figure 18 illustrates the impact of CDI and mass-CDI values on the percentage increase in polymer content, which is directly related to the evaporation of cyclohexane content (volatiles) from the polymer mixture, for all 69 cases. Based on the increase in average polymer content, the cases have been classified and highlighted in various colors. The threshold values for the various evaporation groups are arbitrarily chosen and are as follows: $<13%$ is the low evaporation group, $>22%$ is the high evaporation group, and $13$–$22%$ is the group in the middle. In general, Fig. 18 shows that all the middle-high evaporation group cases lie in the low-CDI, low mass-CDI category, with CDI $≤$ 0.4 and mass-CDI $≤$ 1.5. On the other hand, all the low-evaporation group cases, except for 13, lie either in the high CDI or the high mass-CDI category, with either $CDI>0.4$ or mass-CDI > 1.5. The location of the 13 exceptions is an indication that there might be other factors that are leading to this discrepancy. For instance, both primary and secondary breakup processes of the cement are neglected through the DPM approach that only takes into account the evaporation here. Future studies on devolatilization in steam contactors can focus on the primary sheet breakup and the further secondary breakup of droplets of the polymer to obtain more accurate PSD.

Three-dimensional CFD calculations of multiphase phenomena in polymer devolatilization processes are carried out using the commercial CFD software, ansys fluent. The devolatilization process involves the high-speed flow of steam, treated as a continuous phase, and liquid cement made up of polymer and volatile cyclohexane, treated as a discrete phase, through fluent’s DPM approach. While previous studies [5,6,25,26] have mostly investigated the effect of control parameters such as steam flowrate and temperature on the efficiency of the devolatilization process, the distribution characteristics of the cement particles have not been focused on as a means to evaluate the efficiency of the steam contactor system. The current study attempts to fill this gap by developing quantifying metrics based on cement particle distribution inside the steam contactor, thereby providing insights for further optimizing the devolatilization process toward more efficient solvent removal systems.

The CFD model is verified using a comprehensive grid-convergence analysis conducted by evaluating GCI or grid-convergence index of three different grids. Analysis shows an asymptotically decreasing GCI for multiple variables at multiple locations, thereby justifying the use of the grid in all the calculations. Validation is carried out by comparing PSD to the sieve data from the industrial runs for 24 different cases. Reasonable comparisons were obtained between computations and sieve data with an RRMSE value of 13% and a correlation coefficient of 0.97. Eulerian contours of the continuous phase (steam) show that velocity increases and temperature decreases significantly as the steam impacts the cement injection due to the presence of the shock. The distribution of particles colored by the contours of their diameter and volatile mass fraction indicates that smaller particles move closer to the walls due to inertia and undergo more rapid evaporation.

The interaction between steam and cement within the contactor plays a pivotal role in determining the distribution of particles at the contactor’s exit. To understand the impact of particle position and mass distribution on the exit quality of the polymer, two metrics, namely CDI and mass-CDI, are introduced. They are used to assess the particle position and mass distribution compared to ideal scenarios, where maximum evaporation is possible. The ideal distribution for CDI involves particles spreading out uniformly, while for mass-CDI, an ideal distribution assumes uniform mass distribution among particles. The utility of these metrics is demonstrated through a presentation of 12 cases constituting different cement inlet gap sizes and contactor diameters. CDI for the various cases shows that particle clustering increases with cement inlet gap size and reduction in contactor diameter. Furthermore, mass-CDI also displays similar trends where larger particles (with larger mass) tend to move away from the core of the contactor, as gap size decreases and contactor diameter increases. Based on the behavior of CDI and mass-CDI, the best and worst cases are identified among all the conditions. Finally, these metrics are correlated with the most important deliverable in these devices: the volatile mass fraction at the exit of the contactor. A good majority of the cases, i.e., 56 of 69, showed that both CDI and mass-CDI were good indicators of the devolatilization efficiency.

The current study develops a detailed CFD model for studying the devolatilization in a steam contactor, presenting significant advancements in understanding thermal energy transport within the contactor. By correlating the metrics such as CDI and mass-CDI with the reduction in cyclohexane content in the contactor, the study demonstrates that the particle distribution characteristics can guide the optimization of the devolatilization process through a careful choice of the input conditions in polymer processing. In summary, the guidelines developed through this study can help the polymer industry in optimizing a steam contactor system in primarily two ways: (1) optimizing the amount of steam used to recover the polymer from the cement, which is both a cost and a sustainability issue and (2) optimizing the residual hydrocarbon left in the polymer following the solvent removal operation, which will have to be separated and captured in downstream operations. The model developed in this study has the potential to result in real, quantifiable variable cost savings for the thermoplastics industry. Such a verified and validated model can be utilized for future studies involving scaling exercises, for instance, higher flowrates for higher throughput, while preserving critical phenomena such as heat and mass transfer, shear, solvent removal, and polymer form. In addition, more accurate particle distributions near the inlet can be obtained in the future through detailed studies of the primary breakup of the polymer sheet.

The authors express their gratitude for the financial support extended by Ansys Software Private Limited through the Ansys Ph.D. Fellowship program.

There are no conflicts of interest.

The data that support the findings of this study are available from the corresponding author upon reasonable request.