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