ballBearningModel.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test bearing model with two ball bearings on rigid shaft
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2025-05-17
  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
 18from exudyn.machines import GetBallBearingData, CreateBallBearing
 19
 20
 21useGraphics = True #without test
 22SC = exu.SystemContainer()
 23mbs = SC.AddSystem()
 24
 25
 26#+++++++++++++++++++++++++++++++++++++++++++++++++
 27#simple bearing example
 28#+++++++++++++++++++++++++++++++++++++++++++
 29#6010/C3 example data sheet:
 30outsideDiameter=2*0.040
 31boreDiameter=2*0.025
 32width=0.016
 33## measured:
 34nBalls=14
 35radiusCage=0.0325
 36radiusBalls=0.00873/2 #usually diameter is given
 37outerRingShoulderRadius = 0.035125
 38innerRingShoulderRadius = 0.029875
 39## computed:
 40innerGrooveRadius = radiusBalls*1.25 #assumption
 41outerGrooveRadius = radiusBalls*1.25 #assumption
 42innerEdgeChamfer = 0.001  #in fact, it is round
 43outerEdgeChamfer = 0.001
 44
 45#+++++++++++++++++++++++++++++++++++++++++++
 46#create overall data for bearing, incl. drawing
 47axis=[0,0,1]
 48bearingData = GetBallBearingData(axis, outsideDiameter, boreDiameter, width, nBalls,
 49                                 radiusBalls=radiusBalls, ballsAngleOffset=0,
 50                                 radiusCage=radiusCage, heightCage=radiusBalls*0.2,
 51                                 innerGrooveRadius=innerGrooveRadius, outerGrooveRadius=outerGrooveRadius,
 52                                 innerEdgeChamfer=innerEdgeChamfer, outerEdgeChamfer=outerEdgeChamfer,
 53                                 innerRingShoulderRadius=innerRingShoulderRadius,outerRingShoulderRadius=outerRingShoulderRadius,
 54                                 )
 55bearingGraphics = graphics.BallBearingRings(**bearingData, nTilesRings=64)
 56
 57
 58#+++++++++++++++++++++++++++++++++++++++++++
 59#we need to create markers for outer ring and inner ring:
 60shaftRadius = boreDiameter*0.5
 61shaftLength = 0.120
 62density = 7800
 63gBack = graphics.CheckerBoard([0,0,-shaftLength*0.55], size=outsideDiameter*4)
 64
 65oGround = mbs.CreateGround(graphicsDataList=[gBack])
 66
 67#create body for inner ring:
 68omega = 2*pi*2
 69gravity = [0,-g,0]
 70f = 0.002
 71shaftGraphics = graphics.SolidOfRevolution([0,0,0], axis,
 72                            contour=[[-0.5*shaftLength,0],[-0.5*shaftLength,shaftRadius-f],
 73                                     [-0.5*shaftLength+f,shaftRadius],[0.5*shaftLength-f,shaftRadius],
 74                                     [0.5*shaftLength,shaftRadius-f],[0.5*shaftLength,0],],
 75                            nTiles=64, color=graphics.color.steelblue)
 76
 77bodyInner = mbs.CreateRigidBody(referencePosition=[0,0,0],
 78                                nodeType=exu.NodeType.RotationRxyz,
 79                                inertia=InertiaCylinder(density, shaftLength, shaftRadius, axis=2),
 80                                initialAngularVelocity=[0,0,omega],
 81                                gravity = gravity,
 82                            )
 83
 84#add coordinate constraint to keep velocity constant:
 85mbs.CreateCoordinateConstraint(bodyNumbers=[bodyInner, None],
 86                               coordinates=[5,None],
 87                               velocityLevel=True, offset=omega, show=False)
 88
 89graphicsDataListInner = [shaftGraphics, graphics.Basis(origin=[0,0,width],length=1.5*shaftRadius),]
 90offsetList = [-shaftLength*0.3, shaftLength*0.3, ]
 91for offsetZ in offsetList:
 92    #fixed outer ring:
 93    mOuterRing = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,offsetZ]))
 94    #rotating inner ring:
 95    mInnerRing = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bodyInner, localPosition=[0,0,offsetZ]))
 96
 97    graphicsDataListInner+=[graphics.Move(bearingGraphics['innerRingGraphics'],[0,0,offsetZ])]
 98
 99    #++++++++++++++++++++++++++++++++++
