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])