contactSphereTorusMomentumTest.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  Momentum conservation of ObjectContactSphereTorus (#2127): a free ball in the groove of a free
 5#           ring, both rigid bodies, no gravity; the ball hits the groove and the ring spins. The contact
 6#           forces are internal, so the total linear and angular momentum of the two bodies are conserved,
 7#           with and without friction. Measured: the angular momentum about the origin stays at its initial
 8#           value to 3e-9 (without friction) and 2e-9 (with friction), the linear momentum to 1e-15 - the
 9#           torque of the normal force on the ring, which acts at the contact point and not at the center of
10#           the ring, was missing without friction (angular momentum off by 6e-3 of 0.09).
11#
12# Author:   Johannes Gerstmayr
13# Date:     2026-09-30
14#
15# 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.
16#
17#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
18
19import exudyn as exu
20from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
21
22import numpy as np
23
24testIsActive = exu.sys.get('testIsActive', False)
25
26def Run(mu, stepSize=1e-5, tEnd=0.05):
27    """the ball and the ring; returns the times, the total angular and linear momentum"""
28    SC=exu.SystemContainer(); mbs=SC.AddSystem()
29    iB=InertiaSphere(mass=0.1, radius=0.01)
30    iR=InertiaCylinder(density=7800, length=0.02, outerRadius=0.12, innerRadius=0.1, axis=2)
31    nB=mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[0.1005,0,0]+eulerParameters0, initialVelocities=[0.5,1,0.2]+list(AngularVelocity2EulerParameters_t([0,0,50],eulerParameters0))))
32    nR=mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[0,0,0]+eulerParameters0, initialVelocities=[0,0,0]+list(AngularVelocity2EulerParameters_t([0,2,3],eulerParameters0))))
33    b=mbs.AddObject(ObjectRigidBody(nodeNumber=nB, mass=iB.Mass(), inertia=iB.GetInertia6D()))
34    r=mbs.AddObject(ObjectRigidBody(nodeNumber=nR, mass=iR.Mass(), inertia=iR.GetInertia6D(), centerOfMass=iR.COM()))
35    mB=mbs.AddMarker(MarkerBodyRigid(bodyNumber=b)); mR=mbs.AddMarker(MarkerBodyRigid(bodyNumber=r))
36    nD=mbs.AddNode(NodeGenericData(numberOfDataCoordinates=4, initialCoordinates=[0,0,0,0]))
37    mbs.AddObject(ObjectContactSphereTorus(markerNumbers=[mB,mR], nodeNumber=nD, sphereRadius=0.01, torusMajorRadius=0.1, torusMinorRadius=0.011,
38        torusAxis=[0,0,1], contactStiffness=1e6, contactDamping=1e2, dynamicFriction=mu, frictionProportionalZone=1e-3))
39    sens={}
40    for name,n in [('B',nB),('R',nR)]:
41        for ov in ['Position','Velocity','AngularVelocity','RotationMatrix']:
42            sens[name+ov]=mbs.AddSensor(SensorNode(nodeNumber=n,outputVariableType=getattr(exu.OutputVariableType,ov),storeInternal=True))
43    mbs.Assemble()
44    s=exu.SimulationSettings(); s.timeIntegration.numberOfSteps=int(tEnd/stepSize); s.timeIntegration.endTime=tEnd
45    s.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
46    s.solution.file.write=False; s.solution.sensors.writePeriod=stepSize*10
47    s.timeIntegration.generalizedAlpha.spectralRadius=1
48    mbs.SolveDynamic(s)
49    D={k:mbs.GetSensorStoredData(v)[:,1:] for k,v in sens.items()}
50    t=mbs.GetSensorStoredData(sens['BPosition'])[:,0]
51    def L(name, inertia, mass, com):
52        J=inertia.Inertia()
53        out=[]; P=[]
54        for k in range(len(t)):
55            A=D[name+'RotationMatrix'][k].reshape(3,3); w=D[name+'AngularVelocity'][k]
56            p=D[name+'Position'][k]+A@com; v=D[name+'Velocity'][k]+np.cross(w,A@com)
57            #inertia about COM
58            Jc=J - mass*(np.dot(com,com)*np.eye(3)-np.outer(com,com))
59            out.append(np.cross(p,mass*v)+A@Jc@A.T@w); P.append(mass*v)
60        return np.array(out), np.array(P)
61    LB,PB=L('B',iB,iB.Mass(),np.zeros(3)); LR,PR=L('R',iR,iR.Mass(),np.array(iR.COM()))
62    Lt=LB+LR; Pt=PB+PR
63    return t, Lt, Pt, np.abs(LB).max()
64
65
66testResult = 0
67for mu in [0, 0.3]:
68    t, L, P, _ = Run(mu)
69    errorL = np.abs(L-L[0]).max()
70    errorP = np.abs(P-P[0]).max()
71    exu.Print('friction', mu, ': angular momentum', L[-1].round(8), ', largest change', errorL, ', of the linear momentum', errorP)
72    testResult += sum(L[-1]) + (errorL < 1e-7) + (errorP < 1e-12)
73
74exu.Print('solution of contactSphereTorusMomentumTest=', testResult)
75exu.sys['testResult'] = testResult