The Dynamics of Periodically Forced Prolate Spheroids in a Quiescent Newtonian Fluid with Weak Inertia

The effects of convective and unsteady inertia on the dynamics of periodically forced neutrally buoyant prolate spheroids in a quiescent Newtonian fluid medium, at low Reynolds numbers have been modeled. The resulting nonlinear equations have been solved using appropriate numerical methods. Several tests including a perturbation analysis are performed to validate results. A preferred direction, which is identified as the initial direction of motion, is observed that manifests itself in the properties of the solution. Results of the behaviour of various parameters with respect to the Reynolds number, aspect ratio of the spheroid and the amplitude of the periodic force are presented. The results are technologically important as they may lead to insights in the development of active dampeners and smart fluids.


Introduction
Suspensions of solid particles are encountered both as raw materials and as intermediates in a large number of industries such as printing and paper making, petroleum, pharmaceuticals, and food processing. Suspension rheology leads to insights which may lead to better control of fluid stress deformation behaviour and hence may lead to appropriate changes in processing parameters. In most situations, the particles tend to be non-spherical or even irregularly shaped, the suspension rheology then being sensitive to the orientation distribution of the suspended particles. The motion of non-spherical particles in shear flows at vanishingly small Reynolds numbers has been studied theoretically for a long time and is summarized by Leal 1 . It has, in fact been known since the work of Jeffery 2 and later Bretherton 3 at zero Reynolds numbers, that in the absence of inertia, an axisymmetric particle in a simple shear flow rotates periodically in one of an infinite singleparameter family of closed 'Jeffery' orbits. The particular orbit adopted by the particle, in the absence of hydrodynamic interactions, Brownian motion, etc. depends on the initial conditions, rendering the inertialess limit indeterminate. Subramanian and Koch 4 considered both particle and fluid inertia as a possible mechanism acting to remove this indeterminancy. They developed solutions for aspect ratios close to unity. Hence, their analysis captures the leading order effect of the deviation from sphericity on the particle orientational motion. They found that for the neutrally buoyant case, the inertia of the suspending fluid causes a prolate spheroid to drift toward an axial spin about the vorticity axis of the ambient simple shear.
Ramamohan and coworkers 5,6 have been studying the dynamics and rheology of periodically forced suspensions for a period of about two decades. The class of problems they have studied has fundamental importance and technological potential. This class of problems is one of the simplest physically realizable fluid dynamical systems that can show chaos at the level of the individual particle. It has been shown that there exists a chaotic parametric regime, in the dynamics of periodically forced spheroidal particles in a simple shear flow 7 . This chaotic dynamics can be controlled by controlling the system parameters 8 . These results restricted to zero Reynolds numbers and simple shear flow have been summarized by Asokan et al 5 .
Recently Ramamohan et al 6 , have studied the dynamics and rheology of a dilute suspension of neutrally buoyant periodically forced spherical particles in a quiescent Newtonian fluid at low Reynolds numbers. Since most realistic suspensions have non-spherical shapes, The results reported in this work show that data based on simulations with spheres may need to be modified before they can be applied to realistic situations. Hence, in this paper, these results are extended to prolate spheroidal particle suspensions.

The Hydrodynamic Force Expression for an Arbitrary Shaped Particle
Lovalenti and Brady 9 give the expression for the required hydrodynamic force on an arbitrary shaped particle, in the long time limit at low Reynolds numbers. This is the most rigorous equation which results in an ODE at low but non-zero Re. The reciprocal theorem has been used to obtain the following expression. The details of the derivation can be found in Lovalenti and Brady 9 .
The acceleration reaction term (I 1 =int-xx in Table-2) is computed using the expressions given by Pozrikidis 10 11 . The expression for the Stokes resistance tensor (Φ) in its dimensionless form is given by Here, e is the eccentricity of the spheroid and a is the diagonal matrix mentioned in Chwang and Hu 11 .
In the current work, only one dimensional motion of the particle is considered. The translation is along this major axis of the spheroid and hence is symmetric with respect to the particle. The lift force term is neglected as it contributes only to a force in the perpendicular direction and is zero in this case. This has been verified in these computations.

Solving The Differential Equation
After obtaining suitable expressions and values for different aspect ratios, the dynamics of a spheroidal particle were determined.

Formulation of the Problem
The force equation (1) given by Lovalenti and Brady for an arbitrary shaped particle undergoing an arbitrary time-dependent motion at low Reynolds numbers, in the long time limit is considered. In this case a neutrally buoyant prolate spheroid in an infinite body of quiescent fluid is considered, as well as the effects of an external periodic force acting on the spheroid along the x-axis.