100    contactParametersRingBalls={'contactStiffness':5e6,'contactDamping':50,
101                                'dynamicFriction':0.2,'contactStiffnessExponent':1,
102                                #'restitutionCoefficient':0.5,'impactModel':2
103                                'frictionProportionalZone':1e-2,
104                                }
105
106    listContactSensors = []
107    bearingData['radiusBalls'] *= 1.01 #increase radius for prestressed-configuration
108    bearingItems = CreateBallBearing(mbs, bearingData,
109                                     mInnerRing, mOuterRing, density, density,
110                                     cageInitialAngularVelocity=[0,0,0], #not correct
111                                     ballsInitialAngularVelocity=[0,0,0],
112                                     gravity=gravity,
113                                     springStiffnessCage=1e5, springDampingCage=1e2,
114                                     contactParametersRingBalls=contactParametersRingBalls)
115
116    #add sensor for trace
117    sPosC=mbs.AddSensor(SensorObject(objectNumber=bearingItems['innerRingBallContacts'][0], storeInternal=True,
118                        outputVariableType=exu.OutputVariableType.Position))
119    listContactSensors.append(sPosC)
120    sPosC=mbs.AddSensor(SensorObject(objectNumber=bearingItems['outerRingBallContacts'][0], storeInternal=True,
121                        outputVariableType=exu.OutputVariableType.Position))
122    listContactSensors.append(sPosC)
123
124    #++++++++++++++++++++++++++++++++++
125    #put outer ring graphics here for transparency:
126    oGroundDrawing = mbs.CreateGround(referencePosition=[0,0,offsetZ],
127                                      graphicsDataList=[bearingGraphics['outerRingGraphics'],
128                                                 graphics.Brick(centerPoint=[outsideDiameter*0.6,0,0],
129                                                                size=[outsideDiameter*0.25,outsideDiameter*0.2,width], color=graphics.color.brown[0:3]+[0.4]),
130                                                 graphics.Brick(centerPoint=[-outsideDiameter*0.6,0,0],
131                                                                size=[outsideDiameter*0.25,outsideDiameter*0.2,width], color=graphics.color.brown[0:3]+[0.4]),
132                                                 ])
133
134mbs.SetObjectParameter(bodyInner, 'VgraphicsData', graphicsDataListInner)
135#++++++++++++++++++++++++++++++++++
136
137def UFforce(mbs, t, loadVector):
138    ts = 2
139    f0 = 200
140    force=0 if t < ts else min(f0*(t-ts),f0)
141    if t>2*ts:
142        force=max(-f0, f0-2*f0*(t-2*ts))
143    if t>3*ts:
144        force=min(0, -f0+f0*(t-3*ts))
145
146    return [0,0,force]
147
148def UFtorque(mbs, t, loadVector):
149    ts = 7
150    t0 = 8 #Nm
151    torque=0 if t < ts or t > ts+1 else min(t0*(t-ts),t0)
152    if t > ts+1 and t < ts+2:
153        torque = max(t0*(ts+2-t),0)
154    return [torque,0,0]
155
156mbs.CreateForce(bodyNumber=bodyInner, localPosition=[0,0,0],
157                loadVectorUserFunction=UFforce)
158mbs.CreateTorque(bodyNumber=bodyInner,
159                 loadVectorUserFunction=UFtorque)
160
161
162
163mbs.Assemble()
164
165stepSize = 1e-4
166tEnd = 10
167
168simulationSettings = exu.SimulationSettings()
169simulationSettings.solutionSettings.writeSolutionToFile = True
170simulationSettings.solutionSettings.solutionWritePeriod = 0.004
171simulationSettings.solutionSettings.sensorsWritePeriod = stepSize  #output interval
172simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
173simulationSettings.timeIntegration.endTime = tEnd
174simulationSettings.displayStatistics = True
175simulationSettings.displayComputationTime = True
176#simulationSettings.timeIntegration.simulateInRealtime = True
177#simulationSettings.timeIntegration.realtimeFactor = 0.5
178#simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-2
179# simulationSettings.timeIntegration.discontinuous.maxIterations = 2
180#simulationSettings.timeIntegration.discontinuous.useRecommendedStepSize = False
181simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-8
182simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6
183simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
184simulationSettings.timeIntegration.newton.useModifiedNewton = True
185simulationSettings.parallel.numberOfThreads = 1
186
187simulationSettings.displayStatistics = True
188simulationSettings.timeIntegration.verboseMode = 1
189
190SC.visualizationSettings.window.renderWindowSize=[1600,1400]
191SC.visualizationSettings.openGL.multiSampling=2
192SC.visualizationSettings.openGL.shadow = 0.25
193#SC.visualizationSettings.nodes.showBasis = True
194SC.visualizationSettings.loads.show = False
195SC.visualizationSettings.loads.drawSimplified=False
196SC.visualizationSettings.nodes.basisSize = radiusBalls*1.5
197#SC.visualizationSettings.nodes.tiling = 32
198SC.visualizationSettings.nodes.drawNodesAsPoint = False
199SC.visualizationSettings.connectors.showContact = False
200#SC.visualizationSettings.general.sphereTiling = 32
201
202SC.visualizationSettings.sensors.traces.listOfPositionSensors = listContactSensors
203SC.visualizationSettings.sensors.traces.showPositionTrace = True if len(listContactSensors) else False
204SC.visualizationSettings.sensors.traces.timeSpan = 1.6
205
206SC.renderer.Start()              #start graphics visualization
207SC.renderer.DoIdleTasks()    #wait for pressing SPACE bar to continue
208
209#start solver:q
210mbs.SolveDynamic(simulationSettings)
211
212# SC.renderer.DoIdleTasks()#wait for pressing 'Q' to quit
213SC.renderer.Stop()               #safely close rendering window!
214
215if True:
216    #%%
217    mbs.SolutionViewer()