Lagrangian tracking of fibres in a channel flow

Tracking of fibres in a fluid flow is much more complicated than tracking of spherical particles. In fibre motion, the orientation of fibre against the flow direction plays a very important role. In addition to the standard equation of motion, additional equations for orientation and angular velocity must be solved during the tracking of fibres. A mathematical model describing fibre motion is introduced in this work. Capabilities of this model are demonstrated through simulations of fibre transportation by air in a channel flow. The importance of the terms in the equation of angular velocity are discussed.


Introduction
Microscopic mineral fibres (asbestos) may cause nonmalignant pleural disease, asbestosis, lung cancer and mesothelioma. Man-made vitreous fibres are less dangerous; however, their health effects are poorly understood. One of the first steps of the correct estimation of health risks associated with the inhalation of fibres is to quantify the amount of fibres deposited within the airways. The amount of fibres depositing in the airways depends on many factors, such as the number concentration of fibres in the ambient air, their physical properties (size, aspect ratio, density), type (oral or nasal) and intensity of breathing (depending on the physical activity), among others.
A relatively low number of investigators have reported results of experimental fibre deposition measurements in airway replicas (e.g. Myojo and Takaya, 2001;Su and Cheng, 2009;Cheng, 2007;Belka et al, 2018). Numerical simulations of fibre motion were also performed in the upper (e.g. Inthawong et al, 2013) and bronchial airways (e.g. Feng and Kleinstreuer, 2013) as well as in circular ducts (Tian et al, 2012) and channels (Lo Iacono et al, 2005).
Fibres exhibit more complex aerodynamic behaviour than spherical particles, thus tracking of fibres is a challenging task. The most important difference compared to spherical particles is that as a result of their interaction with the airstreams fibres undergo not only translational, but also rotational motion. Knowledge of fibre orientation is important because the forces which determine the fate of fibre particles within the airways (e.g. drag force) are functions of fibre orientation relative to the direction of the airflow. Many factors may affect the preferential or random orientation of fibres, such as the characteristics of the airflow, the Reynolds number associated with the particle, initial orientation of the fibres, the geometry of the airways or the aspect ratio of the fibres. The aim of this work was to develop and apply numerical techniques in order to track the fibres in channel flow and analyse their time dependent orientation for different initial positions and fibre sizes. The effect of these parameters on fibre behaviour will be discussed. An outlook of future developments will also be provided.

Model of fibre movement
In this study fibres were tracked by using a Lagrangian approach. Superposition principle was applied for the description of fibre motion. The motion of fibre is considered to consist of translation and rotation around its centre. Equations are provided for both translation and rotation parts of the motion. These equations are based on Newton's second law. The fibre is divided into discrete segments. In this study, only two segments were defined. The centre of the fibre was characterized by a position vector, while its orientation was given by a unit vector. Vectors ri and rj are connections of fibre centre with the centres of the two fibre segments (see figure 1).
This study focuses on only one force acting on the fibre which is the drag force from the surrounding fluid. At each segment, the drag force Fi is decomposed into two components: force Si, which is parallel with fibre's axis and force Di perpendicular to this axis. The schematic decomposition of the drag force can be seen in figure 2.

