Preprint
Article

Studies on Structures of Turbulent Swirling Coal Flames by Large-Eddy Simulation

Altmetrics

Downloads

86

Views

14

Comments

0

This version is not peer-reviewed

Submitted:

01 October 2023

Posted:

02 October 2023

You are already at the latest version

Alerts
Abstract
Swirling coal flames are encountered in swirl burners of boiler furnaces and cyclone coal combustors, and were studied using large-eddy simulation (LES) in reported references. However, the structures of swirling coal flames remain to be further studied. In the present paper LES of swirling coal flames is conducted for two combustors using two gas combustion models. The statistical simulation results are accessed by experiments and the instantaneous simulation results are used to analyze the specific features of flame structures of four cases. It is seen that different inlet flow conditions and different secondary-air flow rate lead to different flame structures
Keywords: 
Subject: Engineering  -   Energy and Fuel Technology

1. Introduction

Large-eddy simulation (LES) of combustion, including pulverized-coal combustion, is a powerful tool to study the instantaneous flame structures and can give more accurate statistical results than Reynolds-averaged Navier-Stokes (RANS) modeling,. For LES of coal combustion, many researchers did LES of coal jet flames. Kurose et al. [1] did LES of a pulverized-coal-air jet flame using a Smagorinsky SGS stress model and a conservative scalar-fast chemistry gas combustion model, together with coal-particle devolatilization and char combustion models and Lagrangian particle tracking approach. The LES obtained instantaneous particle dispersion and instantaneous temperature map were obtained. The LES statistically averaged gas velocity, temperature, species concentration and particle number density were also reported. However, no vortex structure was shown, and it is not clear about the particle-turbulence interaction in the coal flame. Also, no experimental validation of the velocity, temperature and species concentration fields was reported. Zhou, H [2] reported the so-called “DEM-LES” of coal combustion in a bubbling fluidized bed. For gas-particle flows the SGS kinetic energy stress model together with a DEM particle model is used. The coal devolatilization and char combustion models are taken into account. It seems that no turbulent combustion model was taken into account. As the modeling results, the particle temperature and heating rate and gas species concentration were reported. No instantaneous flow and flame structures were given and there were no experimental validation. Besides, just as in his case of “DEM-LES” of isothermal gas-particle flows, the grid size is too coarse, it looks like a U-RANS modeling, but not LES. Yamamoto et al. [3] studied the coal jet flame by LES using the code Open-FOAM-1.3. The Smagorinsky SGS stress model, DT radiative heat transfer model, Lagrangian particle tracking approach, and a very simple turbulent combustion model based on the concept of the equilibrium concentration and Kolmogorov time scale are used. The LES predicted time-averaged gas temperature and coal burn-out rate were compared with the measurement results. The agreement is good. However, no radial profiles are given. The flame lift-off height predicted by LES is in much better agreement with experiments than that predicted by RANS modeling using the k-ε model. The temperature maps are shown, where the time-averaged temperature maps obtained using LES and RANS are different from each other. However, there is no quantitative comparison between the LES and RANS modeling temperature distribution. Also, no flow vortex structures are reported and no experimental validation of the predicted species concentration field is given. Muto et al.[4] used LES to study the effect of oxygen concentration on NO formation in a coal jet flame. Smith Philip et al. [5] did LES of a coal jet flame using an Eulerian particle-phase modeling, called < Direct Quadrature Method of Moments (DQMOM). The results show that Eulerian-Eulerian LES can capture the detailed flame structures. Wan et al. [6] studied pilot-assisted pulverized-coal flame jet by LES, focusing on the effect of devolatilizatio model. Wen et al. [7] did LES of a coal jet flame adopting a velocity- scalar joint filtered density function model. The statistical results were validated by experiments.Alternative, an extended flamelet/progress variable (EFPV) model for LES of a coal jet flame was reported by Wen et al. [8] and Xing et al. [9].The results show that it is better than the EBU combustion model.
As regard to LES of swirling coal flame, Edge et al. [10] studied two kinds of IFRF oxygen-coal burners using LES and RANS modeling with k-ε and k-ε-RNG turbulence models. The commercial code ANSYS FLUENT version 12 was used. In LES the Smagorinsky-Lilly SGS stress model, an EDC combustion model, and coal devolatilization and char oxidation models were selected. 1.6 million grid nodes with grid sizes of 1–2 mm were taken for the computation domain of near-burner region. The velocity vectors superposed with temperature distribution obtained by LES and RANS modeling were obtained. The authors like to compare and judge the LES and RANS modeling results, but from the predicted temperature (no experimental validation) and wall heat fluxes it is difficult to say that the LES results are better than the RANS modeling results. Gharebaghi et al. [11] did LES of oxygen-coal combustion in a furnace. The commercial code ANSYS-FLUENT was used by selecting the WALE SGS stress model, the EBU (EDM) volatile and CO combustion model, the DO radiative heat transfer model, together with Lagrangian particle tracking approach, coal-particle devolatilization and char oxidation models. The results give the temperature maps in the near-burner region by both LES and RANS modeling using the k-ε model. The comparison of predicted temperature change along the vertical direction from the burner centreline using both LES and RANS modeling with the measurement results was given. Both modeling results remarkably over-predict the temperature and LES results are not better than RANS results. The author’s opinion is that although LES can more accurately model the turbulence, but for this temperature over-prediction it may be caused by the inaccuracy in radiative heat transfer model and char combustion model. No experimental validation of the predicted gas velocity and species concentration was reported. Also no flow vortex structures are reported. Rabacal et al.[12] studied the coal particle history during combustion in a large-scale furnace using large eddy simulation. A specific feature of the devolatilization model was discussed. Watanabe et al. [13] did large-eddy simulation of coal combustion in a pulverized coal combustion furnace with a complex burner. The simulated gas temperature and NO concentrations are in general agreement with the experiments. Rieth et al [14] introduced the flamelet model into large eddy simulation of realistic coal furnaces. The results of the LES are focussied on instantaneous particle and gas phase data to gain additional insight into the coal conversion process inside the furnace. Shen et al. [15] studied coal combustion in a 660 MW utility boiler furnace using LES. However, the adopted grid sizes for such a large-size furnace are questionable for the requirement of LES. Olenik et al. [16] studied swirl-stabilised pulverised coal combustion in the semi-industrial (2.5 MW) IFRF furnace by means of large eddy simulation with a EBU combustion model and compared to results from the experimental results and RANS predictions. It was found that some temperature underprediction was observed. Franchetti et al. [17] studied a 100 kW(th) swirling oxy-coal furnace using LES. It was reported that in the downstream region the LES overestimated the combustion rates.
It is seen that in most of above-cited studies the focus was paid to the time-averaged LES results and the application of LES in practical coal furnaces. The instantaneous structures of swirling coal flames remain to be further studied. In this paper, swirling coal flames are studied by LES. The specific features of instantaneous strucutures of swirling coal flames are analyzed and discussed.

