Deposition of dry particles on a fin-and-tube heat exchanger by a coupled soft-sphere DEM and CFD

In this study, a novel computational model is utilized for investigating fouling of two commonly encountered heat exchanger fin shapes in an air-conditioning application. The computational method utilizes the discrete element method (DEM) coupled with a largeeddy simulation (LES) framework. The fin-and-tube heat exchangers (FTHE) are investigated for three different Reynolds numbers (ReDh =243, 528, 793), three different particle sizes (Dp = 5, 10, 20 μm) and two different adhesive particle types based on the experimental values in the literature. The code is first benchmarked from the CFD and DEM viewpoints. A comprehensive fouling study of the FTHE’s, consisting of altogether 36 simulations, is then carried out. The major numerical findings of the paper consist of the following four features. First, with low adhesive particles, the plain fin shape has a 3.45 higher volume fouling rate with ReDh =793 than at ReDh =264. With the herringbone fin shape, and the low adhesive particles, the volume fouling rate is 1.76 higher with ReDh =793 than at ReDh =264. Second, for the high adhesive particles, the plain fin has a 5.4 times higher volume fouling rate at ReDh =793 than for ReDh =264. The herringbone fin shape has a 3.92 times higher volume fouling rate with the highest Reynolds number of ReDh =793 compared to ReDh =264. Third, high adhesive particles have 3.0 times higher volume fouling rate than low adhesive particles for both fin shapes, all particle sizes and all Reynolds numbers combined. And finally, herringbone fins have 1.74 times higher volume fouling rate than plain fins for low adhesive particles. For high adhesive particles, herringbone has 1.8 times


