CFD simulations of polymer devolatilization in steam contactors

ABSTRACT Removal of volatiles and other unwanted products in a polymer mix is a critical step in polymer manufacturing. This process serves many purposes, such as improvement in the physical and chemical properties of the polymer, elimination of odours, recovery of monomers and solvents, and fulfilment of environmental regulations. In this study, the mixture of polymer and solvent, or cement, is mixed with superheated steam which causes the unwanted solvent to evaporate. As the cement droplets lose solvent and dry out, they are ejected as cement ‘crumb’ from the mixing device. This process is modelled using computational fluid dynamics (CFD) to gain insight into the complex phenomena. The model simulates the heat transfer and phase changes associated with cement and steam, via 3D axisymmetric calculations. Using an Eulerian–Lagrangian approach, the superheated steam is modelled as the continuous phase and tracked in an Eulerian frame, while the cement droplets are treated using a Lagrangian tracking method, thus as a whole allowing us to track the particle sizes, temperatures, and solvent content. Furthermore, a parametric study is carried out to analyse the effect of initial polymer temperature on final polymer product. The simulation is carried out with three different initial polymer temperatures and the resulting solvent concentration, and the size of the cement particles between the three cases are compared. By increasing the cement operating temperature, the solvent concentration in the cement crumb is significantly lower, and the final cement crumb sizes shows a slight decrease, which indicates better production performance. Indeed, full 3D CFD calculations, presented to verify the axisymmetric assumption, show differences in the particle statistics, which could have implications for design considerations. Such studies can provide further insights into control parameters that can potentially enhance the efficiency of such a manufacturing process.


Introduction
This study focuses on a complex polymer manufacturing process that involves removing unwanted volatiles to refine the final polymer product. When making polymers in a reactor or other device, the polymer may contain some components such as unreacted monomers, solvents, water, and other by-products, which are referred to as volatiles (Albalak, 1996;Gestring & Mewes, 2002). These volatiles can negatively affect the quality of the final polymer product, and therefore volatile removal, or devolatilization, is an important procedure in polymer production. This devolatilization process is similar in certain regards to flash evaporation (McGreavy, 2012). In such problems, the levels of volatiles are in the tens of percentage points. So reductions in volatiles are measured as a percentage.
There are various methods for polymer devolatilization, but this study focuses on one particular method, which is as follows: the polymer is stripped from the solvent by mixing with superheated steam in a CONTACT Abhilash J. Chandy achandy@uakron.edu pressurized vessel, which is referred to as a steam contactor. The steam acts as a carrier to transport the evaporated solvent away from the contactor (Ravindranath & Mashelkar, 1988). Additionally, the steam is particularly stable in the sense that it can strip the polymer from the solvent without significant agglomeration (Albalak, 1996). This part of the devolatilization process lasts only for a few fractions of a second, before being ejected into a coagulator to continue mixing at high temperatures for several hours. As the mixture of polymer and solvent mixes at high temperatures with steam, the solvent begins to evaporate out of the cement, leaving behind a mixture with lower concentrations of solvent. The solvent vapour is then mixed in with the steam, and gets carried out of the contactor (Albalak, 1996).
Research has been carried out previously on polymer processing, specifically on devolatilization techniques like those presented in this study, although their focus was mainly on experiments (see Fleury, Witte, &Schildknecht, 2005, andGestring &Mewes, 2002). For instance, Fleury discussed the effects of using additional energy to drive out solvent via additional steam. With regard to modelling, one of the few studies was that of Yan, Luo, Lu, & Chen (2012), where gas-solid phase change in fluidized-bed polymer reactors was investigated in an Eulerian-Eulerian frame.
The current paper focuses only on the process that occurs in the steam contactor, and not on the coagulation process. This rapid process in the steam contactor involves injecting steam and polymer mixture, at high temperatures and velocities. This cement is a mixture of pure polymer plus known volatiles, and in this case the only volatile is a common solvent, cyclohexane. After the cement undergoes devolatilization in the contactor, the final result is a polymer 'crumb', or small polymer particle of non-uniform shape.
Due to the complexities associated with the multiphysics (fluid flow + heat transfer + multiphase) processes in such manufacturing systems, and also due to the high temperatures and pressures involved, laboratory experiments on such problems are almost impossible to create. On the other hand, computational fluid dynamics (CFD) simulations involving a consideration of most of the sub-processes can be performed at various operating conditions and parameters with relative ease. Hence, CFD simulations can be used as cost effective tools to model the devolatilization process in a steam contactor employed in polymer isolation, and operating parameters can be easily changed to optimize the process.
Several studies in the past have employed various modelling approaches for similar multiphase problems. In spray drying technologies in the pharmaceutical industry, CFD simulations are widely used to study the phase change and evaporation of liquid particles in a continuous gaseous phase. For instance, the focus of Langrish and Fletcher (2001) was on spray drying technologies of foodstuff, where particles were tracked in a Lagrangian manner, and drag effects and heat and mass exchanges with their surroundings were included. Pimentel et al. (2006) discussed the methodology used in randomizing initial particle size in a similar fashion as used here. The use of a Beta distribution of initial droplet sizes was found to be appropriate to use in their scenario. Later, Woo et al. (2008) discussed different spray drying methods using an Eulerian-Lagrangian modelling approach. Dobry et al. (2009) emphasized the importance of mass and heat transfer in studying the drying of active soluble ingredients, and Mezhericher, Levy, and Borde (2009) used a CFD model to show droplet drying using an Eulerian-Lagrangian method. Other studies include those of Gjesing, Hattel, and Fritsching (2009), Ruiz-López, Ruiz-Espinosa, Luna-Guevara, and García-Alvarado (2011), and Sazhin (2006). The work performed by Shalaby, Wozniak, and Wozniak (2008) used numerical simulations to calculate particle cyclone separations by solving for particle trajectories in a Lagrangian approach and continuous phase flow to study cyclone efficiency. More recently, Bonifacio, Maghirang, and Glasgow (2014) discussed methods in numerical simulations of particulate transport and dispersion. The study found the use of the k-Reynolds-averaged Navier-Stokes (RANS) turbulence model appropriate for gas-solid particle trajectory simulations and determining resulting particle concentration profiles. Another study regarding particle dispersion using CFD techniques was performed by Lin and Shen (2014), using RANS modelling to study the turbulent effects in particle flows. Vittorio Messa and Malavasi (2014) did additional research in predicting the multiphase particle distribution in pipe flow. Using a k-RANS turbulence model, the solid phase distribution was calculated. Although these past studies are from different applications, they use simulation techniques that are similar to those used in the current study.

