rollingCoinPenaltyTest.py
You can view and download this file on Github: rollingCoinPenaltyTest.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# In this example, the penalty formulation is used, which additionally treats friction
10#
11# Author: Johannes Gerstmayr
12# Date: 2020-06-19
13#
14# 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.
15#
16#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
17
18import exudyn as exu
19from exudyn.utilities import *
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 = 5./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)+0.01] #origin of disc center point at reference, such that initial contact point is at [0,0,0]
44initialRotation = RotationMatrixY(phi0)
45
46omega0 = [40,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#exu.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 = GraphicsDataOrthoCubePoint(centerPoint=[0,0,0],size=[w*1.1,0.7*r,0.7*r], color=color4lightred)
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 = GraphicsDataOrthoCubePoint(centerPoint=[0,0,-0.001],size=[0.3,0.3,0.002], color=color4lightgrey)
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:
77nGeneric = mbs.AddNode(NodeGenericData(initialCoordinates=[0,0,0], numberOfDataCoordinates=3))
78oRolling=mbs.AddObject(ObjectConnectorRollingDiscPenalty(markerNumbers=[markerGround, markerBody0J0], nodeNumber = nGeneric,
79 discRadius=r, dryFriction=[0.8,0.8], dryFrictionProportionalZone=1e-2,
80 rollingFrictionViscous=0.2,
81 contactStiffness=1e5, contactDamping=1e4,
82 visualization=VObjectConnectorRollingDiscPenalty(discWidth=w, color=color4blue)))
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 sTrailVel=mbs.AddSensor(SensorObject(objectNumber=oRolling, storeInternal=True,#fileName='solution/rollingDiscTrailVel.txt',
91 outputVariableType = exu.OutputVariableType.Velocity))
92
93
94 sAngVel=mbs.AddSensor(SensorBody(bodyNumber=b0, storeInternal=True,#fileName='solution/rollingDiscAngVel.txt',
95 outputVariableType = exu.OutputVariableType.AngularVelocity))
96
97 sPos=mbs.AddSensor(SensorBody(bodyNumber=b0, storeInternal=True,#fileName='solution/rollingDiscPos.txt',
98 outputVariableType = exu.OutputVariableType.Position))
99
100 sForce=mbs.AddSensor(SensorObject(objectNumber=oRolling, storeInternal=True,#fileName='solution/rollingDiscForceLocal.txt',
101 outputVariableType = exu.OutputVariableType.ForceLocal))
102
103mbs.Assemble()
104
105simulationSettings = exu.SimulationSettings() #takes currently set values or default values
106
107tEnd = 0.5
108if useGraphics:
109 tEnd = 0.5
110
111h=0.0001
112
113simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
114simulationSettings.timeIntegration.endTime = tEnd
115#simulationSettings.solutionSettings.solutionWritePeriod = 0.01
116simulationSettings.solutionSettings.sensorsWritePeriod = 0.0005
117#simulationSettings.timeIntegration.verboseMode = 1
118simulationSettings.solutionSettings.writeSolutionToFile = False
119
120simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
121simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
122simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
123simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations=True
124
125
126SC.visualizationSettings.nodes.show = True
127SC.visualizationSettings.nodes.drawNodesAsPoint = False
128SC.visualizationSettings.nodes.showBasis = True
129SC.visualizationSettings.nodes.basisSize = 0.015
130
131if useGraphics:
132 exu.StartRenderer()
133 mbs.WaitForUserToContinue()
134
135mbs.SolveDynamic(simulationSettings)
136
137p0=mbs.GetObjectOutput(oRolling, exu.OutputVariableType.Position)
138exu.Print('solution of rollingCoinPenaltyTest=',p0[0]) #use x-coordinate
139
140exudynTestGlobals.testError = p0[0] - (0.03489603106769764) #2020-06-20: 0.03489603106769764
141exudynTestGlobals.testResult = p0[0]
142
143
144if useGraphics:
145 SC.WaitForRenderEngineStopFlag()
146 exu.StopRenderer() #safely close rendering window!
147
148 ##++++++++++++++++++++++++++++++++++++++++++++++q+++++++
149 #plot results
150 if True:
151
152
153 mbs.PlotSensor(sTrail, componentsX=[0],components=[1], closeAll=True, title='wheel trail')
154
155 mbs.PlotSensor(sPos, components=[0,1,2], title='wheel position')
156 mbs.PlotSensor(sForce, components=[0,1,2], title='wheel force')
157
158 mbs.PlotSensor(sAngVel, components=[0], title='wheel local angular velocity')