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