geometricallyExactBeamOutputTest.py

You can view and download this file on Github: geometricallyExactBeamOutputTest.py

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  The section forces, moments and strains of ObjectBeamGeometricallyExact (#2753): a cantilever
 5#           bent by large tip loads. (1) A bending moment, then a torsional moment at the tip: every element
 6#           carries it as TorqueLocal, and CurvatureLocal is the moment divided by the stiffness, exactly.
 7#           (2) A tip force: in every element, the section force ForceLocal rotated by the rotation matrix
 8#           at the middle of the element gives the tip force of 4.7 - to 5.6e-3 with 8 elements, 3.9e-4 with 32:
 9#           the section forces of an element are constant, the mean of a field that varies along it.
10#           Also the other output variables at the ends of an element, against its nodes.
11#
12# Author:   Johannes Gerstmayr
13# Date:     2026-09-30
14#
15# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
16#
17#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
18
19import exudyn as exu
20from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
21
22import numpy as np
23
24testIsActive = exu.sys.get('testIsActive', False)
25
26L = 1
27EA = 1e5; GA = 4e4; GJ = 2; EIy = 5; EIz = 3
28
29
30def Cantilever(force, torque, nElements=8):
31    """solve statically; return mbs, the elements and the nodes"""
32    lElement = L/nElements
33    SC = exu.SystemContainer(); mbs = SC.AddSystem()
34    nodes = [mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[i*lElement, 0, 0]+list(eulerParameters0)))
35             for i in range(nElements+1)]
36    section = exu.BeamSection()
37    section.stiffnessMatrix = np.diag([EA, GA, GA, GJ, EIy, EIz])
38    section.inertia = np.diag([2e-3, 1e-3, 1e-3])
39    section.massPerLength = 1
40    elements = [mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[nodes[i], nodes[i+1]], length=lElement,
41                                                           sectionData=section)) for i in range(nElements)]
42    mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=mbs.AddNode(NodePointGround()), coordinate=0))
43    for i in [0, 1, 2, 4, 5, 6]:
44        mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mGround, mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nodes[0], coordinate=i))]))
45    mTip = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodes[-1]))
46    mbs.AddLoad(LoadForceVector(markerNumber=mTip, loadVector=force))
47    mbs.AddLoad(LoadTorqueVector(markerNumber=mTip, loadVector=torque, bodyFixed=True))
48    mbs.Assemble()
49    simulationSettings = exu.SimulationSettings()
50    simulationSettings.staticSolver.numberOfLoadSteps = 10
51    simulationSettings.staticSolver.newton.relativeTolerance = 1e-10
52    simulationSettings.staticSolver.newton.absoluteTolerance = 1e-12
53    mbs.SolveStatic(simulationSettings)
54    return (mbs, elements, nodes)
55
56
57testResult = 0
58Out = lambda mbs, e, t, x=0: mbs.GetObjectOutputBody(e, t, localPosition=[x, 0, 0])
59
60#(1) tip moments, one at a time, given in the frame of the tip: constant along the beam
61for torque in [[0, 0, 0.5*EIz/L], [0.4*GJ/L, 0, 0]]:
62    (mbs, elements, nodes) = Cantilever([0, 0, 0], torque)
63    errorMoment = max(np.linalg.norm(Out(mbs, e, exu.OutputVariableType.TorqueLocal) - torque) for e in elements)
64    errorCurvature = max(np.linalg.norm(Out(mbs, e, exu.OutputVariableType.CurvatureLocal) - np.array(torque)/[GJ, EIy, EIz]) for e in elements)
65    exu.Print('tip moment', torque, ': largest error of TorqueLocal', errorMoment, ', of CurvatureLocal', errorCurvature)
66    testResult += sum(Out(mbs, elements[-1], exu.OutputVariableType.CurvatureLocal))
67
68#(2) tip force: every section carries it
69force = [0, -1.5*EIz/L**2, 0.5*EIy/L**2]
70for nElements in [32, 8]:
71    (mbs, elements, nodes) = Cantilever(force, [0, 0, 0], nElements)
72    errorForce = 0
73    for e in elements:
74        rotation = Out(mbs, e, exu.OutputVariableType.RotationMatrix).reshape(3, 3)
75        errorForce = max(errorForce, np.linalg.norm(rotation @ Out(mbs, e, exu.OutputVariableType.ForceLocal) - force))
76    exu.Print(nElements, 'elements, tip force: largest error of the section force, rotated to global', errorForce)
77testResult += sum(Out(mbs, elements[0], exu.OutputVariableType.TorqueLocal))
78testResult += sum(Out(mbs, elements[0], exu.OutputVariableType.StrainLocal))
79
80#the output at the ends of an element is the one of its nodes
81errorEnds = 0
82for (x, node) in [(-0.5*L/nElements, nodes[-2]), (0.5*L/nElements, nodes[-1])]:
83    for outputVariable in [exu.OutputVariableType.Position, exu.OutputVariableType.RotationMatrix,
84                           exu.OutputVariableType.AngularVelocity]:
85        errorEnds = max(errorEnds, np.linalg.norm(Out(mbs, elements[-1], outputVariable, x) - mbs.GetNodeOutput(node, outputVariable)))
86exu.Print('element ends against the nodes: largest difference', errorEnds)
87rotation = Out(mbs, elements[-1], exu.OutputVariableType.Rotation, 0.5*L/nElements)
88testResult += sum(rotation)
89
90exu.Print('solution of geometricallyExactBeamOutputTest=', testResult)
91exu.sys['testResult'] = testResult