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