Vacuum residue coking process simulation using molecular-level kinetic model coupled with vapor-liquid phase separation

Zhengyu Chen,Xinhui Yao,Dong Guan,Suoqi Zhao,Linzhou Zhang,Chunming Xu

State Key Laboratory of Heavy Oil Processing,Petroleum Molecular Engineering Center (PMEC),China University of Petroleum,Beijing 102249,China

Keywords:Kinetic modeling Molecular engineering Vapor–liquid equilibrium Vacuum residue coking

ABSTRACT In this work,a molecular-level kinetic model was built to simulate the vacuum residue(VR) coking process in a semi-batch laboratory-scale reaction kettle.A series of reaction rules for heavy oil coking were summarized and formulated based on the free radical reaction mechanism.Then,a large-scale molecularlevel reaction network was automatically generated by applying the reaction rules on the vacuum residue molecules.In order to accurately describe the physical change of each molecule in the reactor,we coupled the molecular-level kinetic model with a vapor–liquid phase separation model.The vapor–liquid phase separation model adopted the Peng-Robinson equation of state to calculate vapor–liquid equilibrium.A separation efficiency coefficient was introduced to represent the mass transfer during the phase separation.We used six sets of experimental data under various reaction conditions to regress the model parameters.The tuned model showed that there was an excellent agreement between the calculated values and experimental data.Moreover,we investigated the effect of reaction temperature and reaction time on the product yields.After a comprehensive evaluation of the reaction temperature and reaction time,the optimal reaction condition for the vacuum residue coking was also obtained.

1.Introduction

Heavy petroleum has high viscosity,great amounts of metals,sulfur,and carbon residue.It is a significant challenge to process heavy petroleum,and the thermal carbon-rejection process,such as the delayed coking process,is a frequently adapted refining technology.Developing a molecular-level kinetic model to precisely simulate the heavy oil thermal conversion chemistry has become a preferable routine to achieve high-efficiency utilization of unconventional oil.

According to the different representation methods of petroleum molecular structure,the molecular-level kinetic modeling of petroleum refining processes can be divided into four various frameworks:bond-electron matrix (BEM) [1],structure-oriented lumping (SOL) [2,3],molecular type and homologous series(MTHS) matrix [4],and structural unit and bond-electron matrix(SU-BEM) framework [5].In the past three decades,a number of researchers have successfully developed molecular-level kinetic models for some common petroleum refining processes based on these frameworks,including fluid catalytic cracking (FCC) [6,7],hydrotreating [8,9],hydrocracking [10],catalytic reforming[11,12],and steam cracking [13,14].However,there are rarely studies on the detailed kinetic model at the molecular level for vacuum residue(VR)coking process.Tian et al.applied the SOL framework to build a molecular-level kinetic model of VR delayed coking.The developed model can calculate the yield and key bulk properties of the product [15].Moreover,they also utilized the model to investigate the effect of the feedstock composition and reaction conditions on the fraction yield [16].To overcome the problem that the complex reaction system was challenging to solve,Horton et al.proposed an attribute reaction model to describe the molecular conversion of VR based on BEM [17].This model used the composition model built by Zhang et al.as input[18],and accurately calculated the product yield and bulk properties.

Previous studies focused primarily on molecular characterization of VR,molecule-based reaction law for VR,and the solving technology for complex reaction systems.However,the VR coking process involves not only chemical changes but also physical changes,which also have a significant impact on the reaction.The coking reaction is a carbon-rejection process in the liquid phase.The carbon–carbon bonds in hydrocarbon chains break and a large number of light fractions are generated.The generated light components vaporize into the vapor phase and the reaction was soon ceased.The vaporization of vapor phase products changes the molecular composition of the reaction system,causing the change of product yield.Some researchers have developed lumped kinetic models for VR thermal cracking with vapor–liquid equilibrium(VLE)[19].The results showed that the added VLE calculation could increase the overall applicability of the reaction model [20].

At present,there are few reports on the molecular-level kinetic model of VR coking combined with vapor–liquid phase separation.In previous work,we proposed a SU-BEM framework to build the molecular-level kinetic model for petroleum refining processes[21,22].In this work,a kinetic model at the molecular level for VR coking was developed.To describe the physical change in the reactor,we proposed a method of coupling the vapor–liquid phase separation model with the molecular-level kinetic model.The developed model showed a great agreement between the calculated and experimental data under different reaction conditions.

