Constructal Optimization of Rectangular Microchannel Heat Sink with Porous Medium for Entropy Generation Minimization

A model of rectangular microchannel heat sink (MCHS) with porous medium (PM) is developed. Aspect ratio of heat sink (HS) cell and length-width ratio of HS are optimized by numerical simulation method for entropy generation minimization (EGM) according to constructal theory. The effects of inlet Reynolds number (Re) of coolant, heat flux on bottom, porosity and volume proportion of PM on dimensionless entropy generation rate (DEGR) are analyzed. From the results, there are optimal aspect ratios to minimize DEGR. Given the initial condition, DEGR is 33.10% lower than its initial value after the aspect ratio is optimized. With the increase of Re, the optimal aspect ratio declines, and the minimum DEGR drops as well. DEGR gets larger and the optimal aspect ratio remains constant with the increasing of heat flux on bottom. For the different volume proportion of PM, the optimal aspect ratios are diverse, but the minimum DEGR almost stays unchanged. The twice minimized DEGR, which results from aspect ratio and length-width ratio optimized simultaneously, is 10.70% lower than the once minimized DEGR. For a rectangular bottom, a lower DEGR can be reached by choosing the proper direction of fluid flow.


Introduction
Electronic devices are widely used in military and civilian equipment. The components of electronics have been miniaturized and highly integrated in decades. The heat flux of electronics has increased enormously. At present, the heat flux density of some highpower electronic devices has reached 300 W/cm 2 [1], while the local heat flux density of some high-power national defense electronic devices has reached or even exceeded 1000 W/cm 2 [2]. In the foreseeable future, the heat flux density of electronic devices will exceed 2500 w/cm 2 [3]. The heat must be removed in time, otherwise the devices will be destroyed because of the excessive temperature [4,5]. The thermal stress or deformation caused by temperature gradient can also dramatically decrease the reliability of devices [5]. The overheating of electronics has turned into a bottleneck in technological development. The traditional cooling techniques include air cooling, liquid cooling, heat pipe cooling, semiconductor cooling and so on. The cooling effect of air cooling is limited. The heat pipe has good heat dissipation performance, but it has the possibility of natural failure. The semiconductor cooling has the disadvantages of high cost and low efficiency. The traditional cooling techniques are unable to meet with the cooling demand of high heat flux from devices. Therefore, new methods of thermal management must be carried out for electronics. Liquid-cooled MCHSs [6,7] have lots of merits in intensive cooling capacity and are easily packaged with high density. However, due to the small volume and large flow resistance of microchannel, the pressure drop of liquid-cooled microchannel is too high, so it is of great significance to optimize it.
Heat transfer of liquid-cooled MCHS mainly consists of two processes, heat conduction and heat convection, which take place coupled at the same time. Reducing the total irreversibility of the two processes is beneficial for improving the overall thermodynamic performance. Bejan [8] firstly derived the expressions of entropy generation rate (EGR) in the processes of heat transfer and fluid flow, and proposed the principle of EGM [9,10] which enriched the optimization criterion of thermodynamic processes. A lot of researches have made important contributions to the development of EGM [11][12][13][14], and EGM has been widely used for various analyses and optimizations of processes and systems, including heat conduction [15,16], heat convection [17,18], heat exchangers [19], thermodynamic cycles [20,21], and various HSs , and so on. Moghaddam and Saedodin [22] optimized HS with square pin fins using two different algorithms for EGM, and obtained the same optimum parameters of pin fins. Al-Rashed et al. [23] established a cavity HS with circular pin fins, and obtained the numerical influences of length and number of fins on EGM in natural convection. Chauhan et al. [24] numerically compared the relative quantities of thermal EGR and frictional EGR in circular MCHS, and obtained the optimal diameter of HS with the minimum total EGR. According to the principle of bionics, Li et al. [25] established rectangular MCHS models with four kinds of fins in the light of shark-skin, and the optimal shape of fins was obtained for EGM under laminar flow conditions.
Bejan discovered the constructal law and proposed constructal theory [26]. Constructal theory has been extended by many scholars [27][28][29][30][31][32][33][34][35], and been used in heat conduction [35], heat convection [37][38][39], phase transformation [40,41], mass transfer [42], and thermodynamic cycles [43], etc. Constructal theory also has shown distinct superiority in the optimization of liquid-cooled MCHSs. Bello-Ochende et al. [44] numerically conducted the constructal design for rectangular MCHSs under the constraints of a certain volume of HS cell and a range of solid volume proportion. The optimal parameters were achieved for the lowest peak temperature when pressure drop was in the range of 10 kPa to 60 kPa. Xie et al. [45] altered the structure of MCHS by adding bifurcations in channels, which can change the flow conditions of coolant and affect the thermal performance of HS. The results indicated that there was an optimal stage of bifurcations. For triple-layered circular MCHS, constructal design was carried out by Salimpour et al. [46], and the variations of peak temperature with the diameter of every layer were obtained. Lu et al. [47] numerically analyzed the thermodynamic performance of disc HS with Y-shaped microchannel, and obtained the variation laws of peak temperature and pressure drop with branching level.
The combination of EGM and constructal theory is an important research direction in thermal design of electronics. Taking thermoeconomics and EGM into consideration, Badescu et al. [48] optimized heat exchanger according to constructal theory, and obtained the analytic solution between EGR and the parameters of smooth tube heat exchanger. With air as coolant, Chen et al. [49] acquired optimal constructs of cylindrical pin fin HS for EGM under turbulent condition, and Yang et al. [50] obtained optimal constructs with minimum operation cost defined based on EGR. Samal et al. [51] established six HS models embedded circular pipe loops with different shapes, and obtained the optimal shape and channel length of loops. Wang et al. [52] established 3-D heat dissipation model, and obtained the optimal thermal conductivity and heat generating density for EGM under turbulent conditions. Compared with traditional MCHS, PM have two advantages simultaneously when used in MCHS. On the one hand, the specific surface area per volume of PM is larger compared with that of solid fins. On the other hand, inserting PM in the fluid region is equivalent to increase the thermal conductivity of the region. Hence, PM is extensively used in heat convection [53,54]. Bhattacharya et al. [55] analyzed the pressure drop and thermal performance of hybrid HS with PM and solid fins. Compared with the HS with PM only, the experimental results showed that heat transfer coefficient was enhanced significantly with an accessible rise of pressure drop. Hung et al. [56] numerically analyzed the heat transfer coefficient and hydraulic performance of enlarged outlet MCHS, which is fully filled with PM in channels. The results showed that enlarging the outlet can reduce the pressure drop to a certain extent, but the thermal resistance was not necessarily decreased under a low pumping power. In addition, Hung et al. [57] obtained the effect of PM shape on the performance of HS. Leng et al. [58] proposed a new MCHS whose side ribs were replaced by PM entirely. The pressure drop of the new model was nearly 50% lower compared with that of solid side ribs under the same conditions, while the peak temperature increased by 4.7%. Chen et al. [59] altered porosity distribution in the channel of MCHS, and the experimental results showed that thermal resistance was decreased by 62% compared with the conventional MCHS. Gong et al. [60] reduced the thickness of solid side ribs in the MCHSs, and replaced it with PM. The new model has exceedingly thermal and hydraulic performances compared with that of conventional solid rib model. There was an optimal dimensionless thickness of replaced ribs, which can get the best comprehensive performance. Li et al. [61] proposed a new composite HS model with solid pin fins and PM, and the thermal performance was increased by 160% compared with that of pin fin HS. Under the condition of constant wall temperature and fixed volume of PM, Ghorbani et al. [62] analyzed the thermal and hydraulic performances of MCHS partially filled with multi-layer PM, and the optimal positions of PM under different layers were obtained. In numerical assessment of HS filled with PM, Rasam et al. [63] concluded that a lower porosity or a bigger pore density of PM can evidently reduce the total EGR of HS. Wang et al. [64] numerically obtained the optimal parameters for double-layered MCHS with upper porous fins.
Bejan [65] designed the porous media with Darcy seepage flow and obtained the tree-shaped channel structure. Ordonez et al. [66] studied the asymmetric structure of a "volume-point" flow model. Wildi-Tremblay et al. [67] found that the porosity affected the pressure drop absolutely in layered porous media architecture. Lorente [68] studied the eclectrokinetic transfer from the view point of constructal design. And Lorente et al. [69] also investigated the multi-scale heterogeneous porous media. Azoumah et al. [70,71] and Zhou et al. [72] carried out the volume-point optimization for gas-solid reactors, in which heat and mass transfer simultaneously, combining EGM and constructal theory.
According to the constructal law, "For a flow system with structural evolution along the direction of the arrow of time, providing easier and easier paths for the flow flowing through its interior is the fundamental reason for the formation of its structure". The so-called "flow system" and "flow" have broad connotations. The former includes various processes and devices in the engineering field; the latter includes fluid, heat flow, etc. For microchannel heat sink filled with porous medium, fluid passes through the porous medium and heat is also transmitted between the porous medium and the fluid. Obviously, heat sink structure has important influence on fluid flow and heat conduction.
In this paper, the open air cooling model [57] will be transformed into a closed liquid cooling microchannel model, the optimal design with entropy generation minimization for the microchannel heat sink filled with porous medium will be conducted by the constructal design method with wall peak temperature minimization in Ref. [44]. The aspect ratio of HS cell and the length-width ratio of HS will be taken as the design variables under the constraint conditions with fixed cell volume and fixed bottom area. The finite element method will be applied during the numerical simulation. The effects of inlet Reynolds number of coolant, heat flux on bottom, porosity and volume proportion of PM on DEGR will be analyzed. The results can provide theoretical guidelines for the thermal design of MCHS filled with PM.  Figure 1b is selected as the unit computed, and the channel width is Wc. PM is fully filled in the channel. A cell consists of two solid ribs, and the thickness of each rib is W r /2. The volume of cell is:

