**Coal and Biomass Co-Combustion: CFD Prediction of Velocity Field for Multi-Channel Burner in Cement Rotary Kiln** †

## **Zakia Ngadi \* and Mohamed Lhassan Lahlaouti**

Faculty of Sciences of Tétouan, University Abdelmalek Essaadi, BP. 2121, Tetouan 93002, Morocco; hlahlaouti@hotmail.com

**\*** Correspondence: Zakia.ngadi@gmail.com

† Presented at the 14th International Conference INTER-ENG 2020 Interdisciplinarity in Engineering, Mures, , Romania, 8–9 October 2020.

Published: 11 December 2020

**Abstract:** This paper represents the medialization of alternative fuels co-combustion, in a cement rotary kiln, established on the commercial computational fluid dynamic (CFD) software ANSYS FLUENT. The focus is placed on the key issues in the flow field, mainly on how they are affected by turbulence models and co-processing conditions. Real data, from a Moroccan cement plant, are used for model input. The simulation results have shown a potential effect of the physics model on turbulent and gas-solid flow prediction. The CFD results can be taken as a guideline for improved co-processing burner design and reduce the effect of using alternative fuels.

**Keywords:** CFD; co-combustion; turbulence model; cement rotary kiln

## **1. Introduction**

Alternative fuels (AF) have different chemical and physical properties compared to fossil fuels [1]. Thus, AF particles differ from coal in shape and size, their aerodynamics, and their heating up and combustion mechanism. Co-firing coal/OP (olive pomace) modeling in a cement kiln burner is an intricate process that includes gas and particle phases as well as a turbulence effect on the combustion mechanism and heat transfer.

Although much work [2–6] has been done on AF co-combustion, few researches on CFD (computational fluid dynamic) modeling of co-firing AF under rotary cement kiln burner conditions are available in the literature. Even though the efficiency of different K-epsilon models has been improved in recent years on predicting PF combustion, most improvements have been achieved by minimizing the calculation cost. Nonetheless, it is possible to further improve efficiency by comparing these models to choose a well-predicting model for novel alternative fuel. With this goal, this work seeks to develop a comparison between the three K-epsilon models (standard (STD), RNG (re-normalization group), and RKE (realizable K-epsilon)) combined with the eddy dissipation model, in a 2D modeling case of co-firing coal with olive pomace, under real rotary cement kiln burner conditions. The object is to explore how the flow field is affected by the turbulence models and how the co-processing environment can change the turbulence shape in the rotary kiln. The goal is to have a model that could well predict the turbulence behavior in the kiln environment, could take into account the difference in physical properties of AF versus coal, and could be able to give an overview of the kiln process parameters when changing the co-firing rate.

The remainder of the paper is organized as follows: The proposed mathematical models are presented in Section 2. The computational methodology is discussed in Section 3. Section 4 shows the results of different cases. Finally, Section 5 concludes with a summary.

#### **2. Mathematical Models**

The co-combustion in rotary cement kiln burner startup by carried combustibles particles, from the pulverizer to the furnace, using pressurized air. The flow in the rotary kiln is a gas-solid turbulent flow with chemical reactions. The Eulerian-Lagrangian approach is usually used and recommended for modeling dilute gas-solid flows in pulverized fuels (PF) combustion applications. In the Eulerian-Lagrangian approach, the motion of the continuous phase is modeled using the Eulerian framework, and the motion of dispersed-phase particles is modeled using the Lagrangian framework [7]. In this work, the focus was placed on the comparison of K-epsilon model varieties coupled with turbulence-chemistry interaction model (FR/ED (finite rate/eddy dissipation), i.e., how each model impacts the flow field).

#### *2.1. Gas-Phase Governing Equations*

The time-averaged steady-state Navier–Stokes equations, as well as the mass and energy conservation equations, are solved. The governing equations for the conservations of mass, momentum, energy, and species are given as:

$$\nabla \cdot (\rho \overrightarrow{\boldsymbol{\mu}}) = \mathcal{S}\_{\mathfrak{m}\_{\mathcal{I}}} \tag{1}$$

