Magnetic Helical Rods
Based on Sano (2022)
A helical rod magnetized along its screw axis contracts under dipole–dipole interactions. An axial-gradient magnetic field stretches the helix, and quasi-static loading–unloading cycles produce mechanical hysteresis.
Parameters
| Parameter | Value |
|---|---|
| Centerline length \(L\) | 103 mm |
| Rod diameter \(d\) | 2 mm |
| Young's modulus \(E\) | 1.16 MPa |
| Poisson's ratio \(\nu\) | 0.5 |
| Density \(\rho\) | 1.23 g/cm\(^3\) |
| Number of turns | 1.64 |
| Normalized pitch \(\tilde{\psi}\) | 0.96 |
| Number of bonds \(N\) | 103 (bond length \(\ell_0 = 1\) mm) |
| Dipole strength \(\gamma\) | varies |
| Field strength \(\lambda_m\) | varies |
Setup
The helix is initialized in its stress-free configuration and the topmost atom at \(z \approx 0\) is clamped. The centerline is parametrized by arclength \(s\):
\[\mathbf{r}(s) = \langle R\cos(Ks),\; R\sin(Ks),\; s\cos\psi \rangle,\]
where \(K = R^{-1}\sin\psi\) and \(\tilde{\psi} = \psi/(\pi/2) = 0.96\).
All dipoles are initialized along the \(+z\) screw axis. Dipole–dipole interactions are controlled by \(\gamma = \mu_0 M^2 / E\). Since the problem is governed by dimensionless parameters, the value of \(E\) in the input script can be decreased to speed up the simulation. No gravity is applied because the experiments were performed in a glycerol bath.
The external gradient field is defined by the magnetic scalar potential
\[\psi = -\frac{b}{\mu_0}\left[\frac{x^2 + y^2}{4} - \frac{z^2}{2}\right],\]
which gives \(\mathbf{B} = -\mu_0\nabla\psi = b\langle x/2,\; y/2,\; -z\rangle\). The prefactor \(b\) is related to the dimensionless parameters by
\[b = \sqrt{\frac{\mu_0}{\gamma E}}\,\frac{EI}{A}\left(\frac{\lambda_m}{L}\right)^3,\]
where \(\lambda_m = L\,[MbA/(EI)]^{1/3}\) measures the external field strength relative to elastic stiffness.
The field strength is ramped up from \(\lambda_m = 0\) to \(3.5\) in increments of \(0.1\), then ramped back down.
Key LAMMPS settings
- Units:
units cgs - Atom style:
atom_style hybrid bpm/sphere dipole - Contact:
pair_style gran/hertz/history;special_bonds lj 0 1 1 coul 0 1 1disables contact and dipole–dipole forces for bonded (1–2) neighbors - Dipole–dipole interactions:
pair_style lj/cut/dipole/cutwith the LJ energy set to zero, combined with contact throughpair_style hybrid/overlay - Magnetic loading:
fix efield/lepton - Time integration:
fix nve/bpm/sphere update dipole - Damping:
fix langevinfor quasi-static equilibration
Files
Data file generator
The generator writes helix-N103.lam.
It stores zero dipoles in the data file; the LAMMPS input script sets the dipole magnitude from the selected \(\gamma\).
create_helical_rod.py ⬇ Download
import numpy as np
N = 103 # number of bonds
# Units: CGS EMU
# Physical parameters
L = 10.3
d = 0.2
rho = 1.23
# Helix geometry
R = 1.0 # coil radius
psi_tilde = 0.96 # normalized pitch angle
psi = psi_tilde * np.pi / 2
K = np.sin(psi) / R
nturns = K * L / (2 * np.pi)
cospsi = np.cos(psi)
# Discretization
ds = L / N
A = np.pi * d**2 / 4
mass = rho * ds * A
diam = d
m = mass * np.ones(N+1)
dens = m / (np.pi * diam**3 / 6)
# Helix positions: z goes from -H0 to 0, clamped atom at z = 0 (top)
s = np.linspace(0, L, N+1)
x = R * np.cos(-K * (s - L))
y = R * np.sin(-K * (s - L))
z = cospsi * (s - L)
with open('helix-N{}.lam'.format(int(N)), 'w') as file:
file.write('LAMMPS DATA FILE \n\n')
file.write(str(N+1) + ' atoms \n')
file.write(str(N) + ' bonds \n\n')
file.write('1 atom types \n')
file.write('1 bond types \n\n')
file.write('%f %f xlo xhi \n' % (-2*R, 2*R))
file.write('%f %f ylo yhi \n' % (-2*R, 2*R))
file.write('%f %f zlo zhi \n\n' % (z[0]-d, z[-1]+d))
# atom_style hybrid bpm/sphere dipole
file.write('Atoms \n\n')
for k in range(N+1):
# atom-id type x y z molid diameter density q mux muy muz
file.write('%i 1 %f %f %f 1 %f %f 0 0 0 0\n' %
(k+1, x[k], y[k], z[k], diam, dens[k]))
file.write('\n')
file.write('Bonds \n\n')
for k in range(N):
file.write('%i %i %i %i \n' % (k+1, 1, k+1, k+2))
file.write('\n')
LAMMPS input script
deform-helix.in ⬇ Download
# =========================================================================
# Magnetic loading of a helical rod
# Helix magnetized along +z, axial gradient field
#
# Replace GAMMA with the desired value before running.
# =========================================================================
variable gamma equal GAMMA
variable gammastr string GAMMA
# Helix geometry
variable L equal 10.3
variable N equal 103
variable l0 equal $(v_L/v_N)
variable d equal 0.2
# 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
variable E equal 1.16e7
variable nu equal 0.5
variable G equal "v_E/(2*(1+v_nu))"
# Field prefactor: b = bprefactor * lambda_m^3
variable mu0 equal 4*PI
variable bprefactor equal $(sqrt(v_mu0/(v_gamma*v_E))*v_E*v_I/(v_A*v_L^3))
# =========================================================================
# Initialization
# =========================================================================
units cgs
dimension 3
boundary f f f
processors * * 1
atom_style hybrid bpm/sphere dipole
special_bonds lj 0.0 1.0 1.0 coul 0.0 1.0 1.0
newton on off
neighbor 0.1 bin
neigh_modify delay 0 every 10000
atom_modify sort 10000 1
comm_modify mode single cutoff 1.0 vel yes
read_data helix-N${N}.lam
change_box all x scale 3 y scale 3 z final $(-0.8*v_L) $(0.1*v_L)
# =========================================================================
# Groups
# =========================================================================
group fixed id $(v_N+1)
group free subtract all fixed
# =========================================================================
# 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.1 0.1 0.1 0.1
# =========================================================================
# Dipole-dipole interactions and contact
# =========================================================================
variable mu equal sqrt(${gamma}*v_E/v_mu0)*v_A*v_l0
set group all dipole 0 0 ${mu}
pair_style hybrid/overlay gran/hertz/history $(10*v_E) NULL 0.5 NULL 0.1 1 &
lj/cut/dipole/cut 0.0 1.0
pair_coeff * * gran/hertz/history
pair_coeff * * lj/cut/dipole/cut 0.0 0.0
# =========================================================================
# Computes and output
# =========================================================================
timestep 1e-7
compute zmax all reduce max z
compute zmin all reduce min z
variable dz equal "(c_zmax-c_zmin)/v_L"
variable l equal 0
variable b equal v_bprefactor*v_l^3
thermo 10000
thermo_style custom step temp ke pe etotal v_l v_dz
thermo_modify lost ignore flush yes norm no
dump 1 all custom 100000 ./deform-helix-g${gammastr}.dump &
id type x y z mux muy muz
fix prnt all print 100000 &
"${l} ${dz}" &
file deform-helix-g${gammastr}.txt screen no &
title "lambda_m z/L"
run 0
# =========================================================================
# Loop 1: loading (lambda_m = 0 to 3.5)
# =========================================================================
variable i loop 36
label LOOP_UP
variable l equal "(v_i-1)/10"
fix 1 free nve/bpm/sphere update dipole
fix 2 free langevin 0 0 1 83715463 zero yes omega yes
fix 3 all efield/lepton "-v_b*((x^2+y^2)/4-z^2/2);" step 1
fix_modify 3 energy yes
run 10000000
velocity all zero linear
velocity all zero angular
fix 2 free langevin 0 0 1e-1 83715463 zero yes omega yes
run 10000000
next i
jump SELF LOOP_UP
# =========================================================================
# Loop 2: unloading (lambda_m = 3.4 to 0)
# =========================================================================
variable i loop 35
label LOOP_DOWN
variable l equal "(35-v_i)/10"
fix 1 free nve/bpm/sphere update dipole
fix 2 free langevin 0 0 1 83715463 zero yes omega yes
fix 3 all efield/lepton "-v_b*((x^2+y^2)/4-z^2/2);" step 1
fix_modify 3 energy yes
run 10000000
velocity all zero linear
velocity all zero angular
fix 2 free langevin 0 0 1e-1 83715463 zero yes omega yes
run 10000000
next i
jump SELF LOOP_DOWN
References
- Sano, T.G. (2022). Reduced theory for hard magnetic rods with dipole–dipole interactions. J. Phys. A: Math. Theor., 55(10), 104002. doi:10.1088/1751-8121/ac4de2
- Sano, T.G., Pezzulla, M., Reis, P.M. (2022). A Kirchhoff-like theory for hard magnetic rods under geometrically nonlinear deformation in three dimensions. J. Mech. Phys. Solids, 160, 104739. doi:10.1016/j.jmps.2021.104739