Zhaoxiang Zhang ,Fei Guo,Wei Song ,Xiaohong Jia ,Yuming Wang

1 State Key Laboratory of Tribology,Tsinghua University,Beijing 100084,China

2 Joint Research Center for Rubber and Plastic Seals,Tsinghua University,Beijing 100084,China

3 State Key Laboratory of Advanced Forming Technology &Equipment,China Academy of Machinery Science &Technology,Beijing 100044,China

4 Tianding Sealing Technology (Beijing) Co.,Ltd.,Beijing 100072,China

Keywords:Thermal reaction Polymer processing Reaction kinetics Mathematical modeling Empirical correction

ABSTRACT Based on the theory of first-order reaction kinetics,a thermal reaction kinetic model in integral form has been derive.To make the model more applicable,the effects of time and the conversion degree on the reaction rate parameters were considered.Two types of undetermined functions were used to compensate for the intrinsic variation of the reaction rate,and two types of correction methods are provided.The model was explained and verified using published experimental data of different polymer thermal reaction systems,and its effectiveness and wide adaptability were confirmed.For the given kinetic model,only one parameter needs to be determined.The proposed empirical model is expected to be used in the numerical simulation of polymer thermal reaction process.

1.Introduction

Formation of a crosslinked network during thermoforming imparts unique mechanical,chemical and thermodynamic properties to the polymer[1,2],making these polymers widely available in industrial practice[3–5].It is important to accurately describe the reaction process of polymer molding,such as the curing and vulcanization processes,to not only enable researchers to understand the reaction mechanism,but also to enable them to predict the thermal reaction process by numerical simulation,thereby helping the engineer determine the processing parameters for precise molding and energy saving[6,7].This goal usually needs to be achieved using a reaction kinetic model.Here,the reaction kinetics model is a model describing the relationship between the degree of polymerization or cross-linking conversion and the time and temperature.The model can be either theoretical or empirical[8–13].Through experimental design,the researchers measure the conversion degree under different temperatures(or heating rates),substitute the data into the selected model,and then fit the model parameters,thus obtaining a suitable reaction kinetics model.However,there are still some objective problems with the existing reaction kinetics model,which will be described in the following sections.

Existing reaction kinetics models are not simple enough and often containing too many fitting parameters,making them difficult to implement.Moreover,different fitting parameters are often needed at different temperatures[14–18],which makes prediction of the conversion degree and numerical simulation at different temperatures difficult.At present,in order to meet the needs of different applications,researchers use various blending and filling technologies to improve the thermal,mechanical,and chemical properties of polymers [14,19–22].That is to say,more and more reaction systems have been developed,which makes the reaction kinetic models used more and more diverse.To accurately describe the thermal reaction process of the polymer after treatment with these modifications,it is often necessary to use multiple kinetic models.Therefore,a simple,accurate,and powerful kinetic model is required by researchers and engineers.

Besides,the existing models are mostly given in the form of a function of the conversion rate (dα/dt).One of the disadvantages of such models is that when there are fluctuations in the measured data,the value of the dα/dt is difficult to accurately determine,and it can even be negative[23].On the other hand,the use of artificially smooth methods may introduce additional errors and require more test data points.At present,there is a lack of a kinetic model directly given in the form of a function of the degree of conversion.

To compensate for the above shortcomings,a thermal reaction kinetics model in the integral form is derived in this paper.To make the model more applicable,on the basis of the derived model,the effects of time and the conversion degree on the reaction rate parameters (pre-exponential factor and apparent activation energy) are considered,and two types of undetermined functions are used to compensate for the intrinsic variation of the reaction rate.The model is explained and verified using published experimental data for different thermal reaction systems.The proposed kinetic model only requires the activation energy of the thermal reaction to be determined,and not the preexponential factor,which greatly reduces the workload.The proposed model is expected to provide a more powerful method for various thermal reaction processes.

2.Methods

2.1.Theoretical method

For first-order reaction kinetics,the reaction rate can be described as

