beltDriveALE.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Model for belt drive; using ALE ANCF cable and regular cable
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2022-07-08
  8#
  9# 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.
 10#
 11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 12
 13import exudyn as exu
 14from exudyn.itemInterface import *
 15from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 16import exudyn.graphics as graphics
 17from exudyn.beams import *
 18
 19import numpy as np
 20from math import sin, cos, pi, sqrt , asin, acos, atan2, exp
 21import copy 
 22
 23
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26
 27#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 28
 29#improvedBelt = True #True: improved belt model (tEnd ~= 2.5 seconds simulation, more damping, better initial conditions, etc.)
 30
 31#%%
 32#settings:
 33useGraphics = True
 34useContact = True
 35doDynamic = True
 36makeAnimation = False
 37useALE = True
 38useFrictionStiffness = True
 39
 40stepSize = 0.5*0.5*1e-4  #accurate: 2.5e-5 # for frictionVelocityPenalty = 1e7*... it must be not larger than 2.5e-5
 41discontinuousIterations = 0+3 #larger is more accurate, but smaller step size is equivalent
 42
 43if useFrictionStiffness:
 44    stepSize = 0.25*0.5*1e-4  #accurate: 2.5e-5 # for frictionVelocityPenalty = 1e7*... it must be not larger than 2.5e-5
 45    # discontinuousIterations = 6+3 #larger is more accurate, but smaller step size is equivalent
 46
 47#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 48#Parameters for the belt
 49gVec = [0,-9.81*1,0]     # gravity
 50Emodulus=1e7             # Young's modulus of ANCF element in N/m^2
 51b=0.08 #0.002            # width of rectangular ANCF element in m
 52hc = 0.0001              # height (geometric) of rectangular ANCF element in m
 53hcStiff = 0.01           # stiffness relevant height
 54rhoBeam=1036.            # density of ANCF element in kg/m^3
 55A=b*hcStiff              # cross sectional area of ANCF element in m^2
 56I=(b*hcStiff**3)/12      # second moment of area of ANCF element in m^4
 57EI = 0.02*Emodulus*I
 58EA = Emodulus*A
 59rhoA = rhoBeam*A
 60dEI = 0
 61dEA = 1 
 62
 63# EI *= 1000*2
 64# EI  *= 500*5 #for test
 65#%%
 66
 67#h = 1e-3 #step size
 68tAccStart = 0.05
 69tAccEnd = 0.6
 70omegaFinal = 12
 71
 72useFriction = True
 73dryFriction = 0.5#0.5#1.2
 74contactStiffness = 1e8#2e5
 75contactDamping = 0#1e-3*contactStiffness
 76
 77nSegments = 2 #4, for nANCFnodes=60, nSegments = 2 lead to less oscillations inside, but lot of stick-slip...
 78nANCFnodes =2*2*30#2*60#120 works well, 60 leads to oscillatory tangent/normal forces for improvedBelt=True
 79
 80wheelMass = 50#1 the wheel mass is not given in the paper, only the inertia 
 81# for the second pulley
 82wheelInertia = 0.25#0.01
 83rotationDampingWheels = 0 #zero in example in 2013 paper; torque proportional to rotation speed
 84
 85#torque = 1
 86
 87#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 88#create circles
 89#complicated shape:
 90#initialDisplacement = -0.0025 #not used in improvedBelt!
 91radiusPulley = 0.09995
 92positionPulley2x = 0.1*pi
 93#initialDistance = positionPulley2x
 94initialLength = 2*positionPulley2x + 2* pi*(radiusPulley + hcStiff/2)
 95#finalLength = initialLength - 2* initialDisplacement
 96#preStretch = -(finalLength - initialLength)/ initialLength 
 97
 98#factorStriplen = (2*initialDistance+2*pi*radiusPulley)/(2*initialDistance+2*pi*(radiusPulley + hcStiff/2));
 99#print('factorStriplen =', factorStriplen )
