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