energiesFlexibleBodiesTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Kinetic and potential (elastic) energy of flexible bodies as output variables (#2202): free,
  5#           undamped oscillations of two-element models of each beam and plate element, and a kinematic
  6#           tree under gravity with joint springs and constant joint forces. Each model gets sensors of
  7#           KineticEnergy and PotentialEnergy for every element, and the total is kept by the
  8#           generalized-alpha method without numerical damping (spectralRadius = 1) - up to the time
  9#           discretization of the nonlinear oscillation, which the printed relative change shows.
 10#           The elastic energy of an element is integrated with the strains and the integration rules
 11#           of its elastic forces, so that its derivative is exactly the elastic force.
 12#           (1) ANCFCable2D, (2) ANCFCable, (3) BeamGeometricallyExact2D, (4) BeamGeometricallyExact,
 13#           (5) ANCFBeam, (6) ANCFThinPlate, (7) KinematicTree: a revolute and a prismatic joint with
 14#           P control springs, constant joint forces and gravity.
 15#
 16# Author:   Johannes Gerstmayr
 17# Date:     2026-10-01
 18#
 19# 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.
 20#
 21#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 22
 23import exudyn as exu
 24from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 25from exudyn.shells import ShellMesh
 26
 27import numpy as np
 28
 29testIsActive = exu.sys.get('testIsActive', False)
 30
 31SC = exu.SystemContainer()
 32mbs = SC.AddSystem()
 33
 34KE = exu.OutputVariableType.KineticEnergy
 35PE = exu.OutputVariableType.PotentialEnergy
 36models = {}     #name: list of energy sensors; their sum is the total energy of the model
 37
 38def EnergySensors(name, bodies):
 39    models[name] = [mbs.AddSensor(SensorBody(bodyNumber=b, outputVariableType=ovt, storeInternal=True))
 40                    for b in bodies for ovt in [KE, PE]]
 41
 42L = 0.5         #element length
 43section = exu.BeamSection()
 44section.stiffnessMatrix = np.diag([100, 60, 80, 1, 1.5, 2])
 45section.inertia = np.diag([2e-3, 1e-3, 1e-3])
 46section.massPerLength = 2
 47
 48#(1) ANCFCable2D: three nodes, the middle one moving sideways
 49nodes = [mbs.AddNode(NodePoint2DSlope1(referenceCoordinates=[i*L, 0, 1, 0], initialVelocities=[0, 0.2 if i == 1 else 0, 0, 0]))
 50         for i in range(3)]
 51EnergySensors('ANCFCable2D', [mbs.AddObject(ObjectANCFCable2D(nodeNumbers=[nodes[i], nodes[i+1]], length=L,
 52              massPerLength=2, bendingStiffness=1, axialStiffness=100)) for i in range(2)])
 53
 54#(2) ANCFCable (3D)
 55nodes = [mbs.AddNode(NodePointSlope1(referenceCoordinates=[i*L, 1, 0, 1, 0, 0],
 56                                     initialVelocities=[0, 0.2 if i == 1 else 0, 0.1 if i == 2 else 0, 0, 0, 0]))
 57         for i in range(3)]
 58EnergySensors('ANCFCable', [mbs.AddObject(ObjectANCFCable(nodeNumbers=[nodes[i], nodes[i+1]], length=L,
 59              massPerLength=2, bendingStiffness=1, axialStiffness=100)) for i in range(2)])
 60
 61#(3) BeamGeometricallyExact2D
 62nodes = [mbs.AddNode(NodeRigidBody2D(referenceCoordinates=[i*L, 2, 0], initialVelocities=[0, 0.2 if i == 1 else 0, 0.5 if i == 0 else 0]))
 63         for i in range(3)]
 64EnergySensors('BeamGeometricallyExact2D', [mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers=[nodes[i], nodes[i+1]],
 65              length=L, massPerLength=2, crossSectionInertia=2e-3, bendingStiffness=1,
 66              axialStiffness=100, shearStiffness=60)) for i in range(2)])
 67
 68#(4) BeamGeometricallyExact (3D), Euler parameter nodes
 69nodes = [mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[i*L, 3, 0] + list(eulerParameters0),
 70                                     initialVelocities=[0, 0.2 if i == 1 else 0, 0.1 if i == 2 else 0] + list(
 71                                         AngularVelocity2EulerParameters_t([0.3, 0, 0] if i == 0 else [0, 0, 0], eulerParameters0))))
 72         for i in range(3)]
 73EnergySensors('BeamGeometricallyExact', [mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[nodes[i], nodes[i+1]],
 74              length=L, sectionData=section)) for i in range(2)])
 75
 76#(5) ANCFBeam (3D), with the cross-section deformation as a penalty
 77nodes = [mbs.AddNode(NodePointSlope23(referenceCoordinates=[i*L, 4, 0, 0, 1, 0, 0, 0, 1],
 78                                      initialVelocities=[0, 0.2 if i == 1 else 0, 0.1 if i == 2 else 0, 0, 0, 0, 0, 0, 0]))
 79         for i in range(3)]
 80EnergySensors('ANCFBeam', [mbs.AddObject(ObjectANCFBeam(nodeNumbers=[nodes[i], nodes[i+1]], length=L, sectionData=section,
 81              crossSectionPenaltyFactor=[1, 1, 1])) for i in range(2)])
 82
 83#(6) ANCFThinPlate: 2x1 elements, one corner moving out of the plane
 84plate = ShellMesh(vertices=[[0, 5, 0], [1, 5, 0], [1, 5.5, 0], [0, 5.5, 0]], numberOfElementsX=2, numberOfElementsY=1,
 85                  youngsModulus=2e5, poissonsRatio=0.3, density=500, thickness=0.02)
 86plate.CreateANCFThinPlateElements(mbs)
 87mbs.SetNodeParameter(plate.vertexNodeNumbers[2], 'initialVelocities', [0, 0, 0.1] + [0]*6)
 88EnergySensors('ANCFThinPlate', plate.elementNumbers)
 89
 90#(7) KinematicTree: a revolute joint about z with a link along x, then a prismatic joint along x; P control springs
 91#towards offsets, constant joint forces and gravity in -y; no D control, so the total of T and V is kept
 92nKT = mbs.AddNode(NodeGenericODE2(referenceCoordinates=[0, 0.2], initialCoordinates=[0.3, 0], initialCoordinates_t=[0, 0.1],
 93                                  numberOfODE2Coordinates=2))
 94oKT = mbs.AddObject(ObjectKinematicTree(nodeNumber=nKT, baseOffset=[0, 7, 0], linkParents=[-1, 0],
 95                                        jointTypes=[exu.JointType.RevoluteZ, exu.JointType.PrismaticX],
 96                                        jointTransformations=exu.Matrix3DList([np.eye(3)]*2),
 97                                        jointOffsets=exu.Vector3DList([[0, 0, 0], [0.5, 0, 0]]),
 98                                        linkInertiasCOM=exu.Matrix3DList([np.eye(3)*0.01]*2),
 99                                        linkCOMs=exu.Vector3DList([[0.25, 0, 0], [0.1, 0, 0]]),
100                                        linkMasses=[1, 2], gravity=[0, -9.81, 0],
101                                        jointPositionOffsetVector=[0.1, 0.2], jointPControlVector=[50, 200],
102                                        jointVelocityOffsetVector=[0, 0], jointDControlVector=[0, 0], jointForceVector=[2, -1]))
103EnergySensors('KinematicTree', [oKT])
104
105mbs.Assemble()
106
107simulationSettings = exu.SimulationSettings()
108simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
109simulationSettings.timeIntegration.numberOfSteps = 1000
110simulationSettings.timeIntegration.endTime = 1
111simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
112#the cross-section deformation of ObjectANCFBeam carries no mass, so its mass matrix is singular for initial accelerations
113simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = False
114simulationSettings.timeIntegration.verboseMode = 0
115simulationSettings.solution.file.write = False
116simulationSettings.solution.sensors.writePeriod = 0.01
117mbs.SolveDynamic(simulationSettings)
118
119testResult = 0
120for (name, sensors) in models.items():
121    total = sum(mbs.GetSensorStoredData(s)[:, 1] for s in sensors)
122    elastic = sum(mbs.GetSensorStoredData(s)[:, 1] for s in sensors[1::2])
123    exu.Print('%-28s total energy %.8f -> %.8f, relative change %.2e (largest elastic energy %.4f)' % (
124        name, total[0], total[-1], (total[-1]-total[0])/total[0], elastic.max()))
125    testResult += total[-1]
126
127exu.Print('solution of energiesFlexibleBodiesTest=', testResult)
128exu.sys['testResult'] = testResult