Demo: Cantilever Beams

Note

This demo assumes familiarity with running LAMMPS simulations (input scripts, data files, fixes, and computes). It is aimed at users experienced with molecular dynamics who want to learn the BPM workflow for simulating elastic rods.

This demo walks through the dynamic cantilever beam simulation based on the example from Le et al. 2011. By the end, you'll learn:

How to generate a data file with Python
How to define the BPM parameters
How to set up a simple simulation

Additionally, the scripts for the static loading examples are provided at the end.

Cantilever beam under dynamic loading

A rectangular beam with one end clamped has its free end subject to a sinusoidal transverse force \(P(t) = P_0 \sin(\omega t)\), where \(P_0 = 10\) MN and \(\omega = 50\) rad/s. The following table lists the parameters.

Parameter	Value
Beam geometry (\(L \times w \times h\))	\(10\) m \(\times\) \(0.5\) m \(\times\) \(0.25\) m
Young's modulus \(E\)	210 GPa
Poisson's ratio \(\nu\)	0.3
Density \(\rho\)	7850 kg/m\(^3\)
Tip force amplitude \(P_0\)	10 MN
Angular frequency \(\omega\)	50 rad/s

Part 1: Data File Generation

Although it is possible to create the data file within LAMMPS, especially for simple 1D systems, we prefer to do this using a Python script for flexibility.

For a straight beam, atoms are placed at evenly spaced positions along the x-axis with sequential IDs. Bonds connect consecutive atom pairs (1–2, 2–3, ..., N–N+1), giving \(N\) bonds and \(N+1\) atoms. A Python script is provided below.

create_beam.py ⬇ Download

import numpy as np

N = 64  # number of bonds

# Units: SI

# Physical parameters
L = 10.
w = 0.5
h = 0.25
rho = 7850.

# Discretization
x = np.linspace(0, L, N+1)
dx = L/N

diam = h                    # set atom diameter to h; not required
mass = rho * dx * w * h     # mass of each bond...
m = mass * np.ones(N+1)    
m[0] = m[-1] = 0.5 * mass   # ...evenly split between the atoms
dens = m / (np.pi * diam**3 / 6)