where c (mol.L-1) is the concentration of the precursor component that can form the cross-linking network.k(s-1)is the reaction rate,which is solved according to the Arrhenius equation:

where A0(s-1) is the pre-exponential factor,Ea(J.mol-1) is the activation energy constant of the reaction,T(K)is the absolute temperature at which the reaction is located,and R (J.mol-1.K-1) is the ideal gas state constant.These two equations now establish the relationship between temperature and the reaction progress.

It should be noted that the reaction of the polymer thermoforming process may be performed in multiple steps,but there may be one of them determines the overall reaction kinetics.Here,as mentioned in other literature[8,24],the thermal process of polymers is regarded as a total reaction.Assuming that the total concentration of precursors that generates the target substance (network) is c0(mol.L-1).Under the isothermal condition,according to Eq.(1),the remaining precursor concentration in the system after time t can be expressed as

For non-isothermal conditions,the time t (s) needs to be divided into multiple time microelements.In each infinitely small time micro-element,the reaction can be regarded as proceeding at a constant temperature,and then Eq.(4) will be obtained.

where k1,k2,knis the reaction rate at different times,substituting Eq.(3),Eq.(4c) can be expressed as

where T(t)is a function of temperature with respect to time,so the degree of conversion(conversion fraction)can be written as Eq.(6)[25].

Eq.(6) is an integral expression of the degree of thermal reaction conversion.It can both solve the degree of conversion under isothermal conditions and can also be applied to non-isothermal conditions.In addition,a major advantage of this integral form model is that it is not necessary to solve the derivative term of k when solving the conversion rate (dα/dt),so it is not necessary to solve the term of dα/dT.Even under cooling conditions,the change in the degree of conversion calculated using this equation is still positive.Therefore,this type of integral model can be used to solve the changes in the conversion degree during both the heating and cooling processes.This is critical for predicting the postcooling molding phase of polymer materials,which can help reduce energy consumption in the manufacturing process and also make the molding quality of parts more uniform.

2.2.Data acquisition

In order to verify and modify the model,experimental data of thermal conversion curves of polymers at different temperatures are needed.In general,selected polymers are composed of crosslinked networks,and the thermal reaction process is accompanied by the formation or decomposition of these polymers.Specifically,we mainly selected two types of polymers,rubber and resin.The vulcanization process of rubber is widely accepted to be in accordance with the first-order kinetic reaction process,which is consistent with the basic basis of the model derivation process.Resin curing reaction is considered to be not in accordance with the first-order reaction kinetics of thermal reaction process,but the model can be well fitted by the empirical correction method given in the following part.The relevant data are extracted from a series of publishing data.For generality,these conversion degree curves were subordinated to the thermal reaction process of polymers with different types and different reaction systems.In the following,the data were calculated to obtain the required functional relationship curve,such as ln(1-α)～t or -ln[-ln(1-α)]～1/T,for description and verification of the given model.It should be emphasized that the original source of each set of processed data is specified.

3.Results and Discussion

3.1.Verification of the theoretical model

