rollerBearningModel.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for SphereSphereContact
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2024-11-01
  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
 18useGraphics = True #without test
 19SC = exu.SystemContainer()
 20mbs = SC.AddSystem()
 21exu.Print('EXUDYN version='+exu.config.Version())
 22
 23w = 0.12
 24radius=0.1
 25radiusOuter = 6*radius
 26radiusInner = 4*radius
 27
 28vyInit = 0
 29mass = 1.6                              #mass in kg
 30contactStiffness = 1e8                  #stiffness of spring-damper in N/m
 31contactDamping = 5e-4*contactStiffness  #damping constant in N/(m/s)
 32restitutionCoefficient = 0.5
 33impactModel = 2*0
 34
 35tEnd = 2     #end time of simulation
 36stepSize = 5*2e-5 #*10
 37g = 10
 38#print('impact vel=', np.sqrt(2*yInit*g))
 39
 40
 41#ground node
 42nGround1=mbs.AddNode(NodePointGround(referenceCoordinates = [0,0,0]))
 43
 44nSpheres = 8
 45listMasses = []
 46for i in range(nSpheres-0):
 47    rr = radiusOuter-radius
 48    phi = i/nSpheres*2*pi
 49    x = rr * cos(phi)
 50    y = rr * sin(phi)
 51    #add mass point (this is a 3D object with 3 coordinates):
 52    massPoint = mbs.CreateRigidBody(referencePosition=[x,y,0],
 53                                    initialVelocity=[0,0,0],
 54                                    inertia=InertiaCylinder(7800, w, radius, axis=2),
 55                                    gravity = [0,-g,0],
 56                                    graphicsDataList=[graphics.Cylinder(pAxis=[0,0,-w*0.5], vAxis=[0,0,w], radius=radius, color=graphics.color.orange, nTiles=32)],
 57                                    )
 58    listMasses.append(massPoint)
 59    nMassPoint = mbs.GetObject(massPoint)['nodeNumber']
 60    mbs.SetNodeParameter(nMassPoint, 'Vshow', True)
 61
 62gHollow = graphics.Cylinder(pAxis=[0,0,-0.5*w], vAxis=[0,0,w], radius=radiusOuter+radius,
 63                            radiusInner=radiusOuter, color=[0.6,0.6,0.6,0.4], nTiles=2*64)
 64
 65gBack = graphics.CheckerBoard([0,0,-w*2], size=radiusOuter*2.5)
 66oGround = mbs.CreateGround(referencePosition=[0,0,0], graphicsDataList=[gHollow,gBack])
 67
 68bodyInner = mbs.CreateRigidBody(referencePosition=[0,0,0],
 69                            inertia=InertiaCylinder(7800, w, radiusInner, axis=2),
 70                            initialAngularVelocity=[0,0,2*pi*2],
 71                            gravity = [0,-g,0],
 72                            graphicsDataList=[graphics.Cylinder(pAxis=[0,0,-w*2], vAxis=[0,0,4*w], radius=radiusInner, color=graphics.color.dodgerblue, nTiles=64),
 73                                              graphics.Basis(origin=[0,0,2.1*w],length=2*radius),
 74                                              ],
 75                            )
 76mbs.CreateRevoluteJoint(bodyNumbers=[oGround, bodyInner], position=[0,0,0], axis=[0,0,1])
 77
 78mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround1))
 79mBodyInner = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bodyInner))
 80
 81for oMass in listMasses:
 82    mMass = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass))
 83
 84    for k, mRef in enumerate([mGround,mBodyInner]):
 85        isHollowSphere1 = (k==0)
 86        radiusBase = [radiusOuter,radiusInner][k]
 87        nData1 = mbs.AddNode(NodeGenericData(initialCoordinates=[0.1*0,0,0,0],
 88                                            numberOfDataCoordinates=4))
 89        oSSC = mbs.AddObject(ObjectContactSphereSphere(markerNumbers=[mMass,mRef],
 90                                                       nodeNumber=nData1,
 91                                                       spheresRadii=[radius*(1+1e-3),radiusBase],
 92                                                       contactStiffness = contactStiffness,
 93                                                       isHollowSphere1=isHollowSphere1,
 94                                                       dynamicFriction=0.2,
 95                                                       #contactStiffnessExponent=1.,
 96                                                       contactDamping = contactDamping,
 97                                                       impactModel = impactModel,
 98                                                       frictionProportionalZone = 0.001,
 99                                                       # restitutionCoefficient = restitutionCoefficient,
100                                                       minimumImpactVelocity = 1e-2,
101                                                       visualization=VObjectContactSphereSphere(show=True),
102                                                       ))
103
104
105sPos=mbs.AddSensor(SensorBody(bodyNumber=massPoint, storeInternal=True,
106                                outputVariableType=exu.OutputVariableType.Position))
107sVel=mbs.AddSensor(SensorBody(bodyNumber=massPoint, storeInternal=True,
108                                outputVariableType=exu.OutputVariableType.Velocity))
109
110#exu.Print(mbs)
111mbs.Assemble()
112
113simulationSettings = exu.SimulationSettings()
114simulationSettings.solutionSettings.writeSolutionToFile = True
115simulationSettings.solutionSettings.solutionWritePeriod = 0.005
116simulationSettings.solutionSettings.sensorsWritePeriod = stepSize  #output interval
117simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
118simulationSettings.timeIntegration.endTime = tEnd
119#simulationSettings.timeIntegration.simulateInRealtime = True
120#simulationSettings.timeIntegration.realtimeFactor = 0.5
121#simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-2
122#simulationSettings.timeIntegration.discontinuous.maxIterations = 3
123# simulationSettings.timeIntegration.discontinuous.useRecommendedStepSize = False
124simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-8
125simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6
126
127simulationSettings.timeIntegration.newton.useModifiedNewton = True
128
129simulationSettings.displayStatistics = True
130simulationSettings.timeIntegration.verboseMode = 1
131
132SC.visualizationSettings.window.renderWindowSize=[1600,1400]
133SC.visualizationSettings.openGL.multiSampling=4
134SC.visualizationSettings.openGL.shadow = 0.25
135SC.visualizationSettings.nodes.showBasis = True
136SC.visualizationSettings.nodes.basisSize = radius*1.5
137# SC.visualizationSettings.contact.showSpheres = True
138# SC.visualizationSettings.contact.showTori = True
139
140SC.renderer.Start()              #start graphics visualization
141SC.renderer.DoIdleTasks()    #wait for pressing SPACE bar to continue
142
143#start solver:q
144mbs.SolveDynamic(simulationSettings)
145
146SC.renderer.DoIdleTasks()#wait for pressing 'Q' to quit
147SC.renderer.Stop()               #safely close rendering window!
148
149#evaluate final (=current) output values
150# u = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
151# exu.Print('u     =',u)
152uTotal = mbs.GetNodeOutput(nMassPoint, exu.OutputVariableType.CoordinatesTotal)
153exu.Print('uTotal=',uTotal[0])
154
155mbs.SolutionViewer()
156
157
158#plot results:
159if True:
160    mbs.PlotSensor([sPos,sVel], components=[1,1], closeAll=True)
161    import matplotlib.pyplot as plt
162    plt.show(block=True) #for figures to stay open at end of plot routines
163
164#+++++++++++++++++++++++++++++++++++++++++++++++++++++