前兩天在論壇上看到有人問關於正方形四邊形單元（以下簡稱正四邊形單元）的剛度矩陣是否都是一樣的問題，起初乍一想，覺得應該不一樣。但心里想，既然有人提到這個問題，而且以前也沒有見過類似的討論，決定還是自己推導一把。
以下圖所示的正四邊形單元為例，
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))