2.Experimental

2.1.Feedstock

In this work,Venezuelan vacuum residue was feedstock for the coking reaction and the bulk properties of feedstock were characterized in detail,as shown in Table 1 and Fig.1.The density of VR was relatively high (1.0524 g·cm-3) and sulfur,nitrogen,andoxygen content were also relatively high,reaching values of 4.8%(mass),0.98% (mass),and 1.6% (mass),respectively.The boiling point of 70% (mass) cumulative yield obtained by hightemperature simulated distillation was up to 1000 K.According to the characterization data,Venezuelan VR belongs to an inferior unconventional heavy oil and is suitable for the thermal carbonrejection process.

Table 1 Experimental and calculated properties of VR

Fig.1.Boiling point distribution for the VR feed.

2.2.Experimental equipment and reaction condition

The experiment of VR coking was carried out in a semi-batch laboratory-scale reaction kettle.The schematic diagram of the reaction unit is shown in Fig.2.Before the coking reaction,100 g VR feedstock was placed into the 200 ml reaction kettle.The reactor and pipeline were purged with nitrogen gas and kept the reactor pressure at 101.325 kPa.Then,the reactor was heated and the stirring paddle was turned on.When the first drop of liquid appeared in the liquid collector,the timing started and the stirring paddle was turned off.As the reaction proceeded,vapor phase products were continuously generated.Generated light fractions were swept out of the reactor and collected after cooling.When the set reaction time was reached,the heating was stopped immediately.The gas,liquid,and residue were collected and their yields were calculated.Then,the liquid fraction was cut according to boiling point by a micro-distillation device and the gasoline (IBP–200 °C),diesel (200–350 °C),gas oil (350–500 °C) were obtained.To explore the effect of reaction temperature and reaction time on product yields,we set up seven sets of VR coking experiments under various reaction conditions,as shown in Table 2.These experimental data would be used to tune and validate the coking molecular-level kinetic model.

Table 2 Experimental conditions for VR coking

During the coking experiment,we hope that the laboratoryscale reaction unit can reflect the actual industrial scene as much as possible.However,it is worth noting that there are still some differences between the laboratory-scale reaction kettle and the industrial delayed coking process.The main differences are as follows:

Fig.2.Schematic diagram of the laboratory-scale coking unit.

(1) The dimension of the reactor is different and it will affect the reaction product distribution.For the industrial plant,there is a temperature distribution gradient in the equipment,while the temperature distribution is uniform in the laboratory-scale reaction kettle.In addition,due to the high viscosity of residue,the mass transfer resistance of the industrial plant may be greater than the reaction kettle.The heat and mass transfer will have a significant impact on the reaction results.

(2) The reaction time is different.The vacuum residue(VR)was continuously delivered into the reactor for the industrial delayed coking process.The reaction time in the coke drum is about～20 h,which is much longer than that in the semibatch laboratory-scale reaction kettle.

(3) The recycle ratio is different.For the delayed coking process,the unconverted heavy oil from the bottom of the fractionator will be recycled back to the coke drum to continue the coking reaction.But there is no this process in the laboratory-scale reaction kettle.The circulation of the unconverted VR can improve the conversion rate of feedstock and increase the yield of the light fraction.

3.Micro-Kinetic Model for Vacuum Residue Coking

3.1.Feedstock composition model of vacuum residue

The development of a complex molecule-based conversion model depends on accurate molecular characterization.In previous work [5],we proposed a user-friendly method to represent petroleum molecules and named the SU-BEM framework.This framework combined structural units with BEM and used 34 common group fragments in petroleum to represent petroleum molecules.Engineers can flexibly manipulate these SUs to assemble oil molecules.According to the SU-BEM framework,the molecular composition model of petroleum fraction and the molecular-level kinetic model of petroleum refining processes can be conveniently built.The detailed structural units in the SU-BEM framework were listed in Table 3.

