Rotation and orientation of irregular particles in viscous fluids using the gradient smoothed method (GSM)

ABSTRACT The dynamics (rotation and stabilization) of single irregular particles that are constrained to rotate only and until to a quasi-steady state in a two-dimensional channel with viscous flow is numerically investigated using the gradient smoothed method (GSM). The GSM is proved to be computationally stable for arbitrary, irregular geometries discretized with distorted grids and well agreements with others’ work are revealed in validation cases. This work focuses on the influences of irregular particles’ typical shapes, including ellipse, rectangle and triangle, on the relevant surface forces, the flow dynamics, and the response time of the rotation. The effects of aspect ratio and flow velocity are studied in detail for all three typical types of irregular particles. It is found that the imbalance of the total torque on the surface of the particle causes the rotation, and when the particle is approaching the final stable position, the total torque becomes nearly zero with a small fluctuation, which contributes a local oscillation around the stable direction. In our cases, under the constraints, it is found also that the broad side of the elliptical particle always tends to be along the stream. For the rectangular particle, however, the aspect ratio and the flow velocity collectively determine the final orientation which means the major side or the diagonal line is along the flow stream. In addition, the triangular particle is found to behave quite differently in terms of both rotation and stabilization. The ‘response time’ of three types of particles is finally obtained from our GSM simulations. These findings could be helpful for a better understanding on the fluid-particle interactions and maybe advisory for determining the shape factor which is a key parameter for multiple particles motions.


Introduction
Transportations of particles are involved in many natural situations and industrial fields, such as transportation of red blood cell with blood flow in medical science, carrying of drilling cuttings in petroleum drilling industry, scouring of bottom sediment around structures on seafloor (Xia et al., 2009). This work is initially inspired from the cuttings transportation problems in oil drilling process. Figure 1 shows the cuttings from a drilling site, including cuttings shapes and measurement on grain size (Huang, Chen, & Dong, 2008). Excessive simplifications on attributes of cuttings particles, especially on the shapes, make some computational results less trustable. From past investigations and according to the shapes most commonly found in two dimensions, ellipse, rectangle and triangle are believed to be basic. However, above three basic shapes are idiomatically assumed as a circle in engineering simulation and the behavior differences in viscous flow due to the irregular shapes are barely CONTACT Bing Shao upcshaobing@163.com concerned. Our work is investigating the rotation behavior of these three basic shapes of particles in water flow numerically. Since two decades ago, researchers have done lots of work on dynamics of particles' freefalling in fluids both experimentally and numerically (Ardekani, Dabiri, & Rangel, 2008;Feng & Michaelides, 2004;Goyal & Derksen, 2012;Hu, 1996;Huang, Hu, & Joseph, 1998;Pu, Shao, Huang, & Hussain, 2013). Among the numerical methods, finite element method (FEM) (Jotsa & Pennati, 2015;Liu & Quek, 2013), finite difference method (FDM), finite volume method (FVM), and even meshfree methods (Liu, 2009;Pu et al., 2013) are most widely used. Hu, Joseph, and Crochet (1992) used FEM to model the Navier-Stokes equations for liquid and the Newton's equations of motion for solids. They revealed the effect of vortex shedding on the motion of cylinders and classified the falling process as drafting, kissing, and tumbling scenario. Hu, Fortes, and Joseph (1992) directly simulated the falling of an ellipse in a 2D vertical channel, and the pressure at the front stagnation point on the surface of the ellipse was studied to explain the turning mechanism. Huang, Feng, and Joseph (1994) conducted a more detailed study on the turning couples and the stagnation points. They found that, the high pressure at the stagnation point which was identified with the point where the shear stress vanishes, turns the ellipse broadside to be horizontal and the ellipse oscillates due to the vortex shedding. They also did lots of work on the viscous torque and the pressure torque, concluding the pressure torque which is much larger, dominates the ellipse's turning.
Some researchers used the lattice Boltzmann method (LBM) and the immersed boundary method (IBM) to investigate the dynamics about the settling of solids in fluids (Haghani, Rahimian, & Taghilou, 2013;Lee et al., 2008;Peskin, 2002). Lai and Peskin (2000) used IBM to simulate the flow past a circular cylinder and studied the effect of numerical viscosity on the accuracy of the computation. Using a combined immersed boundary-lattice Boltzmann method (IB-LBM), Feng and Michaelides (2004) simulated the sedimentation of a large number of circular particles in an enclosure filled with fluid. More recently, Xia et al. (2009) studied the dynamics of a single two-dimensional elliptical particle settling in a Newtonian fluid using LBM. In Xia's study, the effect of boundaries on the flow patterns during sedimentation was emphasized and the effects of density ratio, aspect ratio and the channel blockage ratio were studied. Based on the IBM and the FEM, Zhang, Liu, and Khoo (2012) proposed a novel method called immersed smoothed FEM (IS-FEM). The IS-FEM can solve the two-dimensional and three-dimensional fluid-structure problems with largely deformable nonlinear solids in viscous flow (Zeng, Liu, Li, & Dong, 2016). Settlings of both single and multiple particles in viscous liquid are investigated as examples by IS-FEM and it shows good stability and wide usability. Other numerical methods, like the discrete element method (DEM), could also be a well tool to study the particle and fluid interactions and the DEM has exhibited its enormous advantage on the granular systems (Lu, Third, & Muller, 2015;Zhong, Yu, Liu, Tong, & Zhang, 2016).
From the literature, work on the dynamics of particles in fluid mainly lies in falling of a chosen particle and flow around a fixed cylinder. Rotation of a single particle is seldom stripped from the motion and the attempt in our work may bring new thinking from another way to the interactions of irregular particles and fluid. Most recently, Sanjeevi, Padding, and Kuipers (2015) and Zastawny, Mallouppas, Zhao, and Wachem (2012) derived a framework to predict the drag, lift and torque coefficients for non-spherical particles. The effects from Reynolds number, the rotation, and the incidence angle are investigated. Relevant data from their work is adopted to validate our results.
The numerical method used in this paper is called the gradient smoothing method (GSM), enlightened by the attractive merits of gradient smoothing operation in Galerkin weak forms (Wang, Khoo, Liu, Xu, & Chen, 2014;Xu, Li, Tan, & Liu, 2012). However, it is a gradient smoothing operation applied for strong-form (differential form) of governing equations. The main line of this method is, all the unknowns are stored at nodes and the gradient smoothing operation is used to approximate the derivatives based on relevant gradient smoothing domains (GSDs). In GSM, both regular and irregular grids can be used with ease through properly designed types of smoothing domains, which makes the GSM applicable for arbitrary geometries with excellent stability.
In the rest of this article, first, the theory of GSM is briefly presented in Section 2. In order to get the turning couples on the particle, we adopt the method used by Huang et al. (1994) for the ellipse and extend the method to the rectangle and triangle. The forces and fluid field are calculated and compared among particles with different shapes. Further, the dynamic rotations of irregular particles are revealed and compared. Finally, the effects of aspect ratio and flow velocity on the rotation are studied.

