Unsteady 3-D RANS simulations of dust explosion in a fan stirred explosion vessel using an open source code

Dust explosion is a constant threat to industries which deal with combustible powders such as woodworking, metal processing, food and feed, pharmaceuticals and additive industries. The current standards regarding dust explosion venting protecting systems, such as EN 14491 (2012) and NFPA 68 (2018), are based on empirical correlations and neglect effects due to complex geometry. Such a simplification may lead to failure in estimating explosion overpressure, thus, increasing risk for injuries and even fatalities at workplaces. Therefore, there is a strong need for a numerical tool for designing explosion protecting systems. This work aims at contributing to the development of such a tool by (i) implementing a premixed turbulent combustion model into OpenFOAM, (ii) verifying the implementation using benchmark analytical solutions, and (iii) validating the numerical platform against experimental data on cornflour dust explosion in a fan-stirred explosion vessel, obtained by Bradley et al. (1989a) under well-controlled laboratory conditions. For this purpose, the so-called Flame Speed Closure model of the influence of turbulence on premixed combustion is adapted and implemented into OpenFOAM. The implementation of the model is verified using exact and approximate analytical solutions for statistically one-dimensional planar and spherical turbulent flames, respectively. The developed numerical platform is applied to unsteady three-dimensional Reynolds Averaged Navier-Stokes simulations of the aforementioned experiments. The results show that the major trends, i.e. (i) a linear increase in an apparent turbulent flame speed St;b with an increase in the root mean square (rms) turbulent velocity u’ and (ii) and an increase in St;b with an increase in the mean flame radius, are qualitatively predicted. Furthermore, the measured and computed dependencies of St;bðu’Þ agree quantitatively under conditions of weak and moderate turbulence.


