Particulate Matter Dispersion Modeling in Agricultural Applications: Investigation of a Transient Open Source Solver

Agriculture is a major emitter of particulate matter (PM), which causes health problems and can act as a carrier of the pathogen material that spreads diseases. The aim of this study was to investigate an open-source solver that simulates the transport and dispersion of PM for typical agricultural applications. We investigated a coupled Eulerian–Lagrangian solver within the open source software package OpenFOAM. The continuous phase was solved using transient large eddy simulations, where four different subgrid-scale turbulence models and an inflow turbulence generator were tested. The discrete phase was simulated using two different Lagrangian solvers. For the validation case of a turbulent flow of a street canyon, the flowfield could be recaptured very well, with errors of around 5% for the non-equilibrium turbulence models (WALE and dynamicKeq) in the main regions. The inflow turbulence generator could create a stable and accurate boundary layer for the mean vertical velocity and vertical profile of the turbulent Reynolds stresses R11. The validation of the Lagrangian solver showed mixed results, with partly good agreements (simulation results within the measurement uncertainty), and partly high deviations of up to 80% for the concentration of particles. The higher deviations were attributed to an insufficient turbulence regime of the used validation case, which was an experimental chamber. For the simulation case of PM dispersion from manure application on a field, the solver could capture the influence of features such as size and density on the dispersion. The investigated solver is especially useful for further investigations into time-dependent processes in the near-source area of PM sources.


Introduction
Particulate matter consists of microscopic particles, which are suspended in the earth's atmosphere. These microscopic particles can enter the human respiratory system or, depending on their size, can even enter the blood stream and cause severe health problems. These are classified as PM1, PM2.5, PM10 based on the particle sizes with all particles less than 1 µm, 2.5 µm and 10 µm, respectively. The effects of particulate matter exposure may result in premature mortality, aggravation of respiratory and cardiovascular diseases, asthma, chronic bronchitis, decreased lung function and various respiratory and cardiovascular health problems [1]. According to the world health organisation, 4.2 million premature deaths each year can be connected to an excessive exposure to PM [2]. and immission processes. Also, the discrete phase is simulated with respective turbulence models to take into account turbulent diffusion. Probably due to the high computational costs, little research in the use of transient open source solvers with a direct coupling of the particle phase and the gaseous phase for near-source agricultural applications was done so far. To close this gap, this study aims to investigate the use of an transient open source solver for the simulation of transport processes of PM in the near field of typical agricultural applications. The goal was to (1) validate the velocity field (2) validate the transport of particles, (3) validate a boundary condition for the cost-effective way to generate a turbulent wind field, and (4) demonstrate the validated solver with an example application, namely the the application of manure on a field.
The paper is organized as follows. Section 2 gives an overview on the theory of large eddy simulations and particle transport. Section 3 describes the validation process of the solver and presents the results from the validation (goals 1 and 2). Section 4 describes the application of the solver for the manure application and presents the simulated results (goals 3 and 4). Section 5 discusses these results and Section 6 presents the drawn conclusions. The workflow of this study is briefly illustrated in Figure 1.

Materials and Methods
The following section provides a brief overview about the solving techniques for the (continuous) flow field and the (dispersed) particle phase. Two different turbulence modelling approaches, namely, RANS and LES, were applied to calculate the turbulent flow fields. Based on these flow fields, the transport of the particulate matter is predicted by two different Eulerian-Langrangian solvers implemented in the open source CFDframework OpenFOAM. Moreover, the complete numerical setup (e.g., computational domain, boundary conditions, solver setups) and mesh convergence study, as well as the validation of the numerical results, will be discussed in the following. To provide a good reproducibility, all OpenFOAM-related terms and specific settings will be shown in typewriter font.

