I’ve done a simple comparison regarding CalculiX buckling capabilities simply regarding differences in element types, solvers settings, and even your layered example. The geometry is a cylinder subjected to torsion and the results are compared to test data (NACA TN Number 1344) and theoretical values from Roark’s book, which can be seen in the next images.
I use the following input data:
For the element type analysis, I run the default settings for 5 eigenvalues and PastiX solver, so I get the following results:
For the solver settings, I compare Pastix, Pardiso, and Spooles for 5 eigenvalues, 1 eigenvalue, and 1 eigenvalue with increased accuracy:
For the shell layering, I use two, four, and eight layers with Pardiso solver and S8R element type:
The results show that linear elements (mainly S4R and S3) show a high discrepancy in the results. The solver used (Pardiso, PastiX, or Spooles) does not show any difference in the eigenvalue extraction, probably because the ARPACK module is doing the same hard work in all of them. Moreover, if one requests only 1 eigenvalue, the solvers fail to find the lower one and probably the first one it finds is the one it outputs (in this case, the 3rd harmonic). However, if you increase the accuracy sufficiently, it finds the lowest eigenvalue with little to no difference in the time spent on calculation. Last of all, the multilayer composite show very good agreement with the single-layer counterpart and the reference values.
I do think it is just a simple validation, and more scenarios are relevant to analyze each parameter’s contribution to the eigenvalue extraction. If you think in some scenarios like this, I would gladly improve this comparison. This way, we can better understand the software and its limitations.