ANCFcantileverTestDyn.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test with ANCF cantilever beam; excitation with a coordinate constraint which is changed by preStepExecute function
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2019-10-25
  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
 16import numpy as np
 17
 18SC = exu.SystemContainer()
 19mbs = SC.AddSystem()
 20
 21
 22#background
 23rect = [-10,-10,10,10] #xmin,ymin,xmax,ymax
 24background = {'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= [background])))
 26
 27#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 28#cable:
 29
 30L=2                    # length of ANCF element in m
 31E=2.07e11             # Young's modulus of ANCF element in N/m^2
 32rho=7800               # density of ANCF element in kg/m^3
 33b=0.1                  # width of rectangular ANCF element in m
 34h=0.1                  # height of rectangular ANCF element in m
 35A=b*h                  # cross sectional area of ANCF element in m^2
 36I=b*h**3/12            # second moment of area of ANCF element in m^4
 37f=3*E*I/L**2           # tip load applied to ANCF element in N
 38
 39print("load f="+str(f))
 40
 41nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 42mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 43
 44mode = 1
 45if mode==0: #treat one element
 46    nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 47    nc1 = mbs.AddNode(Point2DS1(referenceCoordinates=[L,0,1,0]))
 48    o0 = mbs.AddObject(Cable2D(physicsLength=L, physicsMassPerLength=rho*A, physicsBendingStiffness=E*I, physicsAxialStiffness=E*A, nodeNumbers=[nc0,nc1]))
 49    print(mbs.GetObject(o0))
 50
 51    mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))
 52    mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))
 53    mANCF2b = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3))
 54
 55    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
 56    ccy=mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]),offset=1e-6)
 57    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2b]))
 58
 59    mANCFnode = mbs.AddMarker(MarkerNodePosition(nodeNumber=nc1)) #force
 60    mbs.AddLoad(Force(markerNumber = mANCFnode, loadVector = [0, 0, 0]))
 61
 62
 63else: #treat n elements
 64    nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 65    nElements = 16
 66    lElem = L / nElements
 67    for i in range(nElements):
 68        nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
 69        mbs.AddObject(Cable2D(physicsLength=lElem, physicsMassPerLength=rho*A,
 70                              physicsBendingStiffness=E*I, physicsAxialStiffness=E*A, nodeNumbers=[int(nc0)+i,int(nc0)+i+1]))
 71
 72    mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))
 73    mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))
 74    mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3))
 75
 76    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
 77    ccy=mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
 78    mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2]))
 79
 80    #mANCFLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=nLast)) #force
 81    #nl=mbs.AddLoad(Force(markerNumber = mANCFLast, loadVector = [0, -f*0.01, 0])) #will be changed in load steps
 82
 83
 84
 85mbs.Assemble()
 86print(mbs)
 87
 88
 89
 90simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 91
 92def UFexcitation(mbs, t):
 93    mbs.SetObjectParameter(ccy, 'offset', 0.1*np.sin(2*np.pi*20*t))
 94    return True #True, means that everything is alright, False=stop simulation
 95
 96mbs.SetPreStepUserFunction(UFexcitation)
 97
 98
 99fact = 20000
100simulationSettings.timeIntegration.numberOfSteps = 1*fact
101simulationSettings.timeIntegration.endTime = 0.000025*fact
102simulationSettings.solutionSettings.writeSolutionToFile = True
103simulationSettings.solutionSettings.solutionWritePeriod = simulationSettings.timeIntegration.endTime/fact
104simulationSettings.displayComputationTime = True
105simulationSettings.timeIntegration.verboseMode = 1
106
107simulationSettings.timeIntegration.newton.relativeTolerance = 1e-8*1000 #10000
108simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-10*100
109
110simulationSettings.timeIntegration.newton.useModifiedNewton = True
111simulationSettings.timeIntegration.newton.maxModifiedNewtonIterations = 5
112simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1
113simulationSettings.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 6.055454452393343e-06*0.1 #eps^(1/3)
114simulationSettings.timeIntegration.newton.modifiedNewtonContractivity = 1000
115simulationSettings.timeIntegration.generalizedAlpha.useIndex2Constraints = False
116simulationSettings.timeIntegration.generalizedAlpha.useNewmark = False
117simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.6
118simulationSettings.displayStatistics = True
119simulationSettings.displayComputationTime = True
120
121#SC.visualizationSettings.nodes.showNumbers = True
122SC.visualizationSettings.bodies.showNumbers = False
123#SC.visualizationSettings.connectors.showNumbers = True
124SC.visualizationSettings.nodes.defaultSize = 0.01
125
126simulationSettings.solutionSettings.solutionInformation = "nonlinear beam oscillations"
127
128exu.StartRenderer()
129mbs.SolveDynamic(simulationSettings)
130SC.WaitForRenderEngineStopFlag()
131exu.StopRenderer() #safely close rendering window!