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