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 * #includes itemInterface and rigidBodyUtilities
 20import exudyn.graphics as graphics
 21
 22import numpy as np
 23
 24testIsActive = exu.sys.get('testIsActive', False)
 25
 26SC = exu.SystemContainer()
 27mbs = SC.AddSystem()
 28
 29phi0 = 5./180.*np.pi#initial nick angle of disc, 1 degree
 30g = [0,0,-9.81]     #gravity in m/s^2
 31m = 1               #mass in kg
 32r = 0.01            #radius of disc in m
 33w = 0.001           #width of disc in m, just for drawing
 34p0 = [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]
 35initialRotation = RotationMatrixY(phi0)
 36
 37omega0 = [40,0,0*1800/180*np.pi]                   #initial angular velocity around z-axis
 38v0 = Skew(omega0) @ initialRotation @ [0,0,r]   #initial angular velocity of center point
 39#v0 = [0,0,0]                                   #initial translational velocity
 40#exu.Print("v0=",v0)#," = ", [0,10*np.pi*r*np.sin(phi0),0])
 41
 42#inertia for infinitely small ring:
 43inertiaRing = RigidBodyInertia(mass=1, inertiaTensor= np.diag([0.5*m*r**2, 0.25*m*r**2, 0.25*m*r**2]))
 44#print(inertiaRing)
 45
 46#additional graphics for visualization of rotation:
 47graphicsBody = graphics.Brick(centerPoint=[0,0,0],size=[w*1.1,0.7*r,0.7*r], color=graphics.color.lightred)
 48
 49dict0 = mbs.CreateRigidBody(referencePosition=p0,  
 50                            referenceRotationMatrix=initialRotation,  
 51                            initialVelocity=v0,  
 52                            initialAngularVelocity=omega0,  
 53                            inertia=inertiaRing,  
 54                            gravity=g,  
 55                            graphicsDataList=[graphicsBody],  
 56                            returnDict=True)  
 57[n0, b0] = [dict0['nodeNumber'], dict0['bodyNumber']]
 58
 59#ground body and marker
 60gGround = graphics.Brick(centerPoint=[0,0,-0.001],size=[0.3,0.3,0.002], color=graphics.color.lightgrey)
 61oGround = mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=[gGround])))
 62markerGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,0]))
 63
 64#markers for rigid body:
 65markerBody0J0 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=b0, localPosition=[0,0,0]))
 66
 67#rolling disc:
 68nGeneric = mbs.AddNode(NodeGenericData(initialCoordinates=[0,0,0], numberOfDataCoordinates=3))
 69oRolling=mbs.AddObject(ObjectConnectorRollingDiscPenalty(markerNumbers=[markerGround, markerBody0J0], nodeNumber = nGeneric,
 70                                              discRadius=r, dryFriction=[0.8,0.8], dryFrictionProportionalZone=1e-2, 
 71                                              rollingViscousFriction=0.2,
 72                                              contactStiffness=1e5, contactDamping=1e4,
 73                                              visualization=VObjectConnectorRollingDiscPenalty(discWidth=w, color=graphics.color.blue)))
 74
 75
 76#sensor for trace of contact point:
 77if not testIsActive:
 78    sTrail=mbs.AddSensor(SensorObject(objectNumber=oRolling, storeInternal=True,#fileName='solution/rollingDiscTrail.txt', 
 79                               outputVariableType = exu.OutputVariableType.Position))
 80    
 81    sTrailVel=mbs.AddSensor(SensorObject(objectNumber=oRolling, storeInternal=True,#fileName='solution/rollingDiscTrailVel.txt', 
 82                               outputVariableType = exu.OutputVariableType.Velocity))
 83    
 84
 85    sAngVel=mbs.AddSensor(SensorBody(bodyNumber=b0, storeInternal=True,#fileName='solution/rollingDiscAngVel.txt', 
 86                               outputVariableType = exu.OutputVariableType.AngularVelocity))
 87    
 88    sPos=mbs.AddSensor(SensorBody(bodyNumber=b0, storeInternal=True,#fileName='solution/rollingDiscPos.txt', 
 89                               outputVariableType = exu.OutputVariableType.Position))
 90    
 91    sForce=mbs.AddSensor(SensorObject(objectNumber=oRolling, storeInternal=True,#fileName='solution/rollingDiscForceLocal.txt', 
 92                               outputVariableType = exu.OutputVariableType.ForceLocal))
 93
 94mbs.Assemble()
 95
 96simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 97simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
 98
 99tEnd = 0.5
100if not testIsActive:
101    tEnd = 0.5
102
103h=0.0001 
104
105simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
106simulationSettings.timeIntegration.endTime = tEnd
107#simulationSettings.solution.file.writePeriod = 0.01
108simulationSettings.solution.sensors.writePeriod = 0.0005
109#simulationSettings.timeIntegration.verboseMode = 1
110simulationSettings.solution.file.write = False
111
112simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
113simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
114simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
115simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations=True
116
117
118SC.visualizationSettings.nodes.show = True
119SC.visualizationSettings.nodes.drawNodesAsPoint  = False
120SC.visualizationSettings.nodes.showBasis = True
121SC.visualizationSettings.nodes.basisSize = 0.015
122
123if not testIsActive:
124    SC.renderer.Start()
125    SC.renderer.DoIdleTasks()
126
127mbs.SolveDynamic(simulationSettings)
128
129p0=mbs.GetObjectOutput(oRolling, exu.OutputVariableType.Position)
130exu.Print('solution of rollingCoinPenaltyTest=',p0[0]) #use x-coordinate
131
132exu.sys['testResult'] = p0[0]
133
134
135if not testIsActive:
136    SC.renderer.DoIdleTasks()
137    SC.renderer.Stop() #safely close rendering window!
138
139    ##++++++++++++++++++++++++++++++++++++++++++++++q+++++++
140    #plot results
141    if True:
142        
143        
144        mbs.PlotSensor(sTrail, componentsX=[0],components=[1], closeAll=True, title='wheel trail')
145
146        mbs.PlotSensor(sPos, components=[0,1,2], title='wheel position')
147        mbs.PlotSensor(sForce, components=[0,1,2], title='wheel force')
148
149        mbs.PlotSensor(sAngVel, components=[0], title='wheel local angular velocity')