前两天在论坛上看到有人问关于正方形四边形单元（以下简称正四边形单元）的刚度矩阵是否都是一样的问题，起初乍一想，觉得应该不一样。但心里想，既然有人提到这个问题，而且以前也没有见过类似的讨论，决定还是自己推导一把。
以下图所示的正四边形单元为例，
1、按照等参元方法构造插值函数，
位移模式可表示为，
2、构造应变矩阵B，其中，
以B1为例，因为
所以
针对本问题有
代入插值函数N1得到
所以，
同理得到
于是，
3、构造单元刚度矩阵，以平面应力为例，
根据
即可求得单元刚度矩阵。
取a=b=1,利用MABTLAB编程求解上述刚度矩阵为
更改a=b的值，运用MATLAB求解后仍得到相同的刚度矩阵。
为验证上述结果，在ADINA中建立上述不同单元边长的模型，得到的刚度矩阵都为
对比MATLAB与ADINA输出的刚度矩阵可以看出，有个别元素存在一定的出入，这个问题自己还没想明白，若有人知道其中缘由，还请帮忙指出，先谢过。
通过上述推导与求解可以得出，对于各向同性线弹性材料，当矩形单元的长与宽相等时，其刚度矩阵是一样的。
附MATLAB程序（参考SIMWE论坛上的程序）：
ex=2;
ey=2;
E=1;nu=0.3;h=1.;a=ex/2;b=ey/2;
D=E/(1-nu^2)*[1 nu 0;nu 1 0;0 0 (1-nu)/2];
syms s t;
N1=(1-s)*(1-t)/4;
N2=(1+s)*(1-t)/4;
N3=(1+s)*(1+t)/4;
N4=(1-s)*(1+t)/4;
B=1/（a*b）*[b*diff(N1,s) 0 b*diff(N2,s) 0 b*diff(N3,s) 0 b*diff(N4,s) 0;
0 a*diff(N1,t) 0 a*diff(N2,t) 0 a*diff(N3,t) 0 a*diff(N4,t);
a*diff(N1,t) b*diff(N1,s) a*diff(N2,t) b*diff(N2,s) a*diff(N3,t) b*diff(N3,s) a*diff(N4,t) b*diff(N4,s)];
K=a*b*double(int(int((B'*D*B*h),s,-1,1), t,-1,1))