After theoretically deriving Eq.(6),its correctness needs to be verified,because the equation was obtained on the basis of multiple assumptions.There are many studies on polymer thermal processing,which is sufficient to verify the model and determine its shortcomings.If the model is correct,the curve with ln(1-α)～t should be a straight line under the isothermal condition.Here,the conversion degree curve obtained in the literature was transformed to a ln(1-α)～t curve,as shown in Fig.1(a) and (b).It can be seen that some polymer thermoforming processes (mostly rubber vulcanization processes) strictly conform to this law from the beginning to the end of the reaction.However,there are many other thermal reaction processes that do not follow this law(Fig.1(c–f)).The ln(1-α)～t curve of thermal reaction process of these polymers is only a straight line at a specific temperature or in a specific time period,and there are great differences under different temperatures.It is obvious that in the case mentioned later,the ln(1-α)～t function line of the thermal reaction is curved.At this time,Eq.(6) cannot accurately describe the thermal reaction change process.This is because many polymer crosslinking processes are not considered to follow the first-order reaction kinetics,such as the curing process of resin,so the above equation needs to be modified empirically to expand its applicability.It may be because of the effect of diffusion,which can be eliminated by correction of the reaction rate [31].Since the polymer thermal reaction process is closely related to the reaction time.It is believed that at different stages of the reaction,due to the formation or decomposition of cross-linked bonds,the restricted molecular mobility of cross-linked bonds will change accordingly,resulting in the continuous change of the mobility of the reactive groups,which will lead to the constant change of the reaction rate with the progress of the thermal reaction.Alternatively,it may be caused by changes in the kinetic parameters during the reaction,that is,Eaand A0are changed[32,33],which breaks the assumption that Eaand A0are constants in the Arrhenius equation,resulting in the inaccuracy of the theoretical model describing the thermal reaction process.Two empirical correction methods,timedependent and conversion degree dependent,are proposed in here.The corresponding function is introduced to describe the change of Eaand A0with the continuous progress of the reaction (reaction time and conversion degree change constantly),so as to compensate for the variation of these parameters with the reaction process.

Fig.1.ln(1-α)～t curves for the curing/vulcanization of(a)a NR/SBR stoichiometric mixture(the values were calculated from the experimental data reported by[26]),(b)a silicon rubber (the values were calculated from the experimental data reported by [27]),(c) a self-developed resin for autoclave (the values were calculated from the experimental data reported by [28]),(d) a cyanate ester resin (the values were calculated from the experimental data reported by [18]),(e) a rubber (the values were calculated from the experimental data reported by [29]),and (f) a Ciba-Geigy epoxy resin (the values were calculated from the experimental data reported by [30]).

3.2.Empirical correction of the model by considering the time effect

We first used the time-dependent function to comprehensively compensate for the changes of the reaction rate parameters in the model.It is assumed that the following equation is valid on the basis of Eq.(6):

where f(t)is a function of time,assuming that the expression of f(t)is the same at different temperatures for the same material.The following equation is then established to derive the time from both sides of the above equation:

Rearranging Eq.(8) gives

Note that F(t)is the primitive function of f(t),and integrate the above formula to get Eq.(10).

where 0 ≤α<1.From the initial condition α|t=0=0,C=0 is obtained.

Now substitute Eq.(2) into Eq.(10) to get the following equation:

where g(t)is a function of time,and the function g(t)is introduced to eliminate the pre exponential factor.In this way,the function g(t)only needs to be fitted when the parameters are solved,and the value of A0is no longer needed to be determined.Taking the natural logarithm of Eq.(11) gives

Because g(t) is a function of time,the curve fitting with-ln(1-α).eEa/RTor ln[-ln(1-α)]+Ea/RT as ordinate and time as the abscissa at different temperatures should satisfy the same functional relationship.The values of conversion degree and time can be obtained by experiments.Thus,only the value of the activation energy Eain Eqs.(11)and(12)is unknown.It can be found that Eq.(12) can be rearranged to

For a given time t,lnA0F(t) is a definite value.Therefore,the slope of the curve -ln[-ln(1-α)]～1/T is the value of Ea/R.The relevant curve was constructed using data in the literature(Fig.2),and the corresponding value of Ea/R was obtained.

After the value of Ea/R was obtained,the g(t) curve was drawn(Fig.3).It can be seen that the g(t) function of polymers at different temperatures indeed obey the same function law.In some cases,as shown in Fig.3(b–d),although the fitted lines of g(t)～t cannot overlap into the same line,but g(t)～ln(t) or ln[g(t)]～ln(t) can be plotted as the same function.It should be noted that regardless of which form of fitting is adopted,our aim is to obtain the expression of g(t).After fitting the functional relationship of g(t),the corresponding conversion degree model can be obtained by simultaneous solving Eqs.(7),(8) &(11).It should be noted that under non-isothermal conditions,Eq.(14)should be solved numerically and written as an incremental form like Eq.(4c).

