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