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