Developed Numerical Simulation of Falling and Moving Objects in Viscous Fluids under the Action of a Reynolds Lubrication Theory and Low Reynolds Numbers

The development work focuses on the numerical simulations of free body movement in viscous fluid. The aim is to make the simulation of very slow motion of the small body in viscous fluid. We developed bodies’ immersed dynamics simulations in viscous fluid by seeking numerical solutions for appropriate field variables. We developed the methods for vertically and spherically cylindrical objects’ motions, the forces on bodies close to a plane stationary wall are computed from the velocity and pressure fields using the Stokes equation through COMSOL Multiphysics finite element software. The Navier-Stokes equation is reduced to Stokes equation there is independence of time which means object will have an effect only on the motion and the slightly compressible flow assumption is made in order to obtain smooth solution numerically. The forces on an object in slightly compressible Stokes flow have been exerted on the falling objects. The resulting forces have compared with analytical results from the Reynolds Lubrication Theory, and achieved significant results from the development method in Matlab and achieved significant numerical simulations in COMSOL. In addition, an investigation has been made to an object swimming at low Reynolds number. At low Reynolds number moving is possible when object scale is small and flow pattern is slow and sticky. We have developed a system for a thin two-dimensional (2D) worm-like object wiggle that is passing a wave along its centreline and its motion has simulated by the Ordinary Differential Equations (ODE) system and by the Arbitrary Lagrangian-Eulerian (ALE) moving mesh technology. The development method result shows that it is possible for the small object to have a motion from one position to another through small amplitudes and wavelengths in viscous fluid. How to cite this paper: Paul, S. and Oppelstrup, J. (2020) Developed Numerical Simulation of Falling and Moving Objects in Viscous Fluids under the Action of a Reynolds Lubrication Theory and Low Reynolds Numbers. Open Journal of Fluid Dynamics, 10, 8-30. https://doi.org/10.4236/ojfd.2020.101002 Received: December 18, 2019 Accepted: March 2, 2020 Published: March 5, 2020 Copyright © 2020 by author(s) and Scientific Research Publishing Inc. This work is licensed under the Creative Commons Attribution International License (CC BY 4.0). http://creativecommons.org/licenses/by/4.0/ Open Access DOI: 10.4236/ojfd.2020.101002 Mar. 5, 2020 8 Open Journal of Fluid Dynamics S. Paul, J. Oppelstrup


Introduction
In many areas of science and engineering, the flow of viscous fluids at low Reynolds number plays an important role, especially in applied areas of lubrication theory and micro-organism locomotion. This work focuses on the simulation of very slow motion of solid bodies in viscous fluids, inspired by applications in micro-fluidics and bio-mechanics. Viscous flow past in a sphere and in a circular cylinder at very low Reynolds number was first analyzed by Stokes [1].
At low Reynolds number, the inertia of the surrounding fluid becomes unimportant in comparison with viscous forces and thus is generally neglected. Thus the flow is dominated by the viscous forces and the problem is termed as the Stokes-flow problem. It is usually difficult to formulate the exact or approximate solutions for the sphere, ellipsoids and long cylinder body shapes [2]. To solve the Stokes flow problem approximately, two analytical methods are very well known. One is the boundary-value method and the other is the singularity method. The most important special case is that of the exact solution to Stokes' flow problem using the motion of an ellipsoid derived by [3] [4] [5]. Their solution is based on the boundary value problem.
The axi-symmetric slender-body theory is using singularity method at low Reynolds number flow, which has been studied by [2] [6]- [14] and others. More complete formulations of the asymptotic motion of slender axi-symmetric rigid bodies have been done by Geer [15], Keller and Rubinow [16] and non-axially symmetric straight bodies also have been studied by Batchelor [11]. Recently, in the limit of zero Reynolds number, the slender-body theory has been studied by Butler and Shaqfeh [17], Tornberg and Gustavsson [18]. Butler and Shaqfeh [17] have delivered far-field and short-range interaction between particles. For slender bodies of non-spherical shapes at short ranges [17] have used the lubrication theory while for simulation algorithms they have used slender-body theory to model linear and angular velocity and force distribution along the centerline of each thin objects, they have used periodic boundary conditions for non-spherical body shape. The paper by Geer [15] has the treatment of circular cross-section of axi-symmetric rigid body shape and Batchelor [11] has demonstrated the case of non-circular cross-section of non-axial symmetric straight body.
The developed numerical simulation focused on the flow of a viscous fluid at low velocities around a slender body. Considering two different model settings: vertical and spherical shape cylindrical models have been used and theoretical calculation has carried out by Reynolds Lubrication Theory. For the simulation of vertical cylindrical and spherical models, the fluid tank is assumed to be non DOI: 10 [20] surmised that the flagellum would be able to swim provided the average force exerted by the flagella on the fluids is non-zero. Jeffery [5] has been first introduced ellipsoid shape of particle and has given formula for the increase of viscosity.
In addition, the particular aim of the study here is to perform numerical simulations of motion of development microscopic organ in viscous fluids with reference to the stationary boundary wall. For a solitary particle in motion in a viscous fluid, firstly the forces between the boundary walls and the moving body are considered. The worm like 2D body swims past a travelling wave such that its motion is constrained along the centerline of its body. COMSOL3.5a and Matlab software with predefined models for Navier-Stokes and Stokes flow are employed to solve the problem numerically. A moving coordinate system has been used to describe the ellipsoidal particle swimming in the viscous fluid. The resulting force and torque are calculating along with linear and angular velocity of the ellipsoidal body. This is calculated on the basis that net force and torque on the object are zero. The result shows the force and torque that is exerted on the swimming object by the viscous fluid and swimming object is moving slowly in one position to another through small amplitude and wavelengths.