In other cases,it was found that although the fitted curves at different temperatures obeyed the same function law,their intercepts were different and they did not overlap (Fig.4).This may be caused by the difference in the primitive function in Eq.(10)at different temperatures.In this case,the relationship between the temperature and the intercept of the function needs to be determined in order to establish an accurate reaction kinetic model.In addition,another empirical correction method described below can also be used for model fitting.

3.3.Empirical correction of the model by considering the conversion degree effect

Similar to the method described in Section 3.1,it is assumed that the following equation holds on the basis of Eq.(6):

where f(α) is a function of the conversion degree.This assumption is made because many studies have pointed out that the activation energy and pre-exponential factor of the reaction are not fixed values during the process of the thermal reaction,but they are a function of the degree of conversion [32,33].Therefore,a function related to the conversion degree was used to comprehensively compensate for the changes in the reaction rate parameters in the model.Assuming that the expression of f(α)is the same at different temperatures for the same material.Then,the following equation is established to derive the time from both sides of the above equation:

Rearranging Eq.(16) to gives

The integral of the above formula gives

where 0 ≤α<1.From the initial condition α|t=0=0,C=0 is obtained.

Note that F(α) is the primitive function of 1/(1-α)f(α):

Now substitute Eq.(2) into Eq.(18) to get the following equation:

Fig.2.-ln[-ln(1-α)]～1/T curves for(a)the curing behavior of a modified bismaleimide resin(the values were calculated from the experimental data reported by[28]);(b)the crystallization behavior of an isotactic polypropylene (iPP) melts (the values were calculated from the experimental data reported by [34]).

Fig.4.ln[g(t)]～ln(t) curves for (a) the curing behavior of a Ciba-Geigy epoxy resin (the values were calculated from the experimental data reported by [30]) and (b) the crystallization behavior of a PA6 (the values were calculated from the experimental data reported by [36]).

Fig.5.lnt～1/T curves of(a)the curing behavior of a Ciba-Geigy epoxy resin(the values were calculated from the experimental data reported by[28]),(b)the vulcanization of a rubber (the values were calculated from the experimental data reported by [29]),(c) the curing behavior of a cyanate ester resin (the values were calculated from the experimental data reported by [15]),and (d) the curing behavior of a SH3-DER reactions (the values were calculated from the experimental data reported by [14]).

where h(α) is a function of the conversion degree,which is related to the thermal reaction process of the object under study.The function h(α) is introduced to eliminate the pre exponential factor.In this way,the function h(α) only needs to be fitted when the parameters are solved,and the value of A0is no longer needed to be determined.Taking the natural logarithm of the above formula to get Eq.(21).

Fig.6.h(α)～α curves of(a)the curing behavior of a Ciba-Geigy epoxy resin(the values were calculated from the experimental data reported by[28]),(b)the vulcanization of a rubber(the values were calculated from the experimental data reported by[29]),(c)the curing behavior of a diallyl-bearing epoxy resin(the values were calculated from the experimental data reported by [18]),and (d) the curing behavior of a SH3-DER reactions (the values were calculated from the experimental data reported by [14]).

Fig.7.(a)h(α)～lnα curves of the curing behavior of an E51/DDS epoxy resin(the values were calculated from the experimental data reported by[39]),(b)h(α)～lnα curves of the curing behavior of an o-Cresol formaldehyde resin (the values were calculated from the experimental data reported by [16]),(c) ln[h(α)]～lnα curves of the curing behavior of a cyanate ester resin(the values were calculated from the experimental data reported by[33]),and(d)ln[h(α)]～lnα curves of the curing behavior of a DGEBAmPDA (the values were calculated from the experimental data reported by [17]).

Fig.8.h(α)～α curves of the thermal degradation process of(a)a blood orange juice(the values were calculated from the data reported by[40]),and(b)an ascorbic acid(the values were calculated from the data reported by [41]).

From Eq.(21),it can be found that for a given conversion degree α,since ln [F(α)/A0]is a fixed value,and the curve lnt～1/T should be a straight line with a slope of Ea/R.In fact,this method for calculating activation energy is similar to the model free method in many studies [37].The relevant curve was drawn using the existing data in the literature,as shown in Fig.5,and the corresponding value of Ea/R was obtained.

