newtonsCradle.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Several models of Newton's cradle with different coefficients of restitution
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2025-01-22
  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 *
 15import exudyn.graphics as graphics
 16import numpy as np
 17
 18SC = exu.SystemContainer()
 19mbs = SC.AddSystem()
 20exu.Print('EXUDYN version='+exu.config.Version())
 21
 22useRaytracer = True
 23use2contacts = True
 24yInit = 0
 25zInit = 0
 26radius=0.02
 27L = 0.25
 28d = 2*radius*1.02 #x-distance
 29rho = 7800
 30mass = rho * 4/3*pi*radius**3                         #mass in kg
 31
 32contactStiffness = 1e6                      #stiffness of spring-damper in N/m
 33contactDamping = 0    #damping constant in N/(m/s)
 34impactModel = 2
 35
 36tEnd = 2.5     #end time of simulation
 37stepSize = 0.5e-4
 38g = 9.81
 39exu.Print('impact vel=', np.sqrt(2*yInit*g))
 40
 41
 42color1=graphics.color.lightgrey[0:3]+[graphics.material.indexChrome]
 43color2=graphics.color.lightgrey2[0:3]+[graphics.material.indexChrome]
 44
 45gGround = [graphics.CheckerBoard(point=[0,-4*radius,0],
 46                                        normal=[0,1,0], size=1.2,
 47                                        color=color1,
 48                                        alternatingColor=color2)]
 49colorWall = graphics.color.dodgerblue[0:3]+[graphics.material.indexChrome]
 50gGround += [graphics.CheckerBoard(point=[0.6,-4*radius+0.3,0],
 51                                  normal=[-1,0,0],
 52                                  color=colorWall,
 53                                  alternatingColor=colorWall,
 54                                  size = 0.6, size2 = 1.2,
 55                                  nTiles=1)
 56                                  ]
 57gGround += [graphics.CheckerBoard(point=[0,-4*radius+0.3,-0.6],
 58                                  normal=[0,0,1],
 59                                  color=colorWall,
 60                                  alternatingColor=colorWall,
 61                                  size = 1.2, size2 = 0.6,
 62                                  nTiles=1)
 63                                  ]
 64
 65ground = mbs.CreateGround(referencePosition=[0,0,0],
 66                          graphicsDataList=gGround)
 67
 68nSpheres = 5
 69listCoeffRest = [1,0.95,0.8,1e-3] #1 is fully elastic, 0 fully plastic (but 0 is not allowed)
 70#listCoeffRest = [1]
 71
 72nZ = len(listCoeffRest) #number of mechanisms arranged along z
 73
 74for iz in range(nZ):
 75    restitutionCoefficient = listCoeffRest[iz]
 76    oListSpheres = []
 77    nListSpheres = []
 78    mListSpheres = []
 79
 80    if iz == 0:
 81        color = graphics.color.lightgrey[0:3]+[graphics.material.indexChrome]
 82    else:
 83        color = graphics.colorList[iz-1][0:3]+[graphics.material.indexChrome]
 84    zInit = 4*d-iz*4*d
 85    mbs.CreateGround(referencePosition=[0,yInit+L,zInit],
 86                    graphicsDataList=[graphics.Brick(centerPoint = [(nSpheres-1)*0.5*d,0,0],
 87                                                     size = [(1+nSpheres)*d,0.5*radius,0.25*radius],
 88                                                     color = graphics.material.chrome)])
 89
 90    for ix in range(nSpheres):
 91        xInit = ix*d
 92        yOff = 0
 93        angleZ = 0
 94        if ix == nSpheres-1:
 95            xInit += L
 96            yOff += L
 97            angleZ = pi/2
 98
 99        gSphere = graphics.Sphere(radius=radius, color=color, nTiles=32)
