<!-- GENERATED by tools/generators/examplesDocsEmitter.py from python/ - do not edit -->
(testmodels-ancfbeltdrive)=
# ANCFbeltDrive.py

You can view and download this file on Github: [ANCFbeltDrive.py](https://github.com/jgerstmayr/EXUDYN/tree/master/python/TestModels/ANCFbeltDrive.py)

```{code-block} python
:linenos:

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  ANCF cable elements in contact with circles defined by GeneralContact;
#           this is an advanced example for a belt drive according to literature; 
#           note that the belt is very thin in order to avoid bending effects, to make this comparible to analytical results
#
# Author:   Johannes Gerstmayr
# Date:     2022-01-31
#
# 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.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
import exudyn as exu
from exudyn.itemInterface import *
import numpy as np
from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
import exudyn.graphics as graphics
from exudyn.beams import *
from math import sin, cos, sqrt, pi

testIsActive = exu.sys.get('testIsActive', False)

SC = exu.SystemContainer()
mbs = SC.AddSystem()

#background
rect = [-1,-1.5,3,1.5] #xmin,ymin,xmax,ymax
background0 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[rect[0],rect[1],0, rect[2],rect[1],0, rect[2],rect[3],0, rect[0],rect[3],0, rect[0],rect[1],0]} #background
oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], 
                                   visualization=VObjectGround(graphicsData= [background0])))
nGround = mbs.AddNode(NodePointGround())
mCoordinateGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))


#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#contact
doImplicit = True
useContact = True
useFriction = True
dryFriction = 0.5
contactStiffness = 1e5
contactDamping = 1e-3*contactStiffness

if useContact:
    gContact = mbs.AddGeneralContact()
    gContact.verboseMode = 1*0
    gContact.frictionProportionalZone = 1
    gContact.ancfCableUseExactMethod = False
    gContact.ancfCableNumberOfContactSegments = 8
    ssx = 16#32 #search tree size
    ssy = 8#32 #search tree size
    ssz = 1 #search tree size
    gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssz])
    #gContact.SetSearchTreeBox(pMin=np.array([-1,-1,-1]), pMax=np.array([4,1,1]))
torque=10*2
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#wheel:
dampingWheel1 = 1  #5 #add breaking torque, to limit velocity
#cable:

numberOfElements = 16   # per section; 32 recommended
curvedRefConf=False     # this flag could initialize the elements to be produced curved -> not suitable for belt drive!
L=2                     # length of ANCF element in m
E=1e10                  # Young's modulus of ANCF element in N/m^2
rhoBeam=1000            # density of ANCF element in kg/m^3
b=0.002                 # width of rectangular ANCF element in m
h=0.002                 # height of rectangular ANCF element in m
A=b*h                   # cross sectional area of ANCF element in m^2
I=b*h**3/12             # second moment of area of ANCF element in m^4
dEI = 0*1e-3*E*I
dEA = 1e-2*E*A
# f=3*E*I/L**2            # tip load applied to ANCF element in N
g=-9.81
dimZ = b #z.dimension
preStretch=-0.002
# exu.Print("load f="+str(f))
# exu.Print("EI="+str(E*I))

# nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
# mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action

cableTemplate = Cable2D(#length = L / nElements, #set in GenerateBeamElementsAlongLine(...)
                        massPerLength = rhoBeam*A,
                        bendingStiffness = E*I,
                        axialStiffness = E*A,
                        bendingDamping = dEI,
                        axialDamping = dEA,
                        referenceAxialStrain = preStretch, #prestretch
                        #nodeNumbers = [0, 0], #will be filled in GenerateBeamElementsAlongLine(...)
                        visualization=VCable2D(drawHeight=2*h),
                        )
exu.Print("pre-stretch force=", preStretch*E*A)
exu.Print("beam mass per length=", rhoBeam*A)
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#create belt drive:
distanceWheels = 2 #distance of wheel centers
wheelCenter0 = np.array([0,0,0])
wheelCenter1 = np.array([distanceWheels,0,0])

rWheel0 = 0.5
rWheel1 = 0.5
rRoll = 0.1
rollCenter1 = np.array([0.5*distanceWheels,-rWheel0+rRoll,0])

mWheel = 4
mRoll = 1e-6 #small as compared to beam with rhoA=4e-3
kRoll = 1e-2
dRelRoll = 0.005
fRoll = 0.01

velocityControl = True

yAxis = np.array([0,1.,0])
ancfList=[]

if True:
    startAngle = -pi
    arcAngle = -pi
    positionOfNode0 = wheelCenter0-rWheel0*yAxis # starting point of line
    ancf=GenerateCircularArcANCFCable2D(mbs, positionOfNode0, 
                                        rWheel0, startAngle, arcAngle, numberOfElements, 
                                        cableTemplate,
                                        massProportionalLoad = [0,g,0], #optionally add gravity
                                        #fixedConstraintsNode0 = [1,1,1,1], #add constraints for pos and rot (r'_y)
                                        #fixedConstraintsNode1 = [1,1,1,1],
                                        setCurvedReferenceConfiguration=curvedRefConf, 
                                        #verboseMode=True
                                        )
    ancfList+=[ancf]
    beamInfo=GenerateBeamElementsAlongLine(mbs,
                                         ancf[3][-1], wheelCenter1+rWheel1*yAxis,
                                         numberOfElements,
                                         beamTemplate=cableTemplate, #this defines the beam element properties
                                         gravity = [0,g,0], #optionally add gravity
                                          # groundConstraintsStart = [1,1,1,1], #add constraints for pos and rot (r'_y)
                                          # groundConstraintsEnd = [1,1,1,1],
                                         nodeNumberStart=ancf[0][-1]
                                         )
    ancf = [beamInfo['nodes'], beamInfo['elements'], beamInfo['loads'], beamInfo['nodePositions'], beamInfo['constraintObjects']] #as the list of GenerateCircularArcANCFCable2D
    ancfList+=[ancf]

    if True:
        startAngle = 0
        arcAngle = -pi
        ancf=GenerateCircularArcANCFCable2D(mbs, ancf[3][-1], 
                                            rWheel1, startAngle, arcAngle, numberOfElements, 
                                            cableTemplate,
                                            massProportionalLoad = [0,g,0], #optionally add gravity
                                            #fixedConstraintsNode0 = [1,1,1,1], #add constraints for pos and rot (r'_y)
                                            #fixedConstraintsNode1 = [1,1,1,1],
                                            setCurvedReferenceConfiguration=curvedRefConf, 
                                            nodeNumber0=ancf[0][-1]
                                            )
        ancfList+=[ancf]
        beamInfo=GenerateBeamElementsAlongLine(mbs,
                                             ancf[3][-1], ancfList[0][3][0],
                                             numberOfElements,
                                             beamTemplate=cableTemplate, #this defines the beam element properties
                                             gravity = [0,g,0], #optionally add gravity
                                             nodeNumberStart=ancf[0][-1],
                                             nodeNumberEnd=ancfList[0][0][0]
                                             )
        ancf = [beamInfo['nodes'], beamInfo['elements'], beamInfo['loads'], beamInfo['nodePositions'], beamInfo['constraintObjects']] #as the list of GenerateCircularArcANCFCable2D
        ancfList+=[ancf]


def UFvelocityDrive(mbs, t, itemNumber, lOffset): #time derivative of UFoffset
    vMax = 10 #5m/s
    if 2*t < vMax:
        return 2*t
    return vMax

#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
sAngVel=[]
#add contact:
if useContact:


    halfHeight = 0.5*h*0
    wheels = [{'center':wheelCenter0, 'radius':rWheel0-halfHeight, 'mass':mWheel}, 
              {'center':wheelCenter1, 'radius':rWheel1-halfHeight, 'mass':mWheel}, 
              {'center':rollCenter1, 'radius':rRoll-halfHeight, 'mass':mRoll}, #small mass for roll, not to influence beam
              ]

    for i, wheel in enumerate(wheels):
        r = wheel['radius']
        p = wheel['center']
        mass = wheel['mass']
        rot0 = 0 #initial rotation
        pRef = [p[0], p[1], rot0]
        gList = [graphics.Cylinder(vAxis=[0,0,dimZ], radius=r*0.99, #draw smaller to see cable element 
                                      color= graphics.color.dodgerblue, nTiles=32),
                 graphics.Arrow(pAxis=[0,0,0.02*r], vAxis=[r,0,0], radius=0.02*r, color=graphics.color.orange)]

        omega0 = 0 #initial angular velocity
        v0 = np.array([0,0,omega0]) 

        RBinertia = InertiaCylinder(mass/(r**2*np.pi*b), b, r, axis=2)

        nMass = mbs.AddNode(NodeRigidBody2D(referenceCoordinates=pRef, initialVelocities=v0,
                                            visualization=VNodeRigidBody2D(drawSize=dimZ*2)))
        oMass = mbs.AddObject(ObjectRigidBody2D(mass=RBinertia.mass, inertia=RBinertia.GetInertia6D()[2],
                                                nodeNumber=nMass, visualization=
                                                VObjectRigidBody2D(graphicsData=gList)))
        mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
        mGroundWheel = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=p))
        frictionMaterialIndex=0
        if i < 2:
            mbs.AddObject(RevoluteJoint2D(markerNumbers=[mGroundWheel, mNode]))
        else:
            #add measurement roll:
            dRoll = dRelRoll * (2*sqrt(kRoll*mRoll))
            exu.Print('measurement roll resonance=',sqrt(kRoll/mRoll)/(2*pi))
            frictionMaterialIndex=1 #no friction
            mCWheel0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=0))
            mCWheel1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=1))
            mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCWheel0]))
            mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mCoordinateGround, mCWheel1],
                                                 stiffness = 1e-2, damping=dRoll,
                                                 offset=-fRoll/kRoll,
                                                 visualization=VCoordinateSpringDamper(show=False)))
            sMeasureRoll = mbs.AddSensor(SensorNode(nodeNumber=nMass, 
                                                    storeInternal=True, #read by PlotSensor, not from a file (#2492)
                                                    outputVariableType=exu.OutputVariableType.Displacement))
        
        if i == 0:
            if velocityControl:
                mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
                mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheel],
                                                   velocityLevel=True, offsetUserFunction_t=UFvelocityDrive))
                
            else:
                mbs.AddLoad(LoadTorqueVector(markerNumber=mNode, loadVector=[0,0,torque]))
        if i == 1:
            mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
            mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mCoordinateGround, mCoordinateWheel], 
                                                 damping=dampingWheel1, 
                                                 visualization=VCoordinateSpringDamper(show=False)))

        gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=contactStiffness, 
                                     contactDamping=contactDamping, frictionMaterialIndex=frictionMaterialIndex)
        
        sAngVel += [mbs.AddSensor(SensorNode(nodeNumber=nMass, storeInternal=True, #(#2492)
                                  outputVariableType=exu.OutputVariableType.AngularVelocity))]

    allCables = []
    for ancf in ancfList:
        allCables += ancf[1]
        
    for oIndex in allCables:
        gContact.AddANCFCable(objectIndex=oIndex, halfHeight=halfHeight, #halfHeight should be h/2, but then cylinders should be smaller
                              contactStiffness=contactStiffness, contactDamping=contactDamping, 
                              frictionMaterialIndex=0)

    frictionMatrix = np.zeros((2,2))
    frictionMatrix[0,0]=int(useFriction)*dryFriction
    frictionMatrix[0,1]=0 #no friction between some rolls and cable
    frictionMatrix[1,0]=0 #no friction between some rolls and cable
    gContact.SetFrictionPairings(frictionMatrix)


mbs.Assemble()
#exu.Print(mbs)

simulationSettings = exu.SimulationSettings() #takes currently set values or default values

simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
simulationSettings.solution.file.name = 'solution/coordinatesSolution.txt'
simulationSettings.solution.file.write = (not testIsActive) #only the SolutionViewer reads it (#2492)
simulationSettings.solution.file.writePeriod = 0.005
simulationSettings.solution.sensors.writePeriod = 0.001
#simulationSettings.show.computationTime = True
#1 thread, because the multithreaded assembly sums in a non-deterministic order: with 4
#threads the result of this model moved by up to 5.4e-14 between runs, past the 5e-14 test
#tolerance, and at 16 elements per section 1 thread is also FASTER (0.36s vs 0.69s). Use more
#threads for > 100 ANCF elements, but then the value is not a reference any more (#2368)
simulationSettings.parallel.numberOfThreads = 1
simulationSettings.show.statistics = True

doDynamic = True
tEnd = 0.1
h = 0.5e-3
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
simulationSettings.timeIntegration.stepInformation= 3+128+256

simulationSettings.timeIntegration.verboseMode = 1 #otherwise, load steps are shown ...

simulationSettings.timeIntegration.newton.useModifiedNewton = True

#+++++++++++++++++++++++++++++++++++++++++++++
#static:
simulationSettings.staticSolver.verboseMode = 1 #otherwise, load steps are shown ...
simulationSettings.staticSolver.numberOfLoadSteps  = 40
simulationSettings.staticSolver.discontinuous.iterationTolerance = 1
simulationSettings.staticSolver.stabilizerODE2term = 2 #may only act on position degrees of freedom
#+++++++++++++++++++++++++++++++++++++++++++++

SC.visualizationSettings.view0.scene.drawWorldBasis=True
SC.visualizationSettings.nodes.show = True
SC.visualizationSettings.nodes.defaultSize = h*20
SC.visualizationSettings.loads.show = False

SC.visualizationSettings.contour.outputVariableComponent=0
SC.visualizationSettings.contour.outputVariable=exu.OutputVariableType.ForceLocal

#visualize contact:
if False:
    SC.visualizationSettings.contact.showSearchTree =True
    SC.visualizationSettings.contact.showSearchTreeCells =True
    SC.visualizationSettings.contact.showBoundingBoxes = True

if not testIsActive: 
    SC.renderer.Start()
    SC.renderer.DoIdleTasks()

if doDynamic :
    mbs.SolveDynamic(simulationSettings) #183 Newton iterations, 0.114 seconds
else:
    mbs.SolveStatic(simulationSettings) #183 Newton iterations, 0.114 seconds


if not testIsActive:
    SC.visualizationSettings.general.autoFitScene = False
    SC.visualizationSettings.general.graphicsUpdateInterval=0.02
    
    sol = LoadSolutionFile('solution/coordinatesSolution.txt', safeMode=True)#, maxRows=100)
    exu.Print('start SolutionViewer')
    mbs.SolutionViewer(sol)


if not testIsActive: 
    SC.renderer.DoIdleTasks()
    SC.renderer.Stop() #safely close rendering window!
    
    mbs.PlotSensor(sensorNumbers=[sAngVel[0],sAngVel[1]], components=2, closeAll=True)
    mbs.PlotSensor(sensorNumbers=sMeasureRoll, components=1)
        

#print representative result:
#sol = mbs.systemData.GetODE2Coordinates()
#n = len(sol)
measurePos = mbs.GetSensorValues(sMeasureRoll)

exu.Print('solution of ANCFbeltDrive (sensor y-disp): ',measurePos[1]) 
#updated 2026-09-17; the old value was not reproducible - missing simulation parameters.
#Re-measured single-threaded (see numberOfThreads above), which IS reproducible bit for bit;
#the 4-thread value was -0.0011715990134786858 (#2368)
exu.sys['testResult'] = measurePos[1] # use y-coordinate
```