finiteSegmentMethod.py

You can view and download this file on Github: finiteSegmentMethod.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Example of 2D finite segment method compared with ANCF cable elements
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2021-06-16
  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
 16from exudyn.beams import GenerateBeamElementsAlongLine
 17import exudyn.graphics as graphics
 18
 19SC = exu.SystemContainer()
 20mbs = SC.AddSystem()
 21
 22#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 23#beam as finite segment method
 24L = 2           #m
 25EI = 1          #Nm^2
 26nSegments = 8*8   #
 27rhoA = 1        #kg/m
 28mass = rhoA*L
 29massPerSegment = mass/nSegments
 30segmentLength = L/nSegments
 31a = 0.05        #width (for drawing)
 32g = 9.81        #gravity m/s^2
 33offY = 0.2*0      #position offset of ANCF cable
 34
 35#mode='Trap'
 36mode='GA'
 37
 38
 39inertiaSegment = 0*massPerSegment/(12*segmentLength**2) #inertia of segment needs to be zero to agree with Bernoulli-Euler beam
 40
 41#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 42background = graphics.CheckerBoard([0,0,0],[0,0,1], size=L*3)
 43oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], 
 44                                   visualization=VObjectGround(graphicsData= [background])))
 45
 46mPrevious = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition = [0,0,0]))
 47mRotPrevious = -1
 48oSegmentLast = -1 #store last segment for sensor
 49for i in range(nSegments):
 50    graphicsBeam = graphics.Brick([0,0,0],[segmentLength, a, a], graphics.color.red)
 51    nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[(0.5+i)*segmentLength,0,0]))
 52    oRigid = mbs.AddObject(RigidBody2D(mass=massPerSegment, 
 53                                       inertia=inertiaSegment,
 54                                       nodeNumber=nRigid,
 55                                       visualization=VObjectRigidBody2D(graphicsData= [graphicsBeam])))
 56    oSegmentLast = oRigid
 57    mRigidMass = mbs.AddMarker(MarkerBodyMass(bodyNumber=oRigid))
 58    mbs.AddLoad(LoadMassProportional(markerNumber=mRigidMass, loadVector=[0,-g,0]))
 59    
 60    mLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[-0.5*segmentLength,0.,0.])) 
 61    mRight= mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.5*segmentLength,0.,0.])) 
 62    
 63    
 64    oJoint = mbs.AddObject(RevoluteJoint2D(markerNumbers=[mPrevious, mLeft]))
 65    mRot = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nRigid,coordinate=2)) #rotation coordinate
 66    
 67    if mRotPrevious != -1:
 68        mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mRotPrevious,mRot],
 69                                              stiffness=EI/segmentLength, damping=0,
 70                                              visualization=VCoordinateSpringDamper(show=False)))
 71    mRotPrevious = mRot     #for next segment
 72    mPrevious = mRight      #for next segment
 73    
 74
 75#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 76#create ANCF beam as reference
 77useANCF = False
 78if useANCF:
 79    cable = ObjectANCFCable2D(nodeNumbers=[0,0],length=segmentLength, 
 80                              massPerLength=rhoA, bendingStiffness=EI,
 81                              axialStiffness=EI*1e4, useReducedOrderIntegration=True,
 82                              visualization=VCable2D(drawHeight = a, color=graphics.color.steelblue))
 83    
 84    ANCFcable = GenerateBeamElementsAlongLine(mbs, positionStart=[0,offY,0], positionEnd=[L,offY,0], 
 85                                    numberOfElements=nSegments, beamTemplate=cable,
 86                                    gravity=[0,-g,0], groundConstraintsStart=[1,1,0,0], 
 87                                    groundConstraintsEnd=[0,0,0,0])
 88    cableNodeList, cableObjectList, loadList, cableNodePositionList, cableCoordinateConstraintList = ANCFcable['nodes'], ANCFcable['elements'], ANCFcable['loads'], ANCFcable['nodePositions'], ANCFcable['constraintObjects']
 89    oTipCable = cableObjectList[-1] #last cable element
 90
 91#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 92#sensors
 93if useANCF:
 94    sTipCable     = mbs.AddSensor(SensorBody(bodyNumber=oTipCable, localPosition=[segmentLength, 0,0],
 95                                             fileName='solution/sensorTipCable'+mode+'.txt',
 96                                             outputVariableType=exu.OutputVariableType.Position))
 97
 98sTipSegment   = mbs.AddSensor(SensorBody(bodyNumber=oSegmentLast , localPosition=[0.5*segmentLength, 0,0],
 99                                         fileName='solution/sensorTipSegment'+mode+'.txt',
100                                         outputVariableType=exu.OutputVariableType.Position))
101
102mbs.Assemble()
103
104h = 1e-3 #step size
105tEnd = 4
106
107simulationSettings = exu.SimulationSettings() #takes currently set values or default values
108
109simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
110simulationSettings.timeIntegration.endTime = tEnd
111simulationSettings.timeIntegration.verboseMode = 1
112
113simulationSettings.timeIntegration.newton.useModifiedNewton = True
114simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
115simulationSettings.show.statistics = True
116#simulationSettings.linearSolver.solverType  = exu.LinearSolverType.EigenSparse
117
118#SC.visualizationSettings.nodes.defaultSize = 0.05
119
120simulationSettings.solution.file.information = "Finite segment method"
121
122SC.renderer.Start()
123
124if mode == "Trap":
125    mbs.SolveDynamic(simulationSettings, 
126                     solverType=exu.DynamicSolverType.TrapezoidalIndex2)
127else:
128    mbs.SolveDynamic(simulationSettings)
129    
130
131SC.renderer.DoIdleTasks()
132#SC.renderer.DoIdleTasks()
133SC.renderer.Stop() #safely close rendering window!
134
135
136if True and useANCF:
137    
138    mbs.PlotSensor(sensorNumbers=[sTipCable, sTipSegment], components=[1,1]) #plot y components