$$\nabla(\rho \overrightarrow{\boldsymbol{\mu}} \overrightarrow{\boldsymbol{\mu}}) = -\nabla P + \nabla \cdot \left(\overrightarrow{\boldsymbol{\pi}}\right) + \rho \overrightarrow{\boldsymbol{\mathcal{g}}} + \overrightarrow{F}\_{\text{\textquotedblleft}} \tag{2}$$

$$\nabla \cdot \left( \overrightarrow{\boldsymbol{\mu}} \left( \rho \boldsymbol{E} + P \right) \right) = \nabla \cdot \left( \boldsymbol{k}\_{eff} \nabla T - \sum\_{i} \boldsymbol{h}\_{i} \overrightarrow{\boldsymbol{f}}\_{i} + \left( \overrightarrow{\boldsymbol{\tau}} \cdot \overrightarrow{\boldsymbol{\mu}} \right) \right) + \mathsf{S}\_{\mathbf{h}\prime} \tag{3}$$

$$\nabla \cdot (\rho \overrightarrow{\boldsymbol{u}}\_i \boldsymbol{Y}\_i) = -\nabla \cdot \overrightarrow{\boldsymbol{J}}\_i + \mathcal{R}\_i + \mathcal{S}\_{i\boldsymbol{\nu}} \tag{4}$$

The source *Sm* in Equation (1) is the mass added to the continuous phase due to the combustion of particles. Equation (2) is the steady-state momentum equation, where *<sup>P</sup>* is static pressure, <sup>=</sup> τ is the shear forces, and it is called a stress tensor. The third term on the right-hand side (RHS) represents the gravitational body force and the last term is the drag force source term. The first three terms on the RHS of the energy equation in Equation (3) represent the energy transfer due to conduction, species diffusion, and viscous dissipation, respectively. *Sh* is a source term including the heat of chemical reaction and radiation. Radiation heat transfer is calculated using the P1 model, and the absorption coefficient is evaluated by the weighted sum of the gray gas model (WSGG) [8,9]. Equation (4) is the species mass fraction conservation equation for the *i th* species. *Ri* is the chemical reaction production rate of species *i* (this rate will be calculated by the FR/ED turbulence-chemistry interaction models) and *Si* is the source term added from the dispersed phase. <sup>→</sup> *Ji* is the turbulent diffusion flux of species *i* [10].

### *2.2. Particle-Phase Governing Equations*

Alternative fuels are bigger, lighter in density, and are less spherical compared to coal, which affects the motion of particles in the furnace. In the present work, for simplification, we considered AF as a spherical particle; also, the AF and coal particles were modeled separately with two discrete phases following a Rosin-Rammler size distribution [11]. The particles were tracked in Lagrangian frame reference using a stochastic model [12] as shown in Equation (5).

$$\frac{du\_p}{dt} = \frac{3\rho C\_D}{4d\_p\rho\_p} |u\_p - u|(u - u\_p) + g\frac{\rho\_p - \rho}{\rho\_p},\tag{5}$$

where *g* is the gravitational force, and *CD* is the drag coefficient, which is an empirical function of Re as described by Morsi and Alexander [13]. Here Re is the relative Reynolds number, and it's written as,

$$\text{Re} = \frac{\rho d\_p |u\_p - u|}{\mu},\tag{6}$$

In this work, the combustion of solid particle conversion is treated as heating, devolatilizing, and char burning process. For computing volatilization, the single kinetic rate model [14] was used. The kinetic-diffusion limited model [15] was used for char combustion modeling. For both volatilization and char combustion, only convection and heat released by reactions were computed.

#### **3. Numerical Computation and Strategy**

#### *3.1. Geometry and Grid Description*

The rotary kiln used in this paper was assimilated to a real kiln with 46 m in length and 3.8 m in diameter, and it was specially equipped with a multichannel burner (see Figure 1). For simplification, the calculation was started from the fuel injection (i.e., from the burner nozzle end). Therefore, the domain of calculation considered was a 2D axisymmetric rectangle with a multi-entry for fuels and airs. Also, the burner entry geometry was modified for modeling facilities (for more details, see [16]).

**Figure 1.** Scheme of the 2D axisymmetric geometry and the simplified burner details.

Based on the geometric and calculation purposes, a structured non-uniform mesh was used, and a fine grid with 14,040 elements was selected, based on the mesh independence. To optimize convergence, the mesh had a higher density near to the burner and on the axial direction, and was coarse at the rest of the domain.

