Molybdenum Disulphide Precipitation in Jet Reactors: Introduction of Kinetics Model for Computational Fluid Dynamics Calculations

In our previous work, we used the population balance method to develop a molybdenum disulphide kinetics model consisting of a set of differential equations and constants formulated to express the kinetics of complex chemical reactions leading to molybdenum disulphide precipitation. The purpose of the study is to improved the model to describe the occurring phenomena more thoroughly and have introduced computational fluid dynamics (CFD) modelling to conduct calculations for various reactor geometries. CFD simulations supplemented with our nucleation and growth kinetics model can predict the impact of mixing conditions on particle size with good accuracy. This introduces another engineering tool for designing efficient chemical reactors.


Introduction
The main goals of chemical and process engineering are to select the appropriate apparatus and process operating parameters. The designer can efficiently combine time and space scale analyses to determine the proper configuration. Time scale analysis determines how accurate mathematical models one should use for the process description, thereby allowing the designer to choose from classical chemical reactor engineering methods, all the way to the most advanced direct models for all the process stages and then to select a mathematical modelling method based on mass, momentum, energy, and population balances. Numerous works on time scale analysis [1,2] have shown that the course of many processes of practical importance can be influenced by choosing the reactant contact method, the mixing intensity (including the choice of local flow and mixing parameters such as the local shear rate, stresses, and dissipation rate), and the process running method (e.g., continuous, semiperiodic, and periodic). Proper selection of process methods and parameters often results in better product yield and quality and fewer byproducts (including environmentally harmful ones) and often reduces the number of production stages required. Recently, computational models in combination with computational fluid dynamics (CFD) have been increasingly used for this purpose [3][4][5][6][7][8].
CFD was pioneered by Johnson and Richardson [9], who developed an iterative numerical method of calculating Laplace equations to determine variables based on the numbers calculated in the previous model iteration. Because all the numbers were calculated manually, the method applicability was originally considerably limited. The introduction of computers in the early 1950s led to a major development in computational fluid mechanics methods. The first works were related to the aviation, automotive, and armaments industries [10]; it was not until the 1980s that CFD was first applied to chemical engineering [11].
Since then, many scientific papers have used this method to control the final characteristics of chemical products by selecting the process conditions, particularly the reactant mixing conditions [3,[12][13][14][15].
Molybdenum disulphide (MoS 2 ) is a transition-metal dichalcogenide (TMD) and a valuable 2D-nanomaterial used in a wide range of industrial applications. The main reason for such a wide range of applications of molybdenum disulphide in various fields is its 2D terrace structure and its crystalline properties. One promising application is the use of MoS 2 in several catalytic reactions such as hydrogen evolution reactions (HERs), hydrodesulphurization, oxygen reduction reactions (ORRs), and methane conversion [16][17][18][19][20]. In addition, MoS 2 has been used for many years as a dry lubricant and an oil additive. New research also shows the possibility of improving MoS 2 dispersion by synthesizing MoS 2 particles on the surface of graphene materials [21]. These applications could be essential in future environmental projects involving sustainable energy sources, especially because the environmental impact of wet chemical synthesis is acceptable [22]. Owing to the many possible MoS 2 applications, the MoS 2 synthesis reaction kinetics must be accurately described, and the impact of the process conditions on the final product must be accurately modelled [23]. Most of the modern applications of molybdenum disulphide require a chemically pure material with reproducible properties. This wide range of applications makes molybdenum disulphide one of the most frequently studied materials in terms of its synthesis and properties. Particle size is of particular importance for catalysts and lubricants, but it is also important for applications to improve the properties of polymeric materials due to the necessity of achieving suitable dispersions. Recrystallisation or other processes related to the processing of the raw material are more commonly used in the production of electronic materials, but, even then, a substance with a known particle size distribution and morphology is essential. This type of material can be obtained by precipitation, also known as reactive crystallisation, in which the driving force behind the process is the supersaturation resulting from a chemical reaction.
To date, no studies modelling the effects of the reactor geometry and flow conditions on the precipitation of MoS 2 particles have previously been published. Therefore, this study presents the first such attempt and draws conclusions from the obtained results. The study findings shed new light on MoS 2 precipitated using wet chemical synthesis, which is particularly important for optimally designing reactor geometries to carry out the described process. This will be crucial for the application of this method to produce nanoparticles with the desired properties for industrial application.

Materials and Methods
The reaction model presented herein is based on our previous work [24]. The model describes the wet chemical synthesis of MoS 2 in impinging jet reactors (Figure 1), using ammonium heptamolybdate (HMA) and ammonium sulphide (AS) as reaction substrates in citric acid (CA), which acts as a catalyst and allows the reaction pH to be set. To enable reproducibility of the test, information on substrates and their preparation is described in the publications [24][25][26]. It should be noted that the aqueous solution of ammonium sulphide is unstable and degrades in contact with air, eventually leading to precipitation of free sulphur, which can completely distort the results, therefore clean and, above all, properly stored ammonium sulphide should always be used for the process. The model predicts only the particle nucleation and growth kinetics and does not allow for particle aggregation and agglomeration, which also occur extensively during the process. The reaction and particle morphology analysis were described in more detail in previous works [24][25][26]. Previously published experimental data were also used to validate the model [24,25]. Although the sequence of chemical reactions leading to the product has been presented previously [24,25,26], the last reaction given in the sequence was split into two (reaction Equations (8) and (9)). This is a more suitable way to describe the reaction mechanism separating dissociation from precipitation. Therefore, a modified reaction mechanism is given below, which does not affect the computational results or the modelling approach, but only shows a more likely course of the process.
Under ideal mixing conditions, the model given in Ref. [24] sufficiently described the process. However, in finite volume modelling, the model had to be modified by adding more equations because the original model inaccurately described phenomena at low molybdenum-ion concentrations. The process was previously described as being limited by the concentration of sulphide ions because free sulfur must precipitate to reduce molybdenum from the VI to the IV oxidation state. However, at low molybdenum-ion concentrations, the process is no longer limited solely by the sulphide ion concentration. Although the sequence of chemical reactions leading to the product has been presented previously [24][25][26], the last reaction given in the sequence was split into two (reaction Equations (8) and (9)). This is a more suitable way to describe the reaction mechanism separating dissociation from precipitation. Therefore, a modified reaction mechanism is given below, which does not affect the computational results or the modelling approach, but only shows a more likely course of the process. [ Under ideal mixing conditions, the model given in Ref. [24] sufficiently described the process. However, in finite volume modelling, the model had to be modified by adding more equations because the original model inaccurately described phenomena at low molybdenum-ion concentrations. The process was previously described as being limited by the concentration of sulphide ions because free sulfur must precipitate to reduce molybdenum from the VI to the IV oxidation state. However, at low molybdenum-ion concentrations, the process is no longer limited solely by the sulphide ion concentration. The concentration of molybdenum ions in the VI oxidation state then becomes essential to the process, which requires two precipitation models (Equations (2) and (3) in Table 1) combined using share functions depending on the molybdenum concentration (Equations (6) and (7) in Table 1) to obtain an accurate mathematical description of MoS 2 precipitation. The share functions for the HMA concentration are shown in Figure 2.
Rate of formation from MoS 2 precipitation (2) Rate of formation from S precipitation (3) Linear growth rate coefficient from S precipitation (5) u MoS 2 = exp −C HMA a MoS 2 precipitation kinetics share function in model (6) u S = 1 − exp −C HMA a S precipitation kinetics share function in model (7) Substrate consumption rate (10)   The model constants were fitted based on the experimental data described in Ref.
[24] using nonlinear programming algorithms described in the same publication and assuming ideal mixing conditions in a plug-flow reactor. A similar method of determining the model constants was proposed by Al-Tarazi et al. [27,28]. The model constants are listed in Table 2. The resulting model was further implemented in the CFD code using a user-  The model constants were fitted based on the experimental data described in Ref. [24] using nonlinear programming algorithms described in the same publication and assuming ideal mixing conditions in a plug-flow reactor. A similar method of determining the model constants was proposed by Al-Tarazi et al. [27,28]. The model constants are listed in Table 2. The resulting model was further implemented in the CFD code using a user-defined function (UDF) in Ansys Fluent software. Table 2. Constants for the advanced kinetics model.

Constant
Value Unit

CFD Modeling
Ansys Fluent 2020R1 software was used to solve steady-state system transport equations in 3D geometry model since it is impossible to model the V-type reactor directly using two-dimensional domain. Calculations were performed based on steady-state conditions without using the micromixing closure model. To predict the chemical reaction course, the flow must initially be calculated in the computing domain. Therefore, due to transitional conditions of the flow in impinging jet reactors [12,13], the shear stress transport (SST) k-ω turbulence model was employed, and the diffusivity was defined as follows: where µ T is the turbulent viscosity, Sc T is the turbulent Schmidt number, and the molecular diffusivity is constant: The SST k-ω turbulence model contains two transport equations. Because the standard k-ω and modified k-ε models are used near walls and in free flow, respectively, the SST k-ω model performs well for strong pressure gradients [29]. Both models are combined using a mixing function. As described by Bardina et al. [29], the original k-ω and transformed k-ε models are multiplied by the F 1 and 1 − F 1 mixing functions, respectively, and the equations from both models are summed. F 1 = 1 near the wall and 0 elsewhere in the region.
Turbulent viscosity is a function of the turbulent kinetic energy (k) and the specific energy dispersion rate (ω). To better estimate the boundary layer detachment, it is further modified using an auxiliary function defined as follows: Molecules 2022, 27, 3943 For turbulent kinetic energy, the transport equation is given as follows: ∂ρk ∂t For an energy dissipation rate, the transport equation is as follows: ∂ρω ∂t Model constants , β 2 , and γ 2 can be found in Ref. [29]. The boundary conditions are the same as those for the standard k-ω model.
The standard method of moments (SMM) was used to solve the population balance. The basic population balance can be formulated as follows: Using the moments of the number distribution of a crystal size and the equations of nucleation kinetics and crystal growth, the following balance equations for averaged moments can be formulated: where the equations for the average values of nucleation rate R N , growth rate G, and the considered moments m n are given in Table 1 (Equations (8), (9), and (13)-(15), respectively) and B n = D n = 0. The total average value for nucleation and growth are combinations of two kinetics equations, as described in Table 1 (Equations (8) and (9)). To simultaneously consider the effects of the mixture viscosity and density on the flow and the influence of the flow on the precipitation course, the CA concentration was assumed to most strongly influence the fluid properties. Therefore, the parameters entirely depend on the CA concentration. Because the particle concentration never exceeds 3% vol. in the suspension, the particle quantity negligibly affects the mixture rheology [30]. The particle volumetric concentration is represented by the third particle-decomposition moment described in Section 4, and the particles synthesized from the reaction are amorphous composites composed of sulfur and MoS 2 and, therefore, do not exhibit typical MoS 2 lubricating properties. The mixture viscosity was calculated using Kendall and Monroe's formula as follows: The first compound in the formula is responsible for the CA viscosity, while the second represents the remaining solution-which viscosity was assumed constant-independent of concentration and the same as for water. Using the same principle, the mixture density was calculated. The CA viscosity and density data were originally described in [31] and are listed in Table 3.
Using the determined mixture physicochemical parameters, Reynolds number values were calculated for the reactor inlet and outlet, as shown in Figure 3. Although the reactor inlet and outlet flow are both in the laminar regime, violent collisions and changes in the flow direction cause considerable turbulence in the reaction zone [13]. Therefore, a proper turbulence model must be developed to model the flow profile.
As shown in Figure 4, the computational domain was divided into the following zones to achieve the correct time scale: (1) inlet zone of HMA and CA mixture, (2) inlet zone of AS solution, (3) reactor outlet zone, (4) transition zone between zones 1-3 and reaction zone, and (5) reaction zone. Except for the inlets, which both exhibited the same mesh, each zone was characterized by a higher mesh density than the previous one. The reaction zone was filled with poly-hexcore mesh. The mesh parameters used for both computational domains are listed in Table 4. The first compound in the formula is responsible for the CA viscosity, while the second represents the remaining solution-which viscosity was assumed constant-independent of concentration and the same as for water. Using the same principle, the mixture density was calculated. The CA viscosity and density data were originally described in [31] and are listed in Table 3. Using the determined mixture physicochemical parameters, Reynolds number values were calculated for the reactor inlet and outlet, as shown in Figure 3. Although the reactor inlet and outlet flow are both in the laminar regime, violent collisions and changes in the flow direction cause considerable turbulence in the reaction zone [13]. Therefore, a proper turbulence model must be developed to model the flow profile.  As shown in Figure 4, the computational domain was divided into the following zones to achieve the correct time scale: (1) inlet zone of HMA and CA mixture, (2) inlet zone of AS solution, (3) reactor outlet zone, (4) transition zone between zones 1-3 and reaction zone, and (5) reaction zone. Except for the inlets, which both exhibited the same mesh, each zone was characterized by a higher mesh density than the previous one. The reaction zone was filled with poly-hexcore mesh. The mesh parameters used for both computational domains are listed in Table 4.    The mesh of average cell size of 0.02 mm was used in the impingement zone, and the y+ values at the walls in the reaction zone were around 0.1. It was tested that such a grid size was sufficient to obtain mesh independence, i.e., there was no difference in the average values of turbulence quantities (dissipation rate, kinetic energy) in the mixing zone for this and higher mesh densities. However, to reduce the influence of large scalar gradients on the computational results, the grid in the reaction zone was further refined-the particle growth rate was used as the adaptation criterion. An example of the adaptation results is shown in Figure 5. Because the nucleation and growth model constants were chosen assuming ideal reactor mixing conditions [24], the model had to be corrected to improve the prediction accuracy. Because ideal mixing conditions do not occur in real-life reactors, the constants do not fully reflect the phenomena occurring therein. Therefore, an additional function was introduced to correct the growth rate based on the local Reynolds number. We found that under certain flow conditions, the ideal mixing model predicts MoS2 particle sizes very well for higher Reynolds numbers. Those conditions were assumed to be close to ideal Because the nucleation and growth model constants were chosen assuming ideal reactor mixing conditions [24], the model had to be corrected to improve the prediction accuracy. Because ideal mixing conditions do not occur in real-life reactors, the constants do not fully reflect the phenomena occurring therein. Therefore, an additional function was introduced to correct the growth rate based on the local Reynolds number. We found that under certain flow conditions, the ideal mixing model predicts MoS 2 particle sizes very well for higher Reynolds numbers. Those conditions were assumed to be close to ideal mixing-for which the correction factor is equal or close to unity. Therefore, the correction factor was defined as the quotient of the local Reynolds number to the Reynolds number obtained under near-ideal mixing conditions. The correction factor is plotted as functions of the AS inlet Reynolds number in Figure 6, and the maximum correction factor was unity. The average local Reynolds number was calculated for mesh cells in which the energy dissipation rate was 10 times higher than the average in the computational domain, which corresponds to the collision zone area. Therefore, the growth rate is defined as follows: where the correction factor (ζ) is calculated using the formula: Molecules 2022, 27, x FOR PEER REVIEW 10 of 21 Clearly, the correction factor varies as a function of AS concentration, which corresponds to the mixture viscosity and density at different CA concentrations. The curves indicate the systemic transition flow most likely from the laminar to the turbulent regime. Furthermore, the correction factor is clearly closer to unity for higher flows and Reynolds numbers. Simulation convergence was obtained after tens of hours (2-3 days) for a single case, that is, modeling of the flow and precipitation process. Such computing times were obtained on AMD Ryzen 3700X CPU unit (using four cores, up to 4.0 GHz).

Results and Discussion
The modelling results showed that because the different reactor flows depended on the mixture viscosity, the mixing conditions influenced the reaction course. Therefore, the most critical parameters affecting mixing and MoS2 precipitation were compared for both reactor types. Figure 7 shows the influence of the reactor geometry on the inert-tracer concentration distribution, revealing that because of additional axial mixing, the V-type reactor exhibited better mixing conditions than the T-type one. In particular, this can be seen in Figure  7c, which shows the change in concentration of the nonreacting tracer along the radius for two distances from the bottom of the reactor. For the vortex reactor, the concentration is more uniform and closer to the average value than with the T-reactor. Better mixing was achieved at lower substrate concentrations-and, hence, lower mixture viscosities-than at higher concentrations. The velocity profiles shown in Figure 8 are characteristic for the Clearly, the correction factor varies as a function of AS concentration, which corresponds to the mixture viscosity and density at different CA concentrations. The curves indicate the systemic transition flow most likely from the laminar to the turbulent regime. Furthermore, the correction factor is clearly closer to unity for higher flows and Reynolds numbers.
Simulation convergence was obtained after tens of hours (2-3 days) for a single case, that is, modeling of the flow and precipitation process. Such computing times were obtained on AMD Ryzen 3700X CPU unit (using four cores, up to 4.0 GHz).

Results and Discussion
The modelling results showed that because the different reactor flows depended on the mixture viscosity, the mixing conditions influenced the reaction course. Therefore, the most critical parameters affecting mixing and MoS 2 precipitation were compared for both reactor types. Figure 7 shows the influence of the reactor geometry on the inert-tracer concentration distribution, revealing that because of additional axial mixing, the V-type reactor exhibited better mixing conditions than the T-type one. In particular, this can be seen in Figure 7c, which shows the change in concentration of the nonreacting tracer along the radius for two distances from the bottom of the reactor. For the vortex reactor, the concentration is more uniform and closer to the average value than with the T-reactor. Better mixing was achieved at lower substrate concentrations-and, hence, lower mixture viscosities-than at higher concentrations. The velocity profiles shown in Figure 8 are characteristic for the chosen reactor types.  The diffusivity contour plots show the effect of the high Reynolds number-resulting from the increased flow rate-on increasing diffusivity ( Figure 9). Clearly, the increased diffusivity is more pronounced in the V-type reactor and at higher flow rates, and the maximum diffusivity in the T-type reactor is lower than that in the V-type one. The opposite is true for lower flow rates because of the different CA distributions within the reactors. Therefore, the CA concentration critically affects the mixture viscosity and density. At higher CA concentrations, diffusivity decreased because the mixture viscosity increased. Furthermore, the T-type reactor collision axis clearly shifts toward the HMA and CA mixture inlet because the fluids fed into the reactor exhibit different inertias. The diffusivity contour plots show the effect of the high Reynolds number-resulting from the increased flow rate-on increasing diffusivity ( Figure 9). Clearly, the increased diffusivity is more pronounced in the V-type reactor and at higher flow rates, and the maximum diffusivity in the T-type reactor is lower than that in the V-type one. The opposite is true for lower flow rates because of the different CA distributions within the reactors. Therefore, the CA concentration critically affects the mixture viscosity and density. At higher CA concentrations, diffusivity decreased because the mixture viscosity increased. Furthermore, the T-type reactor collision axis clearly shifts toward the HMA and CA mixture inlet because the fluids fed into the reactor exhibit different inertias.
The nucleation rate contours are shown in Figure 10, which compares two types of reactor types at the same flow rate and the same concentrations. Clearly, with increasing AS concentration, the maximum nucleation rate increases over a limited area in the V-type reactor. At the same flow rate and concentration, on the other hand, nucleation occurs over a broader area in T-type than in V-type reactors because less mixing occurs in T-type reactors. Higher nucleation rates were observed at higher flow rates when comparing the same reactor type and concentration.
Comparing the same reactor and concentration at different flow rates, considerably higher growth rates were obtained at the 50 mL/min flow rate (Figure 11). Similarly, higher growth rates were obtained in the V-type reactor than in the T-type one (Figure 12) because the reactants were better mixed in the V-type reactor. Like the nucleation rate contour plots, higher growth rates were also observed at higher concentrations. However, the particles grew over a smaller area.
Although the characteristic particle sizes shown in Figure 13 were in good agreement with those predicted by the ideal mixing model, the fit was slightly worse when the particle size was modelled using CFD. As can be seen, the results of nucleation and growth modelling including the flow field gives a slightly poorer fit than modelling using a pipe reactor model under ideally mixed conditions. However, it allows the effects of changes in geometry and flow conditions on the nucleation and growth process to be taken into account, and therefore, can be used to optimize or scale up the process or study it in different geometries. The model of perfect mixing cannot be used for this purpose, thus its applicability is very limited (or even its not applicable) in real-life engineering applications. The nucleation rate contours are shown in Figure 10, which compares two types of reactor types at the same flow rate and the same concentrations. Clearly, with increasing AS concentration, the maximum nucleation rate increases over a limited area in the Vtype reactor. At the same flow rate and concentration, on the other hand, nucleation occurs over a broader area in T-type than in V-type reactors because less mixing occurs in T-type reactors. Higher nucleation rates were observed at higher flow rates when comparing the same reactor type and concentration. The nucleation rate contours are shown in Figure 10, which compares two types of reactor types at the same flow rate and the same concentrations. Clearly, with increasing AS concentration, the maximum nucleation rate increases over a limited area in the Vtype reactor. At the same flow rate and concentration, on the other hand, nucleation occurs over a broader area in T-type than in V-type reactors because less mixing occurs in T-type reactors. Higher nucleation rates were observed at higher flow rates when comparing the same reactor type and concentration. Comparing the same reactor and concentration at different flow rates, considerably higher growth rates were obtained at the 50 mL/min flow rate (Figure 11). Similarly, higher growth rates were obtained in the V-type reactor than in the T-type one (Figure 12) because the reactants were better mixed in the V-type reactor. Like the nucleation rate contour plots, higher growth rates were also observed at higher concentrations. However, the particles grew over a smaller area.  Although the characteristic particle sizes shown in Figure 13 were in good agreement with those predicted by the ideal mixing model, the fit was slightly worse when the particle size was modelled using CFD. As can be seen, the results of nucleation and growth modelling including the flow field gives a slightly poorer fit than modelling using a pipe reactor model under ideally mixed conditions. However, it allows the effects of changes in geometry and flow conditions on the nucleation and growth process to be taken into account, and therefore, can be used to optimize or scale up the process or study it in different modelling including the flow field gives a slightly poorer fit than modelling using a pipe re-actor model under ideally mixed conditions. However, it allows the effects of changes in geometry and flow conditions on the nucleation and growth process to be taken into account, and therefore, can be used to optimize or scale up the process or study it in different geometries. The model of perfect mixing cannot be used for this purpose, thus its applicability is very limited (or even its not applicable) in real-life engineering applications.  Figure 14 shows the parity plots obtained for the results of the tube reactor calculations assuming ideal mixing conditions and the CFD numerical modelling under nonideal mixing conditions. It can be seen that, regardless of the calculation method, the results are approximately on the perfect fit line. In the case of ideal mixing, the calculations have a clear tendency to overestimate the particle size, whereas in the case of CFD modelling, the average size is undervalued. Figure 13. Characteristic particle sizes for V-type reactor at (a) 20 and (b) 50 mL/min flow rates and for T-type reactor at (c) 20 and (d) 50 mL/min flow rates (ideal mixing refers to the simulation of a tubular reactor with ideal mixing conditions, and non-ideal mixing refers to CFD modelling that takes into account the reactor geometry and flow conditions). Figure 14 shows the parity plots obtained for the results of the tube reactor calculations assuming ideal mixing conditions and the CFD numerical modelling under non-ideal mixing conditions. It can be seen that, regardless of the calculation method, the results are approximately on the perfect fit line. In the case of ideal mixing, the calculations have a clear tendency to overestimate the particle size, whereas in the case of CFD modelling, the average size is undervalued.
Because the determined particle moment distributions show a narrow zero-moment fit (Figure 15), the numerical model exhibits limited ability to accurately predict the number of particles. However, because a better fit was obtained for higher particle moments, the model predicted the resulting particle sizes with good accuracy. Much more accurate numbers of particles were predicted assuming ideal mixing conditions for the model kinetics [24].
The results indicate that the adopted methodology can be successfully applied to predict the precipitation of nanoparticles even in the face of very fast chemical reactions such as the synthesis of metal sulphides. The precipitation of molybdenum sulphide is shown as an example of a very complex process, which is possible to model and predict particle size. Because the determined particle moment distributions show a narrow zero-moment fit (Figure 15), the numerical model exhibits limited ability to accurately predict the number of particles. However, because a better fit was obtained for higher particle moments, the model predicted the resulting particle sizes with good accuracy. Much more accurate numbers of particles were predicted assuming ideal mixing conditions for the model kinetics [24].
The results indicate that the adopted methodology can be successfully applied to predict the precipitation of nanoparticles even in the face of very fast chemical reactions such as the synthesis of metal sulphides. The precipitation of molybdenum sulphide is shown as an example of a very complex process, which is possible to model and predict particle size.  Because the determined particle moment distributions show a narrow zero-moment fit (Figure 15), the numerical model exhibits limited ability to accurately predict the number of particles. However, because a better fit was obtained for higher particle moments, the model predicted the resulting particle sizes with good accuracy. Much more accurate numbers of particles were predicted assuming ideal mixing conditions for the model kinetics [24].
The results indicate that the adopted methodology can be successfully applied to predict the precipitation of nanoparticles even in the face of very fast chemical reactions such as the synthesis of metal sulphides. The precipitation of molybdenum sulphide is shown as an example of a very complex process, which is possible to model and predict particle size.

Conclusions
An improved model of particle nucleation and growth kinetics has been implemented in CFD modelling for MoS2 particles synthesized using confined impinging jet reactors. Because the improved model can predict particle size distributions with satisfactory accuracy, it can provide a valuable engineering tool for designing chemical reactors to produce MoS2 particles. Parity plots of characteristic particle sizes clearly show that the model can accurately predict particle size distributions regardless of the apparatus geometry.


The model enables kinetic constants to be determined for other complex chemical reactions, even with limited knowledge of the reaction mechanism, and can be applied with CFD to various reaction processes;  The effect of the mixing conditions on the chemical reaction was determined using the SST k-ω model combined with the developed kinetic model to calculate results close to the experimental ones. In addition, the modelling and experimental results deviated more markedly at higher concentrations and lower flow rates, which may be because the reactant mixing was worse and, consequently, deviated even further from the ideal mixing conditions. Under these conditions the viscosity is higher due

Conclusions
An improved model of particle nucleation and growth kinetics has been implemented in CFD modelling for MoS 2 particles synthesized using confined impinging jet reactors. Because the improved model can predict particle size distributions with satisfactory accuracy, it can provide a valuable engineering tool for designing chemical reactors to produce MoS 2 particles. Parity plots of characteristic particle sizes clearly show that the model can accurately predict particle size distributions regardless of the apparatus geometry.

•
The model enables kinetic constants to be determined for other complex chemical reactions, even with limited knowledge of the reaction mechanism, and can be applied with CFD to various reaction processes; • The effect of the mixing conditions on the chemical reaction was determined using the SST k-ω model combined with the developed kinetic model to calculate results close to the experimental ones. In addition, the modelling and experimental results deviated more markedly at higher concentrations and lower flow rates, which may be because the reactant mixing was worse and, consequently, deviated even further from the ideal mixing conditions. Under these conditions the viscosity is higher due to the higher concentration of citric acid, what strongly affects the flow parameters, and thus the mixing-limited reaction; • The concentration and velocity contour plots indicate that the V-type reactor exhibited superior fluid mixing than the T-type one at the same flow rate. However, although this is associated with a more marked pressure drop, the particles were better mixed and nucleated over a much larger area. Fluid mixing is also affected by reagent concentration because fluid viscosity increases with increasing reagent concentration, which affects the Reynolds number.
Future studies may elucidate the reaction mechanism of the wet chemical synthesis of MoS 2 , which will enable the development of a more physical nucleation and growth model. Additionally, particle aggregation, agglomeration, and aggregate and agglomerate disintegration should be investigated in the future to determine their effect on the precipitation process and possible separation of particles from the suspension.
Due to the limited amount of rare earth elements, technologies related to their processing will become increasingly important in the future, making these results important for the production of metal sulphides. Additionally, there are valuable materials with catalytic properties that are likely to find application in the global energy transition.