massSpringFrictionInteractive.py
You can view and download this file on Github: massSpringFrictionInteractive.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Example with 1D spring-mass-damper and friction system
5# In renderer window you see a long band, where you need zoom into
6# the mass-spring-damper to see the effect
7#
8# Author: Johannes Gerstmayr
9# Date: 2020-01-10
10#
11# 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.
12#
13#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
14#import sys
15#sys.path.append('C:/DATA/cpp/EXUDYN_git/main/bin/WorkingRelease') #for exudyn, itemInterface and exudynUtilities
16
17import exudyn as exu
18from exudyn.utilities import *
19from exudyn.interactive import InteractiveDialog
20from exudyn.physics import StribeckFunction, RegularizedFriction
21
22import numpy as np
23import matplotlib.pyplot as plt
24import matplotlib.ticker as ticker
25
26from math import exp, sqrt
27from numpy import sign
28
29regVel = 1e-4
30expVel = 0.1
31
32#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33#draw friction curve vs. velocity
34if True:
35 n= 200
36 x = np.linspace(-1,1,n)
37 y = np.zeros(n)
38 for i in range(n):
39 y[i] = StribeckFunction(x[i], muDynamic=0.1, muStaticOffset=0.05, expVel=expVel)
40
41 plt.close('all')
42 plt.plot(x,y,'b-')
43
44 for i in range(n):
45 y[i] = RegularizedFriction(x[i], muDynamic=0.1, muStaticOffset=0.05, velStatic=0.05, velDynamic=expVel)
46
47 plt.plot(x,y,'r-')
48
49 plt.figure()
50
51#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
52plt.rcParams.update({'font.size': 14})
53
54SC = exu.SystemContainer()
55mbs = SC.AddSystem()
56
57print('EXUDYN version='+exu.__version__)
58
59caseName = "movingBand"
60#dRel = 0.05 #relative damping
61
62dRel = 0.002
63tt=0.05
64
65#parameters of mass-spring with friction
66#self-excitation works with (stick-slip); must have enough initial velocity ..., if started from zero, does not work:
67param = [4, 40, 40, 400] #stable oscillation with sufficient excitation; stable until damping of 0.0240
68param = [4, 40, 10, 400] #still stable oscillation until relative damping of 0.005
69param = [4, 40, 5, 400] #still stable oscillation until relative damping of 0.002
70param = [4, 40, 3, 400] #no stable oscillation with relative damping of 0.002
71
72vBand = param[0]
73fFriction = param[1] #force in Newton, only depends on direction of velocity
74fStaticFrictionOffset = param[2]
75stiffness = param[3]
76kSticking = 1e4
77dSticking = 0.01*kSticking
78
79force = 50 #turned on/off for excitation
80L=0.5
81mass = 0.5
82omega0 = sqrt(stiffness/mass)
83damping = 2 * dRel * omega0
84u0 = 0 #initial displacement
85v0 = 4 #initial velocity
86
87#graphics:
88lBand = 200*L
89w = L*0.5
90z=-tt
91gBackground = [GraphicsDataQuad([[-lBand,-w,z],[ L,-w,z],[ L, w,z],[-lBand, w,z]],
92 color=color4lightgrey, alternatingColor=color4grey,
93 nTiles=200, nTilesY=6)]
94
95#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
96
97#node for mass point:
98nMass=mbs.AddNode(Point(referenceCoordinates = [L,0,0], initialCoordinates = [u0,0,0],
99 initialVelocities= [v0,0,0]))
100
101#add mass points and ground object:
102gCube = GraphicsDataOrthoCube(-tt, -tt, -tt, tt, tt, tt, color4steelblue)
103massPoint = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = nMass,
104 visualization=VObjectMassPoint(graphicsData=[gCube])))
105
106#marker for constraint / springDamper
107nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
108groundCoordinateMarker = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
109
110mbs.variables['bandVelocity'] = vBand
111mbs.variables['relDamping'] = dRel
112mbs.variables['dynamicFriction'] = fFriction
113mbs.variables['staticFrictionOffset'] = fFriction
114mbs.variables['stiffness'] = stiffness
115
116bandCoordinateMarker = 0
117if vBand == 0:
118 #ground
119 objectGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0],
120 visualization=VObjectGround(graphicsData=gBackground)))
121 bandCoordinateMarker = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
122else:
123 #moving bar:
124 n0=mbs.AddNode(Point(referenceCoordinates = [0,0,0], initialCoordinates = [0,0,0],
125 initialVelocities= [vBand,0,0]))
126 bandCoordinateMarker = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n0, coordinate = 0))
127 mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = n0,
128 visualization=VObjectMassPoint(graphicsData=gBackground)))
129 #mbs.AddLoad(LoadCoordinate(markerNumber=bandCoordinateMarker, load=1e5))
130 mbs.variables['oCCband'] = mbs.AddObject(CoordinateConstraint(markerNumbers=[groundCoordinateMarker, bandCoordinateMarker],
131 velocityLevel=True, offset=vBand,
132 #offsetUserFunction_t=UFbandVelocity,
133 visualization=VCoordinateConstraint(show=False)))
134
135nodeCoordinateMarker0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nMass, coordinate = 0))
136nodeCoordinateMarker1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nMass, coordinate = 1))
137nodeCoordinateMarker2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nMass, coordinate = 2))
138
139#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
140#Spring-Dampers
141
142nGeneric = mbs.AddNode(NodeGenericData(initialCoordinates=[1,0,0], numberOfDataCoordinates=3))
143mbs.variables['oFriction'] = mbs.AddObject(CoordinateSpringDamperExt(markerNumbers = [bandCoordinateMarker, nodeCoordinateMarker0],
144 nodeNumber=nGeneric,
145 #stiffness = stiffness, damping = damping, #added to separate CSD, otherwise spring is wrong!
146 fDynamicFriction=fFriction,
147 fStaticFrictionOffset=fStaticFrictionOffset,
148 stickingStiffness=kSticking, stickingDamping=dSticking,
149 exponentialDecayStatic=expVel,
150 frictionProportionalZone=regVel,
151 visualization=VObjectConnectorCoordinateSpringDamper(show=False)))
152
153mbs.variables['oSpring'] = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundCoordinateMarker, nodeCoordinateMarker0],
154 stiffness = stiffness, damping = damping,
155 offset = 0)) #damping added
156
157#transverse springs, not needed:
158mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundCoordinateMarker, nodeCoordinateMarker1], stiffness = 4000,
159 visualization=VObjectConnectorCoordinateSpringDamper(show=False)))
160mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundCoordinateMarker, nodeCoordinateMarker2], stiffness = 4000,
161 visualization=VObjectConnectorCoordinateSpringDamper(show=False)))
162
163mbs.variables['loadForce'] = mbs.AddLoad(LoadCoordinate(markerNumber=nodeCoordinateMarker0, load=0)) #for turning on/off
164
165#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
166#add sensors:
167sensPos = mbs.AddSensor(SensorNode(nodeNumber=nMass, fileName='solution/nonlinearPos.txt',
168 outputVariableType=exu.OutputVariableType.Displacement))
169sensVel = mbs.AddSensor(SensorNode(nodeNumber=nMass, fileName='solution/nonlinearVel.txt',
170 outputVariableType=exu.OutputVariableType.Velocity))
171
172
173#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
174#print(mbs)
175mbs.Assemble()
176
177#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
178
179simulationSettings = exu.SimulationSettings()
180
181simulationSettings.solutionSettings.solutionWritePeriod = 1e-1
182simulationSettings.solutionSettings.solutionInformation = "Mass-spring-damper:"+caseName
183
184simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #SHOULD work with 0.9 as well
185#simulationSettings.timeIntegration.simulateInRealtime = True
186#simulationSettings.timeIntegration.realtimeFactor = 0.2
187
188#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
189#data obtained from SC.GetRenderState(); use np.round(d['modelRotation'],4)
190# SC.visualizationSettings.openGL.initialModelRotation = [[ 0.33 , 0.0882, -0.9399],
191# [ 0.0819, 0.9892, 0.1216],
192# [ 0.9404, -0.1171, 0.3192]]
193SC.visualizationSettings.openGL.initialZoom = 0.5#0.28
194SC.visualizationSettings.openGL.initialCenterPoint = [0.297, 0.000318, 0.0]
195SC.visualizationSettings.openGL.initialMaxSceneSize = 0.5
196SC.visualizationSettings.general.autoFitScene = False
197SC.visualizationSettings.general.textSize = 12
198SC.visualizationSettings.general.showSolverInformation = 12
199SC.visualizationSettings.general.graphicsUpdateInterval = 0.02
200SC.visualizationSettings.window.renderWindowSize=[1200,1000]
201#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
202SC.visualizationSettings.general.autoFitScene = False #otherwise, renderState not accepted for zoom
203
204exu.StartRenderer()
205
206#+++++++++++++++++++++++++++++++++++++++++++++++++++++
207
208
209
210#SC.visualizationSettings.general.autoFitScene = False #otherwise, renderState not accepted for zoom
211
212#time.sleep(0.5) #allow window to adjust view
213
214h = 2e-5*0.5 #step size of solver
215deltaT = 0.01 #time period to be simulated between every update
216
217#++++++++++++++++++++++++++++
218#define items for dialog
219dialogItems = [{'type':'label', 'text':'Friction oscillator', 'grid':(0,0,2), 'options':['L']},
220 #{'type':'radio', 'textValueList':[('linear',0),('nonlinear (f=k*u+1000*k*u**3+d*v)',1)], 'value':0, 'variable':'mode', 'grid': [(2,0),(2,1)]},
221 {'type':'label', 'text':'band velocity:', 'grid':(1,0)},
222 {'type':'slider', 'range': (0, 20), 'value':mbs.variables['bandVelocity'], 'steps':401, 'variable':'bandVelocity', 'grid':(1,1)},
223 {'type':'label', 'text':'dynamic friction force:', 'grid':(2,0)},
224 {'type':'slider', 'range': (0, 100), 'value':mbs.variables['dynamicFriction'], 'steps':101, 'variable':'dynamicFriction', 'grid':(2,1)},
225 {'type':'label', 'text':'static friction offset:', 'grid':(3,0)},
226 {'type':'slider', 'range': (0, 100), 'value':mbs.variables['staticFrictionOffset'], 'steps':101, 'variable':'staticFrictionOffset', 'grid':(3,1)},
227 {'type':'label', 'text':'stiffness:', 'grid':(4,0)},
228 {'type':'slider', 'range':(0, 2000), 'value':mbs.variables['stiffness'], 'steps':101, 'variable':'stiffness', 'grid':(4,1)},
229 {'type':'label', 'text':'relative damping:', 'grid':(5,0)},
230 {'type':'slider', 'range':(0, 0.4), 'value':mbs.variables['relDamping'], 'steps':201, 'variable':'relDamping', 'grid':(5,1)},
231 {'type':'label', 'text':'add force:', 'grid':(6,0)},
232 {'type':'radio', 'textValueList':[('on',1),('off',0)], 'value':0, 'variable':'addForce', 'grid': [(6,1),(6,2)]},
233 ]
234
235#++++++++++++++++++++++++++++++++++++++++
236#specify subplots to be shown interactively
237plt.close('all')
238
239plots={'fontSize':16,'sizeInches':(12,12),'nPoints':200,
240 'subplots':(3,1), 'sensors':[[(sensPos,0),(sensPos,1),'time','mass position'],
241 [(sensVel,0),(sensVel,1),'time','mass velocity'],
242 [(sensPos,1),(sensVel,1),'position (phase space)','velocity (phase space)']
243 ],
244 'limitsX':[(),(),()], #omit if time auto-range
245 #'limitsY':[(-0.05,0.05),()]}
246 'limitsY':[(),(),()]}
247
248#++++++++++++++++++++++++++++++++++++++++
249#setup simulation settings and run interactive dialog:
250simulationSettings = exu.SimulationSettings()
251simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
252simulationSettings.solutionSettings.writeSolutionToFile = False
253simulationSettings.solutionSettings.solutionWritePeriod = 0.1 #data not used
254simulationSettings.solutionSettings.sensorsWritePeriod = 0.1 #data not used
255simulationSettings.solutionSettings.solutionInformation = 'Nonlinear oscillations: compare linear / nonlinear case'
256simulationSettings.timeIntegration.verboseMode = 0 #turn off, because of lots of output
257simulationSettings.solutionSettings.writeInitialValues = False
258#simulationSettings.timeIntegration.stepInformation = 2+64+128+8
259# simulationSettings.timeIntegration.newton.relativeTolerance = 1e-3 #reduce a little bit to improve convergence
260
261simulationSettings.timeIntegration.numberOfSteps = int(deltaT/h)
262simulationSettings.timeIntegration.endTime = deltaT
263simulationSettings.timeIntegration.newton.maxIterations = 8
264simulationSettings.timeIntegration.adaptiveStepDecrease = 0.25
265
266#this is an exemplariy simulation function, which adjusts some values for simulation
267def SimulationUF(mbs, dialog):
268 mbs.SetObjectParameter(mbs.variables['oFriction'],'fDynamicFriction',mbs.variables['dynamicFriction'])
269 mbs.SetObjectParameter(mbs.variables['oFriction'],'fStaticFrictionOffset',mbs.variables['staticFrictionOffset'])
270 mbs.SetObjectParameter(mbs.variables['oSpring'],'stiffness',mbs.variables['stiffness'])
271 mbs.SetObjectParameter(mbs.variables['oSpring'],'damping',2*mbs.variables['relDamping']*omega0)
272
273 mbs.SetLoadParameter(mbs.variables['loadForce'],'load',force*mbs.variables['addForce'])
274
275 if 'oCCband' in mbs.variables:
276 mbs.SetObjectParameter(mbs.variables['oCCband'],'offset',mbs.variables['bandVelocity'])
277
278SC.visualizationSettings.general.autoFitScene = False #use loaded render state
279exu.StartRenderer()
280if 'renderState' in exu.sys:
281 SC.SetRenderState(exu.sys[ 'renderState' ])
282
283dialog = InteractiveDialog(mbs=mbs, simulationSettings=simulationSettings,
284 simulationFunction=SimulationUF,
285 title='Interactive window',
286 dialogItems=dialogItems, period=deltaT, realtimeFactor=10,
287 plots=plots, fontSize=12)
288
289# #stop solver and close render window
290exu.StopRenderer() #safely close rendering window!