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