Magnetically Actuated Artificial Cilia: The Effect of Fluid Inertia

Article pubs.acs.org/Langmuir

Magnetically Actuated Artificial Cilia: The Effect of Fluid Inertia S. N. Khaderi,† J. M. J. den Toonder,‡ and P. R. Onck*,† †

Zernike Institute for Advanced Materials, University of Groningen, Groningen, The Netherlands Eindhoven University of Technology, Eindhoven, The Netherlands

‡

ABSTRACT: Natural cilia are hairlike microtubule-based structures that are able to move fluid on the micrometer scale using asymmetric motion. In this article, we follow a biomimetic approach to design artificial cilia lining the inner surfaces of microfluidic channels with the goal of propelling fluid. The artificial cilia consist of polymer films filled with superparamagnetic nanoparticles, which can mimic the motion of natural cilia when subjected to a rotating magnetic field. To obtain the magnetic field and associated magnetization local to the cilia, we solve the Maxwell equations, from which the magnetic body moments and forces can be deduced. To obtain the ciliary motion, we solve the dynamic equations of motion, which are then fully coupled to the Navier−Stokes equations that describe the fluid flow around the cilia, thus taking full account of fluid inertial forces. The dimensionless parameters that govern the deformation behavior of the cilia and the associated fluid flow are arrived at using the principle of virtual work. The physical response of the cilia and the fluid flow for different combinations of elastic, fluid viscous, and inertia forces are identified.

1. INTRODUCTION In the biomedical field, there is an increasing need for miniaturized lab-on-a-chip analysis systems that are able to analyze small quantities of biological samples such as biofluids (e.g., blood, saliva, and urine).1−3 These so-called biosensors are microfabricated total-analysis systems that typically consist of a system of microscopic channels connecting microchambers where dedicated tests are carried out. Classical means of fluid propulsion no longer suffice on these small length scales, which has led to a search for new methods dedicated to fluid propulsion on the micrometer scale, such as micropumps1 and syringe pumps4,5 or by exploiting electromagnetic actuation, as in electro-osmotic6,7 and magnetohydrodynamic devices.8,9 The use of electric fields, however, in transporting biological fluids (which usually have a high conductivity) may induce heating, bubble formation, and pH gradients from electrochemical reactions.10−12 In this work, we explore an alternative way to manipulate fluids in microfluidic systems by mimicking micrometer-scale fluid propulsion mechanisms present in nature. On small length scales, gravity does not play a role and the fluid dynamics is usually dominated by viscosity rather than inertia. An important consequence of this is that a reciprocal motion of an actuator will not lead to net fluid transport. During reciprocal motion, the flow created during the forward motion of the actuator will be canceled by the flow during the backward motion (which exactly traces the path of the forward stroke) even if the rates of forward and backward motion differ. As a result, actuation systems in nature feature nonreciprocal motions that are used by microorganisms for locomotion and fluid propulsion.13 These are often performed using hairlike motile appendages known as cilia (Figure 1a) that beat in an © 2012 American Chemical Society

asymmetric manner with a distinct effective and recovery stroke (Figure 1b). During the effective stroke, the cilia are straight and push a large amount of fluid, whereas during the recovery stroke, they stay closer to the cell surface and pull back a small amount of fluid. The net fluid propelled is in the direction of the effective stroke. The individual asymmetric motion of the cilia leads to a nonreciprocal motion, which, as such, creates flow. In this work, we employ magnetic artificial cilia to propel fluids through lab-on-a-chip microchannels. These magnetic cilia can be realized using thin films consisting of a polymer matrix filled with magnetic nanoparticles, which can be actuated using an external magnetic field14,15 (Figure 1c). Depending on the nature of the particles, the film can be superparamagnetic or ferromagnetic with remanent magnetization. The key challenge is to tune the applied magnetic field to mimic the asymmetrical motion of natural cilia. It was shown in ref 16 that tapered superparamagnetic cilia can exhibit asymmetrical motion when a rotating magnetic field is applied. In addition, also initially curled permanently magnetic cilia can result in nonreciprocal movement because of a magnetic buckling-induced recovery stroke. The optimal channel dimensions and cilia spacing were studied in ref 17, where it was shown that the fluid flow increases and the pressure gradient decreases with increasing channel height. However, both pressure and flow increase when the cilia spacing is decreased. Although the size of natural cilia is on the order of a few micrometers, the size of artificial cilia is on the order of hundreds of micrometers.18−20 Also, the beat Received: January 12, 2012 Revised: March 12, 2012 Published: March 14, 2012 7921

dx.doi.org/10.1021/la300169f | Langmuir 2012, 28, 7921−7937

Langmuir

Article

Figure 1. (a) A microorganism (Opalina) uses hairlike structures called cilia for locomotion (Source: Tamm, S. L.; Horridge, G. A. Proc. R. Soc. London, Ser. B 1970, 175, 219−233). (b) The asymmetric motion of individual cilia consists of an effective and a recovery stroke. (c) Top view of magnetically driven artificial cilia.17 (Left) No magnetic field applied. (Right) With an applied magnetic field. (d) Schematic representation of the numerical model used. The cilia are made of soft polymer films with embedded magnetic nanoparticles. A magnetic field is used to actuate the artificial cilia creating fluid transport.

dynamic equations. In section 3, we derive the governing dimensionless parameters that account for the interaction of elastic, fluidic, and magnetic forces. In section 4, the results are presented where the asymmetric motion of the artificial cilia and the created fluid flow are analyzed in the presence of fluid inertia. Finally, conclusions are drawn in section 5.

frequency of the cilia can be controlled easily with the frequency of the magnetic field. These two aspects allow the artificial cilia to operate beyond the constraint of low Reynolds numbers, which could drastically affect the solid−fluid interaction between the magnetoelastic deformation and imposed flow. This will be explored in this article. Recently, artificial counterparts of natural cilia that are actuated by electrostatic forces,19 magnetic forces,14,15,21−25 base excitation,18 and light actuation have been synthesized. From a theoretical point of view, the mechanics of artificial cilia has predominantly been analyzed in the limit of low Reynolds numbers.16,26−29 Inertial effects have been studied only in the case of mixing30 and for the kinematically constrained motion of cilia.31,32 In this article, we will study the effect of fluid inertia in microfluidic propulsion by means of the actuation of magnetically responsive elastic cilia. The proposed artificial cilia are thin films that we model as a collection of Euler−Bernoulli beam elements. This Lagrangian model of the cilia accounts for elastic and inertial forces in a nonlinear geometry setting. The Lagrangian solid dynamics model is coupled to an Eulerian fluid dynamics model that solves the Navier−Stokes equations using the fictitious domain method.33,34 Maxwell’s equations are solved at each configuration, yielding the local magnetic field from which the magnetic forces can be obtained. This magnetomechanical model allows the study of the overall deformation of the cilia actuated by an external magnetic field and the resulting fluid flow. The dimensionless numbers that quantify the relative contribution of different forces (elastic and inertia forces of the film, magnetic forces, and viscous and inertia forces in the fluid) are established using the virtual work equation. The goal of this article is to explore the interaction of the magnetoelastic deformation of the cilia with the surrounding fluid with a special view to understanding the effect of fluid inertial forces on the flow. We will identify the specific parameter combinations for which flow is maximal and unidirectional. The article is organized as follows. In section 2, we delineate the computational model in which we simultaneously solve the solid dynamics equations, the Maxwell equations, and the fluid

2. COMPUTATIONAL MAGNETOELASTIC SOLID−FLUID INTERACTION MODEL 2.1. Solid Dynamics Model. The artificial cilia are modeled as thin Euler−Bernoulli beams (using plate theory with a single axis of bending) with plane strain conditions in the z direction. It is assumed that the solid material is isotropic and linearly elastic, as specified by the elastic modulus E and Poisson’s ratio v. The principle of virtual work in the solid dynamics problem is

∫V (σδε + ρ(ü δu + v ̈ δv)) dV ⎛ ∂δv ⎞ dV + ∫ (tx δu + t y δv) dA = ∫ fx δu + f y δv + Nz ⎝ ∂x ⎠ ⎜

⎟

(1) where u and v are the axial and transverse displacements of a point along the beam and ρ is the mass density. Furthermore, σ is the axial stress and ε is the corresponding strain, f x and f y are the magnetic body forces in the axial and transverse directions, A = bh is the crosssectional area, Nz is the magnetic body couple in the out-of-plane direction, tx and ty are the surface tractions, and b and h are the out-ofplane width and thickness of the film, respectively. We follow the approach used in ref 35 to linearize and discretize the principal of virtual work to arrive at the discretized form of the virtual work equation at time t +Δt (Appendix A) +Δt δpT (KΔp + Mp̈ t +Δt − Ftext + Ftint) = 0

(2)

where K is the stiffness matrix that combines both material and t+Δt is the external geometric contributions, M is the mass matrix, Fext t is the internal force vector, Δp is the nodal force vector, Fint displacement increment vector, and p̈ is the nodal acceleration vector. 7922


Langmuir

Article

coupling is established using the method of Lagrange multipliers.37 In what follows, we apply the constraint to the fluid dynamics model and then add the virtual work done by the constraint force to the solid dynamics model and finally couple the equations so that the solid and the fluid equations of motion can be solved implicitly. We first add the variation of the Lagrange multiplier times the constraint to the fluid dynamics model,