After the value of Ea/R is obtained,the h(α) curve can be drawn.Thevalues at different temperatures are used as the ordinate and the conversion degree is used as the abscissafor curve fitting,as shown in Fig.6.It can be seen that the h(α)function at different temperatures indeed obeys the same function law.Moreover,it can be well fitted to a variety of polymers and reaction systems,indicating that the method is highly adaptable.

Fig.9.ln[h(α)]～lnα curves of thermal permanent compression deformation of (a) H-8,and (b) H-10 (the values were calculated from the data reported by [42]).

Fig.10.h(α)～α curves of thermal permanent compression deformation of four materials (the values were calculated from the data reported by [43]).

The advantage of this method is that only one fitting parameter(Eaor Ea/R) needs to be determined,which is much easier than needing to fit a series of parameters (such as A,Ea,m,n,k1,and k2) in the literature [38].In addition,for many traditional kinetic models,the thermal reaction conversion degree curves at different temperatures require different parameters to be fitted,which brings a lot of trouble for the use of relevant models for numerical simulation.For example,in the reference literature[14],the Kamal model is used to fit the resin curing process.Table 1 shows the corresponding fitting parameters.It can be seen that there are many fitting parameters to be used,and the fitting parameters at different temperatures are different.However,there are few fitting parameters when using the model in this study,and the established expression is determined,as shown in Table 2,which helps to reduce the workload and is beneficial to the realization of numerical simulation of polymer thermal reaction process.Similar situations can also be found in many studies using traditional models[15–18].Here,a function is introduced to describe the change of internal variables,which correspond to the correction of thermal reaction process that is not suitable for the first-order reaction kinetics hypothesis.Such a function is somewhat equivalent to a series of fitting parameters (such as A,Ea,m,n,k1,and k2).Therefore,the method given in this article only needs to fit the h(α)function in the Fig.6,as the h(α)is the same at different temperatures,which greatly reduce the workload.

Table 1Kinetic parameters for SH3-DER reactions by fitting the experimental data to the Kamal model obtained by Ref.[14].

Table 2Kinetic parameters for SH3-DER reactions by fitting the experimental data to the model in this study

In some cases,as shown in Fig.7(a–d),although the fitted lines of h(α)～α did not overlap,but h(α)～lnα or ln[ h(α)]～lnα can be plotted as the same line.It should be noted that regardless of which form of fitting is adopted,our aim is to obtain the expression of h(α).

After fitting the functional relationship of h(α),the following equation can be obtained by simultaneously solving Eqs.(2),(19)and (21):First,given the degree of conversion,the times required at different thermal reaction temperatures can be obtained according to Eq.(22),which has important application value for the study of thermal aging.

Of course,for a given curing time,the following equation is obtained from Eq.(15) to calculate the degree of conversion:

Fig.11.(a) and (b) the vulcanization behavior of a rubber (the values were calculated from the experimental data reported by [29]);(c) and (d) the curing behavior of a modified bismaleimide resin (the values were calculated from the experimental data reported by [28]).

The following equation can also be used for numerical calculation according to Eq.(16):

Here,it should be noted that the fitting functions of h(α)do not necessarily need to be completely overlapped together,because Eq.(22)will calculate the derivative of the function when establishing the reaction kinetics model,so the constant term in h(α) does not play a role.Therefore,any fitting function with the same change trend at different temperatures can be used to obtain a thermal reaction kinetic model,which is also the difference between this correction method and the method described in Section 3.1.

In addition to the reaction kinetics of polymer crosslinking described previously,it was found that the model is also suitable for some other processes such as thermal degradation reaction(Fig.8),thermal relaxation (Fig.9),and thermal permanent compression deformation(Fig.10).It can be seen that the fitting effect is better,which shows that the method has good adaptability to different thermal reaction types.

3.4.Discussion

