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()