Eigenfrequencies from complex eigenvalues do not match to the eigenfrequencies from "convert2frequencies"

[SimonAndreas]

Hi,

I am performing a modal analysis on a crane structure using the function "ComputeODE2Eigenvalues".
My model consists of two beams, based on 2D beam elements, a rigid body, generic joints, distance constraints, a springDamper element and coordinate constraints.

The function seems to work well when I have "computeComplexEigenvalues = False" and "convert2Frequencies = True". From this I get the natural frequencies directly and sorted.

However, when I want to obtain complex eigenvalues, using "computeComplexEigenvalues = True", "convert2Frequencies = False" and "useAbsoluteValues = False", I get different eigenfrequencies.

I convert the complex eigenvalues to natural frequencies in [Hz] as follows: f_n = np.abs(complexEigenValue)/(2*np.pi). I sort the eigenfrequencies as well.
I use the same approach as in the "complexEigenvaluesTest.py" (TM). In this TestModel, the frequencies DO match between the two different methods.

1. Is there a certain object or constraint in my model that could cause this difference between the complex and real analysis with "ComputeODE2Eigenvalues"?
2. Is my method to calculate the natural frequencies from the complex eigenvalues correct?

[jgerstmayr]

Your description is sound, however, with the limited amount of information, I cannot see the exact problem.

Make sure that you correctly compute the frequency because
f_n = np.abs(complexEigenValue)/(2*np.pi)
does not give the frequency as you should rather use eigenValues[0].imag, see the example.

Another thing is the different sorting you may have to consider.

Provide more output, exudyn version, etc., then it should become more clear. If your system matrices are M,K, D, the two approaches should match well.

In any case, you can try ignoreAlgebraicEquations=True, in order to see if problems also occur in case that you remove constraints.

[SimonAndreas]

Thank you for your answer. I tried to calculate the frequency in different ways, as you proposed, but this does not solve my issue.

I heavily simplified my model and managed to recreate the problem. I am running Exudyn version 1.9.0
I have pasted my code below, my apologies for its length.

I get the following output:
[ 3.716 33.284 36.123 126.193 218.473]
[ 3.783 36.123 58.799 126.193 284.063]
Where the eigenfrequencies in the first list are calculated using "computeComplexEigenvalues = False" and the eigenfrequencies in the second list are calculated using "computeComplexEigenvalues = True".

As you can see, the axial modes are exactly the same (36.123Hz and 126.193Hz) but all the bending modes are different. So it seems like an issue with the bending modes of the "ObjectBeamGeometricallyExact2D". Namely, the issue did not occur in your example "complexEigenvaluesTest.py" (TM) either.

Am I making a mistake or is there an issue with this element type when using the complex modal analysis tool?

import sys
import math
import copy
import exudyn as exu
from exudyn.utilities import *

import matplotlib.pyplot as plt
import matplotlib.ticker as ticker

import numpy as np

IMain=0.01906*0.8     # second moment of area of beam element in m^4 --- stiffCR is a correction factor that RESULTS IN A TIP DEFLECTION OF 1.7mm AT 1kN TIP LOAD
ITB=0.003062*0.4     # second moment of area of beam element in m^4 --- stiffCR is a correction factor that RESULTS IN A TIP DEFLECTION OF 1.7mm AT 1kN TIP LOAD
mMain = 12182       # mass of gangway Main boom
mTB = 4882       # mass of gangway TB
mTip = 676-24          # mass of Tip (modelled as rigid body) (!!!Subtract 24kg due to railing and rubber blocks missing during testing!!!)
kLw = 8.82E6*0.8 # k*L is constant for a wire rope. [N]
Dc_ratio_Main = 0.05      # Damping coefficient of beam elements in main boom as a ratio of the stiffness      Old Dc (not ratio) Main boom: 140E6
Dc_ratio_TB = 0.05      # Damping coefficient of beam elements in T-boom as a ratio of the stiffness      Old Dc (not ratio) T-boom: TB: 26E6
Dc_wLw = 10000      # Damping coefficient * Lw for wire rope (Note, this value changes with Lw, see implementation)     Old Dc (not ratio) wire rope: 850
Dc_ratio_luff = 0.02   # Damping coefficient of both luffing cylinders combined        Old Dc (not ratio) luffing cylinder: 2000000
Dc_ratio_Tel = 0.10        # Damping coefficient RATIO of both combined telescoping wires            Old Dc (not ratio) telescoping system: 2000
AMain=0.02160             # cross sectional area of beam element in m^2
ATB=0.01157             # cross sectional area of beam element in m^2
Eoil = 1.00E9*1.4           # Bulk modulus of hydraulic fluid [Pa]
EAbelt = 19261E3      # EA telescoping belt [N]
## set up mbs
SC = exu.SystemContainer()
mbs = SC.AddSystem()

## define general parameters for model
strokeTel = 8.8      # Stroke length of T-boom
L_tip = strokeTel + 19.9         # gangway length until end of T-boom MEAN: 25.4 UNTIL WINCH
E=200E9              # Young's modulus of beam element in N/m^2
nu = 0.3               # Poisson's ratio
# ks = 10*(1+nu)/(12+11*nu)
G = E/(2*(1+nu))        # shear modulus
g = 9.81
phi = np.deg2rad(0)     # Luffing angle in [rad]

## define main boom parameters
nElementsMain=2           # number of beam finite elements 
LMain=17.875              # total length of beam 
# lElemMain = LMain/nElementsMain
rhoMain = mMain/LMain/AMain     # density of beam element in kg/m^3
EIMain = E*IMain
EAMain = E*AMain
rhoIMain = rhoMain*IMain
GAMain = G*AMain   