## *3.2. Numerical Method and Boundary Condition*

The simulation of the 2D domain was implemented on the commercial CFD tool, ANSYS Fluent, based on the finite volume method. A pressure-based coupled algorithm (Coupled) was used to solve the Lagrangian-Eulerian model for gas and particle equations. To achieve numerical stability, the under-relaxation factors were carefully modified for each k-epsilon model. The limited residuals values were taken as 10−<sup>3</sup> for continuity, 10−<sup>5</sup> for energy and radiation, and 10−<sup>4</sup> for all the rest of the parameters.

To be consistent with the real kiln operating condition, the inlet for temperature and velocity was similar to the real case data (see Table 1). A velocity inlet was specified at all the inlet boundaries with direction normal to boundary except for the swirl inlet that had two velocity components inlet. All the walls were treated as adiabatic and no-slip walls. The outlet was set as a pressure outlet condition. The K and epsilon inlet conditions were computed from the turbulent intensity and hydraulic diameter.



## *3.3. Case Study*

As mentioned earlier, this paper discusses the effect of the different K-epsilon turbulence models on the velocity field prediction in a cement rotary kiln under co-combustion parameters. To this aim, six different cases were investigated and are summarized in Table 2.


**Table 2.** Investigated cases details.

#### **4. Results and Discussion**

To investigate the effect of the turbulence model on the flow field, we started with an initial test without the DPM model and species reaction for the three k-epsilon models. Figure 2 shows the axial velocity streamlines without the DPM model and species reaction simulations. One can observe that the STD and RNG models predict the overall flow field similarly, however, the velocity field in RKE simulation extends further into the rotary kiln than STD and RNG, before taking a uniform distribution (see Figure 3).

The three models differ in the capture of the recirculation zone length. The RNG has a wide inner bubble recirculation zone compared to the STD and RKE at around 6.8–10.6 m. The RKE simulation gives a small inner bubble recirculation zone but a longer recirculation circle at around 9.7–14 m. The recirculation zone given by the STD is around 4.5 m. The difference comes from the epsilon equation forms and the solution method.

**Figure 2.** Streamlines of axial velocity without DPM (STD, RNG, and RKE respectively from the top).

**Figure 3.** Axial velocity contours without DPM (STD, RNG, and RKE respectively from the top).

From Figure 4, for the STD and RNG cases, it can be seen that when the pulverized particles mixe into the flow, the uniform velocity field is deteriorated from 13.24 to 9.21 m/s for STD and from 10.8 to 8.58 m/s for RNG. Also, the recirculation zone is smaller than the cases without DPM as shown in Figure 5; this is due to the interaction of particle-flow.

**Figure 4.** Axial velocity contour with particle injection (STD, RNG, and RKE respectively from the top).

**Figure 5.** Streamlines of axial velocity with particle injection (STD, RNG, and RKE respectively from the top).

Generally, particles interact with the flow through the drag forces and thermodynamics effect, causing a reduction in the gas turbulence and an increase in the dissipation rate of that energy which explain the decay of velocity; for more details, see [17]. For the RKE model case, we remark that the uniform velocity is higher in the case with particle injection than the DPM-off case (7.22 m/s vs. 2.54 m/s).

From the results illustrated, it is clear that the RKE model in both test cases gives the lower velocity compared to the RNG and STD models, which may be due to the dissipative nature of the RKE model. Another important point is that the STD and RNG give the lower reverse velocities, as shown in Figure 5, ensuring a speedy mixture of the hot gases with cold burner stream and the stabilization of the flame in the kiln compared to the RKE cases, which impacts the species and temperature distribution that will be discussed in another paper. Current numerical simulations of AF co-combustion are compared with the numerical results obtained by [18,19] who have investigated the simulation of oxy-coal combustion with RANS (Reynolds averaged Navier stokes) and LES (Large eddy simulation) models. The results obtained for the standard and RNG models show the same trend.

#### **5. Conclusions**

The effect of the turbulence model (i.e., STD, RNG, and RKE K-epsilon models) on the prediction of the velocity field of the real rotary kiln case was investigated.

