Viscosity of molten polyolefins and stirring performanceThe mixing of highly viscous substances (viscosity (μ) ≥ 10 Pa s)23 is mainly characterized by the difficulty in reaching turbulent flows due to high viscous energy dissipation leading to excessive power consumption and local hot spots in the reaction medium. Polymer melts are known to be non-Newtonian fluids with a viscosity dictated by shear rate (a measure of the rate at which parallel internal surfaces slide past one another; Extended Data Table 1) and temperature, potentially leading to local variations in the reactor that could affect the mixing regime during operation32,33. We selected HDPE (Mw = 200 kDa, denoted HDPE200) and PP (Mw = 340 kDa, denoted PP340) grades found in consumer goods such as plastic caps, cars or textiles, with the characteristic Mw used to label the polymers obtained from their melt flow index34.We first conducted a rheological analysis, the results of which are shown in Fig. 1b (see Supplementary Tables 3 and 4 and Methods for details). As expected, increasing the shear rate led to a decrease in the viscosity at all of the tested temperatures, converging toward a distinctive value for each polymer (~500 Pa s for HDPE200 and ∼320 Pa s for PP340). Laminar flows with inherently low mixing capabilities are thus expected in polyolefin chemical recycling (Supplementary Note 1)35. Estimation of the stirring rates (N) required in a typical reactor vessel to reach the equivalent shear rates (secondary horizontal axis in Fig. 1b) revealed that temperature and local variations in viscosity in the melt under typical stirring rates can be disregarded (see Supplementary Note 1 for details).The torque (τ) required to stir a non-Newtonian fluid is proportional to the product of its viscosity and a power law of the stirring rate, that is, τ ∝ μNα (ref. 36), where α is an experimentally determined parameter dependent on the fluid. Taking into account that molten HDPE200 and PP340 have viscosities approximately one million times greater than that of water at room temperature (\(\mu_{{\rm{H}}_2{\rm{O}},298{\rm{K}}}\) ≈ 0.001 Pa s), we verified that the magnetic stirrers commonly available in laboratories are incapable of stirring high-molecular-weight polyolefins and presumably do not allow good control over the stirring rate for low-molecular-weight ones, as illustrated in Fig. 1c and Supplementary Videos 1 and 2. Magnetic stirrers are suitable for viscosities lower than ~1.5 Pa s (ref. 37), whereas mechanical stirrers can be functional up to 105 Pa s (ref. 38).Catalyst evaluation and computational fluid dynamics simulationsExperiments were developed in a four-parallel reactor set-up (see Methods, Supplementary Fig. 1 and also Supplementary Video 3 for the stirring configuration). Ruthenium nanoparticles supported on titania, a state-of-the-art catalyst for the conversion of HDPE200, was used throughout the study (see Supplementary Fig. 2 for the characterization of the Ru/TiO2 catalyst)39. The factors driving the performance (listed in Extended Data Fig. 1) within the experimental limitations described in Supplementary Note 2 were systematically varied and translated into computational fluid dynamics (CFD) simulations, including experimentally obtained viscosity dependencies, to describe the hydrogen–catalyst–melt contact over time (see Methods for a general description, Supplementary Note 3 for scope and Supplementary Fig. 3 for convergence plots).Initial experiments revealed optimum hydrogen pressures for HDPE200 and PP340 in accordance with the literature, with a pressure of 20 bar selected for subsequent tests (Extended Data Fig. 2, Supplementary Fig. 4 and Supplementary Tables 5 and 6)40,41, likely due to competitive adsorption of the polyolefin and hydrogen on the metal surface (Supplementary Fig. 5). Simulations at different temperatures confirmed minimum differences in the distribution of viscosities and Reynolds numbers (Extended Data Fig. 1 and Supplementary Table 7), in line with the observations depicted in Fig. 1b showing viscosity values almost independent of temperature at shear rates equivalent to stirring rates larger than approximately 15 r.p.m. in laboratory reactors (Supplementary Note 1). A commonly reported reaction temperature (498 K) was thus chosen for further analyses after performance tests (Extended Data Fig. 1)42,43.Internal and external mass transport limitationsWe first evaluated the ability of polymer chains to penetrate micropores and mesopores (see Supplementary Note 4 for the case of micrometer-sized pores usually found in shaped catalysts). The Freely Jointed Chain model (Supplementary Note 5) predicts a typical dimension for the folded chain (Λ) of ~ 22 nm for HDPE20044,45,46. For chain lengths below Λ, polymer chains tend to gradually favor the linear conformation47. The relevant pore size that polyolefins, or their liquid products following reaction, may not be able to penetrate thus ranges from Λ ≈ 1 nm (C6) to Λ ≈ 100 nm (high-molecular-weight polyethylene, Mw ≈ 5,000 kDa). These scales are represented in Fig. 2a and suggest that internal mass transport limitations in the case of porous catalysts may increase in relevance even for polyolefins with very high Mw as the reaction progresses toward shorter chain products, and therefore further studies are required. Internal mass transport phenomena were disregarded for simplicity as the Ru/TiO2 used in our study showed an average pore size of 6 nm within a very low specific pore volume of 0.02 cm3 g−1. Experiments at different stirring rates using a particle diameter of 0.6 mm support this assumption (Supplementary Fig. 6).Fig. 2: Characteristic lengths in catalytic polyolefin hydrogenolysis.a, Typical magnitudes of relevant length scales in the catalytic processing of consumer-grade plastics. Chain size refers to the typical dimensions of folded chains in low- and high-molecular-weight polymers. Reaction front refers to the characteristic penetration length of hydrogen in the melt before its concentration drops below 10% of its value at the H2-melt interface. b, Simulated decay of relative hydrogen concentration with distance from the H2–melt interface for different reaction rate constants (kr). \({\rm{C}}_{{{\rm{H}}_2},{\rm{int}}}\) refers to the concentration of H2 at the H2-melt interface. Pseudo-first-order kinetics (\(r=k_{\rm{r}}c_{{\rm{H}}_2}\)) were considered, as explained in more detail in Supplementary Note 5. The area shaded in gray indicates the typical range of kr values calculated in our experiments and reported in the literature. c, Schematic representation of the circulation of catalyst particles in the reaction vessel, with those exposed to hydrogen and polymer melt shown in green as active particles toward hydrogenolysis, as deduced from b, highlighting the fact that the reaction is mostly constrained to the vicinity of the H2–melt interface.Source dataRegarding external mass transport limitations, an earlier study determined negligible external hydrogen gradients to catalyst particles immersed in the melt if equilibrium bulk concentrations of H2 are reached25. However, the simulation of H2 diffusion into molten HDPE200 in the absence of reaction for a range of hydrogen pressures, times and viscosities (Supplementary Fig. 7 and Supplementary Notes 6 and 7) revealed a characteristic time for equilibration beyond typically reported reaction times. In view of this, we computed the decay of the H2 concentration at the H2–melt interface assuming the direct reaction of H2 with the melt after estimating that the observed reaction rate is around five times that of the diffusion rate of H2, as provided by the Hatta number (Extended Data Table 2 and Supplementary Note 7)48. Figure 2b shows the results for different reaction rate constants (kr), defined according to the expression for pseudo-first-order kinetics, \(r=k_{\rm{r}}c_{{\rm{H}}_2}\), where r is the rate of the reaction and \(c_{{\rm{H}}_2}\) is the concentration of H2. A suitable range of kr was estimated from the typical hydrogen consumption and reaction times observed in our study and reported in the literature (see Supplementary Note 7 for details)30. Poorly active catalysts not yielding any liquid products are characterized by kr values of ∼1.5 × 10−3 s−1, whereas highly active systems able to provide 100% conversion into methane are expected to present kr values of ∼0.1 s−1. These values translate into a range of concentration decays, highlighted in gray in Fig. 2b. Typically, we obtained values for kr of ∼0.01 s−1. As observed in Fig. 2b, the concentration of hydrogen drops below 10% of the interface value within a few millimeters in all cases, strongly suggesting that the reaction is mostly confined to the vicinity of the H2–melt interface with a typical length (λ) of ∼10−3 m, leaving most of the melt non-reactive (Fig. 2a). A representation in which only catalyst particles exposed to this region are active toward the reaction (Fig. 2c) is thus a suitable approximation to study the role of agitation in catalytic performance.Impact of catalyst particle circulation on performanceThe previous analysis shows the benefit of stirring configurations maximizing the presence of catalyst particles in the vicinity of the H2–melt interface. A first analysis based on the ratio of gravitational and viscous forces given by the Archimedes number (Ar; Extended Data Table 2) predicted Ar = 10−8–10−7 and therefore that the density of the catalyst is expected to play a negligible role in particle motion (Supplementary Note 8). Nevertheless, the average catalyst particle diameter (dp) is important as it determines the tendency of particles to follow melt streamlines according to the Stokes number (Stk; Extended Data Table 2)49. Considering λ as the characteristic length, Stk ≈ 10−2–100 for dp = 10−4–10−3 m. These values indicate that small catalyst particles will closely follow streamlines, whereas larger ones may deviate from them to an extent comparable to λ.We evaluated the importance of dp by comparing the catalytic performance of three different catalyst sieve fractions (0–0.2, 0.2–0.4 and 0.4–0.6 mm) in the hydrogenolysis of HDPE200 (Fig. 3a and Supplementary Tables 5 and 6). Equivalent experiments with PP340 did not show substantial variation in the total yield due to its low reactivity effectively limiting the performance of the catalyst (Supplementary Fig. 8). However, the same trend was confirmed by using the shorter and thus more reactive PP12 (Supplementary Table 8 and Supplementary Fig. 9). We found the smallest sieve fraction to be beneficial, producing a 40% greater yield of the C1–C45 products compared with the largest sieve fraction under the tested conditions. CFD simulations for dp = 0.2, 0.4 and 0.6 mm, keeping the same stirrer geometry, predicted differences in the particle trajectories. Larger particles on average required longer times to leave the bottom of the reactor and tended toward a more irregular occupancy of the vessel volume (Fig. 3b and Supplementary Fig. 10 for the three modeled particle sizes). Having determined the benefits of smaller sieve fractions, dp = 0.2 mm was used for the rest of the simulations in this study.Fig. 3: Influence of catalyst particle motion on performance.a, Variation in the product distribution for the hydrogenolysis of HDPE200 with catalyst sieve fraction using a propeller stirrer after 4 h. b, Corresponding three-phase CFD simulations using discrete phase modeling showing the trajectories of 0.2 and 0.6 mm catalyst particles under steady-state conditions: the total number of particles (np, with np = 196 for 0.2 mm and np = 58 for 0.6 mm; top) and representative initial trajectories of individual particles (bottom). c, Variation in the product distribution for the hydrogenolysis of HDPE200 with stirrer type after 4 h. d, Corresponding CFD simulations (top views) of catalyst particle trajectories for different stirrers, colored according to the H2 fraction in the vicinity. Simulations for other sieve fractions and parallel analyses for PP340 can be found in Supplementary Figs. 8–10 and Supplementary Tables 5–7. Simulated particle trajectories are presented in Supplementary Video 4. Reaction and simulation conditions: T = 498 K, \({p}_{{{\rm{H}}}_{2}}\) = 20 bar, catalyst/plastic ratio = 0.05 and stirring rate = 750 r.p.m.Source dataThe stirrer imposes the flow pattern that catalyst particles follow (Supplementary Fig. 11). Figure 3c shows the product distribution for the catalytic hydrogenolysis of HDPE200 using three different stirrer geometries under the same conditions (Supplementary Tables 5 and 6 and Supplementary Fig. 12 for PP340). The yield of C1–C45 products was not greatly affected by the stirrer type, whereas the product distribution shifted from gas to liquid fractions, with the amount of gaseous product decreasing in the order impeller > propeller > turbine, highlighting that stirring strategies can tune selectivity, as a consequence of tuning the activity, given that the hydrogenolysis is a series of reactions, and must be reported to facilitate benchmarking. The total number of carbon–carbon bonds followed the same trend, as determined using the recently published procedure for calculating the number of backbone scission, isomerization and demethylation events15 (Supplementary Table 9). Figure 3d and Supplementary Video 4 show critical differences between the stirrers. Propellers tend to split the catalyst particles into two separate zones with either high or low H2 concentration. Impellers tend to keep catalyst particles circulating around the mid plane, where the H2 concentration is high due to the V shape adopted by the H2–melt interface. Impellers are thus better suited to optimizing catalyst use. The turbine is poorly efficient in transferring particles to H2-rich zones, leading to the modest generation of gaseous products (which require more molecules of H2 per molecule of polymer). These effects can be quantitatively understood by considering the maximum value of the vertical component of the particle Reynolds number (Rep,z,max; Extended Data Figs. 3 and 4, Supplementary Notes 9 and 10, Supplementary Tables 7,10 and 11, and Supplementary Fig. 13) as the first performance descriptor. Rep,z,max can be derived from the melt properties, stirring rate, and particle and stirrer geometries and can be linked to activity and therefore changes in selectivity (Extended Data Fig. 2), offering a first tool to predict performance trends.Criterion to maximize the catalyst effectiveness factorThe yield of C1–C45 products did not monotonically increase with stirring rate, as shown in the catalytic tests for both polymers (Fig. 4a). The existence of an optimum rate, in accord with some reports in the literature on the conversion of low-molecular-weight plastics50,51, led us to study the influence of stirring rate on the extent of the H2–melt interface. CFD simulations developed for the three stirrer types (Fig. 4b) predicted the potential of propellers and impellers to increase the interface. The ability of simulations to reproduce the V shape of the H2–melt interface for highly viscous plastics is confirmed in Supplementary Video 3. CFD simulations of the impeller at different stirring rates strongly hint at a relationship between stirring rate and the H2–melt interfacial area (Fig. 4c). Small variations in the distance between the base of the stirrer and the bottom of the vessel also led to small changes in the H2–melt interface. However, an excessive distance (the top of the stirrer at the free melt surface) led to a decrease in the interface (Supplementary Fig. 14). In general, the shear-thinning character of molten plastics makes stirring only effective in the imaginary volume occupied by the stirrer under rotation or slightly beyond, making it advisable to minimize the distance between the stirrer and reactor walls.Fig. 4: Criterion for maximizing the effectiveness factor.a, Variation in the product distribution with stirring rates for HDPE200 and PP340 with the impeller stirrer. b,c, Two-phase CFD simulations of the hydrogen fraction in the mid z–x plane for different stirrer types (b) and different stirring rates for the impeller stirrer (c). d, Correlation between the effectiveness factor, defined as the ratio between the yield of C1–C45 and maximum yield of C1–C45 in a, and the modified power number for HDPE200 and PP340, calculated using the stirring rates in a and the simulated fraction of H2 in Extended Data Fig. 5. Reaction and simulation conditions: T = 498 K, \({p}_{{{\rm{H}}}_{2}}\) = 20 bar and catalyst/plastic ratio = 0.05.Source dataGiven the difficulty of calculating the extent of the H2–melt interface, we defined as a proxy the fraction of hydrogen (\(\chi_{{\rm{H}}_2}\); Extended Data Table 2) in a volume contained between the bottom of the stirrer and the H2–melt interface (Extended Data Fig. 1) when there is no stirring. The average Reynolds number in this region could serve as a descriptor for \(\chi_{{\rm{H}}_2}\) as more turbulence (larger Re values) may lead to more pronounced hills and valleys on the surface of the melt. However, Re is not observable. For the case of Re ≪ 1 as studied here, the power and Reynolds number are inversely linearly correlated. The power number (Np) expresses the relationship between resistance and inertia forces and can be written in terms of \(\chi_{{\rm{H}}_2}\) and observable variables such as the average density of the melt (\({\bar{\rho }}\)), the reactor diameter (D) and the average density of the melt (ρm) (equation (1), Extended Data Table 2 and Supplementary Note 11)24,52.$$\begin{array}{l}{N}_{{\rm{p}}}=\\\displaystyle\frac{2\uppi N\tau }{60\bar{\rho }(N/60)^{3}{D}^{5}}=\displaystyle\frac{7200\uppi \tau }{[\,{\chi }_{{\rm{H}}_2}{\rho }_{{\rm{H}}_2}+(1-{\chi }_{{\rm{H}}_2}){\rho }_{{\rm{m}}}]{N}^{\,2}{D}^{5}}\approx \displaystyle\frac{7200\uppi \tau }{(1-{\chi }_{{\rm{H}}_2}){\rho }_{{\rm{m}}}{N}^{\,2}{D}^{5}}\end{array}$$
(1)
Extended Data Fig. 5 shows the relationship between \(\chi_{{\rm{H}}_2}\) and Np obtained from CFD simulations for HDPE200 and PP340 under the same conditions as used in Fig. 4a. Propeller and impeller stirrers yielded volcano behavior, with a maximum \(\chi_{{\rm{H}}_2}\) ≈ 0.20–0.30 for Np ≈ 104–105, shifted toward slightly lower values in the case for PP340. The difference in the optimal rates in Fig. 4a and Extended Data Fig. 5 can be ascribed to the lower average viscosity in reaction compared with the simulations (Supplementary Note 12), although in practice this has a small impact as \(\chi_{{\rm{H}}_2}\) displays values of around 0.2 for a broad range of stirring rates.Plots of the effectiveness factor (η), defined as the ratio between the yield of C1–C45 products and the maximum yield of C1–C45 products over a series of experiments (equation (2) and Extended Data Table 1), versus the corresponding Np values based on the results presented in Fig. 4a show the optimal Np ranges for the two polymers (\(N_{{\rm{p}},{\rm{HDPE}}_{200}}\) ≈ 2 × 104 to 3 × 104 and \(N_{{\rm{p}},{\rm{PP}}_{340}}\) ≈ 1.5 × 104 to 2.5 × 104) to achieve high η values (Fig. 4d) and serve as a guide for the design of catalytic tests for performance optimization. From equation (1) and Extended Data Fig. 5, it is possible, for a given stirrer geometry (stirrer type and D), to select the stirring rate (N) and torque (τ) to be applied to deliver the desired Np value. Nevertheless, torque control is not a widely available feature of current reactor systems for catalyst evaluation, hindering the applicability of this criterion.$$\eta =\frac{{{\rm{Yield}}\;{\rm{C}}}_{1}-{{\rm{C}}}_{45}}{{({\rm{Yield}}\;{\rm{C}}_{1}-{{\rm{C}}}_{45})}_{\max }}$$
(2)
Use of the concentric cylinders model to describe the stirrer geometry (Supplementary Note 1) gives access to analytical relationships between viscosity, shear rate and torque, leading to an alternative expression for Np (Fig. 5 and Supplementary Note 11) that now includes contributions from the melt properties, stirring rate, fluid dynamics (through \(\chi_{{\rm{H}}_2}\)), and reactor and stirrer geometry (through D, Dr and L; Extended Data Fig. 1). All of the variables are either directly observable or design parameters, except \(\chi_{{\rm{H}}_2}\), which is available from Extended Data Fig. 1 and Supplementary Table 12 and depends on the stirrer type and plastic under treatment. Practitioners of catalysis can thus select appropriate combinations of stirring rate and reactor geometry to achieve the optimal Np ranges for a certain plastic. We note that deviations from the optimal range led to differences of up to 85% in activity and 40% in selectivity (Fig. 4a).Fig. 5: Application of the developed criteria for maximizing the effectiveness factor.Parameters that can be approximated under typical reaction conditions or are known a priori are indicated. Ranges of optimal stirring rates for a given reactor and stirrer geometry can thus be calculated. \({\bar{\mu }}\) refers to the average viscosity of the melt, Dr refers to the diameter of the stirrer and L to the height of the stirrer blades (Extended Data Table 1).In the most common case where the geometries of the stirrers and reactor are given, a first approximation to the optimal ranges of Np can be obtained from the values provided in Fig. 5. Melt densities and average viscosities at typical operation temperatures (Fig. 1b) and a reasonable value for \(\chi_{{\rm{H}}_2}\) of ∼0.2 (Extended Data Fig. 5) allow a straightforward calculation of stirring rate ranges. For example, in the case of D = 2 cm, L = 1 cm and Dr = 2.5 cm, the approximate ranges for high catalyst effectiveness factors would be N = 880–1,300 r.p.m. for HDPE200 and N = 760–1,100 r.p.m. for PP340 This criterion was shown to be valid in the range of most reported operation pressures (20–30 bar), with lower pressures (10 bar) showing behavior compatible with H2 depletion as the limiting factor (Supplementary Table 5 and Supplementary Fig. 15). These results, together with the small variation in viscosity at commonly applied temperatures (Fig. 1), make this criterion pressure- and temperature-independent under most reported conditions.Model scope and future directionsAs the average chain length of the hydrocarbons decreases due to cleavage, so does the viscosity, spanning six orders of magnitude until reaching values close to water (Fig. 1). Thus, the ability of the criterion to predict performance as the reaction progresses was next investigated.We hypothesized that the transition from the initial non-Newtonian character to a Newtonian character, facilitating the creation of turbulence24, may change the structure of the H2–melt interface. The Freely Jointed Chain model predicts a transitioning chain length of around C200 (Supplementary Note 5). With this in mind, we simulated stirring patterns for HDPE100 (non-Newtonian, a proxy for low conversion stages), a hypothetical C200 under Newtonian and non-Newtonian regimes, and eicosane (C21, Newtonian, a proxy for high conversion stages). The results clearly reflect the transition from a single H2–melt interface to an abundance of H2 bubbles populating the melt (Extended Data Fig. 6 and Supplementary Fig. 16), as supported by direct observations when turbulence starts to dominate as viscosity decreases (Supplementary Video 3). We then performed catalytic tests on HDPE100 and eicosane (Supplementary Table 8 and Fig. 6a,b), calculated \(\chi_{{\rm{H}}_2}\) (Supplementary Table 13 and Supplementary Fig. 17) and applied the criterion (Fig. 6c). The non-Newtonian melts of HDPE200 and HDPE100 exhibited very similar trends, with identical optimal stirring rates (although different Np due to different viscosities), whereas eicosane displayed a C-shaped relationship between effectiveness and Np, clearly suggesting the need for a different modeling strategy for the later stages of the reaction (or for the case of catalytic hydrogenolysis of surrogate molecules or the often-used very-low-molecular-weight plastics). The transition seems to occur at Np = 102–103, corresponding to viscosities of around 3–30 Pa s at 1,000 r.p.m., therefore validating the proposed criterion until the later stages of the reaction.Fig. 6: Model scope and influence of thermal gradients.a,b, Variation in the product distribution with stirring rate and two-phase CFD simulations of the hydrogen fraction in the mid z–x plane for HDPE100 (a) and eicosane (b). c, Correlation between the effectiveness factor, defined as the ratio between the yield of C1–C45 and the maximum yield of C1–C45 for HDPE and between the yield of methane and the maximum yield of methane for eicosane, and the modified power number, calculated using the simulated H2 fractions in Extended Data Fig. 5, using an impeller as stirrer. d, Temperature distribution in the mid x–y plane for different stirrer geometries when the thermocouple reaches the operation temperature (498 K, at the position indicated). e, Temporal evolution of the temperature distribution in the x–z plane for different reactor diameters (Dr). Reaction and simulated conditions: T = 498 K, \({p}_{{{\rm{H}}}_{2}}\) = 20 bar and catalyst/plastic ratio = 0.05.Source dataIn addition to the analysis of mass transport limitations, we also investigated heat transport constraints. We simulated the largest possible temperature gradient within the reactor during operation with three different stirrer configurations when the temperature at the thermocouple reaches the set temperature (498 K in our case, equal to that imposed on the reactor walls). Figure 6d shows the temperature distribution in the reactor, which resembles that of the Reynolds number distribution (Extended Data Fig. 2), with gradients of approximately 100 K for the best impeller and propeller geometries. This led us to conduct temporal simulations to predict the time for the gradient to reduce to less than 10 K. Figure 6e (top) shows a time of 15 min for the worst-case scenario of walls at the set temperature and the interior at room temperature at t = 0 with a stagnant and non-reactive melt (see Methods for more details), representing only 6% of the operation time (4 h). From a forward-looking perspective, more refined models able to stepwise predict product distributions and thus recommend optimal operation times will be possible after incorporating kinetic descriptions for the catalyst under study. Alternatively, developing operando tools to track viscosity could also guide optimal reaction times. We also highlight the generality of the applied analysis that could be adapted to future reactor architecture operating in continuous mode. In this direction, processes such as continuous reactive extrusion53 are first steps, which would also enable the online analysis of products.
Assessment of transport phenomena in catalyst effectiveness for chemical polyolefin recycling
