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=['-','',''])