Oscillatory Channel Flow
Based on Pinelli et al. (2017)
Five elastic filaments are immersed in an oscillatory channel flow between no-slip walls.
Lattice Boltzmann and Immersed Boundary Method
The fluid–structure interaction (FSI) is handled by the CoupLB package, a custom, lightweight LBM implementation in LAMMPS. Momentum is exchanged between the Lagrangian BPM particles and the Eulerian fluid lattice through an immersed boundary method (IBM). A penalty force drives each particle toward the local fluid velocity, and the equal and opposite reaction is spread back onto the fluid grid.
Because the stiff BPM springs require a smaller timestep than the LBM, multiple LAMMPS steps are performed per LB step.
This subcycling is controlled by the md_per_lb keyword in fix couplb.
The fluid velocity is interpolated once at the start of each LBM cycle, and the penalty force is recomputed at every substep using each particle's current velocity.
Parameters
| Parameter | Value |
|---|---|
| Filament length \(L\) | 10 mm |
| Filament diameter \(d\) | 1 mm |
| Young's modulus \(E\) | 1.23 MPa |
| Poisson's ratio \(\nu_s\) | 0.3 |
| Filament density \(\rho_s\) | 0.97 g/cm\(^3\) |
| Fluid kinematic viscosity \(\nu_f\) | 0.1 mm\(^2\)/ms |
| Fluid density \(\rho_f\) | 1.26 g/cm\(^3\) |
| Peak velocity \(u_\text{max}\) | 0.6 m/s |
| Oscillation frequency \(f\) | 1 Hz |
| Reynolds number \(\mathrm{Re}\) | 60 |
| Domain (\(x \times y \times z\)) | \(10L \times 6L \times 5L\) |
Setup
Five straight filaments are clamped to the lower wall (\(y = 0\)) of the channel and arranged in a row along \(x\) at intervals of \(L/2\), with the middle filament at the center. The fluid is driven by a sinusoidal body-force acceleration
\[g_x(t) = g_0 \sin(2\pi f\, t),\]
where the amplitude is \(g_0 = u_\text{max} \cdot 2\pi f\), determined from the analytical solution for oscillatory channel flow in the limit where the Stokes layer thickness is much smaller than the half-channel height. No-slip walls bound the channel in \(y\); the remaining directions are periodic.
Discretization
The LBM grid has spacing \(\Delta x_\text{LB} = 0.625\) mm (\(160 \times 96 \times 80\) nodes).
The BPM bond length matches the LBM spacing (\(\mathrm{Gr} = \Delta x_\text{BP}/\Delta x_\text{LB} = 1\)), giving 16 bonds per filament.
The stiff BPM springs require subcycling at 21 LAMMPS steps per LB step (md_per_lb 21).
The data file generator includes a callable function lbm_analysis() that documents how the LBM parameters were selected.
Key LAMMPS settings
- Units:
units ljusing our milli unit system - Atom style:
atom_style hybrid molecular bpm/sphere - Boundary:
p f p— periodic in \(x\) and \(z\), fixed in \(y\) (walls) - Pair style:
pair_style zero— no contact - FSI:
fix couplbwithwall_y 1 1for no-slip walls andgravity v_gx_osc 0 0for the oscillating body force - Time integration:
fix nve/bpm/sphere - Damping: no explicit damping — all dissipation comes from the fluid coupling
Atom types
Three atom types are used: type 1 (clamped base), type 2 (interior), and type 3 (tip). The base atoms are excluded from the dynamics.
Files
Data file generator
The generator writes filaments.lam with 5 filaments (\(\mathrm{Gr} = 1\) by default).
Run lbm_analysis() to print the full LBM parameter analysis.
create_filaments.py ⬇ Download
import numpy as np
Gr = 1.0 # grid ratio dx_BP/dx_LB
# Physical parameters (milli units: mm, ms, mg)
L = 10. # filament length
D = 1. # filament diameter
rho = 0.97 # filament density
E = 1230. # Young's modulus (kPa)
nu_s = 0.3 # Poisson's ratio
G = E / (2 * (1 + nu_s))
nuf = 0.1 # fluid kinematic viscosity (mm^2/ms)
rho_f = 1.26 # fluid density (mg/mm^3)
u_phys = 0.6 # characteristic velocity (mm/ms)
freq = 1e-3 # oscillation frequency (1/ms = 1 Hz)
# Domain
x_length = 10*L
y_length = 6*L
z_length = 5*L
H = y_length # channel height
# Cilia layout
Ncilia = 5
sep = L / 2 # spacing between filaments
# LBM grid spacing
dx_LB = 0.625
# BPM discretization
dx_BP = dx_LB * Gr
Nb = int(np.round(L / dx_BP)) # bonds per filament
Na = Nb + 1 # atoms per filament
l0 = L / Nb
# Cross-section
A = np.pi * D**2 / 4
I = np.pi * D**4 / 64
J = np.pi * D**4 / 32
# =========================================================================
# Write data file
# =========================================================================
x_center = x_length / 2
z_center = z_length / 2
x_positions = x_center + np.linspace(-(Ncilia-1)*sep/2, (Ncilia-1)*sep/2, Ncilia)
Natoms = Ncilia * Na
Nbonds = Ncilia * Nb
# atom_style hybrid molecular bpm/sphere
# Columns: id type x y z molid diameter density
with open('filaments.lam', '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} {x_length:.6f} xlo xhi\n')
f.write(f'{0:.6f} {y_length:.6f} ylo yhi\n')
f.write(f'{0:.6f} {z_length:.6f} zlo zhi\n\n')
f.write('Atoms\n\n')
atom_id = 0
for c in range(Ncilia):
for k in range(Na):
atom_id += 1
x = x_positions[c]
y = D / 2 + k * l0
z = z_center
if k == 0:
atype = 1 # bottom (clamped)
elif k == Na - 1:
atype = 3 # tip
else:
atype = 2 # interior
f.write('%i %i %f %f %f %i %f %f\n' %
(atom_id, atype, x, y, z, c+1, D, rho))
f.write('\n')
f.write('Bonds\n\n')
bond_id = 0
for c in range(Ncilia):
for k in range(Nb):
bond_id += 1
a1 = c * Na + k + 1
a2 = c * Na + k + 2
f.write('%i 1 %i %i\n' % (bond_id, a1, a2))
f.write('\n')
# =========================================================================
# LBM parameter selection and stability analysis
#
# Call lbm_analysis() to print the full report.
# This documents how the parameters in flow.in were chosen.
#
# Option A (Krueger textbook): fix dx_LB, choose U*
# U* = u_phys * dt / dx_LB => dt = U* * dx_LB / u_phys
# nu* = nu * dt / dx_LB^2
# tau* = 0.5 + nu* / cs2 (cs2 = 1/3)
# Ma = U* / cs (cs = 1/sqrt(3))
#
# Subcycling (CoupLB md_per_lb):
# Multiple LAMMPS steps per LB step when dt_LB > dt_critical for BPM
# =========================================================================
def lbm_analysis(U_star=0.08):
"""Print LBM parameter selection, stability analysis, and LAMMPS values."""
cs = 1.0 / np.sqrt(3)
cs2 = 1.0 / 3.0
nx = int(x_length / dx_LB)
ny = int(y_length / dx_LB)
nz = int(z_length / dx_LB)
# Derived LBM parameters
dt_LB = U_star * dx_LB / u_phys
nu_LB = nuf * dt_LB / dx_LB**2
tau = 0.5 + nu_LB / cs2
Ma = U_star / cs
l_star = L / dx_LB
Re_phys = u_phys * L / nuf
Re_LB = U_star * l_star / nu_LB
vel_scale = dx_LB / dt_LB
force_scale = rho_f * dx_LB**4 / dt_LB**2
print(f"Re = {Re_phys:.2f}")
print(f"\nLBM grid: {nx} x {ny} x {nz} (dx_LB = {dx_LB})")
print(f"BPM grid: Nb = {Nb}, Na = {Na}, l0 = {l0:.6f}, Gr = {l0/dx_LB:.4f}")
print(f"\n{'=' * 70}")
print("LATTICE BOLTZMANN PARAMETERS")
print(f"{'=' * 70}")
print(f"\n Chosen: U* = {U_star}")
print(f"\n Derived:")
print(f" dt_LB = {dt_LB:.6f}")
print(f" nu* = {nu_LB:.6f}")
print(f" tau* = {tau:.6f}")
print(f" Ma = {Ma:.6f}")
print(f" l* = {l_star:.1f}")
print(f"\n Re check: Re_phys = {Re_phys:.4f} Re_LB = {Re_LB:.4f} "
f"{'MATCH' if abs(Re_phys - Re_LB)/max(Re_phys, 1e-15) < 1e-6 else 'MISMATCH'}")
print(f"\n Conversion scales:")
print(f" vel_scale = {vel_scale:.6f}")
print(f" force_scale = {force_scale:.6e}")
# Body force amplitude (oscillatory Stokes flow)
#
# For g(t) = g0 sin(wt) between walls separated by H, the centerline
# velocity is u_max ~ g0/w in the thick-channel limit (H >> delta_s).
# Stokes layer: delta_s = sqrt(2*nu/w). Here delta_s ~ 5.6 mm and
# H/2 = 30 mm, so the approximation is valid (error < 1%).
omega = 2 * np.pi * freq
delta_s = np.sqrt(2 * nuf / omega)
G0 = u_phys * omega
print(f"\n{'=' * 70}")
print("BODY FORCE AMPLITUDE (oscillatory Stokes flow)")
print(f"{'=' * 70}")
print(f"\n omega = 2*pi*f = {omega:.6e}")
print(f" delta_s = sqrt(2*nu/omega) = {delta_s:.4f} mm")
print(f" H/(2*delta_s) = {H/(2*delta_s):.2f} (>> 1: thick-channel limit)")
print(f"\n G0 = u_max * omega = {G0:.6e} (mm/ms^2)")
# Stability analysis
print(f"\n{'=' * 70}")
print("STABILITY ANALYSIS")
print(f"{'=' * 70}")
all_ok = True
print(f"\n 1. RELAXATION TIME: tau* = {tau:.6f}", end=" ")
if tau <= 0.5:
print("UNSTABLE"); all_ok = False
elif tau < 0.55:
print("MARGINAL")
elif tau > 2.0:
print("HIGH (over-damped)")
else:
print("GOOD")
print(f" 2. MACH NUMBER: Ma = {Ma:.6f}", end=" ")
if Ma > 0.3:
print("HIGH"); all_ok = False
elif Ma > 0.1:
print("ACCEPTABLE")
else:
print("EXCELLENT")
print(f" 3. GRID RESOLUTION: l* = {l_star:.1f}", end=" ")
if l_star < 10:
print("COARSE")
elif l_star < 20:
print("MODERATE")
else:
print("GOOD")
# Sensitivity: U*
print(f"\n 4. SENSITIVITY (varying U* at dx_LB = {dx_LB}):")
print(f" {'U*':<8} {'dt_LB':<12} {'nu*':<10} {'tau*':<10} {'Ma':<10} {'Status'}")
print(f" {'-'*66}")
for U_t in [0.02, 0.05, 0.08, 0.10, 0.15, 0.20]:
dt_t = U_t * dx_LB / u_phys
nu_t = nuf * dt_t / dx_LB**2
tau_t = 0.5 + 3 * nu_t
Ma_t = U_t / cs
st = ("UNSTABLE" if tau_t <= 0.5 else "marginal tau" if tau_t < 0.55
else "high Ma" if Ma_t > 0.3 else "over-damped" if tau_t > 2.0 else "OK")
mk = " <--" if abs(U_t - U_star) < 1e-10 else ""
print(f" {U_t:<8.3f} {dt_t:<12.6f} {nu_t:<10.6f} {tau_t:<10.6f} {Ma_t:<10.4f} {st}{mk}")
# Sensitivity: dx_LB
print(f"\n 5. SENSITIVITY (varying dx_LB at U* = {U_star}):")
print(f" {'dx_LB':<10} {'dt_LB':<12} {'nu*':<10} {'tau*':<10} {'l*':<8} {'Nx':<6} {'Status'}")
print(f" {'-'*72}")
for dx_t in [0.3125, 0.625, 1.0, 1.25, 2.0]:
dt_t = U_star * dx_t / u_phys
nu_t = nuf * dt_t / dx_t**2
tau_t = 0.5 + 3 * nu_t
l_t = L / dx_t
nx_t = int(x_length / dx_t)
st = ("UNSTABLE" if tau_t <= 0.5 else "marginal tau" if tau_t < 0.55
else "over-damped" if tau_t > 2.0 else "too coarse" if l_t < 10 else "OK")
mk = " <--" if abs(dx_t - dx_LB) < 1e-10 else ""
print(f" {dx_t:<10.4f} {dt_t:<12.6f} {nu_t:<10.6f} {tau_t:<10.6f} {l_t:<8.1f} {nx_t:<6} {st}{mk}")
# Timestep synchronization
Kr = E * A / l0
Ks = 12 * E * I / l0**3
Kt = G * J / l0
Kb = E * I / l0
m_seg = rho * A * l0
I_seg = m_seg * l0**2
omega_r = np.sqrt(Kr / m_seg)
omega_s = np.sqrt(Ks / m_seg)
omega_t = np.sqrt(Kt / I_seg)
omega_b = np.sqrt(Kb / I_seg)
omega_max = max(omega_r, omega_s, omega_t, omega_b)
f_max = omega_max / (2 * np.pi)
dt_crit = np.pi / omega_max
dt_safe = 0.1 / f_max
dt_lammps_target = min(0.5 * dt_safe, 0.5 * dt_crit)
lammps_per_lb = max(1, int(np.ceil(dt_LB / dt_lammps_target)))
dt_lammps = dt_LB / lammps_per_lb
print(f"\n{'=' * 70}")
print("TIMESTEP SYNCHRONIZATION (CoupLB md_per_lb subcycling)")
print(f"{'=' * 70}")
print(f"\n Spring frequencies:")
print(f" omega_r = {omega_r:.4f} omega_s = {omega_s:.4f} "
f"omega_t = {omega_t:.4f} omega_b = {omega_b:.4f}")
print(f" omega_max = {omega_max:.4f} f_max = {f_max:.4f}")
print(f"\n Stability limits:")
print(f" dt_critical = {dt_crit:.6f} dt_safe = {dt_safe:.6f}")
print(f"\n Timesteps:")
print(f" dt_LB = {dt_LB:.6f}")
print(f" dt_LAMMPS = {dt_lammps:.6f}")
print(f" md_per_lb = {lammps_per_lb}")
r_crit = dt_lammps / dt_crit
r_safe = dt_lammps / dt_safe
print(f"\n Safety: dt/dt_crit = {r_crit:.4f} dt/dt_safe = {r_safe:.4f}")
# Simulation time
T_osc = 1.0 / freq # period = 1000 ms = 1 s
n_cycles = 10
sim_time = n_cycles * T_osc # 10 s
n_steps = int(np.ceil(sim_time / dt_lammps))
n_lb = n_steps // lammps_per_lb
print(f"\n{'=' * 70}")
print("SIMULATION TIME")
print(f"{'=' * 70}")
print(f" Target: {sim_time:.0f} ms ({sim_time/1000:.0f} s, {n_cycles} cycles at {1/T_osc*1000:.0f} Hz)")
print(f" LAMMPS: {n_steps:,} steps @ dt = {dt_lammps:.6f}")
print(f" LB: {n_lb:,} steps @ dt = {dt_LB:.6f}")
print(f" md_per_lb = {lammps_per_lb}")
# Summary
print(f"\n{'=' * 70}")
print("VALUES FOR flow.in")
print(f"{'=' * 70}")
print(f" timestep {dt_lammps:.6f}")
print(f" nuf {nuf:.6f}")
print(f" rhof {rho_f:.4f}")
print(f" G0 {G0:.6e}")
print(f" freq {freq:.6e}")
print(f" grid {nx} {ny} {nz}")
print(f" md_per_lb {lammps_per_lb}")
print(f" run {n_steps}")
if all_ok:
print(f"\nALL CHECKS PASSED")
else:
print(f"\nISSUES DETECTED -- REVIEW ABOVE")
if __name__ == '__main__':
lbm_analysis()
LAMMPS input script
The input script is configured for 45 MPI ranks (processors 5 3 3).
flow.in ⬇ Download
# =========================================================================
# Oscillatory channel flow with elastic filaments
# =========================================================================
# Filament parameters (milli units: mm, ms, mg)
variable Nb equal 16
variable L equal 10
variable d equal 1
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 1230
variable G equal 473
# =========================================================================
# Initialization
# =========================================================================
units lj
dimension 3
processors 5 3 3
atom_style hybrid molecular bpm/sphere
special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
newton on off
boundary p f p
neighbor 0.5 bin
neigh_modify delay 0 every 1
comm_modify mode single cutoff $(2*v_l0) vel yes
read_data filaments.lam
# =========================================================================
# Groups
# =========================================================================
group bottom type 1
group interior type 2
group tip type 3
group free union interior tip
# Track the tip of the rightmost filament
variable trackid equal 5*(v_Nb+1)
group track id ${trackid}
variable xt equal x[${trackid}]
variable yt equal y[${trackid}]
variable zt equal z[${trackid}]
# =========================================================================
# 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 * * $(2*v_l0)
# =========================================================================
# LBM parameters
# =========================================================================
variable nuf equal 0.1
variable rhof equal 1.26
# Oscillating body force: g(t) = g0 sin(2 pi f t)
variable freq equal 1e-3
variable g0 equal 3.8e-3
variable time equal step*v_timestep
variable cyc equal v_freq*v_time
variable gx_osc equal v_g0*sin(2*PI*v_cyc)
# =========================================================================
# Dynamics
# =========================================================================
variable timestep equal 0.003968
timestep ${timestep}
compute mytemp all temp/sphere
compute kerot all erotate/sphere
# Initial tip position (reference for displacement)
variable xt0 equal ${xt}
variable yt0 equal ${yt}
variable zt0 equal ${zt}
variable dxt equal v_xt-v_xt0
variable dyt equal v_yt-v_yt0
variable dzt equal v_zt-v_zt0
thermo 2100
thermo_style custom step v_cyc temp ke c_kerot pe v_dxt v_dyt v_dzt
thermo_modify temp mytemp norm no
fix prnt all print 2100 &
"${cyc} ${dxt} ${dyt} ${dzt}" &
file right-tip.txt screen no &
title "cycle dx dy dz"
fix nve free nve/bpm/sphere
fix ibm free couplb 160 96 80 ${nuf} ${rhof} &
md_per_lb 21 &
wall_y 1 1 &
gravity v_gx_osc 0 0 &
output 6300 fluid.dat &
check_every 25200 &
vtk 6300 ./vtk/flow &
vtk_solid id type mass vx vy vz fx fy fz
run 2520000
References
- Pinelli, A., Omidyeganeh, M., Brücker, C., Revell, A., Sarkar, A., Alinovi, E. (2017). The PELskin project: part IV — control of bluff body wakes using hairy filaments. Meccanica, 52, 1503–1514. doi:10.1007/s11012-016-0513-0
- Agrawal, V., Kulachenko, A., Scapin, N., Tammisola, O., Brandt, L. (2024). An efficient isogeometric/finite-difference immersed boundary method for the fluid–structure interactions of slender flexible structures. Comput. Methods Appl. Mech. Engrg., 418, 116495. doi:10.1016/j.cma.2023.116495