COMBIN14 — two springs in series

Three co-linear nodes joined by two COMBIN14 springs of stiffness k1 and k2. With the far end fixed and a tip load F the free-end displacement is F / k_eq where k_eq = k1 · k2 / (k1 + k2).

from __future__ import annotations

import numpy as np
import pyvista as pv

import femorph_solver

Problem data

Two springs. k1 at half the stiffness of k2; their series combination sits in between.

k1 = 1000.0  # N/m
k2 = 2000.0  # N/m
F = 50.0  # N tip load

k_eq = k1 * k2 / (k1 + k2)

Build the model

Three nodes along x. COMBIN14 only carries an axial spring force; the transverse DOFs have zero stiffness and the solver’s zero-pivot guard pins them automatically (equivalent to MAPDL’s auto-constrain warning).

m = femorph_solver.Model()
m.et(1, "COMBIN14")

m.r(1, k1)
m.r(2, k2)

m.n(1, 0.0, 0.0, 0.0)
m.n(2, 1.0, 0.0, 0.0)
m.n(3, 2.0, 0.0, 0.0)

m.real(1)
m.e(1, 2)
m.real(2)
m.e(2, 3)

m.d(1, "ALL")  # clamp node 1
# Restrain transverse translations on the free nodes so the DOFs the
# solver actually keeps are just UX at nodes 2 and 3.
for nn in (2, 3):
    m.d(nn, "UY")
    m.d(nn, "UZ")

m.f(3, "FX", F)

/home/runner/_work/solver/solver/examples/elements/combin14/example_combin14.py:38: 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, "COMBIN14")
/home/runner/_work/solver/solver/examples/elements/combin14/example_combin14.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, k1)
/home/runner/_work/solver/solver/examples/elements/combin14/example_combin14.py:41: 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(2, k2)
/home/runner/_work/solver/solver/examples/elements/combin14/example_combin14.py:43: 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/combin14/example_combin14.py:44: 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, 1.0, 0.0, 0.0)
/home/runner/_work/solver/solver/examples/elements/combin14/example_combin14.py:45: 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(3, 2.0, 0.0, 0.0)
/home/runner/_work/solver/solver/examples/elements/combin14/example_combin14.py:47: DeprecationWarning: Model.real(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).real(real_id)` for line-by-line APDL deck porting, or the native `Model.assign("HEX8", material, real=[...])` for new code.
  m.real(1)
/home/runner/_work/solver/solver/examples/elements/combin14/example_combin14.py:48: 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/combin14/example_combin14.py:49: DeprecationWarning: Model.real(...) is a MAPDL-dialect shortcut and has moved off the Model public surface.  Use `APDL(model).real(real_id)` for line-by-line APDL deck porting, or the native `Model.assign("HEX8", material, real=[...])` for new code.
  m.real(2)
/home/runner/_work/solver/solver/examples/elements/combin14/example_combin14.py:50: 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(2, 3)
/home/runner/_work/solver/solver/examples/elements/combin14/example_combin14.py:52: 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 node 1
/home/runner/_work/solver/solver/examples/elements/combin14/example_combin14.py:56: 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(nn, "UY")
/home/runner/_work/solver/solver/examples/elements/combin14/example_combin14.py:57: 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(nn, "UZ")
/home/runner/_work/solver/solver/examples/elements/combin14/example_combin14.py:59: 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(3, "FX", F)

Static solve + analytical comparison

With both springs in equilibrium the tip displacement must equal F / k_eq. Because the springs are in series the force is the same in each and we can check the intermediate displacement too (u2 = F / k1).

res = m.solve()
dof = m.dof_map()

u2 = res.displacement[np.where((dof[:, 0] == 2) & (dof[:, 1] == 0))[0][0]]
u3 = res.displacement[np.where((dof[:, 0] == 3) & (dof[:, 1] == 0))[0][0]]

print(f"u at node 2 = {u2:.6e} m (expected {F / k1:.6e})")
print(f"u at node 3 = {u3:.6e} m (expected {F / k_eq:.6e})")

assert np.isclose(u2, F / k1, rtol=1e-12)
assert np.isclose(u3, F / k_eq, rtol=1e-12)

u at node 2 = 5.000000e-02 m (expected 5.000000e-02)
u at node 3 = 7.500000e-02 m (expected 7.500000e-02)

Visualise the deformation

Build a point_data displacement vector and warp the mesh. Because the real deflection is ~7.5 cm on a 2 m baseline a modest exaggeration makes the difference between the two springs visible.

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

warped = grid.warp_by_vector("displacement", factor=1.0)
plotter = pv.Plotter(off_screen=True)
plotter.add_mesh(grid, style="wireframe", color="gray", line_width=3)
plotter.add_mesh(
    warped,
    scalars=np.linalg.norm(displacement, axis=1),
    line_width=6,
    render_points_as_spheres=True,
    point_size=14,
    scalar_bar_args={"title": "|u| [m]"},
)
plotter.add_axes()
plotter.show()