It should be noted that the proposed model is not an intrinsic kinetic model.In fact,this research proposes an empirical fitting model and its specific establishment method.The model can characterize the reaction process at different temperatures as a unified expression,which is conducive to the numerical simulation of polymer thermoforming process.The reaction kinetics model is not limited to specific polymer types,and can be well fitted for different polymers and reaction systems.The model styles obtained by this research method are consistent,but the specific parameters will be different,and these specific parameters will be closely related to the corresponding reaction process.

In addition,there are two empirical correction methods to obtain the reaction kinetic model,time-related and conversiondegree-related correction methods.According to the above section,the conversion-degree-related correction method may be more convenient and more adaptable.But on the other hand,from the perspective of polymer molding simulation,the time-related correction method is more suitable,because the transient simulation of polymer thermoforming process is closely related to time,and the model expression established by this correction method is directly given in the form of time,which may be more helpful to use the incremental method to call the model.What’s more,how to choose the appropriate correction method is an interesting issue.In short,when the above correction method is used to calculate the activation energy of the reaction,in the case of a given time or a given conversion degree,if it is found that the activation energy greatly varies at different temperatures,that is to say,it is very difficult to connect these data points into a straight line,as shown in Fig.11(a) &(c),then the adopted method is unable to obtain an accurate model (Fig.11(b) &(d)).This is because the two correction methods are based on the assumption that the function introduced is temperature independent,the difference in activation energy at different temperatures will cause the introduced correction function to be temperature-dependent,that is to say,the correction function at different temperatures is different,so the corresponding correction method is no longer applicable.

He knocked at the door, which was opened by a little old woman who asked, What do you want at this late hour in the midst of this great forest? He said, Haven t you seen a stag about here? Yes, said she, I know the stag well, and as she spoke24 a little dog ran out of the house and began barking and snapping at the stranger

Fig.12 shows the flow chart of the specific steps of using the model in this paper,to facilitate researchers to establish a unified thermal reaction kinetics expression under different thermal reaction temperatures,which is of great significance for the numerical simulation of thermal reaction process.For a certain thermal reaction process,it is first recommended to use Eq.(6)to establish the reaction kinetic expression.Otherwise,it is necessary to determine whether the activation energy is the same at different temperatures for a given time or a given conversion degree,and then select the appropriate empirical correction method for fitting.

Fig.12.Flow chart of the steps to build an expression of thermal reaction kinetics model.

It should be noted that the proposed model is not universal,and it must be established under the conditions that are assumed in the model derivation process.But in any case,the related attempts are not time consuming and the fitting process can be completed in just a few minutes without the aid of computer optimization algorithms.The following research may need to introduce a binary function containing temperature to modify the process,or segment the temperature to reduce the error.

4.Conclusions

Based on the theory of first-order reaction kinetics,we have derived a thermal reaction kinetic model in integral form.On the basis of the model,the effects of the time and conversion degree on the reaction rate parameters(pre-exponential factor and activation energy) were considered,and two types of undetermined functions were used to compensate for the intrinsic variation of the reaction rate.The model was explained and verified using published experimental data of different thermal reaction systems.Compared with the existing thermal reaction kinetic model,the operation method of the model is simple and there are few fitting parameters.It is found that the model shows good fitting and has a very wide adaptability.The model is of great significance to the numerical simulation of polymer thermoforming.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work was supported by the National Key Research and Development Program of China (Grant No.2018YFB2001002).

A0pre-exponential factor,s-1

c concentration of the precursor component,mol.L-1

c0total concentration of precursors,mol.L-1

Eaactivation energy constant,J.mol-1

F(t) primitive function of f(t)

F(α) primitive function of 1/(1 -α)f(α)

f(t) function of time at different temperatures

f(α) function of the conversion degree

g(t) function of time,g(t)=A0F(t)

h(α) function of the conversion degree,h(α)=F(α)/A0

k reaction rate,s-1

R ideal gas state constant,J.mol-1.K-1

T absolute temperature,K

T(t) function of temperature with time,K

t reaction time,s

α conversion degree