Solving the Flow Field Phase
Large Eddy Simulations is a method of turbulence modelling focusing on the behaviour of large eddies by spatially filtering the transient Navier-Stokes equations. LES reduces the computational cost by ignoring the eddies with smallest length scales. In LES, the vortices are subdivided into large resolved vortices, which transport most of the energy in the flow, and small subgrid-scale (SGS) vortices, which are modeled with an SGS turbulence model. The separation of these large and small vortices is done using spatial filtering. Filtering the Navier-Stokes equations results in Equation (1), where τ SGS ij = u i u j −ũ iũj is the SGS stress, which considers the effect of small vortices.
These stresses are calculated using SGS turbulence models. One of the most popular SGS turbulence models is the Smagorinsky model [23]. This uses the idea of Prandtl mixing theory and calculates the SGS stress as According to Smagorinsky, the subgrid viscosity is calculated as whereS ij is the strain rate tensor and C s is the constant of Smagorinsky. The constant of Smagorinsky was theoretically estimated in the range of 0.17-0.21. The Smagorinsky model is a zero-equation model and is based on the local equilibrium assumption. Even though the model is simple and stable, laminar flow is not modelled, damping is necessary for near wall regions and the Smagorinsky constant is problem-dependent.
To improve the simulation results for wall-bounded flows, the Wall-adapting Local Eddy viscosity (WALE) model was proposed considering the square of the velocity gradient tensor as where C w is the model constant in the range of 0.55 ≤ C w ≤ 0.60, which is appropriate for C s = 0.18, and S d ij is the traceless symmetric part of the square of the velocity gradient tensor [24]. The WALE model is able to handle transition and damping is not necessary for the near-wall region.
Smagorinsky, and WALE are zero equation models, whereas the model proposed by Yoshizawa and Horiuti is the one-equation eddy viscosity model, which additionally solves a transport equation to compute subgrid kinetic energy k SGS . According to this model, the subgrid scale viscosity is computed as in terms of the subgrid kinetic energy [25,26]. Similar to the Dynamic Smagorinsky model, a dynamic one equation eddy viscosity model was proposed by using a dynamic procedure to calculate the coefficients. In these one-equation models, no assumption was made of local equilibrium between the subgrid-scale energy production and dissipation rate. These one-equation models tend to produce more accurate results than the zero-equation models [27]. The used software Package OpenFOAM contains a library with a broad range of SGS scale turbulence models. For this study, we investigated the models of Smagorinsky, TKE, dynamic TKE, and the wall-adapting local eddy viscosity (Smagorinsky, oneEqTKE, dynTKE and WALE, respectively).
Another way to solve the flow field is the Reynolds Averaged Navier Stokes (RANS) equations approach. With the RANS equations, only the mean flow is solved, while all turbulence scales are parametrized. For the turbulence parametrization, a turbulence model is applied. A comprehensive overview of the turbulence models within the RANS approach is given in [14].

Solving the Particle Phase
Particle modelling in CFD is considered as multiphase approach since the particles in the fluid becomes a separate phase. In CFD, the multiphase modelling can be done either by Euler-Euler or Euler-Lagrange approach. In the Euler-Euler approach, the particles are considered as a continuous phase whereas in the Euler-Lagrange approach, the particles are considered as a discrete phase and Newton's equations for motion are solved for the particle phase to track the trajectories of the particle. The Euler-Euler approach is suitable for large particles, where the particle size is greater than the cell-size of the mesh. The Euler-Lagrange approach is suitable for smaller particles [28].
As mentioned above, in the Lagrangian framework, the motion of the particles is calculated by solving a set of ordinary differential equations along the particle trajectories. According to Newton's second law of motion, this requires the consideration of all relevant forces acting on the particle, e.g., drag force, gravitational and buoyancy force or pressure gradient force. Since analytical solutions for the majority of the particle forces are only available for small particle Reynolds numbers (i.e., Stokes regime), the determination of the considered particle forces is based on different empirical correlations derived from experiments. In most particle-laden flows, the drag force dominates for the particle motion [29] and is expressed in terms of a drag coefficient C D , which allows for the calculation of the drag force not only for small-particle Reynolds numbers (Re p < 0.5), where viscous effects dominate and no separation is observed, but also for the transition region (i.e., 0.5 < Re p < 1000) and fully turbulent region or Newton regime (above Re p ≈ 1000). The applied drag models are based on the particle Reynolds number, which is defined as: with the velocity u f ,i , density ρ f and the dynamic viscosity µ f of the fluid (or continuous phase), the particle diameter D p and particle velocity u p,i , as well as the magnitude of the relative slip velocity |u f ,i − u p,i |. Assuming rigid and sparsely distributed spherical particles (i.e., dilute dispersion), the drag coefficient is determined using the drag model proposed by Putnam [30]: which is suitable for high particle Reynolds numbers (Re p ≤ 1000) and ensures the correct limiting behavior within the Newton regime (Re p > 1000). Based on the calculated drag coefficient, the general equation of particle motion is used to evaluate the corresponding drag coefficient for spherical particles as follows: Besides the Putnam drag model, various other drag models exist to account for nonspherical particles or different fluid/particle volume fractions. A widely used drag model for higher particle volume fractions φ p (i.e., dense dispersed two-phase flows), which is usually the case in fluidization, is the Ergun-Wen-Yu model. This combines the well-known Ergun equation (see, e.g., [29,31]) for fluid flows in packed beds and the drag model for spherical particles proposed by Wen and Yu [32]. Thus, for fluid volume fractions φ f smaller than 0.8, the particle drag force is directly calculated from the Ergun equation, whereas, for higher fluid volume fractions (i.e., φ f ≥ 0.8), the drag coefficient can be estimated by the Wen-Yu model [31,32]: which is the drag coefficient formulation originally proposed by Schiller and Naumann [29], but accounting for the fluid volume fraction φ f . Analogously to the application of Putnam's correlation, Equation (7), the estimated drag coefficient is then used to calculate the particle drag force from the general equation of particle motion, Equation (8).
One of the most important parameters to characterize fluid-particle flows is the Stokes number St = τ p /τ f , which relates the particle response time τ p to the characteristic time of the flow field τ f . If the Stokes number is less than < 0.1, the particles follow the flow instantaneously [33], and as the Stokes number increases, the particles start to deviate from the flow. In particle modeling, it is also important to consider the phase coupling to decide how the particles interact with each other and the fluid. Elgobashi [34] classified particle-laden turbulent flows based on the volume fraction of particles φ p . For values of φ p ≤ 10 −6 the particles follow a one-way coupling where the momentum transfer from particles to the flow is negligible. For values of 10 −6 < φ p ≤ 10 −3 , the momentum transfer is large and affects the flow, which is termed a two-way coupling. In this region, when the particle response time τ p decreases, the dissipation rate increases, whereas an increase in τ p increases the production of turbulence energy. When φ p > 10 −3 , the flow is dense and particle-particle collisions occur, along with a momentum transfer between the particles and fluid, which is termed a four-way coupling.

