Hello everybody, I'm trying to solve a system of differential equations that describes the behavior of a mechanical system. I have used non-inertial reference frames so I have additional terms such as the Coriolis force.
An example of equation can be
j*ddx1+ct*dx1+kt*x1+g*dx2=0 with g=2*omega*m
m*ddx2+c*dx2+k*x2=0
M=[j 0
0 m];
C=[ct 0
0 c];
K=[kt 0
0 k];
G=[0 g
0 0];
with X=[x1 x2]'
M*ddX+(C+G)*dX+K*X=0
Now if I want to extract the eigen-frequencies of a linear system such M*ddx+C*dx+K*x=0 I have to calculate the impedance matrix
A=M^-1*K
extract the eigenvectors Avet, calculating the modal mass and stiffness matrixes excluding the damping C*dx.
Mp=Avet'*M*Avet
Kp=Avet'*K*Avet
Fp=Avet'*F
obtaining
[Mp]*{ddq}+[Kp]*{q}=[Fp]*sin(omega*t)
and finally I can get the system response.
But if I have coupled terms such as the Coriolis term, how can I get the eigenfrequencies? If I just calculate the eigenvalues of A=M^-1*K I neglect some terms (G).
Thank you in advance.
You are interested in the solution of quadratic eigenvalue problems.
The usual approach for small systems is to use a state-space formulation.
\dot{u} = v
M \dot{v} + (C+G) \dot{u} + K u = 0
A z = \lambda B z
where A and B denote block matrices
A = [M 0, 0 -K]
B = [0 M, M C+G]
For large systems you might use the eigenvectors of the undamped system to reduce the size of the system matrices. Then you can use a state-space form of the reduced system.
BTW, never use M^{-1} K. The result is a non-symmetric matrix.
If M is spd use a Cholesky decomposition M = LL^T
Dear Nils, first of all thanks for the help.
Your answer is clear up to M \dot{v} + (C+G) \dot{u} + K u = 0 (in which I suppose M, C, G and K are matrices (for example size 2x2) and u, du and ddu vectors (for example 2x1)).
Then you wrote
A z = \lambda B z
A and B are defines ass in the attached pdf (I suppose), but what are z and lambda?
Using for example Scilab, should I just write
[d, v] = eigs(A,B)
that, according to the guide (attached) "solves the generalized eigenvalue problem A * v = lambda * B * v and returns a diagonal matrix d containing the six largest magnitude eigenvalues on the diagonal. v is a n by six matrix whose columns are the six eigenvectors corresponding to the returned eigenvalues".
Once I have he eigenvalues and vectors how con I go back to the frequencies of the mechanical system and the frf?
Is there another way? (to get for example all the eigenvalues/vectors and not just the first six)
I'm sorry but this is not my field and I'm quite new. The method I actual use is attached.
Franco
Dear Franco,
You might use matlab or scipy (www.scipy.org) to solve the generalized eigenvalue problem.
Just feed your eigensolver with the matrices A and B, i.e. eig(A,B)
z denotes a complex eigenvector. \lambda is your complex eigenfrequency.
You might use the Ansatz z(t) = z \exp{\lambda t}
M,C,G,K are the matrices.
HTH,
Nils
Use eig instead of eigs. eig is an eigensolver for dense eigenvalue problems, whereas eigs is based on ARPACK for sparse eigenvalue problems.
so the command [d, v] = eigs(A,B) in Scilab is correct and d and v contains directly the eigengfrequency and the eigenvectors of the new linearized system, right?
-----------------------------------------------------------------------------------------------------------
Considering for example
K = [16792.833 24878.166 - 42293.515
24878.166 72652.218 - 1802.9716
- 42293.515 - 1802.9716 209971.19 ]
M = [0.0048794 0. 0.
0. 0.0000091 0.
0. 0. 0.0049968 ]
I have tried to reproduce the results obtained with the A=M^-1*K method whit which I get
[Avet,Aval]=spec(A) //eigenvalues of matrices and pencils
Avet =[- 0.0006395 - 0.9311615 - 0.1983151
- 0.9999998 0.3142683 0.0926311
0.0000462 - 0.1848610 0.9757513]
Aval = [7.977D+09 0 0
0 7.866D-09 0
0 0 43707274]
-----------------------------------------------------------------------------------------------------------
With your method, the matrices A and B are defined as
A =[
0.0048794 0. 0. 0. 0. 0.
0. 0.0000091 0. 0. 0. 0.
0. 0. 0.0049968 0. 0. 0.
0. 0. 0. - 16792.833 - 24878.166 42293.515
0. 0. 0. - 24878.166 - 72652.218 1802.9716
0. 0. 0. 42293.515 1802.9716 - 209971.19 ]
B = [0. 0. 0. 0.0048794 0. 0.
0. 0. 0. 0. 0.0000091 0.
0. 0. 0. 0. 0. 0.0049968
0.0048794 0. 0. 0. 0. 0.
0. 0.0000091 0. 0. 0. 0.
0. 0. 0.0049968 0. 0. 0. ]
right? Then when I use the command
[d,v]=eigs(A,B,3,'SM') //(3 --> number of eigenvalues, 'SM' compute the k smallest in magnitude eigenvalues (same as sigma = 0)
I get
d = 10^(-11) *[- 0.3493791 0. 0.
0. 0.3735524 0.
0. 0. 0.3744012]
v = [0.0010060 0.0000003 0.0000029
- 0.0002139 - 8.818D-08 - 0.0000006
- 0.0005360 - 9.059D-08 - 0.0000015
- 5382207.9 - 1158.5753 - 15295.532
1816502.8 391.02081 5162.2638
- 1068515.1 - 230.0088 - 3036.5804]
-----------------------------------------------------------------------------------------------------------
Eigenvalues and eigenvectors are completely different from values obtained with the A=M^-1*K method. The starting matrices are different, but in the first case I have that Aval=eigenfrequences of the mechanical system and Avet=responce of the system.
Using the linearized equations A * v = lambda * B * v, what represent the eigenvectors and eigenfrequencies (2xn instead of n?) How can I get the real eigenfrequencies and eigenvectors of the original quadratic eigenvalue problem (n and not 2xn)?
Franco
Thanks again
- Badges
- Science method