finiteSegmentMethod.py
You can view and download this file on Github: finiteSegmentMethod.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Example of 2D finite segment method compared with ANCF cable elements
5#
6# Author: Johannes Gerstmayr
7# Date: 2021-06-16
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.itemInterface import *
15from exudyn.utilities import *
16
17SC = exu.SystemContainer()
18mbs = SC.AddSystem()
19
20#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
21#beam as finite segment method
22L = 2 #m
23EI = 1 #Nm^2
24nSegments = 8*8 #
25rhoA = 1 #kg/m
26mass = rhoA*L
27massPerSegment = mass/nSegments
28segmentLength = L/nSegments
29a = 0.05 #width (for drawing)
30g = 9.81 #gravity m/s^2
31offY = 0.2*0 #position offset of ANCF cable
32
33#mode='Trap'
34mode='GA'
35
36
37inertiaSegment = 0*massPerSegment/(12*segmentLength**2) #inertia of segment needs to be zero to agree with Bernoulli-Euler beam
38
39#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
40background = GraphicsDataCheckerBoard([0,0,0],[0,0,1], size=L*3)
41oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
42 visualization=VObjectGround(graphicsData= [background])))
43
44mPrevious = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oGround, localPosition = [0,0,0]))
45mRotPrevious = -1
46oSegmentLast = -1 #store last segment for sensor
47for i in range(nSegments):
48 graphicsBeam = GraphicsDataOrthoCubePoint([0,0,0],[segmentLength, a, a], color4red)
49 nRigid = mbs.AddNode(Rigid2D(referenceCoordinates=[(0.5+i)*segmentLength,0,0]))
50 oRigid = mbs.AddObject(RigidBody2D(physicsMass=massPerSegment,
51 physicsInertia=inertiaSegment,
52 nodeNumber=nRigid,
53 visualization=VObjectRigidBody2D(graphicsData= [graphicsBeam])))
54 oSegmentLast = oRigid
55 mRigidMass = mbs.AddMarker(MarkerBodyMass(bodyNumber=oRigid))
56 mbs.AddLoad(LoadMassProportional(markerNumber=mRigidMass, loadVector=[0,-g,0]))
57
58 mLeft = mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[-0.5*segmentLength,0.,0.]))
59 mRight= mbs.AddMarker(MarkerBodyPosition(bodyNumber=oRigid, localPosition=[ 0.5*segmentLength,0.,0.]))
60
61
62 oJoint = mbs.AddObject(RevoluteJoint2D(markerNumbers=[mPrevious, mLeft]))
63 mRot = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nRigid,coordinate=2)) #rotation coordinate
64
65 if mRotPrevious != -1:
66 mbs.AddObject(CoordinateSpringDamper(markerNumbers=[mRotPrevious,mRot],
67 stiffness=EI/segmentLength, damping=0,
68 visualization=VCoordinateSpringDamper(show=False)))
69 mRotPrevious = mRot #for next segment
70 mPrevious = mRight #for next segment
71
72
73#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
74#create ANCF beam as reference
75useANCF = False
76if useANCF:
77 cable = ObjectANCFCable2D(nodeNumbers=[0,0],physicsLength=segmentLength,
78 physicsMassPerLength=rhoA, physicsBendingStiffness=EI,
79 physicsAxialStiffness=EI*1e4, useReducedOrderIntegration=True,
80 visualization=VCable2D(drawHeight = a, color=color4steelblue))
81
82 ANCFcable = GenerateStraightLineANCFCable2D(mbs, positionOfNode0=[0,offY,0], positionOfNode1=[L,offY,0],
83 numberOfElements=nSegments, cableTemplate=cable,
84 massProportionalLoad=[0,-g,0], fixedConstraintsNode0=[1,1,0,0],
85 fixedConstraintsNode1=[0,0,0,0])
86 [cableNodeList, cableObjectList, loadList, cableNodePositionList, cableCoordinateConstraintList] = ANCFcable
87 oTipCable = cableObjectList[-1] #last cable element
88
89#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
90#sensors
91if useANCF:
92 sTipCable = mbs.AddSensor(SensorBody(bodyNumber=oTipCable, localPosition=[segmentLength, 0,0],
93 fileName='solution/sensorTipCable'+mode+'.txt',
94 outputVariableType=exu.OutputVariableType.Position))
95
96sTipSegment = mbs.AddSensor(SensorBody(bodyNumber=oSegmentLast , localPosition=[0.5*segmentLength, 0,0],
97 fileName='solution/sensorTipSegment'+mode+'.txt',
98 outputVariableType=exu.OutputVariableType.Position))
99
100mbs.Assemble()
101
102h = 1e-3 #step size
103tEnd = 4
104
105simulationSettings = exu.SimulationSettings() #takes currently set values or default values
106
107simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
108simulationSettings.timeIntegration.endTime = tEnd
109simulationSettings.timeIntegration.verboseMode = 1
110
111simulationSettings.timeIntegration.newton.useModifiedNewton = True
112simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
113simulationSettings.displayStatistics = True
114#simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
115
116#SC.visualizationSettings.nodes.defaultSize = 0.05
117
118simulationSettings.solutionSettings.solutionInformation = "Finite segment method"
119
120exu.StartRenderer()
121
122if mode == "Trap":
123 mbs.SolveDynamic(simulationSettings,
124 solverType=exu.DynamicSolverType.TrapezoidalIndex2)
125else:
126 mbs.SolveDynamic(simulationSettings)
127
128
129SC.WaitForRenderEngineStopFlag()
130#SC.WaitForRenderEngineStopFlag()
131exu.StopRenderer() #safely close rendering window!
132
133
134if True and useANCF:
135
136 mbs.PlotSensor(sensorNumbers=[sTipCable, sTipSegment], components=[1,1]) #plot y components