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 #only import if it does not conflict
 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                                          visualization = VNode1D(show = False)))#, color = [0.,0.,0.,1.])
148    mALE = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nALE, coordinate=0,
149                                              visualization = {'show':True})) #ALE velocity
150    nodesANCF = [-1,-1, nALE]
151
152    #Constraint for eulerian coordinate
153    oCCvALE=mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mALE], offset=0,
154                                                velocityLevel = False,
155                                                visualization=VCoordinateConstraint(show=False)))
156
157cableTemplate = ANCFElementType(#physicsLength = L / nElements, #set in GenerateStraightLineANCFCable2D(...)
158                        nodeNumbers = nodesANCF,
159                        physicsMassPerLength = rhoA,
160                        physicsBendingStiffness = EI,
161                        physicsAxialStiffness = EA,
162                        physicsBendingDamping = dEI,
163                        physicsReferenceAxialStrain = preStretch,
164                        physicsReferenceCurvature = 0.,
165                        useReducedOrderIntegration = 2,
166                        strainIsRelativeToReference = False, #would cause reference configuration to be precurved
167                        visualization=VCable2D(drawHeight=hc),
168                        )
169
170if useALE:
171    cableTemplate.physicsAddALEvariation = False
172
173ancf = PointsAndSlopes2ANCFCable2D(mbs, reevingDict['ancfPointsSlopes'],
174                                   reevingDict['elementLengths'],
175                                   cableTemplate, massProportionalLoad=gVec,
176                                   fixedConstraintsNode0=[1*staticEqulibrium,0,0,0], #fixedConstraintsNode1=[1,1,1,1],
177                                   #elementCurvatures  = elementCurvatures, #do not set this, will cause inhomogeneous curvatures
178                                   firstNodeIsLastNode=True, graphicsSizeConstraints=0.01)
179
180lElem = reevingDict['totalLength'] / nANCFnodes
181cFact=b*lElem/nSegments #stiffness shall be per area, but is applied at every segment
182
183contactStiffness*=40*cFact
184contactDamping = 40*2000*cFact #according to Dufva 2008 paper ... seems also to be used in 2013 PEchstein Gerstmayr
185frictionStiffness = 50e8*cFact #1e7 converges good; 1e8 is already quite accurate
186massSegment = rhoA*lElem/nSegments
187frictionVelocityPenalty = 10*sqrt(frictionStiffness*massSegment) #bristle damping; should be adjusted to reduce vibrations induced by bristle model
188
189if useFrictionStiffness:
190    frictionStiffness*=0.1
191else:
192    frictionStiffness*=0
193    frictionVelocityPenalty*= 2
194
195#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
196#create sensors for all nodes
197sMidVel = []
198sAxialForce = []
199sCable0Pos = []
200# sObjectDisp =[]
201
202ancfNodes = ancf[0]
203ancfObjects = ancf[1]
204positionList2Node = [] #axial position at x=0 and x=0.5*lElem
205positionListMid = [] #axial position at midpoint of element
206positionListSegments = [] #axial position at midpoint of segments
207currentPosition = 0 #is increased at every iteration
208for i,obj in enumerate(ancfObjects):
209    lElem = reevingDict['elementLengths'][i]
210    positionList2Node += [currentPosition, currentPosition + 0.5*lElem]
211    positionListMid += [currentPosition + 0.5*lElem]
212
213    for j in range(nSegments):
214        segPos = (j+0.5)*lElem/nSegments + currentPosition
215        positionListSegments += [segPos]
216    currentPosition += lElem
217
218    sAxialForce += [mbs.AddSensor(SensorBody(bodyNumber = obj,
219                                              storeInternal=True,
220                                              localPosition=[0.*lElem,0,0],
221                                              outputVariableType=exu.OutputVariableType.ForceLocal))]
222    sAxialForce += [mbs.AddSensor(SensorBody(bodyNumber = obj,
223                                              storeInternal=True,
224                                              localPosition=[0.5*lElem,0,0],
225                                              outputVariableType=exu.OutputVariableType.ForceLocal))]
226    sMidVel += [mbs.AddSensor(SensorBody(bodyNumber = obj,
227                                          storeInternal=True,
228                                          localPosition=[0.5*lElem,0,0], #0=at left node
229                                          outputVariableType=exu.OutputVariableType.VelocityLocal))]
230    sCable0Pos += [mbs.AddSensor(SensorBody(bodyNumber = obj,
231                                            storeInternal=True,
232                                            localPosition=[0.*lElem,0,0],
233                                            outputVariableType=exu.OutputVariableType.Position))]
234    # sObjectDisp += [mbs.AddSensor(SensorBody(bodyNumber = obj,
235    #                                           storeInternal=True,
236    #                                           localPosition=[0.5*lElem,0,0],
237    #                                           outputVariableType=exu.OutputVariableType.Displacement))]
238
239
240
241#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
242#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
243#add contact:
244if useContact:
245
246    contactObjects = [[],[]] #list of contact objects
247
248    dimZ= 0.01 #for drawing
249    sWheelRot = [] #sensors for angular velocity
250
251    nMassList = []
252    wheelSprings = [] #for static computation
253    for i, wheel in enumerate(circleList):
254        p = [wheel[0][0], wheel[0][1], 0] #position of wheel center
255        r = wheel[1]
256
257        rot0 = 0 #initial rotation
258        pRef = [p[0], p[1], rot0]
259        gList = [graphics.Cylinder(pAxis=[0,0,-dimZ],vAxis=[0,0,-dimZ], radius=r,
260                                      color= graphics.color.dodgerblue, nTiles=64),
261                 graphics.Arrow(pAxis=[0,0,0], vAxis=[-0.9*r,0,0], radius=0.01*r, color=graphics.color.orange),
262                 graphics.Arrow(pAxis=[0,0,0], vAxis=[0.9*r,0,0], radius=0.01*r, color=graphics.color.orange)]
263
264        omega0 = 0 #initial angular velocity
265        v0 = np.array([0,0,omega0])
266
267        nMass = mbs.AddNode(NodeRigidBody2D(referenceCoordinates=pRef, initialVelocities=v0,
268                                            visualization=VNodeRigidBody2D(drawSize=dimZ*2)))
269        nMassList += [nMass]
270        oMass = mbs.AddObject(ObjectRigidBody2D(physicsMass=wheelMass, physicsInertia=wheelInertia,
271                                                nodeNumber=nMass, visualization=
272                                                VObjectRigidBody2D(graphicsData=gList)))
273        mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
274        mGroundWheel = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=p, visualization = VMarkerBodyRigid(show = False)))
275
276        #mbs.AddObject(RevoluteJoint2D(markerNumbers=[mGroundWheel, mNode], visualization=VRevoluteJoint2D(show=False)))
277
278        mCoordinateWheelX = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=0))
279        mCoordinateWheelY = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=1))
280        constraintX = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheelX],
281                                                 visualization=VCoordinateConstraint(show = False)))
282        constraintY = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheelY],
283                                                 visualization=VCoordinateConstraint(show = False)))
284        if i==0:
285            constraintPulleyLeftX = constraintX
286
287        if True:
288
289            sWheelRot += [mbs.AddSensor(SensorNode(nodeNumber=nMass,
290                                                   storeInternal=True,
291                                                   fileName='solutionDelete/wheel'+str(i)+'angVel.txt',
292                                                   outputVariableType=exu.OutputVariableType.AngularVelocity))]
293        tdisplacement = 0.05
294
295
296        def UFvelocityDrive(mbs, t, itemNumber, lOffset): #time derivative of UFoffset
297            if t < tAccStart:
298                v = 0
299            if t >= tAccStart and t < tAccEnd:
300                v = omegaFinal/(tAccEnd-tAccStart)*(t-tAccStart)
301            elif t >= tAccEnd:
302                v = omegaFinal
303            return v
304
305        if doDynamic:
306            if i == 0:
307                mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
308                velControl = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheel],
309                                                    velocityLevel=True, offsetUserFunction_t= UFvelocityDrive,
310                                                    visualization=VCoordinateConstraint(show = False)))#UFvelocityDrive
311            if i == 1:
312                mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
313                mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mCoordinateGround, mCoordinateWheel],
314                                                     damping = rotationDampingWheels,
315                                                     visualization=VCoordinateSpringDamper(show = False)))
316
317                #this is used for times > 1 in order to see influence of torque step in Wheel1
318                def UFforce(mbs, t, load):
319                    tau = 0.
320                    tau +=  25.*(SmoothStep(t, 1., 1.5, 0., 1.) - SmoothStep(t, 3.5, 4., 0., 1.))
321                    #tau += 16.*(SmoothStep(t, 5, 5.5, 0., 1.) - SmoothStep(t, 7.5, 8., 0., 1.))
322                    return -tau
323
324                mbs.AddLoad(LoadCoordinate(markerNumber=mCoordinateWheel,
325                                           load = 0, loadUserFunction = UFforce))
326
327        if staticEqulibrium:
328            mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
329            csd = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround, mCoordinateWheel],
330                                                     visualization=VCoordinateConstraint(show = False)))
331            wheelSprings += [csd]
332
333
334        cableList = ancf[1]
335        mCircleBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass))
336        #mCircleBody = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
337        for k in range(len(cableList)):
338            initialGapList = [0.1]*nSegments + [-2]*(nSegments) + [0]*(nSegments) #initial gap of 0., isStick (0=slip, +-1=stick, -2 undefined initial state), lastStickingPosition (0)
339
340            mCable = mbs.AddMarker(MarkerBodyCable2DShape(bodyNumber=cableList[k],
341                                                          numberOfSegments = nSegments, verticalOffset=-0*hc/2))
342            nodeDataContactCable = mbs.AddNode(NodeGenericData(initialCoordinates=initialGapList,
343                                                               numberOfDataCoordinates=nSegments*(1+2) ))
344
345            co = mbs.AddObject(ObjectContactFrictionCircleCable2D(markerNumbers=[mCircleBody, mCable], nodeNumber = nodeDataContactCable,
346                                                     numberOfContactSegments=nSegments,
347                                                     contactStiffness = contactStiffness,
348                                                     contactDamping=contactDamping,
349                                                     frictionVelocityPenalty = frictionVelocityPenalty,
350                                                     frictionStiffness = frictionStiffness,
351                                                     frictionCoefficient=int(useFriction)*dryFriction,
352                                                     circleRadius = r,
353                                                     visualization=VObjectContactFrictionCircleCable2D(showContactCircle=False)))
354            contactObjects[i] += [co]
355
356
357
358#+++++++++++++++++++++++++++++++++++++++++++
359#create list of sensors for contact
360sContactDisp = [[],[]]
361sContactForce = [[],[]]
362for i in range(len(contactObjects)):
363    for obj in contactObjects[i]:
364        sContactForce[i] += [mbs.AddSensor(SensorObject(objectNumber = obj,
365                                                        storeInternal=True,
366                                                        outputVariableType=exu.OutputVariableType.ForceLocal))]
367        sContactDisp[i] += [mbs.AddSensor(SensorObject(objectNumber = obj,
368                                                        storeInternal=True,
369                                                        outputVariableType=exu.OutputVariableType.Coordinates))]
370
371
372mbs.Assemble()
373
374
375simulationSettings = exu.SimulationSettings() #takes currently set values or default values
376
377simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
378simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/testCoords.txt'
379
380simulationSettings.solutionSettings.writeSolutionToFile = True
381simulationSettings.solutionSettings.solutionWritePeriod = 0.002
382simulationSettings.solutionSettings.sensorsWritePeriod = 0.001
383simulationSettings.parallel.numberOfThreads = 1 #use 4 to speed up for > 100 ANCF elements
384
385simulationSettings.timeIntegration.endTime = tEnd
386simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
387simulationSettings.timeIntegration.stepInformation= 255
388
389simulationSettings.timeIntegration.verboseMode = 1
390
391simulationSettings.timeIntegration.newton.useModifiedNewton = True
392#simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
393#simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
394
395simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1e-3
396simulationSettings.timeIntegration.discontinuous.maxIterations = discontinuousIterations #3
397
398simulationSettings.displayStatistics = True
399simulationSettings.displayComputationTime = False
400
401
402SC.visualizationSettings.general.circleTiling = 24
403SC.visualizationSettings.loads.show=False
404SC.visualizationSettings.sensors.show=False
405SC.visualizationSettings.markers.show=False
406SC.visualizationSettings.nodes.defaultSize = 0.002
407SC.visualizationSettings.openGL.multiSampling = 4
408SC.visualizationSettings.openGL.lineWidth = 2
409SC.visualizationSettings.window.renderWindowSize = [1920,1080]
410
411SC.visualizationSettings.connectors.showContact = True
412SC.visualizationSettings.contact.contactPointsDefaultSize = 0.0002
413SC.visualizationSettings.contact.showContactForces = True
414SC.visualizationSettings.contact.contactForcesFactor = 0.005
415
416
417if True:
418    SC.visualizationSettings.bodies.beams.axialTiling = 1
419    SC.visualizationSettings.bodies.beams.drawVertical = True
420    SC.visualizationSettings.bodies.beams.drawVerticalLines = True
421
422    SC.visualizationSettings.contour.outputVariableComponent=0
423    SC.visualizationSettings.contour.outputVariable=exu.OutputVariableType.ForceLocal
424
425    SC.visualizationSettings.bodies.beams.drawVerticalFactor = 0.0003
426    SC.visualizationSettings.bodies.beams.drawVerticalOffset = -220
427
428    SC.visualizationSettings.bodies.beams.reducedAxialInterploation = True
429
430    # SC.visualizationSettings.contour.outputVariable=exu.OutputVariableType.VelocityLocal
431    # SC.visualizationSettings.bodies.beams.drawVerticalFactor = -0.25
432    # SC.visualizationSettings.bodies.beams.drawVerticalOffset = 0.30
433    # SC.visualizationSettings.bodies.beams.reducedAxialInterploation = False
434
435
436if useGraphics:
437    exu.StartRenderer()
438    #mbs.WaitForUserToContinue()
439
440#simulationSettings.staticSolver.newton.absoluteTolerance = 1e-10
441simulationSettings.staticSolver.adaptiveStep = True
442simulationSettings.staticSolver.loadStepGeometric = False;
443simulationSettings.staticSolver.loadStepGeometricRange=1e4
444simulationSettings.staticSolver.numberOfLoadSteps = 10
445#simulationSettings.staticSolver.useLoadFactor = False
446simulationSettings.staticSolver.stabilizerODE2term = 1e5*10
447simulationSettings.staticSolver.newton.relativeTolerance = 1e-6
448simulationSettings.staticSolver.newton.absoluteTolerance = 1e-6
449
450if staticEqulibrium: #precompute static equilibrium
451    mbs.SetObjectParameter(velControl, 'activeConnector', False)
452
453    for i in range(len(contactObjects)):
454        for obj in contactObjects[i]:
455            mbs.SetObjectParameter(obj, 'frictionCoefficient', 0.)
456            mbs.SetObjectParameter(obj, 'frictionStiffness', 1e-8) #do not set to zero, as it needs to do some initialization...
457
458    # simulationSettings.solutionSettings.appendToFile=False
459    mbs.SolveStatic(simulationSettings, updateInitialValues=True)
460    # simulationSettings.solutionSettings.appendToFile=True
461
462    #check total force on support, expect: supportLeftX \approx 2*preStretch*EA
463    supportLeftX = mbs.GetObjectOutput(constraintPulleyLeftX,variableType=exu.OutputVariableType.Force)
464    print('Force x in support of left pulley = ', supportLeftX)
465    print('Belt pre-tension=', preStretch*EA)
466
467    for i in range(len(contactObjects)):
468        for obj in contactObjects[i]:
469            mbs.SetObjectParameter(obj, 'frictionCoefficient', dryFriction)
470            mbs.SetObjectParameter(obj, 'frictionStiffness', frictionStiffness)
471
472    # useALE = False
473    for coordinateConstraint in ancf[4]: #release constraints for dynamic Solver
474        if not useALE: #except ALE constraint
475            mbs.SetObjectParameter(coordinateConstraint, 'activeConnector', False)
476
477    if useALE:
478        mbs.SetObjectParameter(oCCvALE, 'activeConnector', False) #do not fix ALE coordinate any more
479
480
481    mbs.SetObjectParameter(velControl, 'activeConnector', True)
482    for csd in wheelSprings:
483        mbs.SetObjectParameter(csd, 'activeConnector', False)
484
485if True:
486    mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.TrapezoidalIndex2) #183 Newton iterations, 0.114 seconds
487#mbs.SolveDynamic(simulationSettings)
488
489if useGraphics and True:
490    SC.visualizationSettings.general.autoFitScene = False
491    SC.visualizationSettings.general.graphicsUpdateInterval=0.02
492
493    sol = LoadSolutionFile('solution/testCoords.txt', safeMode=True)#, maxRows=100)
494    mbs.SolutionViewer(sol)
495
496
497if useGraphics:
498    SC.WaitForRenderEngineStopFlag()
499    exu.StopRenderer() #safely close rendering window!
500
501#%%++++++++++++++++++++++++++++++++++++++++
502if True:
503    solDir = 'solutionDelete/'
504    #shift data depending on axial position by subtracting xOff; put negative x values+shiftValue to end of array
505    def ShiftXoff(data, xOff, shiftValue):
506        indOff = 0
507        n = data.shape[0]
508        data[:,0] -= xOff
509        for i in range(n):
510           if data[i,0] < 0:
511               indOff+=1
512               data[i,0] += shiftValue
513        print('indOff=', indOff)
514        data = np.vstack((data[indOff:,:], data[0:indOff,:]))
515        return data
516
517    import matplotlib.pyplot as plt
518    import matplotlib.ticker as ticker
519    from exudyn.plot import DataArrayFromSensorList
520
521    mbs.PlotSensor(closeAll=True)
522
523    #compute axial offset, to normalize results:
524    nodePos0 = mbs.GetSensorValues(sCable0Pos[0])
525    xOff = nodePos0[0]
526    maxXoff = 0.5*positionPulley2x
527    maxYoff = 0.1*r
528    # indOff = 0 #single data per element
529    # indOff2 = 0 #double data per element
530    correctXoffset = True
531    if abs(nodePos0[1]-r) > maxYoff or nodePos0[0] > maxXoff or nodePos0[0] < -0.1*maxXoff:
532        print('*****************')
533        print('warning: final position not at top of belt or too far away')
534        print('nodePos0=',nodePos0)
535        print('*****************')
536        xOff = 0
537        correctXoffset = False
538    else:
539        print('******************\nxOff=', xOff)
540
541
542    dataVel = DataArrayFromSensorList(mbs, sensorNumbers=sMidVel, positionList=positionListMid)
543    if correctXoffset:
544        dataVel=ShiftXoff(dataVel,xOff, reevingDict['totalLength'])
545
546    # mbs.PlotSensor(sensorNumbers=[dataVel], components=0, labels=['axial velocity'],
547    #            xLabel='axial position (m)', yLabel='velocity (m/s)')
548
549    #axial force over beam length:
550    dataForce = DataArrayFromSensorList(mbs, sensorNumbers=sAxialForce, positionList=positionList2Node)
551    if correctXoffset:
552        dataForce = ShiftXoff(dataForce,xOff, reevingDict['totalLength'])
553    # mbs.PlotSensor(sensorNumbers=[dataForce], components=0, labels=['axial force'], colorCodeOffset=2,
554    #            xLabel='axial position (m)', yLabel='axial force (N)')
555
556
557    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
558    #contact forces are stored (x/y) for every segment ==> put into consecutive array
559    contactForces =[[],[]] #these are the contact forces of the whole belt, but from both pulleys!
560    for i in range(len(sContactForce)):
561        contactForces[i] = np.zeros((len(sContactForce[i])*nSegments, 3)) #per row: [position, Fx, Fy]
562        for j, sensor in enumerate(sContactForce[i]):
563            values = mbs.GetSensorValues(sensor)
564            for k in range(nSegments):
565                row = j*nSegments + k
566                contactForces[i][row,0] = positionListSegments[row]
567                contactForces[i][row, 1:] = values[k*2:k*2+2]
568
569    contactForcesTotal = contactForces[0]
570    contactForcesTotal[:,1:] += contactForces[1][:,1:]
571
572    if correctXoffset:
573        contactForcesTotal = ShiftXoff(contactForcesTotal,-xOff, reevingDict['totalLength'])
574    #plot contact forces over beam length:
575    mbs.PlotSensor(sensorNumbers=[contactForcesTotal,contactForcesTotal], components=[0,1], labels=['tangential force','normal force'],
576               xLabel='axial position (m)', yLabel='contact forces (N)', newFigure=True)
577    # mbs.PlotSensor(sensorNumbers=[contactForces[1],contactForces[1]], components=[0,1], labels=['tangential force','normal force'],
578    #            xLabel='axial position (m)', yLabel='contact forces (N)', newFigure=False)
579    mbs.PlotSensor(sensorNumbers=[solDir+'contactForcesh5e-05n2s2cs40ALE1.txt'], components=[0,1],
580               labels=['tangential force ALE','normal force ALE'],
581                xLabel='axial position (m)', yLabel='contact forces (N)', colorCodeOffset=2, newFigure=False)
582    mbs.PlotSensor(sensorNumbers=[solDir+'contactForcesh5e-05n2s2cs40ALE0.txt'], components=[0,1],
583               labels=['tangential force Ref','normal force Ref'],
584                xLabel='axial position (m)', yLabel='contact forces (N)', colorCodeOffset=4, newFigure=False)
585
586    contactDisp =[[],[]] #slip and gap
587    for i in range(len(sContactDisp)):
588        contactDisp[i] = np.zeros((len(sContactDisp[i])*nSegments, 3)) #per row: [position, Fx, Fy]
589        for j, sensor in enumerate(sContactDisp[i]):
590            values = mbs.GetSensorValues(sensor)
591            for k in range(nSegments):
592                row = j*nSegments + k
593                contactDisp[i][row,0] = positionListSegments[row]
594                contactDisp[i][row, 1:] = values[k*2:k*2+2]
595
596    if correctXoffset:
597        contactDisp[0] = ShiftXoff(contactDisp[0],-xOff, reevingDict['totalLength'])
598        contactDisp[1] = ShiftXoff(contactDisp[1],-xOff, reevingDict['totalLength'])
599
600    if False:
601        mbs.PlotSensor(sensorNumbers=[contactDisp[0],contactDisp[0]], components=[0,1], labels=['slip','gap'],
602                   xLabel='axial position (m)', yLabel='slip, gap (m)', newFigure=True)
603        mbs.PlotSensor(sensorNumbers=[contactDisp[1],contactDisp[1]], components=[0,1], labels=['slip','gap'],
604                   xLabel='axial position (m)', yLabel='slip, gap (m)', newFigure=False)
605        mbs.PlotSensor(sensorNumbers=[solDir+'contactDisph5e-05n2s2cs40ALE1.txt'], components=[0,1], labels=['slipALE','gapALE'],
606                   xLabel='axial position (m)', yLabel='slip, gap (m)', colorCodeOffset=2, newFigure=False)
607
608
609
610    header  = ''
611    header += 'endTime='+str(tEnd)+'\n'
612    header += 'stepSize='+str(stepSize)+'\n'
613    header += 'nSegments='+str(nSegments)+'\n'
614    header += 'nANCFnodes='+str(nANCFnodes)+'\n'
615    header += 'contactStiffness='+str(contactStiffness)+'\n'
616    header += 'contactDamping='+str(contactDamping)+'\n'
617    header += 'frictionStiffness='+str(frictionStiffness)+'\n'
618    header += 'frictionVelocityPenalty='+str(frictionVelocityPenalty)+'\n'
619    header += 'dryFriction='+str(dryFriction)+'\n'
620    header += 'useALE='+str(useALE)+'\n'
621    fstr  = 'h'+str(stepSize)+'n'+str(int(nANCFnodes/60))+'s'+str(nSegments)+'cs'+str(int((contactStiffness/41800)))
622    fstr += 'ALE'+str(int(useALE))
623    #fstr += 'fs'+str(int((frictionStiffness/52300)))
624
625    #export solution:
626    contactDispSave = contactDisp[0]
627    contactDispSave[:,1:] += contactDisp[1][:,1:]
628
629    if False: #for saving
630        np.savetxt(solDir+'contactForces'+fstr+'.txt', contactForcesTotal, delimiter=',',
631                   header='Exudyn: solution of belt drive, contact forces over belt length\n'+header, encoding=None)
632        np.savetxt(solDir+'contactDisp'+fstr+'.txt', contactDispSave, delimiter=',',
633                   header='Exudyn: solution of belt drive, slip and gap over belt length\n'+header, encoding=None)
634
635
636    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
637
638    if False:
639        # mbs.PlotSensor(sensorNumbers=[sWheelRot[0], sWheelRot[1]], components=[2,2])#,sWheelRot[1]
640        mbs.PlotSensor(sensorNumbers=[sWheelRot[0], sWheelRot[1],
641                                       solDir+'wheel0angVelALE.txt',solDir+'wheel1angVelALE.txt'], components=[2,2,2,2],
642                   colorCodeOffset=2)#,sWheelRot[1]
643    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++