Guan et al.built a molecular composition model of VR based on the SU-BEM framework [23].This work continued to apply this method to develop the composition model of Venezuelan VR.We chose 95 representative model compounds in heavy oil as the core structures,including paraffins,naphthenes,aromatics,sulfur,nitrogen,and oxygen[24].These representative model compounds were assembled by SUs and some typical core attributes in VR were shown in Fig.3.Next,a molecular library of VR can be constructed by extending the sidechain with the carbon number of sidechain and boiling point as constraints.Then,the density,number-average molecular weight(NMW),element analysis,SARA(saturates,aromatics,resins,and asphaltenes) analysis,and simulated distillation were as input to tune the composition model.The bulk properties of the tuned composition model were shown in Fig.1 and Table 1.The results show that there was a good agreement between calculated and experimental values.The developed composition model can accurately capture the physical and chemical attributes of VR and can be used as the input of the molecularlevel kinetic model.

3.2.Reaction rule and reaction network

The reaction rule is utilized to describe the molecular conversion relationship in the complex reaction system.It can realize the digitization of chemical reactions and automatic inference of reactions and products by reactant selection rule,reaction site recognition rule,and product generation rule.In this work,we summarized and formulated 23 reaction rules at the pathway level based on structural characteristics of molecular groups and free radical reaction mechanism.The main reaction rules contain crack-ing reaction,ring opening,dehydrogenation,ring closure,and condensation reaction,as listed in Table 4.

Table 3 Structural units in SU-BEM framework

Fig.3.Typical core attributes in VR composition model.

The cracking reaction is the primary reaction in the VR coking process.It follows the free radical reaction mechanism and the reaction site of bond breaking is usually between the α carbon and β carbon of the free radical [25].Therefore,theoretically,the C-C bond cleavage for paraffins may take place at any position.However,with the growth of carbon number,the bond dissociation energy of the C-C bond in the center will gradually weaker than that of the two ends.Thus,we assumed that if the carbon number of paraffins was more than 12,the C-C bond breaking only took place in the center.Paraffins with long-chain can occur the β scission to produce small paraffins and olefins.The generated olefins can continue to undergo the β scission reaction and smaller olefins can be generated.For aromatics and naphthenes,the reaction site of sidechain cracking is also related to the bond dissociation energy.Some published work has reported that the aromatics sidechain cracking reaction is most likely to take place in two positions.One is between the α carbon and β carbon in the alkyl chain,and the other is between β carbon and γ carbon.After undergoing the sidechain cracking,besides the aromatics with short side chain can be obtained,paraffins and olefins can also be generated,respectively [26].The sidechain cracking reaction of naphthenes is similar to that of aromatics and also has two pathways.One is to generate naphthenes without sidechains and olefins,and the other is to generate naphthenes and paraffins,as shown in Table 4[27].Moreover,the naphthenic ring can also undergo the ring opening to form the naphthene with a sidechain.

Table 4 Reaction rules for VR coking at pathway level

The cracking reaction produces light fractions,while the dehydrogenation reaction can form molecules with a higher double bound equivalent (DBE).Paraffins can be dehydrogenated to produce olefins,followed by undergoing dehydrogenation reaction to generate dienes.The naphthenic ring can also undergo a dehydrogenation reaction to produce aromatics.In this work,according to the difference of the naphthenic ring structure,we divided its dehydrogenation reaction into three types to make the dehydrogenation reaction match the molecular structure.In contrast to the ring opening,the ring with a sidechain can also undergo ring closure to form molecules with higher degrees of condensation.Furthermore,the aromatics and the corresponding heteroatom compounds in the reaction system can undergo condensation reaction to generate polycyclic aromatic hydrocarbons (PAH) and further react to form residue and coke [28].

In addition to dehydrogenation,H radicals can also transfer between different types of molecules in the thermal cracking process[29,30].For example,the most typical reaction is that tetralin or decalin and their derivatives react with olefins as hydrogen donors to produce naphthalene and its derivatives [31].

After reaction rules were determined and formulated,we used the core structures in the composition model to validate the reaction rule one by one.Taking the coking process of C5-phenanthrene as an example,Fig.4(a) shows the automatically generated reaction network based on the reaction rules.The primary reaction pathways are as follow:the C5-phenanthrene undergoes the sidechain cracking reaction to generate short side chain aromatics and paraffins/olefins.Next,these paraffins/olefins can be converted into cyclo-olefins by Diels-Alder reaction,followed by the dehydrogenation reaction to generate aromatics.Then,aromatics undergo condensation reactions with each other to generate a large amount of high-boiling PAH.In addition,according to the free radical reaction scheme,we inferred the partial reaction network of C5-phenanthrene coking at the mechanism level,as shown in Fig.4(b) [32].As shown in the figure,the reaction network of elementary steps is highly matched with the automatically generated reaction network.After all reaction rules were validated,the deterministic reaction network generation algorithm was called to automatically generate a VR coking molecular-level reaction network with 3322 molecules and 7591 reactions.

