Simulation approach to study kinetic heterogeneity of gadolinium catalytic system in the 1,4-cis-polyisoprene production

1 Introduction

Industrial production of stereoregular polydienes is predominantly carried out using Ziegler–Natta catalysts (1,2). Titanium-based catalytic systems are commonly employed in polyisoprene production, resulting in products with cis-1,4 units of up to 96% (3). On the other hand, neodymium-based catalytic systems from the lanthanide series polymerize dienes with even higher cis unit content (up to 98%) (4). Consequently, the finished product using neodymium catalysts finds successful applications in the tire and rubber industry. However, researchers have been increasingly exploring gadolinium compounds (5,6) to enhance physical and chemical properties and discover alternative catalysts. Gadolinium-based catalysts exhibit the ability to yield products with even higher cis-1,4 unit content (up to 99.1%) and lack gel fractions when compared to neodymium catalysts.

Microheterogeneous Ziegler–Natta-type catalysts result in multiple kinetically nonequivalent active centers (ACs) that produce macromolecules of specific molecular weights (7). A number of factors that can determine the stereospecificity of ACs include defects that appear on the crystal lattice. In addition, the presence of several types of ACs can be explained by the lack of equipotentiality on any solid surface, and any heterogeneous catalytic system has a fairly characteristic distribution of catalyst particles, which is determined by the degree of dispersion. Another factor that determines the nature of kinetic heterogeneity is the nature of the transition metal and its valence state.

Consequently, further empirical investigation of the process, using a model system definition, reveals the necessity to study the kinetic heterogeneity of the catalyst by solving the inverse problem based on primary physicochemical data, represented as a molecular weight distribution (MWD) curve. In the review of Marien et al. (8), the issue of using different representations of molar mass distribution in polymerization processes is addressed, emphasizing the translation and correction of simulated chain length distributions to facilitate comparisons with experimental size exclusion chromatography (SEC) traces and providing guidelines for reporting in kinetic studies. Stochastic simulation algorithm methods have a huge impact on the results obtained (9). Analyzing the MWD through the model system definition enables the identification of kinetically distinct ACs within the system, each capable of forming polymer fractions with varying molecular weights. The dynamics of these AC data can then be utilized to set promising tasks for planning and optimizing polymer production processes. In world science, alternative approaches to assessing kinetic parameters are also successfully used. Salas et al. (10) evaluated and implemented a data-driven approach for online estimation of kinetic parameters in the copolymerization of ethylene with 1,9-decadiene using a unique catalyst, utilizing global sensitivity analysis and introducing the RCMR algorithm for parameter estimation.

This study seeks to create digital tools for solving the inverse problem of MWD and to describe the kinetic heterogeneity of a gadolinium-based catalytic system in the 1,4-cis-polyisoprene synthesis.

2 Materials and methods

The classical approach to solving the inverse problem of MWD formation involves introducing the Frenkel statistical parameter λ and the search for the distribution function of ACs φ ( λ ) based on the expression (11)

where K ( λ , M ) is the kernel of the integral equation reflecting the mechanism of the polymerization process, and q exp ( M ) is the vector determining the product’s MWD. In a scenario where multiple distinct ACs coexist, the initial MWD curve is represented as the superposition of overlapping distributions that are specific to each particular mechanism governing the polymerization process:

The Flory distribution (12) is the most commonly used as the kernel of the integral equation for microheterogenic catalytic systems:

To facilitate calculations, it is common to switch to logarithmic coordinates (x = ln M and s = ln λ) since the calibration in gel permeation chromatography (GPC) determines the molecular mass distribution in q _w(M)–ln M coordinates.

Based on the type φ(λ), appropriate conclusions can be drawn about the contribution of each AC to the total MWD. However, the problem (Eq. 1) is ill-posed. It can be successfully solved using the regularization method proposed by academician Tikhonov and Arsenin (13), which has been previously applied to restore the character of kinetic heterogeneity in catalytic systems based on titanium and neodymium (14).