Six different cases were simulated and predicted results were verified with available literature study. Results of the size and the position of the recirculation zone closer to the wall showed the sensitivity of the velocity field to the particle-flow interaction, especially in the RKE cases. A longer recirculation circle at around 9.7–14 m was remarked in RKE case and the recirculation zone given by the STD was around 4.5 m. Nevertheless, the STD and RNG models gave a good estimation of the flow field and recirculation zone size and length, while the RKE model performed best.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Proceedings* **Thermoelectric Generator Based on CuSO4 and Na2SiO3** †

## **Mihail Chira 1,\*, Andreea Hegyi 1, Henriette Szilagyi <sup>1</sup> and Hora¸tiu Vermes, an <sup>2</sup>**


Published: 18 December 2020

**Abstract:** Thermoelectric generators can operate at small temperature differences providing enough electricity for low-power electronics, sensors in distribution networks, and biomedical devices. The article presents the obtaining of a thermoelectric generator and its electrical characteristics using usual substances. Experimental research was carried out using a mixture consisting of several substances (copper sulfate, calcium hydroxide, silicon dioxide, and sodium silicate) in different proportions. The mixture was inserted between two plates, one graphite (hot plate) and the other aluminum (cold plate), thus obtaining a thermoelectric generator. Electrical voltage, output current, and electrical power were measured at different temperatures.

**Keywords:** thermoelectric generators; output current; electrical power; copper sulfate; sodium silicate

## **1. Introduction**

The technology-based age we are in today has brought great advantages for a more comfortable and intelligent life and leads the way for even faster development in all areas of science and engineering, but at the same time brings a very high demand for energy. Current and future energy consumption requires alternative energy sources with low environmental impact. Easily available solar and wind energy have led to the use of fossil fuels as alternative energy sources [1]. Moreover, due to more recent developments in thermoelectric materials and devices (TE), the possibility of efficient use of heat energy, usually wasted, has become a feasible alternative [2]. More importantly, thermal waste has become a very important and environmentally friendly source of wasted energy [3,4]. There are already many studies explaining the mechanism of conversion of heat into electricity [5–8]. Such studies were the starting point for the development and optimization of new materials using the resources of rising nanotechnologies [9,10].

Thermoelectric devices are an attractive solution, due to their passive operation and simplicity and lack ofmoving parts. A thermoelectric generatoris based on the Seebeck effect that generates tensioninmaterials in the presence of a temperature gradient. Although many materials exhibit the Seebeck effect, the energy generation efficiency is based on a complex interaction between temperature, temperature gradient, Seebeck coefficient, and electrical and thermal conductivity [9]. Often, these properties do not extend in the same direction, which leads to challenges in the efficient generation of thermoelectric energy. Altenkirch derived the thermoelectric efficiency, ZT, better known as the thermoelectric figure of merit [11], to describe the thermoelectric energy conversion capacity of a material. This value is given as

$$ZT = \frac{\alpha^2 \sigma T}{\lambda} \tag{1}$$

where α is the Seebeck coefficient, σ is the electrical conductivity, *T* is the absolute temperature, and λ is the thermal conductivity.

The maximum conversion efficiency of thermoelectric energy is given by the relationship:

$$\eta\_{\text{max}} = \frac{T\_h - T\_c}{T\_h} \left[ \frac{\sqrt{1 + ZT\_{\text{avg}}} - 1}{\sqrt{1 + ZT\_{\text{avg}}} + \frac{T\_c}{T\_h}} \right] \tag{2}$$