2.2. Fluid Dynamics Model. The principle of virtual work in the rate form for the fluid problem is36

∫V σij δDij dV + ρ f ∫V

d ui ∂u δui dV + δp i dV = 0 dt ∂xi V

∫

(3)

where the first and second terms represent the work due to the internal stresses and inertia forces in the fluid, respectively, and the third term imposes the incompressibility condition. In eq 3, σij represents the components of the stress tensor in the fluid, Dij represent the components of the deformation rate tensor in the fluid, ρf is the fluid density, ui represents the components of fluid velocity in the ith direction, d()/dt = ∂t + ui∂()/∂xi represents the total derivative, p is the pressure, and dV = b dx dy. (u represents the axial displacement of a point on the beam, and ui represents the fluid velocity.) The constitutive relation for the fluid is σij = −pδij + 2μDij, where μ is the fluid viscosity. It is to be noted that x represents a spatial point in the fluid domain (Eulerian point of view) and a position occupied by a material point on the beam in the solid domain (Lagrangian point of view). The following shape functions are used to interpolate the velocity and pressure

δUT KUPP + δUT K̂ f U − δUT F f + δPT (KUP)T U Jf

+ δ(λ iJ (ui − pi̇ J )) = 0 uiJf

where and represent the components of the fluid and solid beam velocity, respectively, in the ith direction at the location where the Jth node of the solid beam is present. Discretizing the constraint equation (Appendix B) yields

δUT KUPP + δUT K̂ f U − δUT F f + δPT (KUP)T U + δUT ϕ J λJ + δλJT (ϕ JT U − Aṗ J ) = 0 where A is a matrix that eliminates the rotational degrees of freedom from ṗ. It can be noted that the term δUTϕ Jλ J in the above equation represents the virtual work done by the Lagrange multiplier (i.e., the Lagrange multiplier acts as a force due to the constraint). Hence, the components of the Lagrange multiplier are to be interpreted as fluid drag forces. Invoking the arbitrary nature of the virtual fields and performing the standard finite element assembly yields the following set of equations

⎡ ϕ 0 ϕ 0 ϕ 0 ... ⎤ ⎧ u1 ⎫ 1 2 3 ⎥{ u1, u1, ... }T ⎨ ⎬ =⎢ ⎢⎣ 0 ϕ 0 ϕ 0 ϕ ... ⎥⎦ 1 2 ⎩ u2 ⎭ 1 2 3 = ϕIiUI = ϕT U

p = [ ψ1 ψ2 ψ3 ... ]{ p1 , p2 , p3 ... }T = ψI PI = ψT P

K̂ f U + KUPP + Φλ = F f

(4) J where, ui represents the components of fluid velocity in the ith direction at the Jth node, pJ represents the magnitude of the pressure at the Jth, node and ϕ J and ψ J are the shape functions used to interpolate the nodal velocities and pressure, respectively. The shape functions are chosen so that the resulting element satisfies the Inf-Sup condition. In this work, we use the Taylor−Hood Q2Q1 element, which interpolates the velocities quadratically and the pressures linearly.36 Using eq 4, we evaluate the various terms in the virtual work equation (Appendix B): T

UU

∫V σij δDij dv = ∫V (–pδij + 2μDij)δDij = δU (K

UP

U+K

(KUP)T U = 0 ΦT U − Aṗ = 0

∂ϕ

t +Δt δpT (KΔp + Mp̈ t +Δt − Fext + Ftint) − λT Aδp = 0

λ = { λ1, λ2 }T = { λ11, λ12, λ12, λ 22 }T

P)

is the Lagrange multiplier vector and λiJ is the Lagrange multiplier at the Jth node in the ith direction. In general, the node of the solid beam is present inside a fluid element and does not coincide with a fluid node so that the constraint is satisfied in an approximate sense. Because virtual quantity δp is arbitrary, we can write

∂ϕ

i

i

= δPT (KUP)T U

(6)

+Δt + Ft − AT λ = 0 KΔp + Mp̈ t +Δt − Ftext int

Assuming that the solution at time t is known, the convective term is linearized with respect to the known solution at time t and the time derivative of velocity is discretized using an implicit Euler scheme. Neglecting the higher-order terms, we get (Appendix B)

∫V

duit +Δt δu i d V = ρ f dt V

∫

̂ ̇ t +Δt − AT λ = Fts+Δt Kp

̂ + δUT K1U + δUT K2U − δUT F f = δUT MU

1 K̂ = (Δt 2βK + M) γΔt

(7)

⎡ ⎤ 1 +Δt − Ftint − Δt K⎢ṗ t + Δt(1 − 2β)p̈ t ⎥ Fts+Δt = Ftext ⎣ ⎦ 2 t t t + Δ − γ p p (1 ) ̇ ̈ + (Δt 2βK + M) γΔt

where K̂ f = K +M̂ +K +K . 2.3. Fluid−Structure Interaction. We now couple the Lagrangian beam element formulation of section 2.1 to the Eulerian finite element fluid model of section 2.2. The solid beam is considered to be an internal boundary of the fluid domain. The solid and fluid models are coupled using the constraint that the velocity (or displacement) of the solid is equal to the velocity (or displacement) of the fluid. This 1

(12)

where

Substituting the above expressions into eq 3 yields its discretized form

UU

(11)

The motion of the film with time is obtained by solving eq 11 with appropriate initial and boundary conditions. Newmark’s integration scheme is used to integrate eq 11 in time, on the basis of which eq 11 can be rewritten to solve for the nodal velocities as

⎛ ∂u t +Δt ∂u t +Δt t +Δt ⎞⎟ ⎜δui i uj + δu i i ⎜ ⎟ dV ∂t ∂xj ⎝ ⎠

δUT KUPP + δUT K̂ f U − δUT F f + δPT (KUP)T U = 0

(10)

where

∫V δp ∂xi dV = ∫V δpI ψI ∂xJi UJ dV = δpI ∫V ψI ∂xJi dVUJ i

(9)

The fluid drag forces are now considered to be nodal forces whose virtual work is an additional contribution to eq 2. This leads to

(5) ∂u

(8)

pi̇ J

2

(13)

and γ and β are integration parameters. Equation 12 contains the discretized equations of motion for one beam element. After performing the standard finite element assembly procedure,38 we get 7923


Langmuir

Article

magnetic field far away from the film, and N is the total number of segments. Equation 21 clearly shows that the magnetic field H in the film is the sum of the external field H0 and the demagnetizing field Hself. By rotating Hj back to the local coordinates, we get

the discretized equations of motion for the whole film (dropping the superscript (t + Δt)):

K̂ sṗ − AT λ = Fs

(14)

Combining the equations of motion for the solid (eq 14) and fluid (eq 9) results in

⎤ ⎡ ̂ KUP 0 Φ ⎥ ⎢ Kf ⎧U ⎫ ⎧F f ⎫ ⎥⎪ ⎪ ⎪ ⎪ ⎢ UP T P⎪ ⎪ 0 ⎪ K ( ) 0 0 0 ⎥⎪ ⎢ ⎨ ⎬=⎨ ⎬ ⎥ ⎢ ṗ ⎪ ⎪ F ⎪ s K̂ s − AT ⎥⎪ 0 ⎢ 0 ⎪ ⎪ ⎪ ⎪ ⎩ λ ⎥ ⎭ ⎩0⎭ ⎢ ⎣ ΦT 0 −A 0 ⎦

Ĥ j = RTj H 0 +

For the situation of a permanently magnetized film with magnetization M̂ i, with i = 1,..., N, eq 22 gives the magnetic field in all of the segments. However, in the case of a superparamagnetic film, the magnetization M̂ j is not known a priori but depends on the local magnetic field through

(15)

M̂ j = χ̂ Ĥ j = χ̂ RTj H 0 +

⎡ χ̂ χ̂ ⎤ ⎢ x xy ⎥ χ=⎢ ⎥ ⎣ χ̂xy χ̂y ⎦ There are N similar pairs of equations. In total, there are 2 × N equations for the 2 × N unknown magnetizations. This set of equations is solved to get the magnetization with respect to the local coordinate frame. From the magnetization, the field can be found from eq 22 and the magnetic flux density can be found by using the constitutive equation. The advantage of the proposed method is that we need not model the medium around the magnetic film to determine the magnetic field in the film. The verification of the model is performed by using a superparamagnetic film and is given in Appendix D. Once the magnetization is calculated by solving eq 23, the magnetic body couple and body force per unit volume can be found from N = M × B0 and f = M·∇B0 (with B0 = μ0H0) and is given as input to eq 1.

(17)

Because ∇ × H = 0, a scalar potential ϕ exists such that H = −∇ϕ. Substituting this into eq 17 yields a Poisson equation for ϕ, ∇2ϕ = −∇·M. By taking into consideration the effect of discontinuity in the medium, the general solution of the Poisson equation can be found,39 resulting in

1 4π

∮ n|x′·M−(xx′|′) dS′ +

1 4π

M(x′) dV ′ ∫ ∇′· |x − x′|

(18)

where n′ is the outward normal to the surface of V. The magnetic field H(x) can be found from the gradient of ϕ(x). We now discretize the film into a chain of rectangular segments, within which the magnetization is assumed to be uniform so that ∇′·M = 0 and the volume integral vanishes. The field is now due to only the jump in magnetization across the surface of each segment and is given by the surface integral in eq 18. The magnetic field in coordinates local to segment i (denoted by ∧) due its four surfaces can now be calculated at any position (x ̂ , y )̂ by evaluating the surface integral in eq 18, resulting in Ĥ i = GiM̂ i

