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