Scope of work
The primary objective of the current study is to model the multiphase processes involved in the stripping of a polymer mixture using steam. Three-dimensional (3D) CFD calculations of the problem are conducted using the commercial software ANSYS R Fluent R . The problem itself involves the cement and steam being injected into the contactor separately, but in close proximity, at different temperatures and pressures, such that as they are injected they begin to mix almost immediately. As the materials mix, the temperature of the steam begins to heat the cement mixture, causing solvent evaporation, thereby stripping the mixture of the solvent. Several assumptions were made when creating this model. The cement mixture is initially injected as a randomized set of spherical particles, while in reality they are not spherical due to the cement shearing process upon injection. Additionally, particle collisions are neglected in the current model, in order to keep the latter relatively simple.
To capture all of the complexities of the polymer devolatilization process, and to solve accurately for the heat and mass transfer of the particles while also simulating the particle trajectories, an Eulerian-Lagrangian (E-L) method is used to solve the particle trajectories. This E-L approach also allows use of a multicomponent species transport method to solve for the mass transfer of the cyclohexane solvent evaporating into vapour and being carried away with the steam vapour. The method also solves for the heat transfer occurring as the particles heat up and lose solvent. This 'flashing' process occurs under high pressures and temperatures. Once the initial model is verified, parametric studies of the operating conditions can be performed.

Conservation equations
When modelling fluid flow, the dimensionless parameter Reynolds number (Re) is used to characterize the flow. Re is defined as where ρ is the density, V is the fluid velocity, L is the characteristic length, and μ is the kinematic viscosity of the fluid. The Re calculated in this study was 334,562 and hence the flow would be turbulent. Turbulent flows are characterized by fluctuating velocity fields, which further cause other transported quantities to fluctuate as well, at high frequency. So the instantaneous conservation equations can be time-averaged to remove the small scales, resulting in a modified set of equations that are computationally less expensive to solve. In addition, due to the high steam inlet velocities in the problem studied here, the flow is compressible, with the Mach number, Ma = V/c being around 0.6. In order to account for the density and temperature fluctuations, the compressible Navier-Stokes equations undergo a density-weighted time-averaging or a Favre-averaging procedure given by where the overbear denotes the conventional Reynolds average. So the Favre-averaged Navier-Stokes (FANS) conservation equations are given by In Equations (3), (4), and (5), ' (· · · )' represents the fluctuating quantity (Burns, Frank, Hamill, & Shi, 2006;Walters & Cokljat, 2008), ρ is the density, u is the velocity, P is the pressure, t ji is the viscous stress tensor, e and h are the specific internal energy and specific enthalpy, respectively, and q L j is the laminar mean heat flux. Equation (6) is the Favre-averaged ideal gas equation used to model the continuous steam phase. Additionally, the compressibility factor, Z, is the ratio of the actual volume to the ideal volume of real gas, and correlates how much a gas can deviate from perfect behaviour (Heidaryan, Moghadasi, & Rahimi, 2010). To ensure an ideal gas model is an appropriate assumption to use in this case, the compressibility factor, which is unity for a perfect gas, is calculated for the superheated steam at the operating pressure and temperature given in this case. The compressibility factor is calculated at Z = 0.91, and is within good enough proximity to Z = 1, and hence modelling the steam as an ideal gas is appropriate. The viscous stress tensor is written as where s ij is the strain rate tensor, ζ is the second viscosity coefficient, and δ ij is the Kronecker delta. The strain rate tensor is written as With regard to the turbulence model, the re-normalization group (RNG) k-model is employed in the simulations presented here. The RNG model, developed using methods by Yakhot and Orszag (1986), renormalizes the Navier-Stokes equations to focus on smaller scales of motion over the standard k-, thus leading to improved simulation performance regarding highly strained flows. This model is employed also due to its improved accuracy for capturing the effect of swirl on turbulence, and flows with strong streamline curvatures, vortices, and rotation (ANSYS R Fluent R , 2012; Hou & Zou, 2005;Jawarneh, Tlilan, Al-Shyyab, & Ababneh, 2008;Papageorgakis & Assanis, 1999). This simulation shows regions of highly strained flow near the inlets, thus making the choice of this model appropriate for this simulation. The RNG model, compared to the standard k-model, has an additional term that improves the accuracy for rapidly strained flows. The compressible RNG k-transport equations are written as where G k represents the generation of turbulent kinetic energy due to the mean velocity gradients, G b is the generation of turbulent kinetic energy due to buoyancy, Y M is the contribution of the fluctuation dilation in compressible turbulence to the overall dissipation rate, α k and α are the inverse effective Prandtl number for k and , respectively, and S k and S are source terms to account for the coupling between the continuous and particle phases.

