ANCFCableBeamDampingTest.py

You can view and download this file on Github: ANCFCableBeamDampingTest.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2#
  3# Details:  This is a test script that compares the axial and bending damping
  4#           of the ANCF beam element with the ANCF 2D cable element.
  5# Authors:  Johannes Gerstmayr and Sebastian Weyrer
  6# Date:     2026-03-23
  7#
  8#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  9
 10import exudyn as exu
 11import numpy as np
 12from exudyn.utilities import *
 13import exudyn.graphics as graphics
 14
 15useGraphics = True #without test
 16#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 17#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 18try: #only if called from test suite
 19    from modelUnitTests import exudynTestGlobals #for globally storing test results
 20    useGraphics = exudynTestGlobals.useGraphics
 21except:
 22    class ExudynTestGlobals:
 23        pass
 24    exudynTestGlobals = ExudynTestGlobals()
 25#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 26useGraphics = False
 27
 28exu.Print('EXUDYN says hello with version', exu.config.Version())
 29
 30#++++++++++++++++++++++++++++++++++++++++++++++++++
 31# set parameters and behavior of script
 32useSolutionViewer = False
 33recordImages = False
 34loadVectorList = [[0, -5, 0],
 35                  [-100, 0, 0]
 36                  ]
 37# cable parameters
 38l = 0.5                                             # [m] length of cable
 39massPerLength = 2                                   # [kg/m]
 40cableDiameter = 0.02                                # [m]
 41E = 1e8                                             # [N/m^2] E modulus of cable (in all directions!)
 42
 43nElementsANCF2D = 4                                 # number of elements the ANCF 2D cable is made out of
 44nElementsANCFBeam = 8                               # number of elements the ANCF beam is made out of
 45
 46if useGraphics: #more accurate results
 47    nElementsANCF2D *= 4
 48    nElementsANCFBeam *= 4
 49
 50beta = 0.02*0                                         # Rayleigh damping for bending
 51betaAxial = 0.0005*0                                  # Rayleigh damping for axial deformation
 52
 53# automatically dependent parameters
 54lElement = l/nElementsANCFBeam                      # [m] length of one element (used in ANCF beam)
 55m = massPerLength*l                                 # [kg] full cable mass
 56cableRadius = cableDiameter/2                       # [m] radius of the cable
 57cableArea = np.pi*cableRadius**2                    # [m^2] cross section area of the cable
 58secondMomentOfInertia = np.pi*cableRadius**4/4      # [m^4] second moment of inertia (bending)
 59rho = massPerLength/cableArea                       # [kg/m^3] density of cable
 60rhoA = rho*cableArea                                # [kg/m] mass per unit length
 61rhoI = rho*secondMomentOfInertia                    # [kg*m] inertia per unit length
 62EA = E*cableArea                                    # [N] axial stiffness
 63EI = E*secondMomentOfInertia                        # [N*m^2] bending stiffness
 64
 65# set up cross section data resulting from parameters (assemble matrices where needed) (X --> torsion, Y and Z --> bending)
 66sectionData = exu.BeamSection()
 67kPenalty = EA * 1 # penalty stiffness that is added to eliminate shear deformation
 68sectionData.stiffnessMatrix = np.diag([EA, kPenalty, kPenalty, kPenalty, EI, EI])
 69sectionData.inertia = np.diag([0, rhoI, rhoI])
 70sectionData.massPerLength = rhoA
 71sectionData.dampingMatrix = np.diag([EA*betaAxial, 0, 0, 0, EI*beta, EI*beta])
 72
 73# visualization of cable
 74nTiles = 18
 75sectionGeometry = exu.BeamSectionGeometry()
 76lp = exu.Vector2DList()
 77phi = 2*np.pi/nTiles
 78for i in range(nTiles):
 79    lp.Append([cableRadius*np.cos(i*phi), cableRadius*np.sin(i*phi)])
 80sectionGeometry.polygonalPoints = lp
 81
 82#%% set up system container (once for all tests)
 83SC = exu.SystemContainer()
 84
 85# general visualization settings
 86SC.visualizationSettings.nodes.show = True
 87SC.visualizationSettings.markers.show = True
 88SC.visualizationSettings.markers.drawSimplified = True
 89SC.visualizationSettings.view0.scene.drawCoordinateSystem = True
 90SC.visualizationSettings.bodies.beams.axialTiling = 1
 91# simulation settings
 92simulationSettings = exu.SimulationSettings()
 93simulationSettings.timeIntegration.newton.useModifiedNewton = True
 94simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = False
 95simulationSettings.displayComputationTime = True
 96simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
 97
 98if recordImages:
 99    simulationSettings.solutionSettings.recordImagesInterval = 0.1
