Hi, how to form element stiffness matrix?

[dldhdz]

I have been trying to reproduce element matrix for 4-node plain strain element. I used isoparametic element, Gauss_Legend integration in Matlab.
Unfortunately, I found the element matrix could not match what I extracted from ABAQUS. Actually, the diagonal vaules looked bigger than corresponding ABAQUS values.

I'm confident on my Matlab code, the calculation follow the same routine as in any standard text book. Say,
Ke=sum(B'DB*det|J|*weight(i)*weight(j))*t.

Can anyone tell me what integration scheme used in ABAQUS for 2D 4 nodes elements? Or where the discrepancy come from?

Or how to debug the code?

Thanks in advance!

 

[xerf]

There are different formulations implemented in ABAQUS for plane strain bi-linear elements.

What type of element do you use in ABAQUS, CPE4 or CPE4R ?

You should compare to CPE4. (full integration)

 

[dldhdz]

thanks for your reply, xerf. i did compare with cpe4. just my results can not agree with abaqus element matrix, for the diagonal values, about 15 percent difference. I dont know why; cannot find a way to crack this.

 

[xerf]

1. You should check that your code is error free. I think you can find some textbook examples for Ke to check your results against.
2. What is 't' in your formula ?
3. According to ABAQUS Doc:
"For fully integrated first-order isoparametric elements (4-node elements in two dimensions and 8-node elements in three dimensions) the actual volume changes at the Gauss points are replaced by the average volume change of the element."
Therefore, you might want to try the same, i.e. to replace det(J(x_i,y_i)) in the Gauss formula with 0.25*Sum(det(J(x_i,y_i))), where as far as I remember x_i,y_i should be the coordinates in the parent domain.

 

[dldhdz]

Thanks for your suggestions, xerf.

1. compared plain stress element matrix with ABAQUS cps4; 100% agreement bearing numerical errors.
2. t is the thickness, 1 as default.
3. i tried 0.25*Sum(det(J(x_i,y_i))), seemed not working

4. compare plain strain again, the same discrepancies exist.
I checked the material matrix, the main values read as below:
D(1) = E*(1.0-v)/(1.-2v)/(1.0+v)
D(2) = D(1)*v/(1.0-v)
D(3) = G
I guess this should be OK.

so it seems that abaqus used some different integration for plain strain element cpe4.

I've got no clue now.

 

[xerf]

It's strange. So far I did not need to compare the stiffness matrices , but I used Gauss quadrature to integrate density of plastic dissipation over the entire model (made of CPE8 and CPE8R) and I got exactly ALLPD (i.e. total plastic dissipation).

 

[xerf]

For plane stress and plane strain the material matrix should be different.

The material stiffness matrix for plane strain (linear elastic isotropic material) should be 3x3 with
C11=(L+2*G) C12=L C13=0
C21=L C22=(L+2*G) C23=0
C31=0 C32=0 C33=G


Where L=lambda, G=Greek mu are Lame's constants.
L=E*Nu/((1+Nu)*(1-2*Nu)) , G=(1/2)E/(1+Nu)

Usually the material stiffness matrix for plane stress is modified to enforce zero stress components in the 3rd direction in the constitutive relation S=C*E.

It seems that you recast the matrix in a vector form.

Also, you should check the IP numbering convetion in the ABAQUS Doc for bi-linear elements, it might be different than the one you use.

 

[dldhdz]

Thanks, Xerf.

You are right, it's strange. I compared stiffness matrix, because i want to do some manipulations in global level. So far, I've got exactly cps4, cpe8 and cps8.

Sorry, I did not make myself clear, all constitutive matrices use 3x3. Your formula are essentially the same with mine and any standard text book.

But for cpe4, I tried Matlab alone and abaqus UEL, respectively, both give the same element stiffness matrices, i.e., they are still different from abaqus cpe4 element matrix. I dont know what make it so different for cpe4. Maybe, it was bad luck to compare cpe4 in the first place. According to abaqus documentation, the "selectively reduced integration" used for first-order isoparametric elements (4-node elements in two dimensions and 8-node
elements in three dimensions). As you mentioned in your first reply. But, apparently, cps4 is a first order 2d 4-node element as well. I used full integration for cps4 and get the same stiffness matrix with abaqus. I am curious about that. I will work on that bit to see if I can find anything useful.

 

[xerf]

I think they use reduced integration for volumetric terms and full integration for deviatoric terms. If you are familiar with hybrid formulation for incopressible materials, I think the approach is similar.

You might have to derive B and C by hand.

 

[dldhdz]

I tried a modified C, which gives a result close to Abaqus, not exactly, comparable. This is easy to implement.

I then tried a modified B, so-called B_bar method, which was tedious and result is not any close either. But I could easily misunderstood this formulation.