Briefing on the gradient smoothing method
The key idea in the GSM theory is the approximation of derivatives of approximated functions using various types of properly designed smoothing domains. Here, we brief the process of approximating the first and second order spatial derivatives. The gradient smoothing operation can be applied to approximate derivatives at nodes, centroids of cells and midpoints of cell edges over relevant gradient smoothing domains (GSDs) that can be node-based, cell-based, and edge based (Liu, 2009).

Governing equations
The flow in the fluid domain t with density ρ f and viscosity μ f is governed by the Navier-Stokes equations: (1) where g denotes the gravity vector, u denotes the velocity, t is the time and σ is the stress tensor given by The governing equation for the motion of the particle is Newton's Law: where m denotes the mass of the particle and I is the moment of the inertia. F is the total force exerted from the fluid on the particle and F m is the total torque. U(t) = dX(t)/dt is the particle velocity and (t) denotes the angular velocity.

Gradient smoothing operation
We now use a 2D problem to illustrate the gradient smoothing operation for approximating the gradients of a function. For any point x i in domain i , the gradients of a field variable U can be approximated in the form of convolution (Liu & Liu, 2003;Liu & Xu, 2008;Lucy, 1977). where ∇ is a gradient operator and ω is the smoothing function. Equation (4) can be integrated by parts as in which ∂ i denotes the external boundary of the GSD and n is the unit outward normal vector on ∂ i , as shown in Figure 2. For simplicity, the following piecewise constant function is used as the smoothing function: where V i is the area of i . With Equation (6), the second term on the right side of Equation (5) vanishes, leading to: Equation (7) offers a more robust way to approximate the gradient because numerical integral along the boundary is more cost-effective than over the surface of the smoothing domain (Wang, Khoo, Liu, & Xu, 2013). In addition, Equation (7) can be further extended for functions that may not be continuous (Liu, 2008). Analogously, applying the gradient smoothing technique to the first-order gradient over another smoothing domain gives the approximated second-order gradient (Liu, 2009;Xu, 2009), and the Laplace operator at x i can be approximated as Using Equation (7) and Equation (8), the first and second order spatial derivatives at any point of interest can be approximated with properly defined smoothing domains which will be presented in next part.

