generalContactFrictionTests.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  test friction of spheres and spheres-trigs
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2021-12-06
  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
 16
 17import numpy as np
 18
 19useGraphics = True #without test
 20#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 21#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 22try: #only if called from test suite
 23    from modelUnitTests import exudynTestGlobals #for globally storing test results
 24    useGraphics = exudynTestGlobals.useGraphics
 25except:
 26    class ExudynTestGlobals:
 27        pass
 28    exudynTestGlobals = ExudynTestGlobals()
 29#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 30
 31SC = exu.SystemContainer()
 32mbs = SC.AddSystem()
 33
 34nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
 35
 36#isPerformanceTest = exudynTestGlobals.isPerformanceTest
 37#useGraphics = False
 38# isPerformanceTest = True
 39# tEnd = 0.1
 40# if isPerformanceTest: tEnd *= 0.5
 41
 42#%%+++++++++++++++++++++++++++++++++
 43#sphere-sphere with coordinate constraints, prestressed; fixed torque on one side, linear increasing torque on other side
 44#sphere on ground, rolling
 45#cube on ground, sliding (f=[f, f*mu*t, 0]), tangential force changing
 46#cube on ground with initial velocity
 47#cube-cube contact (meshed)
 48
 49
 50L = 1 #surrounding
 51a = 0.1 #base dimention of objects
 52r = 0.5*a #radius
 53t = 0.25*a #thickness
 54
 55#contact coefficients:
 56mu = 0.8      #dry friction
 57m = 0.025     #mass
 58k = 1e3       #(linear) normal contact stiffness
 59d = 2*1e-4*k  #(linear) contact damping
 60gFact = 10
 61g = [0,0,-gFact]
 62
 63gContact = mbs.AddGeneralContact()
 64gContact.verboseMode = 1
 65#gContact.sphereSphereContact = False
 66gContact.frictionProportionalZone = 1e-3
 67#gContact.excludeDuplicatedTrigSphereContactPoints = False
 68fricMat = mu*np.eye(2)
 69fricMat[0,1] = 0.2
 70fricMat[1,0] = 0.2
 71gContact.SetFrictionPairings(fricMat)
 72gContact.SetSearchTreeCellSize(numberOfCells=[4,4,4])
 73
 74#%% ground
 75p0 = np.array([0,0,-0.5*t])
 76color4wall = [0.9,0.9,0.7,0.5]
 77addNormals = False
 78gFloor = graphics.Brick(p0,[L,L,t],graphics.color.steelblue,addNormals)
 79gFloorAdd = graphics.Brick(p0+[-0.5*L,0,a],[t,L,2*a],color4wall,addNormals)
 80gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
 81gFloorAdd = graphics.Brick(p0+[ 0.5*L,0,a],[t,L,2*a],color4wall,addNormals)
 82gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
 83gFloorAdd = graphics.Brick(p0+[0,-0.5*L,a],[L,t,2*a],color4wall,addNormals)
 84gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
 85gFloorAdd = graphics.Brick(p0+[0, 0.5*L,a],[L,t,2*a],color4wall,addNormals)
 86gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
 87
 88bb = 0.75*a
 89bh = 0.25*a
 90p1 = np.array([0.5*L,0.5*L,0])
 91gFloorAdd = graphics.Brick(p1+[-bb*3, -bb, 0.5*bh],[bb,bb,bh],color4wall,addNormals)
 92gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
 93gFloorAdd = graphics.Brick(p1+[-bb*2, -bb, 1.5*bh],[bb,bb,bh],color4wall,addNormals)
 94gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
 95gFloorAdd = graphics.Brick(p1+[-bb*1, -bb, 2.5*bh],[bb,bb,bh],color4wall,addNormals)
 96gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
 97
 98gDataList = [gFloor]
 99
