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