Validation of the Flow Solver
To validate the flow solver, an urban street canyon study in an atmospheric boundary layer wind tunnel by [35] was used. This study was chosen because it provided very detailed information on the flowfield in a fully turbulent environment with near obstacles, which is similar to a large number of applications for near field particle modeling in agricultural applications. Wind tunnel experiments were carried out in a test section of width 2 m, height 1 m and length 2 m, with a boundary layer flow with mean velocity profile exponent α = 0.30.

Meshing
Three different meshes were generated, all consisting of purely hexahedral cells with the OpenFOAM integrated mesh generators blockMesh and snappyHexMesh. At first, the domain was meshed with more cells near the bottom wall by grading the cells along the height of the domain. The domain was meshed with 260 × 92 × 70 (coarse), 313 × 110 × 70 (medium) and 375 × 130 × 85 cells (fine). Mesh refinement was performed around the buildings and the canyon formed by the buildings. A single-refinement region was created for the coarse mesh and the mesh was refined two times for medium and fine meshes.
The meshes consisted of 2,103,476 computational cells for the coarse mesh, 5,038,022 cells for the medium mesh and 10,559,532 cells for the fine mesh. The y + values for the meshes are described in Table 1. Even with the finest mesh and inflation layers at the walls, a y + < 1 could not be achieved in all regions; hence, a wall function using Spalding's law of the wall was applied, which is known to perform well in all regions of y + [36,37].

Boundary Conditions
The inlet boundary conditions were chosen as suggested by C. Gromke et al. [35]. The inlet velocity profile follows the power law approach as in Equation (10) (11) with boundary layer depth δ, friction velocity u τ = 0.52 m/s and C µ = 0.09.
The pressure at the inlet was defined with a Dirichlet boundary condition (zeroGradient). At the outlet, the boundary condition for velocity was set to zeroGradient and the pressure to fixedValue. The top, front and back boundaries are set as symmetryPlane, meaning a free-slip condition. The boundaries' bottom and buildings are considered walls, where the velocity is set to zero.

Solver Setup and Turbulence Modeling
Incompressible Large Eddy simulations were set up using the PIMPLEFOAM solver, which is a combination of the Pressure-Implicit with Splitting of Operators, [38] (PISO) solver and the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm, [39]. It has the advantage of a fast convergence for transient simulations, as described, e.g., by [40]. An adaptive time stepping was chosen with a Courant number no higher than 1.
The spatial discretization was performed using second-order linearUpwind schemes, and the time variable was discretized with the second-order backward scheme, which are both described in [41].
In this study, the influence of four different subgrid scale turbulence models on the flow field were investigated, namely kEqn, dynamicKEqn, WALE, and Smagorinsky, all described in Section 2.1.
The simulation was carried out for 5 s with the help of 24 processors running for 6 days. Considering the flow velocity at the height H, U H = 4.7 m/s, the flow passes through the entire domain five times in 5 s and, in the region of interest, near the buildings, the flow passes through approximately 40 times. Hence, a period of 5 s was considered sufficient for the entire simulation.

Mesh Convergence
Mesh convergence study was carried out using the dynamicKEqn turbulence model and the experimental results [35] at three different points in the street canyon along a horizontal line on the xz-plane at y = 0 (symmetry plane) and a height of z = 0.7 H. The points are chosen from the experimental results where the normalized vertical velocities U z /U H are −0.05, −0.1 and −0.15. Hence, the normalized vertical velocities are measured from the simulations at x/H = 0.1425, 0.2 and 0.2525 and compared with the experimental results, as in Figure 3. From the figure, the large influence of the grid size can clearly be seen. At x/H = 0.2525, the relative errors for coarse, medium and fine meshes are 37%, 15% and 4%, respectively. Hence, for all the following simulations, the fine mesh was chosen.  Figure 4 shows the flowfields simulated with the four different SGS turbulence models. The normalized vertical velocity is shown in isotachs in the XZ-area between building A and building B at y/H = 0.5 and compared with the experimental values. All four SGS models predicted the general pattern of a clockwise vortex formed in the canyon, with negative vertical velocities on the side of building B and positive velocities at building A. A clear difference can be seen between the equilibrium-assumption SGS models kEqn and Smagorinsky in Figure 4a,b and the non-equilibrium SGS models dynkEqn and WALE in Figure 4c,d. Both equilibrium models failed to predict the spatial distribution of velocity contours accurately, especially for the positive-velocity region near building A. Here, both models predicted the velocitiy maxima near the edges of the building whereas, in the experiment, the velocitiy maxima were high near the middle region of the building. On the other hand, the non-equilibrium models WALE and dynamicKEqn performed very well when reproducing the experimental results. On the side of building B, a very good agreement can be seen for the negative velocities towards the middle region, with no noticeable difference between experimental and simulation results for dynamicKEqn. However, near the wall region, all SGS models showed larger deviations from the experimental data. The highest deviations can be seen at the height of z/H = 0.6; the relative error was, depending on the compared velocity contour, in the range of 5-20%. Possible reasons for this are discussed in Section 5.1.