100mbs = SC.AddSystem()
101
102solution = 0 #accumulated result for test suite
103
104# iterate over the load vectors
105for loadCase, loadVector in enumerate(loadVectorList):
106    #%% set up general things for the test
107    mbs.Reset() # reset mbs since we now make another test
108    oGround = mbs.CreateGround(referencePosition=[0, 0, 0])
109
110    #%% set up the ANCF 2D cable
111    cableTemplate = Cable2D(physicsMassPerLength=rhoA,
112                            physicsBendingStiffness=EI,
113                            physicsAxialStiffness=EA,
114                            physicsBendingDamping=beta*EI,
115                            physicsAxialDamping=betaAxial*EA,
116                            useReducedOrderIntegration=0,
117                            visualization=VCable2D(drawHeight=cableDiameter))
118    nCable2D, oCable2D, lCable2D, _, _ = GenerateStraightLineANCFCable2D(mbs, positionOfNode0=[0, 0, 0],
119                                                                         positionOfNode1=[l, 0, 0],
120                                                                         fixedConstraintsNode0=[1]*4,
121                                                                         numberOfElements=nElementsANCF2D,
122                                                                         cableTemplate=cableTemplate)
123    mCable2D = mbs.AddMarker(MarkerNodePosition(nodeNumber=nCable2D[-1]))
124    # add load to last node
125    mbs.AddLoad(LoadForceVector(markerNumber=mCable2D,
126                                loadVector=loadVector, bodyFixed=False))
127    # add sensor to get position of last node
128    sCable2D = mbs.AddSensor(SensorMarker(markerNumber=mCable2D,
129                                          outputVariableType=exu.OutputVariableType.Displacement,
130                                          storeInternal=True))
131
132    #%% set up ANCF beam cable
133    referenceOffset = [0, 0, 1]
134    initialRotations = [0, 1, 0] + [0, 0, 1]
135    mCableList = [] # this list holds the markers to which the discs can then be attached
136    n0 = mbs.AddNode(NodePointSlope23(referenceCoordinates=referenceOffset + initialRotations))
137    mCableList += [mbs.AddMarker(MarkerNodeRigid(nodeNumber=n0))]
138    for k in range(nElementsANCFBeam):
139        n1 = mbs.AddNode(NodePointSlope23(referenceCoordinates=[lElement*(k+1) + referenceOffset[0], referenceOffset[1], referenceOffset[2]] + initialRotations, visualization=VNodePointSlope23(show=True)))
140        mCableList += [mbs.AddMarker(MarkerNodeRigid(nodeNumber=n1))]
141        oBeam = mbs.AddObject(ObjectANCFBeam(nodeNumbers=[n0, n1],
142                                             physicsLength=lElement,
143                                             sectionData=sectionData, #includes bending stiffness, axial stiffness, damping, etc.
144                                             crossSectionPenaltyFactor=[1]*3,
145                                             visualization=VANCFBeam(sectionGeometry=sectionGeometry,
146                                                                     color=graphics.color.grey)))
147        n0 = n1
148
149    # fix cable to ground
150    mbs.CreateGenericJoint(bodyNumbers=[oGround, mCableList[0]], show=True)
151    # add load to last node
152    mbs.AddLoad(LoadForceVector(markerNumber=mCableList[-1],
153                                loadVector=loadVector, bodyFixed=False))
154    # add sensor to get position of last node
155    sBeam = mbs.AddSensor(SensorMarker(markerNumber=mCableList[-1],
156                                       outputVariableType=exu.OutputVariableType.Displacement,
157                                       storeInternal=True))
158
159    #%% do simulation (implicit dynamic simulation) (WYSWYS - What You See is What You Simulate)
160    mbs.Assemble()
161
162
163    stepSize = 4e-3
164    tEnd = 4
165    if loadCase == 1:
166        stepSize = 0.5e-3
167        tEnd = 0.25
168
169    if not useGraphics: #for test suite
170        tEnd = 0.1
171
172    #tEnd = 0.1
173
174    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
175    simulationSettings.timeIntegration.endTime = tEnd
176    simulationSettings.solutionSettings.writeSolutionToFile = useSolutionViewer
177    simulationSettings.solutionSettings.sensorsWritePeriod = stepSize
178
179    #++++++++++++++++++++++++++++++++++++++++++++++++++
180    if useGraphics:
181        SC.renderer.Start()
182        SC.renderer.DoIdleTasks()
183    mbs.SolveDynamic(simulationSettings=simulationSettings)
184    if useGraphics:
185        SC.renderer.DoIdleTasks()
186        SC.renderer.Stop() #safely close rendering window!
187
188    #++++++++++++++++++++++++++++++++++++++++++++++++++
189    if useSolutionViewer and useGraphics:
190        mbs.SolutionViewer()
191
192    #%% plot displacement to compare results
193    if useGraphics:
194        [plt, fig, ax, line] = mbs.PlotSensor(sensorNumbers=[sCable2D, sBeam],
195                       components=[1-loadCase]*2,
196                       title="Test Damping: load="+str(loadVector),
197                       labels=['ANCF cable 2D','ANCF beam'], fontSize=12,
198                       legendArgs=(0.65,0.78),
199                       )
200
201    #for test suite:
202    solution += np.linalg.norm(mbs.GetSensorValues(sCable2D))
203    solution += np.linalg.norm(mbs.GetSensorValues(sBeam))
204
205exu.Print('ANCFCableBeamDampingTest: solution=', solution)
206
207exudynTestGlobals.testResult = solution