coordinateSpringDamper.py

You can view and download this file on Github: coordinateSpringDamper.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Example with 1D spring-mass-damper system; 
  5#           Compare viscous damping or friction (useFriction=True); includes exact reference solution
  6#
  7# Author:   Johannes Gerstmayr
  8# Date:     2019-08-15
  9#
 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.
 11#
 12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 13
 14import exudyn as exu
 15#a local run must not open a plot window when windows are suppressed (#2477)
 16if exu.special.userInterface.suppressPlots: import matplotlib; matplotlib.use('Agg')
 17from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 18import exudyn.graphics as graphics
 19
 20SC = exu.SystemContainer()
 21mbs = SC.AddSystem()
 22
 23print('EXUDYN version='+exu.config.Version())
 24useGraphics=True
 25useFriction = False
 26
 27
 28L=0.5
 29mass = 1.6
 30k = 4000
 31omega0 = 50 # sqrt(4000/1.6)
 32dRel = 0.05
 33d = dRel * 2 * 80 #80=sqrt(1.6*4000)
 34u0=-0.08
 35v0=1
 36f = 80
 37x0 = f/k
 38print('static value = '+str(x0))
 39fFriction = 20 #force in Newton, only depends on direction of velocity
 40#print('damping='+str(d))
 41
 42#node for mass point:
 43n1=mbs.AddNode(Point(referenceCoordinates = [L,0,0], initialCoordinates = [u0,0,0], initialVelocities= [v0,0,0]))
 44nGround=mbs.AddNode(NodePointGround(referenceCoordinates = [L,0,0]))
 45
 46#add mass points and ground object:
 47objectGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0]))
 48massPoint = mbs.AddObject(MassPoint(mass = mass, nodeNumber = n1))
 49
 50#marker for constraint / springDamper
 51
 52#groundMarker = mbs.AddMarker(MarkerBodyPosition(bodyNumber = objectGround, localPosition= [0, 0, 0]))
 53#bodyMarker = mbs.AddMarker(MarkerBodyPosition(bodyNumber = massPoint, localPosition= [0, 0, 0]))
 54
 55groundCoordinateMarker = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= nGround, coordinate = 0))
 56nodeCoordinateMarker0  = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 0))
 57nodeCoordinateMarker1  = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 1))
 58nodeCoordinateMarker2  = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= n1, coordinate = 2))
 59
 60#Spring-Dampers
 61mbs.AddObject(CoordinateSpringDamperExt(markerNumbers = [groundCoordinateMarker, nodeCoordinateMarker0], 
 62                                     stiffness = k, 
 63                                     damping = d*(1-useFriction), 
 64                                     dynamicFrictionForce=0.2*fFriction*useFriction,
 65                                     frictionProportionalZone=0.01)) #offset must be zero, because coordinates just represent the displacements
 66mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundCoordinateMarker, nodeCoordinateMarker1], stiffness = k)) 
 67mbs.AddObject(CoordinateSpringDamper(markerNumbers = [groundCoordinateMarker, nodeCoordinateMarker2], stiffness = k)) 
 68
 69
 70#add loads:
 71mbs.AddLoad(LoadCoordinate(markerNumber = nodeCoordinateMarker0, load = f))
 72
 73print(mbs)
 74mbs.Assemble()
 75
 76simulationSettings = exu.SimulationSettings()
 77tEnd = 1 #1
 78steps = 1000    #1000
 79simulationSettings.solution.file.writePeriod = 1e-3
 80simulationSettings.timeIntegration.numberOfSteps = steps
 81simulationSettings.timeIntegration.endTime = tEnd
 82
 83simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1 #SHOULD work with 0.9 as well
 84
 85if useGraphics: 
 86    SC.renderer.Start()
 87
 88mbs.SolveDynamic(simulationSettings)
 89
 90if useGraphics: 
 91    SC.renderer.DoIdleTasks()
 92    SC.renderer.Stop() #safely close rendering window!
 93
 94
 95u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
 96uCartesianSpringDamper= u[0] - L
 97errorCartesianSpringDamper = uCartesianSpringDamper - 0.011834933407066539 #for 1000 steps, endtime=1; this is different from CartesianSpringDamper because of offset L (rounding errors around 1e-14)
 98print('error cartesianSpringDamper=',errorCartesianSpringDamper)
 99
100#return abs(errorCartesianSpringDamper)
101
102#+++++++++++++++++++++++++++++++++++++++++++++++++++++
103#exact solution:
104import numpy as np
105import matplotlib.pyplot as plt
106import matplotlib.ticker as ticker
107
108omega0 = np.sqrt(k/mass)     #recompute for safety
109dRel = d/(2*np.sqrt(k*mass)) #recompute for safety
110omega = omega0*np.sqrt(1-dRel**2)
111C1 = u0-x0 #static solution needs to be considered!
112C2 = (v0+omega0*dRel*C1) / omega #C1 used instead of classical solution with u0, because x0 != 0 !!!
113
114refSol = np.zeros((steps+1,2))
115for i in range(0,steps+1):
116    t = tEnd*i/steps
117    refSol[i,0] = t
118    refSol[i,1] = np.exp(-omega0*dRel*t)*(C1*np.cos(omega*t) + C2*np.sin(omega*t))+x0
119
120print('refSol=',refSol[steps,1])
121print('error exact-numerical=',refSol[steps,1] - uCartesianSpringDamper)
122
123data = np.loadtxt(OutputFilePath('solution/coordinatesSolution.txt'), comments='#', delimiter=',')
124plt.plot(data[:,0], data[:,1], 'b-') #numerical solution
125plt.plot(refSol[:,0], refSol[:,1], 'r-') #exact solution
126
127ax=plt.gca() # get current axes
128ax.grid(True, 'major', 'both')
129ax.xaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
130ax.yaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
131plt.tight_layout()
132plt.show()