coordinateSpringDamperExt.py
You can view and download this file on Github: coordinateSpringDamperExt.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test model for CoordinateSpringDamperExt, which allows to model contact, friction and limit stops
5#
6# Author: Johannes Gerstmayr
7# Date: 2022-01-22
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#import sys
13#sys.path.append('C:/DATA/cpp/EXUDYN_git/main/bin/WorkingRelease') #for exudyn, itemInterface and exudynUtilities
14
15import exudyn as exu
16from exudyn.utilities import *
17
18useGraphics = True #without test
19#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
20#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
21try: #only if called from test suite
22 from modelUnitTests import exudynTestGlobals #for globally storing test results
23 useGraphics = exudynTestGlobals.useGraphics
24except:
25 class ExudynTestGlobals:
26 pass
27 exudynTestGlobals = ExudynTestGlobals()
28#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
29#from exudyn.physics import StribeckFunction, RegularizedFriction
30
31SC = exu.SystemContainer()
32mbs = SC.AddSystem()
33
34
35useFrictionReg = True
36useFrictionBristle = True
37useLimitStops = True
38useGears = True
39
40endTime = 2
41stepSize = 1e-3*2
42
43L=2
44mass = 0.5
45g = 9.81
46stiffness = 1e2
47omega0 = sqrt(stiffness/mass)
48dRel = 0.02*5
49damping = 2 * dRel * omega0
50
51kSticking = 1e4
52dSticking = 0.01*kSticking
53frictionProportionalZone = 1e-3
54expVel = 0.2
55muFriction = 0.3
56fDynamicFriction = muFriction * (mass*g)
57fStaticFrictionOffset = 0.5*fDynamicFriction
58exu.Print('fMu=', fDynamicFriction)
59
60kLimits = 1e4
61dLimits = 0.001*kLimits
62
63fLoad=stiffness*0.5
64
65u0 = 0.1*L #initial displacement
66v0 = 10 #initial velocity
67
68w = 0.05*L #drawing
69
70#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
71gBackground = []
72gBackground += [GraphicsDataOrthoCubePoint([-0.5*(L+w),0,0],[w,0.5*L,w], color=color4darkgrey)]
73gBackground += [GraphicsDataOrthoCubePoint([ 0.5*(L+w),0,0],[w,0.5*L,w], color=color4darkgrey)]
74gBackground += [GraphicsDataOrthoCubePoint([0,-w,0],[L,w,w], color=color4grey)]
75
76objectGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0],
77 visualization=VObjectGround(graphicsData=gBackground)))
78
79
80
81#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
82if useFrictionReg:
83 nGround0=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
84 groundCoordinateMarker0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround0, coordinate = 0))
85
86 nMass0 = mbs.AddNode(Point(referenceCoordinates = [0,0,0], initialCoordinates = [0,0,0],
87 initialVelocities= [v0,0,0]))
88
89 #add mass points and ground object:
90 gCube = GraphicsDataOrthoCubePoint(size=[w,w,w], color=color4steelblue)
91 massPoint0 = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = nMass0,
92 visualization=VObjectMassPoint(graphicsData=[gCube])))
93
94 node0CoordinateMarker0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nMass0, coordinate = 0))
95
96 mbs.AddObject(CoordinateSpringDamperExt(markerNumbers = [groundCoordinateMarker0, node0CoordinateMarker0],
97 stiffness = stiffness, damping = damping,
98 offset = 0, velocityOffset=0,
99 fDynamicFriction=fDynamicFriction, fStaticFrictionOffset=fStaticFrictionOffset,
100 frictionProportionalZone=frictionProportionalZone, exponentialDecayStatic=expVel,
101 #springForceUserFunction = UFspring,
102 visualization=VObjectConnectorCoordinateSpringDamperExt(show=True)))
103
104 #add loads:
105 # mbs.AddLoad(LoadCoordinate(markerNumber = node0CoordinateMarker0, load = fLoad))
106
107 sensPos0 = mbs.AddSensor(SensorNode(nodeNumber=nMass0, storeInternal=True,
108 outputVariableType=exu.OutputVariableType.Displacement))
109 sensVel0 = mbs.AddSensor(SensorNode(nodeNumber=nMass0, storeInternal=True,
110 outputVariableType=exu.OutputVariableType.Velocity))
111
112#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
113if useFrictionBristle:
114 nGround0=mbs.AddNode(NodePointGround(referenceCoordinates = [0,1*1.25*w,0]))
115 groundCoordinateMarker0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround0, coordinate = 0))
116
117 nMass0 = mbs.AddNode(Point(referenceCoordinates = [0,1*1.25*w,0], initialCoordinates = [0,0,0],
118 initialVelocities= [v0,0,0]))
119
120 #add mass points and ground object:
121 gCube = GraphicsDataOrthoCubePoint(size=[w,w,w], color=color4steelblue)
122 massPoint0 = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = nMass0,
123 visualization=VObjectMassPoint(graphicsData=[gCube])))
124
125 node0CoordinateMarker0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nMass0, coordinate = 0))
126
127 nGeneric = mbs.AddNode(NodeGenericData(initialCoordinates=[1,0,0], numberOfDataCoordinates=3))
128 mbs.AddObject(CoordinateSpringDamperExt(markerNumbers = [groundCoordinateMarker0, node0CoordinateMarker0], nodeNumber=nGeneric,
129 stiffness = stiffness, damping = damping,
130 offset = 0, velocityOffset=0,
131 fDynamicFriction=fDynamicFriction, fStaticFrictionOffset=fStaticFrictionOffset,
132 stickingStiffness=kSticking, stickingDamping=dSticking,
133 frictionProportionalZone=0, exponentialDecayStatic=expVel,
134 #springForceUserFunction = UFspring,
135 visualization=VObjectConnectorCoordinateSpringDamperExt(show=True)))
136
137 #add loads:
138 # mbs.AddLoad(LoadCoordinate(markerNumber = node0CoordinateMarker0, load = fLoad))
139
140 sensPos0b = mbs.AddSensor(SensorNode(nodeNumber=nMass0, storeInternal=True,
141 outputVariableType=exu.OutputVariableType.Displacement))
142 sensVel0b = mbs.AddSensor(SensorNode(nodeNumber=nMass0, storeInternal=True,
143 outputVariableType=exu.OutputVariableType.Velocity))
144
145#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
146if useLimitStops:
147 nGround1=mbs.AddNode(NodePointGround(referenceCoordinates = [0,2*1.25*w,0]))
148 groundCoordinateMarker1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround1, coordinate = 0))
149
150 nMass1 = mbs.AddNode(Point(referenceCoordinates = [0,2*1.25*w,0], initialCoordinates = [0,0,0],
151 initialVelocities= [v0,0,0]))
152
153 #add mass points and ground object:
154 gCube = GraphicsDataOrthoCubePoint(size=[w,w,w], color=color4steelblue)
155 massPoint1 = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = nMass1,
156 visualization=VObjectMassPoint(graphicsData=[gCube])))
157
158 node1CoordinateMarker0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nMass1, coordinate = 0))
159
160 nGeneric = mbs.AddNode(NodeGenericData(initialCoordinates=[0,0,0], numberOfDataCoordinates=3))
161 oCSD = mbs.AddObject(CoordinateSpringDamperExt(markerNumbers = [groundCoordinateMarker1, node1CoordinateMarker0], nodeNumber=nGeneric,
162 #stiffness = stiffness, damping = damping,
163 limitStopsUpper=0.5*L-0.5*w, limitStopsLower=-(0.5*L-0.5*w),
164 limitStopsStiffness=kLimits,limitStopsDamping=dLimits,useLimitStops=True,
165 #fDynamicFriction=fDynamicFriction, fStaticFrictionOffset=0,
166 #frictionProportionalZone=frictionProportionalZone,
167 #springForceUserFunction = UFspring,
168 #stickingStiffness=kSticking, stickingDamping=dSticking, #DELETE
169 visualization=VObjectConnectorCoordinateSpringDamperExt(show=False)))
170
171 sensPos1 = mbs.AddSensor(SensorNode(nodeNumber=nMass1, storeInternal=True,
172 outputVariableType=exu.OutputVariableType.Displacement))
173 sensVel1 = mbs.AddSensor(SensorNode(nodeNumber=nMass1, storeInternal=True,
174 outputVariableType=exu.OutputVariableType.Velocity))
175
176#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
177if useGears: #show that also transmission / gear ratio works; test for limit stops needed ...
178
179 gStiffness = 1e5
180 gDamping = 0.01*gStiffness
181 rad0 = 0.5*w
182 rad1 = 1.5*w
183 omega0 = 2*pi
184 omega1 = -omega0*rad0/rad1
185
186 nG0 = mbs.AddNode(Node1D(referenceCoordinates = [0], #\psi_0ref
187 initialCoordinates=[0], #\psi_0ini
188 initialVelocities=[omega0])) #\psi_t0ini
189 nG1 = mbs.AddNode(Node1D(referenceCoordinates = [0], #\psi_0ref
190 initialCoordinates=[0], #\psi_0ini
191 initialVelocities=[omega1])) #\psi_t0ini
192
193 #add mass points and ground object:
194 gRotor0 = [GraphicsDataOrthoCubePoint(size=[0.5*w,0.5*w,w], color=color4grey)]
195 gRotor0 += [GraphicsDataCylinder(pAxis=[0,0,-0.25*w],vAxis=[0,0,0.5*w], radius = rad0, color=color4orange, nTiles=32)]
196 gRotor1 = [GraphicsDataOrthoCubePoint(size=[3*0.5*w,3*0.5*w,w], color=color4grey)]
197 gRotor1 += [GraphicsDataCylinder(pAxis=[0,0,-0.25*w],vAxis=[0,0,0.5*w], radius = rad1, color=color4dodgerblue, nTiles=32)]
198 gear0 = mbs.AddObject(Rotor1D(physicsInertia = 1, nodeNumber = nG0,
199 referencePosition = [-rad0,-4*w,0],
200 visualization=VRotor1D(graphicsData=gRotor0)))
201 gear1 = mbs.AddObject(Rotor1D(physicsInertia = 1, nodeNumber = nG1,
202 referencePosition = [ rad1,-4*w,0],
203 visualization=VRotor1D(graphicsData=gRotor1)))
204
205 mGear0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nG0, coordinate = 0))
206 mGear1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nG1, coordinate = 0))
207
208 #nGeneric = mbs.AddNode(NodeGenericData(initialCoordinates=[0,0,0], numberOfDataCoordinates=3))
209 oCSD = mbs.AddObject(CoordinateSpringDamperExt(markerNumbers = [mGear1, mGear0],
210 stiffness = gStiffness, damping = gDamping,
211 factor0 = 1, factor1 = -rad0/rad1,
212 #limitStopsUpper=0.5*L-0.5*w, limitStopsLower=-(0.5*L-0.5*w),
213 #limitStopsStiffness=kLimits,limitStopsDamping=dLimits,
214 frictionProportionalZone=1e-16,#workaround
215 visualization=VObjectConnectorCoordinateSpringDamperExt(show=False)))
216
217 nGroundG=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
218 mGroundG = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGroundG, coordinate = 0))
219 mbs.AddObject(CoordinateConstraint(markerNumbers=[mGear1, mGroundG], offset = omega1, #velocity
220 velocityLevel = True))
221 mbs.AddLoad(LoadCoordinate(markerNumber = mGear0,
222 load = -10, #breaking torque, transmitted to gear1
223 ))
224
225 sensGearPos0 = mbs.AddSensor(SensorNode(nodeNumber=nG0, storeInternal=True,
226 outputVariableType=exu.OutputVariableType.Coordinates))
227 sensGearPos1 = mbs.AddSensor(SensorNode(nodeNumber=nG1, storeInternal=True,
228 outputVariableType=exu.OutputVariableType.Coordinates))
229 sensGearVel0 = mbs.AddSensor(SensorNode(nodeNumber=nG0, storeInternal=True,
230 outputVariableType=exu.OutputVariableType.Coordinates_t))
231 sensGearVel1 = mbs.AddSensor(SensorNode(nodeNumber=nG1, storeInternal=True,
232 outputVariableType=exu.OutputVariableType.Coordinates_t))
233 sensForce = mbs.AddSensor(SensorObject(objectNumber=oCSD, storeInternal=True,
234 outputVariableType=exu.OutputVariableType.Force))
235
236
237#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
238#print(mbs)
239mbs.Assemble()
240
241#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
242
243simulationSettings = exu.SimulationSettings()
244
245simulationSettings.solutionSettings.writeSolutionToFile = False
246simulationSettings.solutionSettings.solutionWritePeriod = 0.1 #data not used
247simulationSettings.solutionSettings.sensorsWritePeriod = 0.002 #data not used
248#simulationSettings.solutionSettings.solutionInformation = 'Nonlinear oscillations: compare linear / nonlinear case'
249simulationSettings.timeIntegration.verboseMode = 1 #turn off, because of lots of output
250#simulationSettings.timeIntegration.stepInformation = 2+64+128+8
251#simulationSettings.timeIntegration.newton.relativeTolerance = 1e-3 #reduce a little bit to improve convergence
252
253simulationSettings.timeIntegration.numberOfSteps = int(endTime/stepSize)
254simulationSettings.timeIntegration.endTime = endTime
255
256if useGraphics:
257 simulationSettings.timeIntegration.simulateInRealtime = True
258 simulationSettings.timeIntegration.realtimeFactor = 2
259
260SC.visualizationSettings.general.graphicsUpdateInterval = 0.02
261SC.visualizationSettings.window.renderWindowSize=[1200,1024]
262#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
263SC.visualizationSettings.general.autoFitScene = False #otherwise, renderState not accepted for zoom
264
265if useGraphics:
266 exu.StartRenderer()
267 mbs.WaitForUserToContinue()
268
269mbs.SolveDynamic(simulationSettings)
270
271if useGraphics:
272 SC.WaitForRenderEngineStopFlag()
273 exu.StopRenderer() #safely close rendering window!
274
275sol = mbs.systemData.GetODE2Coordinates()
276exudynTestGlobals.testResult = np.sum(abs(sol))
277exu.Print('result of coordinateSpringDamperExt=',exudynTestGlobals.testResult) #17.084935539925155
278
279
280if False:
281
282 mbs.PlotSensor(closeAll = True)
283
284 if useLimitStops:
285 mbs.PlotSensor(sensorNumbers=[sensPos1])
286 mbs.PlotSensor(sensorNumbers=[sensVel1])
287 if useFrictionReg:
288 mbs.PlotSensor(sensorNumbers=[sensPos0])
289 if useFrictionBristle:
290 mbs.PlotSensor(sensorNumbers=[sensPos0b], colorCodeOffset=1, newFigure=False, labels=['bristle'])
291 mbs.PlotSensor(sensorNumbers=[sensVel0])
292 if useFrictionBristle:
293 mbs.PlotSensor(sensorNumbers=[sensVel0b], colorCodeOffset=1, newFigure=False, labels=['bristle'])
294
295 if useGears:
296 mbs.PlotSensor(sensorNumbers=[sensGearPos0, sensGearPos1])
297 mbs.PlotSensor(sensorNumbers=[sensGearVel0, sensGearVel1])
298 mbs.PlotSensor(sensorNumbers=[sensForce])