generalContactTriangleMomentumTest.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  Momentum conservation of the sphere-triangle contact of GeneralContact (#2849): a free ball hits
 5#           a free plate off its center, both rigid bodies, no gravity; the plate is a brick meshed into
 6#           triangles. The contact forces are internal, so the total linear and angular momentum of the two
 7#           bodies are conserved, with and without friction. Measured: the angular momentum about the
 8#           origin stays at its initial value to 3e-9 (without friction) and 8e-10 (with friction), the
 9#           linear momentum to 1e-16 - the torque of the normal force on the plate, which acts at the
10#           contact point and not at the reference point of the plate, was missing without friction
11#           (angular momentum off by 7e-3 of 0.015).
12#
13# Author:   Johannes Gerstmayr
14# Date:     2026-10-05
15#
16# 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.
17#
18#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
19
20import exudyn as exu
21from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
22import exudyn.graphics as graphics
23
24import numpy as np
25
26testIsActive = exu.sys.get('testIsActive', False)
27
28def Run(mu, stepSize=1e-5, tEnd=0.02):
29    """the ball and the plate; returns the total angular and linear momentum over time"""
30    SC = exu.SystemContainer()
31    mbs = SC.AddSystem()
32    rBall = 0.01
33    iBall = InertiaSphere(mass=0.1, radius=rBall)
34    plateSize = [0.2, 0.2, 0.02]
35    iPlate = InertiaCuboid(density=1000, sideLengths=plateSize)
36    nBall = mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[0.07, 0.03, 0.0205] + eulerParameters0,
37                                        initialVelocities=[0.3, 0.1, -1] + list(AngularVelocity2EulerParameters_t([0, 0, 20], eulerParameters0))))
38    nPlate = mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[0, 0, 0] + eulerParameters0,
39                                         initialVelocities=[0, 0, 0] + list(AngularVelocity2EulerParameters_t([0, 1, 2], eulerParameters0))))
40    ball = mbs.AddObject(ObjectRigidBody(nodeNumber=nBall, mass=iBall.Mass(), inertia=iBall.GetInertia6D()))
41    plate = mbs.AddObject(ObjectRigidBody(nodeNumber=nPlate, mass=iPlate.Mass(), inertia=iPlate.GetInertia6D()))
42
43    gContact = mbs.AddGeneralContact()
44    gContact.frictionProportionalZone = 1e-3
45    gContact.SetFrictionPairings(mu*np.eye(1))
46    gContact.SetSearchTreeCellSize(numberOfCells=[4, 4, 4])
47    [points, triangles] = graphics.ToPointsAndTrigs(graphics.Brick([0, 0, 0], plateSize))
48    gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mbs.AddMarker(MarkerBodyRigid(bodyNumber=plate)),
49                                        contactStiffness=2e6, contactDamping=2e2, frictionMaterialIndex=0,
50                                        pointList=points, triangleList=triangles)
51    gContact.AddSphereWithMarker(mbs.AddMarker(MarkerNodeRigid(nodeNumber=nBall)), radius=rBall,
52                                 contactStiffness=2e6, contactDamping=2e2, frictionMaterialIndex=0)
53
54    sensors = {}
55    for (name, node) in [('ball', nBall), ('plate', nPlate)]:
56        for variable in ['Position', 'Velocity', 'AngularVelocity', 'RotationMatrix']:
57            sensors[name+variable] = mbs.AddSensor(SensorNode(nodeNumber=node, storeInternal=True,
58                                                              outputVariableType=getattr(exu.OutputVariableType, variable)))
59    mbs.Assemble()
60    simulationSettings = exu.SimulationSettings()
61    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
62    simulationSettings.timeIntegration.endTime = tEnd
63    simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
64    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 1
65    simulationSettings.solution.file.write = False
66    simulationSettings.solution.sensors.writePeriod = stepSize*10
67    mbs.SolveDynamic(simulationSettings)
68
69    data = dict((key, mbs.GetSensorStoredData(sensor)[:, 1:]) for (key, sensor) in sensors.items())
70    def Momentum(name, inertia):
71        J = inertia.Inertia()
72        m = inertia.Mass()
73        angular = []
74        linear = []
75        for k in range(len(data[name+'Position'])):
76            A = data[name+'RotationMatrix'][k].reshape(3, 3)
77            v = data[name+'Velocity'][k]
78            angular.append(np.cross(data[name+'Position'][k], m*v) + A @ J @ A.T @ data[name+'AngularVelocity'][k])
79            linear.append(m*v)
80        return np.array(angular), np.array(linear)
81    (angularBall, linearBall) = Momentum('ball', iBall)
82    (angularPlate, linearPlate) = Momentum('plate', iPlate)
83    return angularBall + angularPlate, linearBall + linearPlate
84
85
86testResult = 0
87for mu in [0, 0.3]:
88    (L, P) = Run(mu)
89    errorL = np.abs(L-L[0]).max()
90    errorP = np.abs(P-P[0]).max()
91    exu.Print('friction', mu, ': angular momentum', L[-1].round(8), ', largest change', errorL, ', of the linear momentum', errorP)
92    testResult += sum(L[-1]) + (errorL < 1e-7) + (errorP < 1e-12)
93
94exu.Print('solution of generalContactTriangleMomentumTest=', testResult)
95exu.sys['testResult'] = testResult