Introduction
Globally around 10 − 20% of the energy consumption in developed countries is used by the Heating Ventilation and Air Conditioning (HVAC) systems of buildings [1]. Almost half of the total energy consumption in buildings is due to air conditioning [1]. Commonly used type of a heat exchanger in the air conditioning unit is a Fin-and-Tube Heat Exchanger (FTHE). Inside the FTHE, heat is exchanged between the flowing air between the fins and the fluid flow in the tubes. The efficiency of a heat exchanger is defined as the ratio between the exchanged heat and the induced pressure drop [2]. Considering the complete heat transfer process, around 60 − 80% of the resistance has shown to be on the air side [3,4]. The most traditional design of a FTHE with plain fins [5,6] can be made smaller by making the fin wavy [7,8,9]. In this way, smaller heat exchangers could be designed. Later, as enabled by louvered [10] and slit [11] fins even more compact heat exchangers have been made for various different air conditioning applications. Quite recently, vortex generators [12,13] including winglet type flow actuators [14] have been proposed for enhancing heat transfer. Most of the previous CFD studies are conducted with ideal air without impurities. Therefore, new designs obtained with CFD are only valid for a new heat exchanger for a short period of time.
After long term operation, air impurities begin to accumulate inside the heat exchanger [15,16,17], commonly denoted fouling. Such fouling processes can be classified into wet and dry particulate fouling. The present work focuses on the dry particulate fouling. Dry particulate fouling process can be further divided into the nucleate and bulk fouling regimes [18]. Every fouling process will initially start with the nucleation regime where particles accumulate on the surface and form a deposition distribution. Later on, once the surfaces have been initially deposited by particles, the nucleation regime transitions typically towards the bulk fouling regime. At this stage the induced pressure drop from fouling is at its maximum in a specific operation point.
In HVAC applications, typical substances that contribute to fouling are dust and fibres that can originate from clothes, furniture and fur of domestic animals [19]. Fouling in dry FTHE environments occurs mainly in the flow stagnation and recirculation regions [20].
Li et al. [21] propose that both particle size, shape and properties are important factors in the fouling process. While fouling of larger particles is dominated by inertial impact, smaller particles deposit by eddy transport and thermophoresis [22]. Due to the rather low temperature gradients in HVAC systems, thermophoresis is often negligible as one of the fouling mechanisms [23].
The challenges of studying FTHE by analytical and empirical means was pointed out by Siegel and Nazaroff [24]. They present a deterministic model for FTHE fouling with particle impaction, gravitational settling and Brownian diffusion. The authors pointed out deficiencies in the model for particles in the range of 5-20 micrometer scale particles having rather large initial velocity [24]. Another computational model by Inamdar et al. [25] was compared against the experimental results. The model was shown to predict the fouling trends in various different heat exchangers to an acceptable precision along with a qualitative estimation of the deposition distribution.
Experimental application field studies have been conducted to study the effect of the real fouling environment and the effect of fouling on the thermal hydraulic efficiency of the FTHE. Ahn et al. [19] showed that, in residential buildings, the mean fouling particle diameters typically range between 1-20 micrometers. In their study, they concluded that after 7 years, in both dry and wet cycle heat exchanger, the fouling can cause the pressure drop over the heat exchanger to increase up to 44% while the cooling capacity decreases 10 − 15% [19]. Another study by Park et al. [26] shows that dry particle fouling causes 4−12% decrease in heat transfer coefficient whereas, the pressure drop increases by 22−37%.
Above, FTHE's in HVAC context was discussed. Another branch of FTHE fouling investigations is related to heat recovery boilers. In such applications, the heat exchangers are exposed to ash and particulates originating from combustion processes. [27,28,29,30].
The properties of the fouling particles in combustion fumes obviously differ drastically from the normal indoor fouling particle properties. The effect of fouling have been experimentally studied in the laboratory environment with various different standard dust types by various authors. An experimental study by Zhang et al. [31] for air conditioning application, con-4 cluded that interrupted louvre fin shapes are more prone to fouling in a microchannel design than a planar fin. One of their main conclusions was that the interior of the heat exchanger with louvre fin is more vulnerable to fouling, whereas with the plain fin, the fouling occurs in the front surfaces. They emphasized that the type of the fin plays a prominent role in the fouling characteristics of the heat exchanger. [31] The microchannel heat exchanger and dry coolers were studied by Bell et al. [32,33], who concluded an increase of up to 50% due to fouling on the pressure drop on the air side. Zhan et al. [34] studied the wavy fin-and-tube heat exchanger and concluded that areas that are most prone to fouling are the leading edge of the fins and the front part of the tubes. A very recent experimental study by Zhan et al. [18] investigated the fouling mechanism of wet-particle deposition and concluded that louvered fins are the most conducive to wet-particle deposition. They also pointed out that dry-particle and wet-particle fouling share several similarities. Yet, the fouling rate in dehumidifying conditions can be a factor of 5 to 7 times higher than the dry particle environment.
The numerical fouling studies conducted on an air conditioning application are very few.
Zhan et al. [35] studied the bulk regime fouling with 50 µm particles with Reynolds-Averaged Navier-Stokes flow simulation and developed a model for the fouling thickness. The effect of particles sizes on the deposition areas on a metal foam heat exchanger has been studied by Sauret and Hooman [36] in laminar conditions. Kuruneru et al. [37] has developed and validated a CFD-DEM model to study fouling of a metal foam heat exhchanger. Kuruneru et al. [38] studied the fouling bulk regime of 25 µm particles in an oscillating unsteady laminar flow and made a comparison study between sawdust and sandstone type particles [39] in the metal foam heat exchanger. Kuruneru et al. [40] also studied fouling with 50 µm and with very low inlet velocity of around 0.1m/s, which is very rarely encountered in HVAC application. All these studies focus on the bulk regime fouling for a very simple tube bank geometry with surrounding metal foam porous area.
Based on the literature survey, fouling of the heat transfer surfaces is an important concern in heat exchanger design. There is a clear research gap for better understanding FTHE fouling processes by simulations which, not only model the fluid flow by CFD, but 5 also model the particle phase and the particle-surface interactions. The main objectives of the present study are listed as follows. First, a CFD-DEM model is benchmarked to model fouling processes in heat exchangers. Second, two different fin shapes with three different particle sizes, 3 typical inlet velocities and two different particle adhesive types are simulated to illustrate and to quantify, via a pioneering research approach, the fouling rate in the nucleation regime of the heat transfer surfaces. The first fin shape is a plain fin without fin forming. The other, very common fin shape used in the industry, is called the herringbone. The overall motivation is to find quantitative results, which can be used to asses the amount of fouling between the different fin shapes. Thereby, in larger scope, the study demonstrates a next level approach to develop enhanced heat transfer surfaces with lower pressure drop and higher heat transfer rates. 6

