Abstract
The recently developed Link-Wise Artificial Compressibility Method (LW-ACM) with double population is presented, validated, and then applied to magnetoconvective flow. The validation step evolves in two stages. First of all, we used the problem of mixed convection over a vertical duct preceded by a sudden expansion in order to validate the multiple relaxation times (MRT) code. Then, we explore the accuracy and the performance of the LW-ACM versus the MRT model and the classic Bhatnagar-Gross-Krook (BGK) model by using a differentially heated cavity problem for Prandtl number . The numerical results for values of Nusselt numbers and process time clearly showed the accuracy and the efficiency of LW-ACM. The adequate accuracy of the proposed model encouraged us to apply it to a magnetoconvection problem. Dynamic and thermal fields are presented and analyzed at and Rayleigh number for Hartmann number up to . Results revealed that the effects of magnetic fields on dynamic and thermal fields are important and that the heat transfer in the cavity is a decreasing function of magnetic field strength.
1. Introduction
The Lattice Boltzmann Method, commonly called LBM, was introduced as a new alternative for the numerical simulation of physical phenomena that can address many problems in physics [1]. Lattice Boltzmann Method sustains today a rapid evolution in terms of physical models, computer implementation, and engineering applications. This method gains an increasing attention and becomes a convenient alternative to classical Computational Fluid Dynamics (CFD) methods. In fact, authors of LBM exceeded the validation step of the method for a wide range of physical problems and many of them are oriented to make the method more efficient than other CFD methods. However, more studies are still needed to make the method useful in engineering practice. The evolutionary process and the formulation of LBM can be highlighted as the following: The Lattice Gas Cellular Automat (LGCA), the continuous Boltzmann-BGK equation, and then the Grad’s Hermite-quadrature expansion. Since its development, LBM demonstrated a successful progress not only in solving viscous flow problems but also in heat transfer area [2]. The first connection of the lattice Boltzmann equation (LBE) to the Boltzmann equation has been established by He and Luo [3]. These formulations can construct models that recover incompressible Navier-Stokes Equations. The Bhatnagar-Gross-Krook (BGK) approximation [4] is the most popular lattice Boltzmann model. This approach is derived from the Enskog equation. In this model the collision operator requires the same relaxation time for each physical quantity. In spite of its simplicity, BGK model suffers from numerical instability at high Reynolds (or Rayleigh number) number. To avoid this restriction, the direct way is to use a large number of grid points; however this will cost large computer resources and lower the computational efficiency. The instability problems of LBM have been extensively studied by many authors through theoretical and numerical investigations. Thus, care was taken to justify the techniques used to stabilize the method without affecting its accuracy.
To remove numerical instability in LBM, some authors [2] used the multiple relaxation times model (MRT). The advantage of this model is that it has adjustable parameters with an additional degree of freedom. These parameters can be determined by optimizing the hydrodynamic properties of the model and linear stability analysis of the LBE evolution operator. Another way to remedy the stability problem in LBM is to return into the entropic formulation of LBM. Some authors [5] used entropic lattice Boltzmann schemes. The derivation of ELB can be performed in many ways and the most popular one is derived from the analog of the discrete Boltzmann function of standard statistical mechanics.
The artificial compressibility method (ACM) introduced by Chorin in 1967 [6] to solve the incompressible Navier-Stokes Equations (NSE) is recently developed with an easy formulation known as Link-Wise Artificial Compressibility Method (LW-ACM) that benefits from some similarities between its classic formulation and Lattice Boltzmann Methods (LBM). Asinari et al. [7] developed the LW-ACM by a finite set of discrete directions (links) on a regular lattice mesh similar to LBM and demonstrated the stability and accuracy of the proposed model. It should be mentioned that the link between LBM and ACM was observed earlier by He et al. [8] and they performed a comparison between the two methods. Their results showed that LBM and ACM are closely related to each other. The similarities between them appear in the continuous form of the macroscopic governing equations, while they differ from each other in their discrete forms. Also they found that difference between LBM and ACM does not affect the momentum results, but it has an impact on the pressure fields. Our recent paper [9] reports numerical results for attached and separated flows by using LW-ACM. The numerical results are relevant for exploiting the applicability of the method to these setups and the computational efficiency of such model is substantiated.
In fact, LW-ACM can be used to solve convective flows. Obrecht and Kuznik [10] used a hybrid thermal LW-ACM to solve the problem of a differentially heated cubic cavity. Their results showed that this scheme remains stable until Rayleigh number . In fact, LW-ACM can be used with double-population model and simulate convective flows and save the memory usage compared to MRT-LBM [11, 12]. Compared to the MRT implementation, the main advantage of the LW-ACM formulation is its simplicity and suitability for parallel implementation. The stability of LBM simulations continues to be the subject of current studies. Recent papers dealing with LBM are oriented towards its efficiency against classical CFD solvers. Thus, much attention has been accorded to the performance and the efficiency of the implementation. The present study deals with the accuracy and the efficiency of the LW-ACM model and then applying it to a magnetoconvective flow problem. The remainder of the present paper is organized as follows. In Section 2, we describe the problem under investigation. The LBM equations with double-population LW-ACM formulation are presented in Section 3. In Section 4, the prescribed LBM model is validated and compared to BGK and MRT results. Then the LW-ACM code is used to compute the magnetoconvection flow in a differentially heated cavity. In particular, the effects of magnetic strength are discussed.
2. Problem Formulation and Boundary Conditions
The problem under investigation is sketched in Figure 1 and detailed as the following. A uniform and horizontal magnetic field is applied to a square cavity filled with an electrically conducting, viscous, and incompressible fluid with negligible viscous dissipation and radiation effects. The cavity of height is heated from below and cooled from above with constant temperatures and , respectively . The upper and lower walls generate a vertical temperature gradient. The problem is considered as two-dimensional. Using the Boussinesq approximation, the fluid density varies as where is the thermal expansion coefficient and is the fluid density at the reference temperature .
The classic thermohydrodynamic equations governing the flow in the cavity with the Boussinesq approximation take the form, , and are the velocity, pressure, and temperature of the fluid, respectively. , , and are the kinematic viscosity, the thermal diffusivity, and the thermal expansion of the fluid. is the unit vector in the vertical direction. is the gravity and is the current density obtained by using Ohm’s Law, and is the electric potential. The effect of the electric field is neglected as it vanishes everywhere in the domain.
The boundary conditions are as follows. For velocities, no-slip boundary conditions are applied on each side of the cavity. For thermal boundary conditions, two types of boundary conditions are used here. The top and bottom walls are considered isothermal and the side walls adiabatic.
In order to express the governing equations in dimensionless form, we introduce the following scaling:By introducing this scaling and previous assumption, we obtain the dimensionless form of (1a)–(1c): and is the Prandtl number . is the Rayleigh number and is the ratio of Lorentz force to the viscous force known as Hartmann number. is the electrical conductivity and is the dynamic viscosity. Hence the problem the problem is fundamentally governed by the prescribed dimensionless parameters and the boundary conditions can be written in dimensionless quantities as follows:
3. Double-Population LW-ACM Method
The artificial compressibility method substitutes the equation of conservation of mass and momentum (1a)-(1b) into the alternative form of the momentum equation (5a) and the artificial compressibility equation (ACE) (5b) as follows: is the Reynolds number based on some reference velocity and reference length . With the dimensionless time , the velocity , the pressure , and the dimensionless external force . is equal to the sum of buoyant and magnetic force. is a reference velocity that can be taken equal to the characteristic velocity in the buoyant flow in our case. Note that by this choice Reynolds number is equal to the square root of Grashof number such that . The artificial density can be determined easily by using the equation of state . The artificial speed of sound is , where is a disposable parameter, analogous to a relaxation parameter, and known as the artificial compressibility. The auxiliary system with (5a)-(5b) converges much more rapidly to the same steady solution as (1a)-(1b). The introduction of the artificial compressibility into the equation of motion does not affect the results and acts in the same manner as the relaxation parameter in LBM.
The explicit time-marching and working on a regular Cartesian grid make some analogies between ACM and LBM. This is the starting point of the link-wise formulation of the ACM. The LW-ACM is developed in the same manner as LBM models starting from the ACM and using analog integration strategy. The particle distributions defined for the finite set of the discrete particle velocity vectors at a site at time are denoted , . In two dimensions, the direction of a single particle probability distribution function is limited to nine directions as follows: for , and for and for . In the LW-ACM, the evolution process of the distribution function can be written aswhere is the discrete speed, is the relaxation time linked to the kinematic viscosity in a lattice unit and given by , and and are the density and velocity in the lattice unit, respectively. is the projection of the body force onto the velocity space and will be detailed later. is the local equilibrium distribution function that has an appropriately prescribed functional dependence on the local hydrodynamic properties as follows:The weighting factors are given by , and . The pressure and the flow velocity in lattice unit are obtained through moment summations in the velocity space as follows:The speed of sound remains constant throughout this study and is equal to in a lattice unit. The asymptotic analysis in [7] showed that the pressure and the velocity in the physical units can be easily computed by taking and , where and is the number of mesh points along the characteristic length of the problem. Quantities in a lattice units and are readily computed in the LW-ACM code.
To solve the advection-diffusion equation for the temperature, we use the D2Q9 model. Note that using the same stencil as for dynamic fields (D2Q9) is more convenient in terms of implementation (see [11]). The thermal evolution equation can be written as follows: is the relaxation time in lattice a unit and is related to the thermal diffusivity and is the equilibrium distribution function defined as follows:Terms of second order in are removed from (10) in order to maintain the linearity of the energy equation with regard to the velocity. The temperature in a lattice unit is computed by conserving only the first moment, i.e., . The relation between the diffusion coefficient in lattice unit and the relaxation time is taken as . With this formulation, the model solves the following energy equation:The introduction of the force term in equation (6a) can be done in various ways. It consists of changing the equilibrium distribution functions by applying the contribution of the force in each direction. The body force giving rise to the external force is defined such thatwhere is the number of nodes in the direction. Note that this assumption is based on the low magnetic Reynolds number (or quasistatic) approximation. Otherwise the magnetic induction model should be taken into consideration to have reliable results [13]. In this case magnetic field should be computed separately using a third distribution function [14].
As for LBM, the treatment of boundary condition in LW-ACM is not applied to the macroscopic quantities but directly applied to the distribution function at the boundary limit. Let us suppose is the fluid node at which we want to compute the unknown distribution function (close to the solid wall ). and are the direction of the discrete velocity and its inverse direction, respectively. The no-slip boundary condition is ensured by reversing the velocity of the particle which wants to enter a solid area:Note that this assumption supposes that the wall is halfway between and . For thermal boundary conditions, isothermal boundaries are applied by using the normal equilibrium condition. To determine unknown densities for straight boundaries, the normal equilibrium condition givesAdiabatic boundary conditions are insured by setting the heat transferred to zero:where the heat flux for insulated surfaces.
By this manner, LW-ACM satisfies the ACM previously introduced by Chorin [6]. For computer implementation, the kinematic viscosity and thermal diffusivity can be easily computed using and and other parameters in simulations:As can be seen from (3a)–(3c), the dynamic and thermal fields depend on the dimensionless parameters, namely, the Prandtl, Rayleigh, and Hartmann numbers. In the classic CFD solvers, this type of problem is fully defined by the prescribed nondimensional parameters . However, for LW-ACM, as for LBM implementation, these three parameters are not enough to compute and since the thermal expansion is known as a physical property but not known in lattice units. For this reason, it is better to add another parameter, namely, the Mach number. Hence, the Mach number is fixed within the incompressible limit . This choice is motivated by the fact that the Mach number adjustment contributes to the instabilities over LBM simulations on the grid scale [15, 16] and it avoids possible velocity overshoots in the ACM calculations as suggested in [6]. In this way, the characteristic velocity in the buoyant flow is taken sufficiently small compared to the speed of sound for all simulations.
The compact form of (6a) and (9a) simplifies the adaptation of an existing BGK-LBM code to LW-ACM. Motivated by this fact, the algorithm of the double-population LW-ACM looks as follows:(i)Setup of dimensionless parameters: , , , and .(ii)Initialization of , , , , and using (7) and (10).(iii)Advection of and .(iv)Compute the equilibrium functions and at fluid nodes using (6b) and (9b).(v)Apply the dynamic and thermal boundary conditions.(vi)Compute and the external force .(vii)Adding the projection of the force term to the distribution function using (12).(viii)Compute the conserved quantities , .(ix)If and reach the convergence criterion, computation terminated.(x)Otherwise, processes (iii)-(viii) are repeated.
Results are assumed to be converged when we reach the following convergence criterion:The subscript sequence represents the space Cartesian coordinates and in the standard basis with the summation is over the entire system and denotes the norm. The volume average Nusselt number which gives an idea about the overall heat transfer in the whole cavity is computed as follows: is the -component of the velocity and the summation is over the entire system.
4. Validation
4.1. Validation of the MRT-LBM Code
As it has been developed in the introduction, LW-ACM is a new formulation of ACM that operates on lattices similar to LBM. Some important points have to be discussed and analyzed before performing simulations using this model, such as its efficiency in terms of speed and accuracy against BGK which is the most widely used model and MRT which is the most stable and accurate one [17].
The MRT-LBM code is validated by two-dimensional simulation of Rayleigh-Bénard convection at different Rayleigh numbers with a Prandtl number of 0.71 (air). The critical Rayleigh number for the onset of the Rayleigh-Bénard convection is which agrees with the theoretical prediction [18]. As the Rayleigh number increased, the steady two-dimensional convection rolls become unstable. The wavy instability and periodic motion observed are in good agreement with the well-known experimental observations and theoretical predictions. Figure 2 shows a snapshot of the temperature distribution for and using MRT-LBM model. At the first level of the simulation, the generation of the mushroom-like isotherms (plumes) is observed. It is caused by the instability and the fluctuation of the instantaneous thickness of temperature boundary layer. This phenomenon can be observed in turbulent convection at high Rayleigh number. The isotherms show a higher level of convective activity and thin thermal boundary layers. It is important to mention here that similar results of snapshot contour plots are found experimentally by Sparrow et al. [19] using electrically heated horizontal copper surfaces situated in a water medium. They used an electrochemical technique to facilitate physical observations. As observed from Figure 2, one can affirm that the MRT-LBM code reproduces the mushroom-like appearance of the isotherms as observed by the reference.
Another validation of the MRT-LBM code was obtained by the mixed convection heat transfer of fluid over a vertical duct preceded with a double-step expansion. The fluid is injected into the channel at a cooler temperature and with a fully developed velocity profile. Channel expansion ratio is taken equal to 30 and the overall length of the computational domain is , where is the step height. The walls downstream of the step are maintained at a constant temperature , while the other walls are treated adiabatically. At the inlet, a fully developed parabolic profile was enforced. With known from the nonequilibrium part of the population normal to the boundary, and can be found as and . At all walls, bounce-back boundary conditions were applied. The artificial open boundary condition is applied in the outlet of the channel by copying the fields from the penultimate column . Both the ‘bounce-back’ and the ‘antibounce-back’ rules [2] were used to implement isothermal and adiabatic boundary conditions, respectively.
The results for mixed convection flow with air have been compared with numerical results reported in using finite-volume method. Computations were performed on grid size. Velocity profiles at the particular axial position Y=5 for fixed Reynolds number and for the Grashof number ranging from to have been compared with those obtained by Tsui and Shu [20] in Figure 3. Here the Reynolds number is defined as , where is the centerline fluid velocity at the inlet. The Grashof number is also based on the step height . It is shown that the numerical results of the MRT-LBM model compare well with the numerical data. Hence, results of MRT-LBM code can be used as a benchmark when testing the LW-ACM. More importantly, it is well known that the MRT-LBM scheme is more accurate, stable, and efficient than the BGK-LBM scheme [17].
Both streamlines and steady state isotherms are shown in Figure 4. As depicted, the asymmetric flow is detected at Grashof number . For , the flow becomes symmetric and the central region is nearly stagnant. However, for Grashof number , flow reversal takes place near the centerline region. Similar results are reported by Tsui and Shu [20] as well.
4.2. Comparison of the LW-AC Method with Other LB Methods
Table 1 shows the computed Nusselt numbers and CPU times for Pr=0.71 obtained by the BGK-LBM, MRT-LBM, and the double LW-ACM method for a square cavity with differentially heated side walls problem at Rayleigh number from 103 to 105. The flow configuration and boundary conditions are as follows. The right and left boundaries are maintained at constant temperatures and , respectively, while the other two walls are kept adiabatic. The no-slip boundary condition is applied at walls. The computing of averaged Nusselt number representing the heat transfer rate at the isothermal walls is defined as follows: where is the cavity height and is the number of lattice nodes in the vertical direction. The results of Wang et al. [2] using MRT-LBM model are used as reference. As it can be seen from all chosen Rayleigh numbers, the LW-ACM model gives results closer to the reference than the BGK-LBM model. On the other hand, the MRT-LBM model gives analogous results as those reported in the reference even we use a coarse mesh. The relative deviation between LW-ACM and the reported results obtained by the coarse 642 mesh is about for , for , and for . As expected, at low Rayleigh number, values of the Nusselt number are close to the reference data while the deviation is more pronounced at higher Rayleigh number. On the other hand, the maximum deviation of the finest grid solution is about at . The double LW-ACM model gives useful results within acceptable errors and it may be used as an effective tool for thermal problems. The table also provides the CPU time required to reach the steady state when using BGK-LBM, MRT-LBM, and LW-ACM models. Clearly, we can see that the CPU time for BGK-LBM model is the shortest. This gives an explanation of why BGK-LBM is the largely used and the most referred model. Except for the first computation using the coarse mesh at , the CPU time required for the MRT model is always high. The comparison between several models with respect to the CPU time shows that the slowness of the MRT model compared to the BGK model is about , , and using the intermediate 1282 grid size and at , , and , respectively. It is true that MRT-LBM model improves the accuracy of simulations, but it requires more time to reach steady state. The double LW-ACM in turn requires about , , and using the same grid size and Rayleigh number.
The second step in the validation process of the LW-ACM code is the comparison of the given results with those reported in the literature using stability analysis. The average Nusselt number for Rayleigh number ranging from to has been compared with those obtained by Clever and Busse [21] in Figure 5. It is shown that our numerical results compare well with both numerical and analytical data obtained from the correlation . At low Rayleigh number, the grid resolution is sufficient to solve the flow and heat transfer. However, for high Rayleigh number, the LW-ACM model slightly underestimates the heat transfer rate.
5. Magnetic Field Effect
In this section, numerical simulations have been performed for square cavity heated from below and cooled from above in the presence of magnetic fields. This can evaluate the double-population LW-ACM in solving more complex flows. Recently, this problem sustains a considerable attention because of a number of its wide variety of applications in engineering and technology such as metallic alloys, magnetoelectrochemistry (MEC), and astrophysical and environmental systems. Thus, it has been shown both experimentally and numerically that magnetic fields can be used to control thermal convection which is important in problems involving metallurgy, microstructure devices, and protein crystals and under reduced gravity conditions. In order to evaluate the potential effects of magnetic fields in heat transfer area, Benos et al. [22] solved analytically and numerically two-dimensional MHD natural convection flow in an internally heated horizontal shallow cavity. Their comparison of the analytical and numerical results showed the validity and the correctness of their analysis. It is of interest to note that the recent progress in superconductivity and electric resistance of metals at very low temperature allowed researchers to reach larger magnetic fields up to 15 T by using super-conducting magnets. Effects of strong magnetic field on two-dimensional natural convection have been studied by Pirmohammadi and Ghassemi [23]. They found that when the Hartmann number is sufficiently large, the convection is suppressed. In fact, in spite of the existence of many papers dealing with the interaction between convection and magnetic fields, there is as yet no complete understanding of magnetoconvection and further research works still needed to provide certain important features. Motivated by this fact, we apply the LW-ACM to the two-dimensional magnetoconvection problem as described in Section 2.
Figure 6 shows the streamlines, the steady state isotherms and the pressure contours at Rayleigh number without magnetic fields and with magnetic fields for Hartmann number 60. We observe that the thermal gradient gives rise to two convection rolls characterizing the instability of Rayleigh-Bénard. The number that governs this type of flow is the Rayleigh number. For all chosen Rayleigh numbers, an adjustment in the flow pattern is observed when magnetic fields are applied. The rotational flow with a single elliptic vortex is observed for high Hartmann number. Isotherms become more uniformly spaced throughout the cavity and the intensity of convection is considerably decreased by the drag induced by the magnetic field, as indicated by a weak distortion of the isothermal lines. The magnetic field reduces the Nusselt number and the fluid velocity. When the Hartmann number is sufficiently large, the mushroom-like isotherms (plumes) caused by the instability of the thermal boundary layer become thinner which tends to homogenize the temperature throughout the cavity. Without magnetic fields, the streamlines form two symmetric vortices with respect to the vertical centerline of the enclosure and when a magnetic field is applied with a sufficiently high value of Hartmann number, the two vortices move from the center of the cavity to the top. The bifurcation of the flow takes place for all cases. The two large vortices appear on the center of the enclosure, but with increasing in Hartmann number, the two vortices are reduced in size located near the top and bottom of the right horizontal top wall. The core of vortex depends on the magnitude of magnetic field. Also, the increase in the value of the magnetic field moves the vortex into the top wall where the temperature is cold. We note also a change in the shape of vortices from the circular shape to elliptical one. One of the application of such result is in melting processes of an electrically conducting fluid metals when the fluid motion presents turbulent convection flow. The applied magnetic field stabilizes or cancels unwanted oscillations in the flow.
Figure 7 presents the evolution of the averaged Nusselt number for different Hartmann numbers. The convergence history indicates that the volume average Nusselt number has a temporal damped oscillation at the beginning of the simulation. This results from incipient instabilities since all the Nusselt numbers have converged to a constant values after a number of iterations. The Nusselt number keeps decreasing with respect of Hartmann number; we can highlight that when , the effect of magnetic fields on the Nusselt number (and so on the heat transferred through the walls) is significant.
For the sake of completeness, computations of the volume-averaged Nusselt number for different magnetic field magnitude give , , , and for , , , and , respectively. In agreement with what has been concluded from the thermal and flow fields in Figure 6, increasing Hartmann number is associated with a decrease in Nusselt number and the heat transfer in the cavity. In fact, increasing magnetic strength leads to an increase the Lorenz force which prevents buoyant force and substantially suppresses the convection. More importantly, the decrease in the heat transfer reaches when increases by compared to the case without magnetic field. When increases from to , the decrease in heat transfer is about .
Results of natural convective flow in the presence of magnetic fields make it possible to judge the ability of the double-population LW-ACM model to predict realistic flows. Thus, the method seems to be able to predict quantities of engineering interest like the volume-averaged Nusselt number accurately. However, as MRT-LBM and LW-ACM methods are comparable, the question is why should double-population LW-ACM be preferred to BGK-LBM with local grid refinement and MRT-LBM?
First, the LW-ACM method is a formulation of ACM that operates on a lattice. Therefore, it does not involve an enforcement of mass conservation at each time step. This causes abrupt changes in particle distributions and breaks the mass conservation especially at the boundary for BGK-LBM. This generates errors near the boundary and destroys uniform convergence to the solution of LBM schemes.
Second, the double-population LW-ACM has the advantage of lower memory usage and computational effort requirement. Thus, if we look carefully at (6a)-(6b) and (9a)-(9b), only the distribution functions and have to be stored whereas, for a MRT-LBM code, it is necessary to add much many quantities at each position since the transformation matrix has eigenvectors connected with the velocity set.
Finally, the implementation of MRT-LBM code is somehow sensitive to the choice of the relaxation parameters which need to be adjusted by a linear stability analysis of the LBE evolution, while only physical quantities governing the problem are needed in order to implement a double LW-ACM code.
6. Conclusions
In this work, a double-population Link-Wise Artificial Compressibility Method (LW-ACM) has been developed to solve convective flows in the presence of a magnetic field. Computation of mixed convective flow over a vertical duct preceded by a sudden expansion was performed in order to validate the MRT model. Then, results of natural convection flows in a differentially heated cavity problem using LW-ACM were compared and validated with those obtained using BGK-LBM and MRT-LBM models and the accuracy of the prescribed model was evaluated. The numerical results for values of Nusselt numbers and process time prove that LW-ACM is an efficient tool for fluid flow and heat transfer computations. Therefore, from a practical point of view, this eliminates the need for using the MRT model which is more expensive and more demanding in terms of memory. Then we applied the numerical model to study magnetoconvection in electrically conducting fluid with , for and for Hartmann number up to . Results showed that the presence of magnetic fields affects directly the buoyant force. The transversal velocity changed considerably which has a direct manifestation on the volume average Nusselt number. This trend demonstrates that LW-ACM is competitive tool and a good alternative to LBM models for engineering applications.
Data Availability
The data used to support the findings of this study are included within the article.
Disclosure
The research was performed as part of the employment of Mohamed Hamdi in the Research and Technology Center of Energy of Borj-Cedria.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
Mohamed Hamdi gratefully acknowledges the assistance of Professor Abdelmajid Jemni from the National Engineering School of Monastir and Professor Samia Abdi-Ben Nasrallah from the Higher Institute of Applied Sciences and Technology of Sousse.