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