3. DIMENSIONAL ANALYSIS In this section, we use the principle of virtual work to identify the dimensionless parameters that govern the deformation behavior of the artificial cilia. Considering only the transverse deformations and magnetic body torques, we have 2

where M̂ = [M̂ x M̂ y]T, Ĥ = [Ĥ x Ĥ y]T, and Gi can be obtained from eq 18 (Appendix C). Note that in this section Einstein’s summation convention is not applied. The field due to segment i with respect to the global coordinates is Hi = RiĤ i, where Ri is defined as

∫ λδvb dx = 0

where the first term represents the virtual elastic work done by the internal moments, the second term represents the virtual work done by the inertial forces of the beam, the third term represents the virtual work done by the magnetic couple, and the last term represents the work done by the fluid drag forces. In the above equation, λ is the traction due to fluid drag on the film in the transverse direction and has units of force per unit area. This is in contrast to the λ used in section 2.3, which has units of force per unit out-of-plane width. We introduce the dimensionless variables V, T, and X such that v = VL, x = XL, and t = Ttref, where L is a characteristic length (taken to be the

(20)

N

∑ R iGijM̂ i = H0 + Hself i=1

2

−

with θi being the orientation of segment i with respect to the global coordinates. The field at any element j because of the magnetization of all of the segments throughout the film is Hj = H 0 +

2

∫ EI ∂∂xv2 ∂∂xδ2v dx + ∫ ρA ∂∂t 2v δv dx − ∫ Nz ∂δ∂xv A dx

(19)

⎡ cosθi −sinθi ⎤ ⎥ Ri = ⎢ ⎢⎣ sinθi cosθi ⎥⎦

(23)

with

with the constitutive relationship B = μ0(M + H), where B is the magnetic flux density (or magnetic induction), H is the magnetic field, M is the magnetization that includes the remnant magnetization, and μ0 is the permeability of free space. Substituting the constitutive equation into eq 16 yields

ϕ(x) = −

N

∑ χ̂ RTj R iGijM̂ i i=1

(16)

∇·H = −∇·M

(22)

i=1

This set of equations is solved to obtain the velocities at the solid and fluid nodal points, the pressure in the fluid, and the Lagrange multipliers at the solid nodal points. It is to be noted that in eq 15 the velocities of the film and the fluid are simultaneously solved for every time increment. This approach is commonly referred to as the monolithic approach. 2.4. Magnetostatics. Maxwell’s equations for the magnetostatic problem with no free currents are39 ∇·B = 0, ∇ × H = 0

N

∑ RTj R iGijM̂ i

(21)

where Gij properly accounts for the relative positioning of segments i and j and the geometry of segment i, H0 is the externally applied 7924


Langmuir

Article

distance a. Introducing these length and time scales into the above equation leads to

length of the cilia) and tref is a characteristic time. Substitution yields

∫

⎛ Ebh3 ∂ 2V ∂ 2δV ρbhL2 ∂ 2V ⎞ ⎜⎜ + 2 δV 2 ⎟⎟ dX 2 2 2 t ref ∂T ⎠ ⎝ 12L ∂X ∂X −

⎛

⎞ − λL δVb⎟ dX = 0 ⎠

∫ ⎝Nzhb ∂δ∂XV ⎜

(24)

from which the elastic (Ebh3/12L2), the inertial (ρbhL2/tref2), the magnetic (Nzhb), and the viscous (λLδVb) terms can be easily identified. By normalizing with the elastic term, we get

∫

⎛

⎛

⎞

⎞

∫ ⎜⎝Mn⎝ ∂δ∂XV ⎠ + FnδV ⎟⎠ dX = 0 ⎜

⎟

(25)

Here, the three governing dimensionless numbers are defined as the inertia number 2 ⎛ ρ ⎞⎛ L ⎞ ⎛ L ⎞ 2 In = 12⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎝ E ⎠⎝ t ref ⎠ ⎝ h ⎠

(the ratio of inertial to elastic force), the magnetic number

⎛ N ⎞⎛ L ⎞ 2 Mn = 12⎜ z ⎟⎜ ⎟ ⎝ E ⎠⎝ h ⎠ (the ratio of magnetic to elastic force), and the fluid number

⎛ λ ⎞⎛ L ⎞3 Fn = 12⎜ ⎟⎜ ⎟ ⎝ E ⎠⎝ h ⎠ (the ratio of fluid to elastic force). From dimensional considerations, λ should scale with μ/tref, leading to Fn = 12(μ/Etref)(L/h)3. To identify the dimensionless parameters that govern the fluid flow, we start from the Navier−Stokes equations for the fluid du −∇·σ + ρ = 0 in A dt

(26)

− σ ·n + λ = 0 on Γs

(27)

⎛ ∂u ⎞ ⎛ ∂ 2u ⎞ DH⎜ x̅ ⎟ ≈ ⎜⎜ 2x̅ ⎟⎟ for L < y < H ⎝ ∂ t ̅ ⎠ ⎝ ∂y ⎠ ̅

(30)

Nz = M̂ x Bŷ − M̂ y Bx̂ = M̂ x Bŷ = M̂ x B0f ̂ (θ)

(31)

where θ is the film orientation and B0 is the amplitude of the applied magnetic field. To identify the form of the magnetic couple for an SPM film, we have to analyze how the magnetic field inside an SPM film depends on the applied magnetic field, the geometry of the film, and the magnetic properties of the film. To this end, let us subject the SPM film to an external magnetic field (Hx0, Hy0). By assuming that the magnetization induced by an applied field is uniform, we can calculate the magnetic field induced by the magnetization Hself (Appendix C)

The ciliary motion creates a flow that has a dominant velocity in the channel direction (i.e., the x direction). Therefore, we consider only the x component of eq 26. Using the constitutive relation for the fluid gives ⎛ ∂ 2u ⎛ ∂u ∂u ⎞ ∂ 2ux ⎞ ∂p ⎟⎟ ρ f ⎜ x + ux x ⎟ = − + μ⎜⎜ 2x + ⎝ ∂t ∂x ⎠ ∂x ∂y 2 ⎠ ⎝ ∂x

(29)

where the terms in parentheses are nondimensional, DH = ρf(H −L)2/μtref = tdiff/tref is the diffusion number, and Re = ρfL2/μtref is the Reynolds number. The diffusion number signifies how long it takes for the momentum to diffuse into the fluid (tdiff) compared to tref, whereas the Reynolds number Re quantifies how large the inertia forces are compared to the viscous forces. (The fluid velocity near the cilia, the velocity gradient, and the viscous energy dissipated per unit time scale with L/tref, 1/tref, and μ/tref2, respectively. The inertia forces scale with ρL/tref2, hence the kinetic energy input per unit time to the system scales with ρ(L/tref2)(L/tref). Their ratio gives the Reynolds number Re.) The question one might ask is, which of these dimensionless numbers is important? We show with the help of a simple analytical model (Appendix E) that DH has two effects; first, it determines the duration of the transient period, and second, it reduces the fluctuating component of the fluid transported. However, the mean propulsion velocity created by the cilia in the steady state is independent of DH, leaving the Reynolds number Re to be the main parameter that governs the fluid transported. Another interesting observation is that the convection terms in the Navier−Stokes equation do not contribute to the fluid momentum (as also verified by our simulations). Before proceeding, we identify the origin of the magnetic couple Nz for the two magnetic material systems considered in this work: (i) permanently magnetic (PM) and (ii) superparamagnetic (SPM) materials. For permanently magnetic materials, having a local remanent magnetization (M̂ x) pointing along the axial direction of the film, (M̂ x, M̂ y) = (M̂ x, 0), we can write

⎛⎛ ∂ 2V ∂ 2δV ⎞ ⎛ ∂ 2V ⎞⎞ ⎜⎜⎜⎜ ⎟⎟ + In⎜⎜ δV ⎟⎟⎟⎟ dX 2 2 ⎝ ∂T 2 ⎠⎠ ⎝⎝ ∂X ∂X ⎠ −

⎛ ∂u ⎞ ⎛ ∂ 2u ⎞ Re⎜ x̅ ⎟ ≈ ⎜⎜ 2x̅ ⎟⎟ for 0 < y < L ⎝ ∂ t ̅ ⎠ ⎝ ∂y ⎠ ̅

(Hx self , Hy self ) = ( −aMx , −bM y)

(28)

where α = 2 tan−1(h/L)/π and β = 2 tan−1(L/h)/π are positive factors that depend on the geometry of the film. The field Hself is nearly in the direction opposite to the magnetization, hence Hself is called a demagnetizing field. For slender films h ≪ L, the factor α approaches zero and the factor β approaches unity. This makes the demagnetizing field negligible along the length (or tangential direction) and equal to the magnetization in magnitude along the thickness (or transverse direction). The magnetic field H in the film is the sum of the external field H0 and the demagnetizing field Hself; see eq 21. Because the

