ansys / pymapdl

Pythonic interface to MAPDL
https://mapdl.docs.pyansys.com
MIT License
423 stars 120 forks source link

problem with the frequency of ansys calculated and scipy calculated are different #288

Closed Villareally closed 3 years ago

Villareally commented 3 years ago

Hello, everyone,

I have a problem with the frequency calculated. when I use ansys calculated I can get [2871,6168,8689,11468] When I try to get the mass and stiffness matrix from a .full file with the scipy.sparse.linalg.eigsh module to calculate, Python give me [15.811 15.812 15.813 15.814]. I think they should be the same.

Here is the code I use for testing:

import pyansys
import os
import scipy
from scipy.sparse import linalg
import numpy as np

fullfile = os.path.join('G:\STIFF', 'CURVE.full')
full = pyansys.read_binary(fullfile)
dof_ref, k, m = full.load_km()   # returns upper triangle only

k += scipy.sparse.triu(k, 1).T
m += scipy.sparse.triu(m, 1).T

k += scipy.sparse.diags(np.random.random(k.shape[0])/1E20, shape=k.shape)

w, v = linalg.eigsh(k, k=20, M=m, sigma=10000)

f = np.real(w)**0.5/(2*np.pi)
print('First four natural frequencies')
for i in range(4):
    print('{:.3f} Hz'.format(f[i]))
Villareally commented 3 years ago

here my full file and result file. which get from ansys calculated. Hope someone can help me. file.zip

akaszynski commented 3 years ago

I'm running into the same issue as you in regards to the python solve. However, when solving this within MAPDL using MAPDL Math, I'm ending up with a different set of eigenvalues than you:


! read in matrices
*SMAT, stiff, D, import, full, file.full, STIFF
*SMAT, mass, D, import, full, file.full, mass

! solve
/SOLU
ANTYPE,MODAL
MODOPT,LANB,4,1
*VEC, MYVEC, D, ALLOC, 0
*eig,STIFF,MASS,,MYVEC 

Which returns:

  FREQUENCIES AT CURRENT LANCZOS CYCLE
     1  0.85280110E+01     2  0.12485649E+02     3  0.14800782E+02
     4  0.17033644E+02     5  0.21613667E+02     6  0.22647486E+02
     7  0.26119110E+02     8  0.29797630E+02     9  0.30706920E+02
    10  0.31378644E+02    11  0.33062955E+02    12  0.34077783E+02
    13  0.36827381E+02    14  0.42039170E+02

     number of steps      :      2
     eigenvalues found    :     14
     total no. eigenvalues:     14

  LANCZOS CYCLE NUMBER =   2

     new shift:  1.4948D+04     modes still needed:        0
  Allocate a [4] Vector : MYVEC

 *** ANSYS - ENGINEERING ANALYSIS SYSTEM  RELEASE Release 18.2.2   18.2.2   ***
 ANSYS Mechanical Enterprise                       
 88888888  VERSION=LINUX x64     12:22:41  OCT 15, 2020 CP=      0.752

 *** FREQUENCIES FROM BLOCK LANCZOS ITERATION ***

  MODE    FREQUENCY (HERTZ)      

   FREQUENCY RANGE REQUESTED=    4 MODES ABOVE   1.00000     HERTZ

    1     8.528010965666    
    2     12.48564886186    
    3     14.80078189537    
    4     17.03364386520    

Please note that I'm running this on ANSYS v18.2, and there might be a fundamental difference in full files between v17 and v18.2. Can you try running the above on v17?

Pinging @FredAns. Is there something we're potentially missing here?

Villareally commented 3 years ago

I'm running into the same issue as you in regards to the python solve. However, when solving this within MAPDL using MAPDL Math, I'm ending up with a different set of eigenvalues than you:


! read in matrices
*SMAT, stiff, D, import, full, file.full, STIFF
*SMAT, mass, D, import, full, file.full, mass

! solve
/SOLU
ANTYPE,MODAL
MODOPT,LANB,4,1
*VEC, MYVEC, D, ALLOC, 0
*eig,STIFF,MASS,,MYVEC 

