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