Smoothing domains
In GSM, the values of the field functions are stored at nodes and by the connecting nodes, the problem domain is divided into a set of regular or irregular cells. Based on primitive cells, a smooth domain can be constructed for any point. As shown in Figure 3, smoothing domains are constituted on primitive unstructured triangular cells, which can be automatically generated. Three different types of smoothing domains are used for approximating spatial derivatives depending on the location of interested point, which includes the node-associated GSD (nGSD), the centroid-associated GSD (cGSD) and the midpoint-associated GSD (mGSD).

Approximation to spatial derivatives
Multiple discretization schemes for spatial differential terms have been developed using different types of quadrature and methods of approximation at any locations including node, midpoint and centroid. Considering the balance of numerical accuracy and computational efficiency, the one-point quadrature (rectangular rule) is chosen. The approximation of the first-and second-order derivations of a nodal field variable is detailed in this section and further implemented in following simulations. The other schemes and analysis of stencils can be found from Xu (2009) and Liu and Xu (2008).

First order derivatives at nodes
By Equation (7), at node i, the first order derivatives of the field variable U are given by where U m denotes the value of the field variable U at midpoints of cell edges; S x and S y are the two components of a domain-edge vector; n x and n y are the two components of the unit normal vector of the domain edge, n = n x i + n y j; j k denotes the kth node surrounding the node i; M ij k denotes the midpoint of the cell-edge of interest, ij k ; C ij k j k+1 and C ij k−1 j k represent the centroids of two triangular cells on the sides of cell edge ij k ; n i denotes the total number of supporting nodes connected to node i.

Second order derivatives at nodes
In the one-point quadrature schemes of GSM, the second-order derivatives in Laplace operator can be approximated by Equation (8) and given as: As shown in Equation (11), only the values of the field variable and its gradients at the midpoint of cell edges are needed to complete the approximations. It should be noted that these two gradients are approximated by the related mGSD associating with the edge ij k in 2.4.1. The vectors for edges connected with the cell edge ij k can be lumped together, which in return reduces the storage space for geometrical parameters. Here, where A M denotes the area of the mGSD. The domain edge vectors S x M and S y M for the mGSD are calculated as The field variable U C at the centroids of the cells is obtained from interpolation of the function values at related nodes, as follows: It is necessary to note that, the principle of gradient smoothing operation is not limited to the triangular mesh in 2D domain and it can also be easily extended to the hybrid mesh both in 2D and 3D cases.