Materials that offer an intrinsically high ZT value are rare and are generally composed of rare earth-based substances, in particular selenes and telluride-based alloys. SiGe alloys [4], skutterudites [12], clathrates [13], and semi-heusler alloys [14] are other promising systems, among them skutterudite, for example. p-Zn4Sb3, P-CeFe3Co0.5Sb12, and N-CoSb3 are actively studied by numerous research groups, providing a high value of merit in the intermediate temperature range from 500 to 973 K with conversion efficiency of 9.5% [9,15]. One approach to improving the ZT is by nano-structuring, by introducing superlatories [16] or nanostructures [17,18], since it has been recognized that it has a reduction in phononic thermal conductivity. Venkatasubramanian et al. [16] reported that the P-type Bi2Te3/Sb2Te3 superstructures exhibit an exceptionally high ZT of 2.4 at room temperature with a remarkably low thermal conductivity of 0.22 Wm−1K−1. There are other references to Si/Ge [19], GaAs/AlAs [20], and PbTe/PbSe [21] superstructures that have lower thermal conductivity than their alloy counterparts. However, many of these materials have low thermal and chemical stability, difficult to produce in useful forms, and are expensive, with low availability and in some cases high toxicity. There are many studies around the world to find thermoelectric potential alloys and oxides and to create multifunctional properties through microstructural engineering. The expansion of thermoelectric applicability will require the development of low-cost, abundant materials as well as the viability or applicability of scalable manufacturing technologies. A potential source of scalable thermoelectric materials is those based on transition metal oxides, which provide reasonable electrical conductivity and a Seebeck coefficient, while being cost-effective and environmentally friendly. Oxides based on Ti, Mg, W, Zn, Cu, V, Co, Rh, and Mo represent a wide range of materials, less investigated than TE materials [22]. Metal oxides can provide a wide range of electronic properties, from insulators to semiconductors.

TEG are divided into two types, large (or bulk) and micro-TEG. The first category has a millimeter size and provides the output power from several to hundreds of watts in a high heat range. This category is usually used for industrial purposes. The second category operates with low dissipated heat and generates electricity in the range of μW up to a few mW [23].

Waste heat temperature, arrangement of TE modules, enhancing heat transfer co-efficient and parasitic loss plays an important role to maximize the net electric power [24].

In many articles, the thermoelectric material properties are assumed to be constant for easy calculation, but in fact the three properties of thermoelectricity, electrical resistivity, thermal conductivity, and Seebeck coefficient, are temperature dependent. The thermal resistance model, whose properties are assumed to be constant, overestimates the output power and efficiency of the TEG. It can be attributed to the underestimation of electrical resistance and does not solve the temperature distributions within the TEG. Assuming the properties of thermoelectric materials as a constant can induce a numerical deviation at high temperatures [25].

Waste heat recovery is very important in terms of reducing energy costs and environmental impacts. The method used for waste heat recovery must be feasible [26].

#### **2. Materials and Methods**

For the realization of the thermoelectric generator, the following materials and substances were used: graphite plate, aluminum strip, copper sulfate, sodium silicate solution of 30% concentration, calcium hydroxide, and silicon dioxide.

1.2 g CuSO4 with 0.18 g CaOH and 0.08 g SiO2 were mixed until homogenized, after which 2 mL of 30% Na2SiO3 solution was added. A homogeneous paste was obtained, which was applied to a graphite plate of 5 mm thickness, and over it was placed an aluminum strip of thickness 0.5 mm and width 30 mm. The graphite/mixture/aluminum overlap area is 50 <sup>×</sup> 30 mm2.

The obtained thermoelectric generator and the schematic representation of the layers that make up it are shown in Figure 1.

**Figure 1.** Thermoelectric generator made: (**a**) photo; (**b**) schematic representation.

Temperature (Tc and Th) and electrical measurements (voltage, current, and resistivity) were performed using the digital multimeter (Peaktech P4090) from room temperature to 393 K to calculate the Seebeck coefficient, electrical conductivity, and electrical power. The Seebeck coefficient was calculated using the measured voltage generated from the temperature gradient. The temperature difference between the hot and cold part of the sample was generated by heating the graphite plate and cooling the aluminum strip, using an 800 W thermostat hob as a heater and cold air gun as a cooler.

Electrical resistivity was measured by the four-point probe technique, from room temperature to 393 K. The Open circuit was made to measure voltage and current as a temperature function to calculate power.

The thermal conductivity of the TE material was measured using a flash laser analyzer from room temperature up to 393 K.

## **3. Results and Discussion**

The dimensions characterizing the thermoelectric generator are presented in Table 1 and Figure 2.


**Table 1.** Parameters of the thermoelectric generator at 393 K.

Analyzing the data in the table, a high Seebeck coefficient is observed due to the possible p–n junctions formed between the elements of the mixture and a reasonable thermoelectric performance. Further experimental research will make it possible to improve TEG parameters by changing the percentages of substances that make up the mixture.

The Seebeck coefficient allows us to calculate the merit factor ZT and determine the efficiency of the thermoelectric generator η.

The heat source touches the graphite plate, which becomes the positive electrode, and the aluminum plate becomes the negative electrode (Figure 3).