Details on the numerical domain and setup
The geometry in the present study is widely used in HVAC industry and has been studied by others as well Chokeman and Wongwises [9], Pirompugd et al. [41]. Two fin types are studied to illustrate the geometry effects on the fouling characteristics.

Computational domain
The computational domain is limited to one channel between the fins. The simulation domain is shown in Fig. 1 for the herringbone fins. For the plain fins, the geometry is otherwise exactly the same. The inlet is placed a distance of P l upstream from the heat exchanger geometry to make sure the deposition process is independent of the inlet position. Similarly, the outlet is placed 2P l downstream of the heat exchanger geometry. Further details on the geometries are provided in Table 1.
where the flow is in a channel with the length L = 2P l = 44 mm, height F h = F p − t = 1.695 mm and width P t = 25.4 mm. In the normal direction of the fin, a staggered tube array is penetrated through the fin pack with the diameter of D c = 9.76 mm, longitudinal tube pitch P l = 22 mm and transverse tube pitch P t = 25.4 mm. In this type of a flow, the flow length scale can be calculated as shown in equation [1]. The geometry details are shown in Table 1. The reference velocity used in this study is the average flow velocity at the minimum cross sectional area calculated as U max,avg = U inlet

A inlet
A min , where U max,avg is the maximum average flow velocity, U inlet is the inlet flow velocity, A inlet is the area of the inlet boundary and A min is the area of the minimum cross section inside the heat exchanger.
In the present study, we use the no-slip boundary condition for velocity at the solid boundaries. Respectively, for pressure the zero gradient boundary condition is utilized at 8 the walls. At the inlet of the domain, constant velocity is prescribed while the zero gradient boundary condition is used at the outlet. At the lateral boundaries, periodic boundary conditions are used. The pressure is fixed at the outlet and assumed zero gradient at the inlet.

Governing equations for fluid flow
The governing equations for three-dimensional, incompressible, transient viscous flow with two-way coupling between the particle phase and the fluid through the heat exchanger are the continuity (2) and the momentum equations (3).
where α f is the local fluid volume fraction, R pf = K pf (u − v ) is the momentum exchanger with the particle phase. For the momentum exchange between the domains an implicit second-order accurate Crank-Nicolson scheme is used. The source of the momentum equation is the sum of the local particle-fluid forces (F i,fluid ) In practice, the studied particle volume concentrations are in the order of φ = 5 × 10 −5 , which corresponds to the 1-way coupling regime [42]. R pf is small compared to the other terms and is neglected in the present study resulting in a 1-way coupling.
The governing equations are discretized by using a finite volume method and spatial terms are discretized by a second order accurate discretization scheme. The PISO (Pressure-Implicit with Splitting of Operators) [43] algorithm is used to couple the pressure and velocity fields. The CFD-DEM simulations are carried out with the open-source CFD-DEM solver called CFDEM (version 3.8.0) [44], which combines the CFD (version 5.0) toolbox OpenFOAM [45] with the DEM solver LIGGGHTS [46].
For accurate prediction of the trajectory calculation for the deposited particles, it is important to solve the flow field and the turbulence quantities correctly [20]. Because the particles are influenced by the largest, energy containing turbulent eddies throughout the flow, the larger eddies are resolved but the smaller ones are modelled. Such an approach is called Large Eddy Simulation (LES), where only the smallest scales of the turbulence structures are modelled while the larger scales are resolved directly. The sub-grid scale viscosity ν sgs in this study is modelled by using a Wall-Adapting Local Eddy-viscosity (WALE) model with default parameters by Nicoud and Ducros [47], which is suited for wall-bounded flows as the eddy viscosity naturally goes to zero at the walls Mirzaei et al. [48].
As Nagaosa [49] concluded in the turbulence-model free study of a plain FTHE, the flow field inside the FTHE is fully laminar for Re D h = 400 (converted from Nagaosa [49] using the reference velocity and length scale as defined in the present study). A transition regime with both laminar and turbulent regions was reported for Re D h = 400-2000 and fully turbulent flow was observed for Re D h = 2400-3200. As a remark, in the present work Re D h = 243-793 which may involve both laminar and turbulent features. For the herringbone fin shape, the wavy fin shape induces a new lateral deviation to the flow field. Therefore, it will shift the spatial location of the transition as well as the critical Reynolds number for the transition.