Validation of the Particle Transport
For the validation of the particle transport solver, high-resolution measurements of particulate matter concentrations (in both time and space) in a fully turbulent flow under well-documented experimental conditions would ideally be found. Unfortunately, these data cannot be found in the literature. As a compromise, the study of [42] was chosen as a reference case to validate against. [42] conducted experiments inside a model room of size 0.8 m × 0.4 m × 0.4 m; see Figure 5. Silver-coated hollow glass spheres of density 1400 kg/m 3 , with a mean diameter of 10 µm and nominal diameter range between 2 and 20 µm, were injected for a time period of 1800 s. Air-flow velocity and particle concentration were measured using Phase Doppler Anemometry (PDA). The particle concentration was measured at 27 different locations in the mid-plane of the model room along three different x positions.

Computational Domain and Meshing
The room was modeled and meshed, as sketched in in Figure 5. The centres of inlet and outlet were located at x = 0 m, y = 0.2 m, z = 0.36 m and x = 0.8 m, y = 0.2 m and z = 0.04 m, respectively, with the same area of 0.04 m × 0.04 m. Two different meshes were created with 221,184 (96 × 48 × 48) and 432,000 (120 × 60 × 60) hexahedral cells. Excepting the inlet and outlet, all other boundaries were considered as walls with a no-slip boundary condition. In both meshes, the y + was less than 1 in almost all regions except the corners, where maximum values of 6 were reached.

Particle Injection and Simulation
At the inlet, a uniform velocity was set according to the experimental inlet velocity of 0.225 m/s. The particles were injected from the inlet patch for 1800 s. Every second, 1000 particles were injected. To vary the diameter of the particles according to the experiment, a normal distribution model was used, with an expectation of 10 µm. The particle diameters ranged from 8 µm to 13 µm. The density of the particles was set to 1400 kg/m 3 . The particles were simulated with Gravity, Drag and Pressure forces. Boundary interactions were set to an escape condition when particles reach the outlet, or rebound condition when particles reached the walls.
For the simulation of the particle phase, the influence of different LES and URANS models on the particle concentration was investigated, as well as the influence of two different lagrangian solvers, namely pimpleLPTFoam, in combination with the drag model of Putnam (sphereDrag) and MPPICFoam using the Ergun-Wen-Yu drag model (ErgunWenYuDrag). The multi-phase particle in Cell method [43] is similar to the Discrete Parcel Method (DPM). In DPM, the particle interactions are directly resolved, whereas in MPPIC, particle-particle interactions are represented by models, which utilise the mean values calculated on the Eulerian mesh. pimpleLPTFoam was originally introduced by the Chair of Modeling and Simulation (LeMoS) at the University of Rostock [44,45] for dilute particle-laden flows by combining a DPM-based particle-tracking approach, including four-way phase-coupling capabilities, with the standard OpenFOAM solver pimpleFoam.
Ithen case of URANS, the converged RANS solution was mapped onto the fine mesh and then the particle phase was simulated with pimpleLPTFoam using the kEpsilon model for 1850 s. The particles were injected from 20 s until 1820 s. In the case of LES, the dynamicKEqn turbulence model was used, since it produced better results in the flowfield validation case, as shown in Section 3.1.6.
Both pimpleLPTFoam and MPPICFoam were set up with backward time schemes; a second-order linear scheme for div(phi,U) term and linearUpwind for div(phi,k) term were used. The simulation was carried out for 1850 s with 24 processors, with an average simulation time of 40 h for LES.

