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