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)
122dictMass = mbs.CreateRigidBody(
123 inertia=RBinertia,
124 nodeType=exu.NodeType.RotationRotationVector,
125 referencePosition=pRef,
126 initialVelocity=-np.cross([0,0,r], omega0),
127 referenceRotationMatrix=RotationMatrixX(0.),
128 initialAngularVelocity=omega0,
129 # gravity=g,
130 graphicsDataList=gList,
131 returnDict=True)
132[nMass, oMass] = [dictMass['nodeNumber'], dictMass['bodyNumber']]
133
134
135nNode0 = nMass
136mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
137mbs.AddLoad(Force(markerNumber=mNode, loadVector= [0,0,-k*r*0.01])) #==> uz = 2*r*0.01
138exu.Print('expect u0z=',2*r*0.01)
139gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
140if useGraphics:
141 sNode0 = mbs.AddSensor(SensorNode(nodeNumber=nNode0, storeInternal=True, #fileName='solution/contactNode0.txt',
142 outputVariableType=exu.OutputVariableType.Displacement))
143 vNode0 = mbs.AddSensor(SensorNode(nodeNumber=nNode0, storeInternal=True, #fileName='solution/contactNode0Vel.txt',
144 outputVariableType=exu.OutputVariableType.Velocity))
145evalNodes += [nMass]
146
147#%%++++++++++++++++++++++++++++++++++++++++++++
148#free rolling sphere at midpoint, many triangles in close contact; slowly go through critical points:
149gList = [graphics.Sphere(point=[0,0,0], radius=r, color= graphics.color.yellow, nTiles=24)]
150omega0 = -1e-12*np.array([1,0.1,0.])
151pRef = [1e-15,-1e-14,r-2*m*gFact/k]
152RBinertia = InertiaSphere(m, r)
153dictMass = mbs.CreateRigidBody(
154 inertia=RBinertia,
155 nodeType=exu.NodeType.RotationRotationVector,
156 referencePosition=pRef,
157 initialVelocity=-np.cross([0,0,r], omega0),
158 referenceRotationMatrix=RotationMatrixX(0.),
159 initialAngularVelocity=omega0,
160 gravity=g,
161 graphicsDataList=gList,
162 returnDict=True)
163[nMass, oMass] = [dictMass['nodeNumber'], dictMass['bodyNumber']]
164
165nNode1 = nMass
166mNode1 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
167#mbs.AddLoad(Force(markerNumber=mNode1, loadVector= [0,0,-k*r*0.01])) #==> uz = 2*r*0.01
168#exu.Print('expect u0z=',2*r*0.01)
169gContact.AddSphereWithMarker(mNode1, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
170
171sNode1 = mbs.AddSensor(SensorNode(nodeNumber=nNode1, storeInternal=True, #fileName='solution/contactNode1.txt',
172 outputVariableType=exu.OutputVariableType.Displacement))
173# vNode1 = mbs.AddSensor(SensorNode(nodeNumber=nNode0, storeInternal=True, #fileName='solution/contactNode0Vel.txt',
174# outputVariableType=exu.OutputVariableType.Velocity))
175
176
177#%%++++++++++++++++++++++++++++++++++++++++++++
178#fixed pressure tests:
179pf = np.array([-1.2*L,0,0])
180nGroundF = mbs.AddNode(NodePointGround(referenceCoordinates=pf ))
181mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGroundF))
182gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
183gDataList += [graphics.Sphere(point=pf, radius=r, color= graphics.color.grey, nTiles=24)]
184
185gList = [graphics.Sphere(point=[0,0,0], radius=r, color= graphics.color.lightgreen, nTiles=24)]
186
187pRef = pf+[0,2*r,0] #[-0.4*L,-0.4*L,r-m*gFact/k]
188RBinertia = InertiaSphere(m, r)
189dictF = mbs.CreateRigidBody(
190 inertia=RBinertia,
191 nodeType=exu.NodeType.RotationRotationVector,
192 referencePosition=pRef,
193 referenceRotationMatrix=RotationMatrixX(0.),
194 graphicsDataList=gList,
195 returnDict=True)
196[nMassF, oMassF] = [dictF['nodeNumber'], dictF['bodyNumber']]
197
198mC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMassF, coordinate=0))
199mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mC]))
200
201mNodeF = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassF))
202mbs.AddLoad(Force(markerNumber=mNodeF, loadVector= [0,-k*r*0.1,0])) #==> u = k*r*0.1/(0.5*k) = 2*r*0.1
203exu.Print('expect uFy=',2*r*0.1)
204gContact.AddSphereWithMarker(mNodeF, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
205if useGraphics:
206 sNodeF = mbs.AddSensor(SensorNode(nodeNumber=nMassF, storeInternal=True, #fileName='solution/contactNodeF.txt',
207 outputVariableType=exu.OutputVariableType.Displacement))
208evalNodes += [nMassF]
209
210#%%++++++++++++++++++++++++++++++++++++++++++++
211# sliding between spheres:
212pr = np.array([-1.2*L,0.5*L,0])
213nGroundF2 = mbs.AddNode(NodePointGround(referenceCoordinates=pr ))
214mNode2 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGroundF2))
215gContact.AddSphereWithMarker(mNode2, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
216gDataList += [graphics.Sphere(point=pr, radius=r, color= graphics.color.lightgrey, nTiles=24)]
217
218gList = [graphics.Sphere(point=[0,0,0], radius=r, color= graphics.color.lightred, nTiles=24)]
219
220dRol = r*0.01
221pRef = pr+[0,2*r-2*dRol,0] #force=k*r*0.01
222fRol = k*dRol
223exu.Print('force rolling=', fRol, ', torque=', fRol*mu*r)
224RBinertia = InertiaSphere(m, r)
225dictR = mbs.CreateRigidBody(
226 inertia=RBinertia,
227 nodeType=exu.NodeType.RotationRotationVector,
228 referencePosition=pRef,
229 referenceRotationMatrix=RotationMatrixX(0.),
230 graphicsDataList=gList,
231 returnDict=True)
232[nMassR, oMassR] = [dictR['nodeNumber'], dictR['bodyNumber']]
233
234
235
236
237mC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMassR, coordinate=0))
238mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mC]))
239mC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMassR, coordinate=1))
240mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mC]))
241mC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMassR, coordinate=2))
242mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mC]))
243
244mNodeR = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassR))
245
246def UFtorque(mbs, t, loadVector):
247 torque = 10*t*fRol*mu*r
248 if t > 0.3:
249 torque = 0
250 return [torque,0,0]
251mbs.AddLoad(Torque(markerNumber=mNodeR, loadVectorUserFunction=UFtorque,
252 loadVector= [1,0,0])) #==> u = k*r*0.1/(0.5*k) = 2*r*0.1
253
254gContact.AddSphereWithMarker(mNodeR, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
255if useGraphics:
256 sNodeR = mbs.AddSensor(SensorNode(nodeNumber=nMassR, storeInternal=True, #fileName='solution/contactNodeR.txt',
257 outputVariableType=exu.OutputVariableType.Rotation))
258 vNodeR = mbs.AddSensor(SensorNode(nodeNumber=nMassR, storeInternal=True, #fileName='solution/contactNodeRvel.txt',
259 outputVariableType=exu.OutputVariableType.AngularVelocity))
260
261evalNodes += [nMassR]
262
263#%%++++++++++++++++++++++++++++++++++++++++++++
264#sphere on stairs
265#%%++++++++++++++++++++++++++++++++++++++++++++
266#free rolling sphere at midpoint, many triangles in close contact; slowly go through critical points:
267gList = [graphics.Sphere(point=[0,0,0], radius=0.5*r, color= graphics.color.yellow, nTiles=24)]
268omega0 = np.array([-0.05,-5,0.])
269pRef = [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
270RBinertia = InertiaSphere(m, 0.5*r)
271dictStair = mbs.CreateRigidBody( #note that this node causes inaccuracies/repeatability issues in test suite!
272 inertia=RBinertia,
273 nodeType=exu.NodeType.RotationRotationVector,
274 referencePosition=pRef,
275 initialVelocity=-np.cross([0,0,0.5*r], omega0),
276 referenceRotationMatrix=RotationMatrixX(0.),
277 initialAngularVelocity=omega0,
278 gravity=g,
279 graphicsDataList=gList,
280 returnDict=True)
281[nMassStair, oMassStair] = [dictStair['nodeNumber'], dictStair['bodyNumber']]
282
283nNode3 = nMassStair
284mNode3 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassStair))
285gContact.AddSphereWithMarker(mNode3, radius=0.5*r, contactStiffness=k, contactDamping=20*d, frictionMaterialIndex=0)
286
287if useGraphics:
288 sNode3 = mbs.AddSensor(SensorNode(nodeNumber=nNode3, storeInternal=True, #fileName='solution/contactNode3.txt',
289 outputVariableType=exu.OutputVariableType.Displacement))
290evalNodes += [nMassStair]
291
292
293#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
294#contact of cube with ground
295tTrig = 0.25*r #size of contact points on mesh ('thickness')
296gCube = graphics.Brick(size=[3*r,2*r,r], color= graphics.color.steelblue,addNormals=addNormals)
297[meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gCube)
298
299#for tests, 1 refinement!
300[meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
301# exu.Print("n points=",len(meshPoints))
302# [meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
303# exu.Print("==> n points refined=",len(meshPoints))
304# refinements give 26,98,386 points!
305
306[meshPoints2, meshTrigs2] = ShrinkMeshNormalToSurface(meshPoints, meshTrigs, tTrig)
307
308#add mesh to visualization
309gCube = graphics.FromPointsAndTrigs(meshPoints, meshTrigs, color=graphics.color.steelblue) #show refined mesh
310gList = [gCube]
311
312#add points for contact to visualization (shrinked)
313for p in meshPoints2:
314 gList += [graphics.Sphere(point=p, radius=tTrig, color=graphics.color.red)]
315
316pRef = [0.5*L-2*r, 0.25*L, 0.5*r+1.5*tTrig]
317v0 = np.array([-2,0,0])
318RBinertia = InertiaCuboid(density=m/(r*2*r*3*r), sideLengths=[3*r,2*r,r])
319dictCube0 = mbs.CreateRigidBody(
320 inertia=RBinertia,
321 nodeType=exu.NodeType.RotationRotationVector,
322 referencePosition=pRef,
323 initialVelocity=v0,
324 initialAngularVelocity=[0,0,0],
325 gravity=g,
326 graphicsDataList=gList,
327 returnDict=True)
328[nMassCube0, oMassCube0] = [dictCube0['nodeNumber'], dictCube0['bodyNumber']]
329
330nCube0 = nMassCube0
331mCube0 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassCube0))
332
333
334gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mCube0, contactStiffness=k, contactDamping=d, frictionMaterialIndex=1,
335 pointList=meshPoints, triangleList=meshTrigs)
336
337for p in meshPoints2:
338 mPoint = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMassCube0, localPosition=p))
339 gContact.AddSphereWithMarker(mPoint, radius=tTrig, contactStiffness=k, contactDamping=d, frictionMaterialIndex=1)
340
341if useGraphics:
342 sCube0 = mbs.AddSensor(SensorNode(nodeNumber=nCube0, storeInternal=True, #fileName='solution/contactCube0.txt',
343 outputVariableType=exu.OutputVariableType.Displacement))
344
345evalNodes += [nMassCube0]
346
347
348#%%++++++++++++++++++++++++++++++++++++++++++++
349
350#add as last because of transparency
351oGround = mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=gDataList)))
352
353#%%+++++++++++++++++++++++++++++++++
354mbs.Assemble()
355
356tEnd = 0.8 #tEnd = 0.8 for test suite
357h= 0.0002 #h= 0.0002 for test suite
358# h*=0.1
359# tEnd*=3
360simulationSettings = exu.SimulationSettings()
361#simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
362simulationSettings.solutionSettings.writeSolutionToFile = False
363if useGraphics:
364 simulationSettings.solutionSettings.solutionWritePeriod = 0.001
365 simulationSettings.solutionSettings.writeSolutionToFile = True
366 simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/coordinatesSolution.txt'
367else:
368 simulationSettings.solutionSettings.exportAccelerations = False
369 simulationSettings.solutionSettings.exportVelocities = False
370
371simulationSettings.solutionSettings.sensorsWritePeriod = h*10
372simulationSettings.solutionSettings.outputPrecision = 8 #make files smaller
373simulationSettings.timeIntegration.verboseMode = 1
374
375simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
376simulationSettings.timeIntegration.newton.useModifiedNewton = False
377
378SC.visualizationSettings.general.graphicsUpdateInterval=0.05
379# SC.visualizationSettings.general.drawWorldBasis = True
380SC.visualizationSettings.general.circleTiling=200
381SC.visualizationSettings.general.drawCoordinateSystem=True
382SC.visualizationSettings.loads.show=False
383SC.visualizationSettings.bodies.show=True
384SC.visualizationSettings.markers.show=False
385
386SC.visualizationSettings.nodes.show=True
387SC.visualizationSettings.nodes.showBasis =True
388SC.visualizationSettings.nodes.drawNodesAsPoint = False
389SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
390SC.visualizationSettings.nodes.tiling = 4
391SC.visualizationSettings.openGL.drawFaceNormals = False
392
393SC.visualizationSettings.openGL.multiSampling = 4
394SC.visualizationSettings.openGL.shadow = 0.25
395SC.visualizationSettings.openGL.light0position = [-3,3,10,0]
396
397if useGraphics:
398 SC.visualizationSettings.general.autoFitScene = False
399 SC.renderer.Start()
400 if 'renderState' in exu.sys:
401 SC.renderer.SetState(exu.sys['renderState'])
402 SC.renderer.DoIdleTasks()
403
404simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
405simulationSettings.timeIntegration.endTime = tEnd
406simulationSettings.timeIntegration.explicitIntegration.computeEndOfStepAccelerations = False #increase performance, accelerations less accurate
407mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
408# mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ODE23)
409
410#compute error:
411uSum=0
412for node in evalNodes:
413 u = mbs.GetNodeOutput(node, exu.OutputVariableType.Coordinates)
414 exu.Print('coords node'+str(node)+' =',u)
415 for c in u:
416 uSum += abs(c) #add up all coordinates for comparison
417
418
419exu.Print('solution of generalContactFrictionTest=',uSum)
420exudynTestGlobals.testError = uSum - (10.132106712933348 )
421
422exudynTestGlobals.testResult = uSum
423
424
425if useGraphics:
426 SC.renderer.DoIdleTasks()
427
428 if True:
429 SC.visualizationSettings.general.autoFitScene = False
430 SC.visualizationSettings.general.graphicsUpdateInterval=0.02
431
432 sol = LoadSolutionFile('solution/coordinatesSolution.txt', safeMode=True)#, maxRows=100)
433 exu.Print('start SolutionViewer')
434 mbs.SolutionViewer(sol)
435
436 SC.renderer.Stop() #safely close rendering window!
437
438if useGraphics:
439
440
441 mbs.PlotSensor([], closeAll=True)
442 mbs.PlotSensor([sNode3]*3, [0,1,2], figureName='node stair')