Numerical Study on Bubble Rising in Complex Channels Saturated with Liquid Using a Phase-Field Lattice-Boltzmann Method

: Packed bed reactors have been widely applied in industrial production, such as for catalytic hydrogenation. Numerical simulations are essential for the design and scale-up of packed beds, especially direct numerical simulation (DNS) methods, such as the lattice-Boltzmann method (LBM), which are the focus of future researches. However, the large density di ﬀ erence between gas and liquid in packed beds often leads to numerical instability near phase interface when using LBM. In this paper, a lattice-Boltzmann (LB) model based on di ﬀ use-interface phase-ﬁeld is employed to simulate bubble rising in complex channels saturated with liquid, while the numerical problems caused by large liquid-to-gas density ratio are solved. Among them, the channel boundaries are constructed with regularly arranged circles and semicircles, and the bubbles pass through the channels accompanied by deformation, breakup, and coalescence behaviors. The phase-ﬁeld LB model is found to exhibit good numerical stability and accuracy in handing the problem of the bubbles rising through the high-density liquid. The e ﬀ ects of channel structures, gas-liquid physical properties, and operating conditions on bubble deformation, motion velocity, and drag coe ﬃ cient are simulated in detail. Moreover, different flow patterns are distinguished according to bubble behavior and are found to be associated with channel structure parameters, gravity Reynolds number ( Re Gr ), and Eötvös number ( Eo ).


Introduction
Packed bed reactors have the advantages of simple structure, convenient operation, low operating cost, good heat, and mass transfer performance, and have thus been widely used in industrial processes such as catalytic hydrogenation [1], oxidation reaction [2], nitrification reaction [3], and wastewater treatment [4], among others. For gas-liquid-solid three-phase packed bed reactors, the dynamic behaviors of the dispersed phase directly affect the mixing, mass transfer, and reaction efficiency between gas and liquid, which are extremely important for the design, optimization, and scale-up of the packed beds. However, the detailed dynamics of the dispersed phase and the complicated topological evolution of the gas-liquid interface are difficult to fully understand through experiments or numerical simulations on the reactor scale. The direct numerical simulation (DNS) provides us with a promising solution, in which the lattice-Boltzmann method (LBM) has attracted more and more attention due to its simple computation and clear physical background. At present, LBM has realized the precise tracking of the gas-liquid interface and the accurate prediction of the bubble (dispersed phase) dynamics from the mesoscopic scale [5]. However, the density ratios of liquid to gas in the

Phase-Field LB Model
The phase-field LB model we currently use is proposed by Fakhari et al. [20], in which the lattice-Boltzmann equations (LBE) combine a modified interface tracking equation and improved hydrodynamic evolution equations to achieve numerical computations of isothermal incompressible multiphase fluid systems.

Macroscopic Governing Equations
In the phase-field LB model, a phase-field variable φ is defined and its value is zero in the light fluid (gas bubbles) while one in the heavy fluid (liquid), which changes smoothly across the phase interface between two different fluids. Unlike those phase-field models based on the Cahn-Hilliard equation [13,17,25], a so-called conservative phase-field equation is used in the current model for interface tracking [26]: where t is the time, u is the macroscopic velocity vector, M is the mobility, ξ is the interface thickness, andn represents the unit vector normal to the gas-liquid interface, whose direction is away from the liquid side and points to the gas bubble side, that iŝ The isothermal incompressible multiphase flows are governed by Navier-Stokes equations, including the continuity equation: and the momentum conservation equation, ρ ∂u ∂t + u · ∇u = −∇p + ∇ · µ ∇u + (∇u) T + F s + F b (4) where ρ and µ represent the fluid density and viscosity respectively, p is the macroscopic pressure, F b is the body force, and F s is the surface tension force due to the presence of interface, which can be calculated by [27]: F s = µ φ ∇φ (5) in which µ φ is the chemical potential of the binary-fluid system, which is defined by the derivative of the volumetric free energy E f with respect to the phase field φ: where the volumetric free energy is given by [27,28]: in which the coefficients β = 12σ/ξ and κ = 3σξ/2 are related to the surface tension σ and the interface thickness ξ, and Ψ(φ) is the bulk free energy and expressed as: It should be noted that in the above equations, the gradient and Laplacian of the phase field φ are computed using a second-order isotropic central difference scheme. This scheme avoids the occurrence of a fourth-order derivative when solving the Laplacian in Equation (6) [13,25], and the combination of central and biased differences that was employed in previous studies [11,12], which weakens the numerical dispersion and maintains the conservation of mass and momentum. Therefore, the phase-field LB model is prone to numerical stability at large density ratios.

