contactCurveCirclesFrictionTest.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  ObjectContactCurveCircles with friction and its output variables (#2867): a circle (a planar rigid body)
 5#           slides along a straight curve under gravity with an initial velocity; the friction force mu*|f_N|
 6#           decelerates it by mu*g, so after the time t its velocity is v0 - mu*g*t. The output variables per
 7#           segment: DisplacementLocal the gap (the static penetration m*g/(k*l) for contactModel 0, with the
 8#           length l of the segment), VelocityLocal its rate, ForceLocal the force on the circle, [-mu*m*g, m*g]
 9#           while it slides.
10#
11# Author:   Johannes Gerstmayr
12# Date:     2026-10-07
13#
14# 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.
15#
16#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
17
18import exudyn as exu
19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
20import numpy as np
21
22testIsActive = exu.sys.get('testIsActive', False)
23
24mass, g, mu, v0 = 1., 10., 0.2, 1.
25k, d, lSegment, radius = 1e5, 400., 4., 0.1
26tEnd = 0.3
27
28SC = exu.SystemContainer()
29mbs = SC.AddSystem()
30oGround = mbs.AddObject(ObjectGround())
31gapStatic = mass*g/(k*lSegment)
32node = mbs.AddNode(NodeRigidBody2D(referenceCoordinates=[0, radius - gapStatic, 0], initialVelocities=[v0, 0, 0]))
33oBody = mbs.AddObject(ObjectRigidBody2D(nodeNumber=node, mass=mass, inertia=0.01))
34mCircle = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oBody))
35mbs.AddLoad(LoadForceVector(markerNumber=mCircle, loadVector=[0, -mass*g, 0]))
36mCurve = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround))
37segments = np.array([[2,0, -2,0]]) #one segment of length 4; the contact side is y > 0
38nData = mbs.AddNode(NodeGenericData(numberOfDataCoordinates=3, initialCoordinates=[-1,0,0]))
39oContact = mbs.AddObject(ObjectContactCurveCircles(markerNumbers=[mCurve, mCircle], nodeNumber=nData, circlesRadii=[radius],
40                                                   segmentsData=exu.MatrixContainer(segments),
41                                                   contactStiffness=k, contactDamping=d,
42                                                   dynamicFriction=mu, frictionProportionalZone=1e-4))
43mbs.Assemble()
44
45simulationSettings = exu.SimulationSettings()
46simulationSettings.timeIntegration.numberOfSteps = 3000
47simulationSettings.timeIntegration.endTime = tEnd
48simulationSettings.timeIntegration.verboseMode = 0
49simulationSettings.solution.file.write = False
50mbs.SolveDynamic(simulationSettings)
51
52velocity = mbs.GetNodeOutput(node, exu.OutputVariableType.Velocity)
53gap = np.atleast_1d(mbs.GetObjectOutput(oContact, exu.OutputVariableType.DisplacementLocal)) #one value per segment
54gap_t = np.atleast_1d(mbs.GetObjectOutput(oContact, exu.OutputVariableType.VelocityLocal))
55force = mbs.GetObjectOutput(oContact, exu.OutputVariableType.ForceLocal)
56exu.Print('velocity', velocity[0], ', expected', v0 - mu*g*tEnd)
57exu.Print('gap', gap, ', expected', -gapStatic, '; gap_t', gap_t)
58exu.Print('force on the circle', force, ', expected', [-mu*mass*g, mass*g])
59
60errors = (int(abs(velocity[0] - (v0 - mu*g*tEnd)) > 1e-3) + int(abs(gap[0] + gapStatic) > 1e-6)
61          + int(abs(gap_t[0]) > 1e-3) + int(np.abs(force - [-mu*mass*g, mass*g]).max() > 1e-2))
62exu.Print('contactCurveCirclesFrictionTest: errors', errors)
63u = errors + velocity[0] + force[0] + force[1]
64exu.Print('solution of contactCurveCirclesFrictionTest=', u)
65exu.sys['testResult'] = u