## visualization parameters
hMain = 0.5 # m Height of main boom


#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## add ground node, markers
nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))    
mNCground = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
mNRground = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))      # define ground marker to constrain moving deck node

## --- Create the T-boom --- (T-boom is created first for visualization purposes)
## create nodes with for loop
nodeListTB=[]
pRefListTB=[]
elementListTB=[]
nodeListWheelsTB = []


## create the Main boom
## create nodes with for loop
nodeListMain=[]
pRefListMain=[]
elementListMain=[]
nodeListWheelsMain = []

nodePositionsMain = np.arange(0,LMain+LMain/nElementsMain,LMain/nElementsMain)

for i in range(len(nodePositionsMain)):
    p1Ref = [cos(phi)*nodePositionsMain[i],sin(phi)*nodePositionsMain[i],phi]
    ni=mbs.AddNode(Rigid2D(referenceCoordinates = p1Ref, initialCoordinates = [0,0,0], 
                        initialVelocities= [0,0,0],visualization=VRigid2D(drawSize=0.3,color=color4dodgerblue)))
    nodeListMain += [ni]
    pRefListMain += [p1Ref]
mDistrMain = mMain/LMain

## create elements:
for i in range(nElementsMain):
    lElemMain = nodePositionsMain[i+1]-nodePositionsMain[i]
    oBeam = mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers = [nodeListMain[i],nodeListMain[i+1]], 
                                                            physicsLength=lElemMain,
                                                            physicsMassPerLength=mDistrMain,
                                                            physicsCrossSectionInertia=rhoIMain,
                                                            physicsBendingStiffness=EIMain,
                                                            physicsBendingDamping=EIMain*Dc_ratio_Main,
                                                            physicsAxialStiffness=EAMain,
                                                            physicsShearStiffness=GAMain,
                                                            visualization=VObjectBeamGeometricallyExact2D(drawHeight = hMain, color=color4lightgrey)
                                                ))
    elementListMain+=[oBeam]

## add constraints for main hinge
n0 = nodeListMain[0]
mC0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=0))
mC1 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=1))
mC2 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=2))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC0]))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC1]))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mNCground, mC2]))

## assemble system and check some quantities
mbs.Assemble()

dof=mbs.ComputeSystemDegreeOfFreedom()
print('dof',dof)

eigenFreq,eigenVec=mbs.ComputeODE2Eigenvalues(
                                            #    useSparseSolver=True, 
                        # numberOfEigenvalues= nModes, 
                        #constrainedCoordinates=constrainedCoordinates,
                        convert2Frequencies= True)

#compute complex eigenvalues: THIS DOES NOT WORK YET
[eigenvaluesComplex, evComplex] = mbs.ComputeODE2Eigenvalues(computeComplexEigenvalues=True,
                                                useAbsoluteValues=False)

nModes = 5

# Direct eigenfrequency calculation
print(np.round(eigenFreq[:nModes],3))

# Eigenfrequencies from Complex eigenvalues
evComplex = np.transpose(evComplex)
eigenvaluesComplex_short = np.delete(eigenvaluesComplex,np.imag(eigenvaluesComplex)<=0)
evComplex_short = np.delete(evComplex,np.imag(eigenvaluesComplex)<=0,0)
eigenFreq_undamped, unique_indices = np.unique(np.sqrt(eigenvaluesComplex_short.imag**2+eigenvaluesComplex_short.real**2)/ (2*np.pi),return_index=True)
evComplex_short = evComplex_short[unique_indices]
print(np.round(eigenFreq_undamped[:nModes],3))
evComplex_short = np.transpose(evComplex_short)

if True: #show modes:
    exu.StartRenderer()
    for i in range(nModes):
        iMode = i
        scaleFactor = 50
        exu.SetWriteToConsole(False)
        mbs.SolveStatic()
        mbs.WaitForUserToContinue()  
        exu.SetWriteToConsole(True)
        mbs.systemData.SetODE2Coordinates(scaleFactor*np.real(eigenVec[:,iMode]), exudyn.ConfigurationType.Visualization)
        mbs.systemData.SetTime(eigenFreq[iMode], exudyn.ConfigurationType.Visualization)
        mbs.SendRedrawSignal()
    
        mbs.WaitForUserToContinue()    
    SC.WaitForRenderEngineStopFlag()
    exu.StopRenderer() #safely close rendering window!

if True: #show modes:
    exu.StartRenderer()
    for i in range(nModes):
        iMode = i
        scaleFactor = 50/20
        exu.SetWriteToConsole(False)
        mbs.SolveStatic()
        mbs.WaitForUserToContinue()  
        exu.SetWriteToConsole(True)
        mbs.systemData.SetODE2Coordinates(scaleFactor*np.real(evComplex_short[:,iMode]), exudyn.ConfigurationType.Visualization)
        mbs.systemData.SetTime(eigenFreq_undamped[iMode], exudyn.ConfigurationType.Visualization)
        mbs.SendRedrawSignal()
    
        mbs.WaitForUserToContinue()    
    SC.WaitForRenderEngineStopFlag()
    exu.StopRenderer() #safely close rendering window! ```

[jgerstmayr]

As you mentioned in another discussion, you were using bending damping but no axial damping. Therefore, for axial modes, the eigenfrequencies are the same for real values (not considering the damping matrix) and complex values (including damping). For the bending modes, the eigenvalues are different, because they include damping.

As you know from the simple one-mass oscillator, damping influences the eigenfrequencies (slightly). This is what you see in your results.

As seen in the complexEigenvaluesTest.py, when distinguishing between eigenValues[0].imag and the abs(eigenValues[0]), you can retrieve the damped and the undamped case.