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