Introduction
Dust explosion is a constant threat to industries which deal with combustible solids, e.g. woodworking, metal processing, food and feed, pharmaceuticals and additive industries. A dust explosion can occur if the concentration of dust-air cloud is within flammability limits and there is an energy supply sufficient for igniting the mixture. Statistics shows that there is one serious dust explosion every day in Europe alone (Beck and Jeske, 1982).
To manage the dust explosion risk, a proper risk assessment must be carried out. As a basis of the risk assessment, one have to establish the explosion characteristics of the dust, which involve Minimum Explosive Concentration (MEC), deflagration index K St , maximum explosion overpressure P max , Minimum Ignition Energy (MIE), Minimum Ignition Temperature (MIT), etc. These data are then used to assess the likelihood and the consequences of an explosion. Subsequently, risk evaluation is performed by comparing results of the risk assessment with established risk criteria in order to determine whether additional protection measures are required. Such measures can be of both organisational and technical nature. The organisational measures involve training the staff to raise awareness of the explosion hazards, proper housekeeping, maintenance, procedures for hot work permit, etc. The technical measures include, e.g., installation of dust extraction systems to limit explosive atmospheres, the use of a proper ATEX-approved equipment in explosive atmospheres to avoid ignition, or installation of explosion protection systems such as venting, suppression, isolation, and containment.
The consequences of a dust explosion are highly dependent on the explosion development, i.e., the rate of the explosion pressure rise, the maximum explosion pressure, the flame propagation speed, etc. Accordingly, the explosion development depends on dust characteristics, operating conditions (pressure, flow rate, turbulence level, initial temperature), geometry of the equipment (scale, interconnection of vessels, congestion, etc.), as well as the location of ignition point and the strength of the ignition source. Consequently, assessing the explosion severity is a huge challenge, especially in facilities with complex processes and geometries. Therefore, there is a great need for new tools to improve predicting explosion patterns. Development of Computational Fluid Dynamics (CFD) tools is a promising approach to meet this request.
The advantages of CFD tools are as follows. First, they can be used to design explosion protecting systems for process plants with complicated geometries where the current standards, e.g. EN14491 (2012), VDI3673 (2002) and NFPA68 (2018), are not applicable. Second, many virtual experiments can be performed using CFD tools, thus, saving a lot of money for avoiding expensive large-scale experiments. Third, physics-based dust explosion models and numerical tools can be used to improve the current standards regarding explosion protection design.
Open Field Operation and Manipulation (OpenFOAM, 2019) software is a suitable platform for developing efficient CFD tools for research into dust explosions. It is a free, open-source general-purpose CFD software package mainly for simulating thermodynamics, fluid dynamics, and chemical reactions. On the technical side, OpenFOAM excels in modern architecture using object-orientated programming language, high parallelization and unstructured grid for dealing with curved geometry. New models and methods can be implemented and tested thanks to the open source. Furthermore, OpenFOAM creates more value for a customer, because it offers an opportunity to create a tailor-made tool that suits the special need of the customer at zero license cost. At the same time, the cost of personal hours in using Open-FOAM may be higher when compared to a commercial CFD code, which benefits from detailed documentation, training, and dedicated technical support. It is also worth noting that OpenFOAM is distributed under version 3 of the GNU General Public License (GPLv3), which gives the users a great freedom in adapting it. For example, users are free (i) to use OpenFOAM for both commercial and non-commercial purposes, (ii) to change the code, (iii) to share the code, and (iv) to share the changes they made. However, this freedom may have negative consequences, e. g., the code developments often suffer from the lack of the detailed documentation readily available to the community.
OpenFOAM was already applied to simulating gas explosions. For instance, Bauwens et al. (2008) studied methane-air deflagration in a vented enclosure experimentally and numerically within the framework of Large Eddy Simulation (LES). Encouraging agreement in terms of explosion overpressure and flame velocities were obtained between experiments and simulations. Later, Bauwens et al. (2011) modelled vented explosion of lean hydrogen-air mixtures by taking into consideration the Darrieus-Landau and Rayleigh-Taylor instabilities.
Based on OpenFOAM, Vendra and Wen (2019) developed a CFD tool for LES of vented deflagration of lean hydrogen-air mixtures. The combustion model used by them is based on (i) a default premixed turbulent combustion model in OpenFOAM (Weller et al., 1998) and (ii) a model of the Darrieus-Landau and Rayleigh-Taylor instabilities, proposed by Bauwens et al. (2011). Sinha et al. (2019) also used the above-mentioned tool for the same application.
Examples of comprehensive benchmark studies on vented hydrogen explosions were reported by Tolias et al. (2018), Vyazmina et al. (2019), and Skjold et al. (2019aSkjold et al. ( , 2019b. The scatter of their results, e.g. explosion overpressure yielded by different CFD tools, calls for more research in this area. As far as dust explosions are concerned, Spijker et al. (2013) modelled lycopodium dust explosion in a tube using a modified solver in OpenFOAM. More often, dust explosions are computed using other numerical tools.
In particular, within the framework of a European Union project called DESC (Dust Explosion Simulation Code), researchers developed the DESC code, which later became a submodule in a commercial CFD software FLACS distributed by Gexcon for simulating gas and dust explosions (Skjold, 2007). The dust explosion module, i.e. FLACS-DustEx, has been developed by comparing simulation results with data measured in various experiments of different scales (Skjold, 2003(Skjold, , 2007Skjold et al., 2005Skjold et al., , 2006. The FLACS-DustEx flame propagation model addresses fine dust particles with high volatile contents (Skjold, 2014) and, accordingly, is based on a model used earlier for gas-air explosions. Tasc� on and Aguado (2015) applied FLACS-DustEx to study dust explosions in industrial scenarios and reported satisfactory agreement between computed results and experimental data or dust explosion venting standards. However, certain problems need to be resolved to improve the model and code, e.g. (i) the use of a single-block structured mesh, which yields over-prediction of momentum dissipation over sloping area (Gant and Hoyes, 2010), (ii) the lack of choice of different turbulence models when compared to a general-purpose CFD software, and (iii) grid-dependence problem associated with combustion model (Skjold, 2014;Tasc� on and Aguado, 2017).
Besides commercial software, inhouse code were also developed by several research groups for studying gas and dust explosions. For example, Ugarte et al. (2016) and Sezer et al. (2017) simulated vented gas explosions using a phenomenological model and a single-zone approximation. Demir et al. (2017Demir et al. ( , 2018 modelled gas and dust explosions in coal mines by studying two flame acceleration mechanisms and compressibility effects. Moreover, numerical studies of induced layered dust explosions were performed by Song et al. (2017), Song and Zhang (2019), and Shimura and Matsuo (2019), but codes used in those studies are not discussed in the cited papers.
Furthermore, hydrogen explosion with suspended inert dust particles were studied by Liberman et al. (2015) in the framework of Direct Numerical Simulation (DNS). However, DNS technique is not realistic for solving engineering problems due to its extremely high computational cost.
Thus, gas and dust explosions are computed adopting various numerical tools that range from (i) engineering models (e.g. FLACS-DustEx), which can be used at industrial scales, but invoke oversimplified phenomenological models of ignition and turbulent combustion, to (ii) detailed models that allow for basic physical phenomena at relatively small spatial scales, but are not yet feasible for simulations at industrial scales. Accordingly, a comparison between the different types of models is not straightforward. Nevertheless, due to rapid development of computer hardware, there is clear trend towards applications of detailed models to problems characterized by larger scales. In line with this trend, there is a growing interest in applying efficient and well-validated detailed models of turbulent combustion to safety problems.
In particular, under certain conditions, the dust explosion process resembles turbulent burning of a gas cloud. This is especially true for very fine organic dust particles with high volatile content (Bradley et al., 1988(Bradley et al., , 1989a1989b). Therefore, a premixed turbulent combustion model is sometimes used for CFD modelling of dust explosions (Spijker et al., 2013;Skjold, 2014). The present work aims at adapting the so-called Flame Speed Closure (FSC) model of the influence of turbulence on premixed burning Chomiak, 1997, 2002a) for this purpose. The main reason is that, as reviewed elsewhere (Lipatnikov and Chomiak, 2002a;Lipatnikov, 2012Lipatnikov, , 2018, the FSC model has been quantitatively validated against a wide set of experimental data obtained by various research groups from various (both expanding and statistically stationary) flames under a wide range of substantially different conditions (various fuels, equivalence ratios, initial temperatures, pressures, rms turbulent velocities, and turbulent length scales). Moreover, the FSC model is an extension of the well-known Turbulent Flame Closure (TFC) model Lipatnikov, 1993, 1995;Karpov et al., 1996), which has been implemented into various commercial CFD codes, e.g., Ansys CFX (Ansys, 2020a), Ansys Fluent (Ansys, 2020b), Converge (2020), FINE (FINE, 2019), AVL FIRE (AVL FIRE, 2020). Accordingly, the TFC model is widely used by automotive and gas turbine industry.
More specifically, this work aims at (i) implementing the FSC model into the open source toolbox OpenFOAM, (ii) verifying the model implementation using exact and approximate analytical solutions for statistically 1-D, planar and spherical, respectively, premixed flames expanding in "frozen" turbulence, and (iii) validating the developed numerical platform against experimental data obtained by Bradley et al. (1989a) utilizing the well-known Leeds fan stirred explosion vessel.
In the next section, the FSC model is briefly summarized. Numerical and experimental setups adopted for verification of the model implementation and for validation of the numerical platform are described in Sect. 3. Computed results are discussed in Sect. 4, followed by conclusions.

Model equations
The FSC model (i) characterizes the thermochemical state of a reacting mixture in a flame using a single combustion progress variable c, which is equal to zero and unity in fresh reactants and equilibrium combustion products, respectively, (ii) invokes the following wellknown Bray-Moss-Libby (BML) equations (Bray and Moss, 1977;Libby and Bray, 1977) ( 1) and (iii) deals with the following transport equation for the Favre-averaged combustion progress variable c. Here, ρ is the density; σ ¼ ρ u =ρ b is the density ratio; t is the time; u is the flow velocity vector; κ is the molecular heat diffusivity of the mixture; D t and U t are the turbulent diffusivity and burning velocity, respectively; Q is a source term discussed later, see Eq. (6); over-lines designate the Reynolds average, while q ¼ ρq=ρ is the Favre-averaged value of q with q} ¼ q À q; subscripts u and b designate unburned and burned gas, respectively. Within the framework of the FSC model, D t and U t are evaluated as follows Chomiak, 1997, 2002a) where D t;∞ is the fully developed turbulent diffusivity, which can be determined using a turbulence model, as will be discussed later; t fd is the flame development time counted starting from end of ignition; τ L ¼ D t;∞ =u' 2 is the Lagrangian time scale of turbulence; u' is the rms turbulent velocity; is an intermediately steady turbulent burning velocity; A ¼ 0.4 (Lipatnikov and Chomiak, 1997) is the sole constant of the FSC model; Da ¼ τ t =τ f is the Damk€ ohler number; τ t ¼ L=u' and τ f ¼ δ L =S L are turbulent and laminar-flame time scales, respectively; L is an integral turbulent length scale; S L and δ L ¼ κ u =S L are the laminar flame speed and thickness, respectively. As discussed in detail elsewhere (Zimont, 1979;Lipatnikov and Chomiak, 2002a;Lipatnikov, 2012), at moderate turbulence, Eq. (5) is qualitatively consistent with various experimental data on the influence of mixture composition, turbulence characteristics, and pressure on turbulent burning velocity or flame speed. Originally, Eq. (5) was analytically derived by Zimont (1979) by assuming that (i) small-scale eddies increase local burning velocity by thickening flamelets and increasing heat and mass transfer within them, with the width of the thickened flamelets being still significantly smaller than L; (ii) large-scale eddies increase local burning velocity by wrinkling the thickened flamelets; and (iii) the mean turbulent flame brush thickness δ t grows by the turbulent diffusion law. Such a regime characterized by apparently stationary turbulent burning velocity given by Eq. (5), but growing δ t was later called Intermediate Steady Propagation (ISP) regime (Zimont, 2000). Numerous experimental data reviewed elsewhere (Prudnikov, 1967;Lipatnikov and Chomiak, 2002a) indicate that such a combustion regime is a widespread regime of premixed turbulent burning.
The original derivation of Eq. (5) was performed under the following constraints: the turbulent Reynolds number (Zimont, 1979). Here, ν u is the kinematic viscosity of unburned gas. Subsequently, Lipatnikov and Chomiak (2002a) argued that the aforementioned assumption (i), i.e. thickening of flamelets by small-scale eddies, could be substituted with a more general assumption that the interaction between small-scale turbulent eddies and flamelets is controlled by the mean dissipation rate ε and chemical time scale τ c . Under this assumption, which is in fact an extension of the well-recognized Kolmogorov hypothesis to the case of premixed turbulent combustion, the constraint of Ka > 1 is substituted with u'=S L > 1 and the model is applicable to moderately turbulent burning also.
(2), the FSC model reduces to the TFC model by Lipatnikov (1993, 1995). The time-dependent terms in square brackets in Eqs. (3) and (4) extend the TFC model and allow us to simulate early stages of premixed turbulent flame development, including the formation of a small flame kernel after ignition, transition to turbulent burning, and development of the turbulent flame. Equation (3) is well known in the turbulence literature (Hinze, 1975) and results from the Taylor (1935) theory of turbulent diffusion. Equation (4) was derived by Lipatnikov and Chomiak (1997) by adapting the Taylor's theory to extend the Zimont model of the intermediately steady turbulent burning velocity.
In order to (i) simulate an early stage of flame kernel growth after spark ignition and (ii) obtain an appropriate balance equation in the limit case of u ' →0, the TFC model was further extended and the following source term Chomiak, 1997, 2002a was incorporated into in Eq.
(2). Here, Θ is the activation temperature for a single reaction that the combustion chemistry is reduced to (Θ ¼ 20000 K in the present work); the Favre-averaged temperature T is evaluated using the simplest form ρT ¼ ρ u T u of the ideal gas state equation; and the reaction time scale t r is set so that, in the case of u ' ¼ 0, the burning velocity yielded by stationary, 1-D Eqs.
(1), (2) and (6) is equal to the laminar burning velocity S L , which is an input parameter of the model. This constraint results in where the non-dimensional function Ψ approximates values of S L ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi t r =κ u p , pre-computed for various ratios of T b =T u and Θ=T u by numerically integrating 1-D Eqs. (1), (2) and (6) with D t ¼ U t ¼ 0. In this case, the FSC Eq. (2) reduces to an equation that models a laminar premixed flame in the case of a single combustion reaction, with the source term Q being introduced into Eq. (2) by Chomiak (1997, 2002a) in order to satisfy this constraint. A polynomial approximation of the function Ψ does not feature any tuning parameter and is reported by Huang et al. (2016).
(1)-(4) for a statistically 1-D planar flame that propagates from left to right.
Equations (8) and (9) describe a flame with a self-similar mean structure. As reviewed elsewhere (Prudnikov, 1967;Lipatnikov and Chomiak, 2002a;Driscoll, 2008;Lipatnikov, 2018) and supported by more recent experimental data (Tamadonfar and Gülder, 2014;Han et al., 2018), various premixed turbulent flames do have such a self-similar mean structure well described by Eqs. (8) and (9). It is worth stressing that a transport equation, which (i) was basically similar to Eq.
(2), (ii) had the exact solution given by Eqs. (8)-(11), but (iii) was written in a different form, was introduced into the combustion literature by Prudnikov (1967) by considering statistically 1-D planar case.
In the statistically 1-D, but spherical case, the solution given by Eqs. (8)-(11) is not exact. Nevertheless, if the mean flame structure is assumed to be self-similar, i.e. cðr; tÞ ¼ c½ðr À R f Þ =δ t � where r is the radial distance and R f is a mean flame radius, the following analytical relation holds (Lipatnikov and Chomiak, 2007) dR f dt for a particular mean flame radius defined elsewhere (Lipatnikov and Chomiak, 2002b), with this result agreeing quantitatively with the Leeds experimental data (Bradley et al., 2003).
Finally, as already noted in Sect. 1 and discussed in detail elsewhere (Lipatnikov and Chomiak, 2002a;Lipatnikov, 2012Lipatnikov, , 2018, the FSC model was quantitatively validated in RANS simulations of a wide set of experiments performed by various research groups with various (both expanding and statistically stationary) flames under a wide range of substantially different conditions (various fuels, equivalence ratios, initial temperatures, pressures, rms turbulent velocities, and turbulent length scales).

Numerical and experimental setups
Equations (1)-(7) were implemented into OpenFOAM, followed by application of the developed numerical platform to simulating experiments by Bradley et al. (1989a). In this section the numerical setups for verification of the model implementation is reported, followed by descriptions of the experimental and numerical setups for the Leeds fan-stirred explosion vessel.

Numerical setup for verification of model implementation
In order to directly verify the implementation of the FSC model into OpenFOAM, statistically 1-D, planar or spherical, premixed turbulent flames propagating in a statistically "frozen" turbulence characterized by stationary and spatially uniform u', L, and D t;∞ were simulated and computed results were compared with the exact or approximate, respectively, solution given by Eqs. (8)-(11) or (12), respectively. The simplification of statistically "frozen" turbulence was necessary at this stage of research, because the solutions were obtained in the case of spatially uniform u', L, and D t;∞ . While the two benchmark problems are statistically 1-D, the RANS simulations were performed in 3-D cases.
Bearing in mind subsequent application of the FSC model to the dust explosion experiments, the simulation input parameters were set based on the experimental conditions. In particular, premixed burning of cornflour dust cloud was simulated. The cornflour chemical equivalent formula is C 6 H 7.88 O 4.98 with a heat of reaction being 15.8 MJ/kg (Bradley et al., 1989a). The corresponding single-step chemical mechanism is as follows To verify the model implementation in the planar benchmark case, the computational domain is set to be a parallelepiped of a length of 0.1 m and a cross section of 0.003 � 0.003 m. The numerical mesh consists of 100 cells in the x direction and three cells in the y or z direction. The flame propagates from left to right. Zero velocity and free entrainment boundary conditions are set on the right (unburned) and the left (burned) boundaries, respectively.
To verify the model implementation in the spherical benchmark case, the computational domain is set to be a cube, which represents one eighth of the total volume filled with reactants or products. The cube side length is 60 mm and the mesh size is 0.25 mm in each direction. Totally 13 824 000 cells are created. The initial conditions correspond to a spherical kernel of a radius of 20 mm, filled with combustion products (c ¼ 1). The rest of the domain is filled with unburned mixture (c ¼ 0). Symmetry boundary conditions are set at three boundaries. Zero velocity and free entrainment boundary conditions are set at other three boundaries.
Other details are provided in Tables 1-3, where W is molecular weight, μ is the molecular dynamic viscosity of a mixture, k is the Favreaveraged turbulent energy, and b ¼ 1 Àc is the combustion regress

Experimental setup for leeds fan stirred explosion vessel
The Leeds experiments (Bradley et al., 1989a) were performed by (i) filling a vessel with air and cornflour dust, (ii) mixing the dust and air by counter-rotating fans, which also generated statistically stationary homogeneous isotropic turbulence in the central part of the vessel, (iii) igniting the mixture by a spark located in the vessel centre, and (iv) recording images of expanding flames using high speed Schlieren technique. The vessel diameter is equal to 305 mm and its volume is about 0.023 m 3 . The vessel has three pairs of orthogonal quartz windows of 150 mm diameter for optical measurements.
Turbulence was generated by four fans, whose rotation speed was changed to vary the rms turbulent velocity u'. In the discussed experiments, the fan speed was varied from 8 to 50 Hz, which corresponded to variations in u' from 0.80 to 5.0 m/s. While the integral length scale was not reported by Bradley et al. (1989a), it was reported in other papers by the Leeds group. In particular, Bradley et al. (2003) stated that the longitudinal integral length scale measured using laser Doppler velocimetry was found "to be 20 mm and independent of fan speed between 1000 and 10 000 rpm", with 1000 rpm corresponding to 16.5 Hz. It is worth noting, however, that, in the experiments with the lean dust-air mixture, the lowest fan speed was less than 16.5 Hz and a decrease in L at low fan speeds was reported in an earlier paper by the Leeds group (Abdel-Gayed et al., 1984). However, those data cannot be used here, because they were obtained using thermo-anemometry, but such a method performs poorly in flows with zero mean velocity. For instance, the earlier Leeds measurements with thermo-anemometry overestimated L at large fan speeds by a factor of about two. Thus, in the Leeds experiments with the dust-air mixture, the turbulence length scale of 20 mm could be overestimated at low fan speeds. Nevertheless, when compared to other experimental data on dust explosions, the Leeds measurements were performed under well-defined laboratory conditions, i.e. the initial and boundary conditions were controlled sufficiently well.
To study dust explosion, a premixed dust-air cloud was ignited by a spark in turbulent medium in the centre of the vessel. Subsequently, turbulent flame kernel growth was recorded using high-speed Schlieren system. By processing Schlieren images, an equivalent mean flame radius R f , i.e. the radius of a circle whose area was equal to the area enveloped by the flame surface on the image, was calculated and turbulent flame speed with respect to combustion products was evaluated by differentiating the measured R f ðtÞ-curves, i.e.
To mitigate an influence of the spark on the flame speed, the measurements were performed in a range of 20 mm� R f ðtÞ �35 mm. For such flame kernels, whose radius was less than the vessel radius by a factor of more than four, an increase in the pressure in the vessel was negligible.
In addition to the values of S t;b obtained at four different R f and five different fan speeds, Bradley et al. (1989a) also reported the values of the laminar flame speed S L and density ratio σ for the studied dust-air mixture. However, methods and precision of evaluation of these quantities are not discussed in the cited paper. Furthermore, the value of the laminar flame thickness δ L , which is required to calculate an important input parameter of the FSC mode such as the chemical time scale τ c ¼ δ L =S L , is not reported either. Thus, even in the considered case of the small-scale well-controlled Leeds experiments, some information important for the model validation is missing. This is a typical problem for testing any model of dust explosion.
In this regard, it is worth noting that one of the biggest challenges associated with dust explosion modelling consists of a limited amount of data on the laminar flame speed. Unlike gas explosions, where S L is a well-studied basic characteristic, the laminar flame speed of a dust-air cloud is rarely available. This parameter depends not only on the dust concentration, but also the dust particle size distribution and moisture content. Indeed, Cloney et al. (2018) studied the laminar flame propagation in a dust-air cloud using a detailed CFD model in OpenFOAM. The model took into account surface reactions, devolitization, and gas phase reactions. The reported results indicate substantial effects of the particle size distribution and initial temperature on the computed S L . A method for modelling the effect of the particle size distribution on the speed of a laminar dust-air flame was recently proposed by Ghaffari et al. (2019). Fortunately, Bradley et al. (1989a) reported the value of S L for the conditions of their experiments and this value is used in the present study.

Simulations of dust explosion in the Leeds fan stirred explosion vessel
To save computational time, one eighth of a cube whose volume is equal to the volume of the vessel is simulated, as sketched in Fig. 1. Accordingly, the length of the computational domain is 0.14 m. A grading method is used to generate the computational mesh with a finer mesh of about 0.125 mm in the centre and coarser mesh of about 3.6 mm in the far-field. The total mesh consists of 2 197 000 cells. A single simulation run takes around 2-4 days on two cores on a computer with 128 GB RAM and totally 28 Xeon Gold cores for a simulation duration of   Fig. 1. Schematic illustration of cornflour dust explosion model. 30 ms using a maximum Courant number of 0.1. Grid sensitivity of numerical results was tested by varying the mesh in the centre from 0.0625 to 0.25 mm, with the computed dependencies of S t;b on R f being close to one another. The thermo-physical properties and boundary conditions are reported in Tables 1 and 3, respectively. It is worth noting that the measured burned temperature of 1500 K (Bradley et al., 1989a) is used here instead of the calculated burned temperature of 1592 K, because neither the method, nor precision of the calculation is discussed in the cited paper. The use of the former temperature yields the density ratio of 5.06, whereas σ ¼ 5.49 reported by Bradley et al. (1989a) corresponds to the latter (higher) temperature.
To mimic spark ignition, the following extra source term is added on the right hand side of Eq.
(2). Here, Hðt 0 À tÞ is Heaviside function, while W 0 , σ r , σ t are parameters of the ignition submodel. More specifically, t 0 is associated with ignition time, σ t characterizes ignition duration, and σ r corresponds to the size of ignition kernel. The factor W 0 is associated with the ignition strength and should be set sufficiently large in order for cðr ¼ 0; t 0 Þ to be close to unity. A similar method was earlier used by Lipatnikov (1993, 1995) and Karpov et al. (1996) to simulate experiments with flames expanding in a fan-stirred bomb, but the Heaviside function was skipped in the cited papers. The point is that if Q ¼ 0 in Eq.
(2), then, any cðx;tÞ ¼ const is a solution to it. Accordingly, due to numerical errors, cðx; tÞ tends often to zero in the entire computational domain. The use of the algebraic source term Q modelled with Eq. (6) resolves the problem (Lipatnikov and Chomiak, 1997), but provided that the magnitude of Q is sufficiently large. Under conditions of the present simulations, this is not so due to a low value of the laminar flame speed for the considered dust-air mixture. Therefore, to retain cðr ¼ 0; t > t 0 Þ � 1, the source term given by Eq. (14) is independent of time if t > t 0 . Nevertheless, due to a low value of σ r , the ignition source term is significant at small r only and weakly affects the speed S t;b of large flame kernels. For instance, S t;b computed at 20 mm� R f ðtÞ �35 mm setting either W 0 ¼ 10 7 s À 1 and σ r ¼ 0:5 mm or W 0 ¼ 10 13 s À 1 and σ r ¼ 0:25 mm are very close to one another. Results reported in Sect. 4.2 were obtained using σ r ¼ 0:25 mm, W 0 ¼ 10 13 s À 1 , and t 0 ¼ 1 ms.
To evaluate turbulence characteristics required by the FSC model, the following equations and the k À ε turbulence model (Launder and Spalding, 1972) implemented into the standard OpenFOAM are used. Here, Pr t is the turbulent Prandtl number, C μ ¼ 0:09 is a constant of the k À ε model, and C d ¼ 1:0. The turbulence model is switched on at t � 2t 0 to avoid unphysically strong generation of turbulence due to rapid density drop in the vessel centre, caused by the ignition source term W ign . To mimic the flux of turbulent energy from the fans to the vessel centre and to simulate statistically stationary turbulence in line with the experiments, the extra source terms ρε 0 and ρε 2 0 =k 0 are added on the right hand sides of the transport equations for k and ε, respectively, following Zimont and Lipatnikov (1995). Here, k 0 and ε 0 are the initial values of the turbulent kinetic energy and its dissipation rate, respectively, reported in Table 4. Due to this modification, the values of u' and L upstream of the flame remain close to their initial values and are weakly affected by C d .

Results and discussions
In this section, results of verification of the model implementation will be presented, followed by validation of the FSC model against the experimental data by Bradley et al. (1989a). Fig. 2a shows that the spatial profiles of the Reynolds-averaged combustion progress variable c, simulated for the statistically 1-D planar flame, change with time. However, the same profiles plotted vs. the normalized distance ξ defined by Eq. (9) are almost the same at all time instants, see broken lines in Fig. 2b, in line with the welldocumented self-similarity of premixed flames (Lipatnikov, 2009). Moreover, this computed self-similar profile agrees very well with the analytical solution given by Eq. (8), cf. solid and broken lines in Fig. 2b.

Verification of model implementation
Comparison between calculated and analytical flame speeds is shown in Fig. 3. In the simulations, the flame speed is evaluated by taking derivative of the mean flame position against time. The mean flame position is defined by the x-coordinate of a surface of cðx; tÞ ¼ 0:5 Chomiak, 2002b, 2007). Since the calculated flame speed exhibits fluctuations, UnivariateSpline function in scipy library of Python is used to smooth the data. Fig. 3 shows that the numerical and analytical results agree well. There is a slight discrepancy in the beginning of the simulations. This is caused by the usage of smoothing function in python. Fig. 4 further verifies the implementation of the FSC model by showing that the mean flame brush thickness evaluated by processing the computed profiles of cðx; tÞ using the following equation agrees very well with the analytical solution given by Eq. (11). Fig. 5 shows that the spatial profiles of the Reynolds-averaged combustion progress variable c, simulated for the statistically 1-D spherical flame, change with time but collapse to the same curve when plotted vs. the normalized distance ξ defined by Eq. (9), see broken lines in Fig. 5b, in line with the well-documented self-similarity of premixed flames. Moreover, this computed self-similar profile agrees well with the analytical solution given by Eq. (8), cf. solid and broken lines in Fig. 5b. It is worth remembering that both the self-similarity of the mean flame structure and Eq. (8) are well supported by experimental data obtained from various flame configurations and different (statistically stationary or expanding) types of flames (Lipatnikov and Chomiak, 2002a;Driscoll, 2008;Lipatnikov, 2012).

Model validation
The discussed combination of the FSC model of the influence of turbulence on combustion and the k À ε model of turbulence involves  (15). In the present study, the standard value of A ¼ 0:4 was used, but Pr t was tuned to get the best agreement with the experimental data. Sensitivity of the computed results to Pr t is illustrated in Fig. 7. In the literature, different values of Pr t can be found, with Pr t being varied between 0.3 and 1.0 when using Eq. (15) jointly with the TFC or FSC model (Lipatnikov, 2018). In the present study, the best agreement with the measured data is obtained for Pr t ¼ 0:4, which is sufficiently close to Pr t ¼ 0:3 recommended in recent papers (Yasari et al., 2015;Verma and Lipatnikov, 2016). Results of the model validation are reported in Figs. 8 and 9. Here, the mean flame position is defined by the x-coordinate of a surface of cðx; tÞ ¼ 0:1, which differs from a value used to verify the model implementation. The reason is that cðx; tÞ ¼ 0:1 is associated with the high speed Schlieren technique used in the experiments, because this technique captures the leading edge of the mean flame brush, whereas Eq.
(12) has been obtained for the radius of a surface of cðx; tÞ ¼ 0:5 (Lipatnikov and Chomiak, 2007). Comparison of open and filled symbols in Fig. 8 shows that the simulations well predict an increase in S t;b by u ' in weak and moderate turbulence. Note that S t;b measured by Bradley et al. (1989a) at u ' ¼ 5 m/s is significantly less than at u ' ¼ 3:31 m/s. Such a decrease in flame speed by the rms turbulent velocity in intense turbulence is well documented in various experiments with flames expanding in a fan-stirred bomb since the seminal study by Karpov et al. (1959). However, available combustion models cannot predict this phenomenon without tuning, as discussed in detail elsewhere Chomiak, 2002a, 2005c). For this reason, the present study is restricted to u ' � 3:31 m/s. It is of interest to note that Fig. 8 shows almost linear correlation between computed (or measured) flame speeds and rms turbulent velocity. For instance, for the mean flame radius of 30 mm, the computed results (open squares) are well fitted with the following linear equation On the contrary, Eq. (5) adapted by the FSC model yields a weaker dependence of intermediately steady turbulent burning velocity on the rms turbulent velocity, i.e. U t;ISP ∝u' 0:75 . The point is that, as shown by Verma and Lipatnikov (2016), transient, curvature and strain effects can significantly change the scaling exponent for simulated turbulent flame speed as a function of the rms turbulent velocity. Fig. 9 shows that the simulations (see lines) reasonably well predict a slow increase in the mean flame speed with the mean flame radius at various u ' , while the effect magnitude is underestimated at u ' ¼ 3:31 m/ s, especially at smaller R f .

Concluding remarks
The Flame Speed Closure (FSC) model of the influence of turbulence on premixed combustion was implemented into the open source platform OpenFOAM. The implementation of the FSC model was numerically verified using benchmark analytical solutions in the cases of statistically 1-D planar and spherical turbulent flames. The FSC model supplemented with the well-known k À ε model of turbulence was applied to unsteady 3-D RANS simulations of small-scale laboratory dust-explosion experiments performed under well-controlled conditions using the Leeds fan-stirred combustion vessel. Numerical results show that the model well predicts (i) an increase in the apparent turbulent   flame speed by the rms turbulent velocity at moderate turbulence and (ii) a slow increase in the flame speed with growth of the mean flame radius.
Thus, the present work indicates that the FSC model could be an appropriate building block for developing an advanced numerical tool for CFD research into large-scale explosions of fine dust particles with high volatile contents. For this purpose, models of other phenomena (e. g. Darrieus-Landau and Rayleigh-Taylor instabilities of initially laminar flame kernel, acceleration of large-scale flames attributed commonly to flame-generated turbulence, radiative heat transfer, etc.) should be invoked, combined with the FSC model and implemented into Open-FOAM. This is the goal for our future work.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.