ANCFgeneralContactCircle.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  ANCF cable element in contact with circles defined by GeneralContact
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2022-01-31
  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.utilities import * #includes itemInterface and rigidBodyUtilities
 15import exudyn.graphics as graphics
 16from exudyn.beams import *
 17import numpy as np
 18from math import sin, cos, sqrt, pi
 19
 20testIsActive = exu.sys.get('testIsActive', False)
 21SC = exu.SystemContainer()
 22mbs = SC.AddSystem()
 23
 24exu.Print('exudyn version=',exu.config.Version())
 25
 26# testIsActive = True
 27#background
 28rect = [-1,-1.5,3,1.5] #xmin,ymin,xmax,ymax
 29background0 = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[rect[0],rect[1],0, rect[2],rect[1],0, rect[2],rect[3],0, rect[0],rect[3],0, rect[0],rect[1],0]} #background
 30oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], 
 31                                   visualization=VObjectGround(graphicsData= [background0])))
 32nGround = mbs.AddNode(NodePointGround())
 33mCoordinateGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
 34
 35#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 36#contact
 37doImplicit = True
 38useContact = True
 39useFriction = True
 40dryFriction = 0.5
 41contactStiffness = 1e5
 42contactDamping = 1e-3*contactStiffness
 43
 44if useContact:
 45    gContact = mbs.AddGeneralContact()
 46    gContact.verboseMode = 1
 47    gContact.frictionProportionalZone = 1
 48    gContact.ancfCableUseExactMethod = False
 49    gContact.ancfCableNumberOfContactSegments = 8
 50    ssx = 16#32 #search tree size
 51    ssy = 8#32 #search tree size
 52    ssz = 1 #search tree size
 53    gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssz])
 54    #gContact.SetSearchTreeBox(pMin=np.array([-1,-1,-1]), pMax=np.array([4,1,1]))
 55
 56torque=-20
 57#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 58#wheel:
 59dampingWheel1 = 1  #5 #add breaking torque, to limit velocity
 60#cable:
 61
 62numberOfElements = 8    # per section
 63curvedRefConf=False     # this flag could initialize the elements to be produced curved -> not suitable for belt drive!
 64L=2                     # length of ANCF element in m
 65E=1e10                  # Young's modulus of ANCF element in N/m^2
 66rhoBeam=1000            # density of ANCF element in kg/m^3
 67b=0.002                 # width of rectangular ANCF element in m
 68h=0.002                 # height of rectangular ANCF element in m
 69A=b*h                   # cross sectional area of ANCF element in m^2
 70I=b*h**3/12             # second moment of area of ANCF element in m^4
 71dEI = 0*1e-3*E*I
 72dEA = 1e-2*E*A
 73# f=3*E*I/L**2            # tip load applied to ANCF element in N
 74g=-9.81
 75dimZ = b #z.dimension
 76preStretch=-0.002
 77# exu.Print("load f="+str(f))
 78# exu.Print("EI="+str(E*I))
 79
 80# nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 81# mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 82
 83cableTemplate = Cable2D(#length = L / nElements, #set in GenerateBeamElementsAlongLine(...)
 84                        massPerLength = rhoBeam*A,
 85                        bendingStiffness = E*I,
 86                        axialStiffness = E*A,
 87                        bendingDamping = dEI,
 88                        axialDamping = dEA,
 89                        referenceAxialStrain = preStretch, #prestretch
 90                        #nodeNumbers = [0, 0], #will be filled in GenerateBeamElementsAlongLine(...)
 91                        visualization=VCable2D(drawHeight=2*h),
 92                        )
 93exu.Print("pre-stretch force=", preStretch*E*A)
 94exu.Print("beam mass per length=", rhoBeam*A)
 95#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 96#create belt drive:
 97distanceWheels = 2 #distance of wheel centers
 98wheelCenter0 = np.array([0,0,0])
 99wheelCenter1 = np.array([distanceWheels,0,0])
