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