A combined CFD-based simulation and experimental study of particle dispersion in a high-concentration airtight space under unorganized airflow

It is significant to obtain the particle distribution data for workpiece processing and industrial protection strategies. To simulate the particle dispersion in a high-concentration airtight space under unorganized airflow, a CFD-based simulation and experimental were utilized in this study. The proposed method works in three steps. The first step is to use Stokes-MRF-Fan-CFD model to simplify the distribution of high-concentration particles calculation to the dispersion state of the gas components. The second step is to compute the airflow patterns and dispersion of particles. The third step is to prove the feasibility of Stokes-MRF-Fan-CFD model. The results of this model in particle dispersion show ±200 mg·m−3 of absolute error, ±20% of relative error, and 88.5747 grade of simulation accuracy. And the results of particle concentration distribution can be mutually confirmed with the streamline and velocity vector distribution in the airtight space, providing a promising and innovative model for the particle dispersion in a high-concentration airtight space under unorganized airflow.


Introduction
With advantages such as high utilization of spraying materials, rapid construction speed and safety, thermal spraying and electrostatic spraying have been widely promoted in various countries in recent years.Spray technology plays an essential role in the engineering field and attracts a large amount of research based on Computational Fluid Dynamics (CFD), including studies of powder transport behavior, [1][2][3] powder concentration distribution 4,5 and spraying process. 6,7It is significant to obtain the particle distribution data for workpiece processing.In addition, the transmission of airborne particles in the work environment threatens human health.Therefore, it is of great importance to know about the dispersion of airborne particles in the work environment for the formulation of industrial protection strategies.CFD has been widely used in simulating airborne particle transport and dispersion. 4,5,8However, the high computing costs of the traditional CFD method restrict its application in simulating the particle distribution of spraying under high-concentration conditions through unorganized airflow disturbances caused by fan rotation.
In recent years, studies have been carried out on the airflow disturbances caused by fan rotation, focusing on the particle-fan interactions [9][10][11] and airflow patterns, [12][13][14][15] etc.However, the airflow induced by fan has different flow patterns under different blowing directions.To date, it is still challenging to capture those complicated airflow fields at room scale.One is that the airflow induced by the ceiling and floor fan with downward and upward blowing directions is nonuniform, producing a strong air movement in its core zone but a weak airflow in the perimeter zone, resulting in unorganized airflow disturbances in the spray environment.Also, the high-concentration could generate a heavy calculation burden of CFD.Reasonable simplification based on Governing equation of fluid flow is needed to reduce the computational complexity of the CFD mode.
This study aims to investigate changes in the unorganized airflow disturbances with fan rotation and simulation the particle distribution of spraying under highconcentration conditions.The airflow in the model spray workshop was established by the spray airflow and recirculated by fans.Theoretical calculations, experiments, and CFD are the main methods used to study the airflow in airtight space.Theoretical calculations are based on the Governing equation of fluid flow and Stokes number, which simplifies the distribution of high-concentration particles calculation to the distribution state of the gas components.The CFD model is established using MRF model, Fan model to speed up the simulation, and validated through field tests.To our best knowledge, this study is the first time to use CFD model to forecast particle concentration distribution of spraying under high-concentration conditions and unorganized airflow disturbances.It is believed that such Stokes-MRF-Fan-CFD model will find its wide application in the fluid mechanics' field.

