geometricallyExactBeamElbowCantilever.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  The right-angle (elbow) cantilever of Simo and Vu-Quoc (1988) with ObjectBeamGeometricallyExact:
 5#           two legs of length 10, the first along y and clamped at the origin, the second along x, joined
 6#           rigidly at the elbow; EA = GA = 1e6, GJ = EI = 1e3, rho*A = 1, rho*I = diag(20, 10, 10). A force in
 7#           z at the elbow rises from 0 to 50 at t = 1 and falls back to 0 at t = 2; then the cantilever
 8#           oscillates freely in combined bending and torsion with amplitudes of the order of the leg length.
 9#           The literature gives this benchmark as curves only (Simo and Vu-Quoc 1988, Fig. 8). Measured with
10#           10 elements per leg and step 0.05 until t = 30: out-of-plane displacement of the tip 8.24 at
11#           t = 7.40 and -9.75 at t = 15.85, of the elbow -4.07 at t = 10.85 - the extrema of the published
12#           curves; 40 elements with step 0.005 give 8.25 at 7.44, -9.73 at 15.87 and -4.06 at 10.87, and
13#           differ from the test by at most 0.24 (tip) and 0.14 (elbow) over the whole time.
14#
15# Author:   Johannes Gerstmayr
16# Date:     2026-09-30
17#
18# 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.
19#
20#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
21
22import exudyn as exu
23from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
24
25import numpy as np
26
27testIsActive = exu.sys.get('testIsActive', False)
28
29L = 10                  #length of each leg
30nElements = 10          #per leg
31stepSize = 0.05
32tEnd = 30
33
34SC = exu.SystemContainer()
35mbs = SC.AddSystem()
36
37section = exu.BeamSection()
38section.stiffnessMatrix = np.diag([1e6, 1e6, 1e6, 1e3, 1e3, 1e3])
39section.inertia = np.diag([20, 10, 10])
40section.massPerLength = 1
41lElement = L/nElements
42
43rotationLeg1 = RotationMatrixZ(np.pi/2)  #local x along global y
44nodesLeg1 = [mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[0, i*lElement, 0]+list(RotationMatrix2EulerParameters(rotationLeg1))))
45             for i in range(nElements+1)]
46nodesLeg2 = [mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[i*lElement, L, 0]+list(eulerParameters0)))
47             for i in range(nElements+1)]
48for nodes in [nodesLeg1, nodesLeg2]:
49    for i in range(nElements):
50        mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[nodes[i], nodes[i+1]], length=lElement, sectionData=section))
51
52mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=mbs.AddNode(NodePointGround())))
53#the clamping in the frame of the first leg, the corner turned back: rotations of the markers (localHT)
54mGroundLeg1 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=mbs.GetMarker(mGround)['nodeNumber'],
55                                            localHT=exu.HT(rotation=rotationLeg1)))
56mbs.AddObject(GenericJoint(markerNumbers=[mGroundLeg1, mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodesLeg1[0]))]))
57mbs.AddObject(GenericJoint(markerNumbers=[mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodesLeg1[-1],
58                                                                        localHT=exu.HT(rotation=rotationLeg1.T))),
59                                          mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodesLeg2[0]))]))
60
61def LoadElbow(mbs, t, loadVector):
62    force = 50*t if t < 1 else (50*(2-t) if t < 2 else 0)
63    return [0, 0, force]
64
65mbs.AddLoad(LoadForceVector(markerNumber=mbs.AddMarker(MarkerNodePosition(nodeNumber=nodesLeg1[-1])), loadVector=[0, 0, 0],
66                            loadVectorUserFunction=LoadElbow))
67sElbow = mbs.AddSensor(SensorNode(nodeNumber=nodesLeg1[-1], outputVariableType=exu.OutputVariableType.Displacement,
68                                  storeInternal=True, writeToFile=False))
69sTip = mbs.AddSensor(SensorNode(nodeNumber=nodesLeg2[-1], outputVariableType=exu.OutputVariableType.Displacement,
70                                storeInternal=True, writeToFile=False))
71mbs.Assemble()
72
73simulationSettings = exu.SimulationSettings()
74simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
75simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
76simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
77simulationSettings.timeIntegration.endTime = tEnd
78simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1   #no numerical damping
79simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8
80simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-8
81simulationSettings.solution.file.write = False
82simulationSettings.solution.sensors.writePeriod = stepSize
83
84mbs.SolveDynamic(simulationSettings)
85
86elbow = mbs.GetSensorStoredData(sElbow)
87tip = mbs.GetSensorStoredData(sTip)
88testResult = 0
89for (name, data, extremum) in [('tip', tip, np.argmax), ('tip', tip, np.argmin), ('elbow', elbow, np.argmin)]:
90    i = extremum(data[:, 3])
91    exu.Print(name, 'out-of-plane displacement', round(data[i, 3], 4), 'at t =', round(data[i, 0], 2))
92    testResult += data[i, 3]
93
94exu.Print('solution of geometricallyExactBeamElbowCantilever=', testResult)
95exu.sys['testResult'] = testResult
96
97if not testIsActive:
98    mbs.PlotSensor(sensorNumbers=[sElbow, sTip], components=[2, 2], labels=['elbow', 'tip'],
99                   yLabel='out-of-plane displacement')