100
101nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0] ))
102mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))
103mGroundC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
104
105[meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gFloor)
106#[meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
107# [meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
108gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
109    pointList=meshPoints,  triangleList=meshTrigs)
110
111if True: #looses color
112    gFloor = graphics.FromPointsAndTrigs(meshPoints, meshTrigs, color=color4wall) #show refined mesh
113    gDataList = [gFloor]
114
115evalNodes = [] #collect nodes that are evaluated for test
116#%%++++++++++++++++++++++++++++++++++++++++++++
117#free rolling sphere:
118gList = [graphics.Sphere(point=[0,0,0], radius=r, color= graphics.color.red, nTiles=24)]
119omega0 = -4.*np.array([5,1.,0.])
120pRef = [-0.4*L,-0.4*L,r-0*m*gFact/k]
121RBinertia = InertiaSphere(m, r)
122[nMass, oMass] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
123                      nodeType=exu.NodeType.RotationRotationVector,
124                      position=pRef,
125                      velocity=-np.cross([0,0,r],omega0),
126                      rotationMatrix=RotationMatrixX(0.),
127                      angularVelocity=omega0,
128                      #gravity=g,
129                      graphicsDataList=gList,
130                      )
131
132nNode0 = nMass
133mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
134mbs.AddLoad(Force(markerNumber=mNode, loadVector= [0,0,-k*r*0.01])) #==> uz = 2*r*0.01
135exu.Print('expect u0z=',2*r*0.01)
136gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
137if useGraphics:
138    sNode0 = mbs.AddSensor(SensorNode(nodeNumber=nNode0, storeInternal=True, #fileName='solution/contactNode0.txt',
139                                      outputVariableType=exu.OutputVariableType.Displacement))
140    vNode0 = mbs.AddSensor(SensorNode(nodeNumber=nNode0, storeInternal=True, #fileName='solution/contactNode0Vel.txt',
141                                      outputVariableType=exu.OutputVariableType.Velocity))
142evalNodes += [nMass]
143
144#%%++++++++++++++++++++++++++++++++++++++++++++
145#free rolling sphere at midpoint, many triangles in close contact; slowly go through critical points:
146gList = [graphics.Sphere(point=[0,0,0], radius=r, color= graphics.color.yellow, nTiles=24)]
147omega0 = -1e-12*np.array([1,0.1,0.])
148pRef = [1e-15,-1e-14,r-2*m*gFact/k]
149RBinertia = InertiaSphere(m, r)
150[nMass, oMass] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
151                      nodeType=exu.NodeType.RotationRotationVector,
152                      position=pRef,
153                      velocity=-np.cross([0,0,r],omega0),
154                      rotationMatrix=RotationMatrixX(0.),
155                      angularVelocity=omega0,
156                      gravity=g,
157                      graphicsDataList=gList,
158                      )
159
160nNode1 = nMass
161mNode1 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
162#mbs.AddLoad(Force(markerNumber=mNode1, loadVector= [0,0,-k*r*0.01])) #==> uz = 2*r*0.01
163#exu.Print('expect u0z=',2*r*0.01)
164gContact.AddSphereWithMarker(mNode1, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
165
166sNode1 = mbs.AddSensor(SensorNode(nodeNumber=nNode1, storeInternal=True, #fileName='solution/contactNode1.txt',
167                                  outputVariableType=exu.OutputVariableType.Displacement))
168# vNode1 = mbs.AddSensor(SensorNode(nodeNumber=nNode0, storeInternal=True, #fileName='solution/contactNode0Vel.txt',
169#                                   outputVariableType=exu.OutputVariableType.Velocity))
170
171
172#%%++++++++++++++++++++++++++++++++++++++++++++
173#fixed pressure tests:
174pf = np.array([-1.2*L,0,0])
175nGroundF = mbs.AddNode(NodePointGround(referenceCoordinates=pf ))
176mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGroundF))
177gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
178gDataList += [graphics.Sphere(point=pf, radius=r, color= graphics.color.grey, nTiles=24)]
179
180gList = [graphics.Sphere(point=[0,0,0], radius=r, color= graphics.color.lightgreen, nTiles=24)]
181
182pRef = pf+[0,2*r,0] #[-0.4*L,-0.4*L,r-m*gFact/k]
183RBinertia = InertiaSphere(m, r)
184[nMassF, oMassF] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
185                      nodeType=exu.NodeType.RotationRotationVector,
186                      position=pRef,
187                      rotationMatrix=RotationMatrixX(0.),
188                      #gravity=g,
189                      graphicsDataList=gList,
190                      )
191
192mC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMassF, coordinate=0))
193mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mC]))
194
195mNodeF = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassF))
196mbs.AddLoad(Force(markerNumber=mNodeF, loadVector= [0,-k*r*0.1,0])) #==> u =  k*r*0.1/(0.5*k) = 2*r*0.1
197exu.Print('expect uFy=',2*r*0.1)
198gContact.AddSphereWithMarker(mNodeF, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
199if useGraphics:
200    sNodeF = mbs.AddSensor(SensorNode(nodeNumber=nMassF, storeInternal=True, #fileName='solution/contactNodeF.txt',
201                                      outputVariableType=exu.OutputVariableType.Displacement))
202evalNodes += [nMassF]
203
204#%%++++++++++++++++++++++++++++++++++++++++++++
205# sliding between spheres:
206pr = np.array([-1.2*L,0.5*L,0])
207nGroundF2 = mbs.AddNode(NodePointGround(referenceCoordinates=pr ))
208mNode2 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGroundF2))
209gContact.AddSphereWithMarker(mNode2, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
210gDataList += [graphics.Sphere(point=pr, radius=r, color= graphics.color.lightgrey, nTiles=24)]
211
212gList = [graphics.Sphere(point=[0,0,0], radius=r, color= graphics.color.lightred, nTiles=24)]
213
214dRol = r*0.01
215pRef = pr+[0,2*r-2*dRol,0] #force=k*r*0.01
216fRol = k*dRol
217exu.Print('force rolling=', fRol, ', torque=', fRol*mu*r)
218RBinertia = InertiaSphere(m, r)
219[nMassR, oMassR] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
220                      nodeType=exu.NodeType.RotationRotationVector,
221                      position=pRef,
222                      rotationMatrix=RotationMatrixX(0.),
223                      #angularVelocity=[10,0,0],
224                      #gravity=g,
225                      graphicsDataList=gList,
226                      )
227
228mC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMassR, coordinate=0))
229mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mC]))
230mC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMassR, coordinate=1))
231mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mC]))
232mC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMassR, coordinate=2))
233mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mC]))
234
235mNodeR = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassR))
236
237def UFtorque(mbs, t, loadVector):
238    torque = 10*t*fRol*mu*r
239    if t > 0.3:
240        torque = 0
241    return [torque,0,0]
242mbs.AddLoad(Torque(markerNumber=mNodeR, loadVectorUserFunction=UFtorque,
243                   loadVector= [1,0,0])) #==> u =  k*r*0.1/(0.5*k) = 2*r*0.1
244
245gContact.AddSphereWithMarker(mNodeR, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
246if useGraphics:
247    sNodeR = mbs.AddSensor(SensorNode(nodeNumber=nMassR, storeInternal=True, #fileName='solution/contactNodeR.txt',
248                                      outputVariableType=exu.OutputVariableType.Rotation))
249    vNodeR = mbs.AddSensor(SensorNode(nodeNumber=nMassR, storeInternal=True, #fileName='solution/contactNodeRvel.txt',
250                                      outputVariableType=exu.OutputVariableType.AngularVelocity))
251
252evalNodes += [nMassR]
253
254#%%++++++++++++++++++++++++++++++++++++++++++++
255#sphere on stairs
256#%%++++++++++++++++++++++++++++++++++++++++++++
257#free rolling sphere at midpoint, many triangles in close contact; slowly go through critical points:
258gList = [graphics.Sphere(point=[0,0,0], radius=0.5*r, color= graphics.color.yellow, nTiles=24)]
259omega0 = np.array([-0.05,-5,0.])
260pRef = [0.5*L-1.45*bb, 0.5*L-1.20*bb, 3*bh+0.5*r-2*m*gFact/k] #[0.5*L-1.45*bb, 0.5*L-1.40*bb, ..] goes to edge
261RBinertia = InertiaSphere(m, 0.5*r)
262[nMassStair, oMassStair] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
263                      nodeType=exu.NodeType.RotationRotationVector,
264                      position=pRef,
265                      velocity=-np.cross([0,0,0.5*r],omega0),
266                      rotationMatrix=RotationMatrixX(0.),
267                      angularVelocity=omega0,
268                      gravity=g,
269                      graphicsDataList=gList,
270                      )
271
272nNode3 = nMassStair
273mNode3 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassStair))
274gContact.AddSphereWithMarker(mNode3, radius=0.5*r, contactStiffness=k, contactDamping=20*d, frictionMaterialIndex=0)
275
276if useGraphics:
277    sNode3 = mbs.AddSensor(SensorNode(nodeNumber=nNode3, storeInternal=True, #fileName='solution/contactNode3.txt',
278                                      outputVariableType=exu.OutputVariableType.Displacement))
279evalNodes += [nMassStair]
280
281
282#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
283#contact of cube with ground
284tTrig = 0.25*r #size of contact points on mesh ('thickness')
285gCube = graphics.Brick(size=[3*r,2*r,r], color= graphics.color.steelblue,addNormals=addNormals)
286[meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gCube)
287
288#for tests, 1 refinement!
289[meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
290# exu.Print("n points=",len(meshPoints))
291# [meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
292# exu.Print("==> n points refined=",len(meshPoints))
293# refinements give 26,98,386 points!
294
295[meshPoints2, meshTrigs2] = ShrinkMeshNormalToSurface(meshPoints, meshTrigs, tTrig)
296
297#add mesh to visualization
298gCube = graphics.FromPointsAndTrigs(meshPoints, meshTrigs, color=graphics.color.steelblue) #show refined mesh
299gList = [gCube]
300
301#add points for contact to visualization (shrinked)
302for p in meshPoints2:
303    gList += [graphics.Sphere(point=p, radius=tTrig, color=graphics.color.red)]
304
305pRef = [0.5*L-2*r, 0.25*L, 0.5*r+1.5*tTrig]
306v0 = np.array([-2,0,0])
307RBinertia = InertiaCuboid(density=m/(r*2*r*3*r), sideLengths=[3*r,2*r,r])
308[nMassCube0, oMassCube0] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
309                      nodeType=exu.NodeType.RotationRotationVector,
310                      position=pRef,
311                      velocity=v0,
312                      #rotationMatrix=RotationMatrixZ(0.),
313                      angularVelocity=[0,0,0],
314                      gravity=g,
315                      graphicsDataList=gList,
316                      )
317
318nCube0 = nMassCube0
319mCube0 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassCube0))
320
321
322gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mCube0, contactStiffness=k, contactDamping=d, frictionMaterialIndex=1,
323    pointList=meshPoints,  triangleList=meshTrigs)
324
325for p in meshPoints2:
326    mPoint = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMassCube0, localPosition=p))
327    gContact.AddSphereWithMarker(mPoint, radius=tTrig, contactStiffness=k, contactDamping=d, frictionMaterialIndex=1)
328
329if useGraphics:
330    sCube0 = mbs.AddSensor(SensorNode(nodeNumber=nCube0, storeInternal=True, #fileName='solution/contactCube0.txt',
331                                      outputVariableType=exu.OutputVariableType.Displacement))
332
333evalNodes += [nMassCube0]
334
335
336#%%++++++++++++++++++++++++++++++++++++++++++++
337
338#add as last because of transparency
339oGround = mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=gDataList)))
340
341#%%+++++++++++++++++++++++++++++++++
342mbs.Assemble()
343
344tEnd = 0.8 #tEnd = 0.8 for test suite
345h= 0.0002  #h= 0.0002 for test suite
346# h*=0.1
347# tEnd*=3
348simulationSettings = exu.SimulationSettings()
349#simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
350simulationSettings.solutionSettings.writeSolutionToFile = False
351if useGraphics:
352    simulationSettings.solutionSettings.solutionWritePeriod = 0.001
353    simulationSettings.solutionSettings.writeSolutionToFile = True
354    simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/coordinatesSolution.txt'
355else:
356    simulationSettings.solutionSettings.exportAccelerations = False
357    simulationSettings.solutionSettings.exportVelocities = False
358
359simulationSettings.solutionSettings.sensorsWritePeriod = h*10
360simulationSettings.solutionSettings.outputPrecision = 8 #make files smaller
361simulationSettings.timeIntegration.verboseMode = 1
362
363simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
364simulationSettings.timeIntegration.newton.useModifiedNewton = False
365
366SC.visualizationSettings.general.graphicsUpdateInterval=0.05
367# SC.visualizationSettings.general.drawWorldBasis = True
368SC.visualizationSettings.general.circleTiling=200
369SC.visualizationSettings.general.drawCoordinateSystem=True
370SC.visualizationSettings.loads.show=False
371SC.visualizationSettings.bodies.show=True
372SC.visualizationSettings.markers.show=False
373
374SC.visualizationSettings.nodes.show=True
375SC.visualizationSettings.nodes.showBasis =True
376SC.visualizationSettings.nodes.drawNodesAsPoint = False
377SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
378SC.visualizationSettings.nodes.tiling = 4
379SC.visualizationSettings.openGL.drawFaceNormals = False
380
381SC.visualizationSettings.openGL.multiSampling = 4
382SC.visualizationSettings.openGL.shadow = 0.25
383SC.visualizationSettings.openGL.light0position = [-3,3,10,0]
384
385if useGraphics:
386    SC.visualizationSettings.general.autoFitScene = False
387    exu.StartRenderer()
388    if 'renderState' in exu.sys:
389        SC.SetRenderState(exu.sys['renderState'])
390    mbs.WaitForUserToContinue()
391
392simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
393simulationSettings.timeIntegration.endTime = tEnd
394simulationSettings.timeIntegration.explicitIntegration.computeEndOfStepAccelerations = False #increase performance, accelerations less accurate
395mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
396# mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ODE23)
397
398#compute error:
399uSum=0
400for node in evalNodes:
401    u = mbs.GetNodeOutput(node, exu.OutputVariableType.Coordinates)
402    exu.Print('coords node'+str(node)+' =',u)
403    for c in u:
404        uSum += abs(c) #add up all coordinates for comparison
405
406
407exu.Print('solution of generalContactFrictionTest=',uSum)
408exudynTestGlobals.testError = uSum - (10.132106712933348 )
409
410exudynTestGlobals.testResult = uSum
411
412
413if useGraphics:
414    SC.WaitForRenderEngineStopFlag()
415
416    if True:
417        SC.visualizationSettings.general.autoFitScene = False
418        SC.visualizationSettings.general.graphicsUpdateInterval=0.02
419
420        sol = LoadSolutionFile('solution/coordinatesSolution.txt', safeMode=True)#, maxRows=100)
421        print('start SolutionViewer')
422        mbs.SolutionViewer(sol)
423
424    exu.StopRenderer() #safely close rendering window!
425
426if useGraphics:
427
428
429    mbs.PlotSensor([], closeAll=True)
430    mbs.PlotSensor([sNode3]*3, [0,1,2], figureName='node stair')