Species transport
To model the mixing and transport of concentrations of different chemical species in the problem studied here, a species transport model is used. When solving the conservation equations for chemical species, the local mass fraction of each species, Y i , is calculated from the solution of a convection-diffusion equation for the ith species (ANSYS R Fluent R , 2012). The conservation equation is described in the general form where J n is the diffusion flux of species n and is given by Sc t is the turbulent Schmidt number, which is set to 0.7, and was found to be appropriate for the circumstances in this simulation. For turbulent flows, the mass diffusion is influenced mostly by the turbulent transport. This is determined by the turbulent Schmidt number, and measures the ratio of momentum and mass diffusivity (Goldberg, Palaniswamy, Batten, & Gupta, 2010;Reynolds, 1975;Tominaga & Stathopoulos, 2007;Yeung, Xu, & Sreenivasan, 2002). The turbulent Schmidt number is relatively insensitive to the molecular fluid properties, thus there is little reason to alter the default value of Sc t = 0.7 for any species (ANSYS R Fluent R , 2012). D n,m is the mass diffusion coefficient for species i in the mixture, μ t is the turbulent viscosity, and D T is the thermal (Soret) diffusion coefficient. Thermal diffusion is enabled here, which causes heavy molecules to diffuse less rapidly, and light molecules to diffuse more rapidly. The thermal (Soret) diffusion coefficient is computed by applying kinetic theory using the empirically-based expression where M w,n is the molecular weight of species n, Y n is the mass fraction of species n, and X n is the mole fraction of species n (ANSYS R Fluent R , 2012). The diffusion flux term on the right-hand side of Equation (8) includes the enthalpy, which accounts for the treatment of species transport in compressible flow.

Two-way coupling
The coupling between the continuous phase and particle phase is taken into account here as a particle's movement and trajectory are influenced by the gas phase. In this instance, this coupling is achieved by calculating source terms in the flow's turbulent equations. A particle's movement is dictated by the turbulent flow of the continuous gas phase (Berlemont, Grancher, & Gouesbet, 1995;Gouesbet & Berlemont, 1999;Kolaitis & Founti, 2006), and thus within Equations (9) and (10) there are two source terms, S k and S , which are calculated in the simulation. Additionally, when modelling the disperse particles, the degree of coupling between the continuous and discrete phase is dictated by the Stokes number (St). The Stokes number is the ratio of the relaxation time of particles to hydrodynamic time. Essentially, St relates a particle's inertia to its viscous drag by comparing the particle's aerodynamic response time to the characteristic flow time, which in turn affects the particle's trajectory. Thus, this dimensionless parameter measures the coupling between the continuous gas and the liquid droplet phases, and the effects of the continuous phase on the particle's trajectory. The Stokes number is defined as where d p is the average particle diameter, L is the hydrodynamic length scale (the size of the computational domain), R is defined as follows: where ρ f is the carrier fluid density, and ρ p is the particle density. When the Stokes number is small, i.e. St 1, the particle's response time is much less than the characteristic time in the flow field. So particles have sufficient time to respond to flow velocity fluctuations, thus interacting with the continuous phase (Apte, Gorokhovski, & Moin, 2003;Crowe, Sommerfeld, & Tsuji, 1998;Raju & Mejburg, 1995). When the Stokes number is large, i.e. St 1, the particles have little or no time to respond to fluid flow changes, and therefore are only barely affected by the fluid flow. This means the particles will pass through the continuous phase as opposed to interacting with it, and hence the trajectories of the particles remain relatively unchanged. The case modelled here has a relatively high Stokes number, specifically, St = 19. This Stokes number is more of a large-scale or global parameter, determined based on the average particle diameter in the practical contactor (≈ 1000 µm).

Problem description
The current work involves modelling a polymer devolatilization process in a steam contactor whereby a distinct polymer is produced. The steam contactor (see Figure 1(a)) injects superheated steam through two annular slits, while the cement mixture (made up of polymer and cyclohexane) is injected in an annular slit alongside the steam (see Figures 1(b) and 1(c)). The superheated steam comes in contact and mixes with the cement. The higher momentum of the steam shears the cement and fragments it into smaller droplets, while the high temperature of the steam flashes the cement. In other words the solvent evaporates and polymer particles, or crumb, are left behind. The polymer crumb continues mixing with the steam through the contactor, until it reaches the outlet (see Figure 1(e)) and is exited into the coagulator of the polymer isolation process. The scope of the current study is limited to modelling the flow and heat transfer process in the contactor and does not involve the coagulator processes.
Presented here is a parametric study to determine the effects of the cement injection temperature on the final polymer crumb. The objective is to optimize the solvent evaporation process by varying the initial temperature of the cement mixture. To carry out this study, three different simulations were carried out by only varying the cement temperature. Each simulation shared the same mesh and boundary conditions. The cement injection temperatures for each case are as follows.
The temperatures where chosen such that they stayed within typical operating parameters seen in this manufacturing problem. The standard operating temperature is around 443 K, and the goal of the study is to obtain similar results specifically with regard to the cyclohexane mass fraction and particle size, using cement at a lower initial preheated temperature, which could save operating costs. Temperatures below 400 K were not considered here mainly because values that low would result in a product with a significantly worse quality, especially with respect to the cyclohexane content.

Geometry details
The full three-dimensional domain for the steam contactor can be seen in Figure 1(a). The total length of the contactor is 713 mm, with a maximum outer diameter of 307 mm. The geometry and mesh for this simulation were generated using ANSYS R Workbench 14.5 and used for all three cases. Due to the symmetry along the axial plane (see Figure 1(d)), to conserve computational cost, only half of the original 3D domain was simulated, with the symmetry boundary shown in Figure 1.

Operating conditions
The only parameter varied in the different cases presented here was the temperature, and the remaining operating conditions remained consistent between the  cases. The operating conditions used in this simulation are shown in Table 1. It is also assumed that the fluid leaving the contactor exits into standard atmospheric conditions of 300 K and 101 kPa. Table 2 shows the conditions of the cement particles as they are injected.

Computational details
Several models within the commercial software ANSYS R Fluent R were enabled and coupled to simulate this particular multiphase problem. When modelling multiphase flows with relatively small dispersed particles, the DPM is the optimal choice, which solves for the continuous steam phase in the steady state, while tracking the particles in a transient manner.