with open('beam-{}.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')

    # Create box
    file.write('%f %f xlo xhi \n'%(-1, 11))
    file.write('%f %f ylo yhi \n'%(-11, 11))
    file.write('%f %f zlo zhi \n\n'%(-0.1, 0.1))

    # atom_style bpm/sphere
    file.write('Atoms \n\n')
    for k in range(N+1):
        # atom-id atom-type molid diameter density x y z
        file.write('%i 1 1 %f %f %f 0 0 \n'%(k+1, diam, dens[k], x[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')

Important details:

atom_style bpm/sphere uses per-atom mass set via density and not per-type mass. The effective density is back-calculated from the target mass and the chosen diameter.
Mass: each bond with rest length \(\ell_0 = L/N\) represents a rectangular segment of volume \(\ell_0 \times w \times h\). The masses are split between the two atoms forming the bond.
Diameter: the atom diameter is a geometric parameter for LAMMPS, not the physical rod diameter. Here we set it to \(h = 0.25\) m.

Part 2: BPM Parameters

Note

Although the BPM assumes cylindrical bond segments, the four spring constants can be defined independently. For planar bending without twisting, \(K_r\), \(K_s\), and \(K_b\) can be consistently defined for any cross-section. The twist stiffness \(K_t\) does not affect the result in 2D.

For a rectangular beam with width \(w\) and depth \(h\):

Area: \(A = wh\)
Second moment of area: \(I = wh^3/12\)
Polar second moment: \(J = wh(w^2 + h^2)/12\)

The four spring constants are as follows:

Stiffness	Formula	Role
\(K_r\)	\(EA / \ell_0\)	Axial stretching
\(K_s\)	\(12EI / \ell_0^3\)	Transverse shear
\(K_t\)	\(GJ / \ell_0\)	Twisting
\(K_b\)	\(EI / \ell_0\)	Bending

The \(K_s\) and \(K_b\) formulas follow from Euler–Bernoulli beam theory (Chen et al. 2022). All damping coefficients are set to zero since this is a fully dynamic, undamped simulation to test the natural frequency response.

Atom diameter and moment of inertia

The atom diameter is set to \(h\) in the data file. LAMMPS computes the moment of inertia of each bpm/sphere atom as \(I_\text{sphere} = \tfrac{1}{10} m d^2\). With \(d = h\), this is close enough to the Timoshenko rotary inertia \(\tfrac{1}{12} m h^2\).

An alternative is to adjust the density and the diameter of the atom to produce the desired mass and moment of inertia.

Note

The fix nve/bpm/sphere also has a disc keyword, which only changes the default moment of inertia formula. Since the diameter and density in the data file assume spherical atoms, do not add the disc keyword.

Part 3: Simulation Setup

The left and right groups are defined using regions that only capture the leftmost and rightmost atoms, respectively (assuming \(N < 100\)). The clamped boundary condition on the left end is imposed by excluding left from time integration. This means the leftmost atom's position and orientation are never updated. The sinusoidal tip force is applied to the rightmost atom via fix addforce. The total simulation time is \(1.5\) s.

Timestep selection

The timestep \(\Delta t\) for explicit time integration must be smaller than the critical timestep set by the highest-frequency mode. For translational modes, \(\omega = \sqrt{k/m}\); for rotational modes, \(\omega = \sqrt{k/I_\text{sphere}}\).

Assuming \(N = 100\) (\(\ell_0 = 0.1\) m, \(m \approx 98\) kg, \(I \approx 0.61\) kg m\(^2\)):

Mode	Stiffness	\(\omega\) [rad/s]	\(\Delta t_c = 2/\omega\) [s]
Radial	\(K_r = EA/\ell_0 \approx 2.6 \times 10^{11}\)	\(5.2 \times 10^{4}\)	\(3.9 \times 10^{-5}\)
Shear	\(K_s = 12EI/\ell_0^3 \approx 1.6 \times 10^{12}\)	\(1.3 \times 10^{5}\)	\(1.5 \times 10^{-5}\)
Twist	\(K_t = GJ/\ell_0 \approx 2.6 \times 10^{9}\)	\(6.5 \times 10^{4}\)	\(3.1 \times 10^{-5}\)
Bend	\(K_b = EI/\ell_0 \approx 1.4 \times 10^{9}\)	\(4.7 \times 10^{4}\)	\(4.2 \times 10^{-5}\)

The shear spring is the stiffest, giving \(\Delta t_c \approx 1.5 \times 10^{-5}\) s. In molecular dynamics, a good rule of thumb is to set \(\Delta t < 0.1 \Delta t_c\). Here, we use \(\Delta t = 10^{-8}\) s, which is well below the threshold.

Notes:

special_bonds, newton and comm_modify settings documented in How to BPM.
Pair style: pair_style zero — no contact needed for this problem
fix addforce to apply the loading
Time integration: fix nve/bpm/sphere (no damping, no dipole update)
Output: tip deflection \(\delta_x / L\) and \(\delta_y / L\) printed to file at regular intervals. The tip displacement is tracked using compute reduce ave on the right group.

The LAMMPS simulation script is provided below.

bend-dynamic-point-load.in ⬇ Download

# =========================================================================
# Cantilever beam under sinusoidal point load on free end
# 
# Replace NBONDS with the desired value before running.
# =========================================================================

variable    N       equal NBONDS

# Beam geometry
variable    L       equal 10.
variable    l0      equal $(v_L/v_N)
variable    w       equal 0.5
variable    h       equal 0.25

# 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 210e9
variable    nu      equal 0.3
variable    G       equal "v_E/(2*(1+v_nu))"

# Loading
variable    P0      equal 10e6
variable    omega   equal 50

# =========================================================================
# Initialization
# =========================================================================

units           si
dimension       3
boundary        f f f
processors      * 1 1

atom_style      bpm/sphere
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    once yes
comm_modify     mode single cutoff $(2*v_l0) vel yes

read_data       beam-${N}.lam

# =========================================================================
# Groups
# =========================================================================

region          left block INF 0.1 INF INF INF INF
region          right block 9.9 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 0 0 0

pair_style      zero $(2*v_l0)
pair_coeff      * * $(2*v_l0)

# =========================================================================
# Simulation time and output
# =========================================================================

variable        timestep equal 1e-8
variable        runtime  equal 1.5
variable        nsteps   equal v_runtime/v_timestep
variable        time     equal step*v_timestep

timestep        ${timestep}

compute         tx right reduce ave x
compute         ty right reduce ave y
variable        dx equal "c_tx/v_L-1"
variable        dy equal "c_ty/v_L"

compute         mytemp all temp/sphere
compute         erot all erotate/sphere

thermo          $(v_nsteps/1000)
thermo_style    custom step temp v_time v_dx v_dy
thermo_modify   temp mytemp

dump            1 all custom $(v_nsteps/100) ./bend-point-N${N}.dump &
                id x y z fx fy fz

fix             prnt all print $(v_nsteps/1000) &
                "${time} ${dx} ${dy}" &
                file deflection-N${N}.txt screen no &
                title "t dx/L dy/L"

# =========================================================================
# Dynamics
# =========================================================================

fix             nve free nve/bpm/sphere

variable        p equal v_P0*sin(v_omega*step*v_timestep)
fix             adf right addforce 0 v_p 0
fix_modify      adf energy no

run             ${nsteps}

Cantilever beam under static loading

These examples use a different beam geometry: a cylindrical rod with \(L = 100\), \(d = 1\), \(E = 1\), \(\nu = 0.5\) in dimensionless units (units lj). The applied loads are chosen such that the maximum tip deflection is 10% of the beam length, and the simulated shapes are compared against the exact Euler–Bernoulli deflection curves (Hibbeler 2018).

Key differences from the dynamic example:

Damping: fix langevin is used to reach quasi-static equilibrium. The load is first ramped from zero, then held constant while the beam relaxes.
Loading: point load (fix addforce on right) or uniform load (fix addforce on all free atoms, with total load distributed evenly).
Energy verification: the scripts compare the LAMMPS bond energy against the theoretical strain energy (\(U = P^2 L^3 / (6EI)\) for point load; \(U = w^2 L^5 / (40EI)\) for uniform load).

create_beam_for_static.py ⬇ Download

import numpy as np

N = 8  # number of bonds

# Dimensionless units
L = 100.
dx = L / N

with open('beam-{}.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, 101))
    file.write('%f %f ylo yhi \n' % (-101, 101))
    file.write('%f %f zlo zhi \n\n' % (-0.1, 0.1))

    # atom_style bpm/sphere
    file.write('Atoms \n\n')
    x = 0.
    for k in range(N+1):
        # atom-id atom-type molid diameter density x y z
        file.write('%i 1 1 1.0 1.0 %f 0 0 \n' % (k+1, x))
        x += dx
    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')

bend-static-point-load.in ⬇ Download

# =========================================================================
# Cantilever beam under static point load on free end
# 
# Replace NBONDS with the desired value before running.
# =========================================================================

variable    N       equal NBONDS

# Beam geometry (dimensionless)
variable    L       equal 100
variable    diam    equal 1
variable    l0      equal $(v_L/v_N)

# Cylindrical cross-section
variable    A       equal (PI/4)*v_diam^2
variable    I       equal (PI/64)*v_diam^4
variable    J       equal (PI/32)*v_diam^4

# Material properties
variable    E       equal 1
variable    nu      equal 0.5
variable    G       equal "v_E/(2*(1+v_nu))"

# Load chosen such that max tip deflection = 10% of L
variable    dymax   equal "-0.1*v_L"
variable    p       equal $(3*v_E*v_I*v_dymax/v_L^3)

# =========================================================================
# Initialization
# =========================================================================

units           lj
dimension       3
boundary        f f f
processors      * 1 1

atom_style      bpm/sphere
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    once yes
comm_modify     mode single cutoff $(2*v_l0) vel yes

read_data       beam-${N}.lam

change_box      all x final -500 500 y final -500 500 z final -500 500

# =========================================================================
# Groups
# =========================================================================

region          left block INF 0.1 INF INF INF INF
region          right block 99.9 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 &
                1e-2 1e-2 1e-2 1e-2

pair_style      zero $(2*v_l0)
pair_coeff      * * $(2*v_l0)

# =========================================================================
# Computes and output
# =========================================================================

timestep        0.01

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"

# Theoretical deflection curve (Euler-Bernoulli)
variable        ytheo atom "v_p*x^2/(6*v_E*v_I)*(3*v_L-x)"
variable        ydiff atom "y-v_ytheo"
compute         maxdiffy all reduce maxabs v_ydiff

# Energy
variable        Wext equal v_p*c_dely
variable        Uf_theo equal $(v_L^3*v_p^2/(6*v_E*v_I))

compute         mytemp all temp/sphere
compute         erot all erotate/sphere
variable        U equal "pe"
variable        Uerr equal abs(v_U-v_Uf_theo)/v_Uf_theo

thermo          1000000
thermo_style    custom step temp ke c_erot pe v_U v_Uf_theo v_Uerr v_Wext v_dy c_maxdiffy
thermo_modify   temp mytemp norm no

dump            1 all custom 1000000 ./bend-point-N${N}.dump &
                id x y z

# =========================================================================
# Stage 1: Ramp load
# =========================================================================

fix             nve free nve/bpm/sphere

variable        fy equal v_p*ramp(0,1)
fix             adf right addforce 0 v_fy 0
fix_modify      adf energy no

run             50000000

# =========================================================================
# Stage 2: Relax under constant load
# =========================================================================

fix             lan free langevin 0 0 1000 83715463 zero yes omega yes

variable        fy equal v_p
fix             adf right addforce 0 v_fy 0
fix_modify      adf energy no

run             50000000

# =========================================================================
# Post-processing
# =========================================================================

variable        Wextf_theo equal v_p*v_dymax

variable        fddy atom "v_p/(v_E*v_I)*(v_L-x)"
variable        Uf_atom atom "0.5*v_E*v_I*(v_fddy)^2*v_l0"
compute         Uf_discrete all reduce sum v_Uf_atom
run             0

print "======================================================================"
print "External work:       W = ${Wext}"
print "Theoretical:         W = ${Wextf_theo}"
print "======================================================================"
print "Elastic energy:      U = ${Uf_theo}"
print "Discrete approx:     U = $(c_Uf_discrete)"
print "======================================================================"

dump            final all custom 1 ./bend-point-N${N}-final.dump &
                id x y z
run             0

bend-static-uniform-load.in ⬇ Download

# =========================================================================
# Cantilever beam under static uniform load
# 
# Replace NBONDS with the desired value before running.
# =========================================================================

variable    N       equal NBONDS

# Beam geometry (dimensionless)
variable    L       equal 100
variable    diam    equal 1
variable    l0      equal $(v_L/v_N)

# Cylindrical cross-section
variable    A       equal (PI/4)*v_diam^2
variable    I       equal (PI/64)*v_diam^4
variable    J       equal (PI/32)*v_diam^4

# Material properties
variable    E       equal 1
variable    nu      equal 0.5
variable    G       equal "v_E/(2*(1+v_nu))"

# Load chosen such that max tip deflection = 10% of L
variable    dymax   equal "-0.1*v_L"
variable    w       equal 8*v_E*v_I*v_dymax/v_L^4

# =========================================================================
# Initialization
# =========================================================================

units           lj
dimension       3
boundary        f f f
processors      * 1 1

atom_style      bpm/sphere
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    once yes
comm_modify     mode single cutoff $(2*v_l0) vel yes

read_data       beam-${N}.lam

change_box      all x final -500 500 y final -500 500 z final -500 500

# =========================================================================
# Groups
# =========================================================================

region          left block INF 0.1 INF INF INF INF
region          right block 99.9 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 &
                1e-2 1e-2 1e-2 1e-2

pair_style      zero $(2*v_l0)
pair_coeff      * * $(2*v_l0)

# =========================================================================
# Computes and output
# =========================================================================

timestep        0.01

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"

# Distribute total load to free atoms
variable        Fy equal v_w*v_L/count(free)

# Theoretical deflection curve (Euler-Bernoulli)
variable        ytheo atom "v_w*x^2/(24*v_E*v_I)*(6*v_L^2-4*v_L*x+x^2)"
variable        ydiff atom "y-v_ytheo"
compute         maxdiffy all reduce maxabs v_ydiff

# Energy
variable        Wexta atom v_w*v_l0*y
compute         Wext all reduce sum v_Wexta
variable        Uf_theo equal $(v_L^5*v_w^2/(40*v_E*v_I))

compute         mytemp all temp/sphere
compute         erot all erotate/sphere
variable        U equal "pe"
variable        Uerr equal abs(v_U-v_Uf_theo)/v_Uf_theo

thermo          1000000
thermo_style    custom step temp ke c_erot pe v_U v_Uf_theo v_Uerr c_Wext v_dy c_maxdiffy
thermo_modify   temp mytemp norm no

dump            1 all custom 1000000 ./bend-uniform-N${N}.dump &
                id x y z

# =========================================================================
# Stage 1: Ramp load
# =========================================================================

fix             nve free nve/bpm/sphere

variable        fy equal v_Fy*ramp(0,1)
fix             adf free addforce 0 v_fy 0
fix_modify      adf energy no

run             50000000

# =========================================================================
# Stage 2: Relax under constant load
# =========================================================================

fix             lan free langevin 0 0 1000 83715463 zero yes omega yes

variable        fy equal v_Fy
fix             adf free addforce 0 v_fy 0
fix_modify      adf energy no

run             50000000

# =========================================================================
# Post-processing
# =========================================================================

variable        Wexta_theo atom v_w*v_l0*v_ytheo
compute         Wextf_theo all reduce sum v_Wexta_theo

variable        fddy atom "v_w/(v_E*v_I)*(0.5*v_L^2-v_L*x+0.5*x^2)"
variable        Uf_atom atom "0.5*v_E*v_I*(v_fddy)^2*v_l0"
compute         Uf_discrete all reduce sum v_Uf_atom
run             0

print "======================================================================"
print "External work:       W = $(c_Wext)"
print "Theoretical:         W = $(c_Wextf_theo)"
print "======================================================================"
print "Elastic energy:      U = ${Uf_theo}"
print "Discrete approx:     U = $(c_Uf_discrete)"
print "======================================================================"

dump            final all custom 1 ./bend-uniform-N${N}-final.dump &
                id x y z
run             0

References

Le, T.-N., Battini, J.-M., Hjiaj, M. (2011). Efficient formulation for dynamics of corotational 2D beams. Computational Mechanics, 48(2), 153–161. doi:10.1007/s00466-011-0585-6
Chen, X., Peng, D., Morrissey, J.P., Ooi, J.Y. (2022). A comparative assessment and unification of bond models in DEM simulations. Granular Matter, 24(1), 29. doi:10.1007/s10035-021-01187-2
Hibbeler, R.C. (2018). Mechanics of Materials. Pearson.

← Back to Examples