.. _testmodels-ancfslidingandalejointtest: ***************************** ANCFslidingAndALEjointTest.py ***************************** You can view and download this file on Github: `ANCFslidingAndALEjointTest.py `_ .. code-block:: python :linenos: #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # This is an EXUDYN example # # Details: Test model for rigid bodies sliding on cables # # Author: Andreas Zwölfer, Johannes Gerstmayr # Date: 2019-12-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. # #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 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() ################################################################################################################################################################## def AddBodyWithSlidingJoints(mbs,xPositionOfFirstNodes=0,referencePositionOfBodyAlongCable=0,gravityFieldConstant=0): a = 0.8*2 #y-dim/2 of Body b = 0.02 #x-dim/2 of Body massRigid = 5000 #moving mass inertiaRigid = massRigid/12*(2*a)**2 #rigid body which slides: graphicsRigid1 = GraphicsDataRectangle(-b,-a,b,a) #drawing of rigid body yCOM = a #COM distance to attachment point on suspension rope; in vertical direction nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[xPositionOfFirstNodes+referencePositionOfBodyAlongCable,-yCOM,0])); oRigid = mbs.AddObject(RigidBody2D(physicsMass=massRigid, physicsInertia=inertiaRigid,nodeNumber=nRigid,visualization=VObjectRigidBody2D(graphicsData= [graphicsRigid1]))) markerRigidTop=mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[0.,yCOM,0.])) #support point markerRigidTopAle = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[0.,yCOM-offset,0.])) #support point if gravityFieldConstant: mR2 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.,0.,0.])) #center of mass (for load) mbs.AddLoad(Force(markerNumber=mR2, loadVector=[0,-massRigid*gravityFieldConstant,0])) #find index of cable element associated with "referencePositionOfBodyAlongCable" cummulativeLength=0 count=0 while referencePositionOfBodyAlongCable >= cummulativeLength: cummulativeLength+=mbs.GetObjectParameter(cable2ObjectList[count],'physicsLength') count+=1 count+=-1 #AleSlidingJoint: oAleSlidingJoint=GenerateAleSlidingJoint(mbs,cable1ObjectList,markerRigidTopAle,AleNode=nALE,localMarkerIndexOfStartCable=count,AleSlidingOffset=referencePositionOfBodyAlongCable)[0] #slidingJoint: oSlidingJoint=GenerateSlidingJoint(mbs,cable2ObjectList,markerRigidTop,localMarkerIndexOfStartCable=count,slidingCoordinateStartPosition=referencePositionOfBodyAlongCable)[0] return [oRigid,nRigid, oAleSlidingJoint, oSlidingJoint] #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=100 #length of ropes gravityFieldConstant=9.81 complianceFactBend = 1 complianceFactAxial = 1 nEl=6 #must be even number vALE=0 offset=0.35 #distance between the two ropes offsetCarrier=-0.515 #distance between suspension rope and slack carrier wheel nGlobalGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) mGlobalGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGlobalGround, coordinate=0)) fixANCFRotation = 0 #######################ROPE2 (Carrier rope)######################################################################################################################## cable2Template=Cable2D(physicsMassPerLength=10, physicsBendingStiffness=50000*complianceFactBend, physicsAxialStiffness=2e8*complianceFactAxial) [cable2NodeList, cable2ObjectList, suspensionLoadList, cable2NodePositionList, dummy]=GenerateStraightLineANCFCable2D(mbs=mbs, positionOfNode0=[0,0,0], positionOfNode1=[L,0,0], numberOfElements=nEl, cableTemplate=cable2Template, massProportionalLoad=[0,-gravityFieldConstant,0], fixedConstraintsNode0=[1,1,0,fixANCFRotation], fixedConstraintsNode1=[1,1,0,fixANCFRotation]) ################################################################################################################################################################## ######################ROPE1####################################################################################################################################### nALE = mbs.AddNode(NodeGenericODE2(numberOfODE2Coordinates=1, referenceCoordinates=[0], initialCoordinates=[0], initialCoordinates_t=[vALE])) mALE = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nALE, coordinate=0)) #ALE velocity marker cable1Template=ALECable2D(physicsMassPerLength=3, physicsBendingStiffness=4000*complianceFactBend, physicsAxialStiffness=5e7*complianceFactAxial,physicsUseCouplingTerms=False, physicsAddALEvariation=False) #for compatibility with test suite results cable1Template.nodeNumbers[2]=nALE [cable1NodeList, cable1ObjectList, haulageLoadList, cable1NodePositionList, dummy]=GenerateStraightLineANCFCable2D(mbs=mbs, positionOfNode0=[0,-offset,0], positionOfNode1=[L,-offset,0], numberOfElements=nEl, cableTemplate=cable1Template, massProportionalLoad=[0,-gravityFieldConstant,0], fixedConstraintsNode0=[1,1,0,fixANCFRotation], fixedConstraintsNode1=[1,1,0,fixANCFRotation]) cAleConstraint=mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mALE])) ################################################################################################################################################################### #Slack carrier: carrierWheelRadius=0.42/2 mSuspensionRopeAttachmentNodeX=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = cable2NodeList[int(nEl/2)], coordinate=0)) mSuspensionRopeAttachmentNodeY=mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = cable2NodeList[int(nEl/2)], coordinate=1)) graphicsCarrier={'type':'Circle', 'color':[.1,0.1,0.8,1], 'position':[0,0,0], 'radius': carrierWheelRadius} nCarrierRigidBody = mbs.AddNode(Rigid2D(referenceCoordinates=[L/2,offsetCarrier-carrierWheelRadius,0])) oCarrierRigidBody = mbs.AddObject(RigidBody2D(physicsMass=200, physicsInertia=1, nodeNumber=nCarrierRigidBody,visualization=VObjectRigidBody2D(graphicsData= [graphicsCarrier]))) mCarrierX = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nCarrierRigidBody,coordinate=0)) mCarrierY = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nCarrierRigidBody,coordinate=1)) mCarrierRot = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nCarrierRigidBody,coordinate=2)) mbs.AddObject(CoordinateConstraint(markerNumbers=[mSuspensionRopeAttachmentNodeX,mCarrierX])) mbs.AddObject(CoordinateConstraint(markerNumbers=[mSuspensionRopeAttachmentNodeY,mCarrierY])) mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mCarrierRot])) #mCarrierWheelLoad = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oCarrierRigidBody, localPosition=[ 0.,0.,0.])) #center of mass (for load) #mbs.AddLoad(Force(markerNumber=mConnectionWheelLoad, loadVector=[0,-CarrierMass*gravityFieldConstant,0])) nSegments = 4 #number of contact segments; must be consistent between nodedata and contact element useFriction = True nFactFriction = 1 if useFriction: nFactFriction = 3 initialGapList = [0.1]*(nSegments*nFactFriction) #initial gap of 0.1 cStiffness = 1e7 mContactCarrier=mbs.AddMarker(MarkerBodyRigid(bodyNumber = oCarrierRigidBody)) for i in cable1ObjectList: mContactCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=i, numberOfSegments = nSegments)) nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,numberOfDataCoordinates=nSegments*nFactFriction)) if useFriction: mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mContactCarrier, mContactCable], nodeNumber = nodeDataContactCable, numberOfContactSegments=nSegments, contactStiffness = cStiffness, circleRadius = carrierWheelRadius)) else: mbs.AddObject(ObjectContactCircleCable2D(markerNumbers=[mContactCarrier, mContactCable], nodeNumber = nodeDataContactCable, numberOfContactSegments=nSegments, contactStiffness = cStiffness, circleRadius = carrierWheelRadius, offset = 0)) #Add Bodys cRigidBodyRotList=[] cRigidBodyTranslationXlist=[] cRigidBodyTranslationYlist=[] distanceBetweenBodys=32 numberOfBodys=2 positionOfBody=0 for i in range(numberOfBodys): #Add Body [oRigid,nRigid,oAleSlidingJoint,oSlidingJoint]=AddBodyWithSlidingJoints(mbs,xPositionOfFirstNodes=0,referencePositionOfBodyAlongCable=positionOfBody,gravityFieldConstant=gravityFieldConstant) mbs.SetObjectParameter(oSlidingJoint, 'classicalFormulation', False) #test model computed with new sliding joint formulation #fix rotation of rigid body mRigidBodyRot = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nRigid, coordinate=2)) #add rigid body marker cRigidBodyRot = mbs.AddObject(CoordinateConstraint(markerNumbers=[mGlobalGround,mRigidBodyRot])) cRigidBodyRotList+=[cRigidBodyRot] #add damper for rotation of rigid body mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mGlobalGround,mRigidBodyRot],damping = 1e6, visualization=VObjectConnectorCoordinateSpringDamper(show=False))) positionOfBody+=distanceBetweenBodys #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++# # Assemble multibody system #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++# mbs.Assemble() #exu.Print(mbs) #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++# # Simualtion settings: #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++# simulationSettings = exu.SimulationSettings() #simulationSettings.staticSolver.loadStepGeometric = True #simulationSettings.staticSolver.adaptiveStep = False simulationSettings.staticSolver.numberOfLoadSteps=10 #simulationSettings.staticSolver.loadStepGeometricRange = 100 simulationSettings.staticSolver.newton.relativeTolerance = 1e-7 #with this error tolerance, the adaptive step selection needs 4 steps #simulationSettings.staticSolver.verboseMode=1 #simulationSettings.staticSolver.verboseModeFile=2 simulationSettings.staticSolver.stabilizerODE2term = 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) #start static calculation mbs.SolveStatic(simulationSettings) #++++++++++++++++++++++++++++++++++++++++ #compute error for test suite: ltgCable = mbs.systemData.GetObjectLTGODE2(cable1ObjectList[int(len(cable1ObjectList)/2)]) nc = ltgCable[1] #vertical displacement exu.Print("select cable coordinate", nc) sol = mbs.systemData.GetODE2Coordinates(); uStatic = sol[nc]; #y-displacement of first node of four bar mechanism exu.Print('static solution of cable1 =',uStatic) exudynTestGlobals.testError = uStatic - (-2.1973218891272532) #before 2023-05-01 (new loads jacobian): -2.1973218869310713 #before 2022-03-09 (old ObjectContactFrictionCircleCable2D): -2.197321886974786 2020-03-05(corrected Cable2DshapeMarker): -2.197321886974786 #2019-12-26: 2.1973218859908146 exudynTestGlobals.testResult = uStatic #++++++++++++++++++++++++++++++++++++++++ #store solution for next computation u = mbs.systemData.GetODE2Coordinates() data = mbs.systemData.GetDataCoordinates() #Add drive via ALE: def userLoadDriveAle(mbs, t, load): #if t < 1: return t*50000*5 #else: return 0 mbs.AddLoad(LoadCoordinate(markerNumber = mALE, loadUserFunction=userLoadDriveAle)) mbs.Assemble() #because of new load #++++++++++++++++++++++++++++++++++++++++ #resuse old solution 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) # remove rigid body motion constraints for dynamic analysis for item in cRigidBodyRotList: mbs.SetObjectParameter(item, 'activeConnector', False) mbs.SetObjectParameter(cAleConstraint, 'activeConnector', False) solveDynamic = True if solveDynamic: # time related settings: steps=400 tend=0.8 h=tend/steps SC.visualizationSettings.openGL.lineWidth = 2 simulationSettings.timeIntegration.numberOfSteps = steps simulationSettings.timeIntegration.endTime = tend simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True simulationSettings.timeIntegration.generalizedAlpha.useNewmark = simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.3 simulationSettings.timeIntegration.verboseMode = 1 simulationSettings.displayStatistics = True mbs.SolveDynamic(simulationSettings) if useGraphics: #SC.WaitForRenderEngineStopFlag() exu.StopRenderer() #safely close rendering window! sol = mbs.systemData.GetODE2Coordinates(); uDynamic = sol[nc]; #y-displacement of first node of four bar mechanism exu.Print('dynamic solution of cable1 =',uDynamic) exudynTestGlobals.testError += uDynamic - (-2.2290865056280076) #before 2023-05-01 (loads jacobian): -2.229086503625397 #before 2022-12-25(resolved BUG 1274): -2.229081157258582; before 2022-03-09 (old ObjectContactFrictionCircleCable2D) : (-2.2290811574753953) #2020-03-05(corrected Cable2DshapeMarker): -2.2290811574753953 #2019-12-26: -2.2290811558815617; 2019-12-18: -2.229126333291627 exudynTestGlobals.testResult += uDynamic exu.Print('result of ANCFslidingAndALEjointTest=',exudynTestGlobals.testResult)