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()