geometricallyExactBeamMassTest.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  The mass matrix of ObjectBeamGeometricallyExact (#1273): the element-consistent mass matrix
 5#           represents a rigid motion exactly, whatever the mesh. (1) One element in a rigid rotation about
 6#           an oblique axis: the kinetic energy from the mass matrix equals the one of the rigid rod,
 7#           1.356 (the lumped mass matrix of the nodes gives 3.056). (2) A stiff pendulum of one element,
 8#           hinged at one end, small amplitude: its period equals the one of the rigid rod,
 9#           2*pi*sqrt(2L/(3g)) = 1.638 (lumped: 2.006 with one element, 1.664 with four).
10#
11# Author:   Johannes Gerstmayr
12# Date:     2026-09-30
13#
14# 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.
15#
16#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
17
18import exudyn as exu
19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
20
21import numpy as np
22
23testIsActive = exu.sys.get('testIsActive', False)
24
25#(1) kinetic energy of a rigid rotation
26L = 2; m = 3; J = np.diag([0.5, 0.2, 0.3])
27omega = np.array([0.3, -0.7, 1.1])
28SC = exu.SystemContainer(); mbs = SC.AddSystem()
29nodes = []
30for x in [0, L]:
31    velocity = np.cross(omega, [x-0.5*L, 0, 0])
32    nodes.append(mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[x, 0, 0]+eulerParameters0,
33                                             initialVelocities=list(velocity)+list(AngularVelocity2EulerParameters_t(omega, eulerParameters0)))))
34section = exu.BeamSection()
35section.stiffnessMatrix = np.diag([1e4]*6)
36section.inertia = J
37section.massPerLength = m/L
38mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=nodes, length=L, sectionData=section))
39mbs.Assemble()
40solver = exu.MainSolverImplicitSecondOrder()
41solver.InitializeSolver(mbs, exu.SimulationSettings())
42solver.ComputeMassMatrix(mbs)
43q_t = mbs.systemData.GetODE2Coordinates_t()
44kineticEnergy = 0.5*q_t @ solver.GetSystemMassMatrix() @ q_t
45inertiaRigid = np.diag([J[0, 0]*L, m*L**2/12+J[1, 1]*L, m*L**2/12+J[2, 2]*L])
46exu.Print('kinetic energy of the rigid rotation', kineticEnergy, ', of the rigid rod', 0.5*omega @ inertiaRigid @ omega)
47testResult = kineticEnergy
48
49#(2) stiff pendulum of one element, released at 0.05 rad from the vertical
50L = 1; m = 1; g = 9.81
51SC = exu.SystemContainer(); mbs = SC.AddSystem()
52rotation = RotationMatrixZ(-np.pi/2+0.05)
53axis = rotation @ np.array([1, 0, 0])
54nodes = [mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=list(x*axis)+list(RotationMatrix2EulerParameters(rotation)))) for x in [0, L]]
55section = exu.BeamSection()
56section.stiffnessMatrix = np.diag([1e9, 1e9, 1e9, 1e7, 1e7, 1e7])
57section.inertia = np.diag([2e-6, 1e-6, 1e-6])
58section.massPerLength = m/L
59mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=nodes, length=L, sectionData=section))
60for n in nodes:
61    mbs.AddLoad(LoadForceVector(markerNumber=mbs.AddMarker(MarkerNodePosition(nodeNumber=n)), loadVector=[0, -0.5*m*g, 0]))
62mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=mbs.AddNode(NodePointGround()), coordinate=0))
63for i in range(3):
64    mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mGround, mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nodes[0], coordinate=i))]))
65sTip = mbs.AddSensor(SensorNode(nodeNumber=nodes[-1], outputVariableType=exu.OutputVariableType.Position,
66                                storeInternal=True, writeToFile=False))
67mbs.Assemble()
68simulationSettings = exu.SimulationSettings()
69simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
70simulationSettings.timeIntegration.numberOfSteps = 4000
71simulationSettings.timeIntegration.endTime = 4
72simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
73simulationSettings.solution.file.write = False
74simulationSettings.solution.sensors.writePeriod = 1e-3
75mbs.SolveDynamic(simulationSettings)
76data = mbs.GetSensorStoredData(sTip)
77t = data[:, 0]; x = data[:, 1]
78crossings = [t[k]-x[k]*(t[k+1]-t[k])/(x[k+1]-x[k]) for k in range(len(x)-1) if x[k] > 0 and x[k+1] <= 0]
79period = crossings[1]-crossings[0]
80exu.Print('period of the stiff pendulum', period, ', of the rigid rod', 2*np.pi*np.sqrt(2*L/(3*g)))
81testResult += period
82
83exu.Print('solution of geometricallyExactBeamMassTest=', testResult)
84exu.sys['testResult'] = testResult