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