前两天在论坛上看到有人问关于正方形四边形单元（以下简称正四边形单元）的刚度矩阵是否都是一样的问题，起初乍一想，觉得应该不一样。但心里想，既然有人提到这个问题，而且以前也没有见过类似的讨论，决定还是自己推导一把。

以下图所示的正四边形单元为例，

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))