ANCFtests2.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  ANCF Cable2D further test file
  5#           Example shows limitations of static solver: larger bending not possible;
  6#           larger number of elements (>16) leads to convergence problems
  7#
  8# Author:   Johannes Gerstmayr
  9# Date:     2019-11-15
 10#
 11# 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.
 12#
 13#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 14
 15import exudyn as exu
 16from exudyn.itemInterface import *
 17
 18SC = exu.SystemContainer()
 19mbs = SC.AddSystem()
 20
 21
 22#background
 23rect = [-2,-2,2,2] #xmin,ymin,xmax,ymax
 24background0 = {'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
 25oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background0])))
 26
 27#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 28#cable:
 29mypi = 3.141592653589793
 30
 31L=2                   # length of ANCF element in m
 32#L=mypi                 # length of ANCF element in m
 33E=2.07e11              # Young's modulus of ANCF element in N/m^2
 34rho=7800*10               # density of ANCF element in kg/m^3
 35b=0.1                  # width of rectangular ANCF element in m
 36h=0.1                  # height of rectangular ANCF element in m
 37A=b*h                  # cross sectional area of ANCF element in m^2
 38I=b*h**3/12            # second moment of area of ANCF element in m^4
 39f=3*E*I/L**2           # tip load applied to ANCF element in N
 40
 41print("load f="+str(f))
 42print("EI="+str(E*I))
 43
 44nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 45mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 46
 47cableList=[]
 48
 49mode = 1
 50if mode==0: #treat one element
 51    #omega = mypi*2
 52    #nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0],initialVelocities=[0,-L/2*omega,0,omega])) #initial velocity
 53    #nc1 = mbs.AddNode(Point2DS1(referenceCoordinates=[L,0,1,0],initialVelocities=[0, L/2*omega,0,omega])) #initial velocity
 54    nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 55    nc1 = mbs.AddNode(Point2DS1(referenceCoordinates=[L,0,1,0]))
 56
 57    mbs.systemData.Info()
 58    o0 = mbs.AddObject(Cable2D(name='FirstCable', physicsLength=L, physicsMassPerLength=rho*A, physicsBendingStiffness=E*I, physicsAxialStiffness=E*A, nodeNumbers=[nc0,nc1]))
 59    cableList+=[o0]
 60
 61    myObject = mbs.GetObject('FirstCable')
 62    print(myObject)
 63    #print(mbs.GetObject(o0))
 64
 65    mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))
 66    mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))
 67    mANCF2b = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3))
 68
 69    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
 70    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
 71    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2b]))
 72
 73    #mANCFnode = mbs.AddMarker(MarkerNodePosition(nodeNumber=nc1)) #force
 74    #mbs.AddLoad(Force(markerNumber = mANCFnode, loadVector = [0, -10000, 0]))
 75    mANCFrigid = mbs.AddMarker(MarkerBodyRigid(bodyNumber=o0, localPosition=[L,0,0])) #local position L = beam tip
 76    mbs.AddLoad(Torque(markerNumber = mANCFrigid, loadVector = [0, 0, E*I*0.25]))
 77
 78    #mbs.systemData.Info()
 79
 80else: #treat n elements
 81    nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 82    nElements = 8 #16
 83    lElem = L / nElements
 84    for i in range(nElements):
 85        nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
 86        elem=mbs.AddObject(Cable2D(physicsLength=lElem,
 87                                   physicsMassPerLength=rho*A,
 88                                   physicsBendingStiffness=E*I,
 89                                   physicsAxialStiffness=E*A,
 90                                   #useReducedOrderIntegration=True,
 91                                   nodeNumbers=[int(nc0)+i,int(nc0)+i+1]))
 92        cableList+=[elem]
 93
 94    mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))
 95    mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))
 96    mANCF3 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3))
 97
 98    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
 99    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
