reevingSystemSpringsTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  A simple 3D reeving system;
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2022-06-16
  8#
  9# 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.
 10#
 11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 12
 13import exudyn as exu
 14from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 15import exudyn.graphics as graphics #only import if it does not conflict
 16
 17useGraphics = True #without test
 18#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 19#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 20try: #only if called from test suite
 21    from modelUnitTests import exudynTestGlobals #for globally storing test results
 22    useGraphics = exudynTestGlobals.useGraphics
 23except:
 24    class ExudynTestGlobals:
 25        pass
 26    exudynTestGlobals = ExudynTestGlobals()
 27#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 28
 29import numpy as np
 30from math import sin, cos, sqrt,pi
 31
 32SC = exu.SystemContainer()
 33mbs = SC.AddSystem()
 34
 35gGround = graphics.CheckerBoard(point=[8,-8,-2], size=32, nTiles=10)
 36oGround = mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=[gGround])))
 37
 38t = 0.1
 39r0 = 0.5
 40g = [0,-9.81,0]
 41
 42posList = [[0,0,0],
 43           [7,-20,1],
 44           [12,0,2],
 45           ]
 46rList = [r0*1.,4*r0,r0,5.9*r0,r0,r0,r0]
 47
 48# posList = [[0,0,0],
 49#            [5,-20,1],
 50#            [8,0,2],
 51#            [12,-10,0],
 52#            [16,0,0]
 53#            ]
 54# rList = [r0*0.,6*r0,r0,5.9*r0,r0,r0,r0]
 55dirList = [-1,1,-1,1,-1,1,1]
 56
 57
 58n = len(posList)
 59nodeList = []
 60bodyList = []
 61markerList = []
 62sheavesRadii = []
 63sheavesAxes = exu.Vector3DList()
 64Lref = 0
 65pLast = [0,0,0]
 66
 67for i, pos in enumerate(posList):
 68    r = rList[i]
 69
 70    graphicsRoll = graphics.Cylinder(pAxis=[0,0,-0.5*t], vAxis=[0,0,t], nTiles=64, radius=r, color=graphics.color.dodgerblue, alternatingColor=graphics.color.darkgrey)
 71
 72    inertiaRoll = InertiaCylinder(density=1000,length=t,outerRadius=r, axis=2)
 73    inertiaRoll = inertiaRoll.Translated([0,-r,0])
 74
 75    dictR = mbs.CreateRigidBody(referencePosition=pos,
 76                                inertia=inertiaRoll,
 77                                gravity=g,
 78                                graphicsDataList=[graphicsRoll, graphics.Basis(inertiaRoll.COM(), length=0.5)],
 79                                returnDict=True)
 80    nR = dictR['nodeNumber']
 81    bR = dictR['bodyNumber']
 82
 83    mR = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nR))
 84    nodeList += [nR]
 85    bodyList += [bR]
 86    markerList += [mR]
 87
 88    sheavesAxes.Append([0,0,dirList[i]*1])
 89    sheavesRadii += [r]
 90
 91    if i != 0:
 92        Lref += NormL2(np.array(pos)-pLast)
 93    if i > 0 and i < len(posList)-1:
 94        #note that in this test example, Lref is slightly too long, leading to negative spring forces (compression) if not treated nonlinearly with default settings in ReevingSystemSprings
 95        Lref += r*pi #0.8*r*pi would always lead to tension
 96
 97    #mR = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bR, localPosition=[0,-r,0]))
 98    #mbs.AddLoad(Force(markerNumber=mR, loadVector=[0,-inertiaRoll.mass*9.81,0]))
 99