The relevant length scales can be identified by looking at the mechanism of fluid flow inside the channel. The cilia of length L are placed periodically at a distance a in a channel of height H. The fluid between the cilia is directly driven by them, and the momentum diffuses from this region upward into the channel. The relevant length scales are the cilia spacing a in the x direction and L and H−L in the y direction, and the relevant time scale is tref. Because the velocity of the fluid after a units will be the same (in the case of uniformly beating cilia), the net velocity and pressure gradients in the x direction vanish over a 7925


Langmuir

Article

h = 2 μm at the fixed end to 1 μm at the free end, ρf = 0.1 kg/m3, μ = 0.1 m Pa s, ρ = 1000 kg/m3, B0 = 18.8 mT, and the cycle time is tref = 2tbeat = 20 ms. Note that the fluid density and viscosity are chosen so that the Reynolds number and fluid number are small. These values can also be obtained for a different set of fluid properties by properly tuning the frequency. The spacing between the cilia a is equal to 2L. (For a discussion onf the effect of cilia spacing, the reader is referred to Khaderi et al.17) The magnetic susceptibilities of the cilia are 4.6 along the length and 0.8 along the thickness, based on a 20% volume fraction of 8.4 nm Fe3O4 superparamagnetic particles.40 For the simulations, each cilium is divided into 40 beam elements and the fluid domain is divided into 28 × 30 elements. The fluid propelled is characterized by two parameters: the net volume of the fluid transported during a ciliary beat cycle and the effectiveness. The velocity field in the fluid at any position x integrated along the channel height gives the instantaneous flux through the channel. This flux, when integrated in the direction of the effective and recovery strokes, gives the positive (Qp) and negative (Qn) flows, respectively (Figure.2). Because of the asymmetric motion, the positive flow is larger than the negative flow, generating a net area flow per cycle (Qp − Qn) in the direction of the effective stroke. The effectiveness, defined as (Qp − Qn)/(Qp + Qn), indicates which part of the totally displaced fluid is effectively converted into a net flow. An effectiveness of unity represents unidirectional flow. Figure 3a shows the deformed geometry of a typical cilium at different instants in time when a rotating magnetic field is applied. The blue arrow heads and the red arrows represent the direction of the magnetization and the applied magnetic field, respectively. The dashed curve represents the trajectory of the tip of the cilium. The components of the magnetic field H in the coordinate system local to the cilia are plotted in Figure 3b,c. The curves with square markers refer to the demagnetizing field or the self-field (Hself, see eq 21), the triangles refer to the applied magnetic field (H0), and the lines without markers refer to the total field. The magnetization and magnetic couple are plotted in Figure 3d,e, respectively. The contours of absolute fluid velocity normalized with L/tbeat are shown in Figure 4a−f. The white arrows represent the applied magnetic field, and the white circles represent the fluid particles. The instantaneous flux and the evolution of total flow during the beat cycle are shown in Figure 4g. When a rotating magnetic field is applied, the free end of the beam nicely follows the applied magnetic field (Figure 4a,b) and is rotated through 135° to form a U bend near the fixed end (instant 3 in Figure 3a). From the magnetic field distribution along the cilia at these instants (Figure 3b,c), it can be seen that the demagnetization is stronger (due to the shape anisotropy) in the transverse or thickness direction (Figure 3c) compared to that in the tangential direction (Figure 3b), which reaches a maximum at a position on the cilia where the magnetic field is orthogonal to the cilia. This occurs near the fixed end at time instant 2 and near the half-span of the cilia at instant 3; see the corresponding instants in Figure 3a. The magnetic field in the transverse direction is significantly lower than the applied magnetic field, whereas this difference in the case of the tangential field is only marginal (see the lines without markers and the lines marked with triangles in Figure 3b,c). Because of the formation of the U bend, the magnetic field in the tangential direction changes sign at time instant 3 in Figure 3a. (See the corresponding instant in Figure 3b.) The magnetic field H in Figure 3b,c leads to a magnetization as shown in Figure 3d. It can be seen that the

magnetization M and H are related through χ, we have two equations that can be solved to find the magnetization: ⎛ Hy 0 ⎞ Hx 0 ⎟ (Mx , M y) = ⎜⎜χ̂x , χ̂y 1 + χ̂yβ ⎟⎠ ⎝ 1 + χ̂xα

The magnetization along the length is enhanced by the shape of the film (through α and β) and also by the magnetic anisotropy (for χ̂y < χ̂x). The magnetic flux density reads B = μ0 (M + H0 + Hself). It can be seen that the demagnetizing field Hself has the effect of decreasing the magnetic flux density. For a slender film, the magnetic flux density in the thickness direction is nearly equal to the applied magnetic field (because Hyself ≈ −My) and along the length it is the sum of the magnetization and the external magnetic field (because Hxself ≈ 0). Because we know the magnetization in terms of the applied magnetic field, we can calculate the magnetic couple (Nz) acting on the film: μ0H02 sin 2θ(χ̂x − χ̂y + χ̂xχ̂y(β − α)) Nz = 1 + αβχ̂xχ̂y + αχ̂x + βχ̂y

(32)

The magnetic couple is directly proportional to the square of the applied field, depending on twice the angle made by the magnetic field vector with the film, the difference in the susceptibility (magnetic anisotropy), and the difference between α and β (geometric anisotropy). The main message from eq 32 is that for a magnetic couple to act on the film we need an anisotropy, either magnetic or geometric. Because the magnetic cilia are thin magnetic films, shape anisotropy is always present.

4. RESULTS 4.1. Cilia Deformation and Fluid Flow at Low Reynolds Numbers. In this section, we study the fluid propulsion created by a 2D array of tapered superparamagnetic artificial cilia, which are actuated by a rotating magnetic field, Bx = B0 cos(ωt ), By = B0 sin(ωt )

(33)

that is uniform in space, where B0 is the magnitude of the applied magnetic field, ω = 2π/tref is the angular frequency, and tref is the time period of rotation of the magnetic field. Because the magnetic couple scales with sin 2ωt (eq 32), the time taken by the cilia for one beat cycle is tbeat = tref/2. To perform the simulations, we take a unit cell containing one cilium. The top and bottom boundaries are no-slip boundaries, and the left and right ends are periodic (Figure 2).

Figure 2. Periodic unit cell used for the analysis shown by dashed lines.

The fluid, Reynolds, magnetic, and inertia numbers are set to 0.015, 0.001, 10.89, and 4.8 × 10−3, respectively, indicating that the fluid forces are small compared to the elastic forces but are large compared to the fluid inertia forces. The values of the dimensionless parameters correspond to L = 100 μm and E = 1 MPa, the thickness of the cilia decreases linearly from 7926


Langmuir

Article

Figure 3. (a) Deformed configurations of a typical cilium at different instants in time. The blue arrowheads represent the direction of the magnetization of the cilium, and the thick, red arrows at the free end represent the applied magnetic field. Instants 1, 2, 3, 5, and 6 correspond to time instants of 0.15tbeat, 0.5tbeat, 0.8tbeat, 0.85tbeat, and 0.89tbeat, respectively. (See also Figure 4.) The dashed line indicates the trajectory of the cilia tip, which encloses an area termed the area swept. (b) Magnetic field tangential to the local coordinates of the cilia. For the legend of different line types, see plot d. (c) Magnetic field transverse to the local coordinates of the cilia. Curves without and with triangular and square markers represent the total, applied, and demagnetizing field, respectively. (d) Magnetization of the cilia in local coordinates of the cilia. (e) Magnetic body couple along the length.

the portion of the film between the free end and the U bend, making the film more curved, and clockwise (negative) moments to act on the portion between the fixed end and the U bend, making the film come back to the initial position (instants 3 and 5 in Figure 3e). Because of the tapered nature of the film, the clockwise moments are larger than the anticlockwise moments. Under the influence of such a system of moments, the U bend propagates to the free end and the cilia come back to the initial position. From Figure 3a, it can be seen that the artificial cilia show an asymmetric motion similar to that of natural cilia.41 The dashed

magnetization in the tangential direction is much larger than that of the transverse direction (due to shape and magnetic anisotropy). Therefore, the magnetization is always tangential to the cilia. The sign reversal of the magnetic field at instant 3 leads to a direction reversal of the magnetization. In the region between the U bend and the free (or fixed) end, the magnetization points toward the free (or fixed) end. The magnetization vector (blue arrowheads) plotted on the cilia in Figure 3a clearly shows this behavior. This magnetization interacts with the applied magnetic field, causing anticlockwise (positive) moments to act on 7927


Langmuir

Article

Figure 4. (a−f) Contours of absolute velocity (normalized with L/tbeat) at different time instants. The direction of the velocity is given by the streamlines, and the white arrows represent the magnetic field at the respective time instants. The circles represent fluid particles. The parameters used are Fn = 0.015, Re = 0.001, and Mn = 10.89. Four unit cells are shown for clarity. (g) Instantaneous flux and accumulated flow as a function of time. The time instants corresponding to a−f are duly marked. Note that the cilia have completed one beat cycle while the magnetic field rotates by π.

