createRollingDiscPenaltyTest.py

You can view and download this file on Github: createRollingDiscPenaltyTest.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test model for CreateRollingDiscPenalty; simple model of a two-wheeler without steering
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2025-02-27
  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)
 20
 21SC = exu.SystemContainer()
 22mbs = SC.AddSystem()
 23
 24g = [0,0,-9.81]     #gravity in m/s^2
 25m = 0.5             #disc mass in kg
 26r = 0.1             #radius of disc in m
 27w = 0.01            #width of disc in m, just for drawing
 28v0 = 2
 29
 30L = 0.8             #length of board
 31
 32omega0 = [0,4*2*np.pi,0]      #initial angular velocity around y-axis
 33v0 = Skew(omega0) @ [0,0,r]   #initial angular velocity of center point
 34
 35#inertia for infinitely small ring:
 36inertiaRing = RigidBodyInertia(mass=1, inertiaTensor= np.diag([0.25*m*r**2, 0.5*m*r**2, 0.25*m*r**2]))
 37
 38#additional graphics for visualization of rotation:
 39graphicsWheel = graphics.Brick(centerPoint=[0,0,0],size=[1.2*r,2*w,1.2*r], color=graphics.color.lightred)
 40
 41oWheel0 = mbs.CreateRigidBody(referencePosition=[-0.5*L,0,r],
 42                             initialVelocity=v0,  
 43                             initialAngularVelocity=omega0,  
 44                             inertia=inertiaRing,  
 45                             gravity=g,  
 46                             graphicsDataList=[graphicsWheel]
 47                             )  
 48
 49oWheel1 = mbs.CreateRigidBody(referencePosition=[0.5*L,0,r],
 50                             initialVelocity=v0,  
 51                             initialAngularVelocity=omega0,  
 52                             inertia=inertiaRing,  
 53                             gravity=g,  
 54                             graphicsDataList=[graphicsWheel]
 55                             )  
 56
 57graphicsBody = graphics.Brick(centerPoint=[0,0,0],size=[L*1.1,0.6*r,0.6*r], color=graphics.color.grey)
 58oBody = mbs.CreateRigidBody(referencePosition=[0,0,r*1.2],
 59                             initialVelocity=v0,  
 60                             inertia=InertiaCuboid(2800, sideLengths=[L,0.5*r,0.5*r]),
 61                             gravity=g,  
 62                             graphicsDataList=[graphicsBody]
 63                             )  
 64
 65#ground body and marker
 66gGround = graphics.CheckerBoard(size=50, nTiles=50)
 67oGround = mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=[gGround])))
 68
 69mbs.CreateRevoluteJoint(itemNumbers=[oWheel0, oBody],
 70                        position = [0,0,0], axis=[0,1,0], useGlobalFrame=False,
 71                        axisRadius=0.1*r, axisLength=0.6*r)
 72
 73mbs.CreateRevoluteJoint(itemNumbers=[oWheel1, oBody],
 74                        position = [0,0,0], axis=[0,1,0], useGlobalFrame=False,
 75                        axisRadius=0.1*r, axisLength=0.6*r)
 76
 77oRolling0 = mbs.CreateRollingDiscPenalty(bodyNumbers=[oGround, oWheel0],
 78                             axisPosition=[0,0,0], axisVector=[0,1,0],
 79                             discRadius=r, 
 80                             contactStiffness=1e5, contactDamping=1e3, dryFriction=[0.4,0.4],
 81                             discWidth=0.1*r, color=graphics.color.blue)
 82
 83oRolling1 = mbs.CreateRollingDiscPenalty(bodyNumbers=[oGround, oWheel1],
 84                             axisPosition=[0,0,0], axisVector=[0,1,0],
 85                             discRadius=r, 
 86                             contactStiffness=1e5, contactDamping=1e3, dryFriction=[0.4,0.4],
 87                             discWidth=0.1*r, color=graphics.color.blue)
 88
 89
 90
 91#sensor for trace of contact point:
 92if not testIsActive:
 93    sTrail=mbs.AddSensor(SensorObject(objectNumber=oRolling0, storeInternal=True,#fileName='solution/rollingDiscTrail.txt', 
 94                               outputVariableType = exu.OutputVariableType.Position))
 95    
 96    
 97
 98
 99mbs.Assemble()
100
101simulationSettings = exu.SimulationSettings() #takes currently set values or default values
102
103
104stepSize=0.002
105tEnd = 1
106
107simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
108simulationSettings.timeIntegration.endTime = tEnd
109#simulationSettings.solution.file.writePeriod = 0.01
110simulationSettings.solution.sensors.writePeriod = 0.005
111simulationSettings.solution.file.write = (not testIsActive)
112simulationSettings.timeIntegration.verboseMode = 1
113simulationSettings.timeIntegration.newton.useModifiedNewton = True
114#simulationSettings.show.computationTime = True
115#simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse #does not work for initial accelerations
116
117SC.visualizationSettings.nodes.show = True
118SC.visualizationSettings.nodes.drawNodesAsPoint  = False
119SC.visualizationSettings.nodes.showBasis = True
120SC.visualizationSettings.nodes.basisSize = 0.015
121SC.visualizationSettings.view0.camera.perspective = 1
122
123if not testIsActive:
124    SC.renderer.Start()
125    SC.renderer.DoIdleTasks()
126
127mbs.SolveDynamic(simulationSettings)
128
129p0=mbs.GetObjectOutput(oRolling0, exu.OutputVariableType.Position)
130
131u = np.linalg.norm(p0)
132exu.Print('solution of createRollingDiscPenalty=',u) 
133
134exu.sys['testResult'] = u
135
136
137if not testIsActive:
138    SC.renderer.DoIdleTasks()
139    SC.renderer.Stop() #safely close rendering window!
140
141    ##++++++++++++++++++++++++++++++++++++++++++++++q+++++++
142    #plot results
143    if True:
144        mbs.PlotSensor(sTrail, componentsX=[0],components=[1], closeAll=True, title='wheel trail')