2. Filtered Governing Equations for LES of Coal Combustion

The filtered gas governing equations for Eulerian-Lagrangian LES of coal combustion can be given as:
ρ t + x i ρ u ¯ i = S ¯ S ¯ = n ¯ k m ˙ k
t ρ u ¯ i + x j ρ u ¯ i u ¯ j = x j μ u ¯ i x j p ¯ x i τ i j x j + ρ k τ r k ( u ¯ k i u ¯ i ) + u ¯ i S
t ( ρ f ¯ ) + x j ( ρ u ¯ j f ¯ ) = x j ( μ S c f ¯ x j ) m j s g s x j + α s S
ρ h ¯ t + x j ρ h ¯ u ¯ j = x j μ Pr h ¯ x j q s g s , j x j + n k Q + k c p T ¯ S
ρ Y ¯ s t + x j ρ u ¯ j Y ¯ s = x j μ S c Y ¯ s x j w ¯ s w s g s g j s g s x j + α ¯ s S ¯

3. Sub-Grid Scale (SGS) Stress Models

For single-phase flows, frequently used SGS stress models are the Smagorinsky eddy-viscosity model [18], Germano´s dynamic eddy-viscosity model [19] and Kim´s SGS energy-equation model [20]. The SGS stress in the momentum equation is defined by:
τ i j ρ u i u j ¯ ρ u ¯ i u ¯ j
The Smagorinsky eddy viscosity model is given by
Preprints 86640 i001
where Cs=0.16 is the Smagorinsky constant, Δ is the filtered scale. The shortcoming of the Smagorinsky model is its too large dissipation rate of the kinetic energy. Germano proposed a dynamic Smagorinsky eddy viscosity model, in which the coefficient Cs is not a constant, but is determined by the following filtration process as
Preprints 86640 i002
M i j = 2 ρ C s Δ ¯ 2 S ¯ ^ S ¯ ^ i j 2 ρ C s Δ ¯ ^ 2 S ¯ ^ S ¯ ^ i j L i j = 2 ρ C s Δ ¯ ^ 2 S ¯ ^ S ¯ ^ i j 2 ρ C s Δ ¯ 2 S ¯ ^ S ¯ ^ i j
It was indicated by many investigators that in many cases the dynamic eddy-viscosity model gives better results than that given by the standard Smagorinsky model. Kim proposed a SGS kinetic energy model as
τ i j = 2 ρ v t ( S i j 1 3 S k k δ i j ) + 2 3 ρ k s g s δ i j
ρ ¯ k s g s t + x i ρ ¯ u i k s g s = P s g s D s g s + x i ρ ¯ ν t Pr t k s g s x i + W s .
k s g s = 1 2 ( u i u j ¯ u ¯ i u ¯ j )             v t = C v ( k s g s ) 1 / 2 Δ
Preprints 86640 i003
The merit of the SGS kinetic energy model is that it avoids the negative value of eddy viscosity; hence it is better than the eddy-viscosity models. In application it was reported that the results of SGS kinetic energy model are near to those obtained using the dynamic eddy-viscosity model. In this paper, a combined
Smagorinsky-Lilly’s SGS stress model [21] is used as
τ i j 1 3 τ k k δ i j = 2 μ t S ¯ i j ,   S ¯ i j 1 2 u ¯ i x j + u ¯ j x i ,   μ t = ρ L s 2 S ¯ ,   S ¯ 2 S ¯ i j S ¯ i j L s = min ( κ d , C s V 1 / 3 )
where: δ i j -Kronecker Delta,V -cell volume, L s -SGS scale, S i j -stress deformation rate, μ t -SGS viscosity coefficient, d -closest distance to the walls, κ -Von Kármán constant, C s -empirical constant, taken as C s =0.1。
The SGS mass flux and heat flux are closed by gradient modeling as
g s g s , j = ρ u j Y s ¯ u ¯ j Y ¯ s = μ t σ Y Y ¯ s x j
q s g s , j = ρ u j T ¯ u ¯ j T ¯ = μ t σ T T ¯ x j
where σ T and σ Y are model constants