Multiphase particle flow with Lagrangian approach
The multiphase problem here is a complex one that involves a liquid-gas phase change, and three different species. The steam is modelled as a continuous gas phase, and remains so throughout the entirety of the simulation. The cement was composed of two different components/species, polymer and cyclohexane, premixed initially, with the latter being a liquid. This simulation uses a species transport mixture model to combine the polymer and cyclohexane into one fluid, cement, and solves for the composition of the new fluid as it changes in the course of the simulation. In this instance, the polymer is a Newtonian fluid, and its viscosity is only a function of temperature, with a piecewise linear function fit, whose values were supplied by the manufacturer. Furthermore, as the liquid cement and gaseous steam mix inside the contactor, the cyclohexane vapour resulting from the evaporation of cyclohexane liquid within the cement mixes with the steam and is vented outside the contactor. This entire process results in cement with a higher composition of polymer and a lesser amount of solvent than those with which it started out. So it should be noted that polymer refinement, not polymer generation, is the physical process modelled in this study. When modelling a secondary phase, there are two general models used to simulate secondary phase flow: Eulerian-Eulerian (E-E) and Eulerian-Lagrangian (E-L) (Kerdouss, Bannari, & Proulx, 2005;Mezhericher et al., 2009;Soltani, 2012). The fluid phase is treated as a continuum by solving the Navier-Stokes equations, while the dispersed phase is solved by tracking a large number of particles or droplets through the calculated flow field. The dispersed phase can exchange momentum, mass, and energy with the fluid phase (Gouesbet & Berlemont, 1999). To use the E-L approach, the secondary phase must be dilute enough. More specifically, the secondary phase volume fraction must be below 10%. In, this case, the volume fraction of the secondary phase is calculated to be approximately 4%. Since the secondary phase in this simulation is relatively dilute and the particles are much smaller, the collisions are not as significant in the flow, thus making the E-L approach a suitable choice, in which the particle-particle interactions can be neglected and the dispersed second phase occupies a low volume fraction (Beishuizen, 2007). The particle trajectories are computed individually at specified intervals during the fluid phase calculation. Because of the requirements of the E-L approach, DPM is used to simulate the flow of the particles. A DPM adopts the E-L approach, uses a set of governing equations (as presented above) to model the continuous phase (steam) in an Eulerian frame, and uses a set of equations to model the discrete phase (polymer and solvent particles) in a Lagrangian frame by tracking the particles in the flow field.

Discrete-phase model
To predict the trajectory of the discrete-phase particles, the model integrates the force balance on the particle, which is written in a Lagrangian reference frame. This force balance equates the particle inertia with the forces acting on the particle, and can be written as where F D ( u − u p ) is the drag force per unit particle mass, given by and u is the continuous phase velocity, u p is the particle velocity, μ is the dynamic viscosity of the fluid, ρ is the density, ρ p is the density of the particle, and d p is the particle diameter.
As the cement particles are injected into the contactor, they are prescribed a certain diameter size, specified by a given size range ( Table 2). The particle size distribution (PSD) though is determined by a beta density function (Devroye, 1986). The beta density function is defined as where D is diameter, a and b are lower and upper bounds ranging from zero to one, p > 0 and q > 0 are shape parameters, and C is the normalizing constant derived to A user-defined function (UDF) coded in C is coupled with ANSYS R Fluent R to model the beta density function distribution. The UDF was programmed to generate random numbers, fit the random numbers to a beta distribution and to set the generated data to fall within the minimum and maximum values for droplet particle size.
For the simulations here a = 0, b = 1, p = 2.0 and q = 1.3 (Devroye, 1986), and are chosen such that the beta density function matches with the inlet particle distribution provided by the manufacturer.
In addition, the simulations use a high-Mach number drag law to account for a particle Mach number greater than 0.4 and particle Reynolds number greater than 20 (Morsi & Alexander, 1972). Using the DPM, reactions of the particles and how they affect the continuous phase can be examined using one of the heat and mass transfer laws in Fluent R (Dobry et al., 2009). The particles in this case may be composed of either pure liquid cyclohexane, pure polymer, or a combination of both. Since this study focuses on particles containing two species, the multicomponent particle heat and mass transfer conditions are used. The multicomponent droplet vaporization rate is calculated as the sum of the vaporization rates of the individual components. Furthermore, to track the cement particles as they flow through the domain, the stochastic DRW model in Fluent R is used (Ajilkumar, Sundararajan, & Shet, 2008). This method uses the velocity fluctuations caused by turbulence to predict the particle trajectory.

Mesh details
A reasonably fine mesh is necessary to calculate flow and heat transfer accurately and also track the particle trajectories. The simulations presented here use a mesh with 518,070 elements and 327,512 nodes (see Figure 2). A hexahedral dominant mesh is preferred in this simulation due to better performance in terms of convergence and accuracy, but the unique shape of the outlet results in a region of tetrahedral elements near the outlet.

Boundary conditions
For both steam inlets, a mass flow inlet boundary condition is used, while the outlet uses a pressure outlet boundary condition set to ambient pressure and a temperature of 273 K. All wall boundary conditions use an adiabatic no-slip condition. The continuous steam phase is solved for in a steady-state simulation to convergence, before injecting particles into the system. The particles are injected using a DPM paired with the beta distribution to randomize the particle diameter sizes as mentioned earlier. A symmetry along the axial plane is utilized in order to reduce the computational effort. This boundary condition is found appropriate here due to the nature of the flow within the contactor. A justification for this boundary condition is also provided below (see Section 6.4). With regards to simulating particle flow in the steam contactor, all surfaces except for the outlet use a reflection type boundary condition. The reflection boundary condition redirects particles away from any boundary other than the outlet, to enable particles to exit the contactor only though the outlet. The reflection in this case only accounts for a perfect elastic collision, such that no momentum is lost after colliding with any surfaces.

