computeODE2EigenvaluesTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for computation of eigenvalues using utility eigensolver functionality based on scipy.linalg
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2020-12-18
  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 *
 15import numpy as np
 16
 17#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 18#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 19try: #only if called from test suite
 20    from modelUnitTests import exudynTestGlobals #for globally storing test results
 21    useGraphics = exudynTestGlobals.useGraphics
 22except:
 23    class ExudynTestGlobals:
 24        pass
 25    exudynTestGlobals = ExudynTestGlobals()
 26#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 27
 28SC = exu.SystemContainer()
 29mbs = SC.AddSystem()
 30
 31
 32#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 33#cable:
 34mypi = 3.141592653589793
 35
 36L=2.                   # length of ANCF element in m
 37#L=mypi                 # length of ANCF element in m
 38E=2.07e11              # Young's modulus of ANCF element in N/m^2
 39rho=7800               # density of ANCF element in kg/m^3
 40b=0.01                  # width of rectangular ANCF element in m
 41h=0.01                  # height of rectangular ANCF element in m
 42A=b*h                  # cross sectional area of ANCF element in m^2
 43I=b*h**3/12            # second moment of area of ANCF element in m^4
 44EI = E*I
 45rhoA = rho*A
 46
 47exu.Print("EI="+str(EI))
 48exu.Print("rhoA="+str(rhoA))
 49
 50nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
 51mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
 52
 53cableList=[]
 54
 55
 56
 57nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
 58nElements = 32 #32
 59lElem = L / nElements
 60for i in range(nElements):
 61    nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
 62    elem=mbs.AddObject(Cable2D(physicsLength=lElem,
 63                               physicsMassPerLength=rho*A,
 64                               physicsBendingStiffness=E*I,
 65                               physicsAxialStiffness=E*A*0.1,
 66                               nodeNumbers=[int(nc0)+i,int(nc0)+i+1],
 67                               useReducedOrderIntegration=True))
 68    cableList+=[elem]
 69
 70
 71mbs.Assemble()
 72
 73simulationSettings = exu.SimulationSettings() #takes currently set values or default values
 74
 75simulationSettings.staticSolver.verboseMode = 1
 76
 77nEig = 3
 78[values, vectors] = mbs.ComputeODE2Eigenvalues(simulationSettings,
 79                                               numberOfEigenvalues = nEig+3)    #3 eigenvalues + 3 rigid body zero eigenvalues
 80
 81
 82omegaNumerical = np.sqrt(values[3:nEig+3])
 83exu.Print("eigenvalues=",omegaNumerical) #exclude 3 rigid body modes
 84#[ 83.17966459 229.28844645 449.50021798]
 85
 86#analytical: bending eigenfrequency of free-free beam:
 87#4.7300, 7.8532, 10.9956, 14.1371, 17.2787 (cosh(beta) * cos(beta) = 1)
 88#find roots beta:
 89#from mpmath import *
 90#mp.dps = 16 #digits
 91#for i in range(10): print(findroot(lambda x: cosh(x) * cos(x) - 1, 3*i+4.7))
 92beta = [4.730040744862704, 7.853204624095838, 10.99560783800167, 14.13716549125746, 17.27875965739948, 20.42035224562606, 23.56194490204046, 26.70353755550819, 29.84513020910325]
 93omega = np.zeros(nEig)
 94for i in range(nEig):
 95    omega[i] = ((beta[i]/L)**4 * (EI/rhoA))**0.5
 96
 97exu.Print('omega analytical =',omega)
 98u = omega[0]-omegaNumerical[0]
 99exu.Print('omega difference=',u)
100
101exudynTestGlobals.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
102exudynTestGlobals.testResult = 1e-6*u