You can view and download this file on Github: beltDriveReevingSystem.py
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: Model for belt drive; According to thread deliverable D2.2 Test case
#
# Author: Johannes Gerstmayr
# Date: 2022-02-27
#
# 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 *
from exudyn.utilities import *
from exudyn.beams import *
import numpy as np
from math import sin, cos, pi, sqrt , asin, acos, atan2, exp
import copy
SC = exu.SystemContainer()
mbs = SC.AddSystem()
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
improvedBelt = True #True: improved belt model (tEnd ~= 2.5 seconds simulation, more damping, better initial conditions, etc.)
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#Parameters for the belt
gVec = [0,-9.81*1,0] # gravity
Emodulus=1e7*1 # Young's modulus of ANCF element in N/m^2
b=0.08 #0.002 # width of rectangular ANCF element in m
hc = 0.01 # height (geometric) of rectangular ANCF element in m
hcStiff = 0.01 # stiffness relevant height
rhoBeam=1036. # density of ANCF element in kg/m^3
A=b*hcStiff # cross sectional area of ANCF element in m^2
I=(b*hcStiff**3)/12 # second moment of area of ANCF element in m^4
EI = Emodulus*I
EA = Emodulus*A
rhoA = rhoBeam*A
dEI = 0*1e-3*Emodulus*I #REMARK: bending proportional damping. Set zero in the 2013 paper there is not. We need the damping for changing the initial configuration.
#dEA = 0.1*1e-2*Emodulus*A #axial strain proportional damping. Same as for the
dEA = 1 #dEA=1 in paper PechsteinGerstmayr 2013, according to HOTINT C++ files ...
# bending damping.
#%%
#%%
#settings:
useGraphics= True
useContact = True
doDynamic = True
makeAnimation = False
velocityControl = True
staticEqulibrium = False #False in 2013 paper; starts from pre-deformed reference configuration
useBristleModel = improvedBelt
#in 2013 paper, reference curvature is set according to initial geometry and released until tAccStart
preCurved = False #uses preCurvature according to straight and curved initial segments
strainIsRelativeToReference = 1. #0: straight reference, 1.: curved reference
useContactCircleFriction = True
movePulley = False #as in 2013 paper, move within first 0.05 seconds; but this does not work with Index 2 solver
tEnd = 1#1*2.25#*0.1 #*5 #end time of dynamic simulation
stepSize = 0.25*1e-4 #accurate: 2.5e-5 # for frictionVelocityPenalty = 1e7*... it must be not larger than 2.5e-5
if improvedBelt:
stepSize = 1e-4
discontinuousIterations = 3 #larger is more accurate, but smaller step size is equivalent
#h = 1e-3 #step size
tAccStart = 0.05
tAccEnd = 0.6
omegaFinal = 12
useFriction = True
dryFriction = 0.5#0.5#1.2
contactStiffness = 1e8#2e5
contactDamping = 0#1e-3*contactStiffness
nSegments = 2 #4, for nANCFnodes=60, nSegments = 2 lead to less oscillations inside, but lot of stick-slip...
nANCFnodes = 8*30#2*60#120 works well, 60 leads to oscillatory tangent/normal forces for improvedBelt=True
wheelMass = 50#1 the wheel mass is not given in the paper, only the inertia
# for the second pulley
wheelInertia = 0.25#0.01
rotationDampingWheels = 0 #zero in example in 2013 paper; torque proportional to rotation speed
#torque = 1
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#create circles
#complicated shape:
initialDisplacement = -0.0025 #not used in improvedBelt!
initialDisplacement0 = initialDisplacement*int(1-movePulley) #this is set at t=0
#h = 0.25e-3
radiusPulley = 0.09995
positionPulley2x = 0.1*pi
#preStretch = -1*(pi*0.4099+0.005)/ pi*0.4099
initialDistance = positionPulley2x - 0
initialLength = 2*initialDistance +2* pi*(radiusPulley + hcStiff/2)
finalLength = initialLength - 2* initialDisplacement0
preStretch = -(finalLength - initialLength)/ initialLength
factorStriplen = (2*initialDistance+2*pi*radiusPulley)/(2*initialDistance+2*pi*(radiusPulley + hcStiff/2));
print('factorStriplen =', factorStriplen )
preStretch += (1-1./factorStriplen) #this is due to an error in the original paper 2013
if improvedBelt:
rotationDampingWheels = 2 #to reduce vibrations of driven pulley
tEnd = 2.45 #at 2.45 node 1 is approximately at initial position!
preStretch = -0.05
initialDisplacement0 = 0
staticEqulibrium = True
strainIsRelativeToReference = False
#dryFriction = 0
hc *= 0.01
EI *= 0.02
print('preStretch=', preStretch)
circleList = [[[initialDisplacement0,0], radiusPulley,'L'],
[[positionPulley2x,0], radiusPulley,'L'],
# [[initialDisplacement0,0], radiusPulley,'L'],
# [[positionPulley2x,0], radiusPulley,'L'],
]
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#create geometry:
reevingDict = CreateReevingCurve(circleList, drawingLinesPerCircle = 64,
radialOffset=0.5*hc, closedCurve=True, #allows closed curve
numberOfANCFnodes=nANCFnodes, graphicsNodeSize= 0.01)
# set precurvature at location of pulleys:
elementCurvatures = [] #no pre-curvatures
if preCurved:
elementCurvatures = reevingDict['elementCurvatures']
gList=[]
if False: #visualize reeving curve, without simulation
gList = reevingDict['graphicsDataLines'] + reevingDict['graphicsDataCircles']
oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(show=False)))#, visualization = {'show : False'}
nGround = mbs.AddNode(NodePointGround())
mCoordinateGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
#mbs.SetObjectParameter(objectNumber = oGround, parameterName = 'Vshow', value=False)
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#create ANCF elements:
dimZ = b #z.dimension
cableTemplate = Cable2D(#physicsLength = L / nElements, #set in GenerateStraightLineANCFCable2D(...)
physicsMassPerLength = rhoA,
physicsBendingStiffness = EI,
physicsAxialStiffness = EA,
physicsBendingDamping = dEI,
physicsAxialDamping = dEA,
physicsReferenceAxialStrain = preStretch*int(improvedBelt), #prestretch
physicsReferenceCurvature = 0.,#-1/(radiusPulley + hc/2),
useReducedOrderIntegration = 2, #2=improved axial strain in postprocessing!
strainIsRelativeToReference = strainIsRelativeToReference,
visualization=VCable2D(drawHeight=hc),
)
ancf = PointsAndSlopes2ANCFCable2D(mbs, reevingDict['ancfPointsSlopes'],
reevingDict['elementLengths'],
cableTemplate, massProportionalLoad=gVec,
fixedConstraintsNode0=[1*staticEqulibrium,0,0,0], #fixedConstraintsNode1=[1,1,1,1],
elementCurvatures = elementCurvatures,
firstNodeIsLastNode=True, graphicsSizeConstraints=0.01)
if useContactCircleFriction:
lElem = reevingDict['totalLength'] / nANCFnodes
cFact=b*lElem/nSegments #stiffness shall be per area, but is applied at every segment
print('cFact=',cFact, ', lElem=', lElem)
contactStiffness*=cFact
contactDamping = 2000*cFact #according to Dufva 2008 paper ... seems also to be used in 2013 PEchstein Gerstmayr
if useBristleModel:
frictionStiffness = 1e8*cFact #1e7 converges good; 1e8 is already quite accurate
massSegment = rhoA*lElem/nSegments
frictionVelocityPenalty = 1*sqrt(frictionStiffness*massSegment) #bristle damping; should be adjusted to reduce vibrations induced by bristle model
else:
frictionVelocityPenalty = 0.1*1e7*cFact #1e7 is original in 2013 paper; requires smaller time step
#frictionVelocityPenalty = 0.25e7*cFact # 0.25e7*cFact with discontinuous.maxIterations = 4
frictionStiffness = 0 #as in 2013 paper
if improvedBelt:
frictionStiffness *= 10*5
frictionVelocityPenalty *= 10
contactStiffness *= 10*4
contactDamping *=10*4
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#create sensors for all nodes
sMidVel = []
sAxialForce = []
sCable0Pos = []
# sObjectDisp =[]
ancfNodes = ancf[0]
ancfObjects = ancf[1]
positionList2Node = [] #axial position at x=0 and x=0.5*lElem
positionListMid = [] #axial position at midpoint of element
positionListSegments = [] #axial position at midpoint of segments
currentPosition = 0 #is increased at every iteration
for i,obj in enumerate(ancfObjects):
lElem = reevingDict['elementLengths'][i]
positionList2Node += [currentPosition, currentPosition + 0.5*lElem]
positionListMid += [currentPosition + 0.5*lElem]
for j in range(nSegments):
segPos = (j+0.5)*lElem/nSegments + currentPosition
positionListSegments += [segPos]
currentPosition += lElem
sAxialForce += [mbs.AddSensor(SensorBody(bodyNumber = obj,
storeInternal=True,
localPosition=[0.*lElem,0,0],
outputVariableType=exu.OutputVariableType.ForceLocal))]
sAxialForce += [mbs.AddSensor(SensorBody(bodyNumber = obj,
storeInternal=True,
localPosition=[0.5*lElem,0,0],
outputVariableType=exu.OutputVariableType.ForceLocal))]
sMidVel += [mbs.AddSensor(SensorBody(bodyNumber = obj,
storeInternal=True,
localPosition=[0.5*lElem,0,0], #0=at left node
outputVariableType=exu.OutputVariableType.VelocityLocal))]
sCable0Pos += [mbs.AddSensor(SensorBody(bodyNumber = obj,
storeInternal=True,
localPosition=[0.*lElem,0,0],
outputVariableType=exu.OutputVariableType.Position))]
# sObjectDisp += [mbs.AddSensor(SensorBody(bodyNumber = obj,
# storeInternal=True,
# localPosition=[0.5*lElem,0,0],
# outputVariableType=exu.OutputVariableType.Displacement))]
#for testing, fix two nodes:
if False:
ii0 = 1
ii1 = 14
for i in range(4):
mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=ancfNodes[ii0], coordinate=i))
mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=ancfNodes[ii1], coordinate=i))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mANCF0],
visualization=VCoordinateConstraint(show = False)))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mANCF1],
visualization=VCoordinateConstraint(show = False)))
#reference solution for clamped-clamped beam:
lElem = reevingDict['elementLengths'][0] #all same
L = lElem * 13 #span is 13 elements long
wMax = rhoA*9.81*L**4 /(384*EI)
print('wMax=',wMax)
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#add contact:
if useContact:
contactObjects = [[],[]] #list of contact objects
gContact = mbs.AddGeneralContact()
gContact.verboseMode = 1
gContact.frictionProportionalZone = 0.005 #limit velocity. I didn't find
#gContact.frictionVelocityPenalty = 0*1e3 #limit velocity. I didn't find
#this in the paper
gContact.ancfCableUseExactMethod = False
gContact.ancfCableNumberOfContactSegments = nSegments
ssx = 16#32 #search tree size
ssy = 16#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])) #automatically computed!
dimZ= 0.01 #for drawing
sWheelRot = [] #sensors for angular velocity
nMassList = []
wheelSprings = [] #for static computation
for i, wheel in enumerate(circleList):
p = [wheel[0][0], wheel[0][1], 0] #position of wheel center
r = wheel[1]
rot0 = 0 #initial rotation
pRef = [p[0], p[1], rot0]
gList = [GraphicsDataCylinder(pAxis=[0,0,-dimZ],vAxis=[0,0,-dimZ], radius=r,
color= color4dodgerblue, nTiles=64),
GraphicsDataArrow(pAxis=[0,0,0], vAxis=[-0.9*r,0,0], radius=0.01*r, color=color4orange),
GraphicsDataArrow(pAxis=[0,0,0], vAxis=[0.9*r,0,0], radius=0.01*r, color=color4orange)]
omega0 = 0 #initial angular velocity
v0 = np.array([0,0,omega0])
nMass = mbs.AddNode(NodeRigidBody2D(referenceCoordinates=pRef, initialVelocities=v0,
visualization=VNodeRigidBody2D(drawSize=dimZ*2)))
nMassList += [nMass]
oMass = mbs.AddObject(ObjectRigidBody2D(physicsMass=wheelMass, physicsInertia=wheelInertia,
nodeNumber=nMass, visualization=
VObjectRigidBody2D(graphicsData=gList)))
mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
mGroundWheel = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=p, visualization = VMarkerBodyRigid(show = False)))
#mbs.AddObject(RevoluteJoint2D(markerNumbers=[mGroundWheel, mNode], visualization=VRevoluteJoint2D(show=False)))
mCoordinateWheelX = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=0))
mCoordinateWheelY = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=1))
constraintX = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheelX],
visualization=VCoordinateConstraint(show = False)))
constraintY = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheelY],
visualization=VCoordinateConstraint(show = False)))
if i==0:
constraintPulleyLeftX = constraintX
if True:
sWheelRot += [mbs.AddSensor(SensorNode(nodeNumber=nMass,
storeInternal=True,
fileName='solutionDelete/wheel'+str(i)+'angVel.txt',
outputVariableType=exu.OutputVariableType.AngularVelocity))]
tdisplacement = 0.05
def UFvelocityDrive(mbs, t, itemNumber, lOffset): #time derivative of UFoffset
if t < tAccStart:
v = 0
if t >= tAccStart and t < tAccEnd:
v = omegaFinal/(tAccEnd-tAccStart)*(t-tAccStart)
elif t >= tAccEnd:
v = omegaFinal
return v
if doDynamic:
if i == 0:
if velocityControl:
mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
velControl = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheel],
velocityLevel=True, offsetUserFunction_t= UFvelocityDrive,
visualization=VCoordinateConstraint(show = False)))#UFvelocityDrive
if i == 1:
mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mCoordinateGround, mCoordinateWheel],
damping = rotationDampingWheels,
visualization=VCoordinateSpringDamper(show = False)))
#this is used for times > 1 in order to see influence of torque step in Wheel1
def UFforce(mbs, t, load):
tau = 0.
tau += 25.*(SmoothStep(t, 1., 1.5, 0., 1.) - SmoothStep(t, 3.5, 4., 0., 1.))
#tau += 16.*(SmoothStep(t, 5, 5.5, 0., 1.) - SmoothStep(t, 7.5, 8., 0., 1.))
return -tau
mbs.AddLoad(LoadCoordinate(markerNumber=mCoordinateWheel,
load = 0, loadUserFunction = UFforce))
if staticEqulibrium:
mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
csd = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheel],
visualization=VCoordinateConstraint(show = False)))
wheelSprings += [csd]
frictionMaterialIndex=0
gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=contactStiffness,
contactDamping=contactDamping, frictionMaterialIndex=frictionMaterialIndex)
if not useContactCircleFriction:
for oIndex in ancf[1]:
gContact.AddANCFCable(objectIndex=oIndex, halfHeight= hc/2, #halfHeight should be h/2, but then cylinders should be smaller
contactStiffness=contactStiffness, contactDamping=contactDamping,
frictionMaterialIndex=0)
else:
cableList = ancf[1]
mCircleBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass))
#mCircleBody = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
for k in range(len(cableList)):
initialGapList = [0.1]*nSegments + [-2]*(nSegments) + [0]*(nSegments) #initial gap of 0., isStick (0=slip, +-1=stick, -2 undefined initial state), lastStickingPosition (0)
mCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=cableList[k],
numberOfSegments = nSegments, verticalOffset=-hc/2))
nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,
numberOfDataCoordinates=nSegments*(1+2) ))
co = mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mCircleBody, mCable], nodeNumber = nodeDataContactCable,
numberOfContactSegments=nSegments,
contactStiffness = contactStiffness,
contactDamping=contactDamping,
frictionVelocityPenalty = frictionVelocityPenalty,
frictionStiffness = frictionStiffness,
frictionCoefficient=int(useFriction)*dryFriction,
circleRadius = r,
visualization=VObjectContactFrictionCircleCable2D(showContactCircle=False)))
contactObjects[i] += [co]
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)
#+++++++++++++++++++++++++++++++++++++++++++
#create list of sensors for contact
sContactDisp = [[],[]]
sContactForce = [[],[]]
for i in range(len(contactObjects)):
for obj in contactObjects[i]:
sContactForce[i] += [mbs.AddSensor(SensorObject(objectNumber = obj,
storeInternal=True,
outputVariableType=exu.OutputVariableType.ForceLocal))]
sContactDisp[i] += [mbs.AddSensor(SensorObject(objectNumber = obj,
storeInternal=True,
outputVariableType=exu.OutputVariableType.Coordinates))]
#user function to smoothly transform from curved to straight reference configuration as
#in paper 2013, Pechstein, Gerstmayr
def PreStepUserFunction(mbs, t):
if True and t <= tAccStart+1e-10:
cableList = ancf[1]
fact = (tAccStart-t)/tAccStart #from 1 to 0
if fact < 1e-12: fact = 0. #for very small values ...
#curvatures = reevingDict['elementCurvatures']
#print('fact=', fact)
for i in range(len(cableList)):
oANCF = cableList[i]
mbs.SetObjectParameter(oANCF, 'strainIsRelativeToReference',
fact)
mbs.SetObjectParameter(oANCF, 'physicsReferenceAxialStrain',
preStretch*(1.-fact))
# if movePulley:
# #WARNING: this does not work for Index2 solver:
# mbs.SetObjectParameter(constraintPulleyLeftX, 'offset', initialDisplacement*(1.-fact))
# #print('offset=', initialDisplacement*(1.-fact))
return True
mbs.Assemble()
simulationSettings = exu.SimulationSettings() #takes currently set values or default values
simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution_nosync/testCoords.txt'
simulationSettings.solutionSettings.writeSolutionToFile = True
simulationSettings.solutionSettings.solutionWritePeriod = 0.002
simulationSettings.solutionSettings.sensorsWritePeriod = 0.001
simulationSettings.displayComputationTime = True
simulationSettings.parallel.numberOfThreads = 1 #use 4 to speed up for > 100 ANCF elements
simulationSettings.displayStatistics = True
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
simulationSettings.timeIntegration.stepInformation= 255
simulationSettings.timeIntegration.verboseMode = 1
simulationSettings.timeIntegration.newton.useModifiedNewton = True
#simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
#simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-3
simulationSettings.timeIntegration.discontinuous.maxIterations = discontinuousIterations #3
simulationSettings.displayStatistics = True
SC.visualizationSettings.general.circleTiling = 24
SC.visualizationSettings.loads.show=False
SC.visualizationSettings.sensors.show=False
SC.visualizationSettings.markers.show=False
SC.visualizationSettings.nodes.defaultSize = 0.002
SC.visualizationSettings.openGL.multiSampling = 4
SC.visualizationSettings.openGL.lineWidth = 2
SC.visualizationSettings.window.renderWindowSize = [1920,1080]
SC.visualizationSettings.connectors.showContact = True
SC.visualizationSettings.contact.contactPointsDefaultSize = 0.0002
SC.visualizationSettings.contact.showContactForces = True
SC.visualizationSettings.contact.contactForcesFactor = 0.005
if makeAnimation == True:
simulationSettings.solutionSettings.recordImagesInterval = 0.02
SC.visualizationSettings.exportImages.saveImageFileName = "animationNew/frame"
if True:
SC.visualizationSettings.bodies.beams.axialTiling = 1
SC.visualizationSettings.bodies.beams.drawVertical = True
SC.visualizationSettings.bodies.beams.drawVerticalLines = True
SC.visualizationSettings.contour.outputVariableComponent=0
SC.visualizationSettings.contour.outputVariable=exu.OutputVariableType.ForceLocal
SC.visualizationSettings.bodies.beams.drawVerticalFactor = 0.001
SC.visualizationSettings.bodies.beams.drawVerticalOffset = -120
if improvedBelt:
SC.visualizationSettings.bodies.beams.drawVerticalFactor = 0.0003
SC.visualizationSettings.bodies.beams.drawVerticalOffset = -220
SC.visualizationSettings.bodies.beams.reducedAxialInterploation = True
# SC.visualizationSettings.contour.outputVariable=exu.OutputVariableType.VelocityLocal
# SC.visualizationSettings.bodies.beams.drawVerticalFactor = -0.25
# SC.visualizationSettings.bodies.beams.drawVerticalOffset = 0.30
# SC.visualizationSettings.bodies.beams.reducedAxialInterploation = False
#visualize contact:
if False:
SC.visualizationSettings.contact.showSearchTree =True
SC.visualizationSettings.contact.showSearchTreeCells =True
SC.visualizationSettings.contact.showBoundingBoxes = True
if useGraphics:
exu.StartRenderer()
mbs.WaitForUserToContinue()
#simulationSettings.staticSolver.newton.absoluteTolerance = 1e-10
simulationSettings.staticSolver.adaptiveStep = False
simulationSettings.staticSolver.loadStepGeometric = True;
simulationSettings.staticSolver.loadStepGeometricRange=1e4
simulationSettings.staticSolver.numberOfLoadSteps = 10
#simulationSettings.staticSolver.useLoadFactor = False
simulationSettings.staticSolver.stabilizerODE2term = 1e5
simulationSettings.staticSolver.newton.relativeTolerance = 1e-6
simulationSettings.staticSolver.newton.absoluteTolerance = 1e-6
if staticEqulibrium: #precompute static equilibrium
mbs.SetObjectParameter(velControl, 'activeConnector', False)
for i in range(len(contactObjects)):
for obj in contactObjects[i]:
mbs.SetObjectParameter(obj, 'frictionCoefficient', 0.)
mbs.SetObjectParameter(obj, 'frictionStiffness', 1e-8) #do not set to zero, as it needs to do some initialization...
# simulationSettings.solutionSettings.appendToFile=False
mbs.SolveStatic(simulationSettings, updateInitialValues=True)
# simulationSettings.solutionSettings.appendToFile=True
#check total force on support, expect: supportLeftX \approx 2*preStretch*EA
supportLeftX = mbs.GetObjectOutput(constraintPulleyLeftX,variableType=exu.OutputVariableType.Force)
print('Force x in support of left pulley = ', supportLeftX)
print('Belt pre-tension=', preStretch*EA)
for i in range(len(contactObjects)):
for obj in contactObjects[i]:
mbs.SetObjectParameter(obj, 'frictionCoefficient', dryFriction)
mbs.SetObjectParameter(obj, 'frictionStiffness', frictionStiffness)
for coordinateConstraint in ancf[4]:
mbs.SetObjectParameter(coordinateConstraint, 'activeConnector', False)
mbs.SetObjectParameter(velControl, 'activeConnector', True)
for csd in wheelSprings:
mbs.SetObjectParameter(csd, 'activeConnector', False)
else:
mbs.SetPreStepUserFunction(PreStepUserFunction)
if True:
mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.TrapezoidalIndex2) #183 Newton iterations, 0.114 seconds
#mbs.SolveDynamic(simulationSettings)
if useGraphics and False:
SC.visualizationSettings.general.autoFitScene = False
SC.visualizationSettings.general.graphicsUpdateInterval=0.02
sol = LoadSolutionFile('solution_nosync/testCoords.txt', safeMode=True)#, maxRows=100)
mbs.SolutionViewer(sol)
if useGraphics:
SC.WaitForRenderEngineStopFlag()
exu.StopRenderer() #safely close rendering window!
#%%++++++++++++++++++++++++++++++++++++++++
if False:
#shift data depending on axial position by subtracting xOff; put negative x values+shiftValue to end of array
def ShiftXoff(data, xOff, shiftValue):
indOff = 0
n = data.shape[0]
data[:,0] -= xOff
for i in range(n):
if data[i,0] < 0:
indOff+=1
data[i,0] += shiftValue
print('indOff=', indOff)
data = np.vstack((data[indOff:,:], data[0:indOff,:]))
return data
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from exudyn.plot import DataArrayFromSensorList
mbs.PlotSensor(closeAll=True)
#compute axial offset, to normalize results:
nodePos0 = mbs.GetSensorValues(sCable0Pos[0])
xOff = nodePos0[0]
maxXoff = 0.5*positionPulley2x
maxYoff = 0.1*r
# indOff = 0 #single data per element
# indOff2 = 0 #double data per element
correctXoffset = True
if abs(nodePos0[1]-r) > maxYoff or nodePos0[0] > maxXoff or nodePos0[0] < -0.1*maxXoff:
print('*****************')
print('warning: final position not at top of belt or too far away')
print('nodePos0=',nodePos0)
print('*****************')
xOff = 0
correctXoffset = False
else:
#compute offset index:
# for (i,s) in enumerate(sCable0Pos):
# p = mbs.GetSensorValues(s)
# print('p'+str(i)+'=', p)
# if p[0] > 0 and i > int(0.8*nANCFnodes):
# indOff+=1
# indOff2+=2
# indOff -= 1
# indOff2 -= 2
print('******************\nxOff=', xOff)
dataVel = DataArrayFromSensorList(mbs, sensorNumbers=sMidVel, positionList=positionListMid)
if correctXoffset:
dataVel=ShiftXoff(dataVel,xOff, reevingDict['totalLength'])
mbs.PlotSensor(sensorNumbers=[dataVel], components=0, labels=['axial velocity'],
xLabel='axial position (m)', yLabel='velocity (m/s)')
#axial force over beam length:
dataForce = DataArrayFromSensorList(mbs, sensorNumbers=sAxialForce, positionList=positionList2Node)
if correctXoffset:
dataForce = ShiftXoff(dataForce,xOff, reevingDict['totalLength'])
mbs.PlotSensor(sensorNumbers=[dataForce], components=0, labels=['axial force'], colorCodeOffset=2,
xLabel='axial position (m)', yLabel='axial force (N)')
if improvedBelt and dryFriction==0.5 and nANCFnodes==120:
#analytical exponential curve, Euler's/Eytelwein's solution:
na = 12 #number of data points
dataExp = np.zeros((na*2, 2))
#f0 = 278.733 #this is at low level, but exp starts later
f0 = 191.0#287.0
x0 = 0.5860 #1.1513 #starting coordinate, drawn in -x direction
d = 0.28 #amount along x drawn
for i in range(na):
x = i/na*d
beta = x/(radiusPulley + hc/2)
val = f0*exp(beta*dryFriction)
#print('x=',x,',exp=',val)
dataExp[i,0] = x0-x
dataExp[i,1] = val
f0 = 193.4#287.0
x0 = 0.984 #1.1513 #starting coordinate, drawn in -x direction
for i in range(na):
x = i/na*d
beta = x/(radiusPulley + hc/2)
val = f0*exp(beta*dryFriction)
dataExp[i+na,0] = x0+x
dataExp[i+na,1] = val
mbs.PlotSensor(sensorNumbers=[dataExp], components=0, labels=['analytical Eytelwein'], colorCodeOffset=3, newFigure=False,
lineStyles=[''], markerStyles=['x '], markerDensity=2*na)
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#contact forces are stored (x/y) for every segment ==> put into consecutive array
contactForces =[[],[]] #these are the contact forces of the whole belt, but from both pulleys!
for i in range(len(sContactForce)):
contactForces[i] = np.zeros((len(sContactForce[i])*nSegments, 3)) #per row: [position, Fx, Fy]
for j, sensor in enumerate(sContactForce[i]):
values = mbs.GetSensorValues(sensor)
for k in range(nSegments):
row = j*nSegments + k
contactForces[i][row,0] = positionListSegments[row]
contactForces[i][row, 1:] = values[k*2:k*2+2]
contactForcesTotal = contactForces[0]
contactForcesTotal[:,1:] += contactForces[1][:,1:]
#plot contact forces over beam length:
mbs.PlotSensor(sensorNumbers=[contactForcesTotal,contactForcesTotal], components=[0,1], labels=['tangential force','normal force'],
xLabel='axial position (m)', yLabel='contact forces (N)', newFigure=True)
# mbs.PlotSensor(sensorNumbers=[contactForces[1],contactForces[1]], components=[0,1], labels=['tangential force','normal force'],
# xLabel='axial position (m)', yLabel='contact forces (N)', newFigure=False)
contactDisp =[[],[]] #slip and gap
for i in range(len(sContactDisp)):
contactDisp[i] = np.zeros((len(sContactDisp[i])*nSegments, 3)) #per row: [position, Fx, Fy]
for j, sensor in enumerate(sContactDisp[i]):
values = mbs.GetSensorValues(sensor)
for k in range(nSegments):
row = j*nSegments + k
contactDisp[i][row,0] = positionListSegments[row]
contactDisp[i][row, 1:] = values[k*2:k*2+2]
mbs.PlotSensor(sensorNumbers=[contactDisp[0],contactDisp[0]], components=[0,1], labels=['slip','gap'],
xLabel='axial position (m)', yLabel='slip, gap (m)', newFigure=True)
mbs.PlotSensor(sensorNumbers=[contactDisp[1],contactDisp[1]], components=[0,1], labels=['slip','gap'],
xLabel='axial position (m)', yLabel='slip, gap (m)', newFigure=False)
header = ''
header += 'endTime='+str(tEnd)+'\n'
header += 'stepSize='+str(stepSize)+'\n'
header += 'nSegments='+str(nSegments)+'\n'
header += 'nANCFnodes='+str(nANCFnodes)+'\n'
header += 'contactStiffness='+str(contactStiffness)+'\n'
header += 'contactDamping='+str(contactDamping)+'\n'
header += 'frictionStiffness='+str(frictionStiffness)+'\n'
header += 'frictionVelocityPenalty='+str(frictionVelocityPenalty)+'\n'
header += 'dryFriction='+str(dryFriction)+'\n'
fstr = 'h'+str(stepSize)+'n'+str(int(nANCFnodes/60))+'s'+str(nSegments)+'cs'+str(int((contactStiffness/41800)))
fstr += 'fs'+str(int((frictionStiffness/52300)))
#export solution:
if improvedBelt:
np.savetxt('solutionDelete/contactForces'+fstr+'.txt', contactForces[0]+contactForces[1], delimiter=',',
header='Exudyn: solution of belt drive, contact forces over belt length\n'+header, encoding=None)
np.savetxt('solutionDelete/contactDisp'+fstr+'.txt', contactDisp[0]+contactDisp[1], delimiter=',',
header='Exudyn: solution of belt drive, slip and gap over belt length\n'+header, encoding=None)
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
mbs.PlotSensor(sensorNumbers=[sWheelRot[0], sWheelRot[1]], components=[2,2])#,sWheelRot[1]
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++