Which returns:

  FREQUENCIES AT CURRENT LANCZOS CYCLE
     1  0.85280110E+01     2  0.12485649E+02     3  0.14800782E+02
     4  0.17033644E+02     5  0.21613667E+02     6  0.22647486E+02
     7  0.26119110E+02     8  0.29797630E+02     9  0.30706920E+02
    10  0.31378644E+02    11  0.33062955E+02    12  0.34077783E+02
    13  0.36827381E+02    14  0.42039170E+02

     number of steps      :      2
     eigenvalues found    :     14
     total no. eigenvalues:     14

  LANCZOS CYCLE NUMBER =   2

     new shift:  1.4948D+04     modes still needed:        0
  Allocate a [4] Vector : MYVEC

 *** ANSYS - ENGINEERING ANALYSIS SYSTEM  RELEASE Release 18.2.2   18.2.2   ***
 ANSYS Mechanical Enterprise                       
 88888888  VERSION=LINUX x64     12:22:41  OCT 15, 2020 CP=      0.752

 *** FREQUENCIES FROM BLOCK LANCZOS ITERATION ***

  MODE    FREQUENCY (HERTZ)      

   FREQUENCY RANGE REQUESTED=    4 MODES ABOVE   1.00000     HERTZ

    1     8.528010965666    
    2     12.48564886186    
    3     14.80078189537    
    4     17.03364386520    

Please note that I'm running this on ANSYS v18.2, and there might be a fundamental difference in full files between v17 and v18.2. Can you try running the above on v17?

Pinging @FredAns. Is there something we're potentially missing here?

Thanks for your help, I can get the frequencies same as you. I think the frequencies your write is good although it's a little different with scipy.sparse.linalg.eigsh. And maybe what the eigsh calculated is not right I don't know why. Another question I want to know, I use the

/solu
antype,1                  
bucopt,lanb,4            
mxpand,4,,,1             
solve
finish

to calculate the modal and I use Gui/General postproc/ read results / by pick,. There are four frequencies, which are what I write [2871,6168,8689,11468]. I want to know if they are not modal's frequencies?

ok,I maybe know what's this:

sqrt(2871)/(2*pi)=8.5280
sqrt(6168)/(2*pi)=12.499
sqrt(8689)/(2*pi)=14.836
sqrt(11468)/(2*pi)=17.044

So the frequencies in gui is circle frequencies.

FredAns commented 3 years ago

Hi, For what I see, this mathematical problem seems to be a bit hard to solve, and I see the Mapdl Lanczos algorithm making a lot of iterations to converge. My feeling is that scipy func gives wrong results, because of the numerical difficulties.

I've checked using an alternative Modal Algorithm for Symmetric Matrices ( MODOPT,SUBSPACE,4,1), an I get very close solutions to Lanczos, but SUBSPACE converge faster because it's more robust in this numerical situation :

FREQUENCY RANGE REQUESTED= 4 MODES ABOVE 1.00000 HERTZ

1     8.528801109020
2     12.49936641695
3     14.83630288045
4     17.04379504882
Villareally commented 3 years ago

Hi, For what I see, this mathematical problem seems to be a bit hard to solve, and I see the Mapdl Lanczos algorithm making a lot of iterations to converge. My feeling is that scipy func gives wrong results, because of the numerical difficulties.

I've checked using an alternative Modal Algorithm for Symmetric Matrices ( MODOPT,SUBSPACE,4,1), an I get very close solutions to Lanczos, but SUBSPACE converge faster because it's more robust in this numerical situation :

FREQUENCY RANGE REQUESTED= 4 MODES ABOVE 1.00000 HERTZ

1     8.528801109020
2     12.49936641695
3     14.83630288045
4     17.04379504882

Thanks for your help, it's so good to have a math genius like you to help my question. I'll select Subspace and Lanczos to solve generalized eigenvalue questions in engineering.

akaszynski commented 3 years ago

Thanks @FredAns!

Villareally commented 3 years ago

Thanks @FredAns!

I have known why the answer is [2871,6168,8689,11468] in the begining. Because

sqrt(2871)/(2*pi)=8.5280
sqrt(6168)/(2*pi)=12.499
sqrt(8689)/(2*pi)=14.836
sqrt(11468)/(2*pi)=17.044

This is similiar to what you calculated in Lanczos

    1     8.528010965666    
    2     12.48564886186    
    3     14.80078189537    
    4     17.03364386520  .  

Because what in Lanczos calculated is circular frequency(ω), and the Ansys's gui calculated is linear frequency (f). So what you calculated and ansys calculated are both right.

That's all, Thanks for your help.

FredAns commented 3 years ago

Good point :)

akaszynski commented 3 years ago

Closing this as it seems to be resolved. In future releases we should be able to send/receive matrices directly to and from Python.