Discrepancy in Element Stiffness Matrix Output using Substructure Approach

1

Dear Calculix Community,

I am reaching out to seek assistance regarding an issue I have encountered while attempting to output the element stiffness matrix from Calculix using the substructure approach. For this purpose, I have utilized the *SUBSTRUCTURE GENERATE, *RETAINED NODAL DOFS, and *SUBSTRUCTURE MATRIX OUTPUT cards in my analysis.

In order to conduct initial testing of the element stiffness matrix from Calculix, I have set up a domain containing a single 4Node Quad element (elem type = CPS4, solving plane stress problem) and have successfully outputted the element stiffness matrix to a .mtx file. For verification purposes, I am comparing the obtained results with the ones generated by SolidsPy, an open-source FEA solver.

Upon comparison, I have noticed that the displacements are in agreement between Calculix and SolidsPy. However, there seems to be a discrepancy in the element stiffness matrix results. To provide you with a more comprehensive overview, I have attached the relevant files for your reference.

I am writing to seek your guidance and assistance in resolving this issue. If there are any insights or information you can provide regarding the observed discrepancy in the element stiffness matrix output, it would be immensely helpful for my analysis.

I appreciate your time and attention to this matter. Thank you in advance for your support. Please let me know if you require any additional information or files to further investigate this issue.

Looking forward to your prompt response.

Inputfile
Inputfile767Ă—752 93.9 KB

Calculix_vs_Solidspy
Calculix_vs_Solidspy1551Ă—886 220 KB

element_stiffness_matrix

Calculix_and_solidspy_solution_from_wolfram
Calculix_and_solidspy_solution_from_wolfram993Ă—868 102 KB

Best regards,
Anoop

2

Quick thinking…might be that the extra stiffness in ccx is due to stress stiffening? As you generate the stiffness matrix after calculation of the static step it might take that effect into account. Recalculate without the initial static step.

3

the stiffness matrix obtained with ccx corresponds to a plane strain analysis…I do not understand very well what’s going on but I think it has to be with the 3d expansion and how afterwards is transformed into a plane stress or plane strain solution.

1 Like
4

Try also extracting the stiffness matrix from CalculiX using *FREQUENCY, SOLVER=MATRIXSTORAGE as discussed here: To output element stiffness matrix to a result file

1 Like
5

I did it. I got the same result.

6

Used *FREQUENCY, SOLVER=MATRIXSTORAGE. getting global stiffness with same results.

7

I did it. Getting same result.

8

I got this from matrix generation procedure in Abaqus (format: row number, column number, matrix entry):

1 1  1.000000000000000e+36
1 2  3.749999999999999e+04
2 1  3.749999999999999e+04
1 3  2.884615384615019e+03
3 1  2.884615384615019e+03
1 4 -2.884615384615389e+03
4 1 -2.884615384615389e+03
1 5 -4.903846153846118e+04
5 1 -4.903846153846118e+04
1 6 -3.749999999999999e+04
6 1 -3.749999999999999e+04
1 7 -5.192307692307727e+04
7 1 -5.192307692307727e+04
1 8  2.884615384615387e+03
8 1  2.884615384615387e+03
2 2  1.000000000000000e+36
2 3  2.884615384615387e+03
3 2  2.884615384615387e+03
2 4  8.538461538461465e+04
4 2  8.538461538461465e+04
2 5 -3.749999999999999e+04
5 2 -3.749999999999999e+04
2 6 -1.015384615384608e+05
6 2 -1.015384615384608e+05
2 7 -2.884615384615387e+03
7 2 -2.884615384615387e+03
2 8 -1.869230769230776e+05
8 2 -1.869230769230776e+05
3 3  9.807692307692343e+04
3 4 -3.749999999999999e+04
4 3 -3.749999999999999e+04
3 5 -5.192307692307727e+04
5 3 -5.192307692307727e+04
3 6 -2.884615384615388e+03
6 3 -2.884615384615388e+03
3 7 -4.903846153846118e+04
7 3 -4.903846153846118e+04
3 8  3.750000000000000e+04
8 3  3.750000000000000e+04
4 4  2.030769230769238e+05
4 5  2.884615384615386e+03
5 4  2.884615384615386e+03
4 6 -1.869230769230776e+05
6 4 -1.869230769230776e+05
4 7  3.750000000000000e+04
7 4  3.750000000000000e+04
4 8 -1.015384615384608e+05
8 4 -1.015384615384608e+05
5 5  9.807692307692344e+04
5 6  3.750000000000000e+04
6 5  3.750000000000000e+04
5 7  2.884615384615019e+03
7 5  2.884615384615019e+03
5 8 -2.884615384615388e+03
8 5 -2.884615384615388e+03
6 6  2.030769230769238e+05
6 7  2.884615384615387e+03
7 6  2.884615384615387e+03
6 8  8.538461538461463e+04
8 6  8.538461538461463e+04
7 7  1.000000000000000e+36
7 8 -3.749999999999999e+04
8 7 -3.749999999999999e+04
8 8  1.000000000000000e+36
9

R u sure the obtained stiffness matrix is for Plane strain analysis?.
I am using element type = CPS4, which is Four-node plane stress element for this analysis.

10

Do you have the input inp?. I would like to try something.

11

SolidsPy_vs_Abaqus
SolidsPy_vs_Abaqus1238Ă—259 30 KB

