Magnetic Beams
Based on Yan et al. (2022)
These scripts simulate the three configurations in Fig. 6(b) of our paper (cf. Fig. 9 in Yan et al.).
For Fig. 6(a) (cf. Fig. 8(a)), add gravity along \(+x\) to Configuration B1 using fix gravity (see the twisting rod example) and center the field at \((L/2, 0, 0)\) by replacing \(x\) with \(x-L/2\) in the potential.
Parameters
| Parameter | Value |
|---|---|
| Beam geometry (\(L \times W \times H\)) | 25.8 mm \(\times\) 1.21 mm \(\times\) 0.49 mm |
| Young's modulus \(E\) | 1.16 MPa |
| Poisson's ratio \(\nu\) | 0.49 |
| Density \(\rho\) | 2.01 g/cm\(^3\) |
| Magnetization magnitude \(M\) | 94.1 kA/m |
| Number of bonds \(N\) | 258 (bond length \(\ell_0 = 0.1\) mm) |
Setup
All three cases use the same straight rectangular beam. The leftmost atom is clamped, motion is restricted to 2D, and the beam is equilibrated quasi-statically with Langevin dynamics after each magnetic loading increment. The quantity of interest is the normalized tip deflection \(\delta_y/L\), measured at the rightmost atom as in the cantilever beam demo.
The field strength per unit \(\lambda_m^C\) is \(B_1 = EI/(MAL^2)\). For gradient-field configs (B1, B2), the gradient is scaled by \(\lambda_m^\nabla \approx 0.36\,\lambda_m^C\) so that all three scripts report loading on the same \(\lambda_m^C\) scale.
Configuration A (config-A.in)
The beam is magnetized in the \(+x\) direction and a uniform magnetic field is applied in the \(+y\) direction:
\[\mathbf{B} = B \hat{\mathbf{y}}.\]
The simulation ramps the dimensionless parameter \(\lambda_m^C = MBAL^2/(EI)\).
Here fix efield is used by analogy with magnetic dipoles;
it applies the dipole torque \(\boldsymbol{\tau} = \mathbf{m} \times \mathbf{B}\).
Configuration B1 (config-B1.in)
The beam is magnetized in the \(+y\) direction and a constant-gradient magnetic field is applied. The magnetic field is defined using the scalar potential
\[\psi = -\frac{b}{\mu_0}\left[\frac{y^2}{2} - \frac{x^2 + z^2}{4}\right].\]
This is applied with fix efield/lepton.
Unlike fix efield, which only applies the dipole torque, fix efield/lepton also captures the force from the field gradient through the scalar potential.
This yields \(\mathbf{B} = -\mu_0 \nabla \psi = b \langle -x/2, y, -z/2 \rangle\) and \(\nabla \mathbf{B} = b\,\mathrm{diag}(-1/2, 1, -1/2)\), where \(\lambda_m^\nabla \approx 0.36 \lambda_m^C\).
Note
fix efield/lepton is used through the electric-field analogy and expects an electric potential expression.
The equivalent potential passed to LAMMPS is therefore written to obtain the desired gradient prefactor dB; it is not the actual magnetic scalar potential, so no explicit \(\mu_0\) appears in the script.
Since these simulations are 2D with \(z=0\), the input scripts only include the \(x\) and \(y\) terms.
Configuration B2 (config-B2.in)
The beam is magnetized in the \(+x\) direction with the same field as B1, with the addition of a small uniform \(+y\) perturbation to kick the beam out of the initially straight equilibrium.
Key LAMMPS settings
- Units:
units ljusing our milli unit system - Atom style:
atom_style hybrid bpm/sphere dipole - Pair style:
pair_style zero; no contact or dipole-dipole interactions - Magnetic loading:
fix efieldfor the uniform andfix efield/leptonfor the non-uniform fields - Time integration:
fix nve/bpm/sphere update dipole - Damping:
fix langevinfor quasi-static equilibration
Files
Data file generator
The generator writes magnetic-beam-N258.lam, which is read by all three input scripts.
Note
The data file generator writes the dipole as zero. Each configuration script sets the dipole direction: \(+x\) for A and B2, \(+y\) for B1.
create_magnetic_beam.py ⬇ Download
import numpy as np
N = 258 # number of bonds
# Physical properties (milli units: mm, ms, mg)
L = 25.8
W = 1.21
H = 0.49
density = 2.01 # mg/mm^3
# Discretization
l0 = L / N
diam = l0 # arbitrary
mass = density * l0 * W * H
dens = mass / (np.pi / 6 * diam**3)
# atom_style hybrid bpm/sphere dipole
with open('magnetic-beam-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' % (-1, 27))
file.write('%f %f ylo yhi \n' % (-27, 27))
file.write('%f %f zlo zhi \n\n' % (-0.1, 0.1))
file.write('Atoms \n\n')
x = 0.
for k in range(N+1):
# atom-id type x y z molid diameter density q mux muy muz
file.write('%i 1 %f 0 0 1 %f %f 0 0 0 0 \n' % (k+1, x, diam, dens))
x += l0
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 scripts
config-A.in ⬇ Download
# =========================================================================
# Magnetized beam — Config A
# Beam magnetized in +x, uniform field in +y
# =========================================================================
variable N equal 258
# Beam geometry
variable L equal 25.8
variable W equal 1.21
variable H equal 0.49
variable l0 equal $(v_L/v_N)
# Rectangular cross-section
variable A equal v_W*v_H
variable I equal v_W*v_H^3/12
variable J equal v_W*v_H/12*(v_W^2+v_H^2)
# Material properties
variable E equal 1.16e3
variable nu equal 0.49
variable G equal "v_E/(2*(1+v_nu))"
# Magnetization
variable M equal 0.941
variable mu equal v_M*v_l0*v_W*v_H
# =========================================================================
# Initialization
# =========================================================================
units lj # milli units
dimension 2
boundary f f p
processors * 1 1
atom_style hybrid bpm/sphere dipole
special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
newton on off
neighbor $(v_l0) bin
neigh_modify delay 0 every 10000
atom_modify sort 10000 1
comm_modify mode single cutoff $(2*v_l0) vel yes
read_data magnetic-beam-N${N}.lam
# Set dipole in +x direction
set group all dipole ${mu} 0 0
# =========================================================================
# Groups
# =========================================================================
region left block INF $(v_l0/2) INF INF INF INF
region right block $(v_L-v_l0/2) INF INF INF INF INF
group left region left
group right region right
group free subtract all left
# =========================================================================
# 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.01 0.01 0.01 0.01
pair_style zero $(2*v_l0)
pair_coeff * * $(2*v_l0)
# =========================================================================
# Computes and output
# =========================================================================
timestep 5e-4
compute delx right reduce ave x
compute dely right reduce ave y
variable dx equal "c_delx/v_L-1"
variable dy equal "c_dely/v_L"
variable l equal 0
variable Bone equal v_E*v_I/(v_M*v_A*v_L^2)
thermo 10000
thermo_style custom step temp ke pe etotal v_l v_dx v_dy
thermo_modify lost ignore flush yes norm no
dump 1 all custom 100000 ./config-A.dump id type x y fx fy mux muy
fix prnt all print 100000 &
"${l} ${dx} ${dy}" &
file config-A.txt screen no
# =========================================================================
# Loop 1: small field increments (lambda_m^C = 0.5 to 5)
# =========================================================================
variable i index 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
label LOOP1
variable By equal "v_Bone*v_i"
variable l equal "v_i"
bond_coeff 1 ${Kr} ${Ks} ${Kt} ${Kb} &
-1 -1 -1 -1 &
0.01 0.01 0.01 0.01
fix 1 free nve/bpm/sphere update dipole
fix 2 free langevin 0 0 1000 83715463 zero yes omega yes
fix 3 all efield 0 v_By 0
fix 4 all enforce2d
fix_modify 3 energy yes
run 500000
velocity all zero linear
velocity all zero angular
bond_coeff 1 ${Kr} ${Ks} ${Kt} ${Kb} &
-1 -1 -1 -1 &
0.05 0.05 0.05 0.05
fix 2 free langevin 0 0 100 83715463 zero yes omega yes
run 2500000
next i
jump SELF LOOP1
# =========================================================================
# Loop 2: large field increments (lambda_m^C = 6 to 45)
# =========================================================================
variable i index 6 7 8 9 10 12.5 15 17.5 20 25 30 35 40 45
label LOOP2
variable By equal "v_Bone*v_i"
variable l equal "v_i"
bond_coeff 1 ${Kr} ${Ks} ${Kt} ${Kb} &
-1 -1 -1 -1 &
0.05 0.05 0.05 0.05
fix 1 free nve/bpm/sphere update dipole
fix 2 free langevin 0 0 100 83715463 zero yes omega yes
fix 3 all efield 0 v_By 0
fix 4 all enforce2d
fix_modify 3 energy yes
run 500000
next i
jump SELF LOOP2
config-B1.in ⬇ Download
# =========================================================================
# Magnetized beam — Config B1
# Beam magnetized in +y, constant-gradient field (dB_yy > 0)
# =========================================================================
variable N equal 258
# Beam geometry
variable L equal 25.8
variable W equal 1.21
variable H equal 0.49
variable l0 equal $(v_L/v_N)
# Rectangular cross-section
variable A equal v_W*v_H
variable I equal v_W*v_H^3/12
variable J equal v_W*v_H/12*(v_W^2+v_H^2)
# Material properties
variable E equal 1.16e3
variable nu equal 0.49
variable G equal "v_E/(2*(1+v_nu))"
# Magnetization
variable M equal 0.941
variable mu equal v_M*v_l0*v_W*v_H
# =========================================================================
# Initialization
# =========================================================================
units lj # milli units
dimension 2
boundary f f p
processors * 1 1
atom_style hybrid bpm/sphere dipole
special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
newton on off
neighbor $(v_l0) bin
neigh_modify delay 0 every 10000
atom_modify sort 10000 1
comm_modify mode single cutoff $(2*v_l0) vel yes
read_data magnetic-beam-N${N}.lam
# Set dipole in +y direction
set group all dipole 0 ${mu} 0
# =========================================================================
# Groups
# =========================================================================
region left block INF $(v_l0/2) INF INF INF INF
region right block $(v_L-v_l0/2) INF INF INF INF INF
group left region left
group right region right
group free subtract all left
# =========================================================================
# 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.01 0.01 0.01 0.01
pair_style zero $(2*v_l0)
pair_coeff * * $(2*v_l0)
# =========================================================================
# Computes and output
# =========================================================================
timestep 5e-4
compute delx right reduce ave x
compute dely right reduce ave y
variable dx equal "c_delx/v_L-1"
variable dy equal "c_dely/v_L"
variable l equal 0
variable dB equal 0
variable dBone equal v_E*v_I/(v_M*v_A*v_L^3)*0.36
thermo 10000
thermo_style custom step temp ke pe etotal v_l v_dx v_dy
thermo_modify lost ignore flush yes norm no
dump 1 all custom 100000 ./config-B1.dump id type x y fx fy mux muy
fix prnt all print 100000 &
"${l} ${dx} ${dy}" &
file config-B1.txt screen no
# =========================================================================
# Loop: constant-gradient field increments (lambda_m^C = 1 to 45)
# =========================================================================
variable i loop 45
label LOOP1
variable dB equal "v_dBone*v_i"
variable l equal "v_i"
bond_coeff 1 ${Kr} ${Ks} ${Kt} ${Kb} &
-1 -1 -1 -1 &
0.01 0.01 0.01 0.01
fix 1 free nve/bpm/sphere update dipole
fix 2 free langevin 0 0 10000 83715463 zero yes omega yes
fix 3 all efield/lepton "-v_dB*(-x^2+2*y^2)/4" step 1
fix 4 all enforce2d
fix_modify 3 energy yes
run 100000
velocity all zero linear
velocity all zero angular
bond_coeff 1 ${Kr} ${Ks} ${Kt} ${Kb} &
-1 -1 -1 -1 &
0.05 0.05 0.05 0.05
fix 2 free langevin 0 0 10 83715463 zero yes omega yes
run 900000
next i
jump SELF LOOP1
config-B2.in ⬇ Download
# =========================================================================
# Magnetized beam — Config B2
# Beam magnetized in +x, constant-gradient field with a small +y perturbation
# =========================================================================
variable N equal 258
# Beam geometry
variable L equal 25.8
variable W equal 1.21
variable H equal 0.49
variable l0 equal $(v_L/v_N)
# Rectangular cross-section
variable A equal v_W*v_H
variable I equal v_W*v_H^3/12
variable J equal v_W*v_H/12*(v_W^2+v_H^2)
# Material properties
variable E equal 1.16e3
variable nu equal 0.49
variable G equal "v_E/(2*(1+v_nu))"
# Magnetization
variable M equal 0.941
variable mu equal v_M*v_l0*v_W*v_H
# =========================================================================
# Initialization
# =========================================================================
units lj # milli units
dimension 2
boundary f f p
processors * 1 1
atom_style hybrid bpm/sphere dipole
special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0
newton on off
neighbor $(v_l0) bin
neigh_modify delay 0 every 10000
atom_modify sort 10000 1
comm_modify mode single cutoff $(2*v_l0) vel yes
read_data magnetic-beam-N${N}.lam
# Set dipole in +x direction
set group all dipole ${mu} 0 0
# =========================================================================
# Groups
# =========================================================================
region left block INF $(v_l0/2) INF INF INF INF
region right block $(v_L-v_l0/2) INF INF INF INF INF
group left region left
group right region right
group free subtract all left
# =========================================================================
# 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.01 0.01 0.01 0.01
pair_style zero $(2*v_l0)
pair_coeff * * $(2*v_l0)
# =========================================================================
# Computes and output
# =========================================================================
timestep 5e-4
compute delx right reduce ave x
compute dely right reduce ave y
variable dx equal "c_delx/v_L-1"
variable dy equal "c_dely/v_L"
variable l equal 0
variable dB equal 0
variable dBone equal v_E*v_I/(v_M*v_A*v_L^3)*0.36
thermo 10000
thermo_style custom step temp ke pe etotal v_l v_dx v_dy
thermo_modify lost ignore flush yes norm no
dump 1 all custom 100000 ./config-B2.dump id type x y fx fy mux muy
fix prnt all print 100000 &
"${l} ${dx} ${dy}" &
file config-B2.txt screen no
# =========================================================================
# Loop 1: small field increments (lambda_m^C = 1 to 4)
# =========================================================================
variable i loop 4
label LOOP1
variable dB equal "v_dBone*v_i"
variable l equal "v_i"
fix 1 free nve/bpm/sphere update dipole
fix 2 free langevin 0 0 1000 83715463 zero yes omega yes
fix 3 all efield/lepton "-v_dB*(-x^2+2*y^2)/4 -0.001*y" step 1
fix 4 all enforce2d
fix_modify 3 energy yes
run 1000000
next i
jump SELF LOOP1
# =========================================================================
# Loop 2: transition increments (lambda_m^C = 5.0 to 5.5)
# =========================================================================
variable i index 5.0 5.1 5.2 5.3 5.4 5.5
label LOOP2
variable dB equal "v_dBone*v_i"
variable l equal "v_i"
bond_coeff 1 ${Kr} ${Ks} ${Kt} ${Kb} &
-1 -1 -1 -1 &
0.01 0.01 0.01 0.01
fix 1 free nve/bpm/sphere update dipole
fix 2 free langevin 0 0 10000 83715463 zero yes omega yes
fix 3 all efield/lepton "-v_dB*(-x^2+2*y^2)/4 -0.001*y" step 1
fix 4 all enforce2d
fix_modify 3 energy yes
run 500000
velocity all zero linear
velocity all zero angular
bond_coeff 1 ${Kr} ${Ks} ${Kt} ${Kb} &
-1 -1 -1 -1 &
0.05 0.05 0.05 0.05
fix 2 free langevin 0 0 1000 83715463 zero yes omega yes
run 1500000
next i
jump SELF LOOP2
# =========================================================================
# Loop 3: large field increments (lambda_m^C = 6 to 45)
# =========================================================================
variable i index 6 6.5 7 7.5 8 8.5 9 9.5 10 12.5 15 20 25 30 35 40 45
label LOOP3
variable dB equal "v_dBone*v_i"
variable l equal "v_i"
bond_coeff 1 ${Kr} ${Ks} ${Kt} ${Kb} &
-1 -1 -1 -1 &
0.05 0.05 0.05 0.05
fix 1 free nve/bpm/sphere update dipole
fix 2 free langevin 0 0 100 83715463 zero yes omega yes
fix 3 all efield/lepton "-v_dB*(-x^2+2*y^2)/4 -0.001*y" step 1
fix 4 all enforce2d
fix_modify 3 energy yes
run 1000000
next i
jump SELF LOOP3
References
- Yan, D., Abbasi, A., Reis, P.M. (2022). A comprehensive framework for hard-magnetic beams: reduced-order theory, 3D simulations, and experiments. Int. J. Solids Struct., 257, 111319. doi:10.1016/j.ijsolstr.2021.111319