Solution methods
A pressure-based solver is used for the three simulations, with a semi-implicit method for pressure-linked equations (SIMPLE) pressure-velocity coupling scheme, the PRESTO! (PREssure STaggering Option) pressure solver scheme, and second-order upwind schemes for momentum, density, energy, and turbulence. The PRESTO! scheme is used due to its robust solving method in which the scheme uses the discrete continuity balance for a staggered control volume about the face to compute the staggered (i.e. face) pressure. It is a better choice than the standard pressure solver regarding multiphase flows, flows with curved domains, and flows with large pressure gradients (ANSYS R Fluent R , 2012; Huang, Hong, & Kim, 2012;Kaya & Karagoz, 2008).
Additionally, the DPM uses the accuracy control method, which enables solutions of equations of motion within a specific tolerance. This is achieved by computing the error of integration steps and reducing integration steps if the error is too large. When accuracy control is activated, the step-length factor will only be used to estimate the time step of the first integration step, and in all subsequent integration steps the particle integration time step is adapted to achieve the specific tolerance (ANSYS R Fluent R , 2012). Since we are interested in the composition of the particles that reach the outlet, the outlet uses an escape boundary condition, and the statistics of escaped particles can be assessed and compared between cases.

DPM time-stepping
A DPM injection time step size of 10 −3 seconds was found to be appropriate here. The time step used when solving for the particle trajectories using the DPM is needed purely for numerical purposes, specifically for advancing the discrete-phase particles implemented in these simulations. The trajectories are calculated in a Lagrangian reference frame, and thus the particles are computed via unsteady integration (ANSYS R Fluent R , 2012). To ensure this time step size does not affect the final outcome, a time sensitivity study was performed with the same operating and boundary conditions, with a cement injection temperature of 443 K, and changing the DPM time step size to 10 −4 seconds, and there were no changes in the results or computational cost.
Particles are assigned a maximum number of 42,000 time steps with a step-length factor of 5 before their path is terminated to ensure that, if a particle gets stuck in the contactor in a re-circulation zone or other obstruction, it is terminated after a certain number of time steps. The step-length factor is the number of time steps required to transverse the current continuous phase control volume, or computational cell. It controls the time step size used to integrate the equations of motion for the particles and computes the time step in terms of the number of computational cells. A step-length factor of five is found appropriate since each cell is approximately 0.5 mm in length, and thus the equations of motion for each particle are solved every 0.1 mm, where the average particle size is 0.08 mm. The maximum number of time steps is also found to be sufficient such that no particles are terminated too early, and also not too large, such that it does not result in a significantly high computational cost.

Experimental validation
To ensure the validity of this model, the particle size distribution is compared with some experimental data measured in the field. Here, the simulation is carried out for a similar steam contactor, and compared to the experimentally gathered results. Although the contactor geometry is different, the physical process remains the same, and thus the model is unchanged. A sample of approximately 100,000 escaped particles from the simulation is analysed and compared with the provided data. The experimental data obtained are averaged over several test batches. Additionally, since collecting cement particles at the exact location of the contactor outlet poses several challenges, the cement particle sizes are physically analysed after they have passed through the contactor and coagulator phases, and thus some differences are expected. Figure 3 shows a histogram plot of the simulation results of cement particle sizes at the outlet region plotted with the averages of the experimental data gathered. The red line in this figure represents the experimental data, while the bars represent the simulation results. Again, since the data are analysed at a different location, values do not match up precisely, but the trend and particle range sizes and averages match appropriately enough to verify the model is sufficiently good at predicting the end results.

Grid-independence
A comparison between the results of Cases I, II, and III requires one to validate the entire numerical approach, including the mesh, solution techniques, and the type of boundary conditions used for the current model. For Figure 3. Histogram plot of simulation data of escaped particle diameter sizes compared to averaged values of experimentally gathered cement particle sizes after the physical process.
this purpose, a grid-convergence study involving three simulations with different meshes, but identical operating conditions, is conducted to ensure that the mesh is sufficient to resolve all the flow and temperature gradients. For simplicity, each simulation is carried out using a steady-state approach and employing only steam. By comparing ranges of temperature and velocity values in the contactor near the inlet when injection speeds are at a maximum, any large differences can be noted and further tested to determine if the solution is indeed mesh independent. So a medium mesh, with 518,070 elements, for which results are presented later, is compared to a finer mesh with 1,378,718 elements (i.e. approximately 2.5 times more cells than the original case) and a coarser mesh with 231,882 elements (i.e. approximately half the cells of the original case). The velocity and temperature values are measured along an axial line on the symmetry plane, shown in Figure 4(a), which corresponds to a position range of 545 to 645 mm. The velocity and temperature values are plotted for all three meshes in Figures 4(b) and 4(c). The highest velocity difference in the contactor between the coarse mesh and the original mesh is 16%, and the highest difference between the original mesh and the fine mesh is 2%. The corresponding value of temperature between the coarse and original mesh is 13%, and between the original mesh and the fine mesh is 2%. These results show that the original mesh is sufficiently refined to yield almost mesh-independent results.

