flexiblePendulumANCF.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  simple flexible pendulum using 2D ANCF elements
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2022-06-28
  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 sys
 14sys.exudynFast = True
 15import exudyn as exu
 16from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 17from exudyn.beams import GenerateBeamElementsAlongLine
 18import exudyn
 19import exudyn.graphics as graphics
 20
 21import numpy as np
 22#from math import sqrt, sin, cos
 23
 24#%%++++++++++++++++++++++++++++++++++++++++
 25useGraphics = True
 26plotResults=False
 27
 28tEnd = 3
 29h= 1e-4
 30
 31SC = exu.SystemContainer()
 32mbs = SC.AddSystem()
 33
 34
 35gravity = 9.81
 36L=1.        #length of ANCF element in m    
 37rhoA=10     #beam + discrete masses
 38hBeam = 0.05
 39wBeam = 0.05
 40Abeam = hBeam*wBeam
 41Ibeam = wBeam*hBeam**3/12
 42Ebeam = 2.1e10
 43nu = 0.3
 44
 45EA=Abeam*Ebeam
 46EI=Ibeam*Ebeam
 47nElements = 25
 48lElem = L/nElements
 49
 50
 51# #additional bending and axial damping
 52bendingDamping=0*0.1*EI # for ALE Element
 53axialDamping=0 # for ALE Element
 54
 55#generate coordinate marker
 56#nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 57#mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 58
 59#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 60#create one beam template
 61cable = Cable2D(#length=L, 
 62                massPerLength=rhoA, 
 63                bendingStiffness=EI, 
 64                axialStiffness=EA, 
 65                bendingDamping=bendingDamping, 
 66                axialDamping=axialDamping, 
 67                # useCouplingTerms = True,
 68                useReducedOrderIntegration = 0, #faster
 69                visualization=VCable2D(drawHeight=hBeam)
 70                )
 71
 72#alternative to mbs.AddObject(ALECable2D(...)) with nodes:
 73yOff = 0*0.5*hBeam
 74ancf=GenerateBeamElementsAlongLine(mbs=mbs,
 75                positionStart=[0,yOff,0], positionEnd=[L,yOff,0],
 76                numberOfElements=nElements,
 77                beamTemplate=cable, #this defines the beam element properties
 78                gravity = [0,-gravity,0], #add larger gravity for larger deformation
 79                groundConstraintsStart = [1,1,0,0], #hinged
 80                #fixedConstraintsNode1 = [0,0,0,0]) #free
 81                )
 82
 83ancfNodes = ancf['nodes']
 84ancfObjects = ancf['elements']
 85
 86oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
 87                                visualization=VObjectGround(graphicsData=[graphics.CheckerBoard(size=2)])))
 88
 89#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 90#sensorFileName = 'solution/beamTip.txt'
 91sTipNode = mbs.AddSensor(SensorNode(nodeNumber=ancfNodes[-1], storeInternal=True,
 92                            outputVariableType=exu.OutputVariableType.Position))
 93sPos = mbs.AddSensor(SensorBody(bodyNumber=ancfObjects[-1], storeInternal=True, localPosition=[lElem,0,0.],
 94                                outputVariableType=exu.OutputVariableType.Position))
 95sVel = mbs.AddSensor(SensorBody(bodyNumber=ancfObjects[-1], storeInternal=True, localPosition=[lElem,0,0.],
 96                                outputVariableType=exu.OutputVariableType.Velocity))
 97
 98
 99mbs.Assemble()
100
101simulationSettings = exu.SimulationSettings() #takes currently set values or default values
102
103simulationSettings.parallel.numberOfThreads = 4 #4 is optimal for 25 elements
104
105simulationSettings.solution.file.write = False
106simulationSettings.solution.sensors.writePeriod = h*100
107simulationSettings.timeIntegration.verboseMode = 1
108simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
109simulationSettings.timeIntegration.endTime = tEnd
110
111simulationSettings.timeIntegration.newton.useModifiedNewton = True
112simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6
113simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-6
114simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.8
115simulationSettings.timeIntegration.adaptiveStep = True #disable adaptive step reduction
116
117simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
118simulationSettings.show.statistics = True
119SC.visualizationSettings.loads.show = False
120SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.StrainLocal
121#SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.CurvatureLocal
122#SC.visualizationSettings.bodies.beams.axialTiling = 500
123#SC.visualizationSettings.bodies.beams.crossSectionTiling = 8
124
125if useGraphics:
126    SC.renderer.Start()
127    SC.renderer.DoIdleTasks()
128
129success = mbs.SolveDynamic(simulationSettings, 
130                           exudyn.DynamicSolverType.TrapezoidalIndex2)
131
132if useGraphics:
133    SC.renderer.DoIdleTasks()
134    #SC.renderer.DoIdleTasks()
135    SC.renderer.Stop() #safely close rendering window!        
136
137
138#%%++++++++++++++++++
139    if True:
140        import matplotlib.pyplot as plt
141        
142        from exudyn.signalProcessing import FilterSensorOutput
143
144        mbs.PlotSensor(sensorNumbers=[sPos,sPos], components=[0,1], 
145                   title='ang vel', closeAll=True,
146                   markerStyles=['','x ','o '], lineStyles=['-','',''])