line represents the trajectory of the cilia tip, which encloses an area termed area swept. The amount of fluid propelled in the Stokes regime scales linearly with the swept area.16 It is to be noted that the effective stroke takes place for a longer duration of time than does the recovery stroke, whereas in natural cilia the opposite is the case,41 but this holds no consequences for the flow generation due to the operation in the Stokes regime. Let us now analyze the fluid flow due to the ciliary motion. When the cilia move, they set the nearby fluid into motion. The momentum diffuses to the rest of the fluid because of the viscous forces. Because the viscous forces are large compared to the inertia forces, this momentum diffuses instantaneously and the fluid velocity nicely follows the cilia velocity at all time instants. The cilia create a positive flux (to the left) during the effective stroke (see the instantaneous flux at the instants corresponding to Figure 4a,b in Figure 4g) and drag back a part of the fluid creating a negative flux (to the right) during the recovery stroke (see the instantaneous flux at the instants corresponding to Figure 4c−e in Figure 4g). The evolution of flow during the beat cycle can be seen from the total flow created, shown using the solid line in Figure 4g. The flow created by the cilia increases during the

effective stroke and decreases during the recovery. Because of the asymmetric motion of the cilia, the positive flow is larger than the negative flow; the cilia create a net flow of magnitude 0.2 (in units of πL2/2) to the left by the end of the cycle (t = tbeat). In the process of creating this flow, the cilia push the fluid back and forth, which leads to a small value of the effectiveness (0.2). The position of the fluid particles (which formed a straight line at the beginning of the beat cycle) at the end of the beat cycle provides a Lagrangian perspective of the fluid transported. The recovery stroke takes place with a whiplike motion of the cilia in a short duration of time, which leads to a larger instantaneous fluid velocity and flux during the recovery stroke than during the effective stroke (compare instants 4 and 5 with 1 and 2 in Figure 4). Because the Reynolds number is low, the kinetic energy input by the cilia is instantaneously dissipated by the large viscous forces. Therefore, even if the cilia impart a large velocity to the fluid during the whiplike recovery stroke (Figure 4d,e), the velocity and flux drop to zero as soon as the cilia have returned to the initial position (Figure 4f and instant 6 in Figure 4g). It is to be noted that because of the small viscous forces (compared to elastic 7928


Langmuir

Article

inertia number, and magnetic number in the limit of low Reynolds number was studied in ref 16. It was shown that the inertia number does not play a major role and when the fluid number is increased (for a given magnetic number) the asymmetry (and therefore the flow) decreases. Hence, in this section we study the flow through the microfluidic channel as a function of the fluid number and the Reynolds number for a given magnetic number. The area flow is plotted as a function of the Reynolds number (Re) for a fluid number of 0.015 and a magnetic number of 10.89 in Figure 5. As the Reynolds number is increased, the net flow changes direction so that the fluid propelled is in the direction of the recovery stroke (dashed box). On increasing the Reynolds number further, the fluid flow reverses direction again and flow occurs in the direction of the effective stroke. It is to be noted that the flow created at high Reynolds number is larger than that in the Stokes regime. To segregate the contributions to the net fluid flow, we also plot the positive and negative flow created in Figure 5. When Re is increased, the positive flow remains nearly constant but the negative flow increases and causes the area flow to decrease. With a further increase in Re, the negative flow reaches a maximum, after which the positive and negative flows start decreasing. At high Re, the positive flow starts increasing again. We now investigate what mechanisms are involved that drive the flow direction reversal at moderate Reynolds number (Re ≈ 1) and what causes a higher flow at high Re compared to the Stokes regime. Before proceeding any further, we restate the importance of fluid inertia. At high Reynolds numbers, the viscous forces are relatively low. The energy input into the fluid during the effective and recovery strokes is not instantly dissipated but is carried over to the next recovery and effective strokes (effect of Re). Also, at high Reynolds numbers, the diffusion from the cilia to the fluid takes place over a finite interval of time (effect of DH). As a result, the fluid velocity field does not instantaneously follow the velocity of the cilia.

forces) the cilia complete the beat cycle in 0.89tbeat. The cilia remain idle in this “dead” position until the start of the next cycle. 4.2. Effect of the Reynolds Number. In this section, we characterize the response of the microfluidic actuator as a function of the dimensionless numbers that govern the behavior of the system. The flow through the microfluidic channel depends on the inertia number In, fluid number Fn, magnetic number Mn, and Reynolds number Re (section 3). The effect of the fluid number,

Figure 5. Flow at the end of one cycle as a function of the Reynolds number (Re) for Fn = 0.015 and Mn = 10.89. Fluid propelled by the cilia in the direction of the effective and recovery strokes is termed positive (Qp) and negative (Qn) flows, respectively (Figure 2). Different regions of significance are marked I−III. In region I, the large kinetic energy input during the recovery stroke leads to a negative flow during the dead position. In region II, the negative flow decreases because of the delayed momentum diffusion, and the positive flow decreases because the effective stroke first has to cancel the negative flow. In region III, the positive flow increases because of the continuous velocity in the direction of the effective stroke.

Figure 6. (a−e) Contours of absolute velocity (normalized with L/tbeat) at different time instants. The direction of the velocity is given by the streamlines, and the white arrows represent the magnetic field at the respective time instants. The circles represent fluid particles. The parameters used are Fn = 0.015, Mn = 10.89, and Re = 1. Four unit cells are shown for clarity. The legend for the contours is the same as that in Figure 4e. Instantaneous flux and accumulated flow as a function of time. The time instants corresponding to a−d are duly marked. Note that the cilia have completed one beat cycle while the magnetic field rotates by π. 7929


Langmuir

Article

At Re = 1, the velocity diffusion time tdiff is finite (effect of DH, see section 3). The velocity of the fluid does not follow the cilia velocity instantaneously. This causes a negative velocity only near the bottom half of the channel during the recovery stroke (Figure 6b−d) and a consequent reduction of the maximum flux transported at Re ≈ 1 (compare instant 3 in Figure 6e against instant 5 in Figure 4g). Because the fluid number is small (Fn = 0.015), the cilia return to the initial position before t = tbeat and remain idle in the dead position (Figure 6c,d) until the start of the next cycle. Because the viscous forces are comparable to the inertia forces, the kinetic energy input into the fluid during the whiplike recovery stroke is not dissipated instantly but creates a large negative flow during the dead position (after the recovery is complete and before the start of the next cycle) and during the initial part of the effective stroke (Figure 6d,e). This causes an increase in the negative flow during the recovery stroke compared to that in the Stokes regime, despite the lower instantaneous flux. Consequently, during the first half of the subsequent effective stroke, the cilia have to spend energy in canceling the negative flow instead of creating a positive flow. Hence, the positive flow is

smaller than that in the Stokes regime. This leads to flow direction reversal at moderate Reynolds numbers. Next, we analyze the flow behavior at Re = 10 (Figure 7). At Re = 10, the viscous forces are so small compared to the inertia forces that they do not dissipate the energy input into the fluid during the effective stroke. This can be seen from the continuous fluid velocity in the direction of the effective stroke near the top of the channel (contour plots of the absolute velocity in Figure 7). This leads to an increase in the positive flow. It is to be noted that because of resistance exerted by the fluid inertia forces on the cilia, they return to the initial position at t ≈ tbeat, reducing the idle time in the dead position to zero (Figure 7f). The effect of the increased tdiff is to restrict the fluid velocity in the negative direction to a narrow region close to the cilia (Figure 7e,f). This reduces the maximum negative flux even further (compare instant 3 in Figure 6e against instant e in Figure 7g). In this situation (low maximum negative flux and the absence of the dead position), the high energy input during the recovery stroke cannot create a substantial flow because the next effective stroke has already started and effectively blocks the negative flow. This leads to a small negative flow compared to

Figure 7. (a−f) Contours of absolute velocity (normalized with L/tbeat) at different time instants. The direction of the velocity is given by the streamlines, and the white arrows represent the magnetic field at the respective time instants. The circles represent fluid particles. The parameters used are Fn = 0.015, Re = 10, and Mn = 10.89. Four unit cells are shown for clarity. (g) Instantaneous flux and accumulated flow as a function of time. The time instants corresponding to a−f are duly marked. Note that the cilia have completed one beat cycle while the magnetic field rotates by π. 7930


Langmuir

Article

Figure 8. Parametric area flow and effectiveness as a function of Fn and Re at a magnetic number of 10.89.

that of Re = 1. The total flow created by the cilia at the end of the cycle is 0.32 in units of πL2/2 (Figure 7g), which is higher than that in the Stokes regime (0.2). This can also be observed from the position of the fluid particles at the end of the beat cycle. The fluid particles have undergone a larger displacement (Figure 7f) compared to those in the Stokes regime (Figure 4f). The observations in this section can be summarized as follows (Figure 5). The quantity of negative flow depends on the competition between the inertia-induced negative flow and its obstruction by the subsequent effective stroke. At small Re, the negative flux is large but occurs only over a short duration of time (Figure 4g). When Re increases, inertia in the fluid causes the flow to be more local to the cilia but to be sustained over a longer period of time so that the total negative flow increases (Figures 5 and 6). The net effect is that the total flow decreases and becomes negative (i.e., in the recovery direction, see region I of Figure 5). When Re further increases (i.e., region II in Figure 5), the positive flow decreases because a part of the effective stroke is spent in overcoming the negative flow caused by the previous recovery stroke. The negative flow also decreases because of the localization of the negative velocity near the cilia during a recovery stroke. Finally, in region III, inertia is so dominant that the flow becomes almost unidirectional with the total flux being positive during most of the cycle, showing only a small dip during recovery (Figure 7). In region III, the flow continues to increase with inertia. To include the effect of the fluid number (Fn), the flow and effectiveness are calculated for a range of Fn and Re) values in Figure 8. The contour plots suggest three regions of significance, namely, A−C. Region A represents the Stokes regime where the fluid viscous forces dominate the inertia forces. In region A, the effectiveness of the fluid propulsion remains constant at ∼0.2 (Figure 8b). The fluid that is propelled decreases when the fluid number is increased; the increase in the fluid number increases the viscous forces, which decrease the area swept by the cilia tip and therefore the flow (Figure 3a). The reader is referred to ref 16 for a more detailed discussion. In region B, the flow is negative and thus takes place in the direction of the recovery stroke. If either the Reynolds number or the fluid number is increased, the fluid flows in the direction of the effective stroke. In region C, the fluid inertia forces are large and the effectiveness of fluid propulsion is high. In this regime, for low fluid numbers the fluid transported is large compared to the Stokes regime. However, when the fluid number is large, the area swept is smaller (because of large viscous forces) and the flow created is comparable to that in the Stokes regime. Interestingly, the effectiveness of the fluid propulsion is large even in