Particle flow results
The simulations are set up to inject a large quantity of particles, or a parcel which is based on size and mass of the particles, into the system every four iterations. As the first parcel of particles is injected, it comes in contact with the steady flowing steam, causing instabilities in the flow. As the simulation progresses, the fluxes caused by the initial injection of particles begin to dissipate, and the flow begins to regain stability. To ensure quality data is gathered, the simulation is carried out until all of the initial fluxes dissipate and that enough escaped particles, approximately 500,000 particles, are tracked to provide sufficient results for statistical analysis.
The choice of injecting a parcel of particles every four continuous phase iterations is based on trial and error. Too few iterations between injections will be more difficult to converge the solution, but too many is computationally expensive. The time step is constant in this study, and is set to 0.1 ms. Since the particle trajectories are computed in a Lagrangian reference (transient) frame and the continuous phase in an Eulerian reference (steady-state) frame, a time step is necessary to solve the particle trajectories. The DPM in ANSYS Fluent R uses an unsteady integration method to compute the particle paths, which requires a time step purely for computational efforts to solve the stepwise integration over discrete time steps.
The number of particles per parcel is determined by where NP is the number of particles per parcel,ṁ s is the cement mass flow rate (which is provided), t is the time step, and m p is the particle mass, which is determined by the size and density. Given that the particles are injected at randomly selected diameters, the number of particles varies, but is approximately 4500 particles per parcel. To ensure that the simulation is converged and has tracked enough particles, the total number of particles in the simulation is also considered. The cement mass flux is monitored at the inlet and outlet, and when they match it is determined that the simulation has reached a steady state. At that point, a sample of approximately 100,000 particles is taken into consideration, and typically at this point a total of about one million particles has passed through the domain. By setting the outlet as an escape boundary for the particles, the particles that are ejected from the contactor are filtered from the particles that remain inside, and several parameters are qualitatively compared between the different cases. The variables that provide the most insight into this parametric study are cement particle material composition, temperature, total time it takes the particles to travel from the inlet boundary to the outlet boundary, or residence time, cement volume fraction, and particle diameter size. Figure 5 shows velocity contour plots for each case as well as the velocity plotted along the axial line on the symmetry plane at the end of the entire simulation. Note that the overall velocity increases as cement injection temperature increases. On observing the differences between Figure 5(a), Figure 5(b), and 5(c) (which correspond to Cases I, II, and III, respectively), the velocity near the inlet decreases at a faster rate for lower temperatures. Near the outlet, the velocity is much higher when inlet cement temperature is higher. This correlates higher cement temperatures with higher velocities, and when velocities are higher near this region, more mixing could occur, which is beneficial to this problem. This is worth noting since the higher temperature cases performed better, thus leading to the conclusion that higher velocity could contributed to better overall performance.
To get a better look at how the cement flows through the contactor, the cement volume fraction is plotted along the symmetry plane. The cement volume fraction shows the areas containing the highest concentrations of cement particles. If the temperature is not high enough to evaporate the solvent properly, then it can affect the contactor design. Shown in Figure 6(a), 6(b), and 6(c) are the cement volume fractions for Cases I, II, and III, respectively. These figures basically show where most of the cement particles are located inside the contactor during operation. It is important to look at this aspect, for instance, if there is a large deposit of cement in a section with very low velocity, from which it can be concluded that there is a re-circulation zone. Figure 6(d) shows a comparison of the cement volume fraction for each case, and it is clear that lower initial cement temperatures cause a higher concentration of cement in the contactor.
In each figure, the higher volume fraction values designate areas with large amounts of cement, while low volume fraction values show the areas that indicate relatively low amounts of cement. Figure 6(a) shows a region near the outlet with a high amount of cement, whereas in Figure 6(b) and 6(c) the concentration of cement in this region is lower, despite both cases showing a similar pattern of high cement concentrations stuck between the outlet and outer wall, where the cement accumulates. This high amount of cement collects due to the shape of the contactor. There is a region between the end wall and the outlet that appears to be an area in the contactor  that the steam flow does not reach properly. Between the three cases, as the temperature increases, the cyclohexane evaporates, causing less cement to remain in the contactor.
One of the most significant parameters in this study is the particle temperature, which varies between the three cases due to an injection temperature difference of 20 K. Because of this initial temperature change, the particles that escape are expected to show significant temperature differences. Figure 7(a), 7(b), and 7(c) show the particle temperatures inside the contactor. As we see for Case I, Figure 7(a) shows the particle temperature inside the main body of the contactor near the outlet region as approximately 360 K, and the average temperature of the escaped particles is 355.8 K. The temperature of the cement particles decreases as they travel through the contactor, due to the heat loss that occurs when the liquid cyclohexane flashes and evaporates from the particle. For Case II, the particle temperature inside the contactor is approximately 370 K, and the average temperature of the escaped particles is 365.9 K (Figure 7(b)). Compared to Case I, the temperature, on average, increased about 10 K throughout the contactor. Case III shows that the approximate particle temperature is 380 K, while the escaped particles have an average temperature of 374.9 K (Figure 7(c)). So while the initial injection temperatures of the three cases vary by 20 K, the difference between the average temperatures inside the contactor is only around 10 K, which is a result of the mixing between the superheated steam (which is at a constant temperature of 446 K) and the cement.
Each simulation is performed with an initial particle composition of 86% cyclohexane and 14% polymer, but with varying initial temperatures. By monitoring the particle's cyclohexane mass fraction, the amount of this solvent that evaporates can be tracked. As the percentage of cyclohexane drops in the particle, it shows that the temperature inside the contactor is high enough to cause the cyclohexane to evaporate into vapour, which eventually mixes with the steam. This leaves a particle with a lower concentration of cyclohexane and a higher concentration of polymer, which is ultimately the purpose of this method of processing the cement. In other words, as the temperature of the particles increases when they come into contact with the higher temperature steam, the cyclohexane is closer to its boiling point. The boiling Figure 7. (a), (b), and (c) Instantaneous view of cement particles coloured by particle temperature and sized by diameter for Case I, II, and III, respectively. point is 353 K at atmospheric pressure, and under the high pressures inside the contactor, this can be as high as 415 K. As the particle temperature increases, the amount of cyclohexane that evaporates also increases. Additionally, with low initial cement temperatures, the particle's temperature is much further away from its boiling point, and hence more superheated steam is needed for the same level of evaporation. Thus it is important to track the cyclohexane mass fraction per particle in the contactor, and more importantly, in the escaped particles. Figures 8(a), 8(b), and 8(c) show the content of liquid cyclohexane per particle inside the contactor for Cases I, II, and III, respectively. Figure 8(a) shows that the particles contain approximately 45% cyclohexane by mass fraction of the particles located in the main body of the contactor. As the particles escape, the average value of cyclohexane remaining in each particle is 44.8%. As the cement injection temperature is increased by 20 K, as in Case II, the cyclohexane content per particle is approximately 40% in the contactor, and on average the escaped particles contain 39.4% (shown in Figure 8(b)). By increasing the initial temperature, the resulting cyclohexane content has decreased by approximately 5%. By further increasing the cement temperature by 20 K as it is injected, the content of cyclohexane per particle decreases by approximately another 5% in the contactor, and on average the escaped particles contain 34.2% cyclohexane (shown in Figure 8(c)).
To get a better look at the particle sizes, temperatures, and composition between the three cases, histogram plots are shown for each case. Figure 9(a) shows the number of particles (as a percentage) at each particle size for all three cases. The overall diameter range stays approximately within 200 − 2500 µm. However, on averaging these values, Case I has larger particles than Case II, followed by Case III. Therefore, the higher temperature cases produce smaller particles. These values can also be seen in Table 3. Figure 9(b) shows the temperature distribution of escaped particles. It is clear from this figure that the lower the initial cement temperature, the lower is the corresponding escaped particle temperature. Figure 9(c) shows the escaped particle cyclohexane mass fraction distribution between the three cases. Obviously, for Case I, since the temperature is the lowest, the cyclohexane content is higher, and as the temperature increases, the cyclohexane content significantly drops. It is apparent from this figure that there is a broader cyclohexane mass fraction range for Cases II and III, while Case I has a narrower range. Both the escaped particle temperature and cyclohexane information are represented in Table 3. , and (c) Instantaneous view of cement particles coloured by cyclohexane mass fraction per particle and sized by particle diameter for Case I, II, and III, respectively. Figure 9. (a) Cement particle size distribution (by diameter in µm) of escaped particles in outlet region (b) particle temperature distribution of escaped particles (c) mass fraction of cyclohexane per particle between Cases I, II, and III. As the particles increase in temperature and more cyclohexane evaporates, the resulting particles are expected to be smaller in size. Table 3 shows a summary of the escaped particle sizes, including minima, maxima, and averages. Note that the residence time in the table indicates the total time that the cement particles took to travel from the inflow boundary to the outflow boundary. From this table, the average particle diameter sizes do not follow a specific pattern, but the size range between the cases shows a decreasing size trend with increasing temperature. The minimum particle sizes vary by 5.1%, while the maximum particle sizes vary by 5.0%.