LBE for Interface Tracking
The following LBE has been recovered into Equation (1) by Geier et al. [29], which can be used to track the gas-liquid interface: where h α is the phase-field distribution function, τ φ is the phase-field relaxation time, e α represents the lattice-related mesoscopic velocity set. For the lattice of D2Q9 model, e α is given by [30]: where c = δx/δt = 1, δx and δt represent the unit lattice length and unit time, respectively, and both δx and δt are set to one in uniform grids. h eq α is the equilibrium phase-field distribution function, which is expressed as: in which In Equations (11) and (12), c s is the lattice sound speed and c s = c/ √ 3, and w α represents the lattice-related weight coefficient set [30], where w 0 = 4/9, w 1-4 = 1/9, w 5-8 = 1/36. Mobility M is positively correlated with the phase-field relaxation time τ φ as: After the calculation of a two-step collision-streaming sequence, the phase field is updated by taking the zeroth moment of the phase-field distribution function: and subsequently the density ρ can be obtained by linear interpolation: where ρ g and ρ l represent the densities of the gas and liquid, respectively.

LBE for Hydrodynamics
An improved hydrodynamic LBE based on HCZ model [9] was developed by Fakhari et al. [17] to update pressure and velocity fields in nearly incompressible multiphase flows: where g α is called the modified hydrodynamic distribution function, and F α is the forcing term, which is calculated by [20]: Ω α is the collision operator, here the collision operator with a multiple-relaxation-time (MRT) model [31] is employed because its performance is more stable in the implementation than the Bhatnagar-Gross-Krook (BGK) model: where g eq α is called the modified equilibrium distribution function and is given by: in which the equilibrium distribution function g eq α is defined as: M is an orthogonal transformation matrix to transform the distribution functions from physical space into moment space [31] and M −1 is its inverse matrix, andŜ is a diagonal relaxation matrix. For the D2Q9 lattice,Ŝ can be selected as: in which s v = 1 τ + 1/2 (22) where τ is the hydrodynamic relaxation time, which is calculated according to the dynamic viscosity derived by linear interpolation: and Among them, µ g and µ l represent the dynamic viscosities of the gas and liquid, respectively. After solving the modified hydrodynamic distribution function according to the two-step calculation sequence of collision-streaming, the hydrodynamic properties are obtained by: Note that the calculation of pressure p requires the updated velocity u, so the update sequence of velocity must precede the pressure.

Discretization
As described in Section 2.1, a second-order isotropic central difference scheme for the phase field is adopted in the phase-field LB model, which is conducive to stable implementation while retaining the second-order accuracy. In detail, the gradient terms used in Equations (5), (17), and (26) and the Laplacian term in Equation (6) are computed according to the following discrete schemes:

Curved Boundary Treatment
If the solid walls are set in the computational domain, the interaction between the fluids and the solid walls can be described by imposing a specified contact angle at the solid boundary [28]: where x w is the position of a point on the solid boundary,n w is a unit vector normal to the solid boundary, with its direction pointing away from the solid wall, and φ w is the phase field value of the point on the solid boundary. Θ is related to the equilibrium contact angle θ and is written as: As shown in Figure 1, the solid (black dots) is presented with a curved solid boundary (black solid line); then Equation (29) is modified to [20]: where φ m is the phase field value of the interpolated point in the fluid (blue dots), and h = |x m − x w | is the distance from the interpolated point to the solid boundary. φ w is estimated by interpolation: Equation (31) can be solved by replacing φ w with Equation (32): As for the neutral wetting conditions (θ = 90 • ), Equation (31) has a trivial solution φ i,j = φ m . To find φ i,j and impose the contact angle conditions on the solid boundary, φ m is the only unknown parameter, which can be obtained by the bidirectional interpolation scheme by utilizing the phase field values of four nearby nodes, and the coordinate positions of nearby nodes are depicted in Figure 1. If the unknown φ i,j is required for interpolation, it is replaced by the φ i,j from the previous time step.
After specifying the wetting boundary conditions, the unknown incoming distribution functions from the solid boundary nodes towards the fluid nodes can be determined by referring to the model of Yu et al. [32]: where The subscript α − denotes the incoming distribution functions, such that e α − = −e α , and the superscript asterisk denotes the pre-streaming or post-collision state of the distribution functions.