**Figure 2.** Seebeck coefficient of graphite/mixture/aluminum system as a function of temperature.

**Figure 3.** Schematic representation of the mode of operation of the thermoelectric generator.

The temperature of the graphite plate was increased from 297 to 393 K, and at the same time, the voltage and the output current were measured. The graphic representation of the voltage and current variation as a function of temperature can be seen in Figure 4.

**Figure 4.** Thermoelectric generator: (**a**) the voltage; (**b**) the output current.

The output voltage ranges from 375 to 850 mV in the range 297–393 K, and the output current ranges from 1 to 658 μA in the same temperature range.

Analyzing Figure 4, it is found that from 330 K, the system is activated so that the output voltage reaches a maximum value around 850 mV and remains constant, and the output current has a linear increase, the slope of the line being high.

A doubling of the amount of calcium hydroxide in the mixture between the graphite and aluminum plates causes a decrease in the current three times, and a doubling of the amount of silicon dioxide causes a decrease in the current two times.

It was found that the lack of calcium hydroxide in the mixture causes a fluctuation of the generated electric current and voltage.

Serial connection of thermoelectric generators causes the Seebeck voltage to increase, while parallel connections allow the output current to increase. The combination of Serial and parallel connections should lead to an increase in both voltage and electric current.

The output voltage and current were represented according to the number of modules connected in series, respectively, in parallel to a temperature gradient of 96 K in Figure 5.

**Figure 5.** (**a**) The output voltage of the thermoelectric generator for serial and parallel connection; (**b**) the output current of the thermoelectric generator for serial and parallel connection.

When the modules were connected in series, the increase in the output voltage was from 780 mV (one module) to 1224 mV (two modules). The simultaneous decrease in current was measured from 658 μA (single module) to 510 μA (two modules), because the serial connection increased the internal resistance of the thermoelectric generator and as a result decreased the power of the electric current.

To reduce the internal resistance and increase the output current of the thermoelectric generator, the modules were connected in parallel. In this case, the output voltage has changed from 780 mV (single module) to 510 mV (two modules). The simultaneous increase in current was measured from 658 μA (one module) to 1800 μA (two modules), because the parallel connection decreased the internal resistance of the thermoelectric generator and as a result increased the power of the electric current.

For a given heat energy supplied, the temperature difference and the values of the output sizes (voltage, current, electric power) are influenced by the area of the common section—graphite/mixture/ aluminum, by the thickness of the layer forming the mixture, and by the number and mode of connection of the thermoelectric generators. As the joint area increases, the output power also increases. An increase in the thickness of the intermediate layer increases the temperature gradient, but also the internal resistance of the thermoelectric generator. The output power of the thermoelectric generator according to the temperature of the hot plate is shown in Figure 6.

Even if the output power is small, the realization of this mixture of substances with thermoelectric properties represents a step for further research on improving the properties of the thermoelectric generator realized.

The materials used have an advantage over those commonly used, namely, they have thermal and chemical stability, are found in abundance, are cheap, have low toxicity, and are easy to process in order to obtain the thermoelement and low cost.

Figure 7 shows a graph from the technical data sheet of the commercial thermogenerator 1MD03-024-04/1, in which the variation of the electrical voltage and the electrical output power according to the temperature of the hot plate is observed [27].

**Figure 6.** Output power of the thermoelectric generator.

**Figure 7.** Voltage and power output for commercial thermoelectric generator 1MD03-024-04/1.

By connecting the thermoelectric modules in series and in parallel, it is possible to obtain values of electrical voltage and output power comparable to commercial ones.

## **4. Conclusions**

This paper demonstrates by experimental data that the thermoelectric generator made is comparable to those produced industrially and found on the market. The originality of this research is given by the novelty of substances used to generate electric current under the action of the thermal factor.

Experimental results indicated:


**Author Contributions:** Conceptualization: M.C., A.H., H.S., and H.V.; methodology: M.C. and H.V.; validation: M.C., A.H., and H.V.; formal analysis: A.H. and H.S.; investigation: M.C.; data curation: M.C.; writing—original draft preparation: M.C.; writing—review and editing: M.C., A.H., H.S., and H.V.; visualization: A.H. and H.S.; supervision: H.V. and H.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

## **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