Governing equations for particles
The collisions between the particles with each other and the fin surface is modelled with the soft-sphere discrete element method (DEM) approach [50]. If a particle i with the mass m i and radius r i , then the mass momentum of inertia can be calculated as I i = (2/5)m i r 2 i . The governing equation for the location x i is given by 5 : where the F con is a contact force upon collision and F fluid is the combined fluid force acting on the particle.

Fluid forces on particles
A Lagrangian approach is used to track the particles as they flow through the heat exchanger. In this study, particles with diameters of d p = 5 µm, d p = 10 µm and d p = 20 µm and with density of ρ p = 2500 kg m −3 are being considered.
Here, we use the particle drag formulation by Benyahia et al. [51], which is based on the simulations by Hill et al. [52] and Koch and Hill [53], where the modified Stokes drag is defined as In the definition, Re p corresponds to the particle Reynolds number while α f corresponds to the particle volume fraction. By using this definition, a larger range of Reynolds numbers and particle volume fractions are covered.

Contact forces
The adhesive force between two spherical particles was originally studied by H.C. Hamaker in 1937 [54]. Hamaker concluded that the dominant forces of adhesion for two materials are the van der Waals and electrostatic forces. They originate from the continuous change of the electrical potential of atoms as the electrons circle around the core. The model that describes these forces is the Johnson-Kendall-Roberts (JKR) model that was originally developed by Johnson et al. [55]. The JKR model is suitable for the specific type of collision for which the Tabor parameter λ T = (4Rγ 2 /E 2 D 3 min ) > 3 [56], where γ is the surface energy density, which is defined as half of the energy required to separate two particles in contact and D min is the minimum separation distance, usually assumed to be 1.65Å [57,58].
Since the surface energy density is defined for a specific material to interact with itself, it is important to notice that the value cannot be used directly in the collision computations between two different materials. For this a new property called adhesion work w = √ γ 1 γ 2 [59] is used.
The contact forces in the normal direction of the surface that are modelled with the JKR model are the spring force F spring,n , the adhesive force F jkr,n and the damping force F damp,n : where n is the surface normal vector and a is the contact area.
F jkr,n = 4 πγEa 3 n The effective Young's modulus is defined as 1 j E j and the effective radius 1 R = 1 r i + 1 r j for particle collision between two materials where E and ν are the Young's modulus and Poisson's ratio and the subscript corresponds to the colliding materials i and j.
In order to model the dissipation of kinetic energy upon collision, a damping force F damp,n is used: where v n is the relative normal velocity, β is a parameter that takes into account the coefficient of restitution e as: S n is parameter that takes into account the material properties as: where δ n is the the overlap distance. In the tangential direction of the contact, the spring force F spring,t is used: where S t = 8G Rδ st is a parameter for the particle properties and the δ st is the tangential overlap. The effective shear modulus is calculated as 1 As was done for the normal direction, a similar damping force F damp,t for the tangential direction is used: where the 1/m = 1/m i + 1/m j is the effective mass and v t is the tangential velocity respect to the surface.

Computational model validation
Next, model validation results are discussed. First, a mesh sensitivity analysis on the pressure drop over both fin shapes is presented. Second, the collection efficiency of the FTHE for low adhesive particles at relevant particle sizes and two fin shapes is considered for different mesh resolutions. Additionally, the overall fluid dynamics and turbulence modelling is validated for a three-dimensional cylinder in crossflow. Third, the flow field inside the FTHE is validated against experimental results available in the literature. Fourth, the contact mechanics model is validated. Finally, the drag experienced by the particles is compared with the analytical Stokes equation and a particle number sensitivity study is carried out.