Laplace Law
Laplace law states that the pressure difference between the inside and outside of a stationary bubble is proportional to the surface tension and inversely proportional to the bubble radius, which is often used to verify the numerical accuracy of the multiphase LB model [8,11,16,17]. The Laplace law for a two-dimensional bubble is written as: where R b is the radius of the bubble. The computational domain of this verification is a square domain with 201 × 201 grids, and periodic boundary conditions are applied at its four boundaries. At the initial moment of the simulation, the center of the single bubble is located at (101,101), and the bubble is immersed in the liquid phase without any external force. The physical properties of the liquid and bubble are fixed as ρ l /ρ g = 1000 and µ l /µ g = 100. Zheng et al. [16] found that the simulation results were in better agreement with the analytical solutions and the spurious current was smaller when the interface thickness was greater than 4.5 lattice units (lu). Thus, the interface thickness is fixed to 6 lu in this Laplace-law test, and the bubble radius and surface tension are adjusted in the range of 10-50 lu and 0.01-0.1, respectively. Figure 2a,b show the phase field contour of a bubble that has reached a stable state and the distribution curve of the phase field in the radial direction through the bubble, respectively. The analytical curve of the phase field distribution in Figure 2b is derived from the equilibrium phase field profile of diffuse-interface models: where x 0 is the location of the bubble interface. By comparison, it can be seen that the numerical results are almost the same as the analytical solutions. The relationship between pressure difference, surface tension, and radius of the bubble is plotted in Figure 3, which is found to be highly consistent with Laplace law with deviations less than 3% and proves the accuracy of the current phase-field LB model quantitatively.

Bubble Deformation
The rising process of a single bubble under buoyancy action in a rectangular liquid-filled channel is often studied and used to test numerical models [33][34][35][36]. The single bubble will reach a relatively stable state with a constant rising velocity and terminal shape after rising for some time. In our simulations, a static circular bubble with a diameter D b = N x /5 is initially placed at (N x /2, N y /4) in a rectangular computational domain discretized with N x × N y = 512 × 1536 grid cells. A volumetric buoyancy force F b = (ρ − ρ l )G yŷ , where G y is the magnitude of the gravitational acceleration andŷ is a unit vector with a vertical downward direction, is continuously imposed on the bubble. The periodic boundary conditions are used at the left and right boundaries of the computational domain, while the bounce-back boundary schemes are applied at the top and bottom boundaries. The density and viscosity ratios of the liquid and bubble are fixed as ρ l /ρ g = 1000 and µ l /µ g = 100. Besides, there are several dimensionless parameters for characterizing bubble dynamics: (1) The gravity Reynolds number (Re Gr ), (2) The Eötvös number (Eo), (4) The above three are not independent, since Mo = Eo 3 /Re 4 Gr . Table 1 lists the terminal bubble shapes under different Re Gr and Eo obtained by present LBM and other experimental or numerical methods [33][34][35][36][37]. By comparison, the terminal bubble shapes in our simulations can be observed to be highly similar to other research results, which achieves numerical stability under large liquid-to-gas density ratio conditions and demonstrates the qualitative accuracy of the current LB model. During the bubble rising, the Reynolds number of the bubble Re = u g D b ρ l /µ l is continuously recorded. The comparison between the Re obtained from the present LBM and the experimental results is listed in Table 2. The acceptable relative errors further confirm the numerical accuracy of the phase-field LB model.  Figure 4 shows the schematics of two computational domains in 2D Cartesian systems for the bubble rising problem in complex channels, in which the lateral boundaries are set as periodic boundaries, while the bounce-back boundary schemes are adopted at the top and bottom. Among them, the blue circular areas represent the gas bubbles, the red areas are the liquid that fills the channels and surrounds the bubbles, and the white circular and semicircular areas are the artificially constructed bounce-back domains, which are referred to as solid particles and not involved in the iterative computations of phase field and hydrodynamic properties. As depicted in Figure 4, two types of symmetric channels are formed by arranging the solid particles in rectangular and triangular matrices, and they are named wavy vertical channel and S-shaped curved channel, respectively. D b and D p marked in the figures are the bubble and particle diameters, the horizontal spacing L of the particles refers to the shortest distance between two horizontally adjacent particles and similarly the vertical spacing H is the shortest distance between two vertically adjacent particles, and S represents the shortest distance between two diagonally adjacent particles. Hereinafter we use the particle spacing normalized by the bubble diameter to characterize the channel width.

