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