Mesh sensitivity assessment for fin-and-tube heat exchanger
In LES, it is of high priority to assess the mesh sensitivity of the results. From the viewpoint of particle transport, the energy containing flow scales should be resolved in order to capture the deposition process reliably [60,61]. In the present study, we aim at resolving the flow field well, in order to capture the particle dispersion adequately. Thereby the subgrid scale effect on particle dispersion is neglected.  Next, the collection efficiency for low adhesive particles for the highest Reynolds number is shown for both fin shapes and three different particle sizes in Fig. 3. The same mesh resolution is used for the plain and herringbone fin shape, which was illustrated in Fig. 2.
More information about the adhesiveness of the particle types is provided in Table. 2 15 0 1 · 10 6 2 · 10 6 3 · 10 6 4 ·  Based on the mesh sensitivity study, the results with the 2 million cell mesh in Fig. 3 deviate only 3% on average from the results with the 4 million cell mesh for the highest Reynolds number and low-adhesive particles. Such a scenario can be considered as the most conservative simulation case. Therefore, we conclude that the scales resolved by the 2 million cell mesh are sufficient to capture the particle transport in the present Reynolds number range. The low-adhesive particles are expected to bounce multiple times more from the surface in contrast to the high-adhesive ones and therefore the low-adhesive particles are influenced by the turbulent scales for a relatively long period of time. Thereby, we assume that the mesh resolution is sufficient for high-adhesive particles as well.

Single cylinder flow validation
In the previous section, the mesh with around 2 million cells was noted to provide acceptable accuracy. Next, using the same spatial resolution, the flow around a single cylinder is confirmed for the range of Reynolds numbers. The boundary layer flow around a cylinder is validated by comparing the drag coefficient and Strouhal number of a cylinder in a cross flow to the numerical results found from the literature for CFD studies [62,63,64] and experimental results [65]. The same y + -value and maximum cell size was chosen as was selected to be sufficient in the mesh sensitivity study in section 3.  [65] CFD in this study   [65] CFD in this study

Flow field validation inside a fin-and-tube heat exchanger
Model validation for the flow field inside the FTHE is carried out in Fig. 6  Experimental herringbone by Chokeman and Wongwises [9] CFD herringbone this study Experimental plain correlation by Wang et al. [5] CFD plain this study The agreement with the experiments Wang et al. [5] and the simulated plain fin is good.
Slightly more deviation can be seen for the herringbone fin shape [9]. The difference can be caused by the non-matching boundary conditions such as the turbulence level of the inlet flow. However, based on the numerical results presented in Fig. 6, we conclude the present model to be quantitatively reliable.

Particle number sensitivity study
A simple comparison was first performed between the analytical equations and the CFD-DEM framework to validate the drag force for individual particles. Settling velocity for all particle sizes and the stopping distance for all particle sizes with the corresponding initial velocities used in present study was compared to the analytical values derived from the Stokes drag law. A mean error of 1% was seen between the analytical and the value calculated with the CFD-DEM model used in present study. To

Contact mechanics validation
As a last demonstration of the model functionality, the contact mechanics between a particle and surface is validated by investigating the bouncing motion of a 6 mm Teflon particle impacting a soda glass surface. The material properties used for the comparison are tabulated in Table 3 in the Appendix.  The height of each individual bounce was seen to be almost identical with the reference values found from the literature [38,66]. Thereby, we conclude that the contact mechanics for the particle impact with a wall surface is correctly implemented.

Particle property selection based on critical velocity measurements
In the context of particle-surface impact, the term 'critical velocity' refers to the incident velocity threshold under which a particle will stick to the surface. In practice, critical velocity depends on the particle properties (size, density, Young's modulus, Poisson's ratio, surface energy density and shape) as well as surface properties (Young's modulus, Poisson's ratio, surface energy density and surface roughness) [59]. Herein, a major effort is carried out to collect literature data on critical velocities in order to deduce and justify a meaningful parameter range for the simulations.