Despite the success of the Tikhonov regularization method, there are challenges related to the high sensitivity of this approach to initial experimental data, as shown in software implementation and computational experiments (15). In some cases, with an error of laboratory measurements exceeding 10%, estimating the proportion of ACs that produce polymer fractions with the lowest molecular weight becomes nearly impossible (16). The correct choice of the Tikhonov regularization parameter and the occurrence of high oscillations in the resulting solution at the ends pose additional challenges and require further analysis (17).

To address these issues, the author proposes an alternative method for solving the inverse problem. The new simulation approach is based on the assumption that the distribution of ACs can be described by one of the model distributions. It automates the process of selecting the number of ACs, their fractions, and kinetic parameters to produce polymers with a given molecular weight on each AC.

The main idea of the simulation approach is to programmatically explore different polymerization process scenarios based on the number of ACs (N), the proportion of each AC in the catalyst’s total composition p i , and the kinetic parameters defining the Frenkel statistical parameter λ i . This exploration leads to different MWD variants, and the best-fitted MWD, reproducing the experimental results, identifies the desired set of parameters p i and λ i ( i = 1 , … , N ). Unlike the classical approach that aims to solve the inverse problem, the new method involves multiple solutions to the direct problem, superposing the MWD distribution characteristic of each type of AC, and selecting a set of parameters p i and λ i that closely approximate the experimental MWD.

The algorithm starts with only one type of active center and determines possible kinetic parameters to match the MWD to the experimental data. If the difference between the experimental and calculated MWD exceeds the permissible error, the number of ACs is increased to two, and kinetic parameters corresponding to each type are determined. This iterative process continues until the difference between the experimental and calculated MWD is smaller than the maximum allowable error size.

To enable multiple searches for different MWD variants, the step-by-step algorithm was optimized, and all calculated dependencies were taken out of the iterative search. The shift step, ∆(λ), was made to functionally depend on the error value (18), ensuring efficient exploration of the solution space.

Following the key concept of the simulation approach to solving the inverse problem, we present the numerical solution algorithm as a sequence of steps that facilitate the accurate identification of the number and fractions of ACs responsible for kinetic heterogeneity.

Step 1. Conducting a laboratory experiment to detect the MWD pattern in a discrete representation. The result yields an array of points q exp [ M i ] , i = 1 , … , N , where N is the number of split points of the original segment.

Step 2. Determine the permissible error value ε > 0 .

Step 3. Set the number of different types of ACs m = 1.

Step 4. Define a set of arrays for all possible combinations of different types of ACs p i = { p 1 , p 2 , … , p m } , where p 1 + p 2 + … + p m = 1 .

Step 5. Set the initial values for the quantities characterizing the Frenkel statistical parameter λ i = λ 0 and determine the final boundaries λ i for each of them in the form of a single quantity λ f in .

Step 6. Calculate the MWD curve as a distribution superposition characteristic of each AC type according to the parameters λ i and p i . The result is an array of points q calc [ M i ] , where i = 1 , … , N .

Step 7. Determine the value of the residual function:

Step 8. Compare the resulting value δ 2 with the permissible error ε . If δ 2 ≤ ε , the calculation is finished and it is necessary to proceed to step 11, otherwise the iterative search continues.

Step 9. Shift the parameter λ i by an amount ∆ ( λ ) : λ i = λ i + ∆ ( λ ) . If for all the elements of the set p i = { p 1 , p 2 , … , p m } the equation λ i > λ f in holds true, then proceed to step 10, otherwise go back to step 6.

Step 10. If all the elements of this set p i have already been considered in steps 5–9 with δ 2 ≥ ε , then increase the number of AC types m = m + 1 and proceed to step 4 to compile a new set p i = { p 1 , p 2 , … , p m } .

Step 11. We write out the resulting solution in the form of a set of parameters m , λ i , and p i , where i = 1 , … , m .

A significant optimization parameter of the algorithm involves identifying the shift step value, ∆ ( λ ) , based on the calculated error value δ 2 . The running speed of the algorithm is affected by the accuracy of determining the parameters of the set p i = { p 1 , p 2 , … , p m } that characterize the proportions of different types of ACs.

All computational experiments were conducted with higher calculation speed when, at the initial stage of the algorithm implementation, the accuracy, determined by the step value of p i , was 5%. After localizing the domain of a possible solution, the accuracy of determining p i was increased.