The 3D axisymmetric assumption (versus full 3D)
There have been previous CFD studies as well, demonstrating the need to carry out full 3D simulations instead of resorting to the axisymmetric assumption. Alexeenko, Levin, Gimelshein, Collins, and Reed (2002) analysed micro nozzle performance and showed that the axisymmetric assumption could lead to significant design errors. In the study by Schlüter (2004), the importance of axisymmetry to large eddy simulation (LES) calculations was assessed by flow modelled with and without swirl over an axisymmetric expansion. On the other hand, Susan-Resiga, Muntean, Tanasa, and Bosioc (2009) assessed the validity of axisymmetric swirling flow simulations, when in reality the flow was unsteady 3D. It was concluded that, for strongly developing helical vortices, the 2D axisymmetric results were in reasonable agreement with 3D time averaged ones, and hence could reliably be used to obtain the base swirling flow for stability analysis.
Since the flow in this study is complex, including a continuous and particulate phase, the justification of the 3D axisymmetric boundary condition is important. There are limitations to using a symmetric boundary condition, especially in simulations involving cylindrically shaped domains, including rotating and swirling flows. In addition, the flow is turbulent. So although the geometry is axisymmetric, the turbulence is not. The axisymmetric assumption means there is no turbulent transport crosswise from the axis. Furthermore, many unsteady non-axisymmetric flow instabilities occur, especially in turbulent flows. Obviously the axisymmetric assumption would result in errors related to these instabilities, and typically, full 3D calculations are necessary. In order to verify this aspect in the simulations here, one fully 3D calculation was carried out for an inlet cement temperature of 423 K, and results were compared to the 3D axisymmetric calculations. This new full 3D calculation utilizes a mesh consisting of approximately 3.5 million elements for flow and temperature solutions, and 1.5 million particles for statistical analysis. With all the other numerical parameters set to be identical with those in the 3D axisymmetric calculation, the full 3D simulation increased the CPU time by a factor of three. Figure 10 shows the profiles of velocity and cement volume fraction comparisons between the axisymmetric and full 3D calculations, plotted along the axial line on the symmetry plane location at the end of the entire simulation. It is interesting to note that there are some significant differences between the velocity and cement volume fraction between the two calculations. With regard to velocity, the axisymmetric case overpredicts the velocity magnitude up to a factor of 1.4 (at the axial position of 0.36 m), and also overpredicts the cement volume fraction almost consistently throughout the symmetry plane location, with the maximum being a factor of three near the outlet. This is an indication of the lower rates of evaporation predicted by the 3D axisymmetric calculations Figure 11. (a) and (b) Instantaneous view of cement particles coloured by (a) particle temperature and (b) cyclohexane mass fraction per particle, and sized by diameter, for fully 3D calculations at cement inlet temperatures of 423 K. and will have implications for the drying of the polymer particles. Further evidence of these aspects is provided through particle distributions. Figure 11 shows the distribution of particles in the contactor coloured by temperature and cyclohexane mass fraction. Qualitatively the temperature and cyclohexane mass faction results look similar to the 3D axisymmetric calculations (Figures 7 and 8, respectively). For instance, the average particle temperature differs by 10 K between the 3D axisymmetric and full 3D calculations, as seen in the escaped particle statistics summarized in Table 4. The average particle diameter and cyclohexane mass fraction between the two calculations differ by 25 µm and 2%, respectively. However, one will be able to provide a better and more quantitative analysis of the particle statistics by assessing the histograms of the particles in these calculations, and more importantly their deviations from the 3D axisymmetric calculations, over the entire range of particle statistics. Figure 12 shows the distribution of particles in terms of particle diameter, temperature, and mass fraction of cyclohexane for both the full 3D and axisymmetric 3D cases. The axisymmetric 3D case corresponds to Case II in Figure 9. Figure 12(a) shows the particle diameter distribution and it is clear that the overall diameter range moves to the left for the fully 3D case, indicating smaller particle sizes. This is also evident in the particle temperature distribution in Figure 12(b), which shows a higher temperature. For instance, in the full 3D case, 50% of particles have a temperature of 378 K compared to 10% for axisymmetric 3D. This obviously is a result of the higher rate of evaporation in the full 3D case, which is again demonstrated in the distribution of the mass fraction of cyclohexane in Figure 12(c).
Again, although not evident in the averages of the particle parameters, the histograms seem to indicate that the full 3D model results in a higher percentage of particles at high temperatures, which in turn causes larger amounts of evaporation leading to smaller sized particles. For instance, a larger amount of particles with diameter less than 1000 µm are generated by the 3D calculations compared to the axisymmetric case. This can have significant implications in the design of the operating  Figure 12. (a) Cement particle size distribution (by diameter in µm) of escaped particles in outlet region; (b) particle temperature distribution of escaped particles; (c) mass fraction of cyclohexane per particle between axisymmetric 3D and fully 3D cases.
conditions for the contactor, and hence needs to be taken into consideration. For instance, according to the 3D CFD calculations, it would require less steam compared to the 3D axisymmetric calculations, in order to strip the cyclohexane from the polymer, thereby making the steam stripping process more efficient.
In summary, the current analysis in this paper demonstrates that although the flow and thermal distribution of velocity and temperature are qualitatively similar between the full 3D and 3D axisymmetric calculations, there are indeed some differences in the particle statistics such as particle diameter and cyclohexane particle mass fraction distribution, which can be critical in determining the design conditions for such contactors. Hence, even though the full 3D computations are performed at three times the cost of those with the axisymmetric assumption, the former might be necessary in carrying out accurate design optimization studies. It should also be noted that despite the deficiencies of the symmetry assumption, the computations presented in this paper are able to provided invaluable insights into this complex process for the first time, relatively efficiently.

