rollingCoinTest.py
You can view and download this file on Github: rollingCoinTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Rolling coin example;
5# examine example of Rill, Schaeffer, Grundlagen und Methodik der Mehrkörpersimulation, 2010, page 59
6# Note that in comparison to the literature, we use the local x-axis for the local axis of the coin, z is the normal to the plane
7# mass and inertia do not influence the results, as long as mass and inertia of a infinitely small ring are used
8# gravity is set to [0,0,-9.81m/s^2] and the radius is 0.01m
9#
10# Author: Johannes Gerstmayr
11# Date: 2020-06-19
12#
13# 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.
14#
15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
16
17import exudyn as exu
18from exudyn.utilities import *
19
20import numpy as np
21
22useGraphics = True #without test
23#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
24#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
25try: #only if called from test suite
26 from modelUnitTests import exudynTestGlobals #for globally storing test results
27 useGraphics = exudynTestGlobals.useGraphics
28except:
29 class ExudynTestGlobals:
30 pass
31 exudynTestGlobals = ExudynTestGlobals()
32#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
33
34SC = exu.SystemContainer()
35mbs = SC.AddSystem()
36
37phi0 = 1./180.*np.pi#initial nick angle of disc, 1 degree
38g = [0,0,-9.81] #gravity in m/s^2
39m = 1 #mass in kg
40r = 0.01 #radius of disc in m
41w = 0.001 #width of disc in m, just for drawing
42p0 = [r*np.sin(phi0),0,r*np.cos(phi0)] #origin of disc center point at reference, such that initial contact point is at [0,0,0]
43initialRotation = RotationMatrixY(phi0)
44
45omega0 = [0,0,1800/180*np.pi] #initial angular velocity around z-axis
46v0 = Skew(omega0) @ initialRotation @ [0,0,r] #initial angular velocity of center point
47#v0 = [0,0,0] #initial translational velocity
48#print("v0=",v0)#," = ", [0,10*np.pi*r*np.sin(phi0),0])
49
50#inertia for infinitely small ring:
51inertiaRing = RigidBodyInertia(mass=1, inertiaTensor= np.diag([0.5*m*r**2, 0.25*m*r**2, 0.25*m*r**2]))
52#print(inertiaRing)
53
54#additional graphics for visualization of rotation:
55graphicsBody = GraphicsDataOrthoCubePoint(centerPoint=[0,0,0],size=[w*1.1,0.7*r,0.7*r], color=color4lightred)
56
57[n0,b0]=AddRigidBody(mainSys = mbs,
58 inertia = inertiaRing,
59 nodeType = str(exu.NodeType.RotationEulerParameters),
60 position = p0,
61 rotationMatrix = initialRotation, #np.diag([1,1,1]),
62 angularVelocity = omega0,
63 velocity=v0,
64 gravity = g,
65 graphicsDataList = [graphicsBody])
66
67#ground body and marker
68gGround = GraphicsDataOrthoCubePoint(centerPoint=[0,0,-0.001],size=[0.12,0.12,0.002], color=color4lightgrey)
69oGround = mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=[gGround])))
70markerGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,0]))
71
72#markers for rigid body:
73markerBody0J0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=b0, localPosition=[0,0,0]))
74
75#rolling disc:
76oRolling=mbs.AddObject(ObjectJointRollingDisc(markerNumbers=[markerGround, markerBody0J0],
77 constrainedAxes=[1,1,1], discRadius=r,
78 visualization=VObjectJointRollingDisc(discWidth=w,color=color4blue)))
79
80sForce=mbs.AddSensor(SensorObject(objectNumber=oRolling, storeInternal=True,#fileName='solution/rollingDiscTrail.txt',
81 outputVariableType = exu.OutputVariableType.ForceLocal))
82
83
84
85#sensor for trace of contact point:
86if useGraphics:
87 sTrail=mbs.AddSensor(SensorObject(objectNumber=oRolling, storeInternal=True,#fileName='solution/rollingDiscTrail.txt',
88 outputVariableType = exu.OutputVariableType.Position))
89
90 sVel=mbs.AddSensor(SensorObject(objectNumber=oRolling, storeInternal=True,#fileName='solution/rollingDiscTrailVel.txt',
91 outputVariableType = exu.OutputVariableType.Velocity))
92
93
94
95mbs.Assemble()
96
97simulationSettings = exu.SimulationSettings() #takes currently set values or default values
98
99tEnd = 0.5
100if useGraphics:
101 tEnd = 0.5
102
103h=0.0005 #no visual differences for step sizes smaller than 0.0005
104
105simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
106simulationSettings.timeIntegration.endTime = tEnd
107#simulationSettings.solutionSettings.solutionWritePeriod = 0.01
108simulationSettings.solutionSettings.sensorsWritePeriod = 0.0005
109simulationSettings.solutionSettings.writeSolutionToFile = False
110simulationSettings.timeIntegration.verboseMode = 1
111# simulationSettings.displayStatistics = True
112
113simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
114simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
115simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
116simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations=True
117
118
119SC.visualizationSettings.nodes.show = True
120SC.visualizationSettings.nodes.drawNodesAsPoint = False
121SC.visualizationSettings.nodes.showBasis = True
122SC.visualizationSettings.nodes.basisSize = 0.015
123
124if useGraphics:
125 exu.StartRenderer()
126 mbs.WaitForUserToContinue()
127
128mbs.SolveDynamic(simulationSettings)
129
130p0=mbs.GetObjectOutput(oRolling, exu.OutputVariableType.Position)
131force=mbs.GetSensorValues(sForce)
132exu.Print('force in rollingCoinTest=',force) #use x-coordinate
133
134u = p0[0] + 0.1*(force[0]+force[1]+force[2])
135exu.Print('solution of rollingCoinTest=',u) #use x-coordinate
136
137exudynTestGlobals.testError = u - (1.0634381189385853) #2024-04-29: added force #2020-06-20: 0.002004099927340136; 2020-06-19: 0.002004099760845168 #4s looks visually similar to Rill, but not exactly ...
138exudynTestGlobals.testResult = u
139
140
141if useGraphics:
142 SC.WaitForRenderEngineStopFlag()
143 exu.StopRenderer() #safely close rendering window!
144
145 ##++++++++++++++++++++++++++++++++++++++++++++++q+++++++
146 #plot results
147 if True:
148
149
150 mbs.PlotSensor(sTrail, componentsX=[0],components=[1], closeAll=True, title='wheel trail')
151
152
153 # import matplotlib.pyplot as plt
154 # import matplotlib.ticker as ticker
155
156 # if True:
157 # data = np.loadtxt('solution/rollingDiscTrail.txt', comments='#', delimiter=',')
158 # plt.plot(data[:,1], data[:,2], 'r-',label='contact point trail') #x/y coordinates of trail
159 # else:
160 # #show trail velocity computed numerically and from sensor:
161 # data = np.loadtxt('solution/rollingDiscTrail.txt', comments='#', delimiter=',')
162
163 # nData = len(data)
164 # vVec = np.zeros((nData,2))
165 # dt = data[1,0]-data[0,0]
166 # for i in range(nData-1):
167 # vVec[i+1,0:2] = 1/dt*(data[i+1,1:3]-data[i,1:3])
168
169 # plt.plot(data[:,0], vVec[:,0], 'r-',label='contact point vel x')
170 # plt.plot(data[:,0], vVec[:,1], 'k-',label='contact point vel y')
171 # plt.plot(data[:,0], (vVec[:,0]**2+vVec[:,1]**2)**0.5, 'g-',label='|contact point vel|')
172
173 # trailVel = np.loadtxt('solution/rollingDiscTrailVel.txt', comments='#', delimiter=',')
174 # plt.plot(data[:,0], trailVel[:,1], 'r--',label='trail vel x')
175 # plt.plot(data[:,0], trailVel[:,2], 'k--',label='trail vel y')
176 # plt.plot(data[:,0], trailVel[:,3], 'y--',label='trail vel z')
177 # plt.plot(data[:,0], (trailVel[:,1]**2+trailVel[:,2]**2)**0.5, 'b--',label='|trail vel|')
178
179 # ax=plt.gca() # get current axes
180 # ax.grid(True, 'major', 'both')
181 # ax.xaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
182 # ax.yaxis.set_major_locator(ticker.MaxNLocator(10)) #use maximum of 8 ticks on y-axis
183 # plt.tight_layout()
184 # plt.legend()
185 # plt.show()