Channel Construction and Numerical Initialization
Initially, a static and stable circular bubble is placed in the middle of the channel and submerged in the liquid with a diameter D b and a volumetric buoyancy force F b = (ρ − ρ l )G yŷ . The phase field in the solid particles is fixed to 0, and in other computational domains the phase field is initialized by the following hyperbolic tangent profile: where R 0 is the initial bubble radius and represents the position of the gas-liquid interface in the radial direction of the bubble, and R(x) is the distance from any position x in the flow domains to the center of the bubble. Unless otherwise specified, in the following simulations, the initial bubble diameter and particle diameter are set as D b = D p = 51 lu. The density and viscosity ratios of the liquid and gas are fixed to ρ l /ρ g = 1000 and µ l /µ g = 100, respectively. Taking into account the balance of numerical stability and accuracy at large density ratios [22], the mobility and interface thickness are taken as M = 0.03 and ξ = 5 lu in the present simulations, respectively. As for the wetting properties, the contact angle of gas-liquid interface on particle surfaces is specified as 40 • . After initialization, the bubble rises through the channel under the action of buoyancy, and simultaneously undergoes a series of evolutions such as deformation, fragmentation, and coalescence. During these processes, as an important dimensionless parameter related to the bubble dynamics, the drag coefficient (C D ) is constantly monitored: where u g is the average velocity of bubbles rising, which is obtained by averaging the instantaneous y-direction velocities over all gas-phase nodes (ρ < ρ l + ρ g /2).

Grid Independence
The channel structure and the effect of the channel on the bubble are greatly affected by the grid, so the case of the bubble rising in a rectangular channel without particles at Re Gr = 100 and Eo = 20 is selected for testing the grid independence. Five different grids with N x × N y = 96 × 1536, 128 × 1536, 176 × 1536, 256 × 1536 and 512 × 1536 grid cells are applied in the bubble rising problem, respectively. The number of grid cells in the y-direction is sufficiently large and constant because its influence on the bubble steady velocity is negligible. The rising velocity variation of the bubble with D b = N x /5 over the dimensionless time is present in Figure 5, and the gravity-based dimensionless time is defined as t * = t G y /D b . It can be observed that when the grid cells are less than 176 × 1536, the bubble terminal velocity is unstable, and the calculation deviation is unacceptable. Thus, the minimum number of grid cells used in the following simulations is set to 176 × 1536. According to the grid independence test results, considering that the channel width is one of the investigated variables affecting the bubble movement, the two computational domains of wavy and S-shaped channels are discretized using 176 × 1536 to 512 × 1536 and 192 × 1536 to 512 × 1536 grid cells, respectively.

Mass Conservation
In order to verify the mass conservation of the current model, the evolution of the total system mass versus the dimensionless time is recorded in Figure 6 during a numerical simulation of bubble rising in a wavy vertical channel with L/D b = 0.719 and H/D b = 0.875 at Re Gr = 100 and Eo = 20. It is observed that the variation of the normalized total system mass M/M 0 is much less than 10 −6 over a long period of time, indicating that the mass of the gas-liquid two-phase system is conserved well using current phase-field LB model.