3.3.Kinetic model

After the reaction network was generated,we can obtain the molecular transformation relationship and the stoichiometric coefficient of each molecule in the reaction.Then,this information can be extracted from the reaction network and the reaction-ratebased expressions of each molecule can be created.It is worth noting that the reversible reactions were divided into two irreversible in this work [15].The obtained reaction rate expressions are as shown in Eq.(1).

VR coking process is a liquid phase reaction system without catalyst.Therefore,the Arrhenius equation can be used to calculate the reaction rate constant,as shown in Eq.(2).

For the reaction system for complex mixtures,there are a large number of reactions.Therefore,a large number of kinetic parameters need to be tuned by the experimental data,which is an extremely challenging task.This work applied the linear free energy relationship (LFER) and quantitative structure–reactivity correlations (QSRCs) to reduce the number of kinetic parameters [33–35].LFER is used to describe the relationship between the reaction enthalpy change and the activation energy,as shown in Eqs.(3)and (4).Eqs.(3) and (4) can calculate the activation energy of exothermic and endothermic,respectively.The reaction enthalpy change in LFER can be estimated by the group contribution method.For the pre-exponential factor,published studies showed that the reaction rate of cracking reaction is related to the carbon number of molecules [27].Thus,the pre-exponential factor can be correlated with the carbon number by the QSRC,as shown in Eq.(5) [22].

Fig.4.Reaction network of C5-phenanthrene coking.(a) automatically generated reaction network;(b) elementary steps in free radical reaction scheme.

Substituting Eqs.(3)–(5) into Eq.(2),we can get the reactionrate-based expressions in this model,as shown in Eqs.(6) and(7).

4.Molecular-Level Kinetic Model Coupled with Separation Model

4.1.Vapor and liquid phase separation

During the coking process,molecules with different boiling points have different reaction time due to the vaporization of light components.To describe the physical change of molecules in the reactor,we developed a vapor–liquid phase separation model for the VR coking process.

For the petroleum system,the equation of state method is usually used to calculate vapor–liquid equilibrium (VLE) [36,37].In this work,we selected the Peng-Robinson (PR) equation of state to calculate the vapor and liquid phase composition.The interaction coefficient of the PR equation is critical to VLE calculation and generally need to be tuned by the experimental data.However,there are thousands of molecules in heavy oil.It is difficult to obtain all interaction coefficients between each molecule.Therefore,this work ignored the interaction coefficients between petroleum molecules and they were assumed to be 0.

The previous research has shown that the VLE model would overestimate the mass vaporization during the heavy oil thermal cracking process[20].It would cause a negative effect on the calculation of the liquid composition.Vacuum residue has high viscosity and density,and mass transfer resistance is relatively large when the light fractions vaporize.Therefore,it is extremely hard for oil molecules to achieve phase equilibrium during the coking reaction.To reduce the effect of mass transfer resistance on the phase separation,we added the separation efficiency coefficient to indicate the deviation between the actual separation result and the VLE.For the vapor phase,the equation of the efficiency coefficient is shown in Eq.(8).

The efficiency coefficient was correlated with the boiling point of the molecule and reaction temperature,as shown in Eqs.(9)–(11).

The parameters in Eqs.(10)and(11)were tuned by experimental data.After that,the efficiency coefficients of each component under various reaction temperature can be calculated by Eqs.(9)–(11).The flowchart of vapor–liquid phase separation for VR is shown in Fig.5.Firstly,the molar fraction and critical properties of feedstock molecules were obtained according to the composition model.Then,the molar fraction and critical properties were input into the VLE model to calculate the molar fraction and bulk properties of the vapor and liquid phase.The obtained vapor phase composition would be substituted into Eq.(8) to calculate the modified vapor phase composition.Finally,according to the mass balance,the modified liquid phase composition can also be obtained.

Fig.5.Flowchart of vapor–liquid phase separation for VR.

4.2.Kinetic model coupled with vapor–liquid phase separation