Calculation of pressure and stress torques on the particle surface
The stress tensor and pressure at every node on the particle surface are needed to calculate the viscous part of the normal stress and the shear stress. To get the torques due to the pressure and the viscous stress, the method proposed by Huang et al. (1994) is adopted here as follows.
In Figure 4, the node of interest on the surface of the particle is displayed and the elements around the nodes as well as the geometrical parameters are defined.
In Figure 5, the x-axis and y-axis are along the walls of the channel. (X i , Y i ) is the center of the ellipse and k is the node on the surface. θ is the turning angle of the ellipse. r is the distance between node k and the center. a 1 and b 1 are the arms of the forces F n and F t .
In Figure 5, the transformation of coordinates from x, y to x , y is given by  The equation for the ellipse surface is F(x , y ) = (x 2 /a 2 ) + (y 2 /b 2 ) − 1 = 0, where a, b are the semilengths of the major and minor axis respectively. The unit outer normal vector is n = (cos α, sin α) and the unit tangential vector is t = (− sin α, cos α) while in x, y system, it should be like n = (cos β, sin β) and t = (− sin β, cos β), where β = α + θ − (π/2). For the node k at (x k , y k ), the angle α can be calculated by The viscous force can be got from: where T is the stress tensor for node k and τ xx , τ xy , τ yy are the components; A is the area of the region around node k, which is composed of one half of the two segments on both sides: Then the shear force can be given by The pressure torque and viscosity torque can be got from: where By Equations (14)-(20), the torques due to both the pressure and the viscous stress can be obtained and analogously, we can calculate the torques on the surface of the rectangular particle and triangular particle by slight modification.
For the rectangular particle, a and b turn to the semi-lengths of the major and minor side. As shown in Figure 6, some angles need to be modified and are given as: , when the interested node is not at the end of any side, when the interested node is at the end of side, The triangle used here is an isosceles triangle and the semi-length of the two equal sides is defined as a and the semi-length of the base is b. Unlike the rectangle, the base angle δ is calculated by δ = arccos (b/2a). Some angles in Figure 7 need to be replaced with: , when the interested node is not at the end of any side, when the interested node is at point A, where h denotes the height of the triangle, h = √ 4a 2 − b 2 ; A, B and C denote three vertices as shown in Figure 7.
With Equation (21) and (22) substituted into Equation (14) ∼ Equation (20), the surface torques for the rectangular and triangular particles can be obtained.

Convergence study
The validation of GSM has been done on numerical errors and robustness to the irregularity of the cells by Liu and Xu (2008). From what they have done, even with highly distorted triangular cells, stable and accurate results still can be obtained. For the simulation in this paper, with the rotation of the particle, the grids around the particle become distorted to some extent. To investigate the numerical results by the GSM, a numerical error indicator is defined using the L 2 -norm of error as follows: where U i andÛ i respectively are predicted and analytical solutions at node i. This type of error is used to compare the accuracy between the GSM and FVM here. The convergence rate with the reducing average nodal spacing (h) is plotted in Figure 10. h is defined as follows: where A and n node denote the area of the whole computational domain and the total number of nodes in the domain, respectively. From Figure 8, the convergence rate of the GSM is higher than that of the FVM. Therefore, considering the distortion of the cells in the process of rotation and the accuracy, for our simulation, the GSM is believed to be more suitable to be adopted.  Zastawny et al. (2012) did the recent work on the prediction of the drag and lift coefficients for non-spherical particles in a flow. Prior to the simulation on the rotation of irregular, non-spherical particle, the prolate ellipsoid with varying incidence angle and the relevant data are used to validate the method adopted in present work. Figure 9 shows the geometry and the boundary conditions. The fluid is water with density 1.0 g/cm 3 and viscosity 0.01 cm 2 /s. The velocity of inlet flow U ∞ is adjusted to vary the Reynolds number and therein, the characteristic dimension used is the equivalent particle diameter, which is the diameter of a circle owning the same area with the ellipse. The dependency of the drag coefficient on the Reynolds number is shown in Figure 10 and is compared with the work of others. It shows the results follow the trends of Hölzer and Sommerfeld (2008), Zastawny et al. (2012) and Sanjeevi et al. (2015). However, for incidence angle ϕ = 90°and Re > 30, the result is comparatively lager than that from literature in a low   degree. Figure 11 presents a well agreement of current results with Zastawny et al. (2012) at low Reynolds numbers (Re = 10). However, at high Reynolds numbers (Re = 300), it tends to over-predict the drag. Acceptable deviation is also shown in Figure 12, where the torque coefficients are compared with the work of Zastawny et al. (2012).

Computing conditions
In our case, restrictions are applied to isolate the rotation from the movements. The displacements of particles are set to zero, which means the particles can only rotate under the flow around their centroids. In order to compare the rotations among particles with different shapes and sizes, the moments of inertia of all particles are consistent to the value of the elliptical particle with a 0.5 aspect ratio where the physical lengths of semiminor axis and semi-major axis are 0.025 cm and 0.05 cm respectively, I = (m/4)(a 2 + b 2 ). The aspect ratio here is defined as the ratio between the minor and major axis for the elliptical particle while for the rectangular particle, it is the one between height and width. And for the isosceles triangular particle, the aspect ratio is the ratio between the length of the base line and the length of the other two equal sides. The effects from both aspect ratio Moreover, the density of the particle is 2.6 g/cm 3 and the fluid is water with density 1.0 g/cm 3 and viscosity 0.01 cm 2 /s. The simulation domain is referred originally from Hu (1996) and some alterations are made. The inlet of the simulation domain is placed 25a behind of the particle and the outlet is placed 25a ahead of the particle, as illustrated by Figure 13. The particle is put at the center of the domain with a width of 6a and both the particle and the fluid are set to be still initially.

