tippeTop.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  example of tippe top toy using GeneralContact, switching up-down after certain time
  5#           NOTE: this is only a demonstration of GeneralContact, but due to soft contact
  6#           this computational model may not fully capture effects of the real contact, which
  7#           would require further adjustment of parameters and comparison with experimental results!
  8#
  9# Author:   Johannes Gerstmayr
 10# Date:     2021-12-12
 11#
 12# 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.
 13#
 14#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 15
 16import exudyn as exu
 17from exudyn.itemInterface import *
 18import numpy as np
 19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 20import exudyn.graphics as graphics
 21
 22SC = exu.SystemContainer()
 23mbs = SC.AddSystem()
 24
 25
 26#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 27#contact
 28
 29L = 5
 30n = 1
 31r = 0.2
 32mu = 0.4
 33m = 0.01
 34k = 1e3
 35d = 0.001*k
 36
 37ssx = 3 #search tree size
 38ssy = 3 #search tree size
 39ssz = 3 #search tree size
 40
 41tt = 0.2*r
 42zz = -0.04
 43markerList = []
 44p0 = np.array([0.,0,-0.5*tt])
 45gFloor = graphics.Brick(p0,[L,L,tt],graphics.color.lightgrey)
 46[meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gFloor)
 47#gDataList = [gFloor]
 48
 49gDataList = [graphics.CheckerBoard([0,0,0],size=L)]
 50
 51useRigidBody = True
 52#ns = 20
 53minP = np.array([1e10,1e10,1e10])
 54maxP = -minP
 55gContact = mbs.AddGeneralContact()
 56height = r*2.4
 57r2 = 0.2*r
 58
 59for i in range(n):
 60
 61    color4node = graphics.color.steelblue
 62    pS0 = [0,0,-zz]
 63    pS1 = [0,0,height-r-r2-zz]
 64    gList = [graphics.Sphere(point=pS0, radius=r, color= color4node, nTiles=24)]
 65    gList += [graphics.Sphere(point=pS1, radius=r2, color= color4node, nTiles=16)]
 66
 67    pRef = [0,0,r+zz]
 68    v0 = [0,-10*0.,0]
 69    omega0 = 20*np.array([0,0.02,10])
 70    rot0 = np.eye(3)
 71
 72    RBinertia = InertiaSphere(m, r*1)
 73
 74    RBinertia = InertiaCuboid(m/(2*r*2*r*height), [2*r,2*r,height])
 75
 76    dictMass = mbs.CreateRigidBody(
 77                  inertia=RBinertia, 
 78                  nodeType=exu.NodeType.RotationRotationVector,
 79                  referencePosition=pRef, 
 80                  initialVelocity=v0,
 81                  referenceRotationMatrix=rot0,
 82                  initialAngularVelocity=omega0,
 83                  gravity=[0., 0., -9.81],
 84                  graphicsDataList=gList,
 85                  returnDict=True)
 86    [nMass, oMass] = [dictMass['nodeNumber'], dictMass['bodyNumber']]
 87
 88
 89
 90    mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
 91    #if useNodeMarker:
 92    #    markerList += [mNode]
 93    mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=pS0))
 94    gContact.AddSphereWithMarker(mBody, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
 95    mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=pS1))
 96    gContact.AddSphereWithMarker(mBody, radius=r2, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
 97
 98    #mbs.AddLoad(Force(markerNumber=mNode, loadVector= [0,forceY,0]))
 99
100oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
101                                    visualization=VObjectGround(graphicsData=gDataList)))
102mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround))
103
104gContact.verboseMode = 1
105#gContact.sphereSphereContact = False
106gContact.frictionProportionalZone = 1e-5
107#[meshPoints,  meshTrigs] = RefineMesh(meshPoints,  meshTrigs)
108gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
109    pointList=meshPoints,  triangleList=meshTrigs)
110
111gContact.SetFrictionPairings(mu*np.eye(1))
112gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssz])
113
114
115
116#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
117
118mbs.Assemble()
119print(mbs)
120simulationSettings = exu.SimulationSettings() #takes currently set values or default values
121
122tEnd = 40
123h = 2*1e-4
124simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
125simulationSettings.timeIntegration.endTime = tEnd
126simulationSettings.solution.file.write = True
127simulationSettings.solution.file.writePeriod = 0.002
128#simulationSettings.show.computationTime = True
129#simulationSettings.show.globalTimers= True
130simulationSettings.timeIntegration.verboseMode = 1
131# simulationSettings.timeIntegration.realtime.active = True #turn this off to run faster!
132
133simulationSettings.show.statistics = True
134
135SC.visualizationSettings.nodes.defaultSize = 0.01
136
137
138simulationSettings.solution.precision = 6
139simulationSettings.solution.file.export.velocities = False
140simulationSettings.solution.file.export.accelerations = False
141
142simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
143simulationSettings.timeIntegration.newton.useModifiedNewton = False
144simulationSettings.timeIntegration.discontinuous.iterationTolerance = 0.1
145
146
147SC.visualizationSettings.general.graphicsUpdateInterval=0.01
148#SC.visualizationSettings.general.circleTiling=50
149SC.visualizationSettings.view0.scene.drawCoordinateSystem=True
150SC.visualizationSettings.loads.show=False
151SC.visualizationSettings.nodes.show=True
152SC.visualizationSettings.nodes.showBasis = True
153SC.visualizationSettings.nodes.basisSize = r*2
154SC.visualizationSettings.nodes.defaultSize = r
155SC.visualizationSettings.nodes.tiling = 2
156SC.visualizationSettings.view0.window.renderWindowSize=[1200,800]
157SC.visualizationSettings.openGL.multiSampling = 4
158
159showContact = False
160SC.visualizationSettings.contact.showSearchTreeCells = showContact
161SC.visualizationSettings.contact.showSearchTree = showContact
162SC.visualizationSettings.contact.showBoundingBoxes = showContact
163
164simulate=True
165if simulate:
166    useGraphics = True
167    if useGraphics:
168        SC.renderer.Start()
169        SC.renderer.RestoreSavedState()
170        SC.renderer.DoIdleTasks()
171
172
173    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
174    simulationSettings.timeIntegration.endTime = tEnd
175    mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
176    # mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitMidpoint)
177    #mbs.SolveDynamic(simulationSettings)
178
179    if useGraphics:
180        SC.renderer.DoIdleTasks()
181        SC.renderer.Stop() #safely close rendering window!
182
183if not simulate or True:
184    SC.visualizationSettings.general.autoFitScene = False
185    
186    mbs.SolutionViewer()