manualExplicitIntegrator.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  ANCF Cable2D cantilever bent with manual explicit integrator
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2020-01-08
  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
 16testIsActive = exu.sys.get('testIsActive', False)
 17
 18SC = exu.SystemContainer()
 19mbs = SC.AddSystem()
 20
 21
 22exu.Print("\n\n++++++++++++++++++++++++++\nStart EXUDYN version "+exu.config.Version()+"\n")
 23
 24#background
 25rect = [-2,-2,2,2] #xmin,ymin,xmax,ymax
 26background0 = {'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
 27background1 = {'type':'Circle', 'radius': 0.1, 'position': [-1.5,0,0]} 
 28background2 = {'type':'Text', 'position': [-1,-1,0], 'text':'Example with text\nin two lines:.=!'} #background
 29oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background0, background1, background2])))
 30
 31#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 32#cable:
 33mypi = 3.141592653589793
 34
 35L=2.                   # length of ANCF element in m
 36#L=mypi                 # length of ANCF element in m
 37E=2.07e11*1e-5              # Young's modulus of ANCF element in N/m^2
 38rho=7800               # density of ANCF element in kg/m^3
 39b=0.1                  # width of rectangular ANCF element in m
 40h=0.1                  # height of rectangular ANCF element in m
 41A=b*h                  # cross sectional area of ANCF element in m^2
 42I=b*h**3/12            # second moment of area of ANCF element in m^4
 43f=3*E*I/L**2           # tip load applied to ANCF element in N
 44
 45exu.Print("load f="+str(f))
 46exu.Print("EI="+str(E*I))
 47
 48nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 49mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 50
 51cableList=[]
 52
 53
 54
 55nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 56nElements = 4 #for tests use nElements = 4
 57lElem = L / nElements
 58nLast = 0
 59for i in range(nElements):
 60    nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
 61    elem=mbs.AddObject(Cable2D(length=lElem, massPerLength=rho*A, 
 62                               bendingStiffness=E*I, axialStiffness=E*A*0.1, 
 63                               nodeNumbers=[int(nc0)+i,int(nc0)+i+1], useReducedOrderIntegration=True))
 64    cableList+=[elem]
 65
 66mANCF0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=0))
 67mANCF1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=1))
 68mANCF2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nc0, coordinate=3))
 69    
 70#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
 71#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
 72#mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2]))
 73
 74mANCFLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=nLast)) #force
 75mbs.AddLoad(Force(markerNumber = mANCFLast, loadVector = [0, -1000, 0])) #will be changed in load steps
 76#mANCFrigid = mbs.AddMarker(MarkerBodyRigid(bodyNumber=elem, localPosition=[lElem,0,0])) #local position L = beam tip
 77#mbs.AddLoad(Torque(markerNumber = mANCFrigid, loadVector = [0, 0, E*I*0.25*mypi]))
 78#mANCFnode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nLast)) #local position L = beam tip
 79#mbs.AddLoad(Torque(markerNumber = mANCFnode, loadVector = [0, 0, 0.4*E*I*0.25*mypi]))
 80#mbs.AddLoad(Force(markerNumber = mANCFnode, loadVector = [0, 0.4*E*I*0.25*mypi,0]))
 81
 82
 83
 84mbs.Assemble()
 85#exu.Print(mbs)
 86
 87simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 88
 89
 90#SC.visualizationSettings.bodies.showNumbers = False
 91SC.visualizationSettings.nodes.defaultSize = 0.025
 92dSize=0.01
 93SC.visualizationSettings.bodies.defaultSize = [dSize, dSize, dSize]
 94
 95#simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-9
 96simulationSettings.staticSolver.verboseMode = 1
 97simulationSettings.staticSolver.verboseModeFile = 2
 98simulationSettings.solution.solverInformationFileName = 'solution/solverInformation.txt'
 99