Conceptual Models for Simulation of Falling and Moving Objects in a Viscous Fluid
The cylindrical and spherical cases are designed to simulate the motion of very thin objects under gravity in a fluid. The collective motion of many of such objects is studied in Linström and Uesaka [21], Tornberg and Gustavsson [18].
How such a thin object is related when it approaches the bottom is important, and it is found that the retarding force is the pressure field resulting form squeezing out the liquid between the thin object and the plane bottom [18]. As illustrated by the analysis and calculations here, the thin object, however, does not sense the bottom until the distance is on the order of the fiber diameter, a length scale  (Figure 1(a) and Figure 1(b)).
The fluid velocity u matches the movement of the body surface, and, in response, the fluid exerts a force f on the surface. The computational domain is bounded by a wall, which is assumed to be stationary. In order to obtain smooth and exact solutions, we introduce slight compressibility. This change gives a unique pressure even in the absence of pressure boundary conditions. Since the flow is bounded by (slip or no-slip) walls, one has to ensure that the velocity specified on the body surface is nearly compatible with incompressible flow i.e., the body volume must be held constant. Continuing from here, we compute the velocity and pressure. Then we observe the forces in the models when they have reached near the bottom wall of the tank.
In second experiment Figure 1(c), consider two-dimensional ellipsoid body that is swimming or moving in a viscous incompressible fluid at low Reynolds number. The locomotion of microscopic organisms [22] [23] [24] is an example of swimming at low Reynolds number. The motion of particles in Stokes flows is important in many biological systems. Blood platelets are spheroidal or ellipsoidal and bacteria swim with the aid of flagella or cilia [25].
Swimming at the small scales relevant to cell swimming that tens of micro-meters and below when it's submerged in a viscous fluid [26]. It is possible for small objects to move when fluid is viscous or at low Reynolds number [27]. If the solid object motion is symmetric in time, then it is called as periodic motion [28].

