solverFunctionsTestEigenvalues.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test of static solver, computing eigenvalues and showing eigenmodes; uses scipy.linalg
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2020-01-06
  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
 20exu.Print("\n\n++++++++++++++++++++++++++\nStart EXUDYN version "+exu.GetVersionString()+"\n")
 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
 25background1 = {'type':'Circle', 'radius': 0.1, 'position': [-1.5,0,0]}
 26background2 = {'type':'Text', 'position': [-1,-1,0], 'text':'Example with text\nin two lines:.=!'} #background
 27oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0], visualization=VObjectGround(graphicsData= [background0, background1, background2])))
 28
 29#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 30#cable:
 31mypi = 3.141592653589793
 32
 33L=2.                   # length of ANCF element in m
 34#L=mypi                 # length of ANCF element in m
 35E=2.07e11              # Young's modulus of ANCF element in N/m^2
 36rho=7800               # density of ANCF element in kg/m^3
 37b=0.01                  # width of rectangular ANCF element in m
 38h=0.01                  # height of rectangular ANCF element in m
 39A=b*h                  # cross sectional area of ANCF element in m^2
 40I=b*h**3/12            # second moment of area of ANCF element in m^4
 41EI = E*I
 42rhoA = rho*A
 43f=3*EI/L**2           # tip load applied to ANCF element in N
 44
 45exu.Print("load f="+str(f))
 46exu.Print("EI="+str(EI))
 47exu.Print("rhoA="+str(rhoA))
 48
 49nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 50mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 51
 52cableList=[]
 53
 54
 55
 56nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 57nElements = 32 #2020-01-03: now works even with 64 elements if relTol=1e-5; did not work with 16 elements (2019-12-07)
 58lElem = L / nElements
 59for i in range(nElements):
 60    nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
 61    elem=mbs.AddObject(Cable2D(physicsLength=lElem, physicsMassPerLength=rho*A,
 62                               physicsBendingStiffness=E*I, physicsAxialStiffness=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
 70mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF0]))
 71mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF1]))
 72mbs.AddObject(CoordinateConstraint(markerNumbers=[mGround,mANCF2]))
 73
 74#mANCFLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=nLast)) #force
 75#nLoad = mbs.AddLoad(Force(markerNumber = mANCFLast, loadVector = [-450, 0, 0])) #will be changed in load steps
 76
 77mANCFrigid = mbs.AddMarker(MarkerBodyRigid(bodyNumber=elem, localPosition=[lElem,0,0])) #local position L = beam tip
 78mbs.AddLoad(Torque(markerNumber = mANCFrigid, loadVector = [0, 0, E*I*0.25*mypi]))
 79
 80#mANCFnode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nLast)) #local position L = beam tip
 81#mbs.AddLoad(Torque(markerNumber = mANCFnode, loadVector = [0, 0, 4*E*I*0.25*mypi]))
 82#mbs.AddLoad(Force(markerNumber = mANCFnode, loadVector = [0, 0.4*E*I*0.25*mypi,0]))
 83
 84
 85
 86mbs.Assemble()
 87#exu.Print(mbs)
 88
 89simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 90
 91
 92#SC.visualizationSettings.bodies.showNumbers = False
 93SC.visualizationSettings.nodes.defaultSize = 0.025
 94
 95#simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-9
 96simulationSettings.staticSolver.verboseMode = 1
 97#simulationSettings.staticSolver.verboseModeFile = 0
 98
 99#simulationSettings.staticSolver.newton.absoluteTolerance = 1e-8
100simulationSettings.staticSolver.newton.relativeTolerance = 1e-6 #1e-5 works for 64 elements
101simulationSettings.staticSolver.newton.maxIterations = 20 #50 for bending into circle
102#simulationSettings.displayComputationTime = True
103
104
105exu.StartRenderer()
106
107simulationSettings.staticSolver.numberOfLoadSteps = 100
108simulationSettings.staticSolver.adaptiveStep = True
109
110staticSolver = exu.MainSolverStatic()
111#staticSolver.SolveSystem(mbs, simulationSettings)
112#print(staticSolver.timer)
113
114import numpy as np
115from scipy.linalg import eigh, eig #eigh for symmetric matrices, positive definite
116
117#+++++++++++++++++++++++++++++++++++++
118#compute eigenvalue problem:
119
120staticSolver.InitializeSolver(mbs, simulationSettings)
121#staticSolver.SolveSteps(mbs, simulationSettings) #if preloaded
122#staticSolver.FinalizeSolver(mbs, simulationSettings)
123
124
125nODE2 = staticSolver.GetODE2size()
126
127#raise ValueError("")
128#compute mass matrix:
129staticSolver.ComputeMassMatrix(mbs, 1)#simulationSettings)
130m = staticSolver.GetSystemMassMatrix()
131#print("m =",m)
132
133#compute stiffness matrix (systemJacobian is larger!)
134staticSolver.ComputeJacobianODE2RHS(mbs, scalarFactor_ODE2=-1,scalarFactor_ODE2_t=0,scalarFactor_ODE1=0)
135staticSolver.ComputeJacobianAE(mbs, 1)
136K = staticSolver.GetSystemJacobian()
137#print("K =",K)
138
139K2 = K[0:nODE2,0:nODE2]
140
141[eigvals, eigvecs] = eigh(K2, m) #this gives omega^2 ... squared eigen frequencies (rad/s)
142ev = np.sort(a=abs(eigvals)) #there may be very small eigenvalues
143print('eigvals=',eigvals)
144
145nEig = 4
146for i in range(len(ev)):
147    ev[i] = ev[i]**0.5
148
149print("omega numerical =",ev[3:3+nEig])
150
151
152#analytical: bending eigenfrequency of free-free beam:
153#4.7300, 7.8532, 10.9956, 14.1371, 17.2787 (cosh(beta) * cos(beta) = 1)
154#find roots beta:
155#from mpmath import *
156#mp.dps = 16 #digits
157#for i in range(10): print(findroot(lambda x: cosh(x) * cos(x) - 1, 3*i+4.7))
158beta = [4.730040744862704, 7.853204624095838, 10.99560783800167, 14.13716549125746, 17.27875965739948, 20.42035224562606, 23.56194490204046, 26.70353755550819, 29.84513020910325]
159omega = np.zeros(nEig)
160for i in range(nEig):
161    omega[i] = ((beta[i]/L)**4 * (EI/rhoA))**0.5
162
163print('omega analytical =',omega)
164
165
166#mbs.SolveStatic(simulationSettings)
167
168
169SC.WaitForRenderEngineStopFlag()
170exu.StopRenderer() #safely close rendering window!