Let us evaluate the given approach to solving the inverse problem of the MWD formation and the algorithm operation for studying the kinetic heterogeneity of a catalytic system based on gadolinium chloride solvate in the 1,4-cis-polyisoprene synthesis.

3 Results

The following chemicals were used in the experiments: isopentane (main substance content 99.6 wt%), isoprene (supreme grade, main substance content not less than 97.0 wt%), piperylene (main substance content not less than 97 wt%), triisobutylaluminum (technical, supreme grade, toluene solution, main substance content not less than 40 wt%), diisobutylaluminum hydride (weight fraction of diisobutylaluminum hydride not less than 40.0 wt%, weight fraction of triisobutylaluminum not more than 20.0 wt%), nitrogen gas (higher purity, grade 2, nitrogen content not less than 99.95 vol%, oxygen content not more than 0.05 vol%, water vapor content not more than 0.004 vol%), toluene for size-exclusion chromatography (main substance content not less than 99.9 wt%, moisture content not more than 0.01 wt%), and S-789 antioxidant (N-2-ethylhexyl-N′-phenyl-p-phenylenediamine (8 PPDA), main substance content not less than 95.5 wt%).

According to the experimental conditions, the initial molar composition of the catalyst was formulated based on the ratio of triisobutyl aluminum/gadolinium (GdCl₃·ni-PrOH·kH₂O)/piperylene = 20/1/2.5, with the catalyst concentration for gadolinium being 0.024 mol·L⁻¹. Diisobutylaluminum hydride was utilized to regulate the molecular weight of the product.

The periodic polymerization process of isoprene was conducted in an autoclave-type reactor (volume V = 2 L) equipped with a mixing device and a cooling jacket. The isoprene concentration in isopentane was 16.7 wt%. Diisobutylaluminum hydride was introduced at the rate of 3.6 × 10⁻⁵ mol per 1 mol of isoprene. The gadolinium catalyst rate was 1 mol of gadolinium per 15.6 × 10³mol of isoprene. The polymerization process was continued until 75% isoprene conversion. The process was deactivated with ethyl alcohol and the resulting polyisoprene solution was stabilized with the antioxidant N-2-ethylhexyl-N′-phenyl-paraphenylenediamine (8 PFDA) and released by evaporating the solvent in the temperature range of 75–80°C (19). The primary analytical method for characterizing molar mass distribution is SEC, also known as GPC, where longer polymer chains, or species with greater molar masses, elute earlier due to their reduced interactions with the smaller pores within the chromatographic columns. The findings of GPC analysis (20) of the resulting polyisoprene are presented in Table 1.

The molar-mass dispersity of the product is M w / M n = 4.44, indicating that the obtained MWD cannot be described by one type of AC, and the third step of the presented algorithm began with the value m = 2. The permissible error of reproduction, ε, is assumed to be 0.03.

To comply with the requirements of Step 5, all possible combinations of parameters were compiled as p i = {(0.05; 0.95), (0.1; 0.9), (0.15; 0.85), (0.2; 0.8); (0.25; 0.75); (0.3; 0.7), (0.35; 0.65), (0.4; 0.6); (0.45; 0.55); (0.5; 0.5)}, determining the fractions of each of the ACs, and values of kinetic parameters were iteratively chosen to find the values λ 1 and λ 2 corresponding to each AC.

At Step 6 of the algorithm, for each of the presented sets { λ i , p i | i = 1 , . . . , m }, the corresponding vectors defining the discrete representation of the calculated MWD curve were obtained, and at the next step, the vector found q calc [ M i ] , i = 1 , . . . , N was compared with the experimental MWD representation to calculate the error δ 2 . Under rough estimation conditions, with an accuracy of parameter identification p i at 5%, the best solution was achieved for the parameters 0.4 ≤ p 1 ≤ 0.5 and 0.5 ≤ p 2 ≤ 0.6 . The refinement of the existence domain of a possible solution led to identifying the values p 1 = 0.44   and p 2 = 0.56 for λ 1 = 3.6 × 10⁻⁶ and λ 2 = 8 × 10⁻⁷. The results of the calculations are presented in Figure 1.

Figure 1

The results of solving the inverse problem of the MWD formation for two types of ACs (the line is the initial curve of MWD; the stroke is the calculated curve).

