You can view and download this file on Github: massSpringFrictionInteractive.py
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: Example with 1D spring-mass-damper and friction system
# In renderer window you see a long band, where you need zoom into
# the mass-spring-damper to see the effect
#
# Author: Johannes Gerstmayr
# Date: 2020-01-10
#
# 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.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#import sys
#sys.path.append('C:/DATA/cpp/EXUDYN_git/main/bin/WorkingRelease') #for exudyn, itemInterface and exudynUtilities
import exudyn as exu
from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
import exudyn.graphics as graphics #only import if it does not conflict
from exudyn.interactive import InteractiveDialog
from exudyn.physics import StribeckFunction, RegularizedFriction
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from math import exp, sqrt
from numpy import sign
regVel = 1e-4
expVel = 0.1
#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#draw friction curve vs. velocity
if True:
n= 200
x = np.linspace(-1,1,n)
y = np.zeros(n)
for i in range(n):
y[i] = StribeckFunction(x[i], muDynamic=0.1, muStaticOffset=0.05, expVel=expVel)
plt.close('all')
plt.plot(x,y,'b-')
for i in range(n):
y[i] = RegularizedFriction(x[i], muDynamic=0.1, muStaticOffset=0.05, velStatic=0.05, velDynamic=expVel)
plt.plot(x,y,'r-')
plt.figure()
#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
plt.rcParams.update({'font.size': 14})
SC = exu.SystemContainer()
mbs = SC.AddSystem()
print('EXUDYN version='+exu.__version__)
caseName = "movingBand"
#dRel = 0.05 #relative damping
dRel = 0.002
tt=0.05
#parameters of mass-spring with friction
#self-excitation works with (stick-slip); must have enough initial velocity ..., if started from zero, does not work:
param = [4, 40, 40, 400] #stable oscillation with sufficient excitation; stable until damping of 0.0240
param = [4, 40, 10, 400] #still stable oscillation until relative damping of 0.005
param = [4, 40, 5, 400] #still stable oscillation until relative damping of 0.002
param = [4, 40, 3, 400] #no stable oscillation with relative damping of 0.002
vBand = param[0]
fFriction = param[1] #force in Newton, only depends on direction of velocity
fStaticFrictionOffset = param[2]
stiffness = param[3]
kSticking = 1e4
dSticking = 0.01*kSticking
force = 50 #turned on/off for excitation
L=0.5
mass = 0.5
omega0 = sqrt(stiffness/mass)
damping = 2 * dRel * omega0
u0 = 0 #initial displacement
v0 = 4 #initial velocity
#graphics:
lBand = 200*L
w = L*0.5
z=-tt
gBackground = [graphics.Quad([[-lBand,-w,z],[ L,-w,z],[ L, w,z],[-lBand, w,z]],
color=graphics.color.lightgrey, alternatingColor=graphics.color.grey,
nTiles=200, nTilesY=6)]
#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#node for mass point:
nMass=mbs.AddNode(Point(referenceCoordinates = [L,0,0], initialCoordinates = [u0,0,0],
initialVelocities= [v0,0,0]))
#add mass points and ground object:
gCube = graphics.BrickXYZ(-tt, -tt, -tt, tt, tt, tt, graphics.color.steelblue)
massPoint = mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = nMass,
visualization=VObjectMassPoint(graphicsData=[gCube])))
#marker for constraint / springDamper
nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
groundCoordinateMarker = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
mbs.variables['bandVelocity'] = vBand
mbs.variables['relDamping'] = dRel
mbs.variables['dynamicFriction'] = fFriction
mbs.variables['staticFrictionOffset'] = fFriction
mbs.variables['stiffness'] = stiffness
bandCoordinateMarker = 0
if vBand == 0:
#ground
objectGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0],
visualization=VObjectGround(graphicsData=gBackground)))
bandCoordinateMarker = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
else:
#moving bar:
n0=mbs.AddNode(Point(referenceCoordinates = [0,0,0], initialCoordinates = [0,0,0],
initialVelocities= [vBand,0,0]))
bandCoordinateMarker = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n0, coordinate = 0))
mbs.AddObject(MassPoint(physicsMass = mass, nodeNumber = n0,
visualization=VObjectMassPoint(graphicsData=gBackground)))
#mbs.AddLoad(LoadCoordinate(markerNumber=bandCoordinateMarker, load=1e5))
mbs.variables['oCCband'] = mbs.AddObject(CoordinateConstraint(markerNumbers=[groundCoordinateMarker, bandCoordinateMarker],
velocityLevel=True, offset=vBand,
#offsetUserFunction_t=UFbandVelocity,
visualization=VCoordinateConstraint(show=False)))
nodeCoordinateMarker0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nMass, coordinate = 0))
nodeCoordinateMarker1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nMass, coordinate = 1))
nodeCoordinateMarker2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nMass, coordinate = 2))
#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#Spring-Dampers
nGeneric = mbs.AddNode(NodeGenericData(initialCoordinates=[1,0,0], numberOfDataCoordinates=3))
mbs.variables['oFriction'] = mbs.AddObject(CoordinateSpringDamperExt(markerNumbers = [bandCoordinateMarker, nodeCoordinateMarker0],
nodeNumber=nGeneric,
#stiffness = stiffness, damping = damping, #added to separate CSD, otherwise spring is wrong!
fDynamicFriction=fFriction,
fStaticFrictionOffset=fStaticFrictionOffset,
stickingStiffness=kSticking, stickingDamping=dSticking,
exponentialDecayStatic=expVel,
frictionProportionalZone=regVel,
visualization=VObjectConnectorCoordinateSpringDamper(show=False)))
mbs.variables['oSpring'] = mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundCoordinateMarker, nodeCoordinateMarker0],
stiffness = stiffness, damping = damping,
offset = 0)) #damping added
#transverse springs, not needed:
mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundCoordinateMarker, nodeCoordinateMarker1], stiffness = 4000,
visualization=VObjectConnectorCoordinateSpringDamper(show=False)))
mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundCoordinateMarker, nodeCoordinateMarker2], stiffness = 4000,
visualization=VObjectConnectorCoordinateSpringDamper(show=False)))
mbs.variables['loadForce'] = mbs.AddLoad(LoadCoordinate(markerNumber=nodeCoordinateMarker0, load=0)) #for turning on/off
#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#add sensors:
sensPos = mbs.AddSensor(SensorNode(nodeNumber=nMass, fileName='solution/nonlinearPos.txt',
outputVariableType=exu.OutputVariableType.Displacement))
sensVel = mbs.AddSensor(SensorNode(nodeNumber=nMass, fileName='solution/nonlinearVel.txt',
outputVariableType=exu.OutputVariableType.Velocity))
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
#print(mbs)
mbs.Assemble()
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
simulationSettings = exu.SimulationSettings()
simulationSettings.solutionSettings.solutionWritePeriod = 1e-1
simulationSettings.solutionSettings.solutionInformation = "Mass-spring-damper:"+caseName
simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #SHOULD work with 0.9 as well
#simulationSettings.timeIntegration.simulateInRealtime = True
#simulationSettings.timeIntegration.realtimeFactor = 0.2
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
#data obtained from SC.GetRenderState(); use np.round(d['modelRotation'],4)
# SC.visualizationSettings.openGL.initialModelRotation = [[ 0.33 , 0.0882, -0.9399],
# [ 0.0819, 0.9892, 0.1216],
# [ 0.9404, -0.1171, 0.3192]]
SC.visualizationSettings.openGL.initialZoom = 0.5#0.28
SC.visualizationSettings.openGL.initialCenterPoint = [0.297, 0.000318, 0.0]
SC.visualizationSettings.openGL.initialMaxSceneSize = 0.5
SC.visualizationSettings.general.autoFitScene = False
SC.visualizationSettings.general.textSize = 12
SC.visualizationSettings.general.showSolverInformation = 12
SC.visualizationSettings.general.graphicsUpdateInterval = 0.02
SC.visualizationSettings.window.renderWindowSize=[1200,1000]
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++
SC.visualizationSettings.general.autoFitScene = False #otherwise, renderState not accepted for zoom
exu.StartRenderer()
#+++++++++++++++++++++++++++++++++++++++++++++++++++++
#SC.visualizationSettings.general.autoFitScene = False #otherwise, renderState not accepted for zoom
#time.sleep(0.5) #allow window to adjust view
h = 2e-5*0.5 #step size of solver
deltaT = 0.01 #time period to be simulated between every update
#++++++++++++++++++++++++++++
#define items for dialog
dialogItems = [{'type':'label', 'text':'Friction oscillator', 'grid':(0,0,2), 'options':['L']},
#{'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)]},
{'type':'label', 'text':'band velocity:', 'grid':(1,0)},
{'type':'slider', 'range': (0, 20), 'value':mbs.variables['bandVelocity'], 'steps':401, 'variable':'bandVelocity', 'grid':(1,1)},
{'type':'label', 'text':'dynamic friction force:', 'grid':(2,0)},
{'type':'slider', 'range': (0, 100), 'value':mbs.variables['dynamicFriction'], 'steps':101, 'variable':'dynamicFriction', 'grid':(2,1)},
{'type':'label', 'text':'static friction offset:', 'grid':(3,0)},
{'type':'slider', 'range': (0, 100), 'value':mbs.variables['staticFrictionOffset'], 'steps':101, 'variable':'staticFrictionOffset', 'grid':(3,1)},
{'type':'label', 'text':'stiffness:', 'grid':(4,0)},
{'type':'slider', 'range':(0, 2000), 'value':mbs.variables['stiffness'], 'steps':101, 'variable':'stiffness', 'grid':(4,1)},
{'type':'label', 'text':'relative damping:', 'grid':(5,0)},
{'type':'slider', 'range':(0, 0.4), 'value':mbs.variables['relDamping'], 'steps':201, 'variable':'relDamping', 'grid':(5,1)},
{'type':'label', 'text':'add force:', 'grid':(6,0)},
{'type':'radio', 'textValueList':[('on',1),('off',0)], 'value':0, 'variable':'addForce', 'grid': [(6,1),(6,2)]},
]
#++++++++++++++++++++++++++++++++++++++++
#specify subplots to be shown interactively
plt.close('all')
plots={'fontSize':16,'sizeInches':(12,12),'nPoints':200,
'subplots':(3,1), 'sensors':[[(sensPos,0),(sensPos,1),'time','mass position'],
[(sensVel,0),(sensVel,1),'time','mass velocity'],
[(sensPos,1),(sensVel,1),'position (phase space)','velocity (phase space)']
],
'limitsX':[(),(),()], #omit if time auto-range
#'limitsY':[(-0.05,0.05),()]}
'limitsY':[(),(),()]}
#++++++++++++++++++++++++++++++++++++++++
#setup simulation settings and run interactive dialog:
simulationSettings = exu.SimulationSettings()
simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
simulationSettings.solutionSettings.writeSolutionToFile = False
simulationSettings.solutionSettings.solutionWritePeriod = 0.1 #data not used
simulationSettings.solutionSettings.sensorsWritePeriod = 0.1 #data not used
simulationSettings.solutionSettings.solutionInformation = 'Nonlinear oscillations: compare linear / nonlinear case'
simulationSettings.timeIntegration.verboseMode = 0 #turn off, because of lots of output
simulationSettings.solutionSettings.writeInitialValues = False
#simulationSettings.timeIntegration.stepInformation = 2+64+128+8
# simulationSettings.timeIntegration.newton.relativeTolerance = 1e-3 #reduce a little bit to improve convergence
simulationSettings.timeIntegration.numberOfSteps = int(deltaT/h)
simulationSettings.timeIntegration.endTime = deltaT
simulationSettings.timeIntegration.newton.maxIterations = 8
simulationSettings.timeIntegration.adaptiveStepDecrease = 0.25
#this is an exemplariy simulation function, which adjusts some values for simulation
def SimulationUF(mbs, dialog):
mbs.SetObjectParameter(mbs.variables['oFriction'],'fDynamicFriction',mbs.variables['dynamicFriction'])
mbs.SetObjectParameter(mbs.variables['oFriction'],'fStaticFrictionOffset',mbs.variables['staticFrictionOffset'])
mbs.SetObjectParameter(mbs.variables['oSpring'],'stiffness',mbs.variables['stiffness'])
mbs.SetObjectParameter(mbs.variables['oSpring'],'damping',2*mbs.variables['relDamping']*omega0)
mbs.SetLoadParameter(mbs.variables['loadForce'],'load',force*mbs.variables['addForce'])
if 'oCCband' in mbs.variables:
mbs.SetObjectParameter(mbs.variables['oCCband'],'offset',mbs.variables['bandVelocity'])
SC.visualizationSettings.general.autoFitScene = False #use loaded render state
exu.StartRenderer()
if 'renderState' in exu.sys:
SC.SetRenderState(exu.sys[ 'renderState' ])
dialog = InteractiveDialog(mbs=mbs, simulationSettings=simulationSettings,
simulationFunction=SimulationUF,
title='Interactive window',
dialogItems=dialogItems, period=deltaT, realtimeFactor=10,
runOnStart=True,
plots=plots, fontSize=12)
# #stop solver and close render window
exu.StopRenderer() #safely close rendering window!