100#simulationSettings.staticSolver.newton.absoluteTolerance = 1e-8
101simulationSettings.staticSolver.newton.relativeTolerance = 1e-6 #1e-5 works for 64 elements
102simulationSettings.staticSolver.newton.maxIterations = 20 #50 for bending into circle
103    
104if not testIsActive: #only start graphics once, but after background is set
105    SC.renderer.Start()
106
107simulationSettings.staticSolver.numberOfLoadSteps = 10
108simulationSettings.staticSolver.adaptiveStep = True
109
110import numpy as np
111
112testRefVal = 0
113#compute eigenvalues manually:
114calcEig = True
115if calcEig:
116    from scipy.linalg import solve, eigh, eig #eigh for symmetric matrices, positive definite
117
118    staticSolver = exu.MainSolverStatic()
119    #staticSolver.SolveSystem(mbs, simulationSettings)
120        
121    staticSolver.InitializeSolver(mbs, simulationSettings)
122
123    staticSolver.ComputeMassMatrix(mbs)
124    m = staticSolver.GetSystemMassMatrix()
125    #exu.Print("m =",m)
126
127    staticSolver.ComputeJacobianODE2RHS(mbs, scalarFactor_ODE2=-1, scalarFactor_ODE2_t=0)
128    staticSolver.ComputeJacobianAE(mbs)
129    K = staticSolver.GetSystemJacobian()
130    #exu.Print("K =",K)
131    nODE2 = staticSolver.GetODE2size()
132
133
134    K2 = K[0:nODE2,0:nODE2]
135
136    [eigvals, eigvecs] = eigh(K2, m) #this gives omega^2 ... squared eigen frequencies (rad/s)
137    ev = np.sort(a=abs(eigvals))
138    #exu.Print("ev =",ev)
139    if (len(ev) >= 7):
140        f6 = np.sqrt(abs(ev[6]))/(2*np.pi)
141        exu.Print("ev=", f6)
142        testRefVal += f6 #first bending eigenmode
143
144    staticSolver.FinalizeSolver(mbs, simulationSettings)
145
146#++++++++++++++++++++++++++++++++++++++++++++++++++
147#TEST
148def UserFunctionInitializeStep(mainSolver, mainSys, sims):
149    #exu.Print("t=", mainSolver.it.currentTime)
150    mainSolver.UpdateCurrentTime(mainSys, sims)
151    mainSys.systemData.SetTime(mainSolver.it.currentTime);
152    return True
153
154#test for explicit integrator:
155def UserFunctionNewton(mainSolver, mainSys, sims):
156
157    nODE2 = mainSolver.GetODE2size()
158    nAE = mainSolver.GetAEsize()
159    #nSys = nODE2+nAE
160    #print("u=", mainSys.systemData.GetODE2Coordinates())
161    dynamicSolver.ComputeODE2RHS(mbs)
162    res = dynamicSolver.GetSystemResidual()
163    Fode2 = res[0:nODE2]
164    #print("res=", Fode2)
165
166    dynamicSolver.ComputeMassMatrix(mbs)
167    M = dynamicSolver.GetSystemMassMatrix()
168    a = np.linalg.solve(M,Fode2) #acceleration
169
170    h = dynamicSolver.it.currentStepSize
171
172    u0 = mainSys.systemData.GetODE2Coordinates()
173    v0 = mainSys.systemData.GetODE2Coordinates_t()
174
175    mainSys.systemData.SetODE2Coordinates(u0+h*v0)
176    mainSys.systemData.SetODE2Coordinates_t(v0+h*a)
177
178    return True
179
180dynamicSolver = exu.MainSolverImplicitSecondOrder()
181
182simulationSettings.timeIntegration.numberOfSteps = 5000 #1000 steps for test suite/error
183simulationSettings.timeIntegration.endTime = 0.05              #1s for test suite / error
184
185simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
186#simulationSettings.show.computationTime = True
187simulationSettings.timeIntegration.verboseMode = 1
188
189#dynamicSolver.SetUserFunctionInitializeStep(mbs, UserFunctionInitializeStep)
190dynamicSolver.SetUserFunctionNewton(mbs, UserFunctionNewton)
191
192dynamicSolver.SolveSystem(mbs, simulationSettings)
193#mbs.SolveDynamic(simulationSettings)
194
195uy=mbs.GetNodeOutput(nLast,exu.OutputVariableType.Position)[1] #y-coordinate of tip
196exu.Print("uy=", uy)
197exu.Print("testResult=", testRefVal + uy)
198exu.sys['testResult'] = testRefVal + uy
199
200if not testIsActive: #only start graphics once, but after background is set
201    SC.renderer.DoIdleTasks()
202    SC.renderer.Stop() #safely close rendering window!