ANCFcantileverTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  ANCF Cable2D cantilever test
  5#
  6# Model:    Cantilever beam with cable elements
  7#
  8# Author:   Johannes Gerstmayr
  9# Date:     2019-11-15
 10# Update:   2022-03-16: get to run static example again, compared to paper!
 11#
 12# 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.
 13#
 14# *clean example*
 15#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 16
 17## import exudyn and utilities
 18import exudyn as exu
 19from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 20from exudyn.beams import GenerateBeamElementsAlongLine
 21import exudyn.graphics as graphics
 22
 23## create container and main system to work with
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26
 27
 28## create graphics background
 29rect = [-0.5,-2,2.5,0.5] #xmin,ymin,xmax,ymax
 30background = {'type':'Line', 'color':[0.1,0.1,0.8,1], 'data':[rect[0],rect[1],0, rect[2],rect[1],0, rect[2],rect[3],0, rect[0],rect[3],0, rect[0],rect[1],0]} #background
 31oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background])))
 32
 33#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 34## define beam dimensions and tip load
 35L=2                    # length of ANCF element in m
 36E=2.07e11             # Young's modulus of ANCF element in N/m^2
 37rho=7800               # density of ANCF element in kg/m^3
 38b=0.1                  # width of rectangular ANCF element in m
 39h=0.1                  # height of rectangular ANCF element in m
 40A=b*h                  # cross sectional area of ANCF element in m^2
 41I=b*h**3/12            # second moment of area of ANCF element in m^4
 42f=3*E*I/L**2           # tip load applied to ANCF element in N
 43
 44print("load f="+str(f))
 45
 46#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 47## generate ANCFCable2D template containing beam parameters
 48cableTemplate = Cable2D(#length = L / nElements, #set in GenerateBeamElementsAlongLine(...)
 49                        massPerLength = rho*A,
 50                        bendingStiffness = E*I,
 51                        axialStiffness = E*A,
 52                        useReducedOrderIntegration = 0,
 53                        #nodeNumbers = [0, 0], #will be filled in GenerateBeamElementsAlongLine(...)
 54                        )
 55
 56## define nodal positions of beam (3D vectors, while cable element is only 2D)
 57positionOfNode0 = [0, 0, 0] # starting point of line
 58positionOfNode1 = [L, 0, 0] # end point of line
 59
 60## number of cable elements for discretization
 61numberOfElements = 64
 62
 63## use utility function to create set of straight cable elements between two positions with options for constraints at supports
 64#alternative to mbs.AddObject(Cable2D(...)) with nodes:
 65ancf=GenerateBeamElementsAlongLine(mbs,
 66                positionOfNode0, positionOfNode1,
 67                numberOfElements,
 68                beamTemplate=cableTemplate, #this defines the beam element properties
 69                gravity = [0,-9.81*0,0], #optionally add gravity
 70                groundConstraintsStart = [1,1,0,1], #add constraints for pos and rot (r'_y)
 71                groundConstraintsEnd = [0,0,0,0])
 72
 73## add load vector on last node in y-direction
 74mANCFLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=ancf['nodes'][-1])) #ancf['nodes'][-1] = last node
 75mbs.AddLoad(Force(markerNumber = mANCFLast, loadVector = [0, -f, 0])) #will be changed in load steps
 76
 77
 78#%%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 79## assemble system and create simulation settings
 80mbs.Assemble()
 81
 82simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 83
 84tEnd = 0.1
 85h = 1e-4
 86simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
 87simulationSettings.timeIntegration.endTime = tEnd
 88simulationSettings.solution.file.write = True
 89simulationSettings.solution.file.writePeriod = simulationSettings.timeIntegration.endTime/1000
 90simulationSettings.show.computationTime = False
 91simulationSettings.timeIntegration.verboseMode = 1
 92
 93simulationSettings.timeIntegration.newton.useModifiedNewton = True
 94
 95simulationSettings.show.statistics = True
 96simulationSettings.show.computationTime = True
 97
 98SC.visualizationSettings.nodes.defaultSize = 0.01
 99simulationSettings.solution.file.information = "ANCF cantilever beam"
100simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
101
102doDynamicSimulation = True #switch between static and dynamic simulation
103
104
105if doDynamicSimulation:
106    ## do dynamic simulation
107    SC.renderer.Start()
108    mbs.SolveDynamic(simulationSettings)
109    SC.renderer.DoIdleTasks()
110    SC.renderer.Stop() #safely close rendering window!
111    ## 
112else:
113    ## perform static simulation with manual load stepping
114    simulationSettings.staticSolver.verboseMode = 0
115    
116    simulationSettings.staticSolver.newton.relativeTolerance = 1e-8
117    simulationSettings.staticSolver.newton.absoluteTolerance = 1e-3 #1 for 256 elements; needs to be larger for larger number of load steps
118    #simulationSettings.staticSolver.numberOfLoadSteps = 1
119    
120    nLoadSteps = 1;
121    for loadSteps in range(nLoadSteps):
122        nLoad = 0
123        loadValue = f**((loadSteps+1)/nLoadSteps) #geometric increment of loads
124        print('load='+str(loadValue))
125        
126        mbs.SetLoadParameter(nLoad, 'loadVector', [0, -loadValue,0])
127        print('load vector=' + str(mbs.GetLoadParameter(nLoad, 'loadVector')) )
128    
129        mbs.SolveStatic(simulationSettings, updateInitialValues=True)
130    
131        sol = mbs.systemData.GetODE2Coordinates()
132        
133        n = len(sol)
134        print('nEL=',numberOfElements, ', tip displacement: x='+str(sol[n-4])+', y='+str(sol[n-3])) 
135        #MATLAB 1 element: x=0.3622447298905063, y=0.9941447587249748 = paper "on the correct ..."
136        #2022-03-16:
137        # nEL= 1 ,  tip displacement: x=-0.36224472989050654,y=-0.9941447587249747
138        # nEL= 2 ,  tip displacement: x=-0.4889263085609102, y=-1.1752228652637502
139        # nEL= 4 ,  tip displacement: x=-0.5074287154557922, y=-1.2055337025602493
140        # nEL= 8 ,  tip displacement: x=-0.5085092365729895, y=-1.207197756093103
141        # nEL= 16 , tip displacement: x=-0.5085365799149556, y=-1.207238895003594
142        # nEL= 32 , tip displacement: x=-0.508537277761696,  y=-1.2072398264650905
143        # nEL= 64 , tip displacement: x=-0.5085373030408489, y=-1.207239853404364
144        # nEL= 128, tip displacement: x=-0.5085373043168473, y=-1.2072398545511795
145        # nEL= 256, tip displacement: x=-0.5085373043916903, y=-1.207239854614031
146        
147        #with second SolveStatic:
148        #nEL= 256 , tip displacement: x=-0.5085373043209366, y=-1.2072398545457574
149        #converged:                   x=-0.508537304326,     y=-1.207239854550
150
151        #here (OLD):
152        #1:  x=-0.36224472989050543, y=-0.994144758724973
153        #2:  x=-0.4889263083414858, y=-1.1752228650551666
154        #4:  x=-0.5074287151188892, y=-1.2055337022335404
155        #8:  x=-0.5085092364970802, y=-1.2071977560198281
156        #64: x=-0.5085373029700947, y=-1.2072398533360738
157        #256:x=-0.5085373043209689, y=-1.2072398545457785
158        
159        
160    
161    
162        #sol = mbs.systemData.GetODE2Coordinates(exu.ConfigurationType.Initial)
163        #print('initial values='+str(sol))