Wrong eigenvalue calculation since v2.16

[Disla]

Hello,

I have been testing with the unicorn file (Spherical Shell) and results depends on the solver version but also on the number of modes requested.?¿?

I have used Mecway v.14 as it integrates CalculiX Version 2.17.

I have removed the gravity effect to have a good symmetry and review carefully all the mesh details to my best.

Up to 6 modes I have been able to identify the parallelism between the negative eigen values and the corresponding ones in the load reversed problem, they are just ordered different. (They curiously have different lenght of significant numbers)

Mode Buckling factor (0,6 MPa)

1 1.026054522

2 1.051731493

3 1.110911414

4 1.166578217

5 1.480403879

6 1.501401857

Mode Buckling factor (-0,6 MPa)

1 -1.50755

2 -1.48203

3 -1.16658

4 -1.11091

5 -1.05173

6 -1.02605

When requesting up to 20 modes things look worst. The values provided by the solver below 1 seems useless.

Negative Buckling Factors .inp and 1 more file

[JuanP74]

Hi, I don’t see any big issue with this case either. As pressure is applied from inside the eigenvalues are negative. The displacements seem physically logical and regarding the number of modes there is a parameter for eigenvalue extraction that is dependent on the number of modes, in the manual:

First line:
• *BUCKLE
Second line:
• Number of buckling factors desired (usually 1).
• Accuracy desired (default: 0.01).
• # Lanczos vectors calculated in each iteration (default: 4 * #eigenvalues).
• Maximum # of iterations (default: 1000).
It is rarely needed to change the defaults.

