SPRING — two springs in series

Three co-linear nodes joined by two SPRING 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
from vtkmodules.util.vtkConstants import VTK_LINE

import femorph_solver
from femorph_solver import ELEMENTS

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. SPRING only carries an axial spring force; the transverse DOFs have zero stiffness and the solver’s zero-pivot guard pins them automatically.

points = np.array(
    [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]],
    dtype=np.float64,
)
cells = np.array([2, 0, 1, 2, 1, 2], dtype=np.int64)
cell_types = np.array([VTK_LINE, VTK_LINE], dtype=np.uint8)
grid = pv.UnstructuredGrid(cells, cell_types, points)
# Two real-constant ids: cell 0 uses real-id 1 (k1), cell 1 uses real-id 2 (k2).
grid.cell_data["ansys_real_constant"] = np.array([1, 2], dtype=np.int32)

m = femorph_solver.Model.from_grid(grid)
m.assign(ELEMENTS.SPRING, real=(k1,), real_id=1)
m.assign(ELEMENTS.SPRING, real=(k2,), real_id=2)

m.fix(nodes=[1], dof="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.fix(nodes=[nn], dof="UY")
    m.fix(nodes=[nn], dof="UZ")

m.apply_force(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()