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