Validation Results of the Particle Transport
The velocity profile U x from the experiment is compared to the RANS and LES simulations along three lines at x = 0.2, 0.4 and 0.6 m in the mid-plane in Figure 6. The velocity profile from the RANS simulation with the kEpsilon turbulence model predicted the experimental velocity profile very precisely along all three lines. For LES simulation, the velocity profile was time-averaged and then compared. Even though the LES simulation began from the converged solution of the RANS simulation, the velocity along the inlet tends to increase further than the expected experimental values. LES was able to perform well near the wall region and in the region below the inlet, but the velocity values are overpredicted along the inlet. The possible reasons for this are discussed in Section 5.2.
The particle distributions in the domain at four different times are shown in Figure 7. For the LES, the particles are thrown towards the wall near the outlet and hit the wall. Particles with a large diameter tend to fall down due to gravity soon after entering the domain, and then, due to the negative-velocity region below the inlet, they move back, as observed in the experiments. A mixing of particles over the whole volume is generally visible, with slightly higher concentrations below the inlet, along the inlet jet, and at the outlet positioned wall. In contrast, for the URANS, the particles always follow the same path as the flow field and no dispersion is visible in the volume of the domain. Even though the particles with a larger diameter tend to fall earlier within the jet, they only start to fall near the outlet wall. This very uneven distribution can be explained by solving the flow field, which converged after the initial RANS. Due to the lack of turbulence, further simulations with URANS result in stagnation of the flow field, meaning that the particle trajectory at 60 s is the same as the particle trajectory at 300, 600, 1800 s and throughout the simulation.

Modelling of Particulate Matter Dispersion from Manure Application
The solver was applied to simulate the dispersion of PM from a manure application process on a field. The experimental settings and detailed information on the source of PM is described in [10,11].

Computational Domain
The geometry is an empty cuboid with flow along the positive x-direction, as in Figure 9. The dimensions of the domain 500 m × 150 m × 80 m. The dimensions and the grid size were chosen after preliminary studies, which are not presented in this study for brevity. The domain height was chosen based on the available wind data. The front, back and top boundaries were considered as symmetryPlane, the bottom was considered as wall, and inlet and outlet as patch. Two meshes were considered for the simulation. The domain was first meshed with a small number of cells, and then then refined using the solver's in-built refineMesh utility. The cells were graded along the z-direction, so that the mesh near the bottom wall was fine, with a nearly uniform aspect-ratio. The mesh consisted of 19,995,000 hexahedral cells.

Turbulent Inlet
To generate fully turbulent inflow conditions in a cost-efficient way, the Divergence-Free Synthetic Eddy Method (DFSEM) was applied to produce a divergence-free turbulence field directly at the inlet of the computational domain. This approach was introduced by Poletto et al. [46] and represents a consequent development of the so-called Synthetic Eddy Method (SEM), which was proposed earlier by Jarrin et al. [47]. Accordingly, based on the Reynolds decomposition of the velocity field, the velocity fluctuations u i are considered as velocity perturbations caused by N turbulent structures randomly placed within a certain volume of space: where r k i = x i −x k i σ k , including x k i as the location of the center of the k-th eddy and the turbulence length-scale σ k calculated at the eddy center. Moreover, Equation (12) consists of a suitable shape function q 0 as well as the eddy intensity α k i , represented by random numbers with zero average. Further information can be found according to the formulation of the DSFEM in [46,48].
The DSFEM requires input parameters to generate the turbulent inflow, namely, the velocity vector, Reynolds stress tensor and integral length scale. To obtain these parameters, wind tunnel experiments were conducted in an atmospheric boundary layer wind tunnel. The inflow profile was described in a previous validation study [49]; the wind tunnel is described in detail in [50].
Two characteristics of the simulated wind field will be assessed. The first is the vertical profile of the time-averaged velocity in x-direction u x and its stability throughout the whole domain. The second is a turbulence characteristic, namely, the Reynolds stress R 11 , which is defined as the time-average of the square of the velocity fluctuations in x-direction: R 11 = u X · u x , with u x = u x − u x . As can be seen, R 11 is closely related to the turbulence intensity in x-direction. U x and R 11 were chosen because they are considered to be two of the most important characteristics for a modeled windfield in the guidelines for Physical modelling of flow and dispersion processes in the atmospheric boundary layer [51].

Solver Setup
To simulate the flow field, the validated dynamicKEqn turbulence model was used, with pimpleLPTFoam solver. The flow field was simulated for a period of 400 s without particles. The period of 400 s was chosen, assuming the average wind speed to be 5 m/s, so that the flow passes through the domain four times. Second-order linearUpwind schemes were used for the div(phi,U)and div(phi,k) terms. The case setup and the solver can be found and downloaded in [52] for further use by the reader. The simulation was performed with 420 processors on the SUPERMUC at Leibniz Supercomputing Centre running for 8 h with a maximum Courant number of 1. After 400 s, the particles were injected and the simulation was continued for further 120 s, which then required another 8 h of simulation time. During the particle phase, the results were written out every 5 s. Even though pimpleLPTFoam and MPPICFoam were able to produce similar results, as discussed in Section 5.2, pimpleLPTFoam was chosen as the solver, as it was convenient to use with the turbulentDFSEMInlet boundary condition.

Particle Properties
The size distribution and density of the particulate matter aerosolized from the manure were obtained from the wind tunnel experiments published in [10,11], where the measured density of the investigated manure was 316 kg/m 3 , and the density of the aerosolized particulate matter was estimated as 130 kg/m 3 . The size distribution of the particulate matter was measured and is shown in Figure 10.