Numerical Procedure
The differential equation was solved using a timestepping finite difference routine. In order to accommodate the nonlinear integral term, the product trapezoidal rule 12 was implemented. The value of ε was chosen to be equal to 0.0001. Further decrease in its value did not result in any significant changes of the velocity and displacement values. Note that there exists a singularity (point at which the intergral term becomes indeterminate) at s=t. Hence, in order to avoid this singularity, the integral in the interval [0, t − �] is evaluated, where � is chosen to be a very small number. Note that in the limit s→t, the integral converges to a finite limit and hence the value of the integral in the range s=t -� to s=t is negligible. Under these conditions, equation (1) reduces to Here, ��� � � � �� �� and I 1 is the computed acceleration term. Now, the equation of motion for a neutrally buoyant particle immersed in a liquid is given by The periodic force � ��� ��� � � � ������ is used, where time has been scaled with respect to the frequency of the external periodic force. The following equations for the displacement and velocity of the particle using Newtons second law of motion were obtained. (6) where, Here, a is the characteristic particle dimension, ρ is the density of the particle and µ is the fluid viscosity.
Two sets of data points, 150000 and 350000, taken at an interval of 0.0001 in both the dimensionless velocity and dimensionless position were generated. Further increase in resolution did not yield any significant difference in the results.

Tests
Several tests were performed in order to validate the results. They are listed below.

Test 1 Perturbation Analysis
Perturbation solutions were obtained in order to validate our results for small Re. The Taylor series expansion for the nonlinear integral term was used. One important aspect to be noted is that the expression given for an arbitrary shaped particle (prolate spheroid in our case) is correct upto O(ReSl).
The perturbation parameter was chosen to be Reynolds number, Re. The hydrodynamic force expression for an arbitrary shaped particle given by Lovalenti and Brady is valid up to O(Re).
Hence, the perturbed solution upto O(Re) was expressed as follows The displacement is calculated by numerically integrating the interpolated data of the velocity.
MATLAB was used for computing the above expressions. We found that for low values of Re, typically upto Re=0.2, both the perturbed and the numerical solutions agreed well. Figure 1 shows a comparison between the phase plots obtained by both the methods for aspect ratios 2 and 10, Re=0.1 and ReF = 0.1.

Test 2
When the initial direction of the motion was reversed, namely by replacing ReF with -ReF , the phase space plot was reflected about the zero velocity axis. That is, a reflection of the phase space attractor about the zero velocity axis when the direction of the first motion is reversed was obtained, which can be considered as an important result which demonstrates the correctness of the results.
The results showed a preferred direction in the solution. Since the only physical direction in this problem is the initial direction of the external force, a reversal of that direction should result in a reversal of direction in the solution, which was indeed the case. The tests performed above provide considerable confidence in the results.

Results and Discussion
We have four variable parameters in the system; the Reynolds number Re, the Strouhal number Sl, the aspect ratio and the amplitude of the periodic force ReF. It is essential to determine the effect of these parameters on the system. Typical phase space plots (plots of particle velocity vs. position) have been generated for different values Kumar Figure 1. The phase plot obtained compares the numerical solution with the perturbed solution. This has been seen for Re=0.1, ReF=0.1, aspect ratio=2(a) and 10(b). The match is very good in the second case owing to reduced inertial effects at higher aspect ratios and hence reduced nonlinearity.
of the Reynolds number, the aspect ratio and the amplitude of the periodic force. In some of the Figures (Figure 3b and Figure 5c), the displacement and velocity values were scaled with their appropriate maximum values obtained in the Stokes' case, i.e. when Re = 0. These values have been termed as displacement and velocity amplitude, respectively. This illuminates the results better so that the individual cases can be compared. One of the parameters was chosen, namely the Strouhal number, a constant and equal to unity. The plots represent an attractor as they are bounded in phase space. Since Sl always occurs in combination with Re, the limit of small Sl number is automatically obtained by reducing Re.
It was observed that the particle tends to move away from the zero displacement axis with time, i.e. after every cycle, we notice that the mean position is shifted along the displacement axis. This is termed as drift. The average displacement of the particle is determined from the zero-velocity axis and is denoted as Y pmean . There exists a definite relation between Y pmean and Re, as well as Y pmean and Re F .

Effect of Re F
The area bounded by the phase space plot which is bounded and hence represents an attractor in phase space, increases with increasing amplitude of the forcing term, ReF, establishing the obvious relation between the attractors and the amplitude of the periodic force. As ReF increases, the particle oscillates with larger amplitude and thus covers a larger surface area in the phase plot. As can be seen from the Figure 2(a), the increase in area is quite significant when ReF is increased from 0.1 to 0.5.