Channel Width Effect
Undoubtedly, one of the variables that has the most obvious impact on bubble movement in the aforementioned channels is the channel width. Figure 7 shows the evolution curves of bubbles rising velocity and drag coefficient versus the dimensionless time at Re Gr = 100 and Eo = 20. In the wavy vertical channel, the bubble is prone to a stabilized state with a periodically fluctuating terminal velocity, along with C D fluctuating within the certain ranges. In Figure 7, compared with the channel without particles, the rising velocity of the bubble is lower, and the drag coefficient is larger in the channel when the particles are arranged, which indicates that the presence of particles has a significant hindrance to the movement of the bubble.
u g and C D under different channel widths are given in Figure 8, which are calculated by averaging the regularly fluctuating u g and C D over a certain period of time in Figure 7. C D is found to be lower as the channel width becomes wider, and the bubble rises faster. The same results are observed by Patel et al. [38] in the studies on the effects of the amplitude of the sinusoidal channel walls on bubble dynamics. In addition, the horizontal spacing between particles has a more severe impact on the bubble movement than the vertical spacing of the particles. Because the horizontal spacing of the particles directly affects the interaction between the bubble and the particles, and determines the difficulty for the bubble to pass through each gap, while the vertical spacing affects the fluctuation frequency of the bubble drag force, which can only indirectly act on the drag force. In view of this, unless otherwise stated, the following numerical results are obtained from the simulations in the channels with H/D b = 0.875.

Surface Tension Effect
Surface tension is an important parameter to characterize the deformation properties of the gas-liquid interface. As for the bubble, a large surface tension means that the bubble is not easy to deform and break. In our study, since the total bubble mass is well conserved, we count the ratio of the total length of the gas-liquid interface to the area occupied by the bubble to measure the deformation rate of the bubble. As shown in Figure 9, the bubble deformation rate is proved to reduce with the increase of surface tension, which confirms that the bubble with large surface tension tends to maintain its original shape and resist deformation. Furthermore, the declining u g and ascending C D are observed in Figure 10 when the surface tension is gradually increasing. This is because the bubble resists deformation more strongly when colliding with the particles, which restricts the bubble rising velocity. In contrast, the bubble with smaller surface tension is more likely to pass through the narrow gaps in the channel by deforming or even breaking. This phenomenon is consistent with the numerical results of Patel et al. [38], and similar to the studies on the effect of surface tension on liquid penetration by Shi et al. [39].

Bubble Diameter Effect
The effects of bubble diameter on bubble dynamics are depicted in Figure 11. The bubble diameter in this part is not equivalent to the particle diameter; it varies in the range of 25-204 lu, but the particle diameter is still fixed at 51 lu. Figure 11a indicates that it is more difficult for a bubble with larger diameter to pass through the channel due to the greater blocking effects of the particles, so the bubble rising velocity is slowed down, resulting in the ascent of drag coefficient, as demonstrated in Figure 11b.

Driving Force Effect
An additional driving force source term F d can be added to the right side of the governing Equation (4) to simulate the changes in the pressure difference of the packed bed in the actual industry. Here we use the dimensionless ratio F d /(G y ·ρ l ) to measure the magnitude of the driving force.
As illustrated in Figure 12, when the incremental additional driving forces are imposed on the bubble based on the presence of buoyancy force, u g increases linearly while C D decreases rapidly. Because the additional driving force pushes the bubble to pass through the channel faster and more smoothly, the rising velocity is significantly increased, so that u g increases and C D decreases.