100    pLast = pos
101
102
103markerGround0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=posList[0]))
104markerGroundMid = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=posList[2]))
105markerGroundLast = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=posList[-1]))
106
107#sMarkerR = mbs.AddSensor(SensorMarker(markerNumber=markerR, outputVariableType=exu.OutputVariableType.Position))
108
109#%%++++++++++++++++++++++++++++++++++++++++++++++++
110#add joints:
111oJoint0 = mbs.AddObject(GenericJoint(markerNumbers=[markerGround0, markerList[0]],
112                                    constrainedAxes=[1,1,1,1,1,1],
113                                    visualization=VGenericJoint(axesRadius=0.5*t, axesLength=1.5*t)))
114oJointLast = mbs.AddObject(GenericJoint(markerNumbers=[markerGroundLast, markerList[-1]],
115                                    constrainedAxes=[1,1,1,1,1,1],
116                                    visualization=VGenericJoint(axesRadius=0.5*t, axesLength=1.5*t)))
117
118if len(posList) > 3:
119    mbs.AddObject(GenericJoint(markerNumbers=[markerGroundMid, markerList[2]],
120                                        constrainedAxes=[1,1,1,1,1,1],
121                                        visualization=VGenericJoint(axesRadius=0.5*t, axesLength=1.5*t)))
122
123#%%++++++++++++++++++++++++++++++++++++++++++++++++
124#add reeving system spring
125stiffness = 1e5 #stiffness per length
126damping = 0.5*stiffness #dampiung per length
127oRS=mbs.AddObject(ReevingSystemSprings(markerNumbers=markerList, hasCoordinateMarkers=False,
128                                   stiffnessPerLength=stiffness, dampingPerLength=damping,
129                                   referenceLength = Lref,
130                                   dampingTorsional = 0.01*damping, dampingShear = 0.1*damping,
131                                   sheavesAxes=sheavesAxes, sheavesRadii=sheavesRadii,
132                                   #regularizationForce = -1, #purely linear system
133                                   visualization=VReevingSystemSprings(ropeRadius=0.05, color=graphics.color.dodgerblue)))
134
135#%%++++++++++++++++++++++++++++++++++++++++++++++++
136#add sensors
137if True:
138    sPos1 = mbs.AddSensor(SensorNode(nodeNumber=nodeList[1], storeInternal=True,
139                                          outputVariableType=exu.OutputVariableType.Position))
140    sOmega1 = mbs.AddSensor(SensorNode(nodeNumber=nodeList[1], storeInternal=True,
141                                          outputVariableType=exu.OutputVariableType.AngularVelocity))
142    sLength= mbs.AddSensor(SensorObject(objectNumber=oRS, storeInternal=True,
143                                          outputVariableType=exu.OutputVariableType.Distance))
144    sLength_t= mbs.AddSensor(SensorObject(objectNumber=oRS, storeInternal=True,
145                                          outputVariableType=exu.OutputVariableType.VelocityLocal))
146    sForce= mbs.AddSensor(SensorObject(objectNumber=oRS, storeInternal=True,
147                                          outputVariableType=exu.OutputVariableType.ForceLocal))
148
149#%%++++++++++++++++++++++++++++++++++++++++++++++++
150#simulate:
151mbs.Assemble()
152
153simulationSettings = exu.SimulationSettings() #takes currently set values or default values
154
155tEnd = 2
156if useGraphics:
157    tEnd = 2 #200
158h=0.01  #use small step size to detext contact switching
159
160simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
161simulationSettings.timeIntegration.endTime = tEnd
162#simulationSettings.solutionSettings.writeSolutionToFile= True #set False for CPU performance measurement
163simulationSettings.solutionSettings.sensorsWritePeriod = 0.01
164
165simulationSettings.timeIntegration.verboseMode = 1
166
167simulationSettings.timeIntegration.newton.useModifiedNewton = True
168
169SC.visualizationSettings.nodes.show = True
170SC.visualizationSettings.nodes.drawNodesAsPoint  = False
171SC.visualizationSettings.nodes.showBasis = True
172SC.visualizationSettings.nodes.basisSize = 0.2
173
174SC.visualizationSettings.openGL.multiSampling = 4
175
176#SC.visualizationSettings.general.autoFitScene = False #use loaded render state
177# useGraphics = True
178if useGraphics:
179    exu.StartRenderer()
180    if 'renderState' in exu.sys:
181        SC.SetRenderState(exu.sys[ 'renderState' ])
182    mbs.WaitForUserToContinue()
183
184
185mbs.SolveDynamic(simulationSettings,
186                 #solverType=exu.DynamicSolverType.TrapezoidalIndex2 #in this case, drift shows up significantly!
187                 )
188
189if useGraphics:
190    SC.WaitForRenderEngineStopFlag()
191    exu.StopRenderer() #safely close rendering window!
192
193if True:
194
195
196    mbs.PlotSensor(sPos1, components=[0,1,2], labels=['pos X','pos Y','pos Z'], closeAll=True)
197    mbs.PlotSensor(sOmega1, components=[0,1,2], labels=['omega X','omega Y','omega Z'])
198    mbs.PlotSensor(sLength, components=[0], labels=['length'])
199    mbs.PlotSensor(sLength_t, components=[0], labels=['vel'])
200    mbs.PlotSensor(sForce, components=[0], labels=['force'])
201
202
203
204
205#compute error for test suite:
206sol2 = mbs.systemData.GetODE2Coordinates();
207u = np.linalg.norm(sol2);
208exu.Print('solution of ReevingSystemSprings=',u)
209
210exudynTestGlobals.testResult = u