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