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
 16import numpy as np
 17
 18testIsActive = exu.sys.get('testIsActive', False)
 19testIsActive = True #do test
 20
 21from exudyn.machines import GetBallBearingData, CreateBallBearing
 22
 23SC = exu.SystemContainer()
 24mbs = SC.AddSystem()
 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, 
 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*2)
 64
 65oGround = mbs.CreateGround(graphicsDataList=[gBack])
 66#fixed outer ring:
 67mOuterRing = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,0]))
 68
 69#create body for inner ring:
 70omega = 2*pi*2
 71gravity = [0,-g,0]
 72f = 0.002
 73shaftGraphics = graphics.SolidOfRevolution([0,0,0], axis, 
 74                            contour=[[-0.5*shaftLength,0],[-0.5*shaftLength,shaftRadius-f],
 75                                     [-0.5*shaftLength+f,shaftRadius],[0.5*shaftLength-f,shaftRadius],
 76                                     [0.5*shaftLength,shaftRadius-f],[0.5*shaftLength,0],], 
 77                            nTiles=32, color=graphics.color.steelblue)
 78
 79bodyInner = mbs.CreateRigidBody(referencePosition=[0,0,0],
 80                                nodeType=exu.NodeType.RotationRxyz,
 81                                inertia=InertiaCylinder(density, shaftLength, shaftRadius, axis=2),
 82                                initialAngularVelocity=[0,0,omega],
 83                                gravity = gravity,
 84                                graphicsDataList=[bearingGraphics['innerRingGraphics'],
 85                                                  graphics.Basis(origin=[0,0,width],length=1.5*shaftRadius),
 86                                                  shaftGraphics
 87                                                  ],
 88                            )
 89#add coordinate constraint to keep velocity constant:
 90mbs.CreateCoordinateConstraint(itemNumbers=[bodyInner, None], 
 91                               coordinates=[5,None],
 92                               velocityLevel=True, offset=omega, show=False)
 93
 94mInnerRing = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bodyInner, localPosition=[0,0,0]))
 95
 96#++++++++++++++++++++++++++++++++++
 97contactParametersRingBalls={'contactStiffness':5e6,'contactDamping':50,
 98                            'dynamicFriction':0.2,'contactStiffnessExponent':1,
 99                            #'restitutionCoefficient':0.5,'impactModel':2
100                            'frictionProportionalZone':1e-2,
101                            }
102
103listContactSensors = []
104bearingData['radiusBalls'] *= 1.01 #increase radius for prestressed-configuration
105bearingItems = CreateBallBearing(mbs, bearingData, 
106                                 mInnerRing, mOuterRing, density, density,
107                                 cageInitialAngularVelocity=[0,0,0], #not correct
108                                 ballsInitialAngularVelocity=[0,0,0],
109                                 gravity=gravity,
110                                 springStiffnessCage=1e5, springDampingCage=1e2,
111                                 contactParametersRingBalls=contactParametersRingBalls)
112
113#add sensor for trace
114sPosC=mbs.AddSensor(SensorObject(objectNumber=bearingItems['innerRingBallContacts'][0], storeInternal=True,
115                    outputVariableType=exu.OutputVariableType.Position))
116listContactSensors.append(sPosC)
117sPosC=mbs.AddSensor(SensorObject(objectNumber=bearingItems['outerRingBallContacts'][0], storeInternal=True,
118                    outputVariableType=exu.OutputVariableType.Position))
119listContactSensors.append(sPosC)
120
121#++++++++++++++++++++++++++++++++++
122#put outer ring graphics here for transparency:
123oGround = mbs.CreateGround(graphicsDataList=[bearingGraphics['outerRingGraphics'],
124                                             graphics.Brick(centerPoint=[outsideDiameter*0.6,0,0],
125                                                            size=[outsideDiameter*0.25,outsideDiameter*0.2,width], color=graphics.color.brown[0:3]+[0.4]),
126                                             graphics.Brick(centerPoint=[-outsideDiameter*0.6,0,0],
127                                                            size=[outsideDiameter*0.25,outsideDiameter*0.2,width], color=graphics.color.brown[0:3]+[0.4]),
128                                             ])
129             
130#++++++++++++++++++++++++++++++++++
131
132timeStartBB = 2 if not testIsActive else 0
133
134def UFforce(mbs, t, loadVector):
135    global timeStartBB
136    ts = timeStartBB
137    f0 = 200
138    force=0 if t < ts else min(f0*(t-ts),f0)
139    if t>2*ts:
140        force=max(-f0, f0-2*f0*(t-2*ts))
141    if t>3*ts:
142        force=min(0, -f0+f0*(t-3*ts))
143    
144    return [0,0,force]
145
146def UFtorque(mbs, t, loadVector):
147    ts = 3.5*timeStartBB
148    t0 = 8 #Nm
149    torque=0 if t < ts or t > ts+1 else min(t0*(t-ts),t0)
150    if t > ts+1 and t < ts+2:
151        torque = max(t0*(ts+2-t),0)
152    return [torque,0,0]
153
154mbs.CreateForce(itemNumber=bodyInner, localPosition=[0,0,0],
155                loadVectorUserFunction=UFforce)
156mbs.CreateTorque(itemNumber=bodyInner, 
157                 loadVectorUserFunction=UFtorque)
158
159
160
161mbs.Assemble()
162
163tEnd = 0.05
164if not testIsActive:
165    tEnd = 10
166
167stepSize = 1e-4
168
169simulationSettings = exu.SimulationSettings()
170simulationSettings.solution.file.write = (not testIsActive) #only the SolutionViewer reads it (#2492)
171simulationSettings.solution.file.writePeriod = 0.004
172simulationSettings.solution.sensors.writePeriod = stepSize  #output interval
173simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
174simulationSettings.timeIntegration.endTime = tEnd
175simulationSettings.show.statistics = True
176simulationSettings.show.computationTime = True
177#simulationSettings.timeIntegration.realtime.active = True
178#simulationSettings.timeIntegration.realtime.factor = 0.5
179#simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-2
180# simulationSettings.timeIntegration.discontinuous.maxIterations = 2
181#simulationSettings.timeIntegration.discontinuous.useRecommendedStepSize = False
182simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-8
183simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6
184simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
185simulationSettings.timeIntegration.newton.useModifiedNewton = True
186simulationSettings.parallel.numberOfThreads = 1
187
188simulationSettings.show.statistics = True
189simulationSettings.timeIntegration.verboseMode = 1
190
191SC.visualizationSettings.view0.window.renderWindowSize=[1600,1400]
192SC.visualizationSettings.openGL.multiSampling=2
193SC.visualizationSettings.openGL.light0.shadow = 0.25
194#SC.visualizationSettings.nodes.showBasis = True
195SC.visualizationSettings.loads.show = False
196SC.visualizationSettings.loads.drawSimplified=False
197SC.visualizationSettings.nodes.basisSize = radiusBalls*1.5
198SC.visualizationSettings.nodes.drawNodesAsPoint = False
199
200SC.visualizationSettings.sensors.traces.listOfPositionSensors = listContactSensors
201SC.visualizationSettings.sensors.traces.showPositionTrace = True if len(listContactSensors) else False
202SC.visualizationSettings.sensors.traces.timeSpan = 1.6
203
204if not testIsActive:
205    SC.renderer.Start()              #start graphics visualization
206    SC.renderer.DoIdleTasks()    #wait for pressing SPACE bar to continue
207
208#start solver:q
209mbs.SolveDynamic(simulationSettings)
210
211if not testIsActive:
212    SC.renderer.Stop()               #safely close rendering window!
213
214#%%++++
215testError = 0.01*np.linalg.norm(mbs.systemData.GetODE2Coordinates())
216exu.Print('solution of ballBearingTest=',testError)
217
218#%%++++++++++++++++++++++++++++++++++++++++++++++++++++
219exu.sys['testResult'] = testError
220
221if not testIsActive:
222    #%%
223    mbs.SolutionViewer()