4. Combustion Models for LES of Coal Combustion

For gas combustion (CO and volatile combustion) during coal combustion, two models were adopted. One is the Magnusen’s modified eddy-break-up (EBU) model [22]
w ¯ s w s g s = c 1 ρ S ¯ min Y f ¯ , Y o x ¯ β , c 2 Y P ¯ 1 + β
where c 1 =4.0, c 2 =0.5
The gas-phase combustion includes two reactions:
C H 4 + 2 O 2 C O 2 + 2 H 2 O       ( reaction 1 )
2 C O + O 2 2 C O 2       ( reaction 2 )
Alternatively, a presumed-PDF fast chemistry gas combustion model [23] is used. In the latter case a mixture fraction equation should be solved
t ( ρ f ¯ ) + x j ( ρ u ¯ j f ¯ ) = x j ( μ S c f ¯ x j ) m j s g s x j
where m j s g s = ρ u j f ¯ u ¯ j f ¯ = μ t σ f f ¯ x j and f f ¯ f ¯ f ¯ = 1 2 L s 2 f ¯ 2
For radiative heat transfer the P1 model is used as:
q ¯ r = a G 4 a σ T ¯ 4
For particle combustion, neglecting moisture evaporation and taking the sub-models of two-equation devolatilization and diffusion-kinetic char oxidation model with total formation of CO are used [22].

5. Particle-Phase Equations

The Lagrangian particle continuity, momentum and energy equations are [23]
A n k v k n d A = N k = c o n s t
d u k i d t k = 1 τ r k ( u i u k i )
d T k d t k = [ Q h Q k Q r k + m ˙ k ( c p T c k T k ) ] / ( m k c k )