The coking experiment of vacuum residue was carried out in a reaction kettle,as shown in Fig.2.The feedstock was the liquid phase under the reaction condition.However,as the reaction proceeded,the light fractions with a low boiling point were generated continuously.They were swept out of the reactor and collected after cooling.In order to model this process,we coupled the kinetic model with the phase separation model,as shown in Fig.6.Initially,the reaction time was divided into n nodes and the reaction time of each node(Δt)can be obtained.The feedstock composition model was input into the kinetic model of VR coking and the products of the current node were calculated.Then,the product composition was input into the vapor–liquid phase separation model and the vapor and liquid composition can be obtained.The liquid phase was delivered to the kinetic model to continue the reaction,while the vapor phase was removed from the reaction model.According to the above method,repeated iteration until the reaction was completed.Finally,the vapor phase product would be blended with the liquid and they would be delivered to the split model.The split model would cut the product into gas,gasoline,diesel,gas oil,and residue based on the boiling point.The ranges of each fraction are as follow:

(1) Gas:C1–C4;

(2) Gasoline:C5–200 °C;

(3) Diesel:200–350 °C;

(4) Gas oil:350–500 °C;

(5) Residue and coking:＞500 °C.

5.Results and Discussion

5.1.Kinetic model parameter regression

After the model was developed,we used seven sets of experimental data under different processing conditions to tune and validate the kinetic model.Among the 7 sets of data,conditions (1),and(3)–(7)were used to tune the kinetic parameters and condition(2) was used to validate the generalization of the model.For parameter optimization process,the six sets of experimental data were input into the model and the model can automatically call the genetic algorithm(GA) to tune the model parameters.According to the lower bound and upper bound of parameters,the GA can find a set of relatively optimized parameters.Then,the molar fraction of product molecules can be calculated by this set of parameters.The molecular composition of the product would be delivered to the split model to calculate the fraction yields.In this work,the absolute error of the experimental and the calculated values was used as the value of the objective function (obj),as shown in Eq.(12),and subsequently delivered to the GA.The GA would continue to tune the parameters and obtained a smaller value of the objective function.According to above steps,repeated iteration until the error was minimized and the best parameters can be obtained.The tuned result of parameters can be observed in Fig.7.As can be seen from the figure,there is an excellent agreement between experimental and calculated values for various reaction conditions.

Fig.7.Parity plots comparing the experimental data and calculated data.

5.2.Model validation and bulk properties calculation

Fig.8.Comparison of predicted and experimental data of product yields and bulk properties.(a) Product yield;(b) gasoline composition;(c) distillation curve of diesel;(d)distillation curve of liquid product.

To validate the generalization of the model parameter,we used the tuned kinetic model to predict the fraction yield under reaction condition (2),as shown in Fig.8(a).The result indicates that the kinetic model can predict the product yield accurately.Moreover,the molecular-level kinetic model can not only calculate the yield of each fraction,but also predict the key properties of the product.We also used the tuned model to predict the gasoline composition and the distillation curve of diesel under reaction condition (2).The predicted values were compared with the experimental measurements,as shown in Fig.8 (b) and (c).As can be seen from the figure,a good agreement between the predicted and experimental data was observed.For the gasoline fraction,the coking gasoline had relatively high olefins,about 30% (mass),while the proportion of n-paraffins and naphthenes were relatively low,resulting in around 13% (mass) and 10% (mass),respectively.Because the sidechain cracking,ring opening,and paraffin dehydrogenation were the primary reactions in the coking process,these reactions can generate a large number of olefins based on the free radical reaction mechanism.According to the composition of the coking gasoline,it is necessary to further hydrofining before entering the gasoline pool.Moreover,the predicted distillation curves of liquid fractions were displayed in Fig.8 (d).The results indicate that the ranges of predicted distillation curves for gasoline,diesel,and gas oil were 50–200 °C,200–350 °C,350–500 °C,respectively.It is proved that the boiling point distribution of each fraction obtained by the split model was reasonable.

5.3.Sensitive analysis

The reaction conditions are critical to fraction yields.Therefore,we investigated the effect of reaction temperature and reaction time on the yield of each fraction,as shown in Fig.9.The line is the simulation value and the dot is the experimental data.As shown in Fig.9 (a),with the reaction temperature increasing,the residue and gas oil yields decreased gradually,the gas yield continued to increase,and the yields of the gasoline and diesel had a maximum.

