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])