Normalized Size Distribution
Particle Diameter (m) Figure 10. Particle size distribution injected in the application example. The size distribution was obtained from experiments with aerosolized manure in a wind tunnel, as described in [10].
In the simulation, the particles were released 50 m downstream the inlet at 1.5 m height. The particles were released from three points at y = −0.25, 0 and 0.25 m. Each second, 400 particles were randomly released from these three points. The particle phase was modelled in one-way coupling since the Stokes number of the particles is less than 1 (2.19 × 10 −6 ) and the volume fraction is of the order of 1 × 10 −15 . The particles were modeled with drag and gravitational force, and without collision between the particles. The particle phase was simulated for a period of 120 s, and the spatial distribution of particle concentration in the domain was estimated afterwards. Particles of diameters ranging from 0.25 to 10 µm, with a total of 24 particle sizes, were analyzed at 1012 positions.
In the x-direction, the particles were analyzed from x = 50 to x = 500 m for every 10 meters; in y-direction, the particles were analyzed from y = −25 m to y = 25 m for every 5 m; in the z-direction, the particles were analyzed at z = 0.5 m and z = 1.5 m. The regions of analysis were spheres of volume 1 m 3 .
The obtained particle concentrations were time-averaged for 120 s and summarized based on the particulate matter size. For example, the particle concentration of all sizes less than 1 µm was summed to obtain the particulate matter concentration of PM1. Similarly, the PM2.5 and PM10 concentrations were obtained and plotted, as shown later in Figure 14.

Results
The results for the particulate matter dispersion from manure application are shown here, with regards to the quality of the simulated wind field, the dispersion of particulate matter, dependent on the diameter and density, and the potential health risks due to air quality. Figure 11a shows the instantaneous velocity field of the domain after a 400 s simulation time (shortly before the particles are injected) with iso-contours of the second invariant of the velocity gradient tensor, which is a useful characteristic to display vortex morphology. Figure 11b shows the same velocity field as a detailed 2D view on the XZ-midplane (y = 0) with a height up to 40 m, and the vorticity displayed. It can be seen that a fully turbulent flowfield has evolved, beginning at approximately x = 50 m. With more distance from the inlet, more vorticity is visible in the higher regions, up to 20 m in the vicinity of the outlet, which is probably due to the additional friction from the floor.  Figure 12a shows the vertical velocity profile in x-direction at three different locations, x = 125, 250 and 375 m, and the values that were initially set at the inlet. In the region from the bottom until a height of approximately 40 m, all simulation results fit very well with the initial inlet data, which means a stable velocity profile is created from the inlet to the outlet. In the region from 40 m height up to the top of the domain, the simulated velocities gradually become lower than the inlet velocity, e.g., at a height of 60 m, the three simulated values are around 10 % slower than the inlet velocity at that height. This is due to the influence of the top boundary and could be improved by increasing the height of the domain.

Atmospheric Boundary Layer
The vertical profile of the Reynolds stress R 11 was also stable throughout the domain with a good match with the experimental results for heights up to approximately 25 m. The good agreement of this turbulent characteristic at the first sampling position at x = 125 m indicates, that vortices are correctly produced by theDFSEM generator and transported without decay from the inlet to the outlet. However, for regions above 25 m, higher deviations from the experimental values occur, with an inverse trend. The reason for this might be the geometry of the computational mesh, which is concentrated towards the ground, resulting in smaller cells in the bottom region, and coarser cells with higher aspect ratios towards the top. The sensitivity of the DFSEM turbulence generator on the cell size and geometry should, therefore, be investigated further.  Figure 13 shows the instantaneous dispersion of particles in the domain, 120 s after their injection, and individually displayed for particle diameters of 1 µm (Figure 13a) and 10 µm (Figure 13b). Both particle sizes follow the flow immediately due to the very low Stokes number. A difference between the dispersion of particles according to their size is visible when accounting for their age (which is the duration of stay of a particle in the domain since its respective injection). It can be seen that old particles (age between 80 and 120 s) of size 10 µm have deposited on the ground in the downstream region, beginning from x = 100 m until the outlet. This deposition did not occur for particles with a 1 µm diameter. Here, the oldest visible particles are at the age of around 60 s, which means that all particles of this size are transported directly through the domain. A similar behavior with a faster deposition of particles with larger diameters was also reported in experimental studies of, e.g., [53].  Considering the concentration of particles at z = 0.5 m and z = 1.5 m in Figure 14, in the lower level, particles are evenly detected after x = 100 m, whereas at the higher level, particles are sparsely distributed. Along the line y = −10 m at the lower level, particles could be located from x = 150 m for PM1 and PM2.5, but, at the higher level, they could only be detected after x = 200 m. This behaviour of the particles was mainly due to the turbulent vortices present in the flow field.