Fig. 2. Forces acting on a fibre during its movement in the fluid
The relative velocity of the fibre to the velocity of the surrounding fluid is required for the evaluation of force components Si and Di. The velocity of the segment i is the sum of fibre translation velocity up and the contribution from rotation × , where is angular velocity of the fibre around its centre. Inclusion of the term × is the major difference compared to the approach of Lo Iacono et al (2005). The effect of this term will be discussed later in this paper. The relative velocity of a fibre is given by (1) where uf,i is fluid velocity at the center of segment i. The parallel relative velocity is obtained as the projection of relative velocity to the direction of fibre axis defined as (2) and perpendicular relative velocity is obtained in a way that sum of parallel and perpendicular components must result in the total relative velocity: . ( The parallel component of force Si is modelled as the friction force acting on a flat plate of equivalent surface area to the cylinder (A): (4) where  denotes the density of the fluid and (5) In equation (5) Rel is the Reynolds number associated with the particle corresponding to the parallel component and  is the viscosity of the fluid.
The perpendicular force component Di is computed according to the classical drag expression for cylinder of infinite length: (6) where CD is the drag coefficient as measured by Munson and Wieselsberger-Lamb (Munson and Young, 2006).
The overall force acting on segment i is then the sum of parallel and perpendicular components: Only the perpendicular force component Di generates a torque: (8) With knowledge of the force acting on the fibre, the equations of motion can be provided. Newton's second law of motion is applied for the computation of the change of fibre's translation velocity: (9) In the above equation F is sum of forces over all the segments and mp is the mass of the fibre. The change of angular velocity is expressed as: where Ip is the moment of inertia of the cylinder rotating around a perpendicular axis traversing its center of mass.
The solution of the system formed by equations (9) and (10) will give the trajectory of the fibre and orientation of the fibre along its trajectory. The two equations are linked together through the orientation of the fibre. The force acting on the fibre in equation (9) is strongly dependant on the orientation of the fibre.

Test case
The performance and abilities of the introduced model of fibre motion were tested in the case of a simple channel flow with parabolic velocity profile. This choice was driven by the intention to study dynamics of fibre motion in more details, to see the difference between various types of fibres and initial conditions. If more complicated cases of fluid flow were chosen, it could happen that these details would be lost or hidden in the complexity of the flow. The mathematical model described in the previous section was implemented in MATLAB.
The height of the channel was 2H = 0.0056 m, while the length of the channel was limited to 0.02 m. The centre of the coordinate system was in the middle of channel height. The fluid flowed in positive x-direction. Fluid d d = velocity depended only on the distance from the wall (ydirection). The centre-line velocity was 1.35 m/s and the velocity profile was parabolic (average velocity was 0.675 m/s). The density of the fluid was considered to be 1 kg/m 3 , while dynamic viscosity was 1.8e-5 Pa·s. Fibres were from a material with density of 2500 kg/m 3 . Other important parameters characterizing the fibre were its length and diameter. Wall distance, length and diameter were varying parameters in this study. Fibres were released into the channel with zero initial and velocity in all the presented cases.

Results
First, fibres with length of 10 µm and diameter of 1 µm were tracked. Trajectories and orientations of fibres with various initial positions are depicted in figure 3. The initial orientation of fibres was traverse to the flow. The length of the fibres in the figures do not reflect reality. Fibres were enlarged for visualisation purpose only. The effect of velocity gradient on fibre orientation is clearly visible in figure 3. A fibre with initial position almost at the centre of the channel experiences only very small velocity difference at its segments and almost no torque is generated. Therefore, the fibre retains its orientation during the motion. Velocity gradient increases with decreasing distance from the wall of the initial position of fibre and this leads to an increased torque.
Fibres located further from the channel axis tend to align with the flow direction more and more rapidly. It takes only small distance for the fibre located near the wall to orient into the flow direction.
All these observations are summarized in figure 4 where streamwise components of orientation are plotted for various fibre positions. Initially, this component is zero because of initial fibre orientation in wall-normal direction. Than it starts to go toward -1 which indicates that the fibre is aligning to flow direction. At the end of the channel, the orientation vector of the fibre is almost (-1 0) and the fibre is oriented parallel to the wall. In figure 5, the streamwise velocity component is depicted. Fibre accelerates almost immediately to the flow velocity. A fibre with these dimensions has very low particle response time and reacts very quickly to the changes in the velocity of the carrier phase. All the results presented so far were obtained while considering term × in expression (1) for the evaluation of fibre relative velocity. Let's call this term as a damping term. The trajectory and orientation of a fibre with length of 10 µm and diameter of 1 µm and initial position y/H = 0.1 can be seen in figure 6. The damping term was not included in this case. Initially, the fibre starts to orient to the flow direction but then keeps rotating for all the rest of the motion. The fibre won't stabilize and align with flow direction. Neglecting of damping term results in a model in which there is no damping for angular velocity and fibre can rotate without any obstruction. Figure 7 better illustrates the effect of omission of the damping term. Both components of orientation vector are depicted in this figure. Fibre orientation changes periodically and fibre orientation is not stabilized, which is not realistic. The inclusion of the damping term is essential for proper prediction of fibre motion. Now let's increase fibre length to 20 µm and diameter to 5 µm.   It can be seen, that this fibre moves not only in streamwise direction, but there is also a slight shift in wall-normal position. This is caused by non-zero lift force when fibre is inclined to the flow. Actually, this behaviour characterizes all the fibres but it was not so clearly noticeable in the previous cases.
The streamwise components of orientation of fibres with various diameters are compared in figure 9. All the fibres behave very similarly. Change in fibre diameter has only limited effect on the ability of fibres to orient into flow direction. This can be explained by equation (10). When the diameter is increased, then higher force D is acting on fibre segments and more torque is generated. On the other hand, the increase of the diameter will result in the increase of the moment of inertia, as well. The higher torque is partially compensated by higher moment of inertia and this results in similar behaviour of fibres. Very similar statement can be done for fibres with various lengths. More torque is generated because force D is acting at longer lever. But increasing length leads to increase of the moment of inertia as well. This can be illustrated by comparing the curves in figure 9 with the curve representing a fibre with initial position at y/H = 0.5 depicted in figure 4. In the near future the model is intended to be further developed for the case of more complex flows in complicated 3D geometry, such as the airflows in human airways.
This work was supported by the Czech Science Foundation under the grant GA18-25618S.