Numerical methods
In the closed space (spray workshop) with airflow, the main concern of this research is the airflow movement.The air is regarded as a single continuous substance without considering the components of the air and the dispersion of the molecular level of each component.After calculation (see Section 2.4), the maximum velocity of airflow at the spout is not more than 71 m/s while the sound velocity is 340 m/s.The airflow velocity is about 0.2 Ma.In fluid mechanics, the compressibility of the fluid may not be considered when the flow rate is less than 0.3 Ma.Air can be regarded as an incompressible fluid. 16The relationship between air density and temperature conforms to the ideal gas equation of state.And the spray workshop is equipped with constant temperature and humidity equipment, so there is no need to calculate the temperature field.
Flow model.CFD simulation analysis is inseparable from the three conservation equations of fluid mechanics (Governing equation of fluid flow) 17 : Continuity equation (equation ( 1)), Momentum equation (equation ( 2)) and Energy equation (equation ( 7)), which correspond to the three basic physical principles: law of mass conservation, Newton's second law (momentum conservation), and law of energy conservation, as shown in Figure S1 in the Supplemental Material for basic governing equations of fluid mechanics.By solving a series of partial differential equations, such as mass conservation, momentum conservation, energy conservation, chemical component concentration conservation, etc., the analysis of physical phenomena such as fluid flow, heat exchange and material transport can be realized.
In the Governing equation of fluid flow in fluid mechanics, the mass conservation by controlling the volume can be calculated by the Continuity equation 9 : Where, r represents the fluid density and v represents the velocity vector of the fluid, respectively.The time change rate of linear momentum is equal to the resultant force acting on the continuum (in fluid mechanics, the fluid in the spray workshop is regarded as a continuum), and the Momentum equation 12 can be derived: Where, s represents the stress tensor (calculated according to equation Where, t represents the surface stress component, the first letter of its subscript represents the surface on which the stress acts, and the second letter represents the action direction of the stress.For fluids, the stress tensor is usually in the form of the sum of the normal stress and the shear stress: Where, I represents the identity matrix, T represents the viscous stress tensor, respectively.The reason for taking the negative sign of I is that the normal direction of the pressure and the action surface is exactly the opposite. Combining equations ( 2) and ( 4), it can be obtained: The conservation of angular momentum requires that the stress tensor be symmetric (s = s T ): When the first law of thermodynamics is applied to control volume, the expression of the Energy equation 19 is: Where, E represents the total energy per unit mass of fluid, S E is the energy dissipation per unit volume (the energy loss in the process of overcoming the viscous force), respectively.Among the three equations ( 1), (2), and (7), the Momentum equation is a vector equation, the Continuity equation and the Energy equation are scalar equations (including equations in three velocity directions).There are seven equations including the Governing equation of fluid flow and the calculation formula of pressure and energy (internal energy).And the equations include seven unknown quantities including density, velocity components in three velocity directions, internal energy, temperature and pressure.The above equations are closed.In this study, density and temperature are constants, and other unknown quantities can be solved.
In the high-concentration conditions, the trajectory of particle will change after collision.In order to reduce the amount of calculation, the parameters between particles was set to not merge after collision.And the subsequent motion of particles after collision can be calculated according to the Momentum equation.
Drag model.The diameter of the coating material (epoxy resin) is greater than 10 mm and will not float in the air for a long time like aerosol. 20When the epoxy resin hits the surfaces (e.g. the ceiling, floor, and wall) of the workshop, it will adhere to the surfaces.And the epoxy resin will no longer drifting into the workshop with the airflow.In the process of calculating the concentration distribution of particles, it is necessary to consider the sedimentation effect of particles caused by gravity.For the sedimentation of particles, the drag force between particles is a linear function of the drag coefficient.For the interaction between the continuous phase and the dispersed phase, 21 the force acting on the dispersed phase i by the drag of the phase j is: Where, A D represents the linear drag coefficient, v r represents the relative velocity between phase i and phase j, respectively.A D is a linear function of relative speed, and the calculation formula of relative speed is: A D is related to the standard engineering definition of particle drag coefficient (C D ), and the calculation formula is: Where, a cd represents the area density on the interface, the coefficient a cd /4 represents the projection area of the equivalent spherical particle, respectively.In CFD calculation, the correction of non-spherical particle shape can be added to the calculation formula of C D .
Moving reference frames.Eulerian mesh is usually used for hydrodynamics calculation. 22Eulerian mesh can deal with large deformation (because the nodes are not moving), 23 but the ability to deal with node motion is insufficient. 24However, there are many examples of grid boundary motion in reality, such as engine cylinder motion, valve opening and closing, wing motion, etc., and also the rotation of fan blades in the spray workshop in this study.It is necessary to introduce grid motion into CFD calculation to simulate the rotation of fan blades.
The basic physical quantities of CFD calculation are velocity, temperature, pressure and composition, and the displacement of grid nodes is not calculated. 25herefore, in order to make the mesh move, the physical constraint usually imposed on the nodes is speed.The moving speed of the grid is defined relative to the reference coordinate system, which can be static (the laboratory coordinate system, in this study is the reference system of the spray workshop), or rotated and translated relative to the reference coordinate system.In steady-state simulation or transient simulation that does not require an accurate time solution, moving reference frames (MRF) provides a method to model rotation and translation as steady-state problems, 26 while keeping the grid stationary.Figure 1(a) shows a motion reference coordinate system considering rotation and translation at a constant speed.
The velocity of the material point P relative to the motion reference coordinate system (relative velocity) is: Where, v represents the absolute speed in the spray workshop, v MRF,t represents the translation speed of the moving reference coordinate system, that is, the speed of its origin relative to the laboratory coordinate system, v MRF represents the angular velocity of the motion reference coordinate system relative to the laboratory coordinate system, r P , MRF represents the position vector of the material point P relative to the moving reference system, respectively.Considering the complexity of fluid movement, the material point P can be divided into many different coordinate systems at different levels, which is called sub coordinate system embedded into the parent coordinate system.The parent-child relationship of this coordinate system defines a hierarchical or tree coordinate system. 12In the study, the rotation of the motion reference coordinate system relative to the parent reference system can be defined.The parent reference system can rotate or translate, as shown in Figure 1(b).
The relative velocity of the material point P is: Where, v p represents the angular velocity of the parent reference system measured in the laboratory coordinate system, v e represents the angular velocity of the embedded reference frame, r p represents the position vector relative to the parent reference system, r e represents the position vector relative to the embedded reference system, respectively.
Revolution model.The stirring fan in the spray workshop has not only rotation motion but also revolution motion (oscillating).The rotation motion mode and parameters of the fan can be directly added during modeling, but the revolution motion is complex and needs to be defined separately.The three fans in the spray workshop revolve around the rotation axis, it is necessary to establish their local coordinate systems.
For simple CFD simulation applications, the laboratory coordinate system can be used.However, there are rotation and revolution motions of multiple fans in the spray workshop, and the airflow motion is complex.It is necessary to define and use the local coordinate system.The laboratory coordinate system is the benchmark for defining other local coordinate systems, and the laboratory coordinate system belongs to the Cartesian coordinate system.The local coordinate system may be defined as a Cartesian coordinate system, a cylindrical coordinate system, or a spherical coordinate system, and the basic laboratory coordinate system or a previously defined Cartesian local coordinate system may be used as a reference system. 27Three Cartesian local coordinate systems are defined in this research, as shown in Figure 2. The revolution of the fan can be defined by defining the rotation motion with the Cartesian local coordinate system of the three fans as the reference system.Stokes number.The Stokes number (Stk), named after George Gabriel Stokes, is a dimensionless number characterizing the behavior of particles suspended in a fluid flow. 23The Stokes number is defined as the ratio of the characteristic time of a particle (or droplet) to a characteristic time of the flow or of an obstacle, or 28 : Where, t 0 represents the relaxation time of the particle (the time constant in the exponential decay of the particle velocity due to drag), u 0 represents the fluid velocity of the flow well away from the obstacle, l 0 represents the characteristic dimension of the obstacle (typically its diameter), respectively.A particle with a low Stokes number follows fluid streamlines (perfect advection), while a particle with a large Stokes number is dominated by its inertia and continues along its initial trajectory. 29,30n the case of Stokes flow, 31 which is when the particle (or droplet) Reynolds number is less than unity, the particle drag coefficient is inversely proportional to the Reynolds number itself. 32In that case, the characteristic time of the particle can be written as 33 : Where, r p represents the particle density, d p represents the particle diameter, and m g represents the fluid dynamic viscosity, respectively.In experimental fluid dynamics, the Stokes number is a measure of flow tracer fidelity in particle image velocimetry (PIV) experiments where very small particles are entrained in turbulent flows and optically observed to determine the speed and direction of fluid movement (also known as the velocity field of the fluid).For acceptable tracing accuracy, the particle response time should be faster than the smallest time scale of the flow.Smaller Stokes numbers represent better tracing accuracy.For Stk .. 1, particles will detach from a flow especially where the flow decelerates abruptly.For Stk \\ 1, particles follow fluid streamlines closely.If Stk \ 0.1, tracing accuracy errors are below 1%.
Romano et al. 34 studied the advection motion of propylene glycol particles (blue dots) with different particle sizes in the stagnation point air flow field (gray streamline) in the PIV experiment.Figure 3 shows the comparison between two different particle sizes for tracking accuracy for PIV.Simulated particles (blue dots) of propyleneglycol advecting in a stagnation point flow field (gray streamlines).Note the 1 mm particles crash onto the stagnation plate whereas the 0.1 mm particles follow the streamlines.In this research, the particle size of epoxy resin powder was only ;30 mm.The medium is air, the density of propyleneglycol is 1.038 gÁcm 23 , and the density of epoxy resin powder is 1.204 gÁcm 23 .Stk epoxy resin ' 0.104Stk propyleneglycol can be obtained by taking equations ( 13) and ( 14).It can be considered that the epoxy resin powder particles follow the fluid streamline movement in the air.
This study was focused on determining the particle distribution of spraying under high-concentration conditions and unorganized airflow disturbances (i.e. a gas-solid two-phase flow problem).As the volume fraction of particles in this study was small (\\ 1%, only 1/8.7564 3 10 23 , see Section 2.4 for the calculation process), we used a reasonable simplification to solve the problem according to the discussion results in this section.The epoxy resin powder particles follow the airflow, and the mixture of injected air and epoxy resin powder is set as a gas component.This gas component and the original gas in the spray workshop flow together with the airflow, and the distribution of this gas component in the workshop is simulated.Thus, the distribution of epoxy resin powder is evaluated according to the distribution state of the gas components, and the concentration distribution of epoxy resin powder in the spray workshop is numerically calculated.

Spray workshop model
The 3D model of the spray workshop was established after a certain simplification based on the actual spray workshop layout, and includes a spout and three fans.The spray workshop consists of a sealed chamber, three fans and particle concentration monitoring system.Table 1 lists the dimensions of the spray workshop, and cooperating with Figure 2 can help build the geometric impression of the spray workshop.The workshop is 3690 mm 3 1490 mm 3 3000 mm in length (L), width (W), and height (H).The volume of this workshop is 16.49 m 3 .Surfaces (e.g. the ceiling, floor, and wall) of the workshop were made of stainless steel plate.In total, three fans numbered 1-3 were used and located under the ceiling and above the floor.They were placed on the center line (length direction) of the workshop, creating an airflow recirculation zone that was also the spraying zone.In addition, the 16 operation holes are distributed 1.1 m away from the bottom of the spray workshop.
There are eight spray operate channels inside the spray workshop to operate the workpiece.The center line of each channel is 145, 400, 600, 400, 600, 400, 600, 400, 145 mm apart from the walls on both sides and the center line of adjacent channels.The height of the center line of the operate channel is 1100 mm.The outside of the spray workshop is covered with thermal insulation materials.The thickness of both ends is 305 mm.At one end of the length direction of the spray workshop, there is a spout for spraying the coating materials.The diameter of the spout is 12 mm and the height is 1100 mm, which is consistent with the height of the center line of the operate channel.Figure 2 shows the size of the operate channel of the spray workshop and the position of the spout.Figure S2 in the Supplemental Material are topological structure and the details of the operate channel in the fluid domain.
In this study, three ceiling and floor fans with four metal blades (Type: KDS-310) of diameter 31 cm was selected and fixed in the center line (length direction) of the workshop.Geometry dimensions of the fan blade are shown in Figure 4(a).The fan was equipped with an AC/DC brush motor (the rated working voltage is 48 V) and can rotate in different directions at 1450 revolutions per minute (r/min), and the revolution speed was 6.11 r/min.
In order to improve the calculation accuracy and prevent calculation divergence, it is necessary to set a very small time step, which requires a very large amount of calculation.Therefore, the Fan model should be simplified.Before the Fan model is simplified, the flow characteristics of the real fan need to be calculated by MRF model to obtain the relationship between the fan flow and pressure at a specific speed.In this research, the fan is modeled and placed in the fluid channel, and the fluid flow calculation of the real  fan is carried out to obtain the pressure difference between the upstream and downstream of the blade at the speed of 1450 r/min (normal working condition of the fan).The pressure difference between the upstream and downstream of the blade at 1450 r/min speed calculated by CFD is introduced into the Fan model as the momentum source of the fan, and a lumped Fan model that combines the flow and pressure of the fan can be established.The fan is modeled as a volume region, and when the airflow passes through this region, the pressure changes upstream and downstream of the fan will be caused.The fan does not move in the calculation, which can greatly reduce the amount of calculation.In addition, Fan model can also simulate the vortex introduced by the fan and calculate the vortex velocity downstream of the fan.

Mesh generation
A polyhedral mesh was used to discretize the volume of the spray workshop model.The spray workshop contains the fan area whose streamlined direction changes rapidly.In order to ensure the calculation accuracy, the area enveloping the revolution range of the fan shall be subject to mesh encryption, 35 and the fan area and the spray workshop area shall be subject to grid overlap processing through overlapping grids. 36n consideration of the calculation speed of models, the fan area and the spray workshop area are set up.Including the real performance calculation of the fan and the flow calculation of the spray workshop.For the high Reynolds number model (k-Epsilon model) used in this study, the value of y + should meet the requirement of 30 ł y + ł 300, and y + = 30 is taken in general research. 37,38The thickness and number of layers of the boundary layer are 2 mm and 5, the value of y + is 30.The number of fan surface mesh is 1.13 million, and the number of volume mesh is 5 million.The number of spray workshop surface mesh is 600,000, and the number of volume mesh is 3.23 million.Since the airflow nearby fan area is more complex, the grid on fan area is denser than in other areas.The final mesh of spray workshop is shown in Figure 5.

Model validation
To ensure the accuracy of the CFD model, the particle concentration is validated.As there were no similar validation data available for particle distribution under high-concentration, we used the experimental datasets from measured data.The airflow in the spray workshop is complex because it varies with the spout and fans' operating conditions, which is changing during the spraying process.Therefore, the whole period of spraying and sedimentation is selected to conduct the validation.
The testing setups included three phases, for example, spraying, dispersing and sedimentation.The setups were as follows.1-15 s: blowing coating material into spray workshop at 0.7 MPa.1-30 s: using fans to keep coating material dispersing uniformly.31-90 s: let stand to stabilize the concentration of coating material.The average concentration of coating material in the spray workshop is 1000 mgÁm 23 .Table 2 shows the physical parameters of the coating material.
The parameters of the air compressor (volume of the air storage tank V = 160 L, initial pressure p = 0.7 MPa, end pressure p = 0.4 MPa, room temperature T = 294.15K) are brought into the ideal gas equation pV = nRT.The total gas injected from the air storage tank of the air compressor a˜n = 19.64 mol and the volume V = 120 L can be obtained.The volume ratio of the injected particle is 1/8.7564 3 10 23 .The spout radius r = 6 mm, the spraying time t = 15 s, and the injected air velocity v = V/(pr 2 3 t) = 70.74m/ s ' 0.2 \ 0.3 Ma, so air can be regarded as incompressible fluid. 16

Particle concentration monitoring
Nowadays, the standards of particle concentration monitoring have been based on the living environment.One common trait among these standards is that the detection limit is low while the sensibility is high (Table S1 in the Supplemental Material for the standard of particle concentration detection).However, the particle is at high-concentration (.1000 mgÁm 23 ) in spray workshop air.Obtaining this information online has been challenging, because it is easy to exceed the highest detection limit of traditional particle concentration monitors.
Recent advances in electrostatic coupling and digital signal processing made it possible to monitor particles of high-concentration workshop in real-time.A KJ-Z particle concentration monitor (Huarui Environmental Technology Co., LTD) was set up in this study to monitor particle concentration in the spray workshop.The principle of KJ-Z is to monitor the AC-induced electromotive force signals when particle access sensor.The induced electromotive force is proportional to particle concentration in the instrument range.It distinguishes the particles concentration by their induced electromotive force signals to classify the intensity and quantity of the concentration signal.This method has a high upper limit for measuring the concentration of particles, but its time resolution is low.
The air inlet of KJ-Z was situated on the operation hole of the spray workshop, and the direction of the air inlet was perpendicular to spray airflow direction.The distance between the air inlet and the wall was 0.745 m.There was no shielding around the air inlet.

Airflow patterns
Calculation results of fan.The fan performance calculation is steady-state calculation, and the MRF model is used to simulate the fan rotation.After processing the fan result file calculated by MRF model, the pressure difference between the upstream and downstream of the blade is obtained.This pressure difference can be used to measure the characteristics of the momentum source of the fan.When the fan is treated as a Fan model, it is considered that the airflow passing through the Fan model is pressurized by the upstream and downstream pressure difference.According to the calculation, the static pressure at upstream of the fan is 216.01Pa, and the static pressure at downstream of the fan is 0.88 Pa.The pressure rise caused by the fan is 0.882(216.01)= 16.89Pa.According to the performance parameters provided by the manufacturer, the pressure difference between the upstream and downstream of the KDS-310 fan under normal operating conditions is 16 Pa (25°C).The results of the MRF model in pressure rise show 20.89 Pa of absolute error, 25.56% of relative error, which means the simulation results of MRF model are good in general.Figure 6(a) is the static pressure and its contours upstream and downstream for fan. Figure 6(b) is the velocity contours of fan. Figure 6(c) is the streamline of fan.
Calculation results of spray workshop.Figure 7 shows the streamlines distribution and velocity vector field in the spray workshop during the three periods of injection, mixing diffusion and sedimentation.It can be seen that in the injection working condition, due to the high injection speed of the coating material particles, the fan flow line is disturbed to make it deviate from the vertical direction and form three vortices concentrated in the middle of the spray workshop.So that the particles cannot be effectively dispersed in the spray workshop.In the mixing diffusion condition, the vortex formed in the injection condition is broken by the fan airflow.
The newly formed large-scale vortex is located at the upper part of the spray workshop, which can effectively disperse the particles without sedimentation.In the settlement condition, the air velocity gradually decreases (\1 mÁs 21 ), and the vortex gradually stabilizes.The particles begin to settle and accumulate in the corner area away from the spout that is not covered by the vortex.After 60 s, the air velocity is stable at about 0.1 mÁs 21 , and the vortex size becomes larger and gradually covers the whole spray workshop.Part of the particles settled in the corner area begins to float toward the upper part of the spray workshop with the streamline, and the particle concentration in the spray workshop tends to be stable at this time.It can be observed that under the different air distribution modes of the three working conditions, the difference in air distribution in the spray workshop is mainly concentrated in the middle area of the spray workshop, the fan area and the corner area far from the spout.The airflow induced by ceiling and floor fan with downward and upward blowing directions is non-uniform, producing a strong air movement in its core zone but a weak airflow in the perimeter zone, resulting in unorganized airflow disturbances in the spray environment.

Dispersion and fates of particles
Particles dispersion with spray in 1-15 s.In 1-15 s, the fan in the spray workshop is turned on, and the air compressor is used to inject the coating material into the spray workshop.According to the above discussion results, this study simplified the following characteristics of the coating material in the model.The gas injected with the coating material particles is taken as a gas phase.Figure 8 shows the velocity contour and particles mass fraction of the spray workshop in 1-15 s.It can be seen that because the injected air velocity (71 mÁs 21 ) is much higher than the fan wind speed (4 mÁs 21 ), the particles are concentrated in the high-concentration area in the middle of the spray workshop and is not effectively dispersed in the spray workshop.
Particles dispersion with airflow of fan in 16-30 s.In 16-30 s, stop injecting the interfering agent and use a fan to stir to disperse the particles evenly.Figure 9 shows the velocity contour and particles mass fraction of the spray workshop in 16-30 s.It can be seen that after the release of the particles is stopped, under the stirring action of the fan, the particles are effectively dispersed in the spray workshop without sedimentation.
Particles sedimentation in 31-60 s.In 31-90 s, the fan in the spray workshop is turned off and the workpiece operation is started.Among them, the workshop was stand for 31-60 s to stabilize the concentration of particles.At this time, the injection of particles and the airflow disturbance stirred by the fan have stopped, and the particles begin to settle by gravity.Figure 10 shows the velocity contour and particles mass fraction of the spray workshop in 31-60 s.It can be seen that the particles settle in the spray workshop after the fan is turned off, and the gas carrying particles begins to settle and accumulate in the area far from the spout.
Particles sedimentation in 61-90 s.In 61-90 s, the working condition of the spray workshop is consistent with that in the period of 31-60 s.The injection of particles and the airflow disturbance caused by the fan have stopped, and the particles are settled by gravity.Figure 11 shows the velocity contour and particles mass fraction of the spray workshop in 61-90 s.It can be seen that the particles settled in the lower area of the spray workshop 70 s ago.After 70 s, the airflow carrying particles began to float toward the upper part of the spray workshop, indicating that there was still wind speed in the spray workshop.
Analysis of particles dispersion.Figure 12 shows the particle concentration-t simulation curves of the spray workshop.Figure 13 shows the simulation value of particle concentration in operate channels.The spraying, dispersion and sedimentation process of particles can be seen intuitively from the figures.When the particles are injected, the injected particles are concentrated and distributed in the middle of the spray workshop.It is reflected in the sudden jump of  12 and 13, and it is concentrated in the 5-8 channels of the spray workshop.After the injection of the particles, the air carrying particles continues to move forward under the inertial action, resulting in a low-pressure area on the injection path of the particles.At this time, the free air flow disturbed by the stirring fan will fill the lowpressure area, resulting in a decrease of particle concentration.Finally, the injection of particles and the airflow disturbance stopped, and the concentration of particles tended to be stable.In 3-8 channels, the particle concentration is three times of the standard working concentration in the first 20 s.The staff should take respiratory and physical protection and strengthen ventilation.

0-15 s particle concentration in Figures
According to the results in Section 3.2.4,before 70 s, the particles settled in the lower area of the spray workshop.After 70 s, the airflow carrying the particles began to float to the upper part of the spray workshop due to the flow in the spray workshop.The staff should take respiratory and eye protection.Figure 13 shows the standard deviation (s) of particle concentration in seventh and eighth operate channels in 31-60 and 61-90 s.The dispersion degree of the particle concentration simulation value in 61-90 s is lower than 31-60 s, which indicates that the dispersion of particles in the spray workshop is more stable and suitable for workpiece operation.

Quantitive validation
Channel 7 and Channel 8 are the most commonly used positions in spraying work, so we chose these two channels for validation of numerical simulation.Figure 14 shows the comparison of experimental value and simulation value of particle concentration in operate channels.In Figure 14, the experiment value is the real observation particle concentration data, simulation value is simulated by CFD model.The simulation value was worked out and compared with the experiment value to evaluate the model's availability.The simulation results of the model were validated by the field test.It can be recognized that in Figure 14 two curves almost overlap each other.
Figure 15 shows the simulated absolute error rate and relative error rate of particle concentration in operate channels.Analysis of the simulation results error shows that most of the absolute error and relative error are in the interval of 6 200 mgÁm 23 and 6 20%, which means the simulation results are good in general.Several relative errors in Figure 15 are big due to the fact that the simulation value can be affected by many aspects during the study.
The accuracy of CFD model also was computed in this study. 20The formula in accuracy of simulation is as follows: Where, S i represents the simulation value and E i represents the experimental value.problem of heavy calculation burden in CFD forecast particle concentration distribution of spraying under high-concentration conditions and unorganized airflow disturbances is preliminarily solved.

Conclusions
Due to the unorganized airflow and high-concentration of particles in the airtight space, the traditional average concentration calculation method cannot truly represent the concentration distribution of particles.The application of Stokes-MRF-Fan-CFD method is an essential and innovative attempt to simulate the particle concentration distribution in the airtight space.The main conclusions of this study are as follows.
(  dispersion degree of the particle concentration simulation value in 61-90 s is lower than 31-60 s, which indicates that the dispersion of particles in the spray workshop is more stable and suitable for workpiece operation.The results of particle concentration distribution can be mutually confirmed with the streamline and velocity vector distribution in the spray workshop. (3) Application of CFD model displays good results.Most of the absolute error and relative error are in the interval of 6 200 mgÁm 23 and 6 20%, and the accuracy of simulation grade is 88.5747, which statistically indicates that the overall results are preferable with error rate in an acceptable range.
(3)), f b represents the resultant force of force per unit volume (such as gravity and wind force) acting on the continuum, F represents the mass force, that is the force on all fluid particles in the fluid; p represents the hydrostatic pressure; m represents the dynamic viscosity coefficient of the fluid; m# is the second viscosity coefficient of the fluid, respectively.In general gas fluid motion, according to Stokes hypothesis, m = 22m#/3.

Figure 1 .
Figure 1.Diagrammatic sketch of: (a) moving reference frame and (b) superposition moving reference frame.

Figure 2 .
Figure 2. The three-dimensional model of the spray workshop.

Figure 3 .
Figure 3.Comparison between two different particle sizes for streamlines tracking accuracy.34

Figure 4 .
Figure 4. (a) The dimension and (b) the three-dimensional model of the KDS-310 fan.

Figure 4 (
b) is the geometric model established according to the blade data of KDS-310 fan.The blade tilt angle of the fan model is 16.7°, the bearing vibration amplitude is 0.15 mm.When the fan performance curve can be obtained and the exact geometry of the fan has no influence on the calculation results, the Fan model can be simplified by using the simplified Fan model.KDS-310 fan belongs to axial flow fan.The inflow and outflow directions of airflow are the same, that is, flow along the axial direction of the fan.This research uses the Fan Momentum Source model provided by the fluid software to establish the Fan model of the axial flow fan.

Figure 5 .
Figure 5.Comparison between (a) fan mesh refinement and (b and c) surrounding mesh size.(d) The surface mesh and (e) volume mesh of the spray workshop.

Figure 6 .
Figure 6.(a) The static pressure and its contours of upstream and downstream, (b) the velocity contours, and (c) the streamline for fan.
) A reasonable simplification was used to solve the problem according to the Governing equation of fluid flow, Stokes number, MRF model and Fan model.Which simplifies the distribution of high-concentration particles calculation to the distribution state of the gas components.The velocity of airflow carrying particles at the spout is about 0.2 Ma, which can be regarded as an incompressible fluid.It can be considered that the epoxy resin powder particles (Stk epoxy resin \ 0.1) follow the fluid streamline movement in the air, and the mixture of injected air and epoxy resin powder is set as a gas component.(2) The particles are concentrated in the highconcentration area in the middle of the spray workshop and are not effectively dispersed in 1-15 s.Under the stirring action of the fan, the particles are effectively dispersed without sedimentation in 16-30 s.The particles settle in the spray workshop after the fan is turned off, and the gas carrying particles begins to settle and accumulate in the area far from the spout in 31-60 s.The airflow carrying particles began to float toward the upper part of the spray workshop in 61-90 s, indicating that there was still wind speed in the spray workshop.The

Figure 12 .
Figure 12.Particle concentration-t simulation curves of center operate channel in the spray workshop.

Figure 13 .
Figure 13.Simulation value of particle concentration in (a) seventh and (b) eighth operate channels.

Figure 14 .Figure 15 .
Figure 14.Comparison of experimental value and simulation value of particle concentration in (a) seventh and (b) eighth operate channels.
18 s = s x t xy t xz t yx s y t yz t zx t zy s z

Table 1 .
The dimensions of the spray workshop.

Table 2 .
Physical parameters of the coating material.
The grade of CFD model is 88.5747 on the accuracy of simulation results.The grade shows that this CFD model is efficient to simu- late particle concentration in a high-concentration spray workshop under unorganized airflow.The absolute error, relative error and the grade of simulation accuracy indicate that the overall results are preferable with error rate within an acceptable range, which shows this CFD model is a desired model for this study.The