S. Paul, J. Oppelstrup
In situations like this, the reciprocal theorem becomes applicable [29]. Single degree of freedom movement is known to be reciprocal motion [19]. If an object is capable of reciprocal motion then it has one degree of freedom configuration. This implies that time is not of significance and that there is no difference between fast and slow movements. The application of periodic boundary conditions results in periodic motion. Consider an object which has been introduced to a tank which is filled with viscous fluid. In the fluid dynamic equations, if the body moves at some average velocity then the average force must be zero [25] [16] [30]. That means calculation of force and torque on the object by the surrounding fluid in terms of the specified object motion. The unknown linear and angular velocity are determined by the condition that the net force and torque on the organism are zero [31].  [32]. The reduction of Navier Stokes equations to Stokes equations brings in time independence and there is really no difference between slow and quick swimming as demonstrated by [33]. It might also be noted that Stokes flow is also called the reciprocal theorem [26]. The governing equations of Stokes flow are Partial Differential Equations (PDE's) of the elliptic type.

Low Reynolds Number Flow
When we apply the velocity boundary condition on the surface of an ellipsoidal, spherical or any other experimental body, the body being in a viscous fluid tends to slow down. This is due to the effect of viscous forces, which impair the motion of the body. The system behavior is strongly viscous, and the so-called creeping motion or Stoke Equations [34] govern their behavior.

Derivation of Stokes Equation
The Navier-Stokes equations describe flow in viscous fluids with momentum balances for each component of the momentum vector in all spatial dimensions.
The study focus on the boundary movement and swimming at low Reynolds number flow that is either very viscous, very slow, or very small. Low Reynolds number flow for water is about 10 to 4 m 2 ⁄s. Flows with small Re are also called creeping flows or Stokes slow and appear on small scales in porous media, in micro-fluidic devices, and flows with small bodies such as bacteria, amoeba, and innate particles such as blood corpuscles, droplets, bubbles, and mineral grains.
However, the inertial forces are negligible also in very large-scale geological Open Journal of Fluid Dynamics S. Paul, J. Oppelstrup movements such as glacier flow and tectonic plate movement. Important thing is that the low Reynolds number flow plays an important role in the study or lubrication theory, moving micro-organism and many area of biophysical and geophysical interest [27].
The Stokes Flow equations for creeping flow are obtained as the 0 Re = limit of the Navier-Stokes equations for momentum balance augmented by the mass balance law [19], and the fluid flow has made slightly compressible, The slightly compressible Stokes equation can be written as, The Stokes Equation (1) with boundary conditions of specified body surface velocities. This is a linear system of PDE, and it follows that the velocity field and hence the forces exerted by the fluid on a moving body are linear and the body surface velocities are proportional to µ .

Reynolds Lubrication Theory (RLT)
For the analysis of the thin layer fluid flow between cylindrical and spherical object and bottom wall of the tank, the Reynolds Lubrication theory is used which is generally called the Lubrication theory. It is a most important approximation of Stokes equation whenever one of the dimensions of the system is small compared to the others [35] [36].
A body with the surface ( ) , , z x y t = is considered close to a solid plane at 0 z = . One can regard the gap as "thin" provided that "h" is small, and the pressure z-gradient can be neglected as well as the momentum transport terms, Integration across the layer of the mass balance results in the Reynolds lubrication equation, If the body is rigid and moving in the z-direction, we have For the vertical cylinder, we have axial symmetric axis, 2) For Spherical Model (RLT) For the sphere there is axi-symmetry, so with r the radial polar coordinate, S. Paul, J. Oppelstrup

Body Movement
A long body performs a wiggling motion, a sinusoidal lateral (η) displacement of every point on its surface, as traveling waves, This motion, combined with translation and rotation of the local coordinate system, keeps body volume constant, but the body perimeter is not constant. The dimension of the object chosen in this simulation is 0.002 m and the shape of the moving object is elliptical. The wiggling amplitude 0 υ is small, on the order of the worm thickness, and the worm length is a few wavelength 2 λ κ = π . In When we differentiate Equation (5) with respect to t then we get the velocity in the moving frame For the forward swimming motion, we calculate the force and torque that is Then, on the boundary of the wiggling swimmer (see Figure 2), where, ξ must be computed from global coordinates ( ) The sum of internal forces of all bodies are zero, the balance of linear and angular momentum for the deforming body take the form, where "m" and "J" are "small". Since inertia forces are neglected in the momentum balance for the fluid, this is an inconsistent approximation.
However, it shows that if the system is stable, i.e. that the forces tend to act against the velocities, then there should be a well-defined limit as "m" and "J" vanish and it should be the solution of the differential-algebraic system, Equa-S. Paul, J. Oppelstrup tion (12). The numerical implementation uses the Ordinary Differential Equations (ODE) system with small "m" and "J".
The ODE is used for differential algebraic equations or to add degrees of freedom to a PDE (Partial Differential Equations) using the introduced states.
First time derivative represents velocity and second time derivative represents force and moment in ODE setting. In ODE system, there is no need to find consistent initial values for 0 U , 0 V and W. "Smallness" was determined by trial and error in choosing smaller and smaller "m" and "J".
The final set of ODE then becomes, Force F and torque M on the object is calculated [22] for a body in slow motion with a known velocity and rotational rate.
At low Reynolds number flow, swimming motion of a small object depends on different wave parameters.

Vertical Cylindrical Model Using RLT (Figure 3)
Consider the axial symmetric axis and the dynamic simulation of thin rigid fiber body, the force distribution is at zero Reynolds number with slow flow. The vertical cylindrical pressure force on the circular area is calculated by using Reynolds In the minimum separation between particle and bottom wall tank, lubrication dominates the short-range interactions. This short-range hydrodynamic force is calculated using lubrication approximation. The lubrications formulas are derived by [37]. Suppose the fiber has a finite length, a short-range interaction between the fiber and fluid, the lubrication equations of [37] are implemented using a dimensional radius of curvature is equal to the particle length. For the vertical cylindrical case our Reynolds Lubrication Equation is To solve the force distribution the pressure is derived in the fiber coordinates.
( ) The far-field hydrodynamic interaction and the end interaction of thin fiber have no effect. However, the details must be included when calculating the lubrication interactions. The equations in Claeys and Brady [37] are handled in the same manner as above using the lubrication analysis. For the bottom wall tank interacting with a cylindrical body, the total pressure force on the circular area in z-direction becomes,