Forces on the particle surface during the rotation
In this section, the pressure and shear stress on the surface of the particle are computed and discussed. In Figure 14, dividing streamlines around the particle and the separation points indicated by the arrows are shown. The point where the shear stress vanishes is so called 'stagnation point', which will be further discussed in next subsections. Cases in this part are computed with the flow velocity of 0.2 m/s and aspect ratios of three irregular particles of 0.5. Six moments are selected to present the variation of pressure and shear stress on the surface of the Figure 14. Dividing streamlines around the particles and the L for the x-coordinate used in following figures. L denotes the percentage of the distance from a chosen node in the particle's circumference.
elliptical particle during the whole process of rotation. The distribution of pressure at six moments is shown in Figure 15. The shear stress and the pressure are non-dimensionalized by (1/2)ρU ∞ 2 . At the very beginning, the pressure is positive on the front surface and negative on the rear. But when it starts to roll (rotation angle < 45°), the positive pressure zone diverts. Due to the rotation and vertex shedding, the positive area shrinks and the negative one extends, as shown by the curves of 0.113π and 0.227π in Figure 15(a). As the rotation continues, the particle turns its broad side to be along the stream and the area of the front face decreases to a minimum value. At same time, the area of positive pressure on the particle surface extends until the pressure on the whole surface becomes positive. Analogously, six moments are chosen for the rectangular particle. Due to the sharp corners of the rectangular particle's surface, there exist break values of the pressure and the shear stress at the corners, as shown in Figure 16. Except the break points, the whole distribution of pressure and its variation with time are similar to the elliptical particle,  but with a lower value. Compared with the elliptical particle, the sharp change of pressure causes more energetic vortexes along the rectangular particle's surface. Figure 17 shows the pressure distribution on the triangular particle's surface. The difference from the other two particles is, most of the time the surface pressure on the triangular particle are positive. The only exception exists at the beginning and certain locations like the minor side, seen from Figure 17(a), or the top vertex, seen from Figure 17(b).
The change of the shear stress on the ellipse's surface is demonstrated in Figure 18. At the initial stage, compared with Figure 15, the point where the shear stress goes through zero, known as the stagnation point, is also where the pressure reaches the maximum. However, it becomes no longer obvious with the ongoing rotation and the unstable vortex shedding. The maximum values of  shear stress occur near the two ends of the major axis when the particle starts to rotate, which can be illustrated by the sharp turning streamlines around the ends. But when the rotation angle grows, as shown in Figure 18(b), the maximum value appears only near the front end. Since the value of shear stress is much smaller than that of pressure, it is not the shear stress that mainly makes the particle rotate but instead, it is the pressure.
For the rectangular particle, there is an obvious feature, namely, piecewise distribution of shear stress. Although the pressure of rectangular particle is smaller than the elliptical particle somewhat, the shear stress is relatively larger, as shown by Figure 19. From Figure 21, it can be visually seen that the red section of the contour with a high value around the rectangular particle is more extensive and that is also what generates a larger shear stress. Additionally, different from the elliptical particle, Figure 19. The distribution of shear stress on the surface of the rectangular particle. the straight outline of the rectangular particle meeting the stream produces a sharper change of the shear stress, as shown in Figure 19(b). The characteristic of the piecewise distribution of the shear stress for the triangular particle is more obvious and the peak values exist on the three vertexes, as shown in Figure 20. The streamlines around the particle during the rotation can be seen from Figure 21, where also the contour of flow velocity and the vortexes on the particle surface are revealed. The particles in a water flow with the centroids fixed, behave like a windmill in the wind. Results show that the elliptical particle is finally stabilized with its broad side along the stream (as shown by the last snapshot), but small oscillation around the stable position is always there. Compared with the elliptical particle, the flow around the rotating rectangular particle is speedier and the vortexes around the rectangular one are bigger at same locations, which gives rise to a more aggressive oscillation. With identical simulation conditions, the rotation of the triangular particle brings about a more violent disturbance to the stream which is due to the faster rotation of the particle.
The total torques for the three irregular particles caused by the pressure and the shear stress are shown in Figure 22. The torques are normalized by T g = (ρ p − ρ)A p ga, where ρ p and A p denote the density and the area of the particle, respectively. At the initial position, both the elliptical and the rectangular particle occupy a relatively high positive torque, creating a quick response to the flow. But the torques drop immediately as the particles rotate, which is trying to slow down the rotation. At time step of 20, the total torques of the elliptical and rectangular particles decrease to zero and the particles reach the stable position and start oscillating. However, it is totally different for the triangular particle. Before time step of 40, the triangular particle maintains a well balance on the torque and keeps itself in a tiny rotation angle. But after, the vortex shedding breaks this balance and meanwhile, the total torque rises and reaches a maximum value at time step of 50 and then begins its decrease. The rotation is slowing down and at time step of 80, the triangular particle arrives at the stable position and starts to oscillate. More detailed discussion is continued in the next section.