Influence of Size of Particulate Matter
Since the source is located at z = 1.5 m, particulate matter could be seen from the source until the end of domain as in Figure 15a, where the PM10 moves straight for a certain distance and then starts to disperse. At z = 1.5 m level, PM10 was dispersed to a width of y = 40 m. From Figure 15, it could be noted that when the height increases, PM10 was detected further away from the source. At z = 5, 10 and 15 m, the particles could be seen only after 125 m, 250 m and 300 m, respectively.

Influence of Density
To investigate the influence of density, the particulate matter with densities 130 kg/m 3 and 316.6 kg/m 3 were simulated with the same initial conditions and mesh as described before. The time averaged particulate matter concentration was spatially averaged over the y-direction and plotted as in Figure 16. The difference between the PM10 concentrations at z = 0.5 m is small. At x = 70 m, the difference is only 0.7 and similar difference is observed in other points, which are not visible. The concentration of PM10 with 316.6 kg/m 3 is high at z = 0.5 m.
A clear trend for the PM10 concentration at z = 1.5 m can be seen for greater distances from the source than approximately 250 m. Here, the PM10 concentrations from the particles with the lower density are always higher than the concentrations from the particles with higher density. On average, after 250 m, at the height of 1.5 m, the concentration of high-density particles was round 50 % less than the concentration of low-density particles. It could be understood that the particles with a lower density tend to spread further than the particles with a higher density, which tend to deposit to the ground.

Risk Estimation
The simulated concentrations of particulate matter from the manure application example was used to assess the risk of particulate matter exposure for humans dependent on their distance to the pollutant source. This was done based on the Air Quality Index (AQI) by the USEPA [54]. According to USEPA, the risk is assessed such that AQI < 50 is Good, AQI 51 to 100 is Moderate, AQI 101 to 150 is Unhealthy for sensitive groups, AQI 151 to 200 is Unhealthy, AQI 201 to 300 is Very Unhealthy and AQI > 301 is Hazardous.The spatially averaged particulate matter concentration was converted to µg/m 3 by using the density and volume of the particulate matter. For computational cost-effective simulations, the injected particle mass was 0.0066 µg per second. To use the AQI under the assumption of a hazardous environment directly at the source, the simulated concentration results were proportionally increased such that the injected particle mass was 100 µg per second.
Based on the simulation results, the PM10 concentration was classified into four types as in Figure 17. It could be observed from the figure, up to 200 m the AQI is greater than 300, and up to 400 m the AQI is in between 100 and 300 and beyond 400 m the AQI is less than 100.

Flowfield
Regarding the turbulence models, the equilibrium-assumption models kEqn and Smagorinsky were not efficient in reproducing the results, which is clearly due to their inability to handle the transition by these models, as discussed in [36]. However, the solver was able to simulate a very well-matching velocity field in the investigated region with the two non-equilibrium SGS turbulence models WALE and dynamicKEqn. Larger deviations were only seen in the near-wall regions. This may be due to the presence of y + values in the buffer region. Even though nutUSpaldingWallFunction was used to compensate the y + values, a better mesh refinement could have resulted in more accurate results and should be investigated in further simulations. Alternatively, so-called detached eddy simulations (DES, a mixed approach of RANS for simulating the near-wall flow, and LES for simulating the non-wall-affected flow) could improve the wall-near results and simultaneously reduce computational effort, which is described, e.g., in [36].
Nevertheless, the LES results obtained in this study were more accurate than the results from the simulations, which were additionally conducted in the experimental study by Gromke et al. [35], who used an RANS approach. At the points that were used for the mesh convergence study, their published results had an error of round 60%, while the LES dynamicKEqn model showed an error of only 4%. This is a good indicator for the advantages of transient LES over steady-state RANS in the context of complex flows, which will be discussed further in Section 5.4. The dynamicKEqn model was also validated in a previous study to simulate the flow inside a naturally ventilated dairy barn and was found to perform very well for the most important characteristics [55]. Future studies on the coupling of indoor and outdoor flow will, therefore, focus on this turbulence model.

Particle Dispersion
The results for validating the particle dispersion showed no differences between the two investigated solvers for the solid phase. The simulations were obviously less sensitive to the Lagrangian solver for the particle transport, and more attention should be paid to the accurate solution of the flow field, especially the turbulent characteristics. Larger discrepancies between experimental an simulation data for the particle concentrations were found in the lower region of the computational domain, which are most likely a result of the velocity overprediction at the inlet jet, when LES was applied. We tested several different meshes, wall functions and SGS models (for brevity, not shown in this study) with no improvement in the results for the velocity field. Probably, the application of well-defined velocity boundary conditions (in terms of time series) at the inlet would improve the velocity field. Unfortunately, the study for the validation of the particle dispersion did not provide information on the inflow characteristics of the test room. This was carried out under nearly laminar conditions with a Reynolds number of ≈1200.
An experiment in a fully turbulent flow, with documented inlet velocity characteristics, would likely have resulted in a more accurate prediction of the flow field and, therefore, of the particle dispersion. Overall, an acceptable match between the experimental results and the simulations results can be stated, with a good recapturing of the general trend and a good match in the upper third of the domain. However, at this stage of validation, the investigated solvers PimpleLPTFoam and MPPICFoam is rather reliable for an assessment of trends (e.g., general effectiveness of flow barriers) than for the derivation of absolute values as needed, e.g., in civil protection applications. A bottleneck here is the lack of highquality datasets under realistic turbulent conditions, which can be used for validatation.