Description of the simulated properties
The particle concentration, size distribution and HVAC working conditions of the FTHE can be arbitrary. The fouling process is very slow and can be measured in months or years.
A fouling process consists of a chain of events, where first, the smaller particles with higher critical velocity will start to deposit on the surface. For a clean surface, the critical velocity of larger particles is too low and therefore they will just rebound on the surface and re-entrain to the flow. This means that in the start of the nucleation regime, the smaller particles start to deposit on the surface and change the adhesion mechanics of the surface for the larger particles. This increases the critical velocity and therefore increases the amount of deposited larger particles. Because of the complexity of the fouling process and the variety in particle properties in the real application, it is essential to simplify the process of selecting the material properties so that the results are repeatable and comparable.
The material properties selection process is based on the critical velocity of a specific particle size. The information found in the literature on different particle and surface pairs and their critical velocity is shown in Fig. 9. Based on this figure, the critical velocities are typically in the range of U crit ∈ [0.1, 5] for D p ∈ [1, 20] µm. A relation between the particle diameter and the critical velocity can be observed from this data. Therefore, even with no information on the material of the particles or the deposition surface, we can assume that by using values that corresponds to the minimum and maximum of the data set, the results will represent both extremes of the particle types encountered in the real application. When the 22 contact mechanics of the chosen properties for a specific particle size are investigated further in Fig. 10, we can observe the restitution curve for each particle-surface pair, here for the D p = 20 µm. Finally, in Fig. 11, the impact kinetics of the collision with the aforementioned properties is illustrated for both adhesive levels with the same incident velocity. All the chosen pairs that lead to the size specific critical velocity are tabulated in the Appendix in Table 4.   [72] Uranine/Al and brass [73] Simulated particles: Experimental results from literate:

Non-dimensional groups
A fouling process of a FTHE can be characterized by six parameters, namely the flow Reynolds number Re D h , the particle Stokes number St τe , the adhesion parameter Ad, the elasticity parameter λ, the density ratio χ and the particle to hydraulic diameter ratio further discussed in what follows.
The Reynolds number is defined as Re D h = U max D h /ν. Note that the used definition is based on the hydraulic diameter calculated in equation (1), the characteristic velocity at the minimum cross sectional area U max,avg and the dynamic viscosity ν = 15 µm 2 / sec. The Stokes number St τe = τ e /τ f . The parameters that define the adhesiveness of the particles are the adhesion coefficient and elasticity parameter that are defined as Ad = γ U 2 max ρp dp and λ = E eff U 2 max ρp [59]. Finally, the density ratio χ = ρ f /ρ p = 0.0005 and the particle diameter = d p /D h .
Although the fouling process is a slow process, the computational resources allow us to perform LES simulations of few-tenth of a second. For this reason, the volume fraction φ =V p /(V p +V f ) of the particles at the inlet is increased to a higher value and kept constant between the different cases. The concentration value is kept low enough so that its effect on the flow field is negligible. The volume fraction being φ = 5 · 10 −5 for all the simulations. As these values are smaller than 10 −4 , their effect on the flow field can be neglected [74]. The mass coupling parameter for the simulations is kept under 0.1 for all simulations to ensure that the mass coupling effects are unimportant [75]. The coefficient of restitution e = 0.5 is kept constant for all the particle types. An overview of the simulations carried out in this study is provided in Table 2.

Results and Discussion
Then the collection efficiency is shown with respect to the Reynolds number as well as Stokes number and adhesion parameter. Finally, an example is provided on the practical relevance and applicability of the present results. It is noted how the laminar inflow impinges on the tube surface (I). Based on the literature, this point corresponds to a major region of particle deposition in various FTHE applications [34,35]. Behind the tubes, a recirculation region (II) is formed. Such regions experience typically rather poor heat transfer characteristics [12]. Then, the flow undergoes transition from laminar to turbulent (III) while the flow is qualitatively relatively turbulent with incoherent features close to the outlet boundary (IV). We note that the channel height and cylinder diameter based Reynolds numbers are respectively 312 and 1800. Hence, the flow is not fully turbulent from the viewpoint of standard channel flow but, instead, the unsteady cylinder wake provokes dynamic and incoherent flow features. Please see Nagaosa [49] for further information.