100        gString = graphics.Cylinder(pAxis=[0,0,0], vAxis=[0,L,0], radius=0.1*radius,
101                                    nTiles=32, color=graphics.material.glass)
102        #add mass point (this is a 3D object with 3 coordinates):
103        massPoint = mbs.CreateRigidBody(referencePosition=[xInit,yInit+yOff,zInit],
104                                        initialVelocity=[0,0,0],
105                                        referenceRotationMatrix=RotationMatrixZ(angleZ),
106                                        inertia=InertiaSphere(mass, radius),
107                                        gravity = [0,-g,0],
108                                        graphicsDataList=[gSphere, gString])
109
110        nMassPoint = mbs.GetObject(massPoint)['nodeNumber']
111        mMassPoint = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassPoint))
112
113        oListSpheres.append(massPoint)
114        nListSpheres.append(nMassPoint)
115        mListSpheres.append(mMassPoint)
116
117        #
118        mbs.CreateSphericalJoint(bodyNumbers=[ground, massPoint],
119                                 position=[ix*d, yInit+L, zInit],
120                                 show=False)
121
122        #distance constraint would work, but spherical joint allows to show string
123        # mbs.CreateDistanceConstraint(bodyNumbers=[ground, massPoint],
124        #                              localPosition0=[ix*d, yInit+L, zInit])
125
126
127    for i in range(nSpheres-1):
128
129        nData1 = mbs.AddNode(NodeGenericData(initialCoordinates=[0.1,0,0,0],
130                                            numberOfDataCoordinates=4))
131        oSSC = mbs.AddObject(ObjectContactSphereSphere(markerNumbers=[mListSpheres[i],mListSpheres[i+1]],
132                                                        nodeNumber=nData1,
133                                                        spheresRadii=[radius,radius],
134                                                        contactStiffness = contactStiffness,
135                                                        contactDamping = contactDamping,
136                                                        impactModel = impactModel,
137                                                        restitutionCoefficient = restitutionCoefficient,
138                                                        #minimumImpactVelocity = 1e-3,
139                                                        ))
140
141
142sPos=mbs.AddSensor(SensorBody(bodyNumber=massPoint, storeInternal=True,
143                                outputVariableType=exu.OutputVariableType.Position))
144sVel=mbs.AddSensor(SensorBody(bodyNumber=massPoint, storeInternal=True,
145                                outputVariableType=exu.OutputVariableType.Velocity))
146
147#exu.Print(mbs)
148mbs.Assemble()
149
150simulationSettings = exu.SimulationSettings()
151simulationSettings.solutionSettings.writeSolutionToFile = True
152simulationSettings.solutionSettings.solutionWritePeriod = 0.02
153simulationSettings.solutionSettings.sensorsWritePeriod = stepSize  #output interval
154simulationSettings.solutionSettings.solutionInformation = 'variation of restitution coefficients: 1, 0.95, 0.8, 1e-3'
155simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
156simulationSettings.timeIntegration.endTime = tEnd
157#simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-3
158#simulationSettings.timeIntegration.discontinuous.useRecommendedStepSize = False
159#simulationSettings.timeIntegration.discontinuous.maxIterations = 2
160
161simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
162
163# simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-2
164# simulationSettings.timeIntegration.discontinuous.useRecommendedStepSize = False
165
166simulationSettings.timeIntegration.newton.useModifiedNewton = True
167simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.9
168
169simulationSettings.displayStatistics = True
170simulationSettings.timeIntegration.verboseMode = 1
171
172SC.visualizationSettings.general.drawCoordinateSystem = False
173SC.visualizationSettings.general.showSolverInformation = False
174SC.visualizationSettings.window.renderWindowSize=[1200,800]
175SC.visualizationSettings.openGL.multiSampling=4
176SC.visualizationSettings.openGL.shadow = 0.2
177SC.visualizationSettings.openGL.lineWidth = 2
178SC.visualizationSettings.openGL.light0position = [-2.0, 4.0, 1.0, 1.0]
179SC.visualizationSettings.openGL.light1position = [-0.2, 0.2, -0.55, 1.0]
180SC.visualizationSettings.openGL.light1specular = 1
181SC.visualizationSettings.openGL.light1diffuse = 0
182SC.visualizationSettings.openGL.light1ambient = 0
183SC.visualizationSettings.openGL.light1specular = 1
184SC.visualizationSettings.openGL.light1quadraticAttenuation = 1
185SC.visualizationSettings.openGL.light1linearAttenuation = 0
186SC.visualizationSettings.loads.show = False
187
188
189SC.renderer.Start()              #start graphics visualization
190if 'renderState' in exu.sys: #reload last view
191    SC.renderer.SetState(exu.sys['renderState'])
192SC.renderer.DoIdleTasks()    #wait for pressing SPACE bar to continue
193
194#start solver:
195mbs.SolveDynamic(simulationSettings)
196
197SC.renderer.DoIdleTasks()#wait for pressing 'Q' to quit
198SC.renderer.Stop()               #safely close rendering window!
199
200#evaluate final (=current) output values
201# u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
202# exu.Print('u     =',u)
203uTotal = mbs.GetNodeOutput(nMassPoint, exu.OutputVariableType.CoordinatesTotal)
204exu.Print('uTotal=',uTotal[1])
205
206if useRaytracer:
207    SC.visualizationSettings.openGL.multiSampling = 1
208    SC.visualizationSettings.openGL.enableLight1 = False
209    SC.visualizationSettings.raytracer.numberOfThreads = 32 #number of threads!
210    SC.visualizationSettings.raytracer.enable = True
211    SC.visualizationSettings.raytracer.searchTreeFactor = 8
212
213
214mbs.SolutionViewer()