ballBearingTest.py

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

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