6. Simulated Pulverized-Coal Swirl Combustor and Numerical Method

The LES of swirling pulverized-coal combustion was made for a combustor, shown in Figure 1.
Figure 1. A coal combustors with swirl burner.
Figure 1. A coal combustors with swirl burner.
Preprints 86640 g001
The swirl combustor was measured at Shanghai Jiaotong University (Figure 1a) [24]. In this combustor, non-swirling primary air laden with coal particles was issued from the inner annular inlet. The secondary swirl air was issued from the outer annular inlet. The flow parameters of Cases 1 and 2 are given in Table 1. The flow parameters of Case 3 are the same as those of Case 2 with only a difference in adding a inlet perturbation using a vortex method.
Table 1. Flow Parameters for Cases 1 and 2.
Table 1. Flow Parameters for Cases 1 and 2.
Primary air flow laden with coal particles Secondary swirling air flow
Axial velocity(m/s) Temperature(K) Coal feeding rate(kg/h) Axial velocity(m/s) Temperature(K) Swirl number
Case 1 23.1 353 14.0 38.3 573 1.0
Case 2 28 353 16.4 45.1 573 1.03
The proximate and element analyses of coal are given in Table 2 and Table 3
Table 2. Coal proximate analysis (Weight%).
Table 2. Coal proximate analysis (Weight%).
moisture volatile fixed carbon ash
6.3 35.8 53.7 4.2
Table 3. Coal element analysis(Weight%).
Table 3. Coal element analysis(Weight%).
carbon hydrogen nitrogen sulfur oxygen
72.6 5.05 1.29 1.55 15.31
The particle size distribution are shown in Table 4
Table 4. Particle size distribution.
Table 4. Particle size distribution.
<200 μ m <75 μ m <40 μ m <25 μ m <10 μ m
100 80 60 40 15
Kinetic constants of coal devolatilization:Bv1=2e+5, Ev1=1.046e+8j/kgmol,Bv2=1.3e+7,Ev2=1.674e+8 j/kg-mol。Kinetic constants of char reaction: Bh=0.86kg/m2.s (N/m2),Eh=1.02 × 108J/kg-mol。Gross Calorific Value:29.29MJ/kg。
The grid sizes are taken as 0.1mm to 1mm. The total grid number is 1943214. The contribution of the resolved part of the total kinetic energy is 75% in most of the computational region. The effect of grid size has been tested and it is found that the statistical results are independent of the grid sizes smaller than the adopted ones. The time step is taken as 1e-3s, which is chosen based on the computation test. Within each time step the convergence can be reached after 30 to 70 iterations. The filtered governing equations are discretized by a finite-volume scheme. For the numerical procedure, the pressure-implicit with splitting operators (PISO) algorithm is used for p-v corrections; the second-order implicit difference scheme is used for the time-dependent term, the convection and diffusion terms. The stochastic component of the flow field at the inlet boundary is taken into account by superposing random perturbations on each velocity component. The boundary condition at the exit is based on pressure outlet boundary conditions. At the near-wall grid nodes the wall function approximation is employed..

7. Simulation Results and Discussion