100#preStretch += (1-1./factorStriplen) #this is due to an error in the original paper 2013
101
102
103rotationDampingWheels = 2 #to reduce vibrations of driven pulley
104tEnd = 2.45 #at 2.45 node 1 is approximately at initial position!
105preStretch = -0.05
106staticEqulibrium = True
107#dryFriction = 0
108
109print('preStretch=', preStretch)
110circleList = [[[0,0], radiusPulley,'L'],
111              [[positionPulley2x,0], radiusPulley,'L'],
112              # [[initialDisplacement0,0], radiusPulley,'L'],              
113              # [[positionPulley2x,0], radiusPulley,'L'],
114              ]
115
116#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
117#create geometry:
118reevingDict = CreateReevingCurve(circleList, drawingLinesPerCircle = 64, 
119                                radialOffset=0.5*hc, closedCurve=True, #allows closed curve
120                                numberOfANCFnodes=nANCFnodes, graphicsNodeSize= 0.01)
121
122
123
124# set precurvature at location of pulleys:
125elementCurvatures = reevingDict['elementCurvatures']
126
127gList=[]
128if False: #visualize reeving curve, without simulation
129    gList = reevingDict['graphicsDataLines'] + reevingDict['graphicsDataCircles']
130
131oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(show=False)))#, visualization = {'show : False'}
132nGround = mbs.AddNode(NodePointGround())
133mCoordinateGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
134
135
136#mbs.SetObjectParameter(objectNumber = oGround, parameterName = 'Vshow', value=False)
137#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
138#create ANCF elements:
139dimZ = b #z.dimension
140
141ANCFElementType = Cable2D
142nodesANCF = [-1,-1]
143if useALE:
144    ANCFElementType = ALECable2D
145    nALE = mbs.AddNode(NodeGenericODE2(numberOfODE2Coordinates=1, referenceCoordinates=[0], 
146                                          initialCoordinates=[0], initialCoordinates_t=[0]))
147    mALE = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nALE, coordinate=0, 
148                                              visualization = {'show':True})) #ALE velocity
149    nodesANCF = [-1,-1, nALE]
150
151    #Constraint for eulerian coordinate
152    oCCvALE=mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mALE], offset=0, 
153                                                velocityLevel = False, 
154                                                visualization=VCoordinateConstraint(show=False)))
155
156cableTemplate = ANCFElementType(#length = L / nElements, #set in GenerateBeamElementsAlongLine(...)
157                        nodeNumbers = nodesANCF,
158                        massPerLength = rhoA,
159                        bendingStiffness = EI,
160                        axialStiffness = EA,
161                        bendingDamping = dEI,
162                        referenceAxialStrain = preStretch,
163                        referenceCurvature = 0.,
164                        useReducedOrderIntegration = 2,
165                        strainIsRelativeToReference = False, #would cause reference configuration to be precurved
166                        visualization=VCable2D(drawHeight=hc),
167                        )
168
169if useALE:
170    cableTemplate.addALEvariation = False
171
172ancf = PointsAndSlopes2ANCFCable2D(mbs, reevingDict['ancfPointsSlopes'], 
173                                   reevingDict['elementLengths'], 
174                                   cableTemplate, massProportionalLoad=gVec, 
175                                   fixedConstraintsNode0=[1*staticEqulibrium,0,0,0], #fixedConstraintsNode1=[1,1,1,1],
176                                   #elementCurvatures  = elementCurvatures, #do not set this, will cause inhomogeneous curvatures
177                                   firstNodeIsLastNode=True, graphicsSizeConstraints=0.01)
178
179lElem = reevingDict['totalLength'] / nANCFnodes
180cFact=b*lElem/nSegments #stiffness shall be per area, but is applied at every segment
181
182contactStiffness*=40*cFact
183contactDamping = 40*2000*cFact #according to Dufva 2008 paper ... seems also to be used in 2013 PEchstein Gerstmayr
184frictionStiffness = 50e8*cFact #1e7 converges good; 1e8 is already quite accurate
185massSegment = rhoA*lElem/nSegments
186frictionVelocityPenalty = 10*sqrt(frictionStiffness*massSegment) #bristle damping; should be adjusted to reduce vibrations induced by bristle model
187
188if useFrictionStiffness:
189    frictionStiffness*=0.1
190else:
191    frictionStiffness*=0
192    frictionVelocityPenalty*= 2
193    
194#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
195#create sensors for all nodes
196sMidVel = []
197sAxialForce = []
198sCable0Pos = []
199# sObjectDisp =[] 
200
201ancfNodes = ancf[0]
202ancfObjects = ancf[1]
203positionList2Node = [] #axial position at x=0 and x=0.5*lElem
204positionListMid = [] #axial position at midpoint of element
205positionListSegments = [] #axial position at midpoint of segments
206currentPosition = 0 #is increased at every iteration
207for i,obj in enumerate(ancfObjects):
208    lElem = reevingDict['elementLengths'][i]
209    positionList2Node += [currentPosition, currentPosition + 0.5*lElem]
210    positionListMid += [currentPosition + 0.5*lElem]
211
212    for j in range(nSegments):
213        segPos = (j+0.5)*lElem/nSegments + currentPosition
214        positionListSegments += [segPos]
215    currentPosition += lElem
216
217    sAxialForce += [mbs.AddSensor(SensorBody(bodyNumber = obj, 
218                                              storeInternal=True,
219                                              localPosition=[0.*lElem,0,0], 
220                                              outputVariableType=exu.OutputVariableType.ForceLocal))]
221    sAxialForce += [mbs.AddSensor(SensorBody(bodyNumber = obj, 
222                                              storeInternal=True,
223                                              localPosition=[0.5*lElem,0,0], 
224                                              outputVariableType=exu.OutputVariableType.ForceLocal))]
225    sMidVel += [mbs.AddSensor(SensorBody(bodyNumber = obj, 
226                                          storeInternal=True,
227                                          localPosition=[0.5*lElem,0,0], #0=at left node
228                                          outputVariableType=exu.OutputVariableType.VelocityLocal))]
229    sCable0Pos += [mbs.AddSensor(SensorBody(bodyNumber = obj, 
230                                            storeInternal=True,
231                                            localPosition=[0.*lElem,0,0],
232                                            outputVariableType=exu.OutputVariableType.Position))]
233    # sObjectDisp += [mbs.AddSensor(SensorBody(bodyNumber = obj, 
234    #                                           storeInternal=True,
235    #                                           localPosition=[0.5*lElem,0,0],
236    #                                           outputVariableType=exu.OutputVariableType.Displacement))]
237
238
239    
240#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
241#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
242#add contact:
243if useContact:
244
245    contactObjects = [[],[]] #list of contact objects
246    
247    dimZ= 0.01 #for drawing
248    sWheelRot = [] #sensors for angular velocity
249
250    nMassList = []
251    wheelSprings = [] #for static computation
252    for i, wheel in enumerate(circleList):
253        p = [wheel[0][0], wheel[0][1], 0] #position of wheel center
254        r = wheel[1]
255    
256        rot0 = 0 #initial rotation
257        pRef = [p[0], p[1], rot0]
258        gList = [graphics.Cylinder(pAxis=[0,0,-dimZ],vAxis=[0,0,-dimZ], radius=r,
259                                      color= graphics.color.dodgerblue, nTiles=32),
260                 graphics.Arrow(pAxis=[0,0,0], vAxis=[-0.9*r,0,0], radius=0.01*r, color=graphics.color.orange),
261                 graphics.Arrow(pAxis=[0,0,0], vAxis=[0.9*r,0,0], radius=0.01*r, color=graphics.color.orange)]
262
263        omega0 = 0 #initial angular velocity
264        v0 = np.array([0,0,omega0]) 
265
266        nMass = mbs.AddNode(NodeRigidBody2D(referenceCoordinates=pRef, initialVelocities=v0,
267                                            visualization=VNodeRigidBody2D(drawSize=dimZ*2)))
268        nMassList += [nMass]
269        oMass = mbs.AddObject(ObjectRigidBody2D(mass=wheelMass, inertia=wheelInertia,
270                                                nodeNumber=nMass, visualization=
271                                                VObjectRigidBody2D(graphicsData=gList)))
272        mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
273        mGroundWheel = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=p, visualization = VMarkerBodyRigid(show = False)))
274    
275        #mbs.AddObject(RevoluteJoint2D(markerNumbers=[mGroundWheel, mNode], visualization=VRevoluteJoint2D(show=False)))
276
277        mCoordinateWheelX = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=0))
278        mCoordinateWheelY = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=1))
279        constraintX = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheelX],
280                                                 visualization=VCoordinateConstraint(show = False)))
281        constraintY = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheelY],
282                                                 visualization=VCoordinateConstraint(show = False)))
283        if i==0:
284            constraintPulleyLeftX = constraintX
285
286        if True:
287        
288            sWheelRot += [mbs.AddSensor(SensorNode(nodeNumber=nMass, 
289                                                   storeInternal=True,
290                                                   fileName='solution/wheel'+str(i)+'angVel.txt',
291                                                   outputVariableType=exu.OutputVariableType.AngularVelocity))]
292        tdisplacement = 0.05
293  
294                         
295        def UFvelocityDrive(mbs, t, itemNumber, lOffset): #time derivative of UFoffset
296            if t < tAccStart:
297                v = 0
298            if t >= tAccStart and t < tAccEnd:
299                v = omegaFinal/(tAccEnd-tAccStart)*(t-tAccStart)
300            elif t >= tAccEnd:
301                v = omegaFinal
302            return v    
303        
304        if doDynamic:    
305            if i == 0:
306                mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
307                velControl = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheel],
308                                                    velocityLevel=True, offsetUserFunction_t= UFvelocityDrive,
309                                                    visualization=VCoordinateConstraint(show = False)))#UFvelocityDrive
310            if i == 1:
311                mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
312                mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mCoordinateGround, mCoordinateWheel],
313                                                     damping = rotationDampingWheels,
314                                                     visualization=VCoordinateSpringDamper(show = False)))
315                
316                #this is used for times > 1 in order to see influence of torque step in Wheel1
317                def UFforce(mbs, t, load):
318                    tau = 0.
319                    tau +=  25.*(SmoothStep(t, 1., 1.5, 0., 1.) - SmoothStep(t, 3.5, 4., 0., 1.))
320                    #tau += 16.*(SmoothStep(t, 5, 5.5, 0., 1.) - SmoothStep(t, 7.5, 8., 0., 1.))
321                    return -tau
322                
323                mbs.AddLoad(LoadCoordinate(markerNumber=mCoordinateWheel,
324                                           load = 0, loadUserFunction = UFforce))
325
326        if staticEqulibrium:
327            mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
328            csd = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheel],
329                                                     visualization=VCoordinateConstraint(show = False)))
330            wheelSprings += [csd]
331        
332
333        cableList = ancf[1]
334        mCircleBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass))
335        #mCircleBody = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
336        for k in range(len(cableList)):
337            initialGapList = [0.1]*nSegments + [-2]*(nSegments) + [0]*(nSegments) #initial gap of 0., isStick (0=slip, +-1=stick, -2 undefined initial state), lastStickingPosition (0)
338
339            mCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=cableList[k], 
340                                                          numberOfSegments = nSegments, verticalOffset=-0*hc/2))
341            nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,
342                                                               numberOfDataCoordinates=nSegments*(1+2) ))
343
344            co = mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mCircleBody, mCable], nodeNumber = nodeDataContactCable, 
345                                                     numberOfContactSegments=nSegments, 
346                                                     contactStiffness = contactStiffness, 
347                                                     contactDamping=contactDamping, 
348                                                     frictionVelocityPenalty = frictionVelocityPenalty, 
349                                                     frictionStiffness = frictionStiffness, 
350                                                     frictionCoefficient=int(useFriction)*dryFriction,
351                                                     circleRadius = r,
352                                                     visualization=VObjectContactFrictionCircleCable2D(showContactCircle=False)))
353            contactObjects[i] += [co]
354
355
356
357#+++++++++++++++++++++++++++++++++++++++++++
358#create list of sensors for contact
359sContactDisp = [[],[]]
360sContactForce = [[],[]]
361for i in range(len(contactObjects)):
362    for obj in contactObjects[i]:
363        sContactForce[i] += [mbs.AddSensor(SensorObject(objectNumber = obj, 
364                                                        storeInternal=True,
365                                                        outputVariableType=exu.OutputVariableType.ForceLocal))]
366        sContactDisp[i] += [mbs.AddSensor(SensorObject(objectNumber = obj, 
367                                                        storeInternal=True,
368                                                        outputVariableType=exu.OutputVariableType.Coordinates))]
369        
370
371mbs.Assemble()
372
373
374simulationSettings = exu.SimulationSettings() #takes currently set values or default values
375
376simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
377simulationSettings.solution.file.name = 'solution/testCoords.txt'
378
379simulationSettings.solution.file.write = True
380simulationSettings.solution.file.writePeriod = 0.002
381simulationSettings.solution.sensors.writePeriod = 0.001
382simulationSettings.parallel.numberOfThreads = 1 #use 4 to speed up for > 100 ANCF elements
383
384simulationSettings.timeIntegration.endTime = tEnd
385simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
386simulationSettings.timeIntegration.stepInformation= 255
387
388simulationSettings.timeIntegration.verboseMode = 1
389
390simulationSettings.timeIntegration.newton.useModifiedNewton = True
391#simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
392#simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
393
394simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-3
395simulationSettings.timeIntegration.discontinuous.maxIterations = discontinuousIterations #3
396
397simulationSettings.show.statistics = True
398simulationSettings.show.computationTime = False
399
400
401SC.visualizationSettings.general.circleTiling = 24
402SC.visualizationSettings.loads.show=False
403SC.visualizationSettings.sensors.show=False
404SC.visualizationSettings.markers.show=False
405SC.visualizationSettings.nodes.defaultSize = 0.002
406SC.visualizationSettings.openGL.multiSampling = 4
407SC.visualizationSettings.openGL.lineWidth = 2
408SC.visualizationSettings.view0.window.renderWindowSize = [1920,1080]
409
410SC.visualizationSettings.connectors.showContact = True
411SC.visualizationSettings.contact.contactPointsDefaultSize = 0.0002
412SC.visualizationSettings.contact.showContactForces = True
413SC.visualizationSettings.contact.contactForcesFactor = 0.005
414
415
416if True:
417    SC.visualizationSettings.bodies.beams.axialTiling = 1
418    SC.visualizationSettings.bodies.beams.drawVertical = True
419    SC.visualizationSettings.bodies.beams.drawVerticalLines = True
420
421    SC.visualizationSettings.contour.outputVariableComponent=0
422    SC.visualizationSettings.contour.outputVariable=exu.OutputVariableType.ForceLocal
423
424    SC.visualizationSettings.bodies.beams.drawVerticalFactor = 0.0003
425    SC.visualizationSettings.bodies.beams.drawVerticalOffset = -220
426        
427    SC.visualizationSettings.bodies.beams.reducedAxialInterploation = True
428    
429    # SC.visualizationSettings.contour.outputVariable=exu.OutputVariableType.VelocityLocal
430    # SC.visualizationSettings.bodies.beams.drawVerticalFactor = -0.25
431    # SC.visualizationSettings.bodies.beams.drawVerticalOffset = 0.30
432    # SC.visualizationSettings.bodies.beams.reducedAxialInterploation = False
433
434
435if useGraphics: 
436    SC.renderer.Start()
437    #SC.renderer.DoIdleTasks()
438
439#simulationSettings.staticSolver.newton.absoluteTolerance = 1e-10
440simulationSettings.staticSolver.adaptiveStep = True
441simulationSettings.staticSolver.loadStepGeometric = False;
442simulationSettings.staticSolver.loadStepGeometricRange=1e4
443simulationSettings.staticSolver.numberOfLoadSteps = 10
444#simulationSettings.staticSolver.useLoadFactor = False
445simulationSettings.staticSolver.stabilizerODE2term = 1e5*10
446simulationSettings.staticSolver.newton.relativeTolerance = 1e-6
447simulationSettings.staticSolver.newton.absoluteTolerance = 1e-6
448
449if staticEqulibrium: #precompute static equilibrium
450    mbs.SetObjectParameter(velControl, 'activeConnector', False)
451
452    for i in range(len(contactObjects)):
453        for obj in contactObjects[i]:
454            mbs.SetObjectParameter(obj, 'frictionCoefficient', 0.)
455            mbs.SetObjectParameter(obj, 'frictionStiffness', 1e-8) #do not set to zero, as it needs to do some initialization...
456            
457    # simulationSettings.solution.file.append=False
458    mbs.SolveStatic(simulationSettings, updateInitialValues=True)
459    # simulationSettings.solution.file.append=True    
460
461    #check total force on support, expect: supportLeftX \approx 2*preStretch*EA
462    supportLeftX = mbs.GetObjectOutput(constraintPulleyLeftX,variableType=exu.OutputVariableType.Force)
463    print('Force x in support of left pulley = ', supportLeftX)
464    print('Belt pre-tension=', preStretch*EA)
465    
466    for i in range(len(contactObjects)):
467        for obj in contactObjects[i]:
468            mbs.SetObjectParameter(obj, 'frictionCoefficient', dryFriction)
469            mbs.SetObjectParameter(obj, 'frictionStiffness', frictionStiffness)
470
471    # useALE = False
472    for coordinateConstraint in ancf[4]: #release constraints for dynamic Solver
473        if not useALE: #except ALE constraint
474            mbs.SetObjectParameter(coordinateConstraint, 'activeConnector', False)
475    
476    if useALE: 
477        mbs.SetObjectParameter(oCCvALE, 'activeConnector', False) #do not fix ALE coordinate any more
478        
479        
480    mbs.SetObjectParameter(velControl, 'activeConnector', True)
481    for csd in wheelSprings:
482        mbs.SetObjectParameter(csd, 'activeConnector', False)
483
484if True:
485    mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.TrapezoidalIndex2) #183 Newton iterations, 0.114 seconds
486#mbs.SolveDynamic(simulationSettings)
487
488if useGraphics and True:
489    SC.visualizationSettings.general.autoFitScene = False
490    SC.visualizationSettings.general.graphicsUpdateInterval=0.02
491    
492    sol = LoadSolutionFile(OutputFilePath('solution/testCoords.txt'), safeMode=True)#, maxRows=100)
493    mbs.SolutionViewer(sol)
494
495
496if useGraphics: 
497    SC.renderer.DoIdleTasks()
498    SC.renderer.Stop() #safely close rendering window!
499
500#%%++++++++++++++++++++++++++++++++++++++++
501if True:
502    solDir = OutputFilePath('solution/')
503    #shift data depending on axial position by subtracting xOff; put negative x values+shiftValue to end of array
504    def ShiftXoff(data, xOff, shiftValue):
505        indOff = 0
506        n = data.shape[0]
507        data[:,0] -= xOff
508        for i in range(n):
509           if data[i,0] < 0:
510               indOff+=1
511               data[i,0] += shiftValue
512        print('indOff=', indOff)
513        data = np.vstack((data[indOff:,:], data[0:indOff,:]))
514        return data
515               
516    import matplotlib.pyplot as plt
517    import matplotlib.ticker as ticker
518    from exudyn.plot import DataArrayFromSensorList
519
520    mbs.PlotSensor(closeAll=True)
521    
522    #compute axial offset, to normalize results:
523    nodePos0 = mbs.GetSensorValues(sCable0Pos[0])
524    xOff = nodePos0[0]
525    maxXoff = 0.5*positionPulley2x
526    maxYoff = 0.1*r
527    # indOff = 0 #single data per element
528    # indOff2 = 0 #double data per element
529    correctXoffset = True
530    if abs(nodePos0[1]-r) > maxYoff or nodePos0[0] > maxXoff or nodePos0[0] < -0.1*maxXoff:
531        print('*****************')
532        print('warning: final position not at top of belt or too far away')
533        print('nodePos0=',nodePos0)
534        print('*****************')
535        xOff = 0
536        correctXoffset = False
537    else:
538        print('******************\nxOff=', xOff)
539            
540    
541    dataVel = DataArrayFromSensorList(mbs, sensorNumbers=sMidVel, positionList=positionListMid)
542    if correctXoffset:
543        dataVel=ShiftXoff(dataVel,xOff, reevingDict['totalLength'])
544    
545    # mbs.PlotSensor(sensorNumbers=[dataVel], components=0, labels=['axial velocity'], 
546    #            xLabel='axial position (m)', yLabel='velocity (m/s)')
547
548    #axial force over beam length:
549    dataForce = DataArrayFromSensorList(mbs, sensorNumbers=sAxialForce, positionList=positionList2Node)
550    if correctXoffset:
551        dataForce = ShiftXoff(dataForce,xOff, reevingDict['totalLength'])
552    # mbs.PlotSensor(sensorNumbers=[dataForce], components=0, labels=['axial force'], colorCodeOffset=2,
553    #            xLabel='axial position (m)', yLabel='axial force (N)')
554
555
556    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
557    #contact forces are stored (x/y) for every segment ==> put into consecutive array
558    contactForces =[[],[]] #these are the contact forces of the whole belt, but from both pulleys!
559    for i in range(len(sContactForce)):
560        contactForces[i] = np.zeros((len(sContactForce[i])*nSegments, 3)) #per row: [position, Fx, Fy]
561        for j, sensor in enumerate(sContactForce[i]):
562            values = mbs.GetSensorValues(sensor)
563            for k in range(nSegments):
564                row = j*nSegments + k
565                contactForces[i][row,0] = positionListSegments[row]
566                contactForces[i][row, 1:] = values[k*2:k*2+2]
567
568    contactForcesTotal = contactForces[0]
569    contactForcesTotal[:,1:] += contactForces[1][:,1:]
570
571    if correctXoffset:
572        contactForcesTotal = ShiftXoff(contactForcesTotal,-xOff, reevingDict['totalLength'])
573    #plot contact forces over beam length:
574    mbs.PlotSensor(sensorNumbers=[contactForcesTotal,contactForcesTotal], components=[0,1], labels=['tangential force','normal force'], 
575               xLabel='axial position (m)', yLabel='contact forces (N)', newFigure=True)
576    # mbs.PlotSensor(sensorNumbers=[contactForces[1],contactForces[1]], components=[0,1], labels=['tangential force','normal force'], 
577    #            xLabel='axial position (m)', yLabel='contact forces (N)', newFigure=False)
578    mbs.PlotSensor(sensorNumbers=[solDir+'contactForcesh5e-05n2s2cs40ALE1.txt'], components=[0,1], 
579               labels=['tangential force ALE','normal force ALE'], 
580                xLabel='axial position (m)', yLabel='contact forces (N)', colorCodeOffset=2, newFigure=False)
581    mbs.PlotSensor(sensorNumbers=[solDir+'contactForcesh5e-05n2s2cs40ALE0.txt'], components=[0,1], 
582               labels=['tangential force Ref','normal force Ref'], 
583                xLabel='axial position (m)', yLabel='contact forces (N)', colorCodeOffset=4, newFigure=False)
584
585    contactDisp =[[],[]] #slip and gap
586    for i in range(len(sContactDisp)):
587        contactDisp[i] = np.zeros((len(sContactDisp[i])*nSegments, 3)) #per row: [position, Fx, Fy]
588        for j, sensor in enumerate(sContactDisp[i]):
589            values = mbs.GetSensorValues(sensor)
590            for k in range(nSegments):
591                row = j*nSegments + k
592                contactDisp[i][row,0] = positionListSegments[row]
593                contactDisp[i][row, 1:] = values[k*2:k*2+2]
594
595    if correctXoffset:
596        contactDisp[0] = ShiftXoff(contactDisp[0],-xOff, reevingDict['totalLength'])
597        contactDisp[1] = ShiftXoff(contactDisp[1],-xOff, reevingDict['totalLength'])
598
599    if False:
600        mbs.PlotSensor(sensorNumbers=[contactDisp[0],contactDisp[0]], components=[0,1], labels=['slip','gap'], 
601                   xLabel='axial position (m)', yLabel='slip, gap (m)', newFigure=True)
602        mbs.PlotSensor(sensorNumbers=[contactDisp[1],contactDisp[1]], components=[0,1], labels=['slip','gap'], 
603                   xLabel='axial position (m)', yLabel='slip, gap (m)', newFigure=False)
604        mbs.PlotSensor(sensorNumbers=[solDir+'contactDisph5e-05n2s2cs40ALE1.txt'], components=[0,1], labels=['slipALE','gapALE'], 
605                   xLabel='axial position (m)', yLabel='slip, gap (m)', colorCodeOffset=2, newFigure=False)
606    
607
608
609    header  = ''
610    header += 'endTime='+str(tEnd)+'\n'
611    header += 'stepSize='+str(stepSize)+'\n'
612    header += 'nSegments='+str(nSegments)+'\n'
613    header += 'nANCFnodes='+str(nANCFnodes)+'\n'
614    header += 'contactStiffness='+str(contactStiffness)+'\n'
615    header += 'contactDamping='+str(contactDamping)+'\n'
616    header += 'frictionStiffness='+str(frictionStiffness)+'\n'
617    header += 'frictionVelocityPenalty='+str(frictionVelocityPenalty)+'\n'
618    header += 'dryFriction='+str(dryFriction)+'\n'
619    header += 'useALE='+str(useALE)+'\n'
620    fstr  = 'h'+str(stepSize)+'n'+str(int(nANCFnodes/60))+'s'+str(nSegments)+'cs'+str(int((contactStiffness/41800)))
621    fstr += 'ALE'+str(int(useALE))
622    #fstr += 'fs'+str(int((frictionStiffness/52300)))
623
624    #export solution:
625    contactDispSave = contactDisp[0]
626    contactDispSave[:,1:] += contactDisp[1][:,1:]
627
628    if False: #for saving
629        np.savetxt(solDir+'contactForces'+fstr+'.txt', contactForcesTotal, delimiter=',', 
630                   header='Exudyn: solution of belt drive, contact forces over belt length\n'+header, encoding=None)
631        np.savetxt(solDir+'contactDisp'+fstr+'.txt', contactDispSave, delimiter=',', 
632                   header='Exudyn: solution of belt drive, slip and gap over belt length\n'+header, encoding=None)
633
634
635    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
636
637    if False:
638        # mbs.PlotSensor(sensorNumbers=[sWheelRot[0], sWheelRot[1]], components=[2,2])#,sWheelRot[1]
639        mbs.PlotSensor(sensorNumbers=[sWheelRot[0], sWheelRot[1], 
640                                       solDir+'wheel0angVelALE.txt',solDir+'wheel1angVelALE.txt'], components=[2,2,2,2],
641                   colorCodeOffset=2)#,sWheelRot[1]
642    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++