flexiblePendulumBeamComparison.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  A highly flexible pendulum, released from horizontal under gravity, modelled three times
  5#           with 10 elements: ObjectBeamGeometricallyExact2D, ObjectBeamGeometricallyExact (3D, moving in
  6#           the x-y plane) and ObjectANCFCable2D (#2730). Gravity acts as nodal forces, half the weight
  7#           of an element at each of its nodes, the same for all three. Measured: with 10 elements the
  8#           three differ from a 40-element GE2D solution by 9 mm (GE2D), 11 mm (GE3D) and 5 mm (ANCF2D)
  9#           in the tip position; with 40 elements GE3D agrees with GE2D to 0.1 mm and ANCF2D to 0.7 mm.
 10#           The result is the sum of the tip positions at the end, so that a change of any of the three
 11#           moves it.
 12#
 13# Author:   Johannes Gerstmayr
 14# Date:     2026-09-29
 15#
 16# 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.
 17#
 18#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 19
 20import exudyn as exu
 21from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 22
 23import numpy as np
 24
 25testIsActive = exu.sys.get('testIsActive', False)
 26
 27L = 0.5                 #length of the pendulum
 28nElements = 10
 29E = 1e8                 #a very soft elastomer
 30rho = 1000
 31h = 0.002               #height of the cross section (in y)
 32b = 0.01                #width (in z)
 33A = b*h
 34Izz = b*h**3/12         #bending in the x-y plane
 35Iyy = h*b**3/12
 36nu = 0.3
 37ks = 10*(1+nu)/(12+11*nu)   #shear correction
 38G = E/(2*(1+nu))
 39g = 9.81
 40lElement = L/nElements
 41tEnd = 1
 42stepSize = 0.0025
 43
 44
 45def Pendulum(kind):
 46    """the tip position over time"""
 47    SC = exu.SystemContainer()
 48    mbs = SC.AddSystem()
 49    nodes = []
 50    for i in range(nElements+1):
 51        x = i*lElement
 52        if kind == 'GE2D':
 53            nodes.append(mbs.AddNode(NodeRigidBody2D(referenceCoordinates=[x, 0, 0])))
 54        elif kind == 'GE3D':
 55            nodes.append(mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[x, 0, 0] + eulerParameters0)))
 56        else:
 57            nodes.append(mbs.AddNode(NodePoint2DSlope1(referenceCoordinates=[x, 0, 1, 0])))
 58    for i in range(nElements):
 59        if kind == 'GE2D':
 60            mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers=[nodes[i], nodes[i+1]], length=lElement,
 61                          massPerLength=rho*A, crossSectionInertia=rho*Izz, bendingStiffness=E*Izz,
 62                          axialStiffness=E*A, shearStiffness=ks*G*A))
 63        elif kind == 'GE3D':
 64            section = exu.BeamSection()
 65            section.stiffnessMatrix = np.diag([E*A, ks*G*A, ks*G*A, G*(Iyy+Izz), E*Iyy, E*Izz])
 66            section.inertia = np.diag([rho*(Iyy+Izz), rho*Iyy, rho*Izz])
 67            section.massPerLength = rho*A
 68            mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[nodes[i], nodes[i+1]], length=lElement,
 69                                                       sectionData=section))
 70        else:
 71            mbs.AddObject(ObjectANCFCable2D(nodeNumbers=[nodes[i], nodes[i+1]], length=lElement,
 72                          massPerLength=rho*A, bendingStiffness=E*Izz, axialStiffness=E*A))
 73    #gravity as nodal forces: half the weight of an element at each of its nodes
 74    for i in range(nElements+1):
 75        weight = rho*A*lElement*(0.5 if i in [0, nElements] else 1)*g
 76        mbs.AddLoad(LoadForceVector(markerNumber=mbs.AddMarker(MarkerNodePosition(nodeNumber=nodes[i])),
 77                                    loadVector=[0, -weight, 0]))
 78    #the hinge at the first node
 79    mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=mbs.AddNode(NodePointGround()), coordinate=0))
 80    for i in range(3 if kind == 'GE3D' else 2):
 81        mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mGround,
 82                      mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nodes[0], coordinate=i))]))
 83    sTip = mbs.AddSensor(SensorNode(nodeNumber=nodes[-1], outputVariableType=exu.OutputVariableType.Position,
 84                                    storeInternal=True, writeToFile=False))
 85    mbs.Assemble()
 86    simulationSettings = exu.SimulationSettings()
 87    simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
 88    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
 89    simulationSettings.timeIntegration.endTime = tEnd
 90    simulationSettings.timeIntegration.verboseMode = 0
 91    simulationSettings.solution.file.write = False
 92    simulationSettings.solution.sensors.writePeriod = stepSize
 93    simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
 94    mbs.SolveDynamic(simulationSettings)
 95    return mbs.GetSensorStoredData(sTip)
 96
 97
 98tips = dict((kind, Pendulum(kind)) for kind in ['GE2D', 'GE3D', 'ANCF2D'])
 99for kind in ['GE3D', 'ANCF2D']:
100    difference = np.abs(tips[kind][:, 1:3] - tips['GE2D'][:, 1:3]).max(axis=1)
101    exu.Print(kind + ' against GE2D: max tip difference', difference.max(), 'until t=0.5:', difference[tips['GE2D'][:, 0] <= 0.5].max())
102exu.Print('tip at t=1:', dict((kind, tips[kind][-1, 1:3].round(6).tolist()) for kind in tips))
103
104testResult = sum(tips[kind][-1, 1] + tips[kind][-1, 2] for kind in tips)
105exu.Print('solution of flexiblePendulumBeamComparison=', testResult)
106exu.sys['testResult'] = testResult