Fig.9.The effect of reaction conditions on the fraction yield.(a) Reaction temperature;(b) reaction time.

For reaction time,a large number of VR molecules were rapidly cracked to generate gas,gasoline,diesel,and gas oil in the first 20 minutes.With the progress of the coking reaction,the reaction rate of the cracking reaction was decreased little by little and the residue yield was also slowly reduced.The results show that residue yield was gradually reduced with the reaction temperature and reaction time rising.However,it is generally believed that the residue and coke are the end cut in the coking process.The yield should be positively correlated with the reaction temperature and reaction time.It is because the coking reaction was carried out in a semi-batch reaction kettle.With the progress of the reaction,the light fractions volatilized into the vapor phase continuously,reducing the liquid phase fraction.However,there were a lot of coking precursors in the light components.The reduction of coke precursors led to the decrease of residue and coke yield.

Besides the residue and coke fraction,the vapor–liquid separation model also affects the yield of the light fraction (such as gas,gasoline,and diesel).According to the theory of parallelsequential reaction,as a middle fraction,the diesel yield should have a maximum,while the yield of gas and gasoline should be rising with the progress of the coking reaction.It is because that the light fraction will continue to undergo the cracking reaction to generate more small molecules.It leads to a continuous increase in gas and gasoline yield and a decrease in diesel yield.When the kinetic model coupled with a vapor–liquid phase separation model,the light fraction vaporized into the vapor phase and the secondary reaction was soon ceased.It causes that the yield of gas,gasoline,and diesel kept a constant value with the progress of the reaction.Therefore,if the model did not combine with the vapor–liquid phase separation,it is difficult to simulate the VR coking process of the semi-batch reactor accurately.

Fig.10.The effect of reaction conditions on the total liquid yield.

In addition to the yield of each fraction,the total liquid yield is also a key index for the heavy oil coking process.The effect of reaction temperature and reaction on the total liquid yield was shown in Fig.10.It can be seen from the figure that the total liquid yield was in the range of 42%–53% (mass) under this range of reaction conditions.With the reaction temperature varied,there was a maximum in the total liquid yield,about 492°C.With the increase of reaction time,the total liquid yield also increased.However,the total liquid yield basically kept a constant value after 60 minutes.Although the increase of reaction time can improve the total liquid yield,the energy consumption of the heater was also growing.Therefore,there is a trade-off between total liquid yield and energy consumption of the heater.After weighing the liquid yield and energy consumption,we recommended that the optimal reaction condition was 492 °C,60 min.

6.Conclusions

In this work,a molecular-level kinetic model was built to simulate the VR coking process in the semi-batch laboratory-scale unit.23 reaction rules for heavy oil coking were summarized and formulated based on the free radical reaction mechanism.Then,A molecular-level reaction network containing 3,322 molecules and 7,591 reactions was automatically generated by applying the reaction rules to VR molecules.For the vapor–liquid two-phase reaction system of the heavy oil coking process,a vapor–liquid phase separation model was proposed.The molecular-level kinetic model for VR coking was coupled with the vapor–liquid phase separation model.We used 6 sets of experimental data under various reaction conditions to tune the model parameters and one set of experimental data to validate the model.The tuned model showed a great agreement between the calculated and experimental data.Furthermore,the sensitivity analysis was carried out to investigate the effect of reaction temperature and reaction time on product yields.The results showed that the total liquid yield had a maximum value with the change of the reaction temperature.After a comprehensive evaluation of the reaction temperature and reaction time,the optimal reaction condition for the VR coking was obtained.

Nomenclature

A QSRC parameters

Ajarrhenius constant of reaction j

CLiconcentrations in the liquid phase of species i,mol·cm-3

Eanactivation energy of reaction family n

E0nactivation energy factor in the Bell-Evans-Polyani LFER of reaction family n

ΔHjenthalpy of reaction of reaction j

kijreaction rate constant of reaction j

R universal gas constant,J·mol-1·K-1

rijreaction rate of reaction j

T reaction temperature,K

X fitted parameter used to calculate efficiency coefficient

Y fitted parameter used to calculate efficiency coefficient

αnreaction index factor in the Bell-Evans-Polyani LFER of reaction family n

φ efficiency coefficient

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 Natural Science Foundation of China (22021004 and U19B2002).