The simulated time-averaged results and their assessment by experiments can be found in Ref. [25]. The following are simulated instantaneous results for structures of swirling coal flames, which has not reported before. Figure 2, Figure 3 and Figure 4 show the instantaneous vorticity isolines of swirling coal flames, The large-eddy coherent structures can be seen for all cases. These structures are produced near the inlet and rapidly damped in the downstream region. It is seen that the vortex structures of Cases 2 and 3 are stronger than those of Case 1, indicating that stronger turbulence exist in Cases 2 and 3.
Figure 2. Instantaneous vorticity isolines (1/s)(Case1).
Figure 2. Instantaneous vorticity isolines (1/s)(Case1).
Preprints 86640 g002
Figure 3. Instantaneous vorticity isolines (1/s)(Case 2)
Figure 3. Instantaneous vorticity isolines (1/s)(Case 2)
Preprints 86640 g003
Figure 4. Instantaneous vorticity isolines (1/s)(Case 3)
Figure 4. Instantaneous vorticity isolines (1/s)(Case 3)
Preprints 86640 g004
Figure 5, Figure 6 and Figure 7 give the instantaneous velocity vectors of Cases 1 to 3. It can also be seen that the vortex structures of Cases 2 and 3 are stronger than that of Case 1.
Figure 5. Instantaneous velocity vectors(Case 1)
Figure 5. Instantaneous velocity vectors(Case 1)
Preprints 86640 g005
Figure 6. Instantaneous velocity vectors(Case 2)
Figure 6. Instantaneous velocity vectors(Case 2)
Preprints 86640 g006
Figure 7. Instantaneous velocity vectors(Case 3)
Figure 7. Instantaneous velocity vectors(Case 3)
Preprints 86640 g007
Figure 8, Figure 9 and Figure 10 show instantaneous particle concentration isolines. In Cases 1 to 3, the particles are concentrated in the near-axis zones. In other regions the particle concentration is low. The particle distribution of Case 2 is more uniform than that of Case 1, and that of Case 3 is more uniform than that of Case 2. This is because of different inlet velocities and coal feeding rates and adding inlet purturbation in Case 3. The near-wall particle concentrtion is higher for Case 3 than that of other cases.
Figure 8. Instantaneous particle concentration isolines(Case1)
Figure 8. Instantaneous particle concentration isolines(Case1)
Preprints 86640 g008
Figure 9. Instantaneous particle concentration isolines(Case 2)
Figure 9. Instantaneous particle concentration isolines(Case 2)
Preprints 86640 g009
Figure 10. Instantaneous particle concentration isolines(Case 3)
Figure 10. Instantaneous particle concentration isolines(Case 3)
Preprints 86640 g010
Figure 11, Figure 12 and Figure 13 give the instantaneous temperature isolines. The flames concentrated in the near-inlet highly shear region and then spreaded in the downstream region. The flame length of Case 1 is longer than that of Case 2 and Case 3 due to different inlet velocities and better mixing of coal particles in Cases 2 and 3.
Figure 11. Instantaneous temperature isolines (K)(Case 1)
Figure 11. Instantaneous temperature isolines (K)(Case 1)
Preprints 86640 g011
Figure 12. Instantaneous temperature isolines (K)(Case 2)
Figure 12. Instantaneous temperature isolines (K)(Case 2)
Preprints 86640 g012
Figure 13. Instantaneous temperature isolines (K)(Case 3)
Figure 13. Instantaneous temperature isolines (K)(Case 3)
Preprints 86640 g013
Conclusions
(1) The structures of swirling coal flames were studied in detail by LES
(2) For all cases coherent vortex structures exist in flames, indicated by vorticity isolines and velocity vectors.
(3) For all cases flames are located in the regions of high particle concentration and high Vorticities.
(4) The annular inlet of primary air flow laden with coal particles leads to higher vorticity, stronger mixing and higher temperature near the inlet, hence intensifies combustion than the central inlet.
(5) Larger secondary swirling air flow enhances mixing and intensifies combustion.

Acknowledgement

This study was supported by the National Natural Science Foundation of China under the Grants 50606026 and 50736006.