the numbers I got are (changing # of Lanzos vector to 120):

    B U C K L I N G   F A C T O R   O U T P U T

 MODE NO       BUCKLING
                FACTOR

      1  -0.1554302E+01
      2  -0.1545839E+01
      3  -0.1538477E+01
      4  -0.1534916E+01
      5  -0.1532558E+01
      6  -0.1522505E+01
      7  -0.1521401E+01
      8  -0.1508375E+01
      9  -0.1508001E+01
     10  -0.1500724E+01
     11  -0.1491560E+01
     12  -0.1491560E+01
     13  -0.1480714E+01
     14  -0.1480714E+01
     15  -0.1479073E+01
     16  -0.1479070E+01
     17  -0.1166578E+01
     18  -0.1110911E+01
     19  -0.1051731E+01
     20  -0.1026054E+01

regarding the accuracy parameter I really don’t understand how it works…

[Disla]

Hello Juan and Victor and thank you very much for your comments and clarifications. They have been very helpful.

This is so far what I get and my conclusions. Sorry if the post is too large.

1-Ccx orders negative eigenvalues starting from biggest absolute value (confusing). If the eigenvalues extraction is done with the same criteria, the result will miss the lowest negative eigenvalues and large negative Eigenvalues may show on the list.

2-Ccx v2.17 from Mecway, v2.15MT and v2.16 Spools from Prepomax compilation have preserved the sign of all the eigenvalues in this particular model.

It helps to identify the load direction of each mode.

-Ccx v2.17 from Mecway match well when both configurations are compared. (Seems the right one to me). If Nastran standard is applied, I think this could be the golden test for the buckling solver. Both configurations ( positive and negative loads should provide the same results ¿isn’t it?) . We probably need to increase the accuracy to compare properly. It seems that the last modes on the list are less accurate. I have noticed some missing or shifted modes when comparing results.

• ccx v2.15MT and ccx v2.16 Spools provides a new set of modes with low values, but Mode shapes shows noise, are not symmetric and some gives varying results between several runs.

-All later solvers tested in my machine gives large negative values under internal pressure. Some of them do not respond or gives varying results between several runs.

Ccx 2.17.1 SPOOLS - (- 0.6/0.6 Mpa)

Ccx 2.18 Dynamic Website - (- 0.6/0.6 Mpa)

Ccx 2.18 Estatic Website - (- 0.6/0.6 Mpa)

Ccx 2.18 Website - (- 0.6/0.6 Mpa)

Ccx 2.19 Dynamic Website- (- 0.6/0.6 Mpa)

Ccx 2.19 Static Website- (- 0.6/0.6 Mpa)

From my point of view, it would be perfect to sort all the modes according to the Nastran standard but pointing to the user which modes belong to the requested load configuration and which to the opposite loading case (as it may be impossible and without interest to the user).

Victor, would be nice to be able to adjust the accuracy. ¿Is it possible with “modify Keyword”?

My calculation abstract is attached.

I have one question. ¿Is the gravity with inverted sign also considered into the possible set of buckling modes?.

Thanks in advance

Buckling Negative values.pdf

[JuanP74]

Hello…I think I now understand your point. At the begining I was just focused on the mathematical issue of eigenvalue extraction and didn’t pay enough attention to the engineering aspect. It is right that some loads do not make sense if reversed. A vessel with internal pressure shouldn’t be calculated with reversed pressure. In fact internal pressure is stabilizing so probably in this problem what makes sense is a static step with internal pressure and then a perturbation buckling step with own weight (gravity) which is the load that can make a spherical shell buckle (gravity s very particular since its magnitude is constant so a buckling factor below 1 will mean there is buckling and above 1 means no bukling but the factor itself is physicaly meaningless since gravity will not change (unless you move to another planet).
The case is different when loads can be reversed, for example in aircraft engineering depending on the manoeuvre you can have positive or negative torsion on the wing, leading to positive or negative shear loads in skin panel, spar/rib webs and then the buckling load is the lowest absolute value disregarding if it is positive o negative since both cases may take place. Hence the reason to have eigenvalues extracted and ordered considering the lowest absolute value.

Both configurations ( positive and negative loads should provide the same results ¿isn’t it?

this is only true if there is material isotropy as well (not only geometrical), for example composite material panels fail at different shear loads if they are applied positive or negative due to material anisotropy.

I have one question. ¿Is the gravity with inverted sign also considered into the possible set of buckling modes?.

for completness, I think I answered above but to be clear: if gravity (own weight) is the load that will make buckle your structure you have to consider it in the perturbation buckle step. If it is present but other type of mechanical loads are applied probably you should run own weight as a static step and on top of that you run the mechanical loads in a buckling perturbation step (this is typical when mechanical loads are significantly higher that own weight). The reason for this is as I said that having gravity multiplied by the buckling factor is wrong, it will be constant anyway.

[JuanP74]

I just modified my previous file to consider only shear loads (removed biaxial loading in the perforated panel) and then I noticed that CCX does not give the lowest positive and lowest negative (absolute value) eigenvalue when 2 are requested. Loads have to be reversed…this is not the expected behaviour and this is surely due to the interface to arpack, i.e. what parameters are passed to Lanczos solver for eigenvalue extraction. In any case the changes introduced in v2.17 are not helpful and do cast wrong results and should be removed. I’ll try to see how the interface with arpack is coded to try to get always the lowest eigenvalues in absolute value.
Attached the files (they require the pp01_ccx.inp file from initial post): pp01run_shear.zip - Google Drive

following picture show results:

[Disla]

Hello Juan,

I understand your point of view.

You would like to solve two problems at once and list the results in absolute value and properly ordered but this is very specific and could give wrong results to other users.

Think about the engineer that is designing beside you the undercarriage. It is also aircraft engineering.

If the user needs to consider both configurations, I will introduce two loading cases.

My example is not so specific. ¿What happens with the wind effects on the tank wall, or the hydrostatic pressure, or the wind friction force of the wing?.

As I comment, the Nastran standard is an option, but the user should know where those eigenmodes come from.

¿Maybe Calculix could consider two options when extracting the modes?. One for the strictly positive values (requested configuration ) and one more specific including the opposite configuration too.

[JuanP74]

Hello,
well, I consider that it is a single problem, if the loads can be reversed. In fact this is the standard in most (all?) structural analysis softwares (nastran, abaqus, samcef, altair optistruct, I think ansys as well, etc.).
I guess what would be best is if the user could select the eigenvalue range of interest for his problem as can be done in other softwares (maybe is not so easy to code this), see picture attached:

[JuanP74]

Well, the above is just a suggestion based on my experience.

At least structural users I guess would like to know if there is a negative eigenvaue with a lower absolute value than the first positive one.

In any case I think some additional information regarding this point in the manual is not only advisable but strongly recommended.

[Disla]

I agree with you.

I guess that -5E16<V1<V2<5E16 available in Lanczos would provide to the users the two options discussed above.

Manual says the range can be defined by the user, at least for frequency calculation .

• *FREQUENCY

Second line:

• Number of eigenfrequencies desired.

• Lower value of requested eigenfrequency range (in cycles/time; default:0 ).

• Upper value of requested eigenfrequency range (in cycles/time; default:∞).

This option don’t seems available in the *BUCKLE card

[vicmw]

Thanks for all this analyis @Disla and @JuanP74. It seems to confirm that the 2.15 code should be reinstated at a minimum.

The option to select a range of modes in *FREQUENCY still calculates the modes below the range but just deletes them from the output, so it wouldn’t help this problem to have that for *BUCKLE. I guess the most power would be an option to specify sigma, but it doesn’t seem to be there.

Everything below only relates to Mecway so normal people can ignore:

The CCX binary included with all versions of Mecway except the short-lived 13.0 include the pre-2.16 code, so as Disla indicated, it makes a difference where you obtained the CCX binary from.

@Disla, “modify keyword” in Mecway doesn’t allow you to set values on the data lines. But you can use it to disable the generated *BUCKLE card and enter you own replacement using “custom step contents”.

[JuanP74]

Hello, I’ve been doing some tests related to this and I would say old code isn’t necessary (I think at this moment but not 100% sure) BUT the shift (sigma) that should be used for buckling is not sigma=1, but something close to zero (NOT zero), for the time being it works well using sigma=0.0001 and very important tol=1.e-8 at least (ARPACK manual recommends using machine epsilon = 1.e-15!! for a PC) this should be changed always. Whith sigma=0.0001 the BE key for eigenvalue extraction gives both negative and positive eigenvalues but algebraically ordered, not by absolute value (like in Nastran or Abaqus) while LM gives only positive ones. I guess the best compromise would be to allow the user to modify these parameters in the input deck maybe keeping current values as default with the sole exception of tol (accuracy in CCX manual), that should be changed to 1.e-8 or 1.e-10 from current value of 1.e-2.

PS: when using BE key at least 2 eigenvalues have to be requested (BE=Both Ends). Anyone can test eigenvalue extraction capabilities of ARPACK directly using Scipy: Sparse eigenvalue problems with ARPACK — SciPy v1.8.1 Manual

[OliviaS]

On a related note - would it be possible to add a reminder to the manual that makes it clear that you want to *increase the # Lanczos vectors requested if you only request the lowest eigenvalues (e.g. BUCKLE \n 1) and do some convergence checks - otherwise the results are likely not going to be what you would expect!

[Alexandru]

Here is a test model that I computed with Calculix 2.15 and 2.19.

I also computed the model with the Permas solver from the INTES company located in Stuttgart.

The eigenvalues computed with Calculix 2.15 match those computed with Permas, so I assume this is the correct result.

The eigenvalues computed with Calculix 2.19 greatly deviate from 2.15.

I’m on OS Windows 10 but I cannot imagine, that the Linux OS would lead to other results.

Here is a picture of the results. Also interesting is that I request 10 Eigenvalues above 500Hz and both CalculiX versions provide only four or three.

[johanngil]

I just ran your file on Windows 10 with ccx 2.20 (PARDISO) solver which results in (vibration) Eigenfrequencies of:
10845.38 Hz
10859.61 Hz
28643.66 Hz
36556.73 Hz

which matches pretty well your Permas reference case.

[Alexandru]

@johanngil

ccx 2.20 (PARDISO)

Cool, many thanks.
How can I get my hands on ccx_2.20.exe?

[Disla]

@Alexandru

How can I get my hands on ccx_2.20.exe?

Same results on 2.19 (PARDISO).

10,845.4
10,859.6
28,643.7
36,556.7
36,963.4

In case you want to keep using ccx 2.19, just add :

*FREQUENCY,SOLVER=PARDISO

[JuanP74]

@Alexandru

The eigenvalues computed with Calculix 2.19 greatly deviate from 2.15.

The issue is related to solver Pastix, use pardiso or spooles

[xyont]

@Disla

Just add :

*FREQUENCY,SOLVER=PARDISO

it seems not mandatory, since versions 2.20 the solver automatically switch to Pardiso when dynamic libraries provided.

[Alexandru]

Thanks, I added the PARADISO opion and I get the error

*ERROR in arpack: the PARDISO library is not linked

I don’t remember where I got the 2.19 version from and internet research only shows the 2.15 version on github. Is there another possibility to get 2.19 including the external dependecies such as PARADISO?

[IdnoB]

Hi… is there a way to reach the user Calc_em ?