{'fraction_non_alphanumeric': 0.025497425366749118, 'fraction_numerical': 0.01043854327095098, 'mean_word_length': 4.581488234783364, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Algorithm design and the development of a discrete element method for simulating particulate flow and heat transfer

An algorithm using the discrete element method (DEM) for simulating the particulate behaviour of flow and heat transfer is developed and described, the reasonable hypothesis and the ingenious design of which have been presented in detail. The organizational structure of the developed algorithm contains an efficient method for determining particle collisions, the status analysis for each particle and the particulate-kinematics analysis during the time step. The reasonability and correctness of the developed DEM algorithm are validated by laboratory experiments: the discharge process of glass beads from a silo; and heating of metal alloy particles in a calciner. Afterwards, a group of validated mechanics parameter values for coal and sand have been tested and verified in the article, preparing for the simulation of the pyrolysis process in a downer or screw reactor in subsequent research projects.


Introduction
Granular materials in nature and particulate technology in process industries are globally researched and developed because the interactions among individual units in particulate systems are so complicated that the behaviours of particulate-laden gases are quite different from those of conventional fluids. Thus, it is important for real interdisciplinary research on particulate systems to understand the mechanism at the microscopic level.
As a well-documented numerical tool and a prime or outstanding example of a dynamic algorithm, the DEM algorithm was originated by Cundall in 1971 [26] and applied to simulate the behaviour of soil particles under dynamic loading conditions by Cundall and Strack in 1979 [27]. The DEM algorithm is more efficient in quasi-static systems because it is capable of handling multiple particle contacts. The particle deformation can be explicitly determined by the interaction force of collisions. The transient contacting surface of particle-to-particle collisions can be confirmed by the particle positions determined by the wellestablished Newton's laws of motion. That is to say, the DEM algorithm is a powerful numerical method, in which the motion of each individual particle is determined by the net force acting upon it. Another advantage of the DEM algorithm is that it can present more detailed information on each particle, such as its position, trajectory, velocity and contact force, which are usually not easy to gain from experiments. In addition, the transient contacting surface can be modified to consider the heat-transfer area of two colliding particles with different temperatures. Thus, the DEM algorithm coupled with thermodynamics and heattransfer models can simulate the transient heat-transfer phenomena in a particulate system. Furthermore, the particulate macro-chemical-reaction kinetics also can be introduced to the kinetics-analysis program to describe the variance of density or size at the particle scale. Therefore, the soft-sphere model of the DEM algorithm can provide an effective numerical method to simulate the particle behaviour of flow and heat transfer in a particulate system.
To utilize the DEM algorithm to serve the simulation of particulate processes extensively and openly, for example, the mixing and heat transfer of coal and sand during the coal-pyrolysis process in a downer reactor or screw reactor, an integrated DEM program package has been developed and presented in detail. For one thing, the detailed structure of the DEM algorithm, such as the searching algorithm for determining the particle collisions, has not been introduced in previous publications. For another, the parameter values are meaningless if the contact-mechanics model or heat-transfer model is unknown or undocumented in some simulation tools. The contact-mechanics model and heat-transfer model should be stated clearly and the physical meaning and values of the parameters should be well understood.
During the algorithm design of the DEM, a new method for cubic-grid formation and connection was developed, especially for application to irregularly shaped simulation domains, which can save the overhead and computation time of systems. The particle-dynamics model in the kinematics analysis of the DEM program has been supplemented by the particle-rotation model controlled by the rigid-body rotation's law to simulate the angle of repose of a particle pile accurately. Most importantly, to simulate coal pyrolysis in a downer reactor, the particle-heattransfer model and particle-thermodynamics model have been coupled to the developed DEM algorithm. Second, the correctness and reasonability of the developed DEM algorithm have been validated and checked through simulating the conditions compared with experiments. In this research, the mechanics model for particle-status analysis and the heat-transfer model for particle-kinematics analysis of the developed DEM algorithm have been validated by an experiment with glass beads discharging from a silo and metal particles heating in a calciner. Finally, to prepare the simulation of the coal-pyrolysis process in a downer reactor using the developed DEM algorithm in the next step, some necessary physical-property parameters of the fuel particles of coal and the heat-carrier particles of sand have been validated and confirmed through comparing the experiments of particle filling in the tube, particle heaping on the platform and particle flowing through an hourglass for timing.
1 Development of the DEM algorithm

Hypothesis and principle
In the DEM algorithm, each element can be considered as one rigid body with a fixed shape, such as a sphere, a billet or a pyramid. A little overlap is allowed to happen when two elements contact each other, which is calculated by a force-interaction model. As a dynamic algorithm, the status of all elements is unchanged in each time step, such as force analysis and thermodynamics analysis. The foundational principle or big picture of the DEM algorithm is described in Fig. 1.
As shown in Fig. 1, the DEM algorithm, as a Lagrangian method for modelling the individual trajectory of each element, can be described as a dynamic algorithm, whose structure consists of three parts: the status analysis, the kinematics analysis and an algorithm for the prediction of particle collisions. In the particle-status analysis, the particle-mechanics model calculates the particle-force analysis and the particle-torque analysis, and the particlethermodynamics model confirms the particle temperature according to the heat balance of one single particle. In the particle-kinematics analysis, the particle-dynamics model takes charge of simulating the particle position and the velocity of the particle in the next time step, and the heattransfer model can calculate the quantity of heat transfer between particles or from a particle to its surroundings. For a large number of particles in the system, it is a major task to detect the particle collisions and to calculate the contact forces acting on the particle. A search algorithm based on

Searching technology
Time driven: t = t + d t Bing Liu | 143 cubic grids for particle collisions is adopted in the whole program to realize the particulate three-dimensional softsphere model of the DEM on the computer more efficiently. The whole calculation process of the dynamic algorithm of the three-dimensional soft-sphere model of the DEM will not be finished until the calculation step approaches the end.

Particle-status analysis
Each particle has its own status parameters, such as the space-position vector x = (x 0 , , particle temperature T p and its enthalpy value H p . All the status parameters of each particle simulated in the DEM algorithm should be determined at the end of each time step.

Mechanics model
Normal-contact-force model. It is very important to calculate the net force and total moment of each particle when particles contact each other, as illustrated in Fig. 2, because accurate particle motion is determined by the exact particle-mechanics model introduced in this section.
To calculate the normal contact force acting on each particle, the Hertz theory [28] is applied to obtain the relationship between the normal contact force and the particle deformation, which is summarized as In the Hertz theory, the normal stiffness coefficient k n is associated with the physical parameters of particles, such as elastic modulus E, Poisson's ratio ν and the particle radius of two colliding spherical particles. e n is the normalunit vector of the contact surface, which can be obtained by the position vectors of two colliding particles as The calculation method of the normal-unit vector can be obtained from Fig. 2. The normal damping coefficient η n can be calculated by [29]: where the reduced mass m * can be calculated by two colliding particles; e represents the recovery coefficient according to the initial velocity and the rebound velocity of one particle; and v n,ij is the normal relative velocity of particle i relative to particle j.
Tangential-contact-force model. Vemuri theory [30] is employed to obtain the tangential interactions between the colliding particles, which can be expressed as follows: However, the calculation method of the tangential force between two colliding particles is complicated because the physical parameters cannot be tested easily in experiments. Therefore, the following methods [31,32] are used to calculate the tangential force between the colliding particles: where μ s is the coefficient of side friction; F cn,ij is the same as mentioned above as the normal contact force of two collision particles; δ t,ij is the relative displacement of two collision particles at contact time Δt; and e t is the unit tangential vector of the contact surface, which can be obtained by The calculating method for the normal-unit vector can be obtained from Fig. 2. η t, ij is the tangential damping coefficient, which can be calculated by where v t,ji is the tangential relative velocity of particle j relative to particle i. Conservative forces, such as gravity or electromagnetic forces, acting on the centre of a homogeneous spherical particle cannot cause the particles to rotate. However, the line of total contact force acting on the particle surface does not go through the centre of a spherical particle, which makes the spherical particle rotate. The total rotation moment is composed of the tangential contact torque and the rotational damping moment: where r i is the vector of the radius of particle i, whose direction is from the centre of the spherical particle to the contact point of collision; F t,ij is the tangential component of the contact force from particle j to particle i; and μ r is the damping coefficient of the rolling friction.

Thermodynamics model
In this research, the enthalpy of a solid particle is assumed to be a function of the temperature, i.e. the enthalpy H of a particle at temperature T can be calculated from that at the specified temperature T 0 , i.e.: (11) where m p is the mass of the particle; and c p is the heat capacity of the particle at constant pressure, which can be calculated as c p = a + bT + cT 2 . Therefore, the enthalpy change ΔH of the particle can be calculated as where parameters a/kJ· kg -1 ·K -1 , b/kJ· kg -1 ·K -2 and c/kJ ·kg -1 ·K -3 are constants collected in the thermodynamics handbooks. Therefore, the change in particle temperature T i can be calculated from the enthalpy difference ΔH i through the heat balance of particle i:

Dynamics model
The translational movement of the particle is driven by the net force acting on it, and accords with Newton's Second Law of Motion. The rotating movement of the particle is driven by the total moment acting upon it and is described by Newton's Law of Rotation: According to the assumption mentioned above, the acceleration of the particle motion in a small-enough time step is a constant and the angle of acceleration of the particle rotation in a small-enough time step is also a constant. The velocity vector of the particle translational movement v i and the angular velocity vector of the particle rotation ω i at the next time step can be expressed by the first-order Taylor series expansion at t = t 0 . The spatial-position vector of particle x i and the phase-position vector of particle ϕ i at the next time step can be derived by the secondorder Taylor series expansion at The particle-kinematics equations of translational motion and rotational motion can be discretized into a Verletleapfrog scheme as in Equations (16) and (17):

Heat-transfer model
The calculation method for the heat transfer between two colliding particles or between a particle and a wall with different temperatures can be obtained by the theoretical analysis of Fourier's law, which can be described as follows: where |ΔT ij | is the absolute temperature difference between two colliding particles or between a particle and a wall standing for the impetus of heat transfer through thermal conduction; A ij is the area of the thermal conduction between two colliding particles or between a particle and a wall; α c, ij is the contact coefficient of thermal conduction; and the factor of 1 α c,ij A c,ij stands for the resistance of heat transfer. It is worth noting that the contact coefficient of thermal conduction is affected by the roughness of the contacting surface, the acting force on the contacting surface, the gas pressure in the gap between two colliding contact surfaces, etc. It cannot be tested by experimentation for most engineering materials except for some highly polished metals. Therefore, the calculation method for the contact coefficient of thermal conduction is the core of the modelling of particulate thermal conduction. The equivalent thermal-resistance model of the heat transfer between contact particles with different temperatures is illustrated in Fig. 3.
As shown in Fig. 3, the particle can be assumed to be surrounded by a thin gas membrane with thickness δ, which is much less than the particle diameter d p . According to the conclusion of Delvosalle from the measurement in experiment [33], the thickness of the assumed gas layer can be considered as 0.1 times the particle diameter, i.e. δ = 0.1d p . Thus, the thermal conduction between two colliding particles with different temperatures is affected by the thin gas layer around the particle. Further assumption should be supplemented to build up the model of particulate thermal conduction. First, the contacting surface is smooth, whose normal vector is parallel with the centre line of the two colliding particles, the direction of the heat transfer through thermal conduction. Second, the area of the contacting surface is equal to that of the overlap circle based on the assumption of DEM modelling. Third, there is a minimum tiny gap interval with a thickness of 4 × 10 -10 m between the contacting surfaces. Last, the thermal conduction can be transferred through the internal minimum tiny gap and the annular overlap of the gas layer around it. Therefore, the heat transfer through thermal conduction between two particles with different temperatures can be discussed from two conditions, as follows.
(i) The heat transfer based on thermal conduction will not occur when two particles are far away from each other, i.e. l 12 ≥ r 1 + r 2 + δ. l 12 is the distance between two particle centres; r 1 and, r 2 are the radii of two particles separately. However, thermal conduction will not begin until the overlap of the gas layers of two particles with different temperatures appears, i.e. r 1 + r 2 < l 12 < r 1 + r 2 + δ, as illustrated in Fig. 3a-c. At this time, the heat transfer through thermal conduction is only controlled by the overlap of the gas layers of the two particles, whose thermal resistance R g can be calculated as follows: where k g is the thermal-conduction coefficient of the gas layer; and Δl is the distance between two relevant particle surfaces along the direction of the centre line in metres. α i is illustrated in Fig. 3b and can be calculated as follows: A similar calculation method can be applied to the thermal resistance between a particle and a wall with different temperatures through overlapped gas layers, as follows: where l pw presents the distance between the particle centre and the wall surface along the direction of the normal vector of the wall surface.
(ii) The thermal conduction will be enhanced when two particles with different temperatures collide with each other, i.e. l 12 ≤ r 1 + r 2 . The heat can be conducted through not only the overlapped gas layers, but also the tiny gas interval d gap between the contacting surfaces of two particles with different temperatures. Thus, the thermal resistance is composed of two parts: the thermal resistance of the tiny gas interval R g1 and the thermal resistance of the annular gas layer R g2 . Therefore, the thermal resistance under this condition can be calculated as follows: where α i and β i (illustrated in Fig. 3e) can be calculated as follows: A similar calculation method can be applied to the thermal resistance between a particle and a wall with Fig. 3: Illustration of the calculation of the thermal resistance of thermal conduction between two colliding particles different temperatures through the gap interval and the annular overlapped gas layers during the contact of a particle with a wall with different temperatures as follows: where l pw presents the distance between the particle centre and the wall surface along the direction of the normal vector of the wall surface.
As known to all, the thermal flux can be also affected by the solid resistance of the particle's interior, which can be calculated as follows: where k p is the thermal-conduction coefficient of the particulate material; r i presents the radius of the particle; and . The spherical surface a radius of r m divides the internal energy of the particle equally.
Above all, the total thermal resistance of heat transfer through thermal conduction between two particles with different temperatures can be calculated as follows:

Formation of cubic grids
Detecting colliding particles is an important component of the algorithm of the three-dimensional soft-sphere model of the DEM, which makes up a large percentage of the computation time. In the particle-status analysis, calculation of the total contact force acting on one particle requires calculation of the contact forces from all the particles colliding with it. Therefore, it is essential to determine all particles that collide with the current one at each time step. To save calculation time, a series of search algorithms have been developed and applied in the soft-sphere model of the DEM, such as the Neighbour List Method [34], Bounding Box Method [35,36] and Boxing Method [36]. In this research, a search algorithm based on cubic grids has been developed to execute the task of determining particle collisions in the three-dimensional soft-sphere model of the DEM.

Cubic-grid formation in a cuboids space.
Before the calculation program of the three-dimensional soft-sphere model of the DEM is executed, the entire cuboids space of the particle motion will be meshed into a series of cubic units, as shown in Fig. 4.
In a cuboids space with a constant volume, the scale of one cubic unit determines the distribution density of the cubic grids. The scale of one cubic unit is determined by its side length. The distribution density of cubic grids in a constant volume can be determined by the particlesize distribution and the density of the particulate phase. Generally, a small scale for the cubic unit is suitable for the dense phase of particle flow and a larger side length of the cubic unit is acceptable with the dilute particle flow. According to the search algorithm based on cubic grids, the side length of the cubic unit is usually larger, but not too much larger, than the maximum particle diameter, i.e. l box > d p, max .
Each cubic unit has its own serial number (ID). As presented in Fig However, it is obvious that half of the neighbours of one cubic unit should be recorded in order to save computer-memory space, because the neighbour relationship is mutual. For example, cubic unit U is the neighbour of the current box unit coloured in pink and the current box unit is also the neighbour of cubic unit U at the same time. When the neighbours of cubic unit U are searched, the current box unit coloured in pink can be omitted, because cubic unit U has been searched when the cubic unit coloured in pink is treated as the current one. Therefore, those in the 26 neighbour units whose ID numbers are larger than (or less than) that of the current cubic unit can be considered as the neighbours recorded in the data structure of the current one. In this research, the uncoloured neighbours whose ID numbers are less than that of the current cubic unit coloured in pink (D, DNW, DN, DNE, DE, DSE, DS, DSW, DW, W, SW, S, SE) are omitted to save calculation time. In other words, the spatial structure consisting of cubic units can be built up and described by the 13 neighbours of the current unit coloured in pink (NW, N, NE, E, U, UNW, UN, UNE, UE, USE, US, USW, UW) in the memory of the computer and the criteria for the formation rule mentioned above have been listed in Table 1.
In Table 1, ID_BOX indicates the series number of the cubic unit; COLUMN shows the column number of the cubic unit in one layer of the whole cuboids space; LAYER represents the number of the cubic number in one layer of the whole cuboids space. If one of the recorded neighbours does not exist, the linker field of the corresponding neighbour of the current cubic will be set as NULL.
Splicing cubic grids between two cuboids spaces. The number of cubic units in the cubic grids seriously affects the memory overhead and the calculation speed of the computer. The spatial domain of the particle motion is usually an irregular one, i.e. the space of particle motion is not usually a cuboids space. If the cubic grids are meshed all at once in a larger cuboids space involving an Bing Liu | 147 irregular shape space of particle motion, a large number of redundant cubic units will be generated in the memory of the computer and a large quantity of computational time will be wasted. Therefore, an irregularly shaped space can be divided into several regularly or quasi-regularly shaped subspaces, in which the cubic grids can be meshed. Therefore, the whole irregularly shaped space of particle motion can be recombined through splicing the cubic grids in the regularly or quasi-regularly shaped space of particle motion together, as described below.
The design of the splicing algorithm for cubic grids is completely determined by the partitioning of the irregularly shaped space of particle motion that can be segmented along the direction vertical to D-U into two parts: the segment part U and the segment part D. Each cubic unit in the fracture surface of segment part D is spliced with UNW, UN, UNE, UW, U, UE, USW, US and USE as neighbours from the segment part U, and the one in the fracture surface of segment part U does not need to be spliced with any cubic unit in segment part D according to the cubic-grid-formation rule mentioned in the 'Cubic grids formation in a cuboids space' section. If the segmentation direction vertical to W-E is adopted, the segment part W and segment part E are the two subspaces of particle motion. Each cubic unit in the fracture surface of segment part W is connected with UNE, UE, E, NE and USE as neighbours from segment part E, and the one in the fracture of segment part E is needed to link to USW, UW, UNW and NW as neighbours from part W. The third cutting direction can be perpendicular to S-N and the subspaces of segment part S and segment part N are separated from the irregularly shaped space of particle motion. All the cubic units in the fracture surface of segment part S are attached to UNW, UN, UNE, NW, N and NE as neighbours from segment part N, and those of segment part N are needed to adjoin USW, US and USE as neighbours from segment part S. It can be concluded from these three segmentation cases that nine neighbour relationships should be connected. However, the first segmentation case is beneficial for the design of the splice algorithm, because the cubic units on the fracture surface of only one segment part D are matched with the neighbours from segment part U. The algorithm can be described as below.
In the fracture surface, the series number of the first cubic unit in segment part D is begin 0 ; step E0 and step N0 are count steps along the direction of E and N in segment part D. Similarly, the series number of the first cubic unit in segment part U is begin1; stepE1 and stepN1 are count steps along the direction of E and N in segment part U. If the fracture is a rectangle with row lengthN and column lengthE, the splice units on the two fracture segment parts are obj 0 and obj 1 , which can be calculated as: For the same value of i and j for the cubic grids of two fracture surfaces, the corresponding series number of cubic units obj 0 and obj 1 in two segment parts are just the two splice units, also a neighbour relationship, along the D-U segmentation direction. Additionally, eight neighbours in segment part N of the current cubic unit in segment part U can be gained according to the criteria list in Table 2. The neighbours will be NULL if the conditions are not satisfied with the criteria listed in Table 2.

Association between particles and cubic grids
The program of the three-dimensional soft-sphere model of the DEM will launch the dynamic mechanism after the formation of cubic grids. In each time step, all particles are associated with the corresponding cubic units according to their spatial-position vectors, in which particle collisions will be determined. In this section, two problems are addressed: (i) how to find the series number of the cubic unit with which the current particle is associated; (ii) how to record or remark the associations between particles and cubic units.
Mapping from the position of a particle to the series number of the cubic unit. It is not difficult to imagine that a particle can be assigned to the corresponding cubic unit according to the position of the particle. The mapping from the position of the particle to the series number of the corresponding cubic unit can be built up based on the cubic grids created, as  The criteria condition of the splice cubic unit between two segment parts of an irregularly shaped space of particle motion Neighbours Criteria described in Section 1.4.1. In a regularly shaped space of particle motion, the series number of the cubic unit is increased from the W direction to the E direction stepped by 1, from the S direction to the N direction stepped by the column number, and from the D direction to the U direction stepped by the layer number. Therefore, the mapping from the position of the particle to the series number of the cubic unit can be built up as: where (x a , x b , x c ) is the position of the vertex of units D, W, and S of the first cubic unit in the cubic grids whose series number is 0; l a , l b and l c are the three side lengths of a cubic unit; n column stands for the column number in each layer of cubic units; and n layer represents the number of cubic units in each layer of the cubic grids.
Algorithm of association between the particles and the cubic grids.
The data structures of the particle and the cubic unit are illustrated in Fig. 5. The data structure of the particle not only records the basis information of the particle, such as the radius, density, mass, rotational inertia, position vector, translational-velocity vector, rotating-velocity vector, force vector and moment vector, but also includes a pointer member of the particle type. When one cubic unit contains more than one particle, the particletype pointer field will link other particle units to form a chain table in a series of increasing particle numbers. The structure of the cubic unit has the pointer members of 13 neighbours of the cubic-unit type and a pointer member of the particle type. The pointer member of the particle type is linked to the particle unit of the minimum series number of all particles matched with this current cubic unit.
The particles are associated with cubic grids all over again in each new time step because of the change in the spatial position of all particles. To avoid influence from the last association between particles and cubic grids, all pointers of a particle type in the data structure of particles and cubic units should be set as NULL before the association mechanism of the next time step is launched. When the series number of the cubic unit is calculated from Equation (29) according to the spatial position of the particle, the particle unit will be linked to the pointer of the particle type in the cubic-unit data structure. If there are some other particles associated with the same cubic unit, a chain table of particle units whose heads are linked to the associated cubic unit's pointer of particle type should be formed in ascending order of the series number of the particle unit. The flow sheet of the pseudo code is presented as Fig. 6.
An example has been prepared for illustrating the association between the particles and the cubic units, as shown in Fig. 7.
When the subroutine that records the association between particle units and cubic grids is launched, a tree structure like that in Fig. 7 will be generated and updated in the computer memory, which is beneficial for the realization and design of the algorithm for detecting particle collisions. As represented in Fig. 7a, the cubic unit of 628 # contains particle units of 428 # , 432 # , 620 # , 656 # and 697 # ; those of 617 # , 751 # and 758 # are associated with the 629 # cubic unit; the 630 # cubic unit involves only the 302 # particle unit. The records of the association between the particle units and cubic grids in the computer memory can be shown as in Fig. 7b. It is clear that the series number of particle units in the chain table is increased from head to tail. It is worth noting that the 629 # cubic box unit is the east (E) neighbouring unit of 628 # , and the 628 # cubic unit is the west (W) neighbouring unit of 629 # . The west (E) neighbouring unit is invalid according to the principle mentioned in Section 1.4.1.Therefore, the 628 # cubic unit is not linked to one of the 13 neighbour pointer fields of 629 # . The same principle acts on 629 # and 630 # .

Algorithm for detecting particle collisions
A particle is associated with the corresponding cubic unit according to the spatial coordinate values of its centre. The cubic grids will be beneficial for detecting particle collisions efficiently. It is not difficult to imagine that two particles close to each other may be associated with the same cubic units or two neighbouring cubic units separately, as shown in Fig. 7a. Therefore, the algorithm for detecting particle collisions has two steps. The particles associated with the cubic unit containing the current particle (in the case of a cubic unit containing more than one particle) should be assessed to determine whether they will collide with the current particle. In addition, the particles associated with the neighbours of the cubic unit containing the current particle should be also assessed as to whether they will collide with the current particle or not. The flow sheet of a pseudo code for detecting particle collisions is presented as Fig. 8 (a) particle (b) cubic unit

Fig. 5:
The data structure of a particle and cubic unit in the threedimensional soft-sphere model of the DEM It will be helpful to take the case as shown in Fig. 7b as an example. When the 428 # particle is the current particle, the particle units associated with the same cubic unit of 628 # will be searched along the chain table after the current particle unit of 428 # first, which are the particle units 432 # , 620 # , 656 # and 697 # . Afterwards, the particle units associated with the 629 # cubic unit, the east neighbour of 628 # containing the current particle unit of 428 # , will be searched along the chain table to assess whether they collide with the 428 # particle or not, which are particle units 617 # , 751 # and 758 # .
However, it is not necessary to search the particles associated with the same cubic unit whose series numbers are less than that of the current particle, or whose positions in the chain table are in front of the current particle. As shown in Fig. 7b, particle units 620 # , 656 # and 697 # will be searched when particle 432 # is the current one. Thus, the particle unit of 428 # will not be searched because particle 432 # has already been searched when 428 # was the current particle.
The conclusion is that the current particle does not have to be judged as to whether or not it will collide with all other particles in the space of particle motion. Therefore, the computation time spent on detecting particle collisions can be saved through introducing the search algorithm based on cubic grids.

Algorithm of the three-dimensional softsphere model of the DEM
As mentioned in the discussion of the hypothesis and principle, the algorithm of the three-dimensional head->next==NULL _pball->next==NULL _pball=_pball->next; //move the pointer to the next particle node _pball->next =&ball[balli]; //link the current particle to the pointer field of the end of the chain  soft-sphere model of the DEM contains three key parts: the particle-status analysis, the particle-kinematics analysis and the search algorithm for particle collisions based on cubic grids. The algorithm flow chart is illustrated in Fig. 9, the calculation process of which can be presented as follows.
(i) Before the main program of the three-dimensional soft-sphere model of the DEM is launched, the initialization information of all particles will be scanned from an original data file, which is formatted by an initialization program or a backup of the results of the main program for the continuing calculation. The initialization information involves all the detailed records of each particle, such as the radius r p , density ρ p , temperature T p , the spatial-position coordinates x = (x 0 , x 1 , x 2 ), the translational-velocity vector v = (v 0 , v 1 , v 2 ) and the rotating angular velocity ω = (ω 0 , ω 1 , ω 2 ).
(ii) The space of particle motion will be meshed to cubic grids. As preparative work, the formation of cubic grids in the space of particle motion is indispensable to the whole algorithm of the three-dimensional soft-sphere model of the DEM.
(iii) All particles at any time step should be associated with the corresponding cubic unit of the cubic grids according to the particle spatial position, i.e. the mapping from the particles to the cubic grids should be found. The tree structure describing the association between particle units and cubic units will be constructed and it will be updated at each time step according to the new mapping between the particles and the cubic grids.
(iv) The particle collisions will be searched based on the constructed tree structure describing the association between the particles and the cubic grids. The particle unit linked to the current cubic unit will be first searched and judged whether or not collisions with current particles will occur along the chain table of particle units. The particle units linked to the neighbours of the one containing the current particle unit will be searched along the chain table of particle units.
(v) If the current particle does not collide with any other particle in the space of particle motion, the status analysis of the current particle will not be changed, i.e. the current particle-kinematics analysis will remain the same as in the last time step. Otherwise, the new status analysis will be launched and the results of dynamic analysis of the current particle will be changed according to the current new status analysis.
(vi) The iteration will be continued until the number of time steps exceeds the maximum value.
_pball->next==NULL _pball=_pball->next; //move the pointer to the next particle node _pball=&ball[balli]; //get the pointer of the current particle node of the chain  The program-development process of the threedimensional soft-sphere model of the DEM has been presented above in detail. However, the reasonability and correctness of the program, especially the particlemechanics model and the particle-thermodynamics model, should be validated against experimental data. In this research, one major goal was to confirm some key parameters in the particle-mechanics model and the particle-thermodynamics model of the threedimensional soft-sphere model of the DEM, which is useful and necessary in the application of the model to research problems.

Validation of the character of the particulate flow
In this section, the character of the particle flow simulated by the developed three-dimensional soft-sphere model of the DEM will be compared with a validation experiment of glass beads discharging from a transparent silo to validate the reasonability and correctness of the particle-mechanics model in the particle-status analysis of the three-dimensional soft-sphere model of the DEM. The experiment [37] has been carried out by González-Montellano et al., whose results will be referenced to compare with the results calculated from the developed three-dimensional soft-sphere model of the DEM program in this research. The illustration of González-Montellano's experimental equipment appears in Fig. 10. As shown in Fig. 10, the transparent silo is made of plastic (PMMA), through which the character of particulate flow can be recorded by a camera (Genie H1400-Monochrome) with a rate of 50 fps and resolution ratio of 280 × 630 pixels. The shape of the transparent silo is composed of a right quadrangular prism and a right frustum, whose size is illustrated in Fig. 10. The glass beads with a diameter of 13.8 mm were chosen as particulate materials filling the transparent silo, whose parameter values of the mechanics model listed in Table 3 are the same as reported by González-Montellano.

Filling height and filling density in the silo
About 14 000 particles were calculated according to the mass and the mean particle diameter of the glass beads. Similarly, 14 000 glass-bead units filling the same size of silo freely were simulated by the developed threedimensional soft-sphere model of the DEM program. In the validation experiment, the filling height of glass beads in the silo is ~86% of the height of the right quadrangular prism, and the filling density of the glass beads in the silo is 1423 kg m -3 [37]. The simulated results from the developed three-dimensional soft-sphere model of the DEM program show that the filling height of glass beads in the silo is 87.7% of the height of the right quadrangular prism, with an relative error of 1.98% compared to the measured value in the validation experiment, and the calculated filling density of glass beads in the silo is 1490.45 kg m -3 , with a relative error of 4.74% compared to the measured value in the validation experiment. Thus, the simulation results from the developed three-dimensional soft-sphere model of the DEM are acceptable.

Discharge time
The discharge time of glass beads T discharge from the silo can be measured by simply opening the gate at the exit of the hopper to complete the discharge process. The measured T discharge and the simulated T discharge are listed in Table 4.
Compared with the validation experiment, the T discharge calculated by the developed three-dimensional soft-sphere model of the DEM has a relative error of 3.24%. Thus, the agreement between the measured T discharge in the check experiment and the simulated T discharge in the DEM simulation is satisfactory.

Flow pattern of glass beads in the silo
The flow pattern of glass beads during the process of discharge is an important result in the validation of the mechanics model of particle-status analysis in the threedimensional soft-sphere model of the DEM program. The simulated results by the program developed in this research not only can be validated by the flow pattern recorded in the experiment, but also can be checked by commercial software, such as EDEM TM . The validation results of the flow pattern of the glass beads in discharge from the silo are shown in Fig. 11.
Each column of Fig. 11 is the validation results of the flow pattern in the silo at the same reduced time, which is defined as the ratio of real time t real to discharge time T discharge . As shown in the recorded experimental graphs in Fig. 11, three colours of the glass-beads bed layer are filled in the silo in the order of grey, white, black and white (the order of blue, white, black and white in the simulation) at the reduced time of 0T discharge . During the discharge process of glass beads in the silo, the flat bed layer of glass beads is curved in the right quadrangular prism part of silo and is curved more seriously from a U shape to a V shape in the right-frustum part of the silo because the particle velocity in the centreline of the silo is larger than that closer to the wall of the silo. It can be concluded from Fig. 11 that the flow pattern of the glass beads in the silo simulated by the developed threedimensional soft-sphere model of the DEM program is similar to that recorded in the experiment and the results simulated by the commercial EDEM TM Academic 2.3 version software.

Particle-velocity distribution in the silo
The mechanics model of particle-status analysis in the three-dimensional soft-sphere model of the DEM program also can be validated through the distribution of the velocity of the glass beads in the silo, which appears in Fig. 12.
As shown in Fig. 12a, the radial-velocity distribution from the exit of the silo 0.02 m at the time of 0.2 sec after opening the gate of the silo is similar to but a little less than the simulated results from the commercial EDEM TM Academic 2.3 version software. In addition, MFI, the ratio of the glass-bead velocity near the wall of the silo to that in the centreline of the silo along the height of the silo at the time of opening the gate of the silo, has been compared with the simulated results from the commercial EDEM TM Academic 2.3 version software, as illustrated in Fig. 12b. Therefore, the particulate flow simulated by the developed three-dimensional soft-sphere model of the DEM can be satisfactorily validated by experimentation and commercial software.

Validation of the character of particulate heat transfer
The particulate-heat-transfer model has been coupled into the developed three-dimensional DEM based on a softsphere model. However, the new coupled particulate heattransfer model and the values of its parameters can be validated by some laboratory-scale experiments. In this research, a rotary calciner built by Chaudhuri [38], as shown in Fig. 13, is adopted to validate the correctness of the particulate-heat-transfer model embedded in the developed 3D DEM program and estimate the effect of various materials on the heat transfer.
The particulate materials consist of two types of spherical particles: alumina particles and copper-alloy particles, whose parameter values in the mechanics model [39] and the thermodynamic model [38] are summarized in Table 5.

Alumina particles heated in a rotating calciner
About half the volume of the calciner was filled with the alumina particles at a room temperature of 25°C and the wall of the calciner was uniformly heated using an industrial heat gun maintained at 100°C. Meanwhile, the calciner with a rotation speed of 20 r min -1 was controlled by a step motor. The calciner was stopped at a certain time interval and the temperature of the particulate bed was tested using 10 inserted thermocouples, as shown in Fig. 4. After   the temperatures of 10 thermocouples were recorded, the thermocouples were extracted and the calciner was rotated again. The temperature profile of the particulate bed in the calciner is simulated and deduced by the new improved 3D DEM program coupled to the thermodynamics model from 0 to 12 sec, which is illustrated in Fig. 14.
It is not difficult to observe that the heat is transferred from the boundary to the centre of the particulate bed gradually. During the uniform rotation of the calciner, the alumina particles near the hot wall are heated directly and rise up to a higher position on the incline of the dynamic angle of repose. Furthermore, the alumina particles at the free surface of the inclined particulate bed are warmed up indirectly by the higher-temperature particles heated by the calciner wall rolling down from the top of the inclined bed layer to the bottom. The curve of the average temperature at the positions of the 10 thermocouples is simulated by the new improved 3D DEM program and compared with that measured by Chaudhuri, which are illustrated in Fig. 15a. The coefficient of heat transfer of this process also can be found as shown in Fig. 15b.
As shown in Fig. 15a, the average temperatures of the particulate bed at the positions of the thermocouples are increased with the rotating process turbulently because the relative standard deviation of the alumina particles in the calciner will be increased at the beginning of the heating process. The heat balance of the particulate bed layer in the calciner can be carried out and the coefficient of heat transfer can be calculated as follows: where m s is the whole mass of the aluminium-particulate bed layer in the calciner, 0.875 kg; c p,s is the heat capacity of the alumina particle, 875 J kg -1 K -1 ; T s and T s, 0 are the particulate layer temperature and the its initial value, K; T wall is wall temperature of the calciner, 398.15 K; C stands for the circle length of the calciner, 0.203π m; e s represents the ratio of the flank area of the calciner occupied by the aluminium-particulate bed layer, ~0.5; α s describes the coefficient of heat transfer of the calciner heating process, W m -2 K -1 . The calculated value of the heat-transfer coefficient is ~170.751 W m -2 K -1 , a relative error of -27.07% with the measured value by Chaudhuri of 124.527 W m -2 K -1 .
According to the results simulated by the new improved 3D DEM program coupled with the heat-transfer model compared with the measurement in the laboratory experiment, the correctness and agreement of the new coupled heat-transfer model are acceptable.

Copper-alloy particles heated in a rotating calciner
Similarly, the calciner heating process of copper-alloy particles, as another type of particulate material for validating the accuracy of the new improved 3D DEM program, is also simulated in this research. About half the volume of the calciner was filled with the copper-alloy particles at a room temperature of 25°C and the wall of the calciner was uniformly heated using an industrial heat gun maintained at 1025°C. Meanwhile, the calciner with a rotation speed of 20 r min -1 was controlled using a step motor. The average-temperature curve of the copper-alloy particles and the relative standard deviation curve of all copperalloy temperatures can be calculated and compared with those simulated by Chaudhuri, as shown in Fig. 16. It can be concluded from Fig. 16a that the average temperature of the copper-alloy particulate bed in the calciner increased with the heating time. However, the amplification of that declined gradually because of the reduced temperature difference between the copper-alloy particles and the calciner wall. Fig. 16b shows that the relative standard deviation of the temperatures of the copper-alloy particles increases with heating time from 0 to 5 sec, and is decreased with heating time after ~5 sec. The uniformity of the original copper-alloy-particle temperature is destroyed by the higher temperature of the particles near the calciner wall as compared to that of the particles in centre of the calciner particulate bed layer. Since the average temperature of the copper-alloy particles is limited to the calciner-wall temperature, the temperature difference  among all particles in the particulate layer is reduced and a new uniformity of copper-alloy-particle temperature appears at the end of the heating process in the rotating calciner.
The temperature profile of the copper-alloy particulate bed in the calciner is simulated and deduced by the new improved 3D DEM program coupled to the thermodynamics model from 0 to 20 sec, which is illustrated in Fig. 17.
The temperature profile of copper-alloy particles in the rotating calciner also can be simulated by the new improved 3D DEM program and compared with the simulation results from Chaudhuri, which are shown in Fig. 17. According to the results simulated by the threedimensional soft-sphere model of the DEM program coupled with the heat-transfer model compared with the measurement in the laboratory experiment, the correctness and the agreement of the new coupled heat-transfer model are acceptable.

Prediction of particulate behaviour of coal and sand
The reasonability and correctness of the developed threedimensional soft-sphere model of the DEM program, especially the particle-mechanics model and the particlethermodynamics model, have been validated through a series of validation experiments. However, to realize the ultimate goal of the simulation of particulate behaviour of flow and heat transfer in the process of coal pyrolysis in a downer reactor by the developed three-dimensional softsphere model of the DEM program, the physical-property parameters of coal particles as fuel and sand particles as heat carriers should be clearly validated and confirmed.
The parameter values of the mechanics model of coal and sand in the particle-status analysis of the threedimensional soft-sphere model of the DEM program for simulating the particulate behaviour of flow are listed in Table 6 [32,40].

Filling tests for coal and sand
The filling of a vertical pipe with coal and sand separately for testing the filling height and filling density of coal and sand individually has been simulated by the developed three-dimensional soft-sphere model of the DEM program and compared with the practical experiment. Experimental results for the filling height and the filling density are compared with the model simulation. The particle diameter with uniform distribution is confined to between 1 and 2 mm through screening. The vertical pipe with an inner diameter of 1.78 cm is made of quartz. Five different initial masses of coal particles and sand particles, listed in Table 7, are weighed using an electronic balance and filled in the same type of vertical quartz pipe, respectively. The particles of coal or sand in the simulation are generated by an initialization program of the three-dimensional softsphere model of the DEM simulation with the uniform distribution of particle diameters ranging from 1 to 2 mm, whose particle densities are considered as the values listed in Table 6, and the generated particle numbers of coal and sand are listed in Table 7, separately. The particle-filling state of coal and sand in the vertical quartz pipe are photographed, as shown in Fig. 18a and b, separately. Meanwhile, the simulations of the particlefilling state for coal and sand using the three-dimensional soft-sphere model of the DEM are illustrated in Fig. 18c and d. The particle-filling height shown in Fig. 18e and the particle-filling density shown in Fig. 18f have been   experimentally measured and calculated in the simulation, the distinction of which is also shown in Fig. 18e and f. The conclusions from Fig. 18e and f can be summarized as follows. (i) The maximum relative errors of the fill height in the vertical pipe between the simulation and the measurements are -9.5% for coal particles and -5.7% for sand particles, respectively. (ii) The effect of the initial filling mass on the particle-filling density in the vertical pipe is not remarkable, although the filling density is slightly increased with the initial filling mass or filling height, which can be interpreted as the pressure acting on the lower bed layer from the upper zone to reduce the bed layer voidage, further increasing the filling density. (iii) The average relative errors of the filling density in the vertical pipe between the simulated and measured values are ~6.7% for coal particles and 6.5% for sand particles. The agreement between the simulation and the experiment is satisfactory.

Voidage distribution along the bed height
The voidage distribution along the bed height can be used to check the reliability of the parameters in the mechanics model of coal and sand in the particlestatus analysis of the three-dimensional soft-sphere model of the DEM program. At the same time, the effect of the particle diameter on the profile of voidage along the bed height can be simulated using the developed three-dimensional soft-sphere model of the DEM program. In this section, 19.083 g of coal particles and 38.013 g of sand particles with five different diameters are respectively placed into five different tubes with an inner diameter of 1.78 cm. The particle numbers corresponding to the different particle diameters at one constant mass are listed in Table 8 for coal and sand separately.
The simulation results for particles of coal and sand with five diameters but with the same constant mass filling in the five vertical pipes are shown in Fig. 19.
The bed-voidage profiles along the bed height for particles with a uniform diameter from 1.0 to 2.0 mm filling in the vertical glass tube with an inner diameter of 1.78 cm for coal and sand separately are shown in Fig. 19c and f. It can be observed that the average bed voidage ranged from 0.3 to 0.4 and slightly increased with the bed height. The reason is that the lower-layer particles are compelled by the pressure from the upper bed layer, and the voidage in the lower bed is slight smaller than that in upper bed layer. It is not difficult to confirm that the voidage of spherical particles under the cubic arrangement is 0.48 and the voidage of spherical particles under the closest arrangement is 0.26. It can be concluded that the parameter values of the mechanics model of coal and sand in Table 8: The relationship between the particle diameter and the particle number with a constant mass of materials

Heaping test for coal and sand
The purpose of the particle-heaping test is to check the reasonability of the sliding-friction coefficient μ s and the rolling-friction coefficient μ r in the mechanics model of coal and sand in the particle-status analysis of the threedimensional soft-sphere model of the DEM program. The pile height and pile-base radius of coal and sand are measured and simulated in this case. The particle diameter with uniform distribution is confined to between 1 and 2 mm through screening. The particles of coal and sand in the simulation are generated by an initialization program of the three-dimensional soft-sphere model of the DEM simulation with a uniform distribution of particle diameter ranging from 1 to 2 mm, and the generated particle numbers of coal and sand that are listed in Table  7, separately. The experimental results and the simulation results for particles of coal and sand heaping on the platform are shown in Fig. 20. Fig. 20a shows the front views of the heaping state in the experiment and simulation for observing the angles of repose of coal and sand. Meanwhile, the simulations of the particle-heaping state for coal and sand using the threedimensional soft-sphere model of the DEM are illustrated  in Fig. 20b. The heaping height and heaping diameter of the coal and sand particles are chosen to describe the comparison of the simulation and the experimental results, which are presented as Fig. 20c. The heaping height and heaping diameter are measured during coal and sand heaping on the platform with five initial masses as listed in Table 7 and plotted with filled points in Fig. 20c. The relationship between the heaping height and the heaping diameter can be fixed through linear regression, which is also plotted with the solid line in Fig. 20c. The filling height is related to the filling diameter as follows: h coal,measured ≈ 0.51 · r coal,measured (32) h sand,measured ≈ 0.66 · r sand,measured The slopes of Equations (32) and (33) are presented as the tangent values of the angles of repose of coal and sand, separately. Therefore, the angles of repose are 27.07° for coal and 33.35° for sand. The heaping height and heaping diameter are simulated, and have been plotted using open symbols in Fig.  20c. The relationship between the heaping height and the heaping diameter can be fixed through linear regression, which is also plotted with a dotted line in Fig. 20c. The filling height is related to the filling diameter as follows: h coal,simulated ≈ 0.45 · r coal,simulated (34) h sand,simulated ≈ 0.68 · r sand,simulated (35) Similarly, the simulated angle of repose of the coal pile is 24.17° with a relative error of -10.71% compared with the experimental value, and the angle of repose of the sand pile in the simulation is 34.17° with a relative error of 2.46% compared with the experimental value.

Hourglass using coal and sand for timing
Particles can be filled into an hourglass for timing. However, the discharge time of the same volume of particles with the same uniform diameter distribution is determined by the physical properties of the solid particles. Therefore, the discharge times of the same volume of coal and sand with the same particle size from the same hourglass must be different. 13 500 coal particles (about 70.65 g) and 13 500 sand particles (about 138.47 g) with the same diameter of 2 mm have been filled into the hopper of the hourglass to measure the time scale of the two hourglasses, which has been also simulated by the developed three-dimensional soft-sphere model of the DEM program. The simulation results of the discharge process of coal and sand from the hopper of the hourglass are shown in Fig. 21.
The time scales of the two hourglasses containing coal and sand separately have been compared between the experiment and simulation, respectively, which has been listed in Table 9.
It can be concluded from Fig. 21 and Table 9 that the discharge time of the coal hourglass is a little shorter than that of the sand hourglass. The relative errors of the simulated results of coal hourglass and sand hourglass are 13.59% and 12.68% compared with the experimental results.
According to the conclusions from the discussion content in the research, the developed three-dimensional soft-sphere model of the DEM program can be capable of predicting the particulate behaviour of static state and flow for coal and sand.

Conclusions
A three-dimensional soft-sphere model of the DEM for simulating the particulate behaviour of flow and heat transfer is developed, the reasonable hypothesis and the ingenious algorithm of which have been presented in as much detail as possible. The organizational structure of the developed three-dimensional particulate soft-sphere model of the DEM contains algorithms for determining particle collisions, for the status analysis of particles and for the kinematics analysis of particles. The search algorithm for particle collisions is based on cubic grids and the splicing of cubic grids can originally solve the problem of searching for particle collisions in a particulate system. The particulate-status analysis is used for solving the causes of particle motion and heat transfer, which involves the mechanics model and the thermodynamics model. The particulate-kinematics analysis is described by some basic laws of physics containing Newton's second law for particle translational motion, the rigid-body law of rotation for rotating particles and the heat-transfer law for two particles with different temperatures.
The correctness and reasonability of the developed three-dimensional soft-sphere model of the DEM program has been validated and checked through simulating the conditions compared with those in the experiment. In this research, the mechanics model in particle-status analysis and the heat-transfer model in particle-kinematics analysis of the three-dimensional soft-sphere model of the DEM have been validated by the experiment of glass-bead discharge from a silo and metal particles heating up in a calciner. The agreement between the experimental and simulated results is consistent and satisfactory, which can indicate that the developed three-dimensional softsphere model of the DEM program is successful and capable of predicting the flow and heat transfer in particulate systems.
In addition, to prepare the simulation of the coalpyrolysis process in a downer reactor using the developed three-dimensional soft-sphere model of the DEM program in the next step, some necessary physical-property parameters of the fuel particles of coal and the heat-carrier particles of sand have been validated and confirmed through comparing the experiments of particle filling in the tube and particle heaping on the platform. The agreement between the simulated results and tested data is satisfactory. Therefore, it can be concluded that the parameter values of the status-analysis model and the kinematics-analysis model for coal and sand particles are acceptable for the next step of simulation of the coal-pyrolysis process in a downer reactor using the developed three-dimensional soft-sphere model of the DEM program.

Scalar:
A (p) = particle surface area, m 2 A (p),c = particle contact area, m 2 A gap = area of tiny gas interval between two contact surfaces, m 2  c p = specific heat capacity, J kg -1 K -1 d gap = width of tiny gas interval between two contact surfaces, m d p = particle diameter, m I (p) = particle rotational inertia, kg m 2 l 12 = distance of two particles, m l a = length of side a of a cubic grid unit, m l b = length of side b of a cubic grid unit, m l c = length of side c of a cubic grid unit, m l pw = distance between particle centre and wall, m E = elastic modulus, Pa e = recovery coefficient, − r (p) = particle radius, m m (p) = particle mass, kg m* = reduced mass of two particles, kg k g = thermal-conduction coefficient of gas phase, W m -1 K -1 k n = normal stiffness, N m -3/2 k t = tangential stiffness, N-m -1 H (p) = particle enthalpy, J Q = thermal quantity, W Q con(,p) = quantity of particulate thermal conduction, W r p = particle radius, m R g = thermal resistance in gas phase, K W -1 R s = thermal resistance in solid particle, K W -1 R total = total thermal resistance, K W -1 T = Kelvin temperature, K T p = particle temperature, K α c = contact coefficient, − δ n = normal overlap, m δ t = tangential overlap, m ε = voidage of bed, − η n = damping coefficient, kg m -1 η t = tangential damping coefficient, kg m -1 μ g = dynamic viscosity of gas phase, Pa s μ r = rolling-friction coefficient, − μ s = sliding-friction coefficient, − ν = Poisson's ratio, − ρ p = particle density, kg m -3 Δl = surface distance of two particles in the contact zone, m Δt = time step, 10 -5 sec ΔT = temperature difference, K Vector: e n = normal-unit vector e t = tangential-unit vector F = net-force vector, N F n = normal-force vector, N F cn = normal-contact-force vector, N F dn = normal-damping-force vector, N F t = tangential-force vector, N F ct = tangential-contact-force vector, N F dt = tangential-damping-force vector, N g = gravity acceleration, m s -2 M = total-moment vector, N m M t = moment of tangential contact force, N m M r = moment of rotation damping, N m v = translational-movement-velocity vector, m s -1 v n = particle normal translational velocity, m s -1 v t = particle tangential translational velocity, m s -1 x = space-position vector, m ϕ = phase-position vector, rad ω = rotating-movement-velocity vector, rad s -1