Fouling locations
The deposition locations for plain fins are shown from the upstream side in Higher initial momentum will lead to larger amount of deposited particles on the fin and lower amount on the tube surface (II and III). In contrast, the high adhesive particles will most likely stick at the first collision, excluding the ones that have accelerated with the flow to achieve an impact velocity high enough to enable the particles to bounce, re-entrain to the flow and even hit the fin surface (IV). In fact, we have noted that the qualitatively, very similar trends on particle deposition would be noted for 5 µm (0.16 ≤ St ≤ 0.47) and 20 µm For all parameter values the qualitative trends are rather similar: most particles will deposit on the front side of the tubes. The deposition location results illustrate that with the ID =13, with the lowest Reynolds number and low adhesiveness, the particles at the center of the flow field will have just enough momentum so that they will bounce once and only deposit if they will hit the surface multiple times. Therefore only few particles will deposit in the middle of the tube. With the higher Reynolds number (ID=14 and ID=15) the particle inertia increases to a level that is enough for the particles to escape from the surface after the collision and therefore the tube surface will have less and less deposited particles. When the ID=15 is compared to ID=14, it can be seen that when the Reynolds number is even higher, the particles will start to deposit on the fin after the collision with the tube surface. With the higher adhesion levels (ID=16-18), the particles will deposit on the tube when they deviate from the streamlines and only rarely bounce back and deposit on the fin. For ID=18, an empty region can be seen on the second tube row where less particles have deposited on the middle of the tube. This is due to the fact that as the flow speed is increased in the middle of the channel between the fins, the particles will have just high enough impact velocity so that they will bounce back, even though there exists a strong adhesive force between the particles and the fin. The low adhesive case with the highest Reynolds number, ID=21, is the only case where no particles are not deposited on the downwind side of the second tube row.
For the herringbone fins, the deposition locations are demonstrated for both upper and lower fins separately. The findings for the herringbone case are consistent with the plain fin cases above. Due to its wavy shape, the number of deposited particles is observed to be higher for the herringbone fin cases. As the particles enter the heat exchanger, the flow is initially guided towards the upper fin while most of the particles first hit the first wave of the lower fin (I). As for the ID=19 and ID=22, the deposition locations are almost identical between the cases. This means that in these cases, when the particle hits the surface, the impact velocity is under the corresponding critical velocity. Therefore, the adhesiveness of the particle does

Collection efficiencies
The fin collection efficiency can be defined as the ratio of the deposited and the inserted particles. It is clearly noted from Figs. 15-16 that, for high adhesive particles the collection efficiency will increase along with Reynolds number. The main reason for this phenomenon is the Stokes number, which increases along with the Reynolds number: the higher the Stokes number the more likely a particle deviates from the streamline. For low adhesive particles, the deposition process is much more complicated since the deposition rarely takes low-adhesive herringbone high-adhesive herringbone low-adhesive plain high-adhesive plain low-adhesive herringbone high-adhesive herringbone low-adhesive plain high-adhesive plain low-adhesive herringbone high-adhesive herringbone low-adhesive plain high-adhesive plain In Fig. 18, the collection efficiency is shown with respect to the Stokes number. A conclusion can be made that a higher Stokes number does not always correlate to a higher collection efficiency and therefore a higher rate of fouling. For high adhesive particles the collection efficiency is shown to be a function of the Stokes number. But when investigated more carefully, different particle size classes will lead to a small difference in the adhesive properties and therefore to a discontinuation in the collection efficiency. Thus, in fouling investigations, one should not only consider St and Re but also take into account the particle adhesiveness (Ad) and the elasticity parameter (λ) which further depend on the adhesion work, Young's modulus, particle size and the density of the particle as seen in Eqs. 6-12.
Finally on the Fig. 19, the collection efficiency is illustrated with respect to the adhesion parameter. It is clearly seen that a higher adhesion parameter does not correlate directly with a higher collection efficiency but in fact the opposite seems to be true. This is because the adhesion parameter is the measure between adhesive force and the particle inertia [59].
In other words, the parameter acts as a relative measure of the impact process between the particle and the fin. The fouling process of a FTHE is not only a function of the adhesion parameter but also a function of the Stokes number, Reynolds number and the elasticity parameter as shown in this study.

