pendulumFriction.py
You can view and download this file on Github: pendulumFriction.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Mathematical pendulum with friction;
5# Remark: uses two friction models: CoordinateSpringDamper and CartesianSpringDamper
6#
7# Author: Johannes Gerstmayr
8# Date: 2019-12-26
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 *
16from exudyn.FEM import *
17
18import numpy as np
19
20useGraphics = True #without test
21#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
22#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
23try: #only if called from test suite
24 from modelUnitTests import exudynTestGlobals #for globally storing test results
25 useGraphics = exudynTestGlobals.useGraphics
26except:
27 class ExudynTestGlobals:
28 pass
29 exudynTestGlobals = ExudynTestGlobals()
30#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31
32SC = exu.SystemContainer()
33mbs = SC.AddSystem()
34
35
36L = 0.8 #length of arm
37mass = 2.5
38g = 9.81
39
40r = 0.05 #just for graphics
41d = r/2
42
43#add ground object and mass point:
44graphicsBackground = GraphicsDataRectangle(-1.2*L,-1.2*L, 1.2*L, 0.2*L, [1,1,1,1]) #for appropriate zoom
45oGround = mbs.AddObject(ObjectGround(referencePosition = [0,0,0],
46 visualization = VObjectGround(graphicsData = [graphicsBackground])))
47
48
49graphicsSphere = GraphicsDataSphere(point=[L/2,0,0], radius=r, color=[1.,0.2,0.2,1], nTiles = 16)
50graphicsSphere2 = GraphicsDataSphere(point=[0,0,0], radius=r, color=color4steelblue, nTiles = 16)
51graphicsLink = GraphicsDataOrthoCube(-L/2,-d/2,-d/2, L/2,d/2, d/2, [0.5,0.5,0.5,0.5])
52
53inertia = InertiaCuboid(density=mass/(L*d*d), sideLengths=[L,d,d])
54exu.Print("mass=",inertia.mass)
55
56nR0 = mbs.AddNode(Rigid2D(referenceCoordinates=[L/2,0,0])) #body goes from [0,0,0] to [L,0,0]
57oR0 = mbs.AddObject(RigidBody2D(nodeNumber=nR0, physicsMass = inertia.mass, physicsInertia=inertia.inertiaTensor[2][2],
58 visualization = VObjectRigidBody2D(graphicsData = [graphicsLink,graphicsSphere])))
59
60#markers:
61mGround0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition = [0,0,0]))
62mR0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oR0, localPosition=[-L/2,0,0]))
63mTip0 = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oR0, localPosition=[L/2,0,0]))
64mNodeR0 = mbs.AddMarker(MarkerNodePosition(nodeNumber=nR0))
65mR0COM = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oR0, localPosition=[0,0,0]))
66
67
68oRJ0 = mbs.AddObject(RevoluteJoint2D(markerNumbers=[mGround0,mR0]))
69#
70mbs.AddLoad(Force(markerNumber = mNodeR0, loadVector = [0, -mass*g, 0]))
71
72zeroZoneFriction = 1e-3 #zero-zone for velocity in friction
73fFriction = 1 #friction force (norm); acts against velocity
74#user function for friction against velocity vector, including zeroZone
75def UserFunctionSpringDamper(mbs, t, itemIndex, u, v, k, d, offset):
76 vNorm = NormL2(v)
77 f=[v[0],v[1],v[2]]
78 if abs(vNorm) < offset[0]:
79 f = ScalarMult(offset[1]/offset[0], f)
80 else:
81 f = ScalarMult(offset[1]/vNorm, f)
82 return f
83
84mbs.AddObject(CartesianSpringDamper(markerNumbers=[mGround0, mTip0],
85 offset=[zeroZoneFriction, fFriction, 0],
86 springForceUserFunction=UserFunctionSpringDamper))
87
88if useGraphics:
89 sRot1 = mbs.AddSensor(SensorBody(bodyNumber = oR0, fileName='solution/pendulumFrictionRotation0.txt',
90 outputVariableType=exu.OutputVariableType.Rotation))
91
92 sRot2 = mbs.AddSensor(SensorMarker(markerNumber = mR0COM, fileName='solution/pendulumFrictionRotation0marker.txt',
93 writeToFile = useGraphics,
94 outputVariableType=exu.OutputVariableType.Rotation))
95
96sPos = mbs.AddSensor(SensorMarker(markerNumber = mR0COM, writeToFile = False,
97 outputVariableType=exu.OutputVariableType.Position))
98
99mbs.Assemble()
100
101simulationSettings = exu.SimulationSettings()
102
103f = 4000
104simulationSettings.timeIntegration.numberOfSteps = int(1*f)
105simulationSettings.timeIntegration.endTime = 0.0001*f
106simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/5000
107simulationSettings.solutionSettings.sensorsWritePeriod = simulationSettings.timeIntegration.endTime/2000
108#simulationSettings.displayComputationTime = True
109simulationSettings.timeIntegration.verboseMode = 1
110
111#simulationSettings.timeIntegration.newton.useModifiedNewton = False
112simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
113simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = simulationSettings.timeIntegration.generalizedAlpha.useNewmark
114simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.60 #0.62 is approx. the limit
115
116simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = True
117simulationSettings.solutionSettings.coordinatesSolutionFileName= "solution/coordinatesSolution.txt"
118simulationSettings.solutionSettings.writeSolutionToFile=False
119
120#simulationSettings.displayStatistics = True
121
122SC.visualizationSettings.nodes.defaultSize = 0.05
123
124if useGraphics:
125 exu.StartRenderer()
126 mbs.WaitForUserToContinue()
127
128mbs.SolveDynamic(simulationSettings)
129
130p0=mbs.GetObjectOutputBody(oR0, exu.OutputVariableType.Position, localPosition=[0,0,0])
131exu.Print("p0=", p0)
132
133p0 = mbs.GetSensorValues(sPos) #obtain values from marker
134exu.Print("p0=", p0, '(marker)')
135u=NormL2(p0)
136exu.Print('solution of pendulumFriction=',u)
137
138exudynTestGlobals.testError = u - (0.3999999877698205) #2020-04-22: 0.3999999877698205
139exudynTestGlobals.testResult = u
140
141
142if useGraphics:
143 SC.WaitForRenderEngineStopFlag()
144 exu.StopRenderer() #safely close rendering window!
145
146
147
148 mbs.PlotSensor([sRot1, sRot2], components=[0,2], closeAll=True, markerStyles=['x','+'])
149
150 # import matplotlib.pyplot as plt
151 # import matplotlib.ticker as ticker
152
153 # data = np.loadtxt('solution/pendulumFrictionRotation0.txt', comments='#', delimiter=',')
154 # plt.plot(data[:,0], data[:,1], 'b-', label='rotation 0') #ccordinate 1 = rotation, scalar for ObjectRigidBody2D
155 # data = np.loadtxt('solution/pendulumFrictionRotation0marker.txt', comments='#', delimiter=',')
156 # plt.plot(data[:,0], data[:,3], 'r-', label='rotation 0') #ccordinate 3 = rotation, Z-coordinate because marker always 3D
157
158 # ax=plt.gca() # get current axes
159 # ax.grid(True, 'major', 'both')
160 # ax.xaxis.set_major_locator(ticker.MaxNLocator(10))
161 # ax.yaxis.set_major_locator(ticker.MaxNLocator(10))
162 # plt.xlabel("time (s)")
163 # plt.ylabel("angle (rad)")
164 # plt.tight_layout() #better arrangement of plot
165 # plt.legend()
166 # plt.show()