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.GetVersionString())
 21
 22use2contacts = True
 23yInit = 0
 24zInit = 0
 25radius=0.02
 26L = 0.25
 27d = 2*radius*1.02 #x-distance
 28rho = 7800
 29mass = rho * 4/3*pi*radius**3                         #mass in kg
 30
 31contactStiffness = 1e6                      #stiffness of spring-damper in N/m
 32contactDamping = 0*1e-4*contactStiffness    #damping constant in N/(m/s)
 33impactModel = 2
 34
 35tEnd = 5     #end time of simulation
 36stepSize = 0.5e-4
 37g = 9.81
 38exu.Print('impact vel=', np.sqrt(2*yInit*g))
 39
 40ground = mbs.CreateGround(referencePosition=[0,0,0],
 41                          graphicsDataList=[graphics.CheckerBoard(point=[0,-4*radius,0], normal=[0,1,0], size=2)])
 42
 43nSpheres = 5
 44listCoeffRest = [1,0.95,0.8,1e-3] #1 is fully elastic, 0 fully plastic (but 0 is not allowed)
 45#listCoeffRest = [1]
 46
 47nZ = len(listCoeffRest) #number of mechanisms arranged along z
 48
 49for iz in range(nZ):
 50    restitutionCoefficient = listCoeffRest[iz]
 51    oListSpheres = []
 52    nListSpheres = []
 53    mListSpheres = []
 54
 55    color = graphics.colorList[iz]
 56    zInit = 4*d-iz*4*d
 57    mbs.CreateGround(referencePosition=[0,yInit+L,zInit],
 58                              graphicsDataList=[graphics.Brick(centerPoint = [(nSpheres-1)*0.5*d,0,0],
 59                                                               size = [(1+nSpheres)*d,0.5*radius,0.25*radius],
 60                                                               color = graphics.color.grey)])
 61
 62    for ix in range(nSpheres):
 63        xInit = ix*d
 64        yOff = 0
 65        if ix == nSpheres-1:
 66            xInit += L
 67            yOff += L
 68        #add mass point (this is a 3D object with 3 coordinates):
 69        massPoint = mbs.CreateRigidBody(referencePosition=[xInit,yInit+yOff,zInit],
 70                                        initialVelocity=[0,0,0],
 71                                        inertia=InertiaSphere(mass, radius),
 72                                        gravity = [0,-g,0],
 73                                        graphicsDataList=[graphics.Sphere(radius=radius,
 74                                                                          color=color,#graphics.color.orange,
 75                                                                          nTiles=64)])
 76        nMassPoint = mbs.GetObject(massPoint)['nodeNumber']
 77        mMassPoint = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassPoint))
 78
 79        oListSpheres.append(massPoint)
 80        nListSpheres.append(nMassPoint)
 81        mListSpheres.append(mMassPoint)
 82
 83        mbs.CreateDistanceConstraint(bodyNumbers=[ground, massPoint],
 84                                     localPosition0=[ix*d, yInit+L, zInit])
 85
 86
 87    for i in range(nSpheres-1):
 88
 89        nData1 = mbs.AddNode(NodeGenericData(initialCoordinates=[0.1,0,0,0],
 90                                            numberOfDataCoordinates=4))
 91        oSSC = mbs.AddObject(ObjectContactSphereSphere(markerNumbers=[mListSpheres[i],mListSpheres[i+1]],
 92                                                        nodeNumber=nData1,
 93                                                        spheresRadii=[radius,radius],
 94                                                        contactStiffness = contactStiffness,
 95                                                        contactDamping = contactDamping,
 96                                                        impactModel = impactModel,
 97                                                        restitutionCoefficient = restitutionCoefficient,
 98                                                        #minimumImpactVelocity = 1e-3,
 99                                                        ))
100
101
102# sPos=mbs.AddSensor(SensorBody(bodyNumber=massPoint, storeInternal=True,
103#                                 outputVariableType=exu.OutputVariableType.Position))
104# sVel=mbs.AddSensor(SensorBody(bodyNumber=massPoint, storeInternal=True,
105#                                 outputVariableType=exu.OutputVariableType.Velocity))
106
107#exu.Print(mbs)
108mbs.Assemble()
109
110simulationSettings = exu.SimulationSettings()
111simulationSettings.solutionSettings.writeSolutionToFile = True
112simulationSettings.solutionSettings.solutionWritePeriod = 0.02
113simulationSettings.solutionSettings.sensorsWritePeriod = stepSize  #output interval
114simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
115simulationSettings.timeIntegration.endTime = tEnd
116#simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-3
117#simulationSettings.timeIntegration.discontinuous.useRecommendedStepSize = False
118#simulationSettings.timeIntegration.discontinuous.maxIterations = 2
119
120simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
121
122# simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-2
123# simulationSettings.timeIntegration.discontinuous.useRecommendedStepSize = False
124
125simulationSettings.timeIntegration.newton.useModifiedNewton = True
126
127simulationSettings.displayStatistics = True
128simulationSettings.timeIntegration.verboseMode = 1
129
130SC.visualizationSettings.window.renderWindowSize=[1600,1200]
131SC.visualizationSettings.openGL.multiSampling=4
132SC.visualizationSettings.openGL.shadow = 0.2
133SC.visualizationSettings.openGL.lineWidth = 2
134SC.visualizationSettings.openGL.light0position = [2,4,-1,0]
135SC.visualizationSettings.loads.show = False
136
137exu.StartRenderer()              #start graphics visualization
138mbs.WaitForUserToContinue()    #wait for pressing SPACE bar to continue
139
140#start solver:
141mbs.SolveDynamic(simulationSettings,
142                 #solverType=exu.DynamicSolverType.TrapezoidalIndex2
143                 )
144
145SC.WaitForRenderEngineStopFlag()#wait for pressing 'Q' to quit
146exu.StopRenderer()               #safely close rendering window!
147
148#evaluate final (=current) output values
149# u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
150# exu.Print('u     =',u)
151uTotal = mbs.GetNodeOutput(nMassPoint, exu.OutputVariableType.CoordinatesTotal)
152exu.Print('uTotal=',uTotal[1])
153
154mbs.SolutionViewer()