Application Example
The DFSEM inlet boundary condition was able to create a vertical velocity profile, which matched the experimental data very well. The profile was stable throughout the domain, which is important for the dispersion studies and indicates a fully developed turbulent boundary layer from the beginning. This is usually achieved only after a longer simulated inflow section with natural turbulence generation or a recycling strategy, where the velocity values from the outlet are mapped and superimposed on the inlet [56][57][58]. These strategies require a considerable amount of computational effort. By using the synthetic turbulence generator DFSEM, this effort was not required; hence, this boundary condition is a promising tool for further efficient investigations of atmospheric dispersion processes. Up to a height of 25 m (which covers the usual region of interest when investigating near source immissions) the Reynolds stresses R 11 were in good agreement with the experimentally derived initial conditions throughout the domain. However, for regions above this height, larger deviations were seen, which could lead to inaccurate dispersion simulations. These deviations are likely due to an interaction between the DFSEM and a non-uniform mesh. Hence, the sensitivity of generated eddies on the grid aspect ratio should be further investigated on unstructured meshes. The simulations of the particle phase in the validated flow field showed two main trends for the dispersion of PM. First, particles with larger diameters tend to sediment earlier on the bottom and are not spread as widely as particles with smaller diameters. This was previously reported in experimental studies e.g., in [53,59]. Second, the density of the PM had a strong influence on the dispersion. The concentrations of particles with a higher density after a certain distance were always lower than the ones with a lower density. This can be explained by the interaction between drag and gravitational forces, where the latter are greater for particles with a higher density; hence, these particles tend to sediment earlier, as described by Sehmel [60]. The ability of the investigated solver to reproduce this behavior is an important feature, because it shows the potential for further simulations of varying sources of PM (e.g., from different types of fertilizers) and their respective immission impact in the near field.

Performance: Cost-Benefit Ratio
It is obvious that the transient large-eddy simulations presented in this study require more computing capacities than the usually applied steady-state RANS simulations. However, we believe the inherent advantages of LES compensate for this drawback. In general, for complex turbulent flows with flow separation and transitional effects (as is typical for the flowfield in the near field of animal housing systems), LES are superior to RANS, because the main portions of the turbulence are directly computed and not modeled, as discussed, e.g., in [19,61]. This good applicability for complex flows can also be confirmed by our results for the street canyon flow discussed in Section 5.1. For an efficient use of LES in agricultural applications, which are characterized by usually large dimensions with several scales of interest, the use of massive parallel processing (use of large computer clusters/supercomputers) is indispensable. One or two decades ago, this would make the practical implementation of LES difficult or nearly impossible. However, in recent years, external cloud computing services became more popular, resulting in a decline in costs for simulation time, so that even with a smaller budget, large-scale simulations can be carried out. The use of open-source solvers saves money for the sometimes exorbitant high license fees for proprietary solvers in parallel computing [55], and these saved resources are usually higher than the fees required to rent external simulation time.

Conclusions and Outlook
The goals of this study, which were formulated in Section 1, could be fulfilled. An open source solver for the simulation of PM dispersion was validated with respect to the velocity field, the transport of particles, and the generation of synthetic turbulence, and this was demonstrated with an example application of manure on a field.
It could be shown that a precise analysis of the spatial distribution of PM concentration is possible, with distinguishing between certain particle properties such as size, density or age. This is very useful for the further use of the solver in risk assessment in the near source emmission area. At this stage of validation, the investigated solver can already be used to investigate the impact of measures, such as, e.g., flow barriers, on the particle concentrations in terms of relative changes. To improve the validation process of the particle transport solver, high-resolution measurements of PM dispersion under realistic conditions are needed for validation datasets. Experiments should be designed and carried out to generate these datasets.
However, the results discussed here are valid for dry conditions with a neutral stratification. Particle dispersion in atmospheric flows is usually influenced by humidity effects (coagulation or evaporation), and also by thermal buoyancy in case of non-neutral stratification. Simulations with the tested solver with additional humidity-particle interaction and solving the temperature field will, therefore, improve the results and be the subject of further investigations. Additionally, the presence of pathogens in the PM could be taken into consideration by giving the particles 'pathogen-like' properties, so that the particles are active for a specific period of time and interact with environmental conditions such as UV-Radiation. These extensions can be added to the model by including already existing libraries, such as the OpenFOAM software package.
Given that, the model has the potential for a powerful tool, especially to accurately simulate near-source pollution, e.g., in case of a disease outbreak with pathogen-laden dust in a certain barn, the risk of contamination of the neighboring barns can be assessed, which might lead to less culling of the unaffected animals.

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

Abbreviations
The following abbreviations are used in this manuscript: