ANCFtestHalfcircle.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  ANCF Cable2D cantilever bent into a half circle; uses multiple static computations
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2019-09-01
  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 *
 15
 16SC = exu.SystemContainer()
 17mbs = SC.AddSystem()
 18
 19
 20#background
 21rect = [-2,-2,2,2] #xmin,ymin,xmax,ymax
 22background0 = {'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
 23background1 = {'type':'Circle', 'radius': 0.1, 'position': [-1.5,0,0]}
 24background2 = {'type':'Text', 'position': [-1,-1,0], 'text':'Example with text\nin two lines:.=!'} #background
 25oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background0, background1, background2])))
 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               # density of ANCF element in kg/m^3
 35b=0.01                 # width of rectangular ANCF element in m; solver has problems with h=0.1 and nElem>8
 36h=0.01                 # 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    o0 = mbs.AddObject(Cable2D(physicsLength=L, physicsMassPerLength=rho*A, physicsBendingStiffness=E*I, physicsAxialStiffness=E*A, nodeNumbers=[nc0,nc1]))
 57    cableList+=[o0]
 58    #print(mbs.GetObject(o0))
 59
 60    mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))
 61    mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))
 62    mANCF2b = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3))
 63
 64    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
 65    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
 66    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2b]))
 67
 68    #mANCFnode = mbs.AddMarker(MarkerNodePosition(nodeNumber=nc1)) #force
 69    #mbs.AddLoad(Force(markerNumber = mANCFnode, loadVector = [0, -10000, 0]))
 70    mANCFrigid = mbs.AddMarker(MarkerBodyRigid(bodyNumber=o0, localPosition=[L,0,0])) #local position L = beam tip
 71    mbs.AddLoad(Torque(markerNumber = mANCFrigid, loadVector = [0, 0, E*I*0.25]))
 72
 73
 74else: #treat n elements
 75    nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 76    nElements = 8*32 #2020-01-02: 64 elements; works now better 2020-01-02 with h=0.01; does not work with 16 elements (2019-12-07)
 77    lElem = L / nElements
 78    for i in range(nElements):
 79        nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
 80        elem=mbs.AddObject(Cable2D(physicsLength=lElem, physicsMassPerLength=rho*A,
 81                                   physicsBendingStiffness=E*I, physicsAxialStiffness=E*A, nodeNumbers=[int(nc0)+i,int(nc0)+i+1]))
 82        cableList+=[elem]
 83
 84    mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))
 85    mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))
 86    mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3))
 87
 88    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
 89    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
 90    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2]))
 91
 92    #mANCFLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=nLast)) #force
 93    #mbs.AddLoad(Force(markerNumber = mANCFLast, loadVector = [0, -100000*0, 0])) #will be changed in load steps
 94    #mANCFrigid = mbs.AddMarker(MarkerBodyRigid(bodyNumber=elem, localPosition=[lElem,0,0])) #local position L = beam tip
 95    #mbs.AddLoad(Torque(markerNumber = mANCFrigid, loadVector = [0, 0, E*I*0.25*mypi]))
 96    mANCFnode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nLast)) #local position L = beam tip
 97    mbs.AddLoad(Torque(markerNumber = mANCFnode, loadVector = [0, 0, E*I*mypi]))
 98
 99
100
101mbs.Assemble()
102print(mbs)
103
104simulationSettings = exu.SimulationSettings() #takes currently set values or default values
105
106fact = 1000
107simulationSettings.timeIntegration.numberOfSteps = 1*fact
108simulationSettings.timeIntegration.endTime = 0.001*fact
109simulationSettings.solutionSettings.writeSolutionToFile = True
110simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/fact
111simulationSettings.displayComputationTime = True
112simulationSettings.timeIntegration.verboseMode = 1
113
114simulationSettings.timeIntegration.newton.relativeTolerance = 1e-7 #10000
115#simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-8*1000
116
117simulationSettings.timeIntegration.newton.useModifiedNewton = False
118#simulationSettings.timeIntegration.newton.maxModifiedNewtonIterations = 5
119#simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
120#simulationSettings.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 6.055454452393343e-06*0.1 #eps^(1/3)
121#simulationSettings.timeIntegration.newton.modifiedNewtonContractivity = 1000
122simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = False
123simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False
124simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6 #0.6 works well
125simulationSettings.displayStatistics = True
126simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
127
128#SC.visualizationSettings.nodes.showNumbers = True
129SC.visualizationSettings.bodies.showNumbers = False
130#SC.visualizationSettings.connectors.showNumbers = True
131SC.visualizationSettings.nodes.defaultSize = 0.025
132
133simulationSettings.solutionSettings.solutionInformation = "ANCF test halfcircle"
134
135solveDynamic = False
136if solveDynamic:
137    exu.StartRenderer()
138    simulationSettings.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 1e-9*0.25
139
140    def UFchangeLoad(mbs, t):
141        mbs.SetLoadParameter(0,'loadVector',[0, 0, E*I*3.141592653589793*t])
142        return True #True, means that everything is alright, False=stop simulation
143
144    mbs.SetPreStepUserFunction(UFchangeLoad)
145
146    mbs.SolveDynamic(simulationSettings)
147    #v = mbs.CallObjectFunction(1,'GetAngularVelocity',{'localPosition':[L/2,0,0],'configuration':'Current'})
148    #print('angular vel='+str(v))
149    SC.WaitForRenderEngineStopFlag()
150    exu.StopRenderer() #safely close rendering window!
151
152else:
153    simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-9*0.25
154#    simulationSettings.staticSolver.verboseMode = 1
155#
156#    simulationSettings.staticSolver.newton.absoluteTolerance = 1e-10
157    simulationSettings.staticSolver.newton.maxIterations = 50 #for bending into circle
158
159    exu.StartRenderer()
160
161    doLoadStepping = False
162    if doLoadStepping:
163        nLoadSteps = 80 #80
164        for loadSteps in range(nLoadSteps):
165            #nLoad = 0
166            #loadValue = f**((loadSteps+1)/nLoadSteps) #geometric increment of loads
167            #print('load='+str(loadValue))
168            #loadDict = mbs.GetLoad(nLoad)
169            #loadDict['loadVector'] = [0, -loadValue,0]
170            #mbs.ModifyLoad(nLoad, loadDict)
171            loadFact = ((loadSteps+1)/nLoadSteps)
172            simulationSettings.staticSolver.currentTime = loadFact
173            simulationSettings.staticSolver.newton.relativeTolerance = 1e-8*loadFact #10000
174
175            loadDict = mbs.GetLoad(0)
176            loadDict['loadVector'] = [0, 0, E*I/L*2*mypi*loadFact]
177            mbs.ModifyLoad(0, loadDict)
178
179            #curvatureValue = 2*((loadSteps+1)/nLoadSteps) #geometric increment of loads
180            #print('curvature='+str(curvatureValue))
181
182            #for nCable in cableList:
183            #    cableDict = mbs.GetObject(nCable)
184            #    cableDict['physicsReferenceCurvature'] = curvatureValue
185            #    cableDict['physicsReferenceAxialStrain'] = 0.1*curvatureValue
186            #    mbs.ModifyObject(nCable, cableDict)
187
188            mbs.SolveStatic(simulationSettings)
189
190            sol = mbs.systemData.GetODE2Coords()
191            mbs.systemData.SetODE2Coords(coords=sol, configurationType=exu.ConfigurationType.Initial) #set initial conditions for next step
192
193            print('sol step  ' + str(loadSteps) + ':')
194            n = len(sol)
195            print('tip displacement: x='+str(sol[n-4])+', y='+str(sol[n-3]))
196            n2 = int(len(sol)/8)
197            print('mid displacement: x='+str(sol[n2*4])+', y='+str(sol[n2*4+1]))
198
199
200            #sol = mbs.systemData.GetODE2Coords(exu.ConfigurationType.Initial)
201            #print('initial values='+str(sol))
202    else:
203        simulationSettings.staticSolver.numberOfLoadSteps = 1
204        simulationSettings.staticSolver.adaptiveStep = True
205        mbs.SolveStatic(simulationSettings)
206
207
208    SC.WaitForRenderEngineStopFlag()
209    exu.StopRenderer() #safely close rendering window!