The matrix value except the highlighted values from Abaqus matches with solidspy.

12

**
** Heading +++++++++++++++++++++++++++++++++++++++++++++++++
**
*Heading
Hash: 5oS75DSa, Date: 12/27/2023, Unit system: MM_TON_S_C
**
** Nodes +++++++++++++++++++++++++++++++++++++++++++++++++++
**
*Node
1, 0.00000000E+000, 0.00000000E+000
2, 5.00000000E+000, 0.00000000E+000
3, 5.00000000E+000, 2.00000000E+000
4, 0.00000000E+000, 2.00000000E+000
**
** Elements ++++++++++++++++++++++++++++++++++++++++++++++++
**
*Element, Type=CPS4, Elset=Shell_part-1
1, 1, 2, 3, 4
**
** Node sets +++++++++++++++++++++++++++++++++++++++++++++++
**
*Nset, Nset=Internal_Selection-1_Concentrated_Force-1
2, 3
*Nset, Nset=Internal_Selection-1_Displacement_Rotation-1
1, 4
**
** Element sets ++++++++++++++++++++++++++++++++++++++++++++
**
*Elset, Elset=Internal_Selection-1_Pshell
Shell_part-1
**
** Surfaces ++++++++++++++++++++++++++++++++++++++++++++++++
**
**
** Physical constants ++++++++++++++++++++++++++++++++++++++
**
**
** Materials +++++++++++++++++++++++++++++++++++++++++++++++
**
*Material, Name=Steel
*Elastic
210000, 0.3
**
** Sections ++++++++++++++++++++++++++++++++++++++++++++++++
**
*Solid section, Elset=Internal_Selection-1_Pshell, Material=Steel
1
**
** Pre-tension sections ++++++++++++++++++++++++++++++++++++
**
**
** Constraints +++++++++++++++++++++++++++++++++++++++++++++
**
**
** Surface interactions ++++++++++++++++++++++++++++++++++++
**
**
** Contact pairs +++++++++++++++++++++++++++++++++++++++++++
**
**
** Amplitudes ++++++++++++++++++++++++++++++++++++++++++++++
**
**
** Initial conditions ++++++++++++++++++++++++++++++++++++++
**
**
** Steps +++++++++++++++++++++++++++++++++++++++++++++++++++
**
**
** Step-1 ++++++++++++++++++++++++++++++++++++++++++++++++++
**
*Step
*Static
**
** Output frequency ++++++++++++++++++++++++++++++++++++++++
**
*Output, Frequency=1
**
** Boundary conditions +++++++++++++++++++++++++++++++++++++
**
*Boundary, op=New
** Name: Fixed
*Boundary
Internal_Selection-1_Displacement_Rotation-1, 1, 1, 0
Internal_Selection-1_Displacement_Rotation-1, 2, 2, 0
**
** Loads +++++++++++++++++++++++++++++++++++++++++++++++++++
**
*Cload, op=New
*Dload, op=New
** Name: Load
*Cload
Internal_Selection-1_Concentrated_Force-1, 1, 1000000
**
** Defined fields ++++++++++++++++++++++++++++++++++++++++++
**
**
** History outputs +++++++++++++++++++++++++++++++++++++++++
**
**
** Field outputs +++++++++++++++++++++++++++++++++++++++++++
**
*Node file
RF, U
*El file
S, E
**
** End step ++++++++++++++++++++++++++++++++++++++++++++++++
**
*End step

*NSET,NSET=N1
1, 2, 3, 4
*STEP
*SUBSTRUCTURE GENERATE
*RETAINED NODAL DOFS,SORTED=NO
N1,1,2
*SUBSTRUCTURE MATRIX OUTPUT,STIFFNESS=YES,OUTPUT FILE=USER DEFINED,FILE NAME=2D_Rectangle_prepomax_concentrated_load
*END STEP

13

Those are the fixed DOFs, apparently your matrix doesn’t account for BCs.

Actually, why are you using so large force for such a small model (2x5 mm) ?

1 Like
14

Is this is a good example for comparing stiffness matrices?.
Both programs are far from satisfying the assumptions established by the element selection.
The boundary conditions and element thickness do not satisfy Plane Stress and the extreme loads together with the Poissson ratio deliver a solution in wich the element is inverted.

screenshot.12
screenshot.121356Ă—230 23 KB

I would start with some suitable BC and load conditions.

1 Like
15

The reason for using large force is to get displacements in the order greater than 10^0

16

Then maybe displacement control would be better. Because currently, you can’t e.g. enable NLGEOM in this preload step (it won’t converge).

17

The force applied won’t have effect on stiffness calculation right?, just it is necessary to calculate displacement. Please correct me if I am wrong.

18

just change it to CPE4 and then you’ll see it…

19

At least in Abaqus, matrix generation procedure can account for preload:

Matrix generation:
…

  • includes initial stress and load stiffness effects due to preloads and initial conditions if nonlinear geometric effects are included in the analysis

I’m pretty sure that CalculiX does the same.

Unfortunately, such comparisons can be tricky due to expansion of 2D elements to 3D ones in CalculiX.

20

How about start with a simpler case like a 3D element with no constraints? That eliminates the changes caused by 3D expansion and whether constraints are included in the output matrix or not. It looks like a lot of struggle to work out all those complicating factors in one go.