The rotation and stabilization of the irregular particles
In this section, the rotation and the stabilization laws of these three particles are presented and the effects from the aspect ratio and flow velocity on the response time are discussed. The response time defined here is the time range from the particular time step after which the angular velocity initially breaks zero, to another time    step when the angular velocity returns to zero first. The response time represents the first numerical rotation cycle and is intended to reflect the reaction sensitivity of the particle to the flow.

Elliptical particle
In order to observe the effect from flow velocity on the rotation, a series of flow velocities are used, including U i /U im 1.0, 0.75, 0.5, 0.25 and 0.1 (U im = 0.2 m/s). The elliptical particle is placed initially with its major axis perpendicular to the flowing direction, as shown by the lower icon in Figure 23. It can be easily seen from Figure 23 that, when the stream flows faster, less time is needed for the inside particle to reach the stable position. The particle rotates anticlockwise and finally stabilizes with its major axis along the stream, but the wave-shaped curve indicates a continuous oscillation around the finial position. A set of aspect ratios are used to study its effect and the result curves are shown in Figure 24, where the dash lines denote certain rotation angles indicated by the icons. It can be seen that the aspect ratio does not affect the final stable direction. However, with a smaller aspect ratio which means a slender body, less time is needed to reach the stable position. On the contrary, a more circular particle needs more time to be stabilized. The combination curves can be found in Figure 25.
By Figure 25, for the response time, there is a power function law with the flow velocity while a linear relationship with the aspect ratio. Further comparisons are made in sections below.

Rectangular particle
Rotation curves of the rectangular particle with the aspect ratio of 1/4 under a set of flow velocities are shown  in Figure 26. The three dash lines with icons in each subgraph represent three positions: the initial position with the major side perpendicular to the flow direction, the middle position with the diagonal line along the flow direction and the other position with the major side along the flow direction. Obviously, in a faster flowing fluid which means a higher Reynold number, the rotation becomes quicker. The rectangular particle in a 1/4 aspect ratio is going to stabilize with its long edge parallel to the flowing direction. But some difference can be obviously found in Figure 27.
By Figure 27, it is clearer that the stable position of the rectangular particle is bound up with its aspect ratio. When the particle is more slender with a smaller aspect ratio, it tends to stabilize in a direction with its long edge along the stream. Oppositely, with a larger aspect ratio, the particle is more like a square and it is more likely to turn a smaller angle. Moreover, this stable direction is also related to the flow velocity. A faster stream can make the rectangular particle rotate a bit more from where its long edge is along the stream to where its diagonal line is.
Similarly, here are also given the relation curves between the response time with both the flow velocity and aspect ratio. Same to the elliptical particle, a power function law exists between the response time and flow velocity. The relationship between the response time and the aspect ratio is affected by the stable position. With the aspect ratio smaller than 0.5, the response time changes not much due to the common final direction. Nevertheless, when the aspect ratio is over 0.5, as shown in Figure 28(b), a drop of the response time occurs due to the shift of the stable direction.

