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