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
18SC = exu.SystemContainer()
19mbs = SC.AddSystem()
20
21
22#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
23#cable:
24mypi = 3.141592653589793
25
26L=2.                   # length of ANCF element in m
27#L=mypi                 # length of ANCF element in m
28E=2.07e11              # Young's modulus of ANCF element in N/m^2
29rho=7800               # density of ANCF element in kg/m^3
30b=0.01                  # width of rectangular ANCF element in m
31h=0.01                  # height of rectangular ANCF element in m
32A=b*h                  # cross sectional area of ANCF element in m^2
33I=b*h**3/12            # second moment of area of ANCF element in m^4
34EI = E*I
35rhoA = rho*A
36
37exu.Print("EI="+str(EI))
38exu.Print("rhoA="+str(rhoA))
39
40nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
41mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber = nGround, coordinate=0)) #Ground node ==> no action
42
43cableList=[]
44
45
46
47nc0 = mbs.AddNode(Point2DS1(referenceCoordinates=[0,0,1,0]))
48nElements = 32 #32
49lElem = L / nElements
50for i in range(nElements):
51    nLast = mbs.AddNode(Point2DS1(referenceCoordinates=[lElem*(i+1),0,1,0]))
52    elem=mbs.AddObject(Cable2D(length=lElem, 
53                               massPerLength=rho*A, 
54                               bendingStiffness=E*I, 
55                               axialStiffness=E*A*0.1, 
56                               nodeNumbers=[int(nc0)+i,int(nc0)+i+1], 
57                               useReducedOrderIntegration=True))
58    cableList+=[elem]
59
60
61mbs.Assemble()
62
63simulationSettings = exu.SimulationSettings() #takes currently set values or default values
64
65simulationSettings.staticSolver.verboseMode = 1
66
67nEig = 3
68[values, vectors] = mbs.ComputeODE2Eigenvalues(simulationSettings, 
69                                               numberOfEigenvalues = nEig+3)    #3 eigenvalues + 3 rigid body zero eigenvalues
70
71
72omegaNumerical = np.sqrt(values[3:nEig+3])
73exu.Print("eigenvalues=",omegaNumerical) #exclude 3 rigid body modes
74#[ 83.17966459 229.28844645 449.50021798] 
75
76#analytical: bending eigenfrequency of free-free beam:
77#4.7300, 7.8532, 10.9956, 14.1371, 17.2787 (cosh(beta) * cos(beta) = 1)
78#find roots beta:
79#from mpmath import *
80#mp.dps = 16 #digits
81#for i in range(10): exu.Print(findroot(lambda x: cosh(x) * cos(x) - 1, 3*i+4.7))
82beta = [4.730040744862704, 7.853204624095838, 10.99560783800167, 14.13716549125746, 17.27875965739948, 20.42035224562606, 23.56194490204046, 26.70353755550819, 29.84513020910325]
83omega = np.zeros(nEig)
84for i in range(nEig):
85    omega[i] = ((beta[i]/L)**4 * (EI/rhoA))**0.5
86
87exu.Print('omega analytical =',omega)
88u = omega[0]-omegaNumerical[0]
89exu.Print('omega difference=',u)
90
91exu.sys['testResult'] = 1e-6*u