Triangular particle
As what is done for the other two particles, Figure 29 shows the rotation curves of the triangular particle in stream with a set of velocities. The dash lines and icons are denoting the stable angle (or stable direction) and rotating direction respectively. Obviously, a faster flow prompts the triangular particle to rotate quicker to reach the stable position. The particle with aspect ratio of 1/4 starts to turn firstly anticlockwise until to the angle of 3π/4 and then stops to turn oppositely. Finally it rotates clockwise to the position where its major height is along the flow with the vertex meeting the stream. More about this can be found in Figure 30.
By Figure 30, when the aspect ratio is under 2/5, the triangular particle possesses another kind of action. That is, the particle firstly turns anticlockwise for 3π /4 and then stops to turn oppositely until its major height is parallel to the flowing direction and the vertex is meeting the stream. A local oscillation finally takes over without any effect from the flow velocity. It is the total torque on the surface of the particle that produces this, seen in Figure 31. A more slender triangular particle (aspect ratio of 1/3 in Figure 31) owns a larger total torque, which gives the particle a quicker response.
Similarly, the response time of the triangular particle is in a power function law with the flow velocity and has a linear relation with the aspect ratio, as shown in Figure 32.

Comparison among irregular particles
In this section, the rotations of three irregular particles are compared and the effecting law from the shapes is expected.
Curves in Figure 33 show the elliptical and rectangular particles share something in common. Quick responses are made both by the elliptical and the rectangular particles meaning that they start to rotate quickly in the stream. But less time is taken by the elliptical particle to    reach the stable position. Nevertheless, it is different for the triangular particle as its response is much more lagging. With the flow velocity of 0.25U im , there is almost no reaction at time step of 100. Referring to Figure 22, the initial balance of total torque is going to be broken at time step of 40. However, once the triangle particle starts to roll, it rotates much faster and reaches the stable position at almost same time with the other two. The reason is the positive total torque on the surface of the triangular particle lasts longer, which provides a longer rotating acceleration.
From Figure 25, Figure 28 and Figure 32, the response time always possesses a power function law with the flow velocity and a linear relationship with the aspect ratio. So, an equation can be derived from the fitted curves as T = (cr/V), where the T denotes the response time, r is the aspect ratio, V is the flow velocity and c is a constant.   Combining the time ranges covered by the response time curves of three irregular particles, we get Figure 34. The green region, red region and the purple region denote the response time coverages for the rectangular particle, the elliptical particle and the triangular particle respectively. In Figure 34, the green region and the purple region, which denote the response times of the rectangular particle and the elliptical particle respectively, visibly overlap in a large scale. Additionally, the red region is completely covered by the green one. As a whole, the upper boundary is owned by the rectangular particle while the lower boundary is hold by the elliptical particle. From what we got, for the simulation with massive particles, the shape of the particles can be simplistically considered as rectangle or ellipse.

Conclusion
In this paper, the GSM is adopted to study the rotation of the irregular particles in a channel filled with moving viscous fluid. Three irregular particles are studied, including the elliptical particle, the rectangular particle and the triangular particle.
At the very beginning of rotation, the total torque on the surface is relatively much bigger, but with the ongoing rotation, the total torque decreases. When the particle is approaching the final stable position, the total torque is nearly zero with a small fluctuation, which contributes a local oscillation around the stable direction in a small scale.
The forces for the elliptical particle are close to the results found by other researchers and something novel is obtained. The much larger value of pressure is the main factor making the particle rotate, rather than the shear stress; in our cases, the broad side of the elliptical particle is finally turned to be along the stream. For the rectangular particle, the final stable direction is affected by its own aspect ratio and the flow velocity. With a more slender body and a slower flow, the major edge of the rectangular particle tends to be parallel to the flowing direction. Otherwise, the rectangular particle has more tendency to turn its diagonal line to be along the flow direction. Different from the other two, the triangular particle reveals a lagging response in the same stream. The top end of the triangular particle can disturb more vortexes during the rotation and finally the particle stabilizes with its height parallel to the flowing direction and with the top end meeting the stream.
The response time ranges of three particles are combined and the inclusion relation is illustrated. The work reveals the differences between irregular particles on their rotation in viscous flow. What we got from the work in this paper could be of a reference value to the simulation concerning massive particles where some simplifications on the particle shapes are necessary. But some limitations still exist for our work. The rotation of the particles is studied independently from the comprehensive motion which is most often happening in real cases and moreover, some parameters used by the author maybe not very appropriate. So, further study is being conducted on the actions of irregular particles in the viscous flow and more accurate results are being produced.