this region, thus creating unidirectional flow. This suggests that we can use the inertia forces in the fluid to generate a unidirectional flow with periodically beating artificial cilia. In natural ciliary systems, such unidirectional fluid flow is achieved by the hydrodynamic interaction between individual cilia, which leads to a metachronal wave.42 In artificial ciliary systems, however, the metachronal wave can be established by magnetically forcing the cilia to beat out of phase (using unit cells that include multiple cilia).43 The study performed in the current section also suggests that a given ciliary system can be made to switch the flow direction depending on the operating frequency. (See the dashed arrow in Figure 8a, which shows the direction of increasing frequency or decreasing tbeat.) A ciliary system in region B creates a flow in the direction of the recovery stroke. If the frequency of magnetic actuation is increased, both the fluid number and the Reynolds number will increase, resulting in a strong increase in the effectiveness and large unidirectional flow, that is maximal for the range of Fn and Re studied (Figure 8).

5. CONCLUSIONS In this work, we have analyzed the fluid flow created by magnetically driven artificial cilia in microfluidic channels when the Reynolds number is finite. Such magnetic actuators can mimic the asymmetric motion of natural cilia when subjected to a rotating magnetic field. The competition among the fluid viscous forces, inertia forces, and elastic forces in the cilia was captured using Fn, which is the ratio between the fluid drag forces and elastic forces, and Re, the ratio of fluid inertia to viscous forces. The fluid inertia has two effects. First, the momentum diffusion from the cilia at the boundary into the channel is delayed. Second, the energy input into the fluid in the effective stroke is retained beyond its duration. When the Reynolds number is high, this leads to unidirectional flow throughout a beat cycle. At low fluid numbers (Fn < 1), the flow created at high Reynolds numbers (Re > 1) is much larger than that in the Stokes regime. However, at moderate Reynolds numbers, the fluid flow changes direction and occurs in the direction of the recovery stroke. In these cases, because the fluid number is low, the cilia complete the recovery stroke quickly and are idle for a short period of time before the start of the next effective stroke. During this idle position, the high energy input into the fluid by the fast recovery stroke creates a large flow in the direction of the recovery stroke. This recovery flow overcomes the effective flow and creates a net transport in the direction of the recovery stroke. By exploring the parameter space, we found that the flow created by a given ciliary system can switch direction simply by tuning the frequency. 7931


Langmuir

Article

From a design point of view, it is of interest to explore how the geometry of the cilia (length and thickness), magnetic field, cilia spacing, and frequency should be tuned for a given set of cilia material properties, channel height, and fluid viscosity. Because we want the cilia to operate under moderate pressure heads, we take the length of the cilia to be half that of the channel height.17 In section 4.2, we find that the cilia create large unidirectional flows for Re > 1 and Fn < 1. Using the condition that the Reynolds number should be larger than 1, we can find the frequency of operation. Using the condition on Fn, we can arrive at the ratio of the cilia length to its thickness (aspect ratio). Once we know the aspect ratio, we can find the magnetic couple and magnetic field to be applied using the magnetic number (Mn ≈ 10).

■

APPENDIX A: SOLID DYNAMICS MODEL The artificial cilia are modelled as a discrete assemblage of Euler−Bernoulli beam elements. The beam element has two nodes at the ends and three degree of freedom (DOF) at each node: u is the axial displacement, v is the transverse displacement, and φ = ∂v/∂x is the rotation, where x is the axial coordinate along the beam. The axial displacement along the beam (or film) is interpolated linearly, and the transverse displacement is interpolated cubically, u = Nup, v = Nvp where Nu and Nv are the standard interpolation matrices38 and p = {u1 v1 ϕ1 l0 u2 v2 ϕ2 l0}T, where subscripts 1 and 2 refer to the node numbers and l0 is some reference length. The nonlinear axial strain ε in the beam is 2

ε=

∂u ∂ 2v 1 ⎛ ∂v ⎞ + ⎜ ⎟ − y 2 ≡ ε̅ − y χ ∂x 2 ⎝ ∂x ⎠ ∂x

The principle of virtual work is used as a condition to establish equilibrium:44 when a consistent virtual displacement field is applied to a body, the body will be in equilibrium when the virtual work done by the internal forces equals the virtual work done by the external forces, t +Δt t +Δt δWint = δWext

with t +Δt = δWint

∫V

(σδε + ρ(u δ̈ u + v δ̈ v)) dV

0

where δ(·) represents the variation of a quantity, σ is the axial stress, ρ is the density of the film and (·)̈ implies a second derivative with respect to time. By substituting the strains and defining ∫ σ dA = P and −∫ σy dA = M (A = bh is the area of the cross section, h is the thickness, and b is the out-of-plane width of the film), the internal virtual work at time t + Δt can be written as the sum of an elastic part and an inertia part t +Δt = δWint

∫x

(Pt +Δt δε̅ t +Δt + M t +Δt δχ + ρA(uẗ +Δt δu + v ẗ +Δt δv)) dx 0

We now expand the elastic part of the internal work linearly in time by substituting (Qt+Δt = Qt + ΔQ) for the field parameters in the above equation, which gives t +Δt = δWint

∫x

[(Pt δε̅ t + M t δχ) + (ΔP δε̅ t + ΔM δχ) + Pt Δδε̅ + ρA(uẗ +Δt δu + v ẗ +Δt δv)] dx 0

in which terms of order higher than 1 are neglected. The following notations are introduced for convenience: ∂u/∂x = Bup, ∂v/∂x = Bvp, ∂2u/∂x2 = Cvp. The constitutive relations are ΔP = EAΔε̅ and ΔM = EIΔχ, with E = E̅/(1 − ν2) being the effective elastic modulus, ν being the Poisson ratio, I being the second moment of area defined as I = bh3/12, and E̅ being the elastic modulus. By choosing the domain of integration to be the current configuration (i.e., using an updated Lagrangian framework), the total displacements are zero, p = 0, and we get t +Δt = δpT f t + δpT KΔp + δpT Mpt +Δt δWint ̈ int

where f tint =

∫ [Pt BuT + Mt CvT ] dx

is the nodal internal force vector, K=

∫ EA BTu Bu dx + ∫ EI CTv Cv dx + ∫ Pt BTv Bv dx

is the stiffness matrix, the first two terms of which represent the material stiffness and the third term represents the geometric stiffness, and M=

∫ ρA(NTu Nu + NTv Nv) dx 7932


Langmuir

Article

is the mass matrix. The corresponding external virtual work is t +Δt δWext =

⎛

⎞

∫ ⎝f xt+Δt δu + f yt+Δt δv + Nzt+Δt ∂δ∂xv ⎠A dx + ∫ (txt+Δt δu + t yt+Δt δv)b dx ⎜

= δpT

⎟

∫ [(f xt+Δt NTu + f yt+Δt NTv + Nzt+Δt BTv )A + b(txt+Δt NTu + t yt+Δt NTv )] dx

t +Δt = δpT f ext

where f x and f y are the body forces in the axial and transverse directions (for example, magnetic body forces), Nz is the body couple in the out-of-plane direction (for example, the magnetic body couple), and tx and ty are the surface tractions (for example, fluid drag). By equating the internal and external virtual work and noting that the resulting equation holds for arbitrary δp, we get + Δt − f t KΔp + Mp̈ t +Δt = f text int The motion of the film with time is obtained by solving the above equation using Newmark’s algorithm with appropriate initial and boundary conditions.

■

APPENDIX B: DISCRETIZATION OF VARIOUS TERMS USED IN SECTION 2.2

∫V σijδDij dV

=

∫V (−pδij + 2μDij)δDij

=

∫V ⎜⎜−pδij + μ⎜⎜ ∂xi

⎛

⎛ ∂u ⎝

⎝

= μδUI

∫V

+

j

∂uj ⎞⎞ ∂δui ⎟⎟ dV ∂xi ⎟⎠⎟⎠ ∂xj

∂ϕJj ⎞ ∂ϕIi ⎛ ∂ϕJi ∂ϕIi ⎜ ⎟ dV UJ − δUI + ψ dV PJ ⎜ ⎟ ∂xj ⎝ ∂xj ∂xi ⎠ V ∂xj J

∫