The analysis of Figure 1 reveals that the calculated MWD, resulting from the distribution superposition characteristic of each of the 2 AC types, does not agree with the experiment findings, and the error value δ 2 , according to the calculations, is 0.277. It is evident that two ACs in the system are not sufficient to describe the initial MWD. In this case, the value of m increases to 3, and the iterative search continues for a more complex set p i = {(0.05; 0.05; 0.9), (0.05; 0.1; 0.85), (0.05; 0.15; 0.8), (0.25; 0.25; 0.5)}.

The first approximation showed that the best solution was achieved for the parameters 0.1 ≤ p 1 ≤ 0.2 , 0.4 ≤ p 2 ≤ 0.5 , and 0.4 ≤ p 3 ≤ 0.5 . The refinement of the existence domain of the possible solution led to identifying the values p = (0.14; 0.44; 0.42) for λ = (1.34 × 10⁻⁵; 2.45 × 10⁻⁶; 7 × 10⁻⁷). The resulting calculations are presented in Figure 2, demonstrating that the calculated MWD agrees well with the experimental data. The error value δ 2 in this case is 0.006, indicating the end of the iterative search and there was no need to increase the parameter m.

Figure 2

The results of solving the inverse problem of the MWD formation for three types of ACs (the line is the initial curve of MWD; the stroke is the calculated curve).

Thus, a new simulation approach to solving the inverse problem of the MWD formation successfully identified the presence of three ACs: type A _Gd – ln M = 11.2, type B _Gd – ln M = 12.9, and type C _Gd – ln M = 14.2, with the proportion of ACs of type A _Gd is 0.14; type B _Gd – 0.44; and type C _Gd – 0.42.

Further mathematical description of the polymerization process under different external factors often requires refining a number of kinetic parameters. This refinement is accomplished through the inverse problem of chemical kinetics.

4 Discussion

This study presents a novel simulation approach to address the inverse problem of kinetic heterogeneity in polymerization processes. The approach is founded on a fundamental assumption that the distribution of ACs can be described by specific model distributions. By utilizing primary physicochemical data, such as the polymerization rate and the MWD of the resulting polymer product, the simulation method automatically identifies the kinetic parameters and effectively solves the problem of kinetic heterogeneity.

The proposed simulation approach was rigorously tested and validated in the context of producing polyisoprene using a catalytic system based on gadolinium chloride solvate, with a molar composition of TIBA/GdCl3·ni-PrOH·kH₂O/Pip = 20/1/2.5. Through the developed methodology, the study successfully determined the existence of at least three distinct types of ACs responsible for initiating the polymerization process. These ACs were categorized as follows:

1. Type A _Gd with a statistical parameter ln M = 11.2, contributing to approximately 14% of the polymerization process.

2. Type B _Gd with a statistical parameter ln M = 12.9, accounting for around 44% of the polymerization process.

3. Type C _Gd with a statistical parameter ln M = 14.2, constituting about 42% of the polymerization process.

The results obtained from the simulation approach showcased the capability to unravel the complex kinetic heterogeneity of the catalytic system and enabled the accurate identification of the fractions and kinetic parameters of each AC involved in the polymerization process. This valuable insight will undoubtedly facilitate further investigations into optimizing reaction conditions, product properties, and potentially lead to advancements in industrial polymer synthesis processes.

In conclusion, this pioneering simulation approach marks a significant step forward in solving the inverse problem of kinetic heterogeneity, shedding light on the complexities of catalytic systems and their impact on the polymerization process. The study’s findings open up new possibilities for understanding and controlling polymerization reactions, ultimately contributing to the development of innovative materials and industrial applications.

5 Conclusions

Approbation of the simulation approach for the production of polyisoprene in the presence of a catalytic system based on gadolinium chloride solvate with a molar composition of TIBA/GdCl3·ni-PrOH·kH₂O/Pip = 20/1/2.5 based on the developed methodology made it possible to identify the presence of at least three types of ACs initiating the polymerization process: type A _Gd – ln M = 11.2, fraction – 0.14; type B _Gd – ln M = 12.9, fraction – 0.44; and type C _Gd – ln M = 14.2, fraction – 0.42.