bungeeJump.py
You can view and download this file on Github: bungeeJump.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Example to simulate a bungee jumper
5#
6# Author: Johannes Gerstmayr
7# Date: 2024-04-21
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.itemInterface import *
15from exudyn.utilities import *
16from exudyn.beams import *
17
18import numpy as np
19from math import sin, cos, pi, sqrt , asin, acos, atan2
20import copy
21
22SC = exu.SystemContainer()
23mbs = SC.AddSystem()
24
25
26
27#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
28#settings:
29useGraphics= True
30tEnd = 30 #end time of dynamic simulation
31stepSize = 2e-3 #step size
32
33
34#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
35#create circles
36#complicated shape:
37nANCFnodes = 80
38h = 0.25e-3
39preStretch=0
40circleList = [[[-2,0],1,'L'],
41 [[-0.5,-40],0.5,'R'],
42 [[1,0],1,'L'],
43 ]
44
45
46#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
47#create geometry:
48reevingDict = CreateReevingCurve(circleList, drawingLinesPerCircle = 32,
49 removeLastLine = False, #True allows closed curve
50 closedCurve = False,
51 numberOfANCFnodes=nANCFnodes, graphicsNodeSize= 0.01)
52
53
54gList=[]
55if False: #visualize reeving curve, without simulation
56 gList = reevingDict['graphicsDataLines'] + reevingDict['graphicsDataCircles']
57
58#bridge graphics
59bH = 4
60bW = 30
61bT = 80
62tH = 190
63tW = 8
64gList += [GraphicsDataOrthoCubePoint([-bW*0.5,-bH*0.5,0],size=[bW-2,bH,bT],color=color4grey)]
65gList += [GraphicsDataOrthoCubePoint([-1,-0.1,0],size=[2,0.2,2],color=color4grey)]
66gList += [GraphicsDataOrthoCubePoint([-bW*0.5,-tH*0.5,0],size=[tW,tH,tW],color=[0.6,0.6,0.6,0.2])]
67gList += [GraphicsDataOrthoCubePoint([-25+10,-tH-2,0],size=[50,4,bT],color=color4steelblue)]
68
69
70#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
71#create ANCF elements:
72#values are adjusted to have reasonable stiffness and damping!!! no measured values!
73gVec = [0,-9.81,0] # gravity
74E=1e7 # Young's modulus of ANCF element in N/m^2
75rhoBeam=1000 # density of ANCF element in kg/m^3
76b=0.020 # width of rectangular ANCF element in m
77h=0.020 # 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 = 200e-3*E*I #bending proportional damping
81dEA = 100e-3*E*A #axial strain proportional damping
82
83# dimZ = b #z.dimension
84
85cableTemplate = Cable2D(#physicsLength = L / nElements, #set in GenerateStraightLineANCFCable2D(...)
86 physicsMassPerLength = rhoBeam*A,
87 physicsBendingStiffness = E*I*10, #increase bending stiffness to avoid buckling and numerical issues
88 physicsAxialStiffness = E*A,
89 physicsBendingDamping = dEI,
90 physicsAxialDamping = dEA,
91 physicsReferenceAxialStrain = preStretch, #prestretch
92 visualization=VCable2D(drawHeight=0.05),
93 )
94
95ancf = PointsAndSlopes2ANCFCable2D(mbs, reevingDict['ancfPointsSlopes'], reevingDict['elementLengths'],
96 cableTemplate, massProportionalLoad=gVec,
97 fixedConstraintsNode0=[1,1,1,1], #fixedConstraintsNode1=[1,1,1,1],
98 firstNodeIsLastNode=False, graphicsSizeConstraints=0.01)
99
100nLast = ancf[0][-1]
101mCable = mbs.AddMarker(MarkerNodePosition(nodeNumber=nLast))
102
103#add jumper as rigid body
104gJumper = []
105hJumper = 1.8
106#gJumper += [GraphicsDataOrthoCubePoint([0,0,0],size=[0.4,1.8,0.5],color=color4blue)]
107gJumper += [GraphicsDataOrthoCubePoint([0,0.3,0],size=[0.5,0.64,0.5],color=color4blue)]
108gJumper += [GraphicsDataOrthoCubePoint([0,-0.25*hJumper,0],size=[0.3,0.5*hJumper,0.5],color=color4darkgrey)]
109gJumper += [GraphicsDataSphere([0,0.75,0],radius=0.15,color=color4orange)]
110
111bJumper = mbs.CreateRigidBody(referencePosition=[0,0.5*hJumper,0],
112 inertia=InertiaCuboid(250, sideLengths=[0.4,1.8,0.5]), #90kg
113 initialVelocity=[0.25,0,0],
114 initialAngularVelocity=[0,0,-pi*0.2],
115 gravity=gVec,
116 graphicsDataList=gJumper)
117
118bGround = mbs.CreateGround()
119fixJumper = mbs.CreateGenericJoint(bodyNumbers=[bJumper, bGround],position=[0,0,0],useGlobalFrame=False,
120 constrainedAxes=[1,1,0, 1,1,1])
121
122mJumper = mbs.AddMarker(MarkerBodyPosition(bodyNumber=bJumper, localPosition=[0,-0.5*hJumper,0]))
123mbs.AddObject(SphericalJoint(markerNumbers=[mCable, mJumper]))
124
125sPosJumper = mbs.AddSensor(SensorBody(bodyNumber=bJumper, storeInternal=True,
126 outputVariableType=exu.OutputVariableType.Position))
127sVelJumper = mbs.AddSensor(SensorBody(bodyNumber=bJumper, storeInternal=True,
128 outputVariableType=exu.OutputVariableType.Velocity))
129sAccJumper = mbs.AddSensor(SensorBody(bodyNumber=bJumper, storeInternal=True,
130 outputVariableType=exu.OutputVariableType.Acceleration))
131
132#transparent
133oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
134 visualization=VObjectGround(graphicsData= gList)))
135
136mbs.Assemble()
137
138simulationSettings = exu.SimulationSettings() #takes currently set values or default values
139
140simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
141simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/coordinatesSolution.txt'
142simulationSettings.solutionSettings.writeSolutionToFile = True
143# simulationSettings.displayComputationTime = True
144simulationSettings.parallel.numberOfThreads = 1 #use 4 to speed up for > 100 ANCF elements
145# simulationSettings.displayStatistics = True
146
147simulationSettings.timeIntegration.endTime = tEnd
148simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
149# simulationSettings.timeIntegration.stepInformation= 3+128+256
150simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = True
151simulationSettings.timeIntegration.newton.useModifiedNewton = True
152simulationSettings.timeIntegration.verboseMode = 1
153
154
155SC.visualizationSettings.general.circleTiling = 24
156SC.visualizationSettings.loads.show=False
157SC.visualizationSettings.nodes.defaultSize = 0.01
158SC.visualizationSettings.openGL.multiSampling = 4
159
160
161if useGraphics:
162 exu.StartRenderer()
163 # mbs.WaitForUserToContinue()
164
165simulationSettings.staticSolver.numberOfLoadSteps = 10
166simulationSettings.staticSolver.stabilizerODE2term = 1
167#compute initial static solution
168mbs.SolveStatic(simulationSettings, updateInitialValues=False)
169ode2 = mbs.systemData.GetODE2Coordinates()
170mbs.systemData.SetODE2Coordinates(ode2, configuration=exu.ConfigurationType.Initial)
171
172#turn of constraint of jumper
173mbs.SetObjectParameter(fixJumper[0], parameterName='activeConnector', value=False)
174
175mbs.WaitForUserToContinue()
176
177mbs.SolveDynamic(simulationSettings) #183 Newton iterations, 0.114 seconds
178
179
180# mbs.SolutionViewer()
181
182
183if useGraphics:
184 SC.WaitForRenderEngineStopFlag()
185 exu.StopRenderer() #safely close rendering window!
186
187#%%
188if True:
189 mbs.PlotSensor([sPosJumper],components=[1],closeAll=True)
190 mbs.PlotSensor([sVelJumper],components=[1])
191 mbs.PlotSensor([sAccJumper],components=[1])