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])