LINK180 — axial truss under end load

Single LINK180 spar fixed at node 1 and pulled along the x-axis at node 2. The tip displacement is compared to the closed-form rod solution u = P L / (E A) and the axial stress is plotted on the deformed mesh.

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver

Problem data

Steel rod, 1 m long, 100 mm² cross-section, pulled with a 1 kN tip load.

E = 2.1e11  # Pa
A = 1.0e-4  # m² (100 mm²)
L = 1.0  # m
P = 1.0e3  # N (tensile)

Build the model

femorph_solver.Model accepts MAPDL preprocessor verbs one-for-one: et, mp, r, n, e, d, f. A LINK180 has three translational DOFs per node, so we only fix UX at node 2’s support and zero-out the transverse DOFs to kill the transverse rigid modes.

m = femorph_solver.Model()
m.et(1, "LINK180")
m.mp("EX", 1, E)
m.mp("DENS", 1, 7850.0)
m.r(1, A)
m.n(1, 0.0, 0.0, 0.0)
m.n(2, L, 0.0, 0.0)
m.e(1, 2)

m.d(1, "ALL")  # clamp all translational DOFs at node 1
m.d(2, "UY")  # suppress transverse rigid-body motion
m.d(2, "UZ")
m.f(2, "FX", P)

/home/runner/_work/solver/solver/examples/elements/link180/example_link180.py:37: DeprecationWarning: Model.et(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).et(et_id, name)` for line-by-line APDL deck porting, or the native `Model.assign("HEX8", material)` for new code.
  m.et(1, "LINK180")
/home/runner/_work/solver/solver/examples/elements/link180/example_link180.py:38: DeprecationWarning: Model.mp(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).mp(prop, mat_id, value)` for line-by-line APDL deck porting, or the native `Model.assign(element, {prop: value, ...})` for new code.
  m.mp("EX", 1, E)
/home/runner/_work/solver/solver/examples/elements/link180/example_link180.py:39: DeprecationWarning: Model.mp(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).mp(prop, mat_id, value)` for line-by-line APDL deck porting, or the native `Model.assign(element, {prop: value, ...})` for new code.
  m.mp("DENS", 1, 7850.0)
/home/runner/_work/solver/solver/examples/elements/link180/example_link180.py:40: DeprecationWarning: Model.r(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).r(real_id, *values)` for line-by-line APDL deck porting, or the native `Model.assign(element, material, real=[...])` for new code.
  m.r(1, A)
/home/runner/_work/solver/solver/examples/elements/link180/example_link180.py:41: DeprecationWarning: Model.n(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).n(num, x, y, z)` for line-by-line APDL deck porting, or the native `Model.from_grid(pv_grid)` for new code.
  m.n(1, 0.0, 0.0, 0.0)
/home/runner/_work/solver/solver/examples/elements/link180/example_link180.py:42: DeprecationWarning: Model.n(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).n(num, x, y, z)` for line-by-line APDL deck porting, or the native `Model.from_grid(pv_grid)` for new code.
  m.n(2, L, 0.0, 0.0)
/home/runner/_work/solver/solver/examples/elements/link180/example_link180.py:43: DeprecationWarning: Model.e(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).e(*node_nums)` for line-by-line APDL deck porting, or the native `Model.from_grid(pv_grid)` for new code.
  m.e(1, 2)
/home/runner/_work/solver/solver/examples/elements/link180/example_link180.py:45: DeprecationWarning: Model.d(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).d(node, label, value)` for line-by-line APDL deck porting, or the native `Model.fix(nodes=..., where=..., dof=...)` for new code.
  m.d(1, "ALL")  # clamp all translational DOFs at node 1
/home/runner/_work/solver/solver/examples/elements/link180/example_link180.py:46: DeprecationWarning: Model.d(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).d(node, label, value)` for line-by-line APDL deck porting, or the native `Model.fix(nodes=..., where=..., dof=...)` for new code.
  m.d(2, "UY")  # suppress transverse rigid-body motion
/home/runner/_work/solver/solver/examples/elements/link180/example_link180.py:47: DeprecationWarning: Model.d(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).d(node, label, value)` for line-by-line APDL deck porting, or the native `Model.fix(nodes=..., where=..., dof=...)` for new code.
  m.d(2, "UZ")
/home/runner/_work/solver/solver/examples/elements/link180/example_link180.py:48: DeprecationWarning: Model.f(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).f(node, label, value)` for line-by-line APDL deck porting, or the native `Model.apply_force(node, fx=..., fy=..., fz=...)` for new code.
  m.f(2, "FX", P)

Static solve + analytical comparison

The rod equation gives u_tip = P L / (E A). With a single LINK180 this is exact (linear shape functions are sufficient for a prismatic bar under end load).

res = m.solve()

dof = m.dof_map()
tip_ux_row = np.where((dof[:, 0] == 2) & (dof[:, 1] == 0))[0][0]
u_tip = res.displacement[tip_ux_row]
u_expected = P * L / (E * A)

print(f"LINK180 tip UX        = {u_tip:.6e} m")
print(f"Analytical PL/(EA)    = {u_expected:.6e} m")
assert np.isclose(u_tip, u_expected, rtol=1e-10)

LINK180 tip UX        = 4.761905e-05 m
Analytical PL/(EA)    = 4.761905e-05 m

Post-processing: axial stress

For a single LINK180 the axial stress is uniform: σ = E · (Δu / L) where Δu is the elongation. Carry it as cell data on the mesh for plotting.

sigma_axial = E * (u_tip / L)
print(f"Axial stress          = {sigma_axial:.3e} Pa (= P/A = {P / A:.3e})")

grid = m.grid.copy()
displacement = np.zeros((grid.n_points, 3), dtype=np.float64)
for i, node_num in enumerate(grid.point_data["ansys_node_num"]):
    rows = np.where(dof[:, 0] == int(node_num))[0]
    for r in rows:
        displacement[i, int(dof[r, 1])] = res.displacement[r]
grid.point_data["displacement"] = displacement
grid.cell_data["sigma_axial"] = np.array([sigma_axial])

Axial stress          = 1.000e+07 Pa (= P/A = 1.000e+07)

Plot the deformed truss coloured by axial stress

warped = grid.warp_by_vector("displacement", factor=1.0e5)
plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(
    grid,
    style="wireframe",
    color="gray",
    line_width=3,
    label="undeformed",
)
plotter.add_mesh(
    warped,
    scalars="sigma_axial",
    line_width=6,
    scalar_bar_args={"title": "axial stress [Pa]"},
    label="deformed (×1e5)",
)
plotter.add_legend()
plotter.add_axes()
plotter.show()