UU U + K UPP ) = δUT (KUU U + KUPP) = δUI (KIJ J IJ J

∫V

⎛ ∂u t +Δt duit +Δt ∂u t +Δt t +Δt ⎞⎟ δui dV = ρ δui ⎜⎜ i + i uj ⎟ dV dt ∂ ∂xj t V ⎝ ⎠

∫

⎛ u t +Δt − u t ⎞ ∂(uit + Δui) t i i + (uj + Δuj)⎟⎟ dV Δ ∂ t x j ⎝ ⎠

=ρ

∫V δui⎜⎜

=ρ

∫V δui⎜⎜

=ρ

∫V δui⎜⎜

=ρ

∫V δui⎜⎜

=ρ

∫V δUI ϕIi(U Jt+Δt ΔJit

⎛ u t +Δt − u t ∂u t ∂u t ∂Δui t ⎞⎟ i i + i Δuj + i utj + uj V Δt ∂xj ∂xj ∂xj ⎟⎠ ⎝ ⎛ u t +Δt − u t ∂u t ∂u t ∂u t +Δt t ⎞⎟ i i − i utj + i utj +Δt + i u j ⎟ dV Δt ∂xj ∂xj ∂xj ⎝ ⎠ ⎛ u t +Δt i

⎝ Δt

∂u t ∂u t ∂u t +Δt t ⎞⎟ ut − i − i utj + i utj +Δt + i u j ⎟ dV Δt ∂xj ∂xj ∂xj ⎠ ϕ

= ρδUI

⎛ ut ∂ϕJi ∂u t ⎞ ∂u t − ⎜⎜ i + i utj ⎟⎟ + i ϕJjU Jt +Δt + U Jt +Δt utj ) dV ∂xj ⎠ ∂xj ∂xj ⎝ Δt

∫V Δ1t ϕIiϕJi dV U Jt+Δt + ρδUI ∫V

− ρδUI

⎛ ut

∫V ϕIi⎜⎜ Δit ⎝

+

∂ϕJi ∂uit ϕIiϕJj dV U Jt +Δt + ρδUI utj ϕIi dV U Jt +Δt ∂xj ∂xj V

∫

∂uit t ⎞⎟ u j dV ∂xj ⎟⎠

1 t +Δt + δ = δUIM̂ IJ U Jt +Δt + δUIKIJ UJ UIKIJ2 U Jt +Δt − δUIFI

̂ t +Δt + δUT K1Ut +Δt + δUT K2Ut +Δt − δUT F = δUT MU ̂ + δUT K1U + δUT K2U − δUT F = δUT MU

7933


Langmuir

∑ δ(λiJ(uiJ f

Article

− pi̇ J )) = λ iJ δuiJ + δλ iJ (uiJ − pi̇ J )

J

= λ iJ ϕ IiJδUI + δλ iJ (ϕ IiJUI − pi̇ J ) = δUT ∑ ϕ J λJ + δλJT (ϕ JT U − AJ ṗ J ) J

■

APPENDIX C: MAGNETIC FIELD CAUSED BY A MAGNETIC SEGMENT The magnetic field caused by a rectangular magnetic segment is ⎛ ⎡ h − M̂ x ⎜ 1⎢ 2 − ̂ Hx = ⎜ −tan ⎢ l 2π ⎜ ⎢⎣ − ⎝ 2 ⎛ ⎡ ̂ My 1 ⎜ ln⎢⎛⎜ − l + 2π ⎜⎝ 2 ⎣⎢⎝ 2 M̂ y ⎛ 1 ⎡⎛ l ⎜ ln⎢⎜ − − 2π ⎜⎝ 2 ⎢⎣⎝ 2

⎛ ⎡ h ⎤ ⎡h ⎤⎞ ⎤ ⎡ h ⎤⎞ − y ̂ ⎥⎟ − − ŷ ⎥ − y ̂ ⎥⎟ − ŷ ⎥ ̂ ⎜ ⎢ ⎢ ⎢ M x −1 2 −1 −1 2 2 ⎜ −tan ⎢ l ⎥ + tan ⎢ l ⎥⎟ ⎥ + tan ⎢ l ⎥⎟ − π 2 ⎟ ⎜ + x ̂ ⎥⎦ + x ̂ ⎥⎦ ⎢⎣ ⎢⎣ + x ̂ ⎥⎦⎟ ⎢⎣ − + x ̂ ⎥⎦ ⎠ ⎠ ⎝ 2 2 2 2 2⎤ 2 2 ⎤⎞ ⎡ ⎞ ⎛ h ⎞ ⎛l ⎞ ⎛ h ⎞ 1 − x⎟̂ + ⎜ − + y ̂⎟ ⎥ − ln⎢⎜ − x⎟̂ + ⎜ − + y ̂⎟ ⎥⎟⎟ ⎠ ⎝ 2 ⎠ ⎥⎦ ⎝ ⎠ ⎝ ⎠ ⎥⎦⎠ 2 ⎢⎣ 2 2 ⎞2 ⎛ h ⎞ 2 ⎤ 1 ⎡⎛ l ⎞2 ⎛ h ⎞ 2 ⎤⎞ − x⎟̂ + ⎜ + y ̂⎟ ⎥ − ln⎢⎜ − x⎟̂ + ⎜ + y ̂⎟ ⎥⎟⎟ ⎠ ⎝2 ⎠ ⎥⎦ ⎠ ⎝2 ⎠ ⎥⎦⎠ 2 ⎢⎣⎝ 2

⎛ ⎛ ⎡ l ⎤ ⎡ l ⎤⎞ ⎡ l ⎤ ⎡l ⎤⎞ M̂ y ⎜ M̂ y ⎜ ⎢ − 2 − x̂ ⎥ ⎢ 2 − x ̂ ⎥⎟ ⎢ − 2 − x̂ ⎥ ⎢ 2 − x ̂ ⎥⎟ 1 1 1 1 − − − − Ĥ y = ⎜ − tan ⎢ h ⎜ − tan ⎢ h ⎥ + tan ⎢ h ⎥⎟ − ⎥ + tan ⎢ h ⎥⎟ 2π ⎜ 2π ⎜ + y ̂ ⎥⎦ ⎢⎣ − + y ̂ ⎥⎦ ⎢⎣ − + y ̂ ⎥⎦⎟ ⎢⎣ ⎢⎣ + y ⎥⎦⎟ ⎝ ⎠ ⎝ ⎠ 2 2 2 2 ⎛ 2 2⎤ 2 2 ⎤⎞ ⎡ ⎡ ̂ ⎞ ⎛ h ⎞ ⎛ l ⎞ ⎛h ⎞ M 1 ⎛ l 1 + x ⎜⎜ ln⎢⎜ − + x⎟̂ + ⎜ − − y ̂⎟ ⎥ − ln⎢⎜ − + x⎟̂ + ⎜ − y ̂⎟ ⎥⎟⎟ ⎠ ⎝ 2 ⎠ ⎥⎦ ⎠ ⎝2 ⎠ ⎥⎦⎠ 2π ⎝ 2 ⎢⎣⎝ 2 2 ⎢⎣⎝ 2 ⎛ ⎡⎛ ⎞2 ⎛ ⎞2 ⎤ ⎞2 ⎛ ⎞ 2 ⎤⎞ − M̂ x ⎜ 1 ⎡⎛ l ⎢⎜ + x⎟̂ + ⎜ − h − y ̂⎟ ⎥ − 1 ln⎢⎜ l + x⎟̂ + ⎜ h − y ̂⎟ ⎥⎟ + ln ⎠ ⎝ 2 ⎠ ⎥⎦ ⎠ ⎝2 ⎠ ⎥⎦⎟⎠ 2π ⎜⎝ 2 ⎢⎣⎝ 2 2 ⎢⎣⎝ 2

where M̂ x and M̂ y are magnetizations in the tangential (or length) and normal (or thickness) directions, h is the thickness, and l is the length of the segment. Here, x̂ and ŷ are the local coordinates having their origin in the center of the segment.

■

APPENDIX D: VALIDATION OF THE MAGNETOSTATIC MODEL In this Appendix, we validate the magnetostatic model used in this article on the basis of a superparamagnetic film whose magnetic properties are anisotropic. The contour length of the film is 100 μm, and its thickness is 2 μm. The film has susceptibilities of 4.6 and 0.8 in the tangential and normal directions, respectively. Six configurations of the film are chosen to validate the magnetostatic model. See Figure 1A. The configurations chosen reflect the deformed geometry of the film during one cycle of ciliary motion.16 The solution for the magnetic flux density B is compared to the solution from a commercial multiphysics software program (COMSOL). The field values are

Figure 1A. Positions of the film used to test the magnetostatic model. The hollow arrows show the direction of the applied field for the respective positions. 7934


Langmuir

Article

Figure 2A. Flux density comparison. The symbols represent the solution from COMSOL, and the lines represent the field calculated from the present model.

normalized with the absolute value of the applied field. The comparison is shown only for positions 2−6 in Figure 2. The solid lines represent the solution from our magnetostatic model, and the + symbols represent the solution from COMSOL. The agreement is excellent. 7935


Langmuir

■

Article

APPENDIX E: EFFECT OF THE DIFFUSION NUMBER WITH RESPECT TO THE ANALYTICAL MODEL In this Appendix, we study the effect of the diffusion number DH on the propelled fluid. To this end,45 we analyse the case when the height of the channel is larger than the cilia length. This allows us to assume that the cilia create an oscillating Couette flow (shear flow) such that the velocity at the bottom boundary of the channel (f(t)) has an oscillating component (u0) and a steady component (umean) f (t ) = u0 sin wt + u mean(1 − e−t / t ref ) where ω = 2π/tref is the angular frequency. The solution to the problem of oscillating shear flow can be obtained by solving the 1D Navier− Stokes equation. The 1D Navier−Stokes equation (where in a shear flow the pressure is an arbitrary constant) is ∂u ∂ 2u =μ 2 ∂t ∂y where u is the fluid velocity. The initial and boundary conditions for the Couette flow are ρ

u(0, t ) = f (t ), u(H , t ) = 0, (y , 0) = 0 The solution to the 1D Navier−Stokes equation is given by46 ∞

u(y , t ) =

μt 2n π ref2 ρH n=1

∑

∫0

t

⎛ μt n2 π 2(t − τ) ⎞ dτ ⎛ n πy ⎞ ref ⎟ exp⎜ − f (τ) sin⎜ ⎜ ⎟⎟ ⎝ H ⎠ t ref ⎝ ρH 2 ⎠ t ref

Equation 1 shows that the velocity of the fluid at any time and the location depend on the dimensionless parameter DH = ρH2/μtref. The velocity u(y, t) exhibits transient behavior for two reasons: (i) the inertia of the fluid and (ii) the transient part of the boundary velocity associated with umean. The inertial transience is represented by the exponential term in eq 1, whose duration scales with ρH2/μ. The transient period due to the boundary condition is represented by the exponential term in f(t), and its duration is 10tref. Because the boundary velocity is oscillatory, the steady state will also be oscillatory in nature. The instantaneous flux Q(t) at any time t is Q(t) = ∫ 0H u(y) dy, and the fluid transported is q̅ = ∫ cycle [Q(t) dt]/tref. The fluid transported (q)̅ after the system has reached a steady state (t > ρH2/μ) can be calculated analytically and is given by q̅ = Humean/2. Thus, the mean flow rate depends only on umean and scales linearly with the channel height but is independent of DH. DH governs only the duration of the transient period.

■

(12) Bruss, H.; Brask, A.; Kutter, J. P. Nanofluidic components for Electrokinetic micropumps. 2004. (13) Sleigh, M. A., Ed. Cilia and Flagella; Academic Press: London, 1974. (14) Fahrni, F.; Prins, M. W. J.; van IJzendoorn, L. J. Lab Chip 2009, 9, 3413−3421. (15) Belardi, J.; Schorr, N.; Prucker, O.; Ruhe, J. Adv. Funct. Mater. 2011, 21, 3314−3320. (16) Khaderi, S. N.; Baltussen, M. G. H. M.; Anderson, P. D.; Ioan, D.; den Toonder, J. M. J.; Onck, P. R. Phys. Rev. E 2009, 79, 046304. (17) Khaderi, S. N.; Craus, C. B.; Hussong, J.; Schorr, N.; Belardi, J.; Westerweel, J.; Prucker, O.; Ruhe, J.; den Toonder, J. M. J.; Onck, P. R. Lab Chip 2011, 11, 2002−2010. (18) Oh, K.; Chung, J.-H.; Devasia, S.; Riley, J. J. Lab Chip 2009, 9, 1561−1566. (19) den Toonder, J.; Bos, F.; Broer, D.; Filippini, L.; Gillies, M.; de Goede, J.; Mol, T.; Reijme, M.; Talen, W.; Wilderbeek, H.; Khatavkar, V.; Anderson, P. Lab Chip 2008, 8, 533−541. (20) van Oosten, C. L.; Bastiaansen, C. W. M.; Broer, D. J. Nat. Mater. 2009, 8, 677−682. (21) Hussong, J.; Schorr, N.; Belardi, J.; Prucker, O.; Ruhe, J.; Westerweel, J. Lab Chip 2011, 11, 2017−2022. (22) Shields, A. R.; Fiser, B. L.; Evans, B. A.; Falvo, M. R.; Washburn, S.; Superfine, R. Proc. Natl. Acad. Sci. U.S.A. 2010. (23) Evans, B. A.; Shields, A. R.; Carroll, R. L.; Washburn, S.; Falvo, M. R.; Superfine, R. Nano Lett. 2007, 7, 1428−1434. (24) Vilfan, M.; Potocnik, A.; Kavcic, B.; Osterman, N.; Poberaj, I.; Vilfan, A.; Babic, D. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 1844− 1847. (25) Timonen, J. V. I; Johans, C.; Kontturi, K.; Walther, A.; Ikkala, O.; Ras, R. H. A. ACS Appl. Mater. Interfaces 2010, 2, 2226−2230. (26) Kim, Y. W.; Netz, R. R. Phys. Rev. Lett. 2006, 96, 158101. (27) Gauger, E. M.; Downton, M. T.; Stark, H. Eur. Phys. J. E 2009, 28, 231−242.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is a part of the 6th Framework European Project “Artic” under contract STRP 033274. We acknowledge fruitful discussions with Michiel Baltussen and Patrick Anderson.

■

REFERENCES

(1) Laser, D. J.; Santiago, J. G. J. Micromech. Microeng. 2004, 14, R35−R64. (2) Whitesides, G. M. Nature 2006, 442, 368−373. (3) Chang, S. T.; Paunov, V. N.; Petsev, D. N.; Velev, O. D. Nat. Mater. 2007, 6, 235−240. (4) Schilling, E.; Kamholz, A.; Yager, P. Anal. Chem. 2002, 74, 1798− 1804. (5) Jeon, N.; Dertinger, S.; Chiu, D.; Choi, I.; Stroock, A.; Whitesides, G. Langmuir 2000, 16, 8311−8316. (6) Zeng, S.; Chen, C.; Santiago, J. G.; Chen, J.; Zare, R. N.; Tripp, J. A.; S., F.; Frechet, J. M. J. Sens. Actuators, B 2002, 82, 209−212. (7) Chen, L.; Ma, J.; Tan, F.; Guan, Y. Sens. Actuators, B 2003, 88, 260−265. (8) Lemoff, A. V.; Lee, A. P. Sens. Actuators, B 2000, 63, 178−185. (9) Homsy, A.; Koster, S.; Eijkel, J. C. T.; van den Berg, A.; Lucklum, F.; Verpoorte, E.; de Rooij, N. F. Lab Chip 2000, 5, 466−471. (10) Studer, V.; Pepin, A.; Chen, Y.; Ajdari, A. Analyst 2004, 129, 944−949. (11) Wu, J.; Lian, M.; Yang, K. Appl. Phys. Lett. 2007, 90, 234103. 7936


Langmuir

Article

(28) Alexeev, A.; Yeomans, J. M.; Balazs, A. C. Langmuir 2008, 24, 12102−12106. (29) Downton, M. T.; Stark, H. Europhys. Lett. 2009, 85, 44002. (30) Baltussen, M. G. H. M.; Anderson, P. D.; Bos, F. M.; den Toonder, J. M. J. Lab Chip 2009, 9, 2326−2331. (31) Hussong, J.; Breugem, W. P.; Westerweel, J. J. Fluid Mech. 2011, accepted for publication. (32) Khaderi, S. N.; Baltussen, M. G. H. M.; Anderson, P. D.; den Toonder, J. M. J.; Onck, P. R. Phys. Rev. E 2010, 82, 027302. (33) Baaijens, F. P. T Int. J. Numer. Methods Fluids 2001, 35, 743− 761. (34) van Loon, R.; Anderson, P. D.; van de Vosse, F. N. J. Comput. Phys. 2006, 217, 806−823. (35) Annabattula, R. K.; Huck, W. T. S.; Onck, P. J. Mech. Phys. Solids 2010, 58, 447−465. (36) Bathe, K.-J. Finite Element Procedures; Prentice-Hall: Englewood Cliffs, NJ, 1996. (37) Lanczos, C. The Variational Principles of Mechanics; University of Toronto Press: Toronto, 1952. (38) Cook, R. D.; Malkus, D.; Plesha, M.; Malkus, D. S.; Plesha, M. E. Concepts and Applications of Finite Element Analysis, 4th ed.; John Wiley and Sons: New York, 2002. (39) Jackson, J. D. Classical Electrodynamics; John Wiley and Sons: New York, 1974. (40) van Rijsewijk, L. Electrostatic and Magnetic Microactuation of Polymer Structures for Fluid Transport. M.Sc. Thesis, Eindhoven University of Technology, 2006 (41) Childress, S. Mechanics of Swimming and Flying; Cambridge University Press: Cambridge, U.K., 1981. (42) Satir, P.; Sleigh, M. A. Annu. Rev. Physiol. 1990, 52, 137−155, PMID 2184754 (43) Khaderi, S. N.; den Toonder, J. M. J.; Onck, P. R. J. Fluid Mech. 2011, 688, 44−65. (44) Malvern, L. E. Introduction to the Mechanics of a Continuous Medium; Prentice-Hall: Englewood Cliffs, NJ, 1977. (45) Blake, J. R.; Liron, N.; Aldis, G. K. J. Theor. Biol. 1982, 98, 127− 141. (46) Polyanin, A. D. Handbook of Linear Partial Differential Equations for Engineers and Scientists; Chapman & Hall/CRC: Boca Raton, FL, 2002; section 1.1.2−1.1.5.

7937


Magnetically Actuated Artificial Cilia: The Effect of Fluid Inertia

Recommend Documents