Spherical Model Using RLT (Figure 4)
The spherical model has a two-dimensional (2D) axi-symetric axis. The object is The spherical polar coordinates are, In addition, the Reynolds lubrication equation becomes, , and total force,

Moving Coordinates (Figure 5)
Consider the two-dimensional (2D) numerical simulations of Stokes flow for the small swimming object at low Reynolds number. We have chosen two reference frames, one is fixed (Global Coordinates), and another is the moving frames of reference (Local Coordinates). The local coordinates moves with fluid particle.
This motion is realised when we apply the moving boundary condition on fluid particle. For the swimming object with moving mesh, we choose the 2D plane axis. For the moving mesh, we have allowed for re-meshing using the arbitrary

Lagrangian-Eulerian (ALE) technique and have added Stokes equation from
MEMS module in COMSOL Multiphysics software.
When we allow for re-meshing, we get three frames: The spatial frame, the reference frame, the mesh frame. The spatial frame coordinates are x and y by default, the reference frame coordinates are X and Y by default and the mesh frame coordinates are xm and ym by default. In the moving case, we have chosen Figure 5. Picture of two frames where particle is moving with local frame. Open Journal of Fluid Dynamics S. Paul, J. Oppelstrup a large tank that is filled by viscous fluid and elliptical object shape that is moving or swimming in the fluid. In Subdomain Setting we use free displacement which means mesh displacement is controlled only by the boundary conditions on the surrounding boundaries. In boundary setting we use two kinds of setting. For tank boundary we add mesh displacement module and in x and y direction add zero because our global coordinate is fixed and there is no mesh displacement in here, and for object boundary we add mesh velocity u and v, which is the free boundary.

Velocity and Pressure Model (Figure 6)
There are some necessary steps to be taken to get a complete solution of vertical cylindrical model in COMSOL. Select geometry, constant parameters, set-up equation system, boundary conditions (applied at the global and subdomain level), and refine mesh or non-uniform mesh is to avoid error in solution between two bottom wall. One tries to choose different orders of elements for finite element solutions like the Lagrange quadratic or the Lagrange cubic. If the difference in solution between using quadratic or cubic elements is negligible, one can retain the lower order element. It is seen in this case that there is practically no difference between solutions arising from the choice of quadratic (2 nd order) or cubic (3 rd order) Lagrange shape functions. Hence, quadratic Lagrange elements are used. Finally, the solution is sought and then post processing is carried out to visualize the solution fields of interest.
The velocity and pressure plots are near the bottom wall of the tank. The pressure plots show some spikes at the cylinder edges. This is owing to the singularities that arise at sharp corners or edges of the cylinders [38]. An option to avoid singularity is to refine the meshes at these locations. This could be done either by increasing the order (p-convergence) or by increasing the number of elements (h-convergence). However, either of these options didn't yield the result and one is led to think that the alternative would only be to smoothen out the sharp edges.

Force: Compare RLT with PDE (Figure 7)
The numerical (PDE) and the theoretical forces (LRT) are calculated, the x-axis corresponds to the ratio R δ describing the motion of the cylinder towards the bottom wall of the tank and plotted the along y-coordinate are the forces on the cylinder calculated numerically and theoretically. The force grows quickly when R β δ = is below 0.05. For values of β greater than this the value of force asymptotically approaches the value of zero (0).
Logarithmic plots are used to find the differences in growth rates for the force in the numerical and the theoretical case. The logarithmic plot also shows that the force in the low ranges of β is proportional to 3 β − . The force plot seems to show a favorable agreement between the numerical and theoretical calculations for the vertical cylinder. It must however be kept in mind that we have introduced the slightly compressibility as opposed to the incompressibility that is otherwise assumed.

Velocity and Pressure Model
Consider the sphere is close to the bottom wall boundary. The pressure field results shown in Figure 8(b) shows no spikes and this seems to indicate that the horizontal and spherical model are best used to avoid unsavoury spikes in the pressure fields.

Force: Compare RLT with PDE
The forces have been computed for the cases of sphere close to the bottom wall of the tank. We examine the force when the spherical model reaches near the bottom wall of the tank. The plots in Figure 9 seem to produce a good agreement for the force fields between the theoretical and the numerical simulations.