Bubble Flow Pattern
During the numerical simulations, we found that the fluid properties and operating conditions set in the numerical initializations have significant impacts on the evolutions of the bubble. In order to classify the different bubble evolution processes, the gravity Reynolds number and Eötvös number that include fluid property and operating condition parameters are used as the classification criteria.
For the wavy vertical channel, mainly four types of bubble flow patterns are identified, and their detailed bubble evolution diagrams are listed in Figure 13. The four flow patterns are named Aw, Bw, Cw, and Dw, respectively; their subscript "w" represents the wavy vertical channel, and the capital letters are used to distinguish different flow patterns.
By observing Figure 13, the flow pattern Aw describes that the bubble is completely blocked by the two particles at the entrance of the channel, and even cannot enter the channel. It is inferred from the small Eo that the bubble surface tension in this flow pattern is relatively large, which makes it difficult for the bubble to enter the channel through deformation when encountering obstacles. In contrast, the bubble in the flow pattern Bw will become elongated when it flows into the narrow gaps thanks to the reduced surface tension, thus passing through the channel integrally. In flow pattern Cw, whenever the bubble crosses the gaps and collides with the particles, some secondary bubbles are split from the tail of the parent bubble. Then these secondary bubbles will gradually merge and expand following the parent bubble, and the coalesced secondary bubble with smaller rising resistance will catch up with and merge into the parent bubble finally. As for the flow pattern Dw, compared to the flow pattern Cw, more secondary bubbles are generated after the parent bubble interacts with the particles due to the extremely small surface tension, and the parent bubble will eventually be split into multiple secondary bubbles and dispersed in the liquid. As a summary, a lower Eo leads to a steady rise of the bubble, while a higher Eo instigates the dynamics instability and causes a multiple breakup in the form of secondary bubbles, which is in good agreement with the numerical results of Patel et al. [38] on bubble dynamics in sinusoidal channels using the level set method.  For the S-shaped curved channel, according to different bubble evolutions, five flow patterns As, Bs, Cs, Ds, and Es are distinguished, as illustrated in Figure 14. Among them, flow pattern As is described as the bubble cannot enter and pass through the channel due to the obstruction of the narrow channel. In flow pattern Bs, the bubble enters the channel by being squeezed and elongated, and it shuttles left and right to rise through the channel. For flow patterns As and Bs, the bubbles retain good integrity because of the higher surface tension. As for the flow pattern Cs, the bubble first splits into two after colliding with the particle, and then the two secondary bubbles will coalesce under their interaction. In another case, the two separated secondary bubbles will rise simultaneously with negligible mutual interference due to their far distance, which is named flow pattern Es. Similar to the flow pattern Dw, the bubble in flow pattern Ds is gradually broken into multiple small secondary bubbles and they are dispersed in the liquid phase, continuously breaking and coalescing. It is worth mentioning that the bubbles in the flow patterns Dw and Ds have better dispersibility and higher specific surface area, which are more conducive to gas-liquid mixing and mass transfer in actual industrial production.  Eo is found to have more significant effects on the flow patterns compared to Re Gr , and lower Eo numbers often correspond to the bubbles with higher integrity because the surface tension that is related to Eo number plays a leading role in bubble deformation. In the S-shaped channel, the flow pattern Ds occupies a wider range than the flow pattern Dw in the wavy channel, indicating that the bubbles are more likely to reach the state of fragmentation and dispersion because of more frequent collisions with particles. It is worth reminding that the trend of the boundaries between different flow patterns in an S-shaped curved channel is similar to the shape regime curve of the single bubble under gravitational motion in the work of Clift et al. [40].
It can be qualitatively summarized from Figure 15 that the bubble becomes more broken and dispersed as Re Gr and Eo increase, and this is attributed to the changes in liquid viscosity and bubble surface tension that have significant impacts on the gas-liquid interface. The increase of Re Gr and Eo leads to the reduction of liquid viscosity and surface tension, respectively. The reduced liquid viscosity causes the weakening of the viscous effects, which results in the increase of bubble wobbling [38], and the lowered surface tension further promotes bubble deformation and breakup.
In addition to the Re Gr and Eo, channel width also has an important influence on the bubble flow patterns. The flow pattern divisions according to the channel widths of two types of channels are demonstrated in Figure 16. The horizontal spacing between particles is found to have more obvious impacts on the flow of bubbles compared to the vertical spacing. The horizontal spacing of the particles directly affects the interaction between the bubble and the particles, thereby affecting the bubble shapes and flow patterns. On the other hand, in the wavy channel, the changes in channel width do not contribute much to the transitions of bubble flow patterns while maintaining the same Re Gr and Eo. Conversely, in the relatively complex S-shaped channel, the transitions of the flow patterns are more sensitive to the variations of channel width; four different flow patterns are revealed by changing the channel width. Under this Re Gr and Eo condition, the flow pattern Cs is more conducive to the even distribution of the dispersed phase in practice because too-small channel width is not prone to the disintegration and flow of the dispersed phase, and too large channel width reduces the contact between the dispersed phase and the particles.

Conclusions
In this paper, the numerical studies on two-dimensional bubble rising in complex channels saturated with liquid at large density ratios within a wide range of gravity Reynolds numbers and Eötvös numbers have been implemented using phase-field LB model. The main research conclusions of this work are as follows: (1) The present LB model is tested through three aspects of Laplace law, bubble deformation, and mass conservation, and it has been proven to have good stability, accuracy, and conservation from both qualitative and quantitative perspectives.
(2) In the simulations of bubble rising in complex channels, the effects of channel width, surface tension, bubble diameter and additional driving force on bubble motion are investigated in detail. The larger channel width and additional driving force as well as smaller bubble diameter and surface tension lead to lower drag coefficients, which are conducive to smooth passage through the channels for the bubble.

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