At present, the treatment of municipal solid waste has gradually been upgraded from the primary treatment methods such as simple removal, field stacking and open-air simple landfill to advanced treatment methods such as composting fermentation, sanitary filling, biodegradation, incineration power generation and resource utilization. How to effectively and effectively control the garbage and prevent pollution has increasingly become a hot issue that needs to be solved in the pursuit of global concern. Among them, garbage drying technology is a key part of garbage disposal. This paper introduces the theoretical research and numerical simulation of the new two-phase flow garbage drying model. The main content is the setting of the simulation conditions of the two-phase flow model of rotary heat transfer. The maximum volume fraction of liquid phase and the maximum heat transfer efficiency of the drying device can be achieved under different parameter conditions by finite element analysis software, and the correctness of simulation results is verified by physical equipment.
Recently, the main structure of roller dryer is that the feeding end is slightly higher than the discharging end. The wet material enters from the high end and is dried by effective contact with hot air in the rotary kiln. With the slow rotation of the drum, the dried material is discharged from the low end under gravity. In Zheng Guangqiang and Zhou Yan  study, the cylinder and rotary kiln of the dryer were modeled and solved by the finite element method, and their mechanical analysis was studied. Considering that the wet material is heavier than the dried in both sides, a variable load is applied to the model to reflect the stress of the cylinder and rotary kiln more truly. In abroad, Canada’s A.S. Mujumdar  and others have studied roto-louver dryer. Roller dryer is another type of rotary dryer, in which the horizontal drum is equipped with longitudinal shutters, forming a conical drum in the external drum. The hot air passes through the shutters and liquefies at the bottom of the inner cylinder. Compared with the traditional rotary dryer, the dryer produces better heat and mass transfer rate and has a more compact mechanism. Because of the mild treatment of the product, the shutter dryer produces less friction, so it is suitable for fragile materials, but their cost is higher.
For the study of rotary drum motion and material movement, based on the mechanism of material space throwing, Li Miaoling  et al. established the numerical analysis model of material movement, deduced the formula of material detention time with the weight of lifting plate holding capacity, studied the law of material space movement in rotary dryer, established the analytical model of material movement in cylinder and the numerical analysis model of detention time; Zhu Guihua, Xu Hongwei , et al. used roller dryer to dry materials. When studying the influence of the internal structure of drum on material movement path and drying effect in the drum, the collision effect of materials and the traction effect of fluid on materials are considered comprehensively. F. T. Ademiluyi  et al. of Nigeria have studied the effect of drying parameters on the heat transfer of rotary dryers. The material is dried in a laboratory-scale rotary dryer at different inlet air temperature, inlet air speed, relative humidity, feed rate, driving speed and feeding speed. The results showed that the inlet temperature, inlet velocity, and feed rate have significant effects on the specific heat transfer coefficient and heat load of the material. A model for predicting the specific heat transfer coefficient of the heat load function is as a function of the inlet air temperature and velocity. Meanwhile the inlet air temperature, the inlet air velocity, and the feed rate, is also given. The predictions are compared with the experimental data, and good consistency between them is obtained.
The above-mentioned literature focuses on structural design and material handling for the research of rotary dryers, but does not study the heat transfer conditions and heat transfer medium and temperature. For example, when the drying material is a low-burning substance or an organic substance, an excessively high temperature causes combustion or carbonization; the drying heat source uses hot air or high-temperature flue gas, and the energy utilization rate is low. If the working fluid is changed into steam and water, on the one hand, the temperature inside the dryer can be effectively controlled, and on the other hand, the latent heat of vaporization can be fully utilized to improve the energy utilization rate.
Wang, W.  et al. applied the enhanced two-fluid model to study the influence factors of the coupled flow field between particles and fluids under various operating conditions, and considered the influence of virtual mass force. The mixture model regards fluid and particle phase as a continuous medium in the whole flow field, and particle phase is defined as pseudo-fluid. This method provides a reference for the study of the coupled flow field and its influencing factors on two-phase flow. Overseas, Guibo Li, Yongsheng Lian, Matt Mersch, Chris Omalley and Adam Hofmann  of the United States used the fluid transient interface capture to simulate the water sheet emitted from the nozzle of a fat blower and obtained the simulation data. Maria Adela, Puscas , et al. in France used the hybrid multiscale finite element to analyze two-phase flow. This method can perform the calculation of the basic functions including the small-scale variations of the permeability field in parallel and reduce the size of the global linear system. Two-stage MPI is used to adjust the computational resources of the algorithm to achieve better scalability of the method. This can bring significant acceleration to the two-phase model and hybrid multi-scale method. The above research has an excellent reference for the numerical simulation method of this paper. In the previous research, we carried out the following supplementary research: 1) Combining the two-phase flow model with the heat transfer model, using water and water vapor as heat transfer Working quality. 2) The temperature changes and heat transfer in the dryer under different inlet flow rates and inlet structures were investigated.
Finally, we draw the following conclusions: 1) The heat transfer coefficient in the dryer increases as the inlet flow rate increases. 2) The phase change rate in the dryer decreases as the inlet flow rate increases. 3) Increasing the number of inlets in the dryer structure can effectively increase heat transfer efficiency. 4) Within a certain speed, the dryer speed has little effect on the heat transfer results.
2. Establishment of Dryer Numerical Simulation
Creating an accurate 3D model of the drying device is the key to solving this numerical simulation problem. After the accurate measurement and data collection of the physical object, the geometric parameters without influence on the simulation results were corrected and a simplified three-dimensional model was established.
This paper uses UG to reconstruct the data features of the 3D model. The model is mainly divided into two three-dimensional entity parts:
Part one: The dryer drum is approximately a long cylinder as a rotating part.
Part two: Air inlet, water outlet and the heating housing.
The heating housing is tightly assembled with the cylinder, and the annular cavity is a fixed portion. The outer surface of the cylinder forms a contact surface with the inner surface of the housing. The heating fluid flows in from the inlet at one end of the housing and flows out at the other end. In order to make the simulation close to the real situation, two small cylinders are created as models for the inlet and the outlet. As shown in Figure 1.
An unsteady compressible flow equation is used depending on the flow properties of the fluid in the chamber and the flow of the two-phase flow. The K-epsilon turbulence model is used. Turning on the mixed multiphase flow model and the energy equation in FLUENT is to simulate the two-phase flow field in the cavity model and the process of fluid phase transition in the cavity.
2.1. The Standard k-ε Turbulence Model
The standard k-ε model is a semi-empirical model based on model transport equations for the turbulence kinetic energy (k) and its dissipation rate (ε). The model transport equation for k is derived from the exact equation, while the model transport equation for ε was obtained using physical reasoning and bears little resemblance to its mathematically exact counterpart.
Transport Equations for the Standard k-ε Model .
The turbulence kinetic energy, k, and its rate of dissipation, ε, are obtained from the following transport equations:
Figure 1. Simplified physical model.
The turbulent (or eddy) viscosity, , is computed by combining k and ε as follows
2.2. Mixture Model of Multiphase Flow
We use a hybrid model. The hybrid model is a simplified multiphase flow model that is used to simulate multiphase flows with different velocities in each phase, but assumes local equilibrium on a short spatial scale. The coupling between the phases should be very strong. It is also used to simulate isotropic multiphase flows with strong coupling and multiphase flows with phases moving at the same speed.
The hybrid model uses a single fluid approach:
1) The hybrid model allows the phases to run through each other (interpenetrating). Therefore, the volume fractions and for a control volume can be any value between 0 and 1, depending on the space occupied by phase q and phase p.
2) The hybrid model uses the concept of slip velocity, allowing the phases to move at different speeds. (Note: The phase can also be assumed to move at the same speed, and the hybrid model is reduced to a uniform multiphase flow model).
The hybrid model solves the continuity equation of the mixed phase, the mixed momentum equation, the mixed energy equation, the volume fraction equation of the second phase, and the algebraic expression of the relative velocity (if the phases move at different speeds).
2.2.1. Continuity Equation for the Mixture    
where is the mass-averaged velocity:
and is mixed density:
2.2.2. Momentum Equation for the Mixture  
The momentum equation for the mixture can be obtained by summing the momentum equations for all the phases.
is the mixed viscosity.
is the drift velocity for the second phase k:
2.2.3. Energy Equation for the Mixture 
The energy equation for the mixture takes the following form:
The first term on the right-hand side of Equation (2-11) represents energy transfer due to conduction. includes any other volumetric heat sources.
In Equations (2-14),
for a compressible phase, and for an incompressible phase.
2.2.4. Relative (Slip) Velocity and the Drift Velocity
The relative velocity (also referred to as the slip velocity) is defined as the velocity of a secondary phase (p) relative to the velocity of the main phase (q):
The mass fraction for any phase (k) is defined as
The drift velocity and the relative velocity ( ) are connected by the following expression:
FLUENT’s mixture model makes use of an algebraic slip formulation. The basic assumption of the algebraic slip mixture model is that to prescribe an algebraic relation for the relative velocity, a local equilibrium between the phases should be reached over short length scale. The form of the relative velocity is given by:
where is the particle relaxation time, according to Manninen et al. :
The default drag function is taken from Schiller and Naumann :
and the acceleration is of the form
The simplest algebraic slip formulation is the so-called drift flux model, in which the acceleration of the particle is given by gravity or a centrifugal force and the particulate relaxation time is modified to take into account the presence of other particles.
In turbulent flows the relative velocity should contain a diffusion term due to the dispersion appearing in the momentum equation for the dispersed phase. FLUENT adds this dispersion to the relative velocity:
2.2.5. Volume Fraction Equation for the Secondary Phases
From the continuous equation for secondary phase p, the volume fraction equation for secondary phase p can be obtained:
2.3. Temperature Boundary Conditions of the Coupling Wall
When a fixed temperature condition is applied at the wall, the heat transfer to the wall from a fluid cell is computed as
When you specify a convective heat transfer coefficient boundary condition is specified at a wall, FLUENT uses your inputs of the external heat transfer coefficient and external heat sink temperature are used by FLUENT to compute the heat flux to the wall as
Heat transfer to the wall boundary from a solid cell is computed as
2.4. Dynamic Mesh Update Method of the Rotation Area
For the mesh re-division of the cylindrical rotation model, we use the spring-based smoothing method. In the spring-based smoothing method, the edges between any two mesh nodes are idealized as a network of interconnected springs. The initial spacings of the edges before any boundary motion constitute the equilibrium state of the mesh. A displacement at a given boundary node will generate a force proportional to the displacement along all the springs connected to the node. Using Hook’s Law, the force on a mesh node can be written as
The spring constant for the edge connecting nodes i and j is defined as:
At equilibrium, the net force on a node due to all the springs connected to the node must be zero. This condition results in an iterative equation such that
Since displacements are known at the boundaries (after boundary node positions have been updated), Equation (2-27) is solved using a Jacobi sweep on all interior nodes. At convergence, the positions are updated such that
3. Simulation Results and Analysis
We use the fluid simulation software fluent to simulate the conditions of different inlet flow rates, which are 0.25 m/s - 0.75 m/s, and the calculation of different structure models, which are 1 entrance, 5 entrances, 7 entrances and 9 entries, the output data is imported into CFD-Post for viewing and analysis.
3.1. Temperature Conditions at Different Inlet Flow Rates
For the treatment of numerical simulation results, we first analyze the temperature changes and distributions of various parts of the fluid region at different inlet velocities. As shown in Figure 2, the temperature nephogram is used to analyze the temperature distribution at the horizontal plane near the center of rotation, i.e. the location of the inlet and outlet. The initial inlet temperature is 113 degrees (industrial high-temperature steam). When steam enters the chamber, the diffusion of the stream varies with the inlet velocity. It can be seen that when the inlet velocity is 0.25 m/s, the heat transfer is concentrated in the first half of the chamber, which has been completely cooled at 1.8 m. The outlet temperature
Figure 2. Temperature distribution nephogram in the central plane.
is normal; when the inlet velocity is 0.5 m/s, the fluid has basically cooled completely at 2.6 meters, and the outlet temperature is close to a normal temperature; when the velocity is 0.75, 1.0 and 1.5 m/s respectively, the fluid cannot be completely cooled in the chamber. The outlet temperature is higher than the room temperature, ranging from 310 K to 338 K.
Figure 3 draws a straight line with the center position in the 0.65 planes, thus temperature distribution curve is achieved at each inlet velocity. It can be seen that when the fluid enters the cavity, the stream directly contacts the surface with lower temperature, and the temperature decreases rapidly. From the contrast Figure 4, it can be seen that the turning point of the temperature drop slope occurs, which is just the peak value of the volume fraction of the liquid phase at that point. It shows that the latent heat of vaporization generated by the phase transformation is the main factor leading to the temperature decrease. After that, the rate of temperature change becomes slower, which is due to the fact that most of the vapor phase changes into condensate water after phase transformation, which is stored under the cavity, and the mass of vapor participating in wall heat transfer decreases greatly. Meanwhile, after 0.6 m, the curves of each inlet velocity begin to separate. When the inlet velocity is 0.25 m/s, the inlet velocity is the smallest, and the flow rate at the same time is the smallest. The heat transfer between the gas phase and the wall is more sufficient. At 1.76 m, the heat transfer process of the stream is completed, and the temperature reaches room temperature. When the inlet condition is 0.5 - 1.5 m/s, the temperature change rate decreases with the increase of the inlet velocity. The outlet temperature increases accordingly.
3.2. Phase Transitions at Different Inlet Flow Rates
Looking at the phase transition curve Figure 4, at 0.294 m, the volume fraction of the liquid phase is the peak value. It can be seen that the slower the inlet velocity is, the lower the peak value. Since the slower inlet velocity, the better the
Figure 3. Temperature variation curve of center position at the z-axis.
Figure 4. Liquid volume fraction curve at z-axis at the center position.
transfer near the inlet, the larger amount of liquid phase is stored in the first half, such as at the inlet velocity of 0.25 m/s, before 1.76 m, a large amount of liquid phase is accumulated at the bottom of the wall and is fully completed liquefaction; in contrast, the faster the inlet velocity is, the less liquid phase can be fully liquefied in the first half. For example, in the case of 1 m/s inlet velocity, the volume fraction of the liquid phase in the first half is much less than that of 0.25 m/s, but till the outlet, liquefaction always exists in the gas phase, and at 1.2 m, the volume fraction of liquid phase exceeds that of 0.25 m/s.
Figure 5 is the average volume fraction of liquid at each inlet. It can be seen that when the inlet flow rate is less or equal to 0.5, the volume fraction is small and the phase transition is sufficient. Within the range of 0.5 to 1 m/s, the average volume fraction at the outlet rises sharply. When the inlet flow rate is more than 1 meter per second, the volume fraction the outlet rises slowly and the slope of the ascending curve is lower.
Figure 5. The average outlet liquid volume fraction curve for each inlet velocity.
3.3. Analysis of Wall Heat Transfer
As shown in Figure 6, the average heat flux at each inlet flow rate increases linearly with the increase of the inlet flow rate (the turning point at 1 and 1.5 m/s is caused by coordinate changes).
By calculating the Nusselt number on the heat transfer wall and analyzing the variation curve for the simulation results of different inlet velocities, it can be concluded that with the increase of the inlet velocity, shown in Figure 7, the Nusselt number on the heat transfer wall increases linearly, from 62 to 312, and the comprehensive heat transfer level is not high.
3.4. Separate Analysis of Phase Transitions under Individual Conditions
Since the large difference in the liquid volume fraction nephogram, it is difficult to show in a single diagram. Now the inlet velocity is set with 0.25, 0.75, 1.5 m/s respectively for separate analysis.
Figure 8 shows the liquid volume fraction nephogram of the bottom surface at an inlet velocity of 0.25 m/s. It can be seen that phase transition occurs in the first half of the cavity due to the slow inlet velocity, and a large number of liquids are stored in the central part of the first half. The phase transition is fully completed and the gas phase cannot reach the outlet.
As shown in Figure 9, the liquid volume fraction cloud at the inlet velocity of 0.75 m/s shows that the fluid diffuses rapidly after entering the chamber. The phase transformation process mainly occurs at 2.1 m. The gas phase cannot be completely transformed into a liquid phase in the chamber, and the outlet contains the gas phase.
Figure 10 is a liquid volume fraction nephogram with an inlet velocity of 1.5 m/s. It can be seen that the gas phase directly reaches the back surface of the cavity due to the fast inlet velocity and weak diffusion of the fluids going into the inlet, and then diffuses and flows going back to both sides. Phase transformation
Figure 6. Average wall heat flux curves at various inlet velocities.
Figure 7. The Nusselt number of heat transfer walls at different inlet velocities.
Figure 8. Liquid volume fraction nephogram of bottom surface at inlet velocity of 0.25 m/s.
Figure 9. A volume fraction nephogram of the liquid phase at an inlet velocity of 0.75 m/s.
occurs at a higher position of the cavity, and the phase transformation concentrates near the outlet. The overall liquid volume fraction is low and the outlet is low. It contains a lot of gas phase.
After calculating the average liquid volume fraction and the Nusselt number of the heat transfer wall surface, the phase transition rate of the non-rotating model is 99.051% of that of the rotating model, and the Nusselt number of the model without rotation is 98.8 percent of the rotational model.
3.5. Structural Optimization of the Dryer Model
Aiming at the problem of low overall heat transfer coefficient and low comprehensive heat transfer level, this paper optimizes the structure of the model in order to achieve higher heat transfer efficiency. After analyzing the simulation results of the model fluid, we found that the main reason for the low Nusselt number of the heat transfer wall is that the entrance is smaller than the whole model. The entrance is concentrated in the middle of the model, and the space on both sides of the cavity is not available. To make full use of heat exchange, we have changed the entrance model to 5, 7, 9, shown in Figure 11, and arranged them evenly. As shown in Figure 12, the numerical simulation results show that the Nusselt number of the heat transfer wall is greatly increased at 1054, 15,741, 19,587 respectively.
3.6. A new Number to Analyze the Comprehensive Benefits
In summary, although the heat and heat transfer coefficient increase with the increase of inlet velocity, when the velocity is too high, the phase transformation of fluid is insufficient, and the utilization rate of heat is not high. Therefore, it is
Figure 10. A volume fraction nephogram of the liquid phase at an inlet velocity of 1.5 m/s.
Figure 11. One of the models of structural optimization.
necessary to find a balance between the two conditions. So we define the value CED (Comprehensive heat transfer efficiency of dryer model) to represent the comprehensive heat transfer efficiency of the dryer model. Its expression is:
The results are shown in Figure 13, and it can be concluded that the value is between 0.75 meters per second to 1.0 meter per second.
In this study, we established a 3D model which is consistent with the dryer and calculated the temperature distribution, phase transformation and heat transfer of the fluid under different boundary conditions by using the fluid dynamics simulation software.
The results show that the phase transformation mainly occurs at the inlet and the first half of the cavity. When the inlet velocity is 0.5 m/s or less, the water
Figure 12. Nusselt number of heat transfer walls with different structures.
Figure 13. The CED number in different inter flow velocity.
vapor in the inlet basically completely transforms into a liquid. The outlet temperature is close to normal temperature, which is not conducive to the overall heat transfer. However, the liquid phase transformation rate is high and the utilization rate of heat is high. When the inlet velocity reaches 0.75 m/s or above, heat transfer occurs in the middle and rear parts of the cavity. The gas phase of the inlet cannot be converted into a liquid. The outlet temperature is obviously higher than the normal temperature. Between 50 and 60 degrees Celsius, there are a large number of gas phases at the outlet. When the inlet velocity is too high, at 1.5 m/s or above, a large number of fully heat-exchanged gases impinge on the back surface of the cavity and move towards the cavity. There are relatively more phase transitions near the exit and in the back half of the cavity, and the gas diffusion in the first half is weak, which is not conducive to the overall heat dissipation of the model.
The results show that the phase change rate of the fluid can be obtained by detailed analysis and conversion of phase change parameters and conversion of volume fraction to mass fraction. However, since most of the liquid is stored on the surface of the cavity, the phase change rate of the fluid cannot be obtained by the volume fraction of the liquid at the outlet. We estimate the phase change rate of liquid by the volume fraction of the liquid at the outlet. The inlet condition is 0.25 m/s. The phase transition rate is 100% at the same time. From the curves in the results, it is found that the phase transition rate decreases significantly from 0.5 to 1 m/s at the inlet velocity.
The results show that for the wall heat flux, with the increase of the inlet velocity, the heat flux increases linearly and tends to be stable when the velocity reaches a certain value. Correspondingly, the wall heat transfer coefficient of the cylinder also increases.
The results show that the heat transfer rate and heat flux of the model decrease, but the influence is small. In practice, objects with different internal shapes change direction and position with the rotation. This part of the situation is difficult to be reflected in the simulation software, which makes it different from the actual situation and belongs to a reasonable error.
After optimizing the inlet structure, the model shows a high heat transfer efficiency in the numerical simulation results. Especially when the number of inlets reaches 7, the heat transfer wall Nusselt number exceeds ten thousand. The inlet structure shows a different in heat transfer efficiency of the model.
This work is financially supported by Natural Science Foundation of Shandong Province (Grant No.ZR2017BEE012) and Shenzhen science and technology project (JCYJ20170818103826176) and Focus on research and development plan in Shandong province (Grant No.2018GGX104011).
: A constant. .
: The generation of turbulence kinetic energy due to the mean velocity gradients.
: The generation of turbulence kinetic energy due to buoyancy.
: The contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate.
: Constant. , .
v: The component of the flow velocity parallel to the gravitational vector.
u: The component of the flow velocity perpendicular to the gravitational vector.
: The turbulent Prandtl number for K equation. .
: The turbulent Prandtl number for ε equation. .
: Fluid density.
k: Turbulent kinetic energy (J).
ε: Dissipation rate of turbulent pulsating kinetic energy.
: Speed in the Xi direction.
: User-defined source terms.
p: Pressure on fluid micro-element.
: Average speed of mass.
: Mixed density.
: The mass transfer due to cavitation or user-defined mass sources.
: The volume fraction of the kth phase.
n: The number of phases.
: The volume force.
: The drift velocity for the second phase k.
: The mixed viscosity.
: The effective conductivity ( , where is the turbulent thermal conductivity, defined according to the turbulence model being used).
: The particle relaxation time.
: The diameter of the particles (or droplets or bubbles) of secondary phase p.
: The secondary-phase particle’s acceleration.
: The mixture turbulent viscosity.
: A Prandtl dispersion coefficient.
: The sensible enthalpy for phase k.
: Fluid-side local heat transfer coefficient.
: Wall surface temperature.
: Local fluid temperature.
q: Convective heat flow from the wall.
: Radiative heat flux.
: External heat-sink temperature defined by you.
: External heat transfer coefficient defined by you.
: Thermal conductivity of the solid.
: Local solid temperature.
: Distance between wall surface and the solid cell center.
: The displacements of node i.
: ’s neighbor j.
: The spring constant (or stiffness) between node i and its neighbor node j.
n: The positions at the current time step.
Nu: Nusselt number.
: Liquid phase volume fraction at the outlet.
: Liquid flow rate at the inlet.
 Zhu, G., Xu, H., Wang, J., Liu, Y. and Peng, N. (2018) Characteristics Simulation and Structure Optimization of Drum Sludge Dryer. Journal of Jiangsu University (Natural Science Edition), 39, 556-562.
 Ademiluyi, F.T., Abowei, M.F.N., Puyate, Y.T. and Achinewhu, S.C. (2010) Effects of Drying Parameters on Heat Transfer during Drying of Fermented Ground Cassava in a Rotary Dryer. Drying Technology, 28. https://doi.org/10.1080/07373931003622669
 Wang, W., Li, Z., Zhang, J., et al. (2018) Pap. Simulation Study of Particle-Fluid Two-Phase Coupling Flow Field and Its Influencing Factors of Crystallization Process. Chemical Papers, 72, 3105-3117. https://doi.org/10.1007/s11696-018-0537-0
 Li, G., Lian, Y., Mersch, M., Omalley, C. and Hofmann, A. (2014) Liquid-Gas Two-Phase Flow Simulation for Flat Fan Nozzles. 4th ASME Joint US-European Fluids Engineering Division Summer Meeting. https://doi.org/10.1115/FEDSM2014-21170
 Puscas, M.A., Enchéry, G. and Desroziers, S. (2018) Application of the Mixed Multiscale Finite Element Method to Parallel Simulations of Two-Phase Flows in Porous Media. Oil & Gas Science and Technology. Rev. IFP Energies Nouvelles, 73, 14. https://doi.org/10.2516/ogst/2018022
 Shao, N., Salman, W. and Gavilidis, A. (2008) CFD Simulations of the Effect of Inlet Conditions on Taylor Flow Formation. International Journal of Heat and Fluid Flow, 29, 1603-1611. https://doi.org/10.1016/j.ijheatfluidflow.2008.06.010
 Altimira, M., Rivas, A. and Larraona, G.S. (2009) Characterization of Fan Spray Atomizers through Numerical Simulation. International Journal of Heat and Fluid Flow, 30, 339-355. https://doi.org/10.1016/j.ijheatfluidflow.2008.12.006
 Li, Z.D., Zhang, G.Q. and Li, Z. (2008) Simulation of Gas Flow Field in Laval Nozzle and Straight Nozzle for Powder Metalurgy and Spray Forming. Journal of Iron and Steel Research, International, 15, 44-47. https://doi.org/10.1016/S1006-706X(08)60264-2
 Tong, M. and Browne, D.J. (2009) Modelling Compressible Gas Flow near the Nozzle of a Gas Atomiser Using a New Unified Model. Computers & Fluids, 38, 1183-1190. https://doi.org/10.1016/j.compfluid.2008.11.014