100
101rWheel0 = 0.5
102rWheel1 = 0.5
103mWheel = 2
104
105yAxis = np.array([0,1.,0])
106ancfList=[]
107
108if True: #add ANCF cable elements
109    startAngle = -pi
110    arcAngle = -pi
111    positionOfNode0 = wheelCenter0-rWheel0*yAxis # starting point of line
112    ancf=GenerateCircularArcANCFCable2D(mbs, positionOfNode0, 
113                                        rWheel0, startAngle, arcAngle, numberOfElements, 
114                                        cableTemplate,
115                                        massProportionalLoad = [0,g,0], #optionally add gravity
116                                        #fixedConstraintsNode0 = [1,1,1,1], #add constraints for pos and rot (r'_y)
117                                        #fixedConstraintsNode1 = [1,1,1,1],
118                                        setCurvedReferenceConfiguration=curvedRefConf, 
119                                        )
120    ancfList+=[ancf]
121    beamInfo=GenerateBeamElementsAlongLine(mbs,
122                                         ancf[3][-1], wheelCenter1+rWheel1*yAxis,
123                                         numberOfElements,
124                                         beamTemplate=cableTemplate, #this defines the beam element properties
125                                         gravity = [0,g,0], #optionally add gravity
126                                         nodeNumberStart=ancf[0][-1]
127                                         )
128    ancf = [beamInfo['nodes'], beamInfo['elements'], beamInfo['loads'], beamInfo['nodePositions'], beamInfo['constraintObjects']] #as the list of GenerateCircularArcANCFCable2D
129    ancfList+=[ancf]
130
131    startAngle = 0
132    arcAngle = -pi
133    ancf=GenerateCircularArcANCFCable2D(mbs, ancf[3][-1], 
134                                        rWheel1, startAngle, arcAngle, numberOfElements, 
135                                        cableTemplate,
136                                        massProportionalLoad = [0,g,0], #optionally add gravity
137                                        setCurvedReferenceConfiguration=curvedRefConf, 
138                                        nodeNumber0=ancf[0][-1]
139                                        )
140    ancfList+=[ancf]
141    beamInfo=GenerateBeamElementsAlongLine(mbs,
142                                         ancf[3][-1], ancfList[0][3][0],
143                                         numberOfElements,
144                                         beamTemplate=cableTemplate, #this defines the beam element properties
145                                         gravity = [0,g,0], #optionally add gravity
146                                         nodeNumberStart=ancf[0][-1],
147                                         nodeNumberEnd=ancfList[0][0][0]
148                                         )
149    ancf = [beamInfo['nodes'], beamInfo['elements'], beamInfo['loads'], beamInfo['nodePositions'], beamInfo['constraintObjects']] #as the list of GenerateCircularArcANCFCable2D
150    ancfList+=[ancf]
151
152if not testIsActive: 
153    #add sensor for one node, showing moving coordinates
154    sensorsNode = []
155    for i, aList in enumerate(ancfList):
156        sensorsNode += [mbs.AddSensor(SensorNode(nodeNumber=aList[0][0], #fileName='solutionNode'+str(i)+'.txt', 
157                                                 storeInternal=True,outputVariableType=exu.OutputVariableType.Position))]
158
159
160#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
161sAngVel=[]
162#add contact:
163if useContact:
164
165
166    halfHeight = 0.5*h*0
167    wheels = [{'center':wheelCenter0, 'radius':rWheel0-halfHeight, 'mass':mWheel}, 
168              {'center':wheelCenter1, 'radius':rWheel1-halfHeight, 'mass':mWheel}, ]
169
170    for i, wheel in enumerate(wheels):
171        r = wheel['radius']
172        p = wheel['center']
173        mass = wheel['mass']
174        rot0 = 0 #initial rotation
175        pRef = [p[0], p[1], rot0]
176        gList = [graphics.Cylinder(vAxis=[0,0,dimZ], radius=r*0.99, #draw smaller to see cable element 
177                                      color= graphics.color.dodgerblue, nTiles=32),
178                 graphics.Arrow(pAxis=[0,0,0.02*r], vAxis=[r,0,0], radius=0.02*r, color=graphics.color.orange)]
179
180        omega0 = 0 #initial angular velocity
181        v0 = np.array([0,0,omega0]) 
182
183        RBinertia = InertiaCylinder(mass/(r**2*np.pi*b), b, r, axis=2)
184
185        nMass = mbs.AddNode(NodeRigidBody2D(referenceCoordinates=pRef, initialVelocities=v0,
186                                            visualization=VNodeRigidBody2D(drawSize=dimZ*2)))
187        oMass = mbs.AddObject(ObjectRigidBody2D(mass=RBinertia.mass, inertia=RBinertia.GetInertia6D()[2],
188                                                nodeNumber=nMass, visualization=
189                                                VObjectRigidBody2D(graphicsData=gList)))
190        mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
191        mGroundWheel = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=p))
192        frictionMaterialIndex=0
193
194        mbs.AddObject(RevoluteJoint2D(markerNumbers=[mGroundWheel, mNode]))
195        
196        if i == 0:
197            mbs.AddLoad(LoadTorqueVector(markerNumber=mNode, loadVector=[0,0,torque]))
198        if i == 1:
199            mCoordinateWheel = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMass, coordinate=2))
200            mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mCoordinateGround, mCoordinateWheel], 
201                                                 damping=dampingWheel1, 
202                                                 visualization=VCoordinateSpringDamper(show=False)))
203
204        gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=contactStiffness, 
205                                     contactDamping=contactDamping, frictionMaterialIndex=frictionMaterialIndex)
206        
207        if not testIsActive: 
208            sAngVel += [mbs.AddSensor(SensorNode(nodeNumber=nMass, #fileName='solution/wheel'+str(i)+'angVel.txt',
209                                                 storeInternal=True, outputVariableType=exu.OutputVariableType.AngularVelocity))]
210
211    #generate list of all cable elements:
212    allCables = []
213    for ancf in ancfList:
214        allCables += ancf[1]
215        
216    #add all cable elements to contact
217    for oIndex in allCables:
218        gContact.AddANCFCable(objectIndex=oIndex, halfHeight=halfHeight, #halfHeight should be h/2, but then cylinders should be smaller
219                              contactStiffness=contactStiffness, contactDamping=contactDamping, 
220                              frictionMaterialIndex=0)
221
222    #create matrix of material interaction (in this case, only 1x1):
223    frictionMatrix = np.zeros((1,1))
224    frictionMatrix[0,0]=int(useFriction)*dryFriction
225    gContact.SetFrictionPairings(frictionMatrix)
226    #gContact.verboseMode=2
227
228mbs.Assemble()
229
230simulationSettings = exu.SimulationSettings() #takes currently set values or default values
231
232tEnd = 0.1
233h = 1e-3
234
235simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
236simulationSettings.solution.file.name = 'solution/coordinatesSolution.txt'
237
238if not testIsActive:
239    tEnd = 0.75
240    simulationSettings.solution.file.write = True
241    simulationSettings.solution.file.writePeriod = 0.005
242else:
243    simulationSettings.solution.file.write = False
244
245simulationSettings.solution.sensors.writePeriod = 0.001
246#simulationSettings.show.computationTime = True
247simulationSettings.parallel.numberOfThreads = 1 #use 4 to speed up for > 100 ANCF elements
248simulationSettings.show.statistics = True
249
250doDynamic = True
251simulationSettings.timeIntegration.endTime = tEnd
252simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
253simulationSettings.timeIntegration.stepInformation= 3+128+256 #show step reduction and increase
254
255simulationSettings.timeIntegration.verboseMode = 1 #otherwise, load steps are shown ...
256simulationSettings.timeIntegration.newton.useModifiedNewton = True
257
258SC.visualizationSettings.view0.scene.drawWorldBasis=True
259SC.visualizationSettings.nodes.show = True
260SC.visualizationSettings.nodes.defaultSize = h*20
261SC.visualizationSettings.loads.show = False
262
263SC.visualizationSettings.contour.outputVariableComponent=0
264SC.visualizationSettings.contour.outputVariable=exu.OutputVariableType.ForceLocal
265
266#visualize contact:
267if False:
268    SC.visualizationSettings.contact.showSearchTree =True
269    SC.visualizationSettings.contact.showSearchTreeCells =True
270    SC.visualizationSettings.contact.showBoundingBoxes = True
271
272if not testIsActive: 
273    SC.renderer.Start()
274    SC.renderer.DoIdleTasks()
275
276mbs.SolveDynamic(simulationSettings) #183 Newton iterations, 0.114 seconds
277
278
279if not testIsActive and False:
280    SC.visualizationSettings.general.autoFitScene = False
281    SC.visualizationSettings.general.graphicsUpdateInterval=0.02
282    
283    sol = LoadSolutionFile('solution/coordinatesSolution.txt', safeMode=True)#, maxRows=100)
284    exu.Print('start SolutionViewer')
285    mbs.SolutionViewer(sol)
286
287
288if not testIsActive: 
289    SC.renderer.DoIdleTasks()
290    SC.renderer.Stop() #safely close rendering window!
291    
292    if len(sAngVel) != 0:
293        
294        mbs.PlotSensor(sensorNumbers=[sAngVel[0],sAngVel[1]], components=2, closeAll=True)
295        mbs.PlotSensor(sensorNumbers=sensorsNode, componentsX=0, components=1, 
296                   xLabel='PositionX', newFigure=True, title='trajectories of 4 nodes')
297        
298#print representative result:
299posNode0 = mbs.GetNodeOutput(ancfList[0][0][0], variableType=exu.OutputVariableType.Position)
300exu.Print('node0 pos: ',posNode0) #[-0.0922746  -0.48937754  0.        ] 
301sol = posNode0[0] + posNode0[1]
302exu.Print('ANCFgeneralContactCircle sol=',sol)
303
304exu.sys['testResult'] = sol