You can view and download this file on Github: computeODE2EigenvaluesTest.py
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: Test for computation of eigenvalues using utility eigensolver functionality based on scipy.linalg
#
# Author: Johannes Gerstmayr
# Date: 2020-12-18
#
# 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.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
import exudyn as exu
from exudyn.itemInterface import *
import numpy as np
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
try: #only if called from test suite
from modelUnitTests import exudynTestGlobals #for globally storing test results
useGraphics = exudynTestGlobals.useGraphics
except:
class ExudynTestGlobals:
pass
exudynTestGlobals = ExudynTestGlobals()
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SC = exu.SystemContainer()
mbs = SC.AddSystem()
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#cable:
mypi = 3.141592653589793
L=2. # length of ANCF element in m
#L=mypi # length of ANCF element in m
E=2.07e11 # Young's modulus of ANCF element in N/m^2
rho=7800 # density of ANCF element in kg/m^3
b=0.01 # width of rectangular ANCF element in m
h=0.01 # height of rectangular ANCF element in m
A=b*h # cross sectional area of ANCF element in m^2
I=b*h**3/12 # second moment of area of ANCF element in m^4
EI = E*I
rhoA = rho*A
exu.Print("EI="+str(EI))
exu.Print("rhoA="+str(rhoA))
nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
cableList=[]
nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
nElements = 32 #32
lElem = L / nElements
for i in range(nElements):
nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
elem=mbs.AddObject(Cable2D(physicsLength=lElem,
physicsMassPerLength=rho*A,
physicsBendingStiffness=E*I,
physicsAxialStiffness=E*A*0.1,
nodeNumbers=[int(nc0)+i,int(nc0)+i+1],
useReducedOrderIntegration=True))
cableList+=[elem]
mbs.Assemble()
simulationSettings = exu.SimulationSettings() #takes currently set values or default values
simulationSettings.staticSolver.verboseMode = 1
nEig = 3
[values, vectors] = mbs.ComputeODE2Eigenvalues(simulationSettings,
numberOfEigenvalues = nEig+3) #3 eigenvalues + 3 rigid body zero eigenvalues
omegaNumerical = np.sqrt(values[3:nEig+3])
exu.Print("eigenvalues=",omegaNumerical) #exclude 3 rigid body modes
#[ 83.17966459 229.28844645 449.50021798]
#analytical: bending eigenfrequency of free-free beam:
#4.7300, 7.8532, 10.9956, 14.1371, 17.2787 (cosh(beta) * cos(beta) = 1)
#find roots beta:
#from mpmath import *
#mp.dps = 16 #digits
#for i in range(10): print(findroot(lambda x: cosh(x) * cos(x) - 1, 3*i+4.7))
beta = [4.730040744862704, 7.853204624095838, 10.99560783800167, 14.13716549125746, 17.27875965739948, 20.42035224562606, 23.56194490204046, 26.70353755550819, 29.84513020910325]
omega = np.zeros(nEig)
for i in range(nEig):
omega[i] = ((beta[i]/L)**4 * (EI/rhoA))**0.5
exu.Print('omega analytical =',omega)
u = omega[0]-omegaNumerical[0]
exu.Print('omega difference=',u)
exudynTestGlobals.testError = 1e-6*(u - (-2.7613614363986017e-05)) #2021-01-04: added factor 1e-6, because of larger errors/differences in 32/64bit eigenvalue solvers; 2020-12-18: (nElements=32) -2.7613614363986017e-05
exudynTestGlobals.testResult = 1e-6*u