Effect of Re
The effect of increasing the Reynolds number can be seen from Figure 3(a). The effect is opposite compared to that of increasing the amplitude of the forcing term. Increasing Re results in a smaller attractor plot. This shows the effect of inertia on the motion of the particle. Inertial effects dominate at higher Reynolds number and the mean  position of the particle is seen to shift in the direction of initial motion on increasing Re.

Effect of Aspect Ratio
The values obtained for the acceleration reaction term for different aspect ratios are given in Table 2. The values of the second and third diagonal elements of the tensor are quite similar to one another. This is expected as they both are symmetric to the direction of motion of particle, in the current work. These values provide an idea about the reaction to particle motion and hence the term acceleration reaction. As can be seen, the values decrease with increasing aspect ratio. These values appear in the term Re' in the equation (6). This term additionally contains a factor which is the square of the inverse of aspect ratio. Hence, the effect of this term decreases with increasing aspect ratio. The contribution from the Pseudo-steady Stokes drag also decreases with increasing aspect ratio which again contributes to lower resistance. The lowering of resistance with aspect ratio could be due to the body becoming more streamlined and hence being able to move through the fluid medium easily. Thus, the overall effect of increasing aspect ratio can be seen as an increase in the area bounded by the attractor. The effect of aspect ratio on Y pmean was studied and it was found that Y pmean increases with aspect ratio. This result can be seen in Figure 4.

Effect of Inertia on Steady State
In order to determine the influence of the inertial effects with every cycle, the simulation was run for longer times, typically about 350000 iterations. The phase plots obtained are presented in Figure 5(a,b).
It can be seen that in the initial cycles, the drift is significant when compared to the drift in the later stages of the simulation. Thus, we see that the inertial effects are dominant during the initial stages and later on the particle tends towards an oscillatory steady state. We observe that the inertial effects dominate at lower aspect ratios and their dominance reduces with increasing aspect ratio. The attainment of an oscillatory steady state is found to be quicker when the aspect ratio is larger. This happens due to weaker inertial effects. Figure 5(c) shows this effect clearly. Here, we scale Y pmean obtained after every cycle of the periodic forcing term with respect to the maximum displacement obtained in Stokes' case. A steady value of this Y pmean amplitude is attained more quickly when aspect ratio increases.

Effect of Nonlinearity
The behaviour of the system was analysed by neglecting the unsteady Oseen correction term which contains the non-linear term and compared it with the attractors obtained from the full differential expression ( Figure 6).
From this, the effects of the non-linear and inertial terms are clear. The phase plots excluding the Oseen correction cover a larger surface area. Another important observation is the absence of drift of the particle from the zero velocity axis. The particle is found to oscillate about a single point. The attainment of steady state is thus strongly dependent on inertia. A slight variation in the parameters considered could vary the inertia and thus approach to steady state.

Conclusion
In this paper, an attempt has been made to determine the dynamics of a prolate spheroid under periodic forcing in a quiescent Newtonian fluid medium at low Reynolds numbers. Inertial effects have been included to study the behaviour more realistically. The numerical values for the acceleration reaction term for different aspect ratios were presented. It is observed that these values decrease with increasing aspect ratio.

1.
The particle is seen to oscillate under periodic forcing. A preferred direction of motion is observed and it is seen that the particle shows a net displacement along this direction with time.

2.
The effect of system variables is studied in detail and it is found that increasing Re restricts the particle motion and hence the size of the attractor.

3.
Increasing the periodic force amplitude is found to increase the size of the attractor.

4.
The effect of the shape of the particle is studied by varying the aspect ratio. The size of the attractor increases with increasing aspect ratio due to weaker inertial effects.
The results were supplimented with detailed physical arguments and wherever possible, various tests have been conducted to justify the results. The ultimate goal is the rheology, which can be further obtained using the results in this paper, since this will determine the stress deformation behaviour and will determine processing parameters for these suspensions. It is hoped that this work excites further research in this area. Future work could possibly cover 2-D and 3-D aspects of such a motion. It would also be interesting to study the effects of coupling rotation with translation. Figure 6. The phase portrait comparing the plots obtained with and without the inclusion of the Oseen correction term and thereby the non-linear term. The effect of inertia on the system can be clearly seen. The non-linear term contribution at higher aspect ratio is reduced and we can see it approaching the curve obtained by neglecting the non-linear term. Here, Re=0.3, ReF=0.1, aspect ratio=2 (a) and 10 (b).