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