Magnetically Actuated Cilia Array
Based on Gu et al. (2020)
An \(8 \times 8\) array of magnetically actuated cilia beats in an antiplectic metachronal wave pattern, pumping fluid along a channel. The cilia magnetization varies along the array, with cilia at the same \(x\) position sharing the same dipole orientation, and the array is driven by a rotating external magnetic field. The phase offset between neighboring cilia produces a traveling wave whose wavelength \(\lambda\) controls the pumping efficiency.
Parameters
| Parameter | Value |
|---|---|
| Cilium length \(L\) | 4 mm |
| Cilium diameter \(d\) | 0.8 mm |
| Young's modulus \(E\) | 185 kPa |
| Shear modulus \(G\) | 61.6 kPa |
| Solid density \(\rho_s\) | 2.39 g/cm\(^3\) |
| Number of bonds \(N\) | 10 (bond length \(\ell_0 = 0.4\) mm) |
| Magnetization \(M\) | 20 kA/m |
| Applied field \(B\) | 80 mT |
| Field angular velocity \(\omega_B\) | −30 deg/s |
| Fluid kinematic viscosity \(\nu_f\) | 0.913 mm\(^2\)/ms |
| Fluid density \(\rho_f\) | 1.26 mg/mm\(^3\) |
| Domain (\(x \times y \times z\)) | 90 mm \(\times\) 36 mm \(\times\) 45 mm |
| Metachronal wavelength \(\lambda\) | varies |
Setup
The cilia are immersed in a viscous fluid (99% glycerol) and driven by a magnetic field rotating in the \(x\)–\(z\) plane:
\[\mathbf{B}(t) = B\bigl[\cos\phi(t)\,\hat{\mathbf{x}} + \sin\phi(t)\,\hat{\mathbf{z}}\bigr], \qquad \phi(t) = \frac{\pi}{2} + \omega t.\]
The magnetization pattern varies with position and is defined in the generated data file. One thousand tracer particles are randomly distributed in the channel to visualize the flow and measure the net fluid transport per cycle.
CoupLB
The fluid is periodic in \(x\), bounded by no-slip walls in \(y\), and no-slip at the bottom (\(z = 0\)) with free-slip at the top.
The LBM grid has spacing \(\Delta x_\text{LB} = 1.0\) mm (\(90 \times 36 \times 45\) nodes).
The BPM bond length is shorter than the grid spacing (\(\mathrm{Gr} = \Delta x_\text{BP}/\Delta x_\text{LB} = 0.4\)).
No subcycling is used; the LAMMPS and LBM timesteps are both 0.008 ms.
The data file generator includes a callable function lbm_analysis(dt) to check LBM and BPM stability.
Key LAMMPS settings
- Units:
units ljusing our milli unit system - Atom style:
atom_style hybrid molecular bpm/sphere dipole - Boundary:
p f f— periodic in \(x\), fixed walls in \(y\) and \(z\) - Pair style:
pair_style zero— no contact - Magnetic actuation:
fix efield v_bx 0 v_bzapplies the rotating field - FSI:
fix couplbwithwall_y 1 1andwall_z 1 3 - Time integration:
fix nve/bpm/sphere update dipolefor cilia,fix nvefor tracers - Confinement:
fix wall/regionto prevent lost atoms - Damping: no explicit damping — all dissipation comes from the fluid coupling
Atom types and molecule IDs
The atom types are: 1 (plate), 2 (cilia), and 3 (tracers). The molecule IDs are: 1 (floor = plate + cilia base), 2 (free cilia atoms). Each cilium has \(N+1\) atoms, where the base atom is excluded from the dynamics. The plate atoms are only for visualization.
Files
Data file generator
The generator writes a data file for the \(8 \times 8\) cilia array using the normalized wavelength set by the lamda variable.
Call lbm_analysis(dt) to check LBM and BPM stability for a given timestep.
create_magnetic_cilia_array.py ⬇ Download
import numpy as np
lamda = 1 # metachronal wavelength (in units of Ncx * cilia_spacing)
# Physical parameters (milli units: mm, ms, mg, mu0 = 4pi)
L_plate = 36. # plate side length
W_plate = 36.
t_plate = 0.8 # plate thickness
L = 4. # cilium length
D = 0.8 # cilium diameter
rho_s = 2.39 # solid density
E = 185. # Young's modulus (kPa)
G = 61.6 # shear modulus (kPa)
nu_s = E / (2 * G) - 1
# Magnetic properties
M_mag = 0.2 # magnetization (1e5 A/m)
B_ext = 8.0 # applied field strength
# Fluid (99% glycerol)
nuf = 0.913 # kinematic viscosity (mm^2/ms)
rho_f = 1.26 # fluid density (mg/mm^3)
# Actuation
omega_B = -30 * np.pi / 180 / 1000 # -30 deg/s in rad/ms
u_phys = L * np.abs(omega_B) # tip speed estimate
# Domain
xbox = 90.
ybox = 36.
zbox = 45.
# Cilia array
Ncx = 8 # cilia in x
Ncy = 8 # cilia in y
Na = 11 # atoms per cilium
cilia_spacing = 4.0 # mm between cilia centers
# Plate discretization
Nx_plate = 45
Ny_plate = 45
# LBM grid
dx_LB = 1.0
# BPM discretization
A_cs = np.pi * D**2 / 4
I_cs = np.pi * D**4 / 64
J_cs = np.pi * D**4 / 32
l0 = L / (Na - 1)
# =========================================================================
# Metachronal wave
# =========================================================================
mu_mag = M_mag * A_cs * l0
def theta_dipole(x_pos, lamda, theta0=0, sign=1):
"""Dipole phase angle as function of x position."""
return sign * (2 * np.pi * x_pos / (lamda * Ncx * cilia_spacing) + theta0)
# =========================================================================
# Create geometry
# =========================================================================
Natoms = 0
V_xyz, V_type, V_molid, V_diam, V_rho, V_mu = [], [], [], [], [], []
bond_atoms = []
# --- Plate atoms (type 1, molecule 1) ---
dx_plate = L_plate / Nx_plate
dy_plate = W_plate / Ny_plate
plate_xyz = []
for i in range(Nx_plate):
for j in range(Ny_plate):
x = dx_plate * (i + 0.5)
y = dy_plate * (j + 0.5)
z = D / 2
plate_xyz.append([x, y, z])
plate_xyz = np.array(plate_xyz)
n_plate = len(plate_xyz)
V_xyz.append(plate_xyz)
V_type.append(np.ones((n_plate, 1)))
V_molid.append(np.ones((n_plate, 1)))
V_diam.append(t_plate * np.ones((n_plate, 1)))
V_rho.append(rho_s * np.ones((n_plate, 1)))
V_mu.append(np.zeros((n_plate, 3)))
Natoms += n_plate
# --- Cilia (type 2, molecule 2) ---
def create_cilium(x0, y0, mol_id=2):
global Natoms
xyz = np.zeros((Na, 3))
xyz[:, 0] = x0
xyz[:, 1] = y0
xyz[:, 2] = np.linspace(0, L, Na) + D / 2
types = 2 * np.ones((Na, 1))
molids = mol_id * np.ones((Na, 1))
molids[0] = 1 # base atom belongs to plate molecule
th = theta_dipole(x0, lamda=lamda, theta0=np.pi/2, sign=1) # antiplectic
mu = np.zeros((Na, 3))
mu[:, 0] = mu_mag * np.cos(th)
mu[:, 2] = mu_mag * np.sin(th)
diams = D * np.ones((Na, 1))
rhos = rho_s * np.ones((Na, 1))
bonds = np.zeros((Na - 1, 2))
for k in range(Na - 1):
bonds[k, 0] = Natoms + k + 1
bonds[k, 1] = Natoms + k + 2
V_xyz.append(xyz)
V_type.append(types)
V_molid.append(molids)
V_diam.append(diams)
V_rho.append(rhos)
V_mu.append(mu)
bond_atoms.append(bonds)
Natoms += Na
for i in range(1, Ncx + 1):
for j in range(1, Ncy + 1):
create_cilium(cilia_spacing * i, cilia_spacing * j)
# Concatenate
V_xyz = np.concatenate(V_xyz, axis=0)
V_type = np.concatenate(V_type, axis=0)
V_molid = np.concatenate(V_molid, axis=0)
V_diam = np.concatenate(V_diam, axis=0)
V_rho = np.concatenate(V_rho, axis=0)
V_mu = np.concatenate(V_mu, axis=0)
bond_atoms = np.concatenate(bond_atoms, axis=0)
Nbonds = bond_atoms.shape[0]
# Center structure in box
x_center = xbox / 2
y_center = ybox / 2
x_shift = x_center - V_xyz[:n_plate, 0].mean()
y_shift = y_center - V_xyz[:n_plate, 1].mean()
V_xyz[:, 0] += x_shift
V_xyz[:, 1] += y_shift
# =========================================================================
# Write data file
# =========================================================================
filename = f'magnetic-cilia-array-8x8-antiplectic-lambda-{lamda}.lam'
# atom_style hybrid molecular bpm/sphere dipole
# Columns: id type x y z molid diameter density q mux muy muz
with open(filename, 'w') as f:
f.write('LAMMPS DATA FILE\n\n')
f.write(f'{Natoms} atoms\n')
f.write(f'{Nbonds} bonds\n\n')
f.write('3 atom types\n')
f.write('1 bond types\n\n')
f.write(f'{0:.6f} {xbox:.6f} xlo xhi\n')
f.write(f'{0:.6f} {ybox:.6f} ylo yhi\n')
f.write(f'{0:.6f} {zbox:.6f} zlo zhi\n\n')
f.write('Atoms\n\n')
for i in range(Natoms):
f.write('%i %i %f %f %f %i %f %f 0 %f %f %f\n' %
(i+1, int(V_type[i, 0]),
V_xyz[i, 0], V_xyz[i, 1], V_xyz[i, 2],
int(V_molid[i, 0]), V_diam[i, 0], V_rho[i, 0],
V_mu[i, 0], V_mu[i, 1], V_mu[i, 2]))
f.write('\n')
f.write('Bonds\n\n')
for i in range(Nbonds):
f.write(f'{i+1} 1 {int(bond_atoms[i, 0])} {int(bond_atoms[i, 1])}\n')
print(f"Written: {filename}")
print(f" {Natoms} atoms ({n_plate} plate + {Ncx*Ncy*Na} cilia)")
print(f" {Nbonds} bonds")
def lbm_analysis(dt):
"""Check LBM and BPM stability for the given timestep (no subcycling)."""
cs2 = 1.0 / 3.0
nx_lb = int(xbox / dx_LB)
ny_lb = int(ybox / dx_LB)
nz_lb = int(zbox / dx_LB)
U_star = u_phys * dt / dx_LB
nu_LB = nuf * dt / dx_LB**2
tau = 0.5 + nu_LB / cs2
Ma = U_star / (1.0 / np.sqrt(3))
Re = u_phys * L / nuf
# BPM stability
Kr = E * A_cs / l0
Ks = 12 * E * I_cs / l0**3
Kt = G * J_cs / l0
Kb = E * I_cs / l0
m_seg = rho_s * A_cs * l0
I_seg = m_seg * l0**2
omega_max = max(np.sqrt(Kr/m_seg), np.sqrt(Ks/m_seg),
np.sqrt(Kt/I_seg), np.sqrt(Kb/I_seg))
dt_crit = np.pi / omega_max
print(f"\nRe ≈ {Re:.4f}")
print(f"LBM grid: {nx_lb} x {ny_lb} x {nz_lb} (dx = {dx_LB})")
print(f"BPM grid: Na = {Na}, l0 = {l0:.4f}")
print(f"\n dt = {dt}")
print(f" U* = {U_star:.6e}")
print(f" nu* = {nu_LB:.6f}")
print(f" tau* = {tau:.6f}", end="")
if tau <= 0.5: print(" *** UNSTABLE ***")
elif tau < 0.55: print(" *** WARNING: marginal (close to 0.5) ***")
elif tau > 2.0: print(" *** WARNING: over-damped ***")
else: print(" OK")
print(f" Ma = {Ma:.6e}", end="")
if Ma > 0.3: print(" *** WARNING: high ***")
else: print(" OK")
print(f" dt_crit = {dt_crit:.6f}")
print(f" dt/dt_crit = {dt/dt_crit:.4f}", end="")
if dt > dt_crit: print(" *** UNSTABLE ***")
else: print(" OK")
if __name__ == '__main__':
lbm_analysis(dt=0.008)
LAMMPS input script
magnetic-cilia-array.in ⬇ Download
# =========================================================================
# Magnetically actuated 8x8 cilia array with antiplectic wave pattern
#
# Replace LAMBDA with the desired wavelength before running.
# =========================================================================
# Cilium parameters (milli units: mm, ms, mg, mu0=4pi)
variable Nb equal 10
variable L equal 4
variable d equal 0.8
variable l0 equal $(v_L/v_Nb)
# Cross-section properties
variable A equal (PI/4)*v_d^2
variable I equal (PI/64)*v_d^4
variable J equal (PI/32)*v_d^4
# Elastic moduli (kPa)
variable E equal 185
variable G equal 61.6
# =========================================================================
# Initialization
# =========================================================================
units lj
dimension 3
boundary p f f
processors 4 3 3
atom_style hybrid molecular bpm/sphere dipole
special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
newton on off
neighbor 0.5 bin
neigh_modify delay 0 every 1
comm_modify mode single cutoff $(2*v_l0) vel yes
read_data magnetic-cilia-array-8x8-antiplectic-lambda-LAMBDA.lam
# =========================================================================
# Walls and tracers
# =========================================================================
# Confinement walls (y and z)
region walls block INF INF 0 36 0 45 side in
# Tracer particles for fluid transport measurement
region dropzone block INF INF INF INF 1 44 side in
create_atoms 3 random 1000 12345 dropzone overlap 1.0
set type 3 diameter 0.2
set type 3 density 1.
# =========================================================================
# Groups
# =========================================================================
group floor molecule 1
group cilia molecule 2
group tracers type 3
group mobile union cilia tracers
# =========================================================================
# BPM bond style
# =========================================================================
variable Kr equal v_E*v_A/v_l0
variable Ks equal 12*v_E*v_I/v_l0^3
variable Kt equal v_G*v_J/v_l0
variable Kb equal v_E*v_I/v_l0
bond_style bpm/rotational smooth no break no
bond_coeff 1 ${Kr} ${Ks} ${Kt} ${Kb} &
-1 -1 -1 -1 &
0 0 0 0
pair_style zero $(2*v_l0)
pair_coeff * *
# =========================================================================
# Magnetic actuation
# =========================================================================
# Rotating field in x-z plane
variable B equal 8
variable omega equal -30*PI/180/1000
variable phi equal PI/2+v_omega*step*v_timestep
variable bx equal v_B*cos(v_phi)
variable bz equal v_B*sin(v_phi)
# =========================================================================
# Dynamics
# =========================================================================
variable timestep equal 0.008
timestep ${timestep}
variable nstepspercycle equal "abs(2*PI/(v_omega*v_timestep))"
variable nsteps equal $(10*v_nstepspercycle)
compute mytemp all temp/sphere
compute erot all erotate/sphere
variable cycle equal step/v_nstepspercycle
thermo $(v_nstepspercycle/100)
thermo_style custom step temp ke c_erot pe v_cycle
thermo_modify temp mytemp norm no
dump 1 all custom $(v_nstepspercycle/50) &
./LAMBDA-array.dump &
id type mol x y z mass diameter mux muy muz
dump 2 tracers custom ${nstepspercycle} &
./LAMBDA-tracers.dump id xu yu zu
dump_modify 2 sort id
# =========================================================================
# Fixes
# =========================================================================
fix nve cilia nve/bpm/sphere update dipole
fix nvetr tracers nve
fix wall mobile wall/region walls harmonic 100 0 1.0
fix bext cilia efield v_bx 0 v_bz
fix lb mobile couplb 90 36 45 0.913 1.26 &
dx 1.0 &
wall_y 1 1 &
wall_z 1 3 &
vtk $(v_nstepspercycle/50) vtk/fluid
run ${nsteps}
References
- Gu, H., Boehler, Q., Cui, H., Secchi, E., Savorana, G., De Marco, C., Gervasoni, S., Peyron, Q., Huang, T.-Y., Pane, S., Hirt, A.M., Ahmed, D., Nelson, B.J. (2020). Magnetic cilia carpets with programmable metachronal waves. Nat. Commun., 11, 2637. doi:10.1038/s41467-020-16458-4