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