Hi, how to form element stiffness matrix?
Hi, how to form element stiffness matrix?
(OP)
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!





RE: Hi, how to form element stiffness matrix?
What type of element do you use in ABAQUS, CPE4 or CPE4R ?
You should compare to CPE4. (full integration)
RE: Hi, how to form element stiffness matrix?
RE: Hi, how to form element stiffness matrix?
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.
RE: Hi, how to form element stiffness matrix?
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.
RE: Hi, how to form element stiffness matrix?
RE: Hi, how to form element stiffness matrix?
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.
RE: Hi, how to form element stiffness matrix?
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.
RE: Hi, how to form element stiffness matrix?
You might have to derive B and C by hand.
RE: Hi, how to form element stiffness matrix?
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.