ANCFrotatingCable2D.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  ANCF Cable2D cantilever test
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2023-11-07
  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.utilities import * #includes itemInterface and rigidBodyUtilities
 15from exudyn.beams import GenerateBeamElementsAlongLine
 16from math import sin
 17import exudyn.graphics as graphics
 18
 19SC = exu.SystemContainer()
 20mbs = SC.AddSystem()
 21
 22
 23#background
 24background = graphics.CheckerBoard(point=[0,0,-0.1],size = 5)
 25oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
 26
 27#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 28#cable:
 29
 30L=2                    # length of ANCF element in m
 31E=2e11                  # Young's modulus of ANCF element in N/m^2
 32rho=7800               # density of ANCF element in kg/m^3
 33b=0.01                 # width of rectangular ANCF element in m
 34h=0.01                 # height of rectangular ANCF element in m
 35A=b*h                  # cross sectional area of ANCF element in m^2
 36I=b*h**3/12            # second moment of area of ANCF element in m^4
 37f=3*E*I/L**2           # tip load applied to ANCF element in N
 38
 39print("load f="+str(f))
 40
 41#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 42#generate ANCF beams with utilities function
 43cableTemplate = Cable2D(#length = L / nElements, #set in GenerateBeamElementsAlongLine(...)
 44                        massPerLength = rho*A,
 45                        bendingStiffness = E*I,
 46                        axialStiffness = E*A,
 47                        bendingDamping = 0.02*E*I,
 48                        useReducedOrderIntegration = 0,
 49                        visualization=VCable2D(drawHeight=h),
 50                        #nodeNumbers = [0, 0], #will be filled in GenerateBeamElementsAlongLine(...)
 51                        )
 52
 53positionOfNode0 = [0, 0, 0] # starting point of line
 54positionOfNode1 = [L, 0, 0] # end point of line
 55numberOfElements = 16
 56
 57#alternative to mbs.AddObject(Cable2D(...)) with nodes:
 58ancf=GenerateBeamElementsAlongLine(mbs,
 59                positionOfNode0, positionOfNode1,
 60                numberOfElements,
 61                beamTemplate=cableTemplate, #this defines the beam element properties
 62                gravity = [0,-9.81,0], #optionally add gravity
 63                #fixedConstraintsNode0 = [1,1,0,1], #add constraints for pos and rot (r'_y)
 64                #fixedConstraintsNode1 = [0,0,0,0]
 65                )
 66
 67ancfNodes = ancf['nodes']
 68# #force applied to last node:
 69# mANCFLast = mbs.AddMarker(MarkerNodeRigid(nodeNumber=ancfNodes[1])) #ancf['nodes'][-1] = last node
 70# mbs.AddLoad(Force(markerNumber = mANCFLast, loadVector = [0, -f, 0])) #will be changed in load steps
 71
 72#torque and clamping of first node:
 73mANCFFirst = mbs.AddMarker(MarkerNodeRigid(nodeNumber=ancfNodes[0])) #ancf['nodes'][-1] = last node
 74
 75if True:
 76    #create rigid body:
 77    gBody = graphics.Brick(size = [h,h,h], color=graphics.color.red)
 78    dictBody = mbs.CreateRigidBody(referencePosition=[0,0,0],
 79                                inertia = InertiaCuboid(1000, [h,h,h]),
 80                                graphicsDataList=[gBody],
 81                                create2D = True, returnDict=True)
 82    
 83    #connect rigid body with ANCF
 84    mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=dictBody['bodyNumber'], localPosition=[0,0,0]))
 85    mbs.AddObject(GenericJoint(markerNumbers=[mANCFFirst,mBody], constrainedAxes=[1,1,0, 0,0,1],
 86                               visualization=VGenericJoint(axesRadius=h*0.5,axesLength=h)))
 87
 88    #connect rigid body with ground
 89    mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround,localPosition=[0,0,0]))
 90    mbs.AddObject(RevoluteJoint2D(markerNumbers=[mBody,mGround],
 91                               visualization=VRevoluteJoint2D(drawSize=h*0.5)))
 92
 93    #prescribe rotation of rigid body    
 94    nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 95    mcGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 96    mBodyPhi = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber= dictBody['nodeNumber'], coordinate = 2))
 97    
 98    def UFoffset(mbs, t, itemNumber, lOffset):
 99        if t<2:
100            return 0.
101        elif t<6:
102            return pi*sin(pi*t)
103        else:
104            return 0.
105            
106
107    mbs.AddObject(CoordinateConstraint(markerNumbers = [mcGround, mBodyPhi],
108                                   offset = 0.,
109                                   offsetUserFunction = UFoffset))
110
111else:
112    #possibility to fix to ground:
113    mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround,localPosition=[0,0,0]))
114    mbs.AddObject(GenericJoint(markerNumbers=[mANCFLast,mGround], constrainedAxes=[1,1,0, 0,0,1],
115                               visualization=VGenericJoint(axesRadius=h*0.5,axesLength=h)))
116    
117#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
118mbs.Assemble()
119# print(mbs)
120simulationSettings = exu.SimulationSettings() #takes currently set values or default values
121
122tEnd = 10
123h = 2e-3
124simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
125simulationSettings.timeIntegration.endTime = tEnd
126simulationSettings.solution.file.write = True
127simulationSettings.solution.file.writePeriod = simulationSettings.timeIntegration.endTime/1000
128simulationSettings.show.computationTime = False
129simulationSettings.timeIntegration.verboseMode = 1
130
131simulationSettings.timeIntegration.newton.useModifiedNewton = True
132simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6
133#simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
134#simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
135
136
137SC.visualizationSettings.nodes.defaultSize = 0.01
138
139simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
140
141mbs.SolveDynamic(simulationSettings)
142
143mbs.SolutionViewer()