Conclusion
The current study presents a CFD model based on a hybrid Eulerian/Lagrangian approach for the process of polymer devolatilization in a steam contactor. ANSYS R Fluent R simulations are carried out to model the flow, heat transfer, and phase changes in the contactor, in order to investigate the final polymer particle characteristics at the outlet. The polymer mixture originally has excess hydrocarbon solvent that needs removal. To remove this excess solvent, the mixture of polymer and solvent, or cement, is injected alongside with superheated steam, resulting in an instantaneous mixing process, and causing part of the solvent to evaporate rapidly. The solvent vapour and steam are vented away, leaving the remaining cement mixture with a lower concentration of solvent. As the cement droplets lose solvent and dry out, they are ejected as dry 'crumb' from the mixing device. A discrete-phase model (DPM) is used to inject and track cement droplets as they mix with superheated steam. The DPM allows us to monitor final cement particle sizes, the amount of solvent evaporation, and the change in particle temperatures that occur in the steam contactor. Additionally, the current study investigates the effect of initial cement temperature on the physical properties of the final particles, such as solvent concentration and particle size. These results show that, as the initial cement temperature increases, the particle size decreases slightly, while the cyclohexane mass fraction significantly decreases. More specifically, an increase in the initial cement temperature by 20 K corresponds to a 4.5% temperature increase. With this temperature increase from 423 to 443 K, the minimum diameter size decreased by 2.8% and the maximum diameter size decreased by 1.4%. In addition, the solvent content decreased by 13.9%. While the diameter does decrease in size, this change is not significant, which is due to the cement particle's already small size when entering the contactor. However, the resulting cyclohexane significantly decreased. Such a change in the amount of cyclohexane content present in the exiting polymer crumb is a critical factor in the refining process to produce quality polymers.
In addition, full 3D calculations were also presented in order to assess the validity of the axisymmetric assumption employed here in most of the analysis. Full 3D calculations showed differences in the evaporation rates thereby indicating drier and smaller polymer particles overall. For instance, 43% of the particles in full 3D calculations had diameters below a 1000 µm compared to 36% in axisymmetric 3D. The results demonstrated that even though detailed 3D computations are more computerintensive than 3D axisymmetric by a factor of three, they might be necessary for determining the design operating conditions for the contactor. However, the authors note that the analyses of the polymer devolatilization process presented here is the first of its kind. So despite the assumption and associated deficiencies, they can provide invaluable insights and future directions, especially for parametric studies, especially because they are computationally more efficient.
On the one hand, increasing the cement inlet temperature does not alter the particle size significantly, but results in a drier crumb exiting the contactor, which is favorable for the post-processing of the polymer in the coagulator. However, an increase in the cement temperature can increase production costs. Hence, it is important to converge to an optimal temperature that would increase neither the costs nor the particle size significantly and also keep the cyclohexane content to a reasonable level. Furthermore, the main goal was to understand the changes the particles undergo during the devolatilization process, and hence the model did not consider breakup, collision, or coalescence of the particles in the contactor. A lot of information about the correlations between operating conditions such as cement injection temperature and particle size or the outlet cyclohexane content can still be extracted from the current investigations. It is also worth noting that, among all the physical processes, the devolatilization process is the largest contributing factor to the change in particle size from the inlet to the outlet of the contactor. Such studies allow production companies quickly and cost-effectively to understand how significant is the impact of changing the flow conditions on the final product. The results for the escaped particles have been compared to the final product that is being physically manufactured, and is well within the expected range for particle size and cyclohexane content. With this information, the model can be used as an effective tool to study not only the effects of changes in operating conditions, but also contactor geometry alterations. The model presented here will pave the way for identifying control parameters in polymer isolation processes, thereby improving the design of equipment in such polymer manufacturing techniques.

Disclosure statement
No potential conflict of interest was reported by the authors.