References

  1. Kurose, R., Makino, H., Large eddy simulation of a solid-fuel jet flame, Combustion and Flame, 2003, 135:1–16.
  2. Zhou, H., Flamant, G., & Gauthier, D. (2004). DEM-LES of coal combustion in a bubbling fluidized bed. Part II: Coal combustion at the particle level. Chemical Engineering Science, 59, 4205–4215.
  3. Yamamoto, K., Murota, T., Okazaki,T., Taniguchi, M., Large eddy simulation of a pulverized coal jet flame ignited by a preheated gas flow, Proceedings of the Combustion Institute, 2011, 33:1771–1778.
  4. Muto, M., Watanabe, H. et al., Large-eddy simulation of pulverized coal jet flame - Effect of oxygen concentration on NOx formation, Fuel, 2015,142:152-163.
  5. Smith Philip, J.,Shen Haoshu et al., Large eddy simulation of a pulverized coal jet flame and its validation, Journal of Engineering Thermophysics, 2019:5.
  6. Wan, K.D., Xia, J. et al., Large-eddy simulation of pilot-assisted pulverized-coal combustion in a weakly turbulent jet, Flow, Turbulence and Combustion, 2017, 99: 531-550.
  7. Wen, X., Jin, H.H. et al., Large eddy simulation of piloted pulverized coal combustion using the velocity-scalar joint filtered density function model, Fuel, 2015, 158: 494-502.
  8. Wen, X., Luo, K. et al., Large eddy simulation of piloted pulverised coal combustion using extended flamelet/progress variable model, Combustion Theory and Modelling, 2017, 5:925-953.
  9. Xing, J.K., Luo, K. et al., Large eddy simulation of Cambridge bluff-body coal (CCB2) flames with a flamelet progress variable model, Proceedings of the Combustion Institute, 2021, 38:5347-5354.
  10. Edge, P., Gubba, S.R., Ma, L., Porter, R. Pourkashanian, M., Williams, A., LES modelling of air and oxy-fuel pulverised coal combustion—impact on flame properties, Proceedings of the Combustion Institute, 2011,33: 2709–2716.
  11. Gharebaghi, M., Irons, R.M.A., Ma, L., Pourkashanian, M., Pranzitelli, A., Large eddy simulation of oxy-coal combustion in an industrial combustion test facility, International Journal of Greenhouse Gas Control, 2011, doi:10.1016.
  12. Rabacal, M , Costa, M , Rieth, M , Kempf, A.M., Particle history from massively parallel large eddy simulations of pulverised coal combustion in a large-scale laboratory furnace, Fuel, 2020, 271, 117587DOI10.1016/j.fuel.2020.117587. [CrossRef]
  13. Watanabe, H., Tanno, K. et al., Large-Eddy Simulation of coal combustion in a pulverized coal combustion furnace with a complex burner, 6th International Symposium on Turbulence, Heat and Mass Transfer, 2009, pp:1027-1030.
  14. Rieth, M., Proch, F. et al., Flamelet LES of a semi-industrial pulverized coal furnace, Combustion and Flame, 2016, 173: 39-56.
  15. Shen, H.S., Wu, Y.X. et al., Large eddy simulation of a 660 MW utility boiler under variable load conditions, Frontiers in Energy, 2021, 15:124-131.
  16. Olenik, G., Stein, O.T., Kronenburg, A., LES of swirl-stabilised pulverised coal combustion in IFRF furnace No. 1, Proceedings of the Combustion Institute, 2015, 35:2819-2828.
  17. Franchetti, B.M., Marincola, F.C. et al., Large eddy simulation of a 100 kW(th) swirling oxy-coal furnace, Fuel, 2016, 181: 491-502.
  18. Smagorinsky, J., General circulation experiments with the primitive equation (I): The basic experiment. Monthly Weather Rev, 1963, 91 (3): 99-164.
  19. Germano, M., Piomelli, U., Moin, P., et al., A dynamic sub-grid-scale eddy viscosity model, Physics of Fluids, 1991, A 3:1760-1765.
  20. Kim, W. W., Menon, S. S., A new dynamic one-equation sub-grid-scale model for large eddy simulation, AIAA 95-0356, 1995, 33th Aerospcace Sciences Meeting and Exhibition, USA.
  21. Lilly, D.K., On the application of the eddy viscosity concept in the inertial subrange of turbulence, NCAR Manuscript, 123, 1966.
  22. Magnussen, B. F., Hjertager, B. H., On mathematical modeling of turbulent combustion with special emphasis on soot formation and combustion, In: Combustion Inst. Pittsburgh PA, eds. 16th Symposium (International) on Combustion. Pittsburgh: The Combustion Institute, 1976, pp. 719-729.
  23. Zhou, L.X., Theory and Numerical Modeling of Turbulent Gas-Particle Flows and Combustion, CRC Press, 1993.
  24. Hu, L.Y., Zhou, L.X. et al., Measurements of swirling coal combustion and NO formation, 7th International Symposium on Measurement Techniques of Multiphase Flows, 2011.
  25. Hu, L.Y., Zhou, L.X. et al., Measurement and simulation of swirling coal combustion, Particuology, 2013, 11: 189-197, http://dx.doi.org/10.1016/j.partic.2012.05.009. [CrossRef]
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

© 2024 MDPI (Basel, Switzerland) unless otherwise stated