Swimming Object Simulation
We have focussed on visualizing the resultsfor swimming object, continuously re-meshingthe object to reach final time step. The geometry of deformations has shown at different time steps in Figure 10. Figure 10 depicts the object in its wiggling motion starting from an initial shape and returning back to its initial shape in a certain time. The surface plots on the moving object depict the velocity at different time steps. The deformation of this object depends upon its boundary conditions. It is seen in the initial time step that the velocity at the edges of the moving object is the highest.
The linear and angular velocity for moving object is calculated on the basis that net force and torque on the object are zero. That means these forces are calculated from the global coordinate system and the velocities are applied on the moving object boundary. Using these conditions, allows one to visualise the motion of the object up till its final time step. The convergence of solution for each of the time step is rather slow and this results in a longer duration for the entire solution to be ready, sine one needs to re-mesh at every time step.
As opposed to the previous case, one sees in this Figure 11 that the object doesn't come back to its initial shape. In every wiggling motion, one sees that the velocity is higher than the previous motion. It shows that that the deviation of the motion from an average straight line is sinusoidal in nature.

Mesh Displacement
The six different images of mesh displacement is at different time steps ( Figure   12). The first image corresponds to a time step of t = 0. In this the mesh is equally distributed around the corners of the swimming object. However in the third, fourth and fifth images one can see the mesh in its deformed form. The object deformation is identical to the way the mesh has deformed. On the left Open Journal of Fluid Dynamics S. Paul, J. Oppelstrup hand side there is compression in the negative x-direction and the right mesh elements are stretched. The meshes in the first image are structured while the remaining meshes are seen to be unstructured.

Velocity Profile
Force and torque are applied on the moving object from the global coordinate system and we get the linear and angular velocity of local variables. The slow motion of the worm is noticed by the change in velocities or motion of the local co-ordinate system.   The result shows that the x-direction mean displacement profile for global coordinate system and its derivatives displacement with respect to time profile of same coordinates system. The x-direction displacement refers to the moving object displacement per time period. Imagine the object to be an ellipsoid moving by combinations of horizontal wave pattern travelling from one place to another with small amplitude and wave length. By using a combination of such oscillations for the motion, the swimmer can move in the opposite direction to that of the surface waves [34]. Figure 13 is a plot of the fluid velocity field directions near the swimming object. Here we see that the domain velocity direction has changed periodically with the object motion, except that they would be translated in the direction of the object wave motion.

Velocity profile observation for different motion cases
It has already been stated that the body motion is periodic if amplitude of motion is small [24] [39]. It's a good approximation for periodic motion of ellipsoid when amplitude and wave length are both small. We use identical idea to observe the mean velocity profile and mean displacement in a periodic motion of the ellipsoid body.
The periodic motion is observed for different amplitudes and wave length to see mean velocity profile for different waveforms. Results from the comparison of these with the calculated velocity profile are shown in below. Velocity profile amplitude 0 0.01 υ = and wave length 0.2

λ =
are default values in this observation. There seems to be a good agreement with the calculated values. The results are summarised in Table 1 below. Observation: From this velocity profile observation it is seen that the default mean displacement and mean velocity values agree very well. One, two and four time periods mean velocity values are almost the same implying that the wiggle motion is small and its moving very slowly in every period.
In the first observation, amplitude is reduced by 10% to 0 0.009   . This experiment results shows only for one and two time periods respectively. But the fourth time period motion could not be obtained with consistent initial values.
This simulation result is less efficient from the first observation mean velocity result and the mean velocity value decreases at small amplitudes. But at amplitude below 0.003 then there is an error suggesting a failure to find consistent initial values. In the third observation, we take the default amplitude value, 0 0.01 υ = and increase 75% the wave length to 1.5 λ = . One gets inefficient wiggling at the change of wave length alone and when the mean velocity value is very small. When we increase wave length over 1.5 then the solvers returns an error.
In fourth observation, we choose a smaller value for amplitude, 0 0.009 υ = and a larger value for wavelength 1 λ = . The velocity profiles remain insufficient and the mean displacement and mean velocity of wiggling object is not smooth. In this case we have changed both wave length and amplitude, and it is seen that in every period of this observation mean velocity value is very smaller than previous observations.

Summary and Conclusions
In this study, we have developed low Reynolds number flow to simulate the ef- for the efficient swimming of small objects, it is seen that wave length should be smaller than the size of the moving object (i.e., 1 λ < ). The development methods in COMSOL 3.5 Multiphysics with Matlab has worked sufficiently well in getting a good enough solution as required in our study.