rollingDiscTangentialForces.py
You can view and download this file on Github: rollingDiscTangentialForces.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test for tangential forces of rolling disc
5#
6# Author: Johannes Gerstmayr, Michael Pieber
7# Date: 2024-04-29
8#
9# 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.
10#
11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12
13import exudyn as exu
14from exudyn.utilities import *
15
16import numpy as np
17
18useGraphics = True #without test
19#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
20#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
21try: #only if called from test suite
22 from modelUnitTests import exudynTestGlobals #for globally storing test results
23 useGraphics = exudynTestGlobals.useGraphics
24except:
25 class ExudynTestGlobals:
26 pass
27 exudynTestGlobals = ExudynTestGlobals()
28#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
29SC = exu.SystemContainer()
30mbs = SC.AddSystem()
31
32
33#Dimensions
34mass = 30 #kg
35r = 1.8 #m
36d = 3 #m
37kAB=1.2 #inertia radius
38
39JradAB = mass*kAB**2
40#wheel
41w=0.1 #for drawing
42
43
44
45
46#####################################
47Fn=400 #N
48g=9.81
49
50omegaY = np.sqrt(r*(Fn-mass*g)/JradAB)
51omegaZ = d/r*omegaY #rad/s # = angular velocity omega_y
52
53exu.Print('omegaY=',omegaY)
54exu.Print('omegaZ=',omegaZ)
55
56
57theta = np.array([[mass*r**2,0,0], [0,mass*r**2,0], [0,0,JradAB]])
58wV= np.array([0,omegaY*0,-omegaZ])
59
60
61wL=np.array([0,omegaY,0])
62M=Skew(wL)@theta@wV
63
64exu.Print('necessary torque=',d*(Fn-mass*g), ', gyro torque=', -M[0])
65exu.Print(d*(mass*g-Fn))
66#centrifugal force in simulation: -396.37
67exu.Print('centrifugal force:',d * omegaY**2 * mass)
68
69#%%
70#####################################
71#####################################
72omegaZ=omegaY*1
73omegaX=-d/r*omegaZ #rad/s angular velocity of the wheel
74g = [0,0,-9.81] #gravity
75bodyDim = [d,0.1,0.1] #body dimensions
76
77
78v0 = omegaZ * d
79v0Vec = np.array([0,v0,0])
80p0 = [0,0,0] #origin of pendulum
81# p1 = [bodyDim[0],0,0] #origin of wheel
82
83pMid0 = [d*0.5,0,r] #center of mass, body0
84pMid1 = [d,0,r] #center of mass, body1
85
86theta = np.array([[JradAB,0,0], [0,mass*r**2,0], [0,0,mass*r**2]])
87iWheel = RigidBodyInertia(mass,theta,p0)
88# exu.Print(iWheel)
89iCube = InertiaCuboid(1, sideLengths=[0.1,0.1,0.1])
90# exu.Print(iCube)
91
92
93w1=np.array([-omegaZ,0,omegaX]) #wheel + frame
94w2=np.array([0,0,omegaX]) #frame
95
96M=Skew(w2)@theta@w1 #possible for symmetric rotor, but dangerous!
97
98Mgyro = w2[2]*theta[0,0]*w1[0] #rotor solution; condition: rotor is symmetric
99
100exu.Print('omega1 x (Theta * omega2) = ',M[1])
101exu.Print('Mgyro = ',Mgyro)
102exu.Print('omega1 x (Theta * omega1) = ',Skew(w1)@theta@w1)
103exu.Print(d*(-Fn+mass*9.81))
104
105
106###############################################################################
107#
108
109
110#initial acceleration:
111angAcc = np.array([0., -14.28472222, 0. ])
112exu.Print('Theta * angAcc = ',theta@angAcc)
113exu.Print('Theta * angAcc + omega1 x (Theta * omega1) = ',theta@angAcc + Skew(w1)@theta@w1)
114
115planeNormal = RotationMatrixX(0.*pi)@np.array([0,0,1])
116vOff = -r*planeNormal+[0,0,r]
117
118#graphics for body
119graphicsBody = GraphicsDataRigidLink(p0=[-0.5*bodyDim[0],0,0], p1=[0.5*bodyDim[0],0,0],
120 axis1=[0,0,1], radius=[0.01,0.01],
121 thickness = 0.2, width = [0.2,0.2], color=color4lightred)
122
123#%%
124[n0,b0]=AddRigidBody(mainSys = mbs,
125 inertia = iCube,
126 nodeType = str(exu.NodeType.RotationEulerParameters),
127 position = pMid0,
128 rotationMatrix = np.diag([1,1,1]),
129 velocity = list(0.5*v0Vec),
130 angularVelocity = [0,0,omegaZ],
131 #gravity = g,
132 graphicsDataList = [graphicsBody])
133
134
135graphicsBodyWheel = GraphicsDataOrthoCubePoint(centerPoint=[0,0,0],size=[w*2,1.4*r,1.4*r], color=color4lightred)
136[n1,bWheel]=AddRigidBody(mainSys = mbs,
137 inertia = iWheel,
138 nodeType = str(exu.NodeType.RotationEulerParameters),
139 position = pMid1,
140 velocity = list(v0Vec),
141 rotationMatrix = np.diag([1,1,1]),
142 angularVelocity = [omegaX,0,omegaZ],
143 gravity = g,
144 graphicsDataList = [graphicsBodyWheel])
145
146
147
148#ground body and marker
149#graphicsPlane = GraphicsDataOrthoCubePoint(centerPoint=[0,0,-0.1],size=[3*d,3*d,0.2], color=color4grey)
150graphicsPlane = GraphicsDataCheckerBoard(point=vOff, normal=planeNormal, size=3*d)
151
152oGround = mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=[graphicsPlane])))
153markerSupportGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,r]))
154
155#markers for rigid body:
156markerBody0J0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=b0, localPosition=[-0.5*bodyDim[0],0,0]))
157markerBody0J1 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=b0, localPosition=[0.5*bodyDim[0],0,0]))
158
159markerBodyWheel = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bWheel, localPosition=[0,0,0]))
160
161mbs.AddObject(SphericalJoint(markerNumbers=[markerSupportGround,markerBody0J0],
162 constrainedAxes=[1,1,1],
163 visualization=VObjectJointSpherical(jointRadius=0.025)))
164
165mbs.AddObject(GenericJoint(markerNumbers=[markerBody0J1, markerBodyWheel],
166 constrainedAxes=[0*1,1,1,0,1,1],
167 visualization=VObjectJointGeneric(axesRadius=0.01, axesLength=0.1)))
168
169# rolling disc joint:
170markerRollingPlane = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=vOff))
171#oRolling=mbs.AddObject(ObjectJointRollingDisc(markerNumbers=[markerRollingPlane,markerBody0J1],
172
173if True:
174 oRolling=mbs.AddObject(ObjectJointRollingDisc(markerNumbers=[markerRollingPlane,markerBodyWheel],
175 constrainedAxes=[1,1,1], #note that tangential constraints lead to additional forces on ground ==> changes force on ground!
176 discRadius=r,
177 #planeNormal = planeNormal,
178 visualization=VObjectJointRollingDisc(discWidth=w,color=color4blue)))
179else:
180 nGeneric = mbs.AddNode(NodeGenericData(initialCoordinates=[0,0,0], numberOfDataCoordinates=3))
181 oRolling=mbs.AddObject(ObjectConnectorRollingDiscPenalty(markerNumbers=[markerRollingPlane,markerBodyWheel],
182 nodeNumber=nGeneric,
183 discRadius=r, discAxis=[1,0,0],
184 planeNormal = planeNormal,
185 contactStiffness=1e5, contactDamping=1e3,
186 dryFriction=[1,1],
187 visualization=VObjectConnectorRollingDiscPenalty(discWidth=w,color=color4blue)))
188
189sForce = mbs.AddSensor(SensorObject(objectNumber=oRolling, storeInternal=True,
190 outputVariableType = exu.OutputVariableType.ForceLocal))
191
192sTrailVel = mbs.AddSensor(SensorObject(objectNumber=oRolling, storeInternal=True,
193 outputVariableType = exu.OutputVariableType.Velocity))
194
195sAngVel = mbs.AddSensor(SensorBody(bodyNumber=bWheel, storeInternal=True,
196 outputVariableType = exu.OutputVariableType.AngularVelocity))
197sAngVelLocal = mbs.AddSensor(SensorBody(bodyNumber=bWheel, storeInternal=True,
198 outputVariableType = exu.OutputVariableType.AngularVelocityLocal))
199sAngAcc = mbs.AddSensor(SensorBody(bodyNumber=bWheel, storeInternal=True,
200 outputVariableType = exu.OutputVariableType.AngularAcceleration))
201
202mbs.Assemble()
203
204
205simulationSettings = exu.SimulationSettings() #takes currently set values or default values
206useGraphics=False
207tEnd = 0.1
208if useGraphics:
209 tEnd = 2
210
211h = 0.001
212simulationSettings.timeIntegration.endTime = tEnd #0.2 for testing
213simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
214#simulationSettings.solutionSettings.solutionWritePeriod = 0.01
215#simulationSettings.solutionSettings.sensorsWritePeriod = 0.01
216simulationSettings.timeIntegration.verboseMode = 1
217
218simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations=True
219
220simulationSettings.timeIntegration.simulateInRealtime = True
221
222SC.visualizationSettings.connectors.showJointAxes = True
223SC.visualizationSettings.connectors.jointAxesLength = 0.3
224SC.visualizationSettings.connectors.jointAxesRadius = 0.08
225SC.visualizationSettings.openGL.lineWidth=2 #maximum
226SC.visualizationSettings.openGL.shadow=0.2
227SC.visualizationSettings.openGL.multiSampling = 4
228
229if useGraphics:
230 exu.StartRenderer()
231 if 'renderState' in exu.sys:
232 SC.SetRenderState(exu.sys['renderState'])
233 mbs.WaitForUserToContinue()
234
235exu.SolveDynamic(mbs, simulationSettings)
236
237if useGraphics:
238 SC.WaitForRenderEngineStopFlag()
239 exu.StopRenderer() #safely close rendering window!
240
241force = mbs.GetSensorValues(sForce)
242u = 1e-3*(abs(force[0]) + abs(force[1]) + abs(force[2]))
243exu.Print('rollingDiscTangentialForces: F=',force) #use x-coordinate
244exu.Print('solution of rollingDiscTangentialForces=',u) #use x-coordinate
245
246exudynTestGlobals.testError = u - (1.0342017388721547)
247exudynTestGlobals.testResult = u
248
249
250if True:
251 from exudyn.plot import PlotSensor
252 PlotSensor(mbs,closeAll=True)
253 PlotSensor(mbs, sensorNumbers=[sForce], components=[0,1,2])
254 # PlotSensor(mbs, sensorNumbers=sAngVel, components=[0,1,2])
255 # PlotSensor(mbs, sensorNumbers=sAngVelLocal, components=[0,1,2])
256 # PlotSensor(mbs, sensorNumbers=sAngAcc, components=[0,1,2])
257
258 # PlotSensor(mbs, sensorNumbers=sTrailVel, components=[0,1])