100    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF3]))
101
102    #mANCFLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=nLast)) #force
103    #mbs.AddLoad(Force(markerNumber = mANCFLast, loadVector = [0, -1e8, 0])) #will be changed in load steps
104    #mANCFrigid = mbs.AddMarker(MarkerBodyRigid(bodyNumber=elem, localPosition=[lElem,0,0])) #local position L = beam tip
105    #mbs.AddLoad(Torque(markerNumber = mANCFrigid, loadVector = [0, 0, E*I*0.25*mypi]))
106    mANCFnode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nLast)) #local position L = beam tip
107    mbs.AddLoad(Torque(markerNumber = mANCFnode, loadVector = [0, 0, 0.5*E*I*mypi]))
108    #mbs.AddLoad(Force(markerNumber = mANCFnode, loadVector = [0, 1e3, 0]))
109
110
111mbs.Assemble()
112#print(mbs)
113
114simulationSettings = exu.SimulationSettings() #takes currently set values or default values
115#simulationSettings.solutionSettings.coordinatesSolutionFileName = 'ANCFCable2Dbending' + str(nElements) + '.txt'
116
117fact = 1000
118simulationSettings.timeIntegration.numberOfSteps = 1*fact
119simulationSettings.timeIntegration.endTime = 0.002*fact
120simulationSettings.solutionSettings.writeSolutionToFile = True
121simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/fact
122simulationSettings.displayComputationTime = False
123simulationSettings.displayStatistics = True
124simulationSettings.timeIntegration.verboseMode = 1
125
126simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = True
127simulationSettings.timeIntegration.generalizedAlpha.useNewmark = True
128simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6 #0.6 works well
129simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = False
130
131#SC.visualizationSettings.nodes.showNumbers = True
132SC.visualizationSettings.bodies.showNumbers = False
133#SC.visualizationSettings.connectors.showNumbers = True
134SC.visualizationSettings.nodes.defaultSize = 0.05
135
136simulationSettings.solutionSettings.solutionInformation = "ANCF cable with imposed curvature or applied tip force/torque"
137
138solveDynamic = True
139if solveDynamic:
140    exu.StartRenderer()
141
142    def UFchangeLoad(mbs, t):
143        tt=t
144        if tt > 1:
145            tt=1
146        #mbs.SetLoadParameter(0, 'loadVector', [0, 1e6*tt, 0]) #for force
147        mbs.SetLoadParameter(0, 'loadVector', [0, 0, 2*0.5*E*I*mypi*tt])
148
149        #print('t=',tt,'p=',mbs.GetNodeOutput(nLast, exu.OutputVariableType.Position))
150        return True #True, means that everything is alright, False=stop simulation
151
152    mbs.SetPreStepUserFunction(UFchangeLoad)
153
154
155
156    mbs.SolveDynamic(simulationSettings)
157
158    SC.WaitForRenderEngineStopFlag()
159    exu.StopRenderer() #safely close rendering window!
160
161else:
162    simulationSettings.staticSolver.verboseMode = 1
163    #simulationSettings.staticSolver.loadStepGeometric = True;
164    #.staticSolver.loadStepGeometricRange = 1e2;
165
166    exu.StartRenderer()
167
168    #manual load stepping
169    doLoadStepping = False
170    if doLoadStepping:
171        nLoadSteps = 40;
172        for loadSteps in range(nLoadSteps):
173            loadFact = ((loadSteps+1)/nLoadSteps)
174            simulationSettings.staticSolver.loadStepStart = loadFact
175            simulationSettings.staticSolver.newton.relativeTolerance = 1e-8*loadFact #10000
176
177            loadDict = mbs.GetLoad(0)
178            loadDict['loadVector'] = [0, 0, E*I/L*2*mypi*loadFact]
179            mbs.ModifyLoad(0, loadDict)
180
181            #prescribe curvature:
182            #curvatureValue = 2*((loadSteps+1)/nLoadSteps)
183            #print('curvature='+str(curvatureValue))
184
185            #for nCable in cableList:
186            #    cableDict = mbs.GetObject(nCable)
187            #    cableDict['physicsReferenceCurvature'] = curvatureValue
188            #    cableDict['physicsReferenceAxialStrain'] = 0.1*curvatureValue
189            #    mbs.ModifyObject(nCable, cableDict)
190
191            mbs.SolveStatic(simulationSettings)
192
193            sol = mbs.systemData.GetODE2Coordinates()
194            mbs.systemData.SetODE2Coordinates(coordinates=sol, configurationType=exu.ConfigurationType.Initial) #set initial conditions for next step
195
196            print('sol step  ' + str(loadSteps) + ':')
197            n = len(sol)
198            print('tip displacement: x='+str(sol[n-4])+', y='+str(sol[n-3]))
199            n2 = int(len(sol)/8)
200            print('mid displacement: x='+str(sol[n2*4])+', y='+str(sol[n2*4+1]))
201
202    else:
203        simulationSettings.staticSolver.numberOfLoadSteps  = 8
204        simulationSettings.staticSolver.newton.relativeTolerance = 1e-7
205        simulationSettings.staticSolver.verboseMode = 1
206        simulationSettings.displayStatistics = True
207        mbs.SolveStatic(simulationSettings)
208
209
210    SC.WaitForRenderEngineStopFlag()
211    exu.StopRenderer() #safely close rendering window!