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 *
 18from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 19import exudyn.graphics as graphics #only import if it does not conflict
 20
 21SC = exu.SystemContainer()
 22mbs = SC.AddSystem()
 23
 24
 25#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 26#contact
 27
 28L = 5
 29n = 1
 30r = 0.2
 31mu = 0.4
 32m = 0.01
 33k = 1e3
 34d = 0.001*k
 35
 36ssx = 3 #search tree size
 37ssy = 3 #search tree size
 38ssz = 3 #search tree size
 39
 40tt = 0.2*r
 41zz = -0.04
 42markerList = []
 43p0 = np.array([0.,0,-0.5*tt])
 44gFloor = graphics.Brick(p0,[L,L,tt],graphics.color.lightgrey)
 45[meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gFloor)
 46#gDataList = [gFloor]
 47
 48gDataList = [graphics.CheckerBoard([0,0,0],size=L)]
 49
 50useRigidBody = True
 51#ns = 20
 52minP = np.array([1e10,1e10,1e10])
 53maxP = -minP
 54gContact = mbs.AddGeneralContact()
 55height = r*2.4
 56r2 = 0.2*r
 57
 58for i in range(n):
 59
 60    color4node = graphics.color.steelblue
 61    pS0 = [0,0,-zz]
 62    pS1 = [0,0,height-r-r2-zz]
 63    gList = [graphics.Sphere(point=pS0, radius=r, color= color4node, nTiles=24)]
 64    gList += [graphics.Sphere(point=pS1, radius=r2, color= color4node, nTiles=16)]
 65
 66    pRef = [0,0,r+zz]
 67    v0 = [0,-10*0.,0]
 68    omega0 = 20*np.array([0,0.02,10])
 69    rot0 = np.eye(3)
 70
 71    RBinertia = InertiaSphere(m, r*1)
 72
 73    RBinertia = InertiaCuboid(m/(2*r*2*r*height), [2*r,2*r,height])
 74
 75    [nMass, oMass] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
 76                            #nodeType=exu.NodeType.RotationRxyz,
 77                            nodeType=exu.NodeType.RotationRotationVector,
 78                            position=pRef, velocity=v0,
 79                            rotationMatrix=rot0,
 80                            angularVelocity=omega0,
 81                            gravity=[0., 0.,-9.81],
 82                            graphicsDataList=gList,
 83                            )
 84    #mbs.SetNodeParameter(nMass, 'VdrawSize', r*2)
 85    #mbs.SetNodeParameter(nMass, 'Vcolor', color4node)
 86
 87
 88
 89    mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
 90    #if useNodeMarker:
 91    #    markerList += [mNode]
 92    mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=pS0))
 93    gContact.AddSphereWithMarker(mBody, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
 94    mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=pS1))
 95    gContact.AddSphereWithMarker(mBody, radius=r2, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
 96
 97    #mbs.AddLoad(Force(markerNumber=mNode, loadVector= [0,forceY,0]))
 98
 99oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
100                                    visualization=VObjectGround(graphicsData=gDataList)))
101mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround))
102
103gContact.verboseMode = 1
104#gContact.sphereSphereContact = False
105gContact.frictionProportionalZone = 1e-5
106#[meshPoints,  meshTrigs] = RefineMesh(meshPoints,  meshTrigs)
107gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
108    pointList=meshPoints,  triangleList=meshTrigs)
109
110gContact.SetFrictionPairings(mu*np.eye(1))
111gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssz])
112
113
114
115#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
116
117mbs.Assemble()
118print(mbs)
119simulationSettings = exu.SimulationSettings() #takes currently set values or default values
120
121tEnd = 40
122h = 2*1e-4
123simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
124simulationSettings.timeIntegration.endTime = tEnd
125simulationSettings.solutionSettings.writeSolutionToFile = True
126simulationSettings.solutionSettings.solutionWritePeriod = 0.002
127#simulationSettings.displayComputationTime = True
128#simulationSettings.displayGlobalTimers= True
129simulationSettings.timeIntegration.verboseMode = 1
130# simulationSettings.timeIntegration.simulateInRealtime = True #turn this off to run faster!
131
132simulationSettings.displayStatistics = True
133
134SC.visualizationSettings.nodes.defaultSize = 0.01
135
136
137simulationSettings.solutionSettings.outputPrecision = 6
138simulationSettings.solutionSettings.exportVelocities = False
139simulationSettings.solutionSettings.exportAccelerations = False
140
141simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
142simulationSettings.timeIntegration.newton.useModifiedNewton = False
143simulationSettings.timeIntegration.discontinuous.iterationTolerance = 0.1
144
145
146SC.visualizationSettings.general.graphicsUpdateInterval=0.01
147#SC.visualizationSettings.general.circleTiling=50
148SC.visualizationSettings.general.drawCoordinateSystem=True
149SC.visualizationSettings.loads.show=False
150SC.visualizationSettings.nodes.show=True
151SC.visualizationSettings.nodes.showBasis = True
152SC.visualizationSettings.nodes.basisSize = r*2
153SC.visualizationSettings.nodes.defaultSize = r
154SC.visualizationSettings.nodes.tiling = 2
155SC.visualizationSettings.window.renderWindowSize=[1200,800]
156SC.visualizationSettings.openGL.multiSampling = 4
157
158showContact = False
159SC.visualizationSettings.contact.showSearchTreeCells = showContact
160SC.visualizationSettings.contact.showSearchTree = showContact
161SC.visualizationSettings.contact.showBoundingBoxes = showContact
162
163simulate=True
164if simulate:
165    useGraphics = True
166    if useGraphics:
167        exu.StartRenderer()
168        if 'renderState' in exu.sys:
169            SC.SetRenderState(exu.sys['renderState'])
170        mbs.WaitForUserToContinue()
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.WaitForRenderEngineStopFlag()
181        exu.StopRenderer() #safely close rendering window!
182
183if not simulate or True:
184    SC.visualizationSettings.general.autoFitScene = False
185
186    mbs.SolutionViewer()