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)