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