Geometric Model
The volume of PM in cell is: The aspect ratio of cell is: The volume proportion of PM is: The total bottom area of HS is: The length-width ratio of HS is: The number of cell is: The detail geometric parameters of HS are shown in Table 1. Under the geometric constraints of Equations (1), (2), (4) and (5), the aspect ratio of cell and the length-width ratio of HS will be optimized. Meanwhile, the effects of inlet Reynolds number of coolant, heat flux on bottom, porosity and volume proportion of PM on DEGR will be analyzed, respectively. Table 1. Geometrical parameters of MCHS with PM.

Geometric Parameters Expressions Remarks
Volume V c of the HS cell GL z L x Constant, 0.9 mm 3 , Geometric constraint Bottom area S of the HS L x L y Constant, 100 mm 2 , Geometric constraint

Heat Transfer Model
The following conditions are assumed about the model. The solid ribs and PM are made of isotropic copper with k s = 400 W·m −1 ·K −1 . Deionized water, which is incompressible with constant properties, is selected as the coolant. The flow condition of coolant is single-phase, steady and laminar.
The governing equations of PM regions are: The energy equation for the solid fins: where, V(m · s −1 ) is the velocity vector, p(Pa) is the pressure, and µ(Pa · s) is the dynamic viscosity coefficient, k(W · m −1 · K −1 ) is the thermal conductivity, T( K) is the temperature, c(J · kg −1 · K −1 is the constant pressure specific heat capacity. ε, K(m 2 ) are the porosity and permeability of PM, respectively. C F is the Forchheimer's constant. The subscript eff is the effective value in the simulation. The subscripts s and f represent solid and fluid, respectively. k eff and (ρc) e f f are defined as follows: The boundary conditions are set as follows: Inlet Reynolds number of coolant where, D h (m) is the hydraulic diameter, Re is given as 100, 200, 300, respectively.
Outlet pressure p out = 1atm (17) The interface of solid and liquid in HS where, n is the normal direction of the interface. A uniform heat flux q is applied on the bottom plate of HS The other surfaces of HS are adiabatic: The total pumping power of HS is: The total EGR of heat transfer process is the sum of liquid and solid parts, i.e.,: The EGR of solid region is: where, the subscript thermal represents EGR caused by heat transfer. The EGR of liquid region is: .
where, the subscript frictional represents EGR caused by fluid flow resistance. Φ is the viscous dissipation function per unit volume.
Therefore, the total EGR of the entire HS system is: Bejan number is defined to illustrate the contribution of EGR caused by heat transfer to the total EGR: DEGR are defined as: where, Q = q L x L y . The smaller˜. S g represents the less irreversibility of HS.
In this paper, COMSOL Multiphysics Server 5.5 is used for numerical simulation of the model. In the calculation, the heat sink unit is divided into tetrahedral unstructured grids. The convergence criteria of the continuity equation, momentum and energy equation are 10 −6 . The initial values of the model parameters are set as the following: volume ratio of porous area φ p = 0.4, porosity ε = 0.8, permeability K = 4.73 × 10 −10 m 2 , coolant inlet Re = 100, heat flux density at the bottom of the heat sink q = 100W/cm 2 , initial value of the aspect ratio of the end face of the heat sink α= 1, overall length of the heat sink Aspect ratio β = 1. Table 2 shows the heat sink entropy production rate˜. S g and relative error, calculated in this paper under different grid numbers. In order to obtain more precise results, the third grid strategy will be adopted in this paper.  Figure 2 illustrates the variations of maximum temperature T max , pumping power P with aspect ratio α under the conditions of Re = 100, q = 100W · cm −2 , ε = 0.8, φ p = 0.4, α = 1 and β = 1.  Figure 2 shows that as α increases, T max decreases gradually, but the amount of reduction becomes smaller and smaller. P increases with the increase of α. This is because increasing α increases the convective heat transfer area and strengthens the convective heat transfer of heat sink. However, on the contact surface between the porous region and the solid rib, the flow resistance also increases significantly. That is, the increase of the aspect ratio of the end face of the heat sink unit can strengthen the heat exchange and reduce the maximum temperature of the heat sink, but the price is the significant increase of the input pump power. Figure 3 illustrates the variations of˜. S g ,˜. S g,s ,˜. S g,p with aspect ratio α under the conditions of Re = 100, q = 100W · cm −2 , ε = 0.8, φ p = 0.4, α = 1 and β = 1. Figure 3 shows that there exists a critical value α cr when the value of˜. S g,s is equal to˜. S g,p .˜. S g,s is bigger than˜. S g,p when α < α cr and vice versa. The variations of˜. S g,s and˜. S g,p with α are different. With the increase of α,˜. S g,s diminishes at first and then turns up in a very slight amount. This proves that DEGR of solid fins can not be abated continuously by increasing α. With the increase of α,˜. S g,p first reduces and then goes up, and the increasing tendency becomes larger. The variation rule of˜. S g with α results from the combining of˜. S g,s and . S g,p .˜. S g decreases firstly and then increases significantly with the enlargement of α. There is an optimal α opt = 3 that makes˜. S g reach the minimum value˜. S g,m = 0.003656, and the amount of decrease is 0.001809, which is 33.10% lower than the initial value˜. S g = 0.005464.

Basic Analysis and Optimization
The reason is that˜. S g,s consists of heat transfer across a finite temperature difference. T max decreases as α increases, resulting in˜. S g,s decreasing firstly, and then keeping stable.˜. S g,p consists of heat transfer across a finite temperature difference and fluid flow with viscous effects. T max decreases monotonically as α increases, and EGR of heat transfer decreases firstly and then keeps stable. The channel becomes narrow when α is large, and the flow resistance increases greatly, resulting in a substantial increase in EGR by fluid flow with viscous effects.  Figure 4 shows that Be decreases continually from 0.99 to 0.38 along with the enlargement of α from 1 to 13. The frictional EGR is so pretty that it can be neglected when α = 1. But with the increase of α, the frictional EGR increases significantly and gradually exceeds the value of thermal EGR. This can be interpreted that increasing α benefits the decrease of thermal EGR under the condition with a lower α, but continues to increase frictional EGR. In the practical design of HS, thermal EGR and frictional EGR should both be considered. The thermal EGR should be reduced in every possible way when Be > 0.5, and frictional EGR should be reduced when Be < 0.5. The optimal aspect ratio of cell should be chosen in order to minimize the total DEGR.  Figure 5 illustrates the variations of˜. S g with aspect ratio α for different Re under the conditions of q = 100W · cm −2 , ε = 0.8, φ p = 0.4 and β = 1. Figure 5 shows that˜. S g decreases firstly and then increases with the enlargement of α, and different optimal α opt can be achieved that make˜. S g reach the minimum values˜. S g,m for Re = 100, 200, 300, respectively.

Effects of Reynolds Number on DEGR
The optimal aspect ratio α opt and˜. S g,m all decrease with the increase of Re. For different α, the variation of˜. S g with Re is distinct.˜. S g diminishes with the enlargement of Re in the lower region of α and vice versa. There exists a certain range of α, which makes˜. S g decrease first and then increase with the increase of Re; i.e., there is an optimal Re opt to minimize˜. S g . This is because the increase of Re diminishes thermal EGR and increases frictional EGR in the meantime. In the practical application, the relative variation of DEGRs of heat transfer and fluid friction with inlet Reynolds number of coolant should be comprehensively considered. The geometric parameter design and inlet Reynolds number for HS should be optimized simultaneously to obtain the minimum or near minimum DEGR.  Table 3 shows the optimal aspect ratio α opt of HS cell, the minimum DEGR˜. S g,m , the optimal number N opt of cell and the hydraulic diameter D h for different Re. One can see that as the value of Re increases from 100 to 300, α opt and N opt decrease, and the decreasing proportion of the minimum DEGR is 28.3%. The optimal D h increases, which indicates that a larger hydraulic diameter of HS cell is of advantage to minimize DEGR.  Figure 6 illustrates the temperature and temperature gradient profiles of optimal cell geome tries for different Re. Figure 6a shows that the maximum temperature of optimal cell decreases with the increase of Re. It indicates that increasing Re not only reduces the minimum DEGR, but also reduces the maximum temperature of HS. Figure 6b shows that increasing Re increases the maximum temperature gradient, and temperature gradient of PM region is bigger than that of solid region.  Figure 7 illustrates the variations of˜. S g with aspect ratio α for different q under the conditions of Re = 100, ε = 0.8, φ p = 0.4, β = 1. Figure 7 shows that˜. S g increases with the increases of q for a definite α.˜. S g decreases firstly then increases significantly with the enlargement of α for a given q , and an optimal α opt can be achieved that makes˜. S g get the minimum value˜. S g,m . α opt stay the same value for different q . This is because the increase of q increases thermal EGR directly but has little influence on frictional EGR. When q takes 100W/cm 2 , 150W/cm 2 , and 200W/cm 2 , relative to the initial value under this condition,˜. S g,m decreases by 0.001809, 0.003898, and 0.006399, respectively, and the decrease ranges are 33.10%, 33.71%, and 32.94%, respectively. When α is given,˜. S g increases as q increases. This is because the increase of q will directly cause the increase of heat transfer entropy production rate, and the change law of flow entropy production rate with α is almost not affected by q , so the law of˜. S g changing with the increase of α remains unchanged, and the value of α opt is also constant.  Figure 8 illustrates the 3-D relationship between˜. S g and aspect ratio α, φ p under the conditions of Re = 100, q = 100W · cm −2 , ε = 0.8, β = 1. Figure 8 shows that the variations of˜. S g with α are distinct for different φ p . There exists a critical value φ p,cr that alters the changing rule of˜. S g with α.˜. S g has an ascent after a decline trend with the enlargement of α when φ p < φ p,cr , and an optimal α opt can be achieved that makes˜. S g get the minimum value˜. S g,m .˜. S g monotonically enlarges with the enlargement of α when φ p > φ p,cr .

Effects of Volume Proportion of PM on DEGR
The variations of˜. S g with φ p are distinct for different α. There exists a critical value α cr , which alters the changing rule of˜. S g with φ p .˜. S g monotonically decreases with the increase of φ p when α < α cr . But when α > α cr , optimal φ p,opt can be achieved that make˜. S g get the minimum value˜. S g,m for a certain larger α. The maximum relative difference of the minimum DEGR is 2.2%, which illustrates that the volume proportion of PM has little effect on the minimum DEGR. Figure 9 illustrates the temperature and temperature gradient profiles of optimal cell geometry when φ p are 0.4, 0.6, 0.8 respectively. Figure 9a shows that the maximum temperature of optimal cell increases with the increase of φ p . The reason is that a larger volume proportion of PM directly decreases the effective thermal conductivity of HS, which has negative influence on thermal conduction. Figure 9b shows that increasing volume proportion of PM increases the maximum temperature gradient.   Figure 10 illustrates the variations of˜. S g with aspect ratio α for different ε under the conditions of Re = 100, q = 100W · cm −2 , φ p = 0.4, β = 1. Figure 10 shows that the variation of˜. S g decreases firstly then increases with the enlargement of α for any fixed ε, and an optimal α opt can be achieved that makes˜. S g get the minimum value˜. S g,m .˜. S g,m enlarges and α opt goes down with ε increases. For different α, the variation of˜. S g with ε is distinct.˜. S g diminishes with the enlargement of ε in the lower region of α and vice versa. There exists a certain range of α, which makes˜. S g decrease first and then increases with the increase of ε, i.e., an optimal ε opt can be obtained for minimizing˜. S g . This is because the increase of ε has positive effects on diminishing fluid friction and effective conductivity of PM region.  Table 4 shows the optimal aspect ratio α opt , the optimal number N opt of cell and the minimum DEGR˜. S g,m for different ε. It can be seen from Table 4 that a smaller ε makes a smaller minimum˜. S g,m .  Figure 11 illustrates that the variations of˜. S g ,˜. S g,s ,˜. S g,p with length-width ratio β for different Re under the conditions of q = 100W · cm −2 , ε = 0.8, φ p = 0.4 and α = 3. It can be seen that the variations of˜. S g are distinct for different Re.˜. S g decreases with the increase of β when Re = 100.˜. S g stays almost unchanged with the increase of β when Re = 200.˜. S g increases with the increase of β when Re = 300. This is because that the values of frictional EGR and thermal EGR are diverse for different Re. Figure 11a shows that when Re = 100, . S g,s is bigger than˜. S g,p , and with the increase of β,˜. S g,s diminishes but˜. S g,p stays unchanged. Figure 11b shows that when Re = 200, the variations of˜. S g,s and˜. S g,p with the increase of β are absolutely opposite, i.e.,˜. S g,s increases while˜. S g,p decreases. Figure 11c shows that when Re = 300, the amount of˜. S g,s decreased is less than that of˜. S g,p increased with the increase of β. The reason for the intersection of the˜. S g,s and˜. S g,p curves in Figure 11b,c is that when Re increases from 100 to 300, the maximum temperature of the heat sink decreases that making˜. S g,s decrease, while the flow resistance increases that making˜. S g,p increase.

Aspect Ratio and Length-Width Ratio Are Optimized Simultaneously for EGM
The aspect ratio α and length-width ratio β are optimized simultaneously for EGM. Figure 12 illustrates the variations of˜. S g versus α and β under the conditions of Re = 100, Figure 12 shows that˜. S g descends first and then ascends evidently with the enlargement of α for any fixed β, and an optimal α opt can be achieved that makes˜. S g get the minimum value˜. S g,m for a certain β. For different α, the variations of˜. S g with β are distinct.˜. S g diminishes with the enlargement of β for a certain lower α and vice versa. There exists a certain range of α, which makes˜. S g decrease first and then increase with the increase of β. The twice minimized DEGR˜. S g,mm = 0.003266 is acquired when α = 2.5 and β = 2, and˜. S g,mm is 10.70% lower than the once minimized DEGR (˜. S g,m = 0.003656). This is because that thermal EGR will be increased by increasing α or β alone. The comprehensive influence of α and β on DEGR should be considered at the same time to obtain the minimum DEGR.

Conclusions
In this paper, a constructal optimization model for MCHS fully filled with PM is developed. The aspect ratio of HS cell and the length-width ratio of HS are optimized by numerical simulation methods combining constructal theory and EGM. The effects of inlet Reynolds number of coolant, heat flux on bottom of HS, porosity and volume proportion of PM on DEGR are analyzed, respectively. The following conclusions can be summarized.
(1) There are different aspect ratios to minimize DEGR for different Reynolds numbers. DEGR is 33.10% lower than the initial value through aspect ratio optimization.
With the increase of Reynolds number, the optimal aspect ratio and the minimum DEGR all decrease.
(2) DEGR gets larger with the increasing of heat flux, but the optimal aspect ratios remain constant. The optimal aspect ratios are diverse for different volume proportions of PM, but the maximum relative difference of the minimum DEGR is 2.2%.
(3) The twice minimized DEGR is 10.70% lower than the once minimized DEGR. For a rectangular bottom, the direction of fluid flow has a great effect on DEGR, and a lower DEGR can be reached by setting the right direction of fluid flow.
The model herein is designed for local high heat flux of high-performance computing, smart devices and high-performance CPU, and so on. In practical engineering, it can also be nested with large-scale radiators, which is the next step to design multi-structure integrated radiators with different heat dissipation load areas.
Funding: This work was supported by the National Natural Science Foundation of China (Grant Nos. 51979278 and 51579244).

Acknowledgments:
The authors wish to thank the reviewers for their careful, unbiased and constructive suggestions, which led to this revised manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.