finiteSegmentMethod.rst

finiteSegmentMethod.py

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

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details:  Example of 2D finite segment method compared with ANCF cable elements
#
# Author:   Johannes Gerstmayr
# Date:     2021-06-16
#
# 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.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

import exudyn as exu
from exudyn.itemInterface import *
from exudyn.utilities import *

SC = exu.SystemContainer()
mbs = SC.AddSystem()

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#beam as finite segment method
L = 2           #m
EI = 1          #Nm^2
nSegments = 8*8   #
rhoA = 1        #kg/m
mass = rhoA*L
massPerSegment = mass/nSegments
segmentLength = L/nSegments
a = 0.05        #width (for drawing)
g = 9.81        #gravity m/s^2
offY = 0.2*0      #position offset of ANCF cable

#mode='Trap'
mode='GA'


inertiaSegment = 0*massPerSegment/(12*segmentLength**2) #inertia of segment needs to be zero to agree with Bernoulli-Euler beam

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
background = GraphicsDataCheckerBoard([0,0,0],[0,0,1], size=L*3)
oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
                                   visualization=VObjectGround(graphicsData= [background])))

mPrevious = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition = [0,0,0]))
mRotPrevious = -1
oSegmentLast = -1 #store last segment for sensor
for i in range(nSegments):
    graphicsBeam = GraphicsDataOrthoCubePoint([0,0,0],[segmentLength, a, a], color4red)
    nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[(0.5+i)*segmentLength,0,0]))
    oRigid = mbs.AddObject(RigidBody2D(physicsMass=massPerSegment,
                                       physicsInertia=inertiaSegment,
                                       nodeNumber=nRigid,
                                       visualization=VObjectRigidBody2D(graphicsData= [graphicsBeam])))
    oSegmentLast = oRigid
    mRigidMass = mbs.AddMarker(MarkerBodyMass(bodyNumber=oRigid))
    mbs.AddLoad(LoadMassProportional(markerNumber=mRigidMass, loadVector=[0,-g,0]))

    mLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[-0.5*segmentLength,0.,0.]))
    mRight= mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.5*segmentLength,0.,0.]))


    oJoint = mbs.AddObject(RevoluteJoint2D(markerNumbers=[mPrevious, mLeft]))
    mRot = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nRigid,coordinate=2)) #rotation coordinate

    if mRotPrevious != -1:
        mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mRotPrevious,mRot],
                                              stiffness=EI/segmentLength, damping=0,
                                              visualization=VCoordinateSpringDamper(show=False)))
    mRotPrevious = mRot     #for next segment
    mPrevious = mRight      #for next segment


#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#create ANCF beam as reference
useANCF = False
if useANCF:
    cable = ObjectANCFCable2D(nodeNumbers=[0,0],physicsLength=segmentLength,
                              physicsMassPerLength=rhoA, physicsBendingStiffness=EI,
                              physicsAxialStiffness=EI*1e4, useReducedOrderIntegration=True,
                              visualization=VCable2D(drawHeight = a, color=color4steelblue))

    ANCFcable = GenerateStraightLineANCFCable2D(mbs, positionOfNode0=[0,offY,0], positionOfNode1=[L,offY,0],
                                    numberOfElements=nSegments, cableTemplate=cable,
                                    massProportionalLoad=[0,-g,0], fixedConstraintsNode0=[1,1,0,0],
                                    fixedConstraintsNode1=[0,0,0,0])
    [cableNodeList, cableObjectList, loadList, cableNodePositionList, cableCoordinateConstraintList] = ANCFcable
    oTipCable = cableObjectList[-1] #last cable element

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#sensors
if useANCF:
    sTipCable     = mbs.AddSensor(SensorBody(bodyNumber=oTipCable, localPosition=[segmentLength, 0,0],
                                             fileName='solution/sensorTipCable'+mode+'.txt',
                                             outputVariableType=exu.OutputVariableType.Position))

sTipSegment   = mbs.AddSensor(SensorBody(bodyNumber=oSegmentLast , localPosition=[0.5*segmentLength, 0,0],
                                         fileName='solution/sensorTipSegment'+mode+'.txt',
                                         outputVariableType=exu.OutputVariableType.Position))

mbs.Assemble()

h = 1e-3 #step size
tEnd = 4

simulationSettings = exu.SimulationSettings() #takes currently set values or default values

simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.timeIntegration.verboseMode = 1

simulationSettings.timeIntegration.newton.useModifiedNewton = True
simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
simulationSettings.displayStatistics = True
#simulationSettings.linearSolverType  = exu.LinearSolverType.EigenSparse

#SC.visualizationSettings.nodes.defaultSize = 0.05

simulationSettings.solutionSettings.solutionInformation = "Finite segment method"

exu.StartRenderer()

if mode == "Trap":
    mbs.SolveDynamic(simulationSettings,
                     solverType=exu.DynamicSolverType.TrapezoidalIndex2)
else:
    mbs.SolveDynamic(simulationSettings)


SC.WaitForRenderEngineStopFlag()
#SC.WaitForRenderEngineStopFlag()
exu.StopRenderer() #safely close rendering window!


if True and useANCF:

    mbs.PlotSensor(sensorNumbers=[sTipCable, sTipSegment], components=[1,1]) #plot y components