Fouling rate prediction from air quality measurements
In the earlier sections, the presented numerical results indicated that the fouling process of a FTHE depends monotonically on the Reynolds and Stokes numbers for high-adhesive particles. In contrast, for the low adhesive particles, a non-monotonic relationship was observed. As expected, it was noted that the adhesion and elasticity parameters are equally important in the fouling process. Next, the results will be applied to estimate the fouling rate in a certain FTHE application. The air quality and the particle size distribution of the ambient environment can vary drastically between different applications. For example, the outdoor air during a sandstorm most likely has a much higher concentration of low adhesive particles compared to the indoor air of an average school with high adhesive clothing fibers and other organic substances that are circulated through the heat exchanger. Therefore, when a fouling rate comparison between different fin shapes is performed, it is important to define the particle size distribution and the adhesion properties of the particles in the air.
Next, the present numerical results are applied in the context of the data provided by Cheng and Lin [76]. The data includes the average particle mass size distribution of the air at the Taipei  The thickness of the fouling layer decreases the efficiency of the FTHE in two different ways. First, the fouling layer acts as an insulation in the heat transfer process between the air and the fluid in the tubes. Second, the thickness of the fouling layer will decrease the minimum cross sectional area of the flow and therefore increase the pressure drop over the FTHE. By knowing the volume of each particle size (V dp ), the volume fouling rate can be calculated as F volume = V dp F number . The volume of the deposited particles in unit time for a frontal area of the FTHE can be seen in Fig. 21.  First, it is noted that low adhesive particles have a lower volume fouling rate than the high adhesive particles. Second, for a given adhesive type, the plain fin has a lower volume fouling rate than the herringbone fin. It is clear that the volume fouling rate increases with the Reynolds number. If the Reynolds increases by factor 3, the fouling rate is not desired to increase more than by factor 3. The volume fouling rate for low adhesive particles is noted to be 3.45 times higher for Re D h =793 than at Re D h =264 in the case of the plain fin shape.
In contrast, the herringbone fin shape at Re D h =793, has a 1. 76  Last, we note that in practice FTHEs may need to operate in different ambient conditions. Therefore, based on the numerical findings, the type of particles can affect the choice of the optimal fin type. In fact, such design guidelines are already followed in the industry by practical experiences without detailed knowledge on the phenomena. By summing over the fin types and the Reynolds numbers, the fouling volume rate was found to be approximately 3 times higher for the high adhesive particles than for the low adhesive ones.
As discussed earlier, enhanced fin shapes such as herringbone are used to increase the heat transfer on the air side. In an environment with low adhesive particles the volume fouling rate in the nucleation regime is 1.74 times higher for the herringbone fin shape when compared to the plain fin. For the high adhesive particles, the volume fouling rate is 1.8 times higher with herringbone. Together, when both adhesion levels are summed together, the volume fouling rate for the herringbone fin shape is 1.78 times higher when compared to the plain fin FTHE. The presented numerical results indicate the relevance of combined CFD-DEM studies regarding the fouling rate and the performance of the fin during the lifetime of the FTHE.

Conclusions
In this study, the deposition of dry particles on a fin-and-tube heat exchanger by a coupled soft-sphere DEM and CFD was performed. The novel selection of material properties for the fouling particles is done such that a representative range of Reynolds numbers, Stokes numbers, elasticity and adhesion parameters is covered. The material selection in this study will lead to critical velocities between the particles and the fin surface that corresponds to the measurements performed for various material combinations and therefore different fouling characteristics. All the different models used for the calculation of the flow field and particle drag were validated. A comparison was carried out between low adhesive and high adhesive particle environments, with particle sizes of D p = 5, 10, 20µm, with three typical FTHE Reynolds numbers Re D h =243, 528, 793 and two different fin shapes found in the contemporary HVAC industry. The major findings of this study are summarized below.
1. Novel method for the selection of adhesion properties such as particle size, particle Re D h =793 than at Re D h =264. The herringbone fin shape has a volume fouling rate of 3.92 times higher with Re D h =793 than at Re D h =264. 39 5. High adhesive particles will have 3.0 times higher volume fouling rate than low adhesive particles for both fin shapes, particle sizes and all Reynolds numbers combined.
6. Herringbone fins have 1.74 higher volume fouling rate than plain fin shape for low adhesive type particles. For high adhesive particles, herringbone has 1.8 times higher volume fouling rate and when both particle types are summed together, herringbone has 1.78 times higher volume fouling rate than the plain fin shape.
In the future, the investigated CFD-DEM method could be used for further studies on different fin shapes under different operating conditions and particle properties. Such results would be of high value in design of FTHE's. Novel topics for future research on fouling are different flow control strategies and fin shape optimization. Consequently, it is of constant interest to better manage the costs of FTHE life cycle.

Acknowledgments
The authors wish to acknowledge CSC -IT Center for Science, Finland, for computational resources and Koja Oy, for financial support.