contactFrictionCircleCable2DnoFrictionTest.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  ObjectContactFrictionCircleCable2D without a friction model - frictionStiffness and
 5#           frictionVelocityPenalty zero - has no tangential force, whatever the friction coefficient and
 6#           whatever slip state the data node starts with (#1290). A cable lies on a fixed circle, in contact
 7#           from the start, its data node initialized to slip; the cable is pulled along its axis. The motion
 8#           with frictionCoefficient 0.5 equals the one with 0, and the tangential force is zero in every step.
 9#
10# Author:   Johannes Gerstmayr
11# Date:     2026-10-05
12#
13# 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.
14#
15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
16
17import exudyn as exu
18from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
19import numpy as np
20
21testIsActive = exu.sys.get('testIsActive', False)
22
23def Run(frictionCoefficient):
24    SC = exu.SystemContainer()
25    mbs = SC.AddSystem()
26    oGround = mbs.AddObject(ObjectGround())
27    length, nElements = 1., 4
28    (b, h, E, rho) = (0.01, 0.01, 2e9, 1000)
29    nodes = [mbs.AddNode(Point2DS1(referenceCoordinates=[length/nElements*i,0,1,0])) for i in range(nElements+1)]
30    cables = [mbs.AddObject(Cable2D(length=length/nElements, massPerLength=rho*b*h, bendingStiffness=E*b*h**3/12,
31                                    axialStiffness=E*b*h, nodeNumbers=[nodes[i], nodes[i+1]])) for i in range(nElements)]
32    #gravity, and a pull along the cable at its end
33    for (i, node) in enumerate(nodes):
34        mbs.AddLoad(Force(markerNumber=mbs.AddMarker(MarkerNodePosition(nodeNumber=node)),
35                          loadVector=[1.*(i == nElements), -9.81*rho*b*h*length/nElements, 0]))
36
37    radius = 0.3
38    mCircle = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0.5,-radius-0.001,0]))
39    nSegments = 4
40    for cable in cables:
41        mCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=cable, numberOfSegments=nSegments))
42        #in contact from the start, and in slip: a state that has no meaning without a friction model
43        data = mbs.AddNode(NodeGenericData(initialCoordinates=[-0.001]*nSegments + [1]*nSegments + [0]*nSegments,
44                                           numberOfDataCoordinates=3*nSegments))
45        mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mCircle, mCable], nodeNumber=data,
46                                                         numberOfContactSegments=nSegments, contactStiffness=1e4,
47                                                         contactDamping=10, frictionStiffness=0, frictionVelocityPenalty=0,
48                                                         frictionCoefficient=frictionCoefficient, circleRadius=radius))
49    mbs.Assemble()
50
51    mbs.variables['tangential'] = 0.
52    def PostStep(mbs, t):
53        for i in range(mbs.systemData.NumberOfObjects()):
54            if mbs.GetObject(i)['objectType'] == 'ContactFrictionCircleCable2D':
55                force = mbs.GetObjectOutput(i, exu.OutputVariableType.ForceLocal)
56                mbs.variables['tangential'] = max(mbs.variables['tangential'], np.max(np.abs(force[0::2])))
57        return True
58    mbs.SetPostStepUserFunction(PostStep)
59
60    simulationSettings = exu.SimulationSettings()
61    simulationSettings.timeIntegration.numberOfSteps = 1000 #100 steps did not resolve the contact: Newton failed and the result depended on the rounding (#2912)
62    simulationSettings.timeIntegration.endTime = 0.1
63    simulationSettings.timeIntegration.verboseMode = 0
64    simulationSettings.solution.file.write = False
65    mbs.SolveDynamic(simulationSettings)
66    return (mbs.systemData.GetODE2Coordinates(), mbs.variables['tangential'])
67
68(q0, tangential0) = Run(0.)
69(q1, tangential1) = Run(0.5)
70difference = np.linalg.norm(q1 - q0)
71exu.Print('contact without friction: tangential forces', tangential0, tangential1, ', difference of the motion', difference)
72
73u = np.sum(q0) + 1e3*(difference + tangential0 + tangential1)
74exu.Print('solution of contactFrictionCircleCable2DnoFrictionTest=', u)
75
76exu.sys['testResult'] = u