You can view and download this file on Github: ANCFmovingRigidBodyTest.py
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: Test model for moving rigid body on two cables
#
# Author: Andreas Zwölfer, Johannes Gerstmayr
# Date: 2019-12-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.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
import exudyn as exu
from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
import exudyn.graphics as graphics #only import if it does not conflict
useGraphics = True #without test
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
try: #only if called from test suite
from modelUnitTests import exudynTestGlobals #for globally storing test results
useGraphics = exudynTestGlobals.useGraphics
except:
class ExudynTestGlobals:
pass
exudynTestGlobals = ExudynTestGlobals()
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SC = exu.SystemContainer()
mbs = SC.AddSystem()
#Import function
#import GenerateStraightLineANCFCable2D as func
#background
rect = [-2,-2,4,2] #xmin,ymin,xmax,ymax
background0 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[rect[0],rect[1],0, rect[2],rect[1],0, rect[2],rect[3],0, rect[0],rect[3],0, rect[0],rect[1],0]} #background
background1 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[0,-1,0, 2,-1,0]} #background
oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0]))
L=10
gravityFieldConstant=9.81 #9000
complianceFactBend = 0.1
complianceFactAxial = 0.1
nEl=10 #must be even number
vALE=2*5
offset=-1
nGlobalGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
mGlobalGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGlobalGround, coordinate=0))
fixANCFRotation = 1
#######################SUSPENSION ROPE##################################################################################################################################################################
suspensionCableTemplate=Cable2D(physicsMassPerLength=20.87,
physicsBendingStiffness=78878*complianceFactBend,
physicsAxialStiffness=398240000*complianceFactAxial)
[suspensionCableNodeList, suspensionCableObjectList, suspensionLoadList, suspensionCableNodePositionList, dummy]=GenerateStraightLineANCFCable2D(mbs=mbs, positionOfNode0=[0,0,0], positionOfNode1=[L,0,0], numberOfElements=nEl, cableTemplate=suspensionCableTemplate,
massProportionalLoad=[0,-gravityFieldConstant,0], fixedConstraintsNode0=[1,1,0,fixANCFRotation], fixedConstraintsNode1=[1,1,0,fixANCFRotation])
##################################################################################################################################################################
######################Haulage ROPE##################################################################################################################################################################
nALE = mbs.AddNode(NodeGenericODE2(numberOfODE2Coordinates=1, referenceCoordinates=[0], initialCoordinates=[0], initialCoordinates_t=[vALE]))
mALE = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nALE, coordinate=0)) #ALE velocity marker
haulageCableTemplate=ALECable2D(physicsMassPerLength=6.96,
physicsBendingStiffness=5956*complianceFactBend,
physicsAxialStiffness=96725000*complianceFactAxial,
physicsAddALEvariation=False) #for compatibility with test suite results
haulageCableTemplate.nodeNumbers[2]=nALE #this will not be overwritten!
[haulageCableNodeList, haulageCableObjectList, haulageLoadList, haulageCableNodePositionList, dummy]=GenerateStraightLineANCFCable2D(mbs=mbs,
positionOfNode0=[0,offset,0], positionOfNode1=[L,offset,0], numberOfElements=nEl, cableTemplate=haulageCableTemplate,
massProportionalLoad=[0,-gravityFieldConstant,0], fixedConstraintsNode0=[1,1,0,fixANCFRotation], fixedConstraintsNode1=[1,1,0,fixANCFRotation])
cAleConstraint=mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mALE]))
##################################################################################################################################################################################################################################################
#slack carrier test
slackCarrierWheelRadius=0.75
mSuspensionRopeAttachmentNodeX=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = suspensionCableNodeList[int(nEl/2)], coordinate=0))
mSuspensionRopeAttachmentNodeY=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = suspensionCableNodeList[int(nEl/2)], coordinate=1))
graphicsSlackCarrier={'type':'Circle', 'color':[.1,0.1,0.8,1], 'position':[0,0,0], 'radius': slackCarrierWheelRadius}
nSlackCarrierRigidBody = mbs.AddNode(Rigid2D(referenceCoordinates=[5,offset-slackCarrierWheelRadius,0]))
oSlackCarrierRigidBody = mbs.AddObject(RigidBody2D(physicsMass=1, physicsInertia=1, nodeNumber=nSlackCarrierRigidBody,visualization=VObjectRigidBody2D(graphicsData= [graphicsSlackCarrier]))) #, visualization=VObjectRigidBody2D(graphicsData= [graphicsSupportWheels])
mSlackCarrierX = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nSlackCarrierRigidBody,coordinate=0))
mSlackCarrierY = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nSlackCarrierRigidBody,coordinate=1))
mSlackCarrierRot = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nSlackCarrierRigidBody,coordinate=2))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mSuspensionRopeAttachmentNodeX,mSlackCarrierX]))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mSuspensionRopeAttachmentNodeY,mSlackCarrierY]))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mSlackCarrierRot]))
nSegments = 4 #number of contact segments; must be consistent between nodedata and contact element
useFriction = False
nFactFriction = 1
if useFriction: nFactFriction = 3
initialGapList = [0.1]*(nSegments*nFactFriction) #initial gap of 0.1
cStiffness = 1e7
mContactSlackCarrier=mbs.AddMarker(MarkerBodyRigid(bodyNumber = oSlackCarrierRigidBody))
for i in haulageCableObjectList:
mContactCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=i, numberOfSegments = nSegments))
nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments*nFactFriction))
if useFriction:
mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mContactSlackCarrier, mContactCable], nodeNumber = nodeDataContactCable, numberOfContactSegments=nSegments, contactStiffness = cStiffness, circleRadius = slackCarrierWheelRadius, offset = 0))
else:
mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mContactSlackCarrier, mContactCable], nodeNumber = nodeDataContactCable, numberOfContactSegments=nSegments, contactStiffness = cStiffness, circleRadius = slackCarrierWheelRadius, offset = 0))
##################################################################################
a = 0.8 #y-dim/2 of gondula
b = 0.02 #x-dim/2 of gondula
yCOM = a #COM distance to attachment point; in vertical direction
massRigid = 60 #12
#vInit=40
inertiaRigid = massRigid/12*(2*a)**2
g = 9.81 # gravity
refPos = [0,offset,0]
# refPos = [fieldData['stationData'][0]['referencePointCoordinates'][1][0],fieldData['maxVerticalPositionSuspensionRopeShoes'][0],0]
#rigid body which slides:
graphicsRigid1 = GraphicsDataRectangle(-b,0,b,a) #drawing of rigid body
graphicsRigid2 = GraphicsDataRectangle(-a,-a,a,0) #drawing of rigid body
nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[refPos[0],refPos[1]-yCOM,0], initialVelocities=[vALE,0,0]));
oRigid = mbs.AddObject(RigidBody2D(physicsMass=massRigid, physicsInertia=inertiaRigid,nodeNumber=nRigid,visualization=VObjectRigidBody2D(graphicsData= [graphicsRigid1,graphicsRigid2])))
markerRigidTopAle = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[0.,yCOM,0.])) #support point
mR2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.,0.,0.])) #center of mass (for load)
mbs.AddLoad(Force(markerNumber=mR2, loadVector=[0,-massRigid*g,0]))
aleCableMarkerList = []#list of Cable2DCoordinates markers
aleOffsetList = [] #list of offsets counted from first cable element; needed in sliding joint
aleOffset = 0 #first cable element has offset 0
aleSlidingCoordinateInit = 0
aleInitialLocalMarker = 0
for i in range(len(haulageCableObjectList)): #create markers for cable elements
m = mbs.AddMarker(MarkerBodyCable2DCoordinates(bodyNumber = haulageCableObjectList[i]))
aleCableMarkerList += [m] #list containing 'MarkerBodyCable2DCoordinates' marker for sliding joint
aleOffsetList += [aleOffset] #list of relative (arclength) coordinates of the starting point of a cable
aleOffset += L/nEl
nodeDataAleSlidingJoint = mbs.AddNode(NodeGenericData(initialCoordinates=[aleInitialLocalMarker],numberOfDataCoordinates=1)) #initial index in cable list
aleSlidingJoint = mbs.AddObject(ObjectJointALEMoving2D(name='aleSlider', markerNumbers=[markerRigidTopAle,aleCableMarkerList[aleInitialLocalMarker]],
slidingMarkerNumbers=aleCableMarkerList, slidingMarkerOffsets=aleOffsetList,nodeNumbers=[nodeDataAleSlidingJoint, nALE], activeConnector = False))
mnRigid0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid, coordinate=0)) #add rigid body marker
mnRigid1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid, coordinate=1)) #add rigid body marker
mnRigid2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid, coordinate=2)) #add rigid body marker
nCCRigid0 = mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mnRigid0]))
nCCRigid1 = mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mnRigid1]))
nCCRigid2 = mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mnRigid2]))
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# Assemble multibody system
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
mbs.Assemble()
#exu.Print(mbs)
#mbs.WaitForUserToContinue()
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# Simualtion settings:
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
simulationSettings = exu.SimulationSettings()
simulationSettings.staticSolver.numberOfLoadSteps = 2
SC.visualizationSettings.general.circleTiling = 64
SC.visualizationSettings.nodes.defaultSize=0.125
SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.Displacement
SC.visualizationSettings.contour.outputVariableComponent = 1 # plot y-component
SC.visualizationSettings.contact.contactPointsDefaultSize = .005
SC.visualizationSettings.connectors.showContact = True
if useGraphics:
exu.StartRenderer()
#get initial velocities
vInit = mbs.systemData.GetODE2Coordinates_t(configuration = exu.ConfigurationType.Initial)
#mbs.WaitForUserToContinue()
mbs.SolveStatic(simulationSettings)
#prolong solution for next computation
u = mbs.systemData.GetODE2Coordinates()
#v = mbs.systemData.GetODE2Coordinates_t()
data = mbs.systemData.GetDataCoordinates()
mbs.systemData.SetODE2Coordinates(u,configuration = exu.ConfigurationType.Initial)
mbs.systemData.SetODE2Coordinates_t(vInit,configuration = exu.ConfigurationType.Initial)
mbs.systemData.SetDataCoordinates(data,configuration = exu.ConfigurationType.Initial)
#store some reference value:
ncables = len(suspensionCableNodeList)
sol = mbs.systemData.GetODE2Coordinates();
u = sol[int(ncables/4)*4+1]; #y-displacement of node at midpoint of rope
#mbs.WaitForUserToContinue()
mbs.SetObjectParameter(aleSlidingJoint, 'activeConnector', True)
mbs.SetObjectParameter(nCCRigid0, 'activeConnector', False)
mbs.SetObjectParameter(nCCRigid1, 'activeConnector', False)
mbs.SetObjectParameter(nCCRigid2, 'activeConnector', False)
mbs.SetObjectParameter(cAleConstraint, 'activeConnector', False)
solveDynamic = True
if solveDynamic:
# time related settings:
steps=200
tend=0.1
h=tend/steps
#fact = 15000
simulationSettings.timeIntegration.numberOfSteps = steps #1*fact
simulationSettings.timeIntegration.endTime = tend #0.002*fact
# Integrator related settings:
# simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = False
# simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False
simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.3
simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = False
simulationSettings.timeIntegration.verboseMode = 1
simulationSettings.timeIntegration.verboseModeFile = 0
mbs.SolveDynamic(simulationSettings)
if useGraphics:
#SC.WaitForRenderEngineStopFlag()
exu.StopRenderer() #safely close rendering window!
#compute error for test suite:
ncables = len(suspensionCableNodeList)
sol2 = mbs.systemData.GetODE2Coordinates();
u2 = sol2[int(ncables/4)*4+1]; #y-displacement of node in first quater of rope
exu.Print('static deflection =',u) #2020-03-05(corrected Cable2DshapeMarker): -0.06446474690480661 2019-12-17(new static solver): -0.06446474690512931; 2019-12-16: -0.06446474679809994
exu.Print('dynamic deflection =',u2) #2020-03-05(corrected Cable2DshapeMarker):0.06446627698400298; 2020-01-09: -0.06446627698121662(computeInitialAccelerations = False) 2020-01-09: -0.06446627843202835; 2019-12-26: -0.06446627698104967; 2019-12-17(update residual): -0.06446627698121662; 2019-12-16 (late): -0.06446627699890756; 2019-12-16: -0.06446610364603222
#exudynTestGlobals.testError = u + u2 - (-0.06446474690480661-0.06446627698400298)
exu.Print('ANCFmovingRigidBodyTest=',u+u2)
exudynTestGlobals.testError = u + u2 - (-0.06446474690612931 - 0.06446622244370685) #updated 2022-12-25
exudynTestGlobals.testResult = u + u2