2-2图示悬臂板，属于平面应力问题，其网格图及单元、节点编号见图2-1，E=×1011，u=，演算其单刚阵到总刚阵的组集过程，并用MATLAB软件计算总刚阵。
图2-1答：根据图2-1所示列出单元节点列表：i j k1354225332654162（1）计算单元刚度阵单元1的刚度矩阵：[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=15,514,513,515,414,413,415,314,313,31kkkkkkkkkk，[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=14,514,513,515,414,413,415,314,313,31kkkkkkkkkk；单元2的刚度矩阵：[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=25,523,522,525,323,322,325,223,222,22kkkkkkkkkk，[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=24,523,522,525,323,322,325,223,222,22kkkkkkkkkk；节点单元单元3的刚度矩阵：[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=36,635,632,636,535,532,536,235,232,23k kk k k kk k k k ，[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=36,635,632,636,535,532,536,235.232,2300000000000000000000000000k k k k k k k k k k ； 单元4的刚度矩阵：[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=46,642,641,646,242,241,246,142,141,14k k k k k k k k k k ，[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=46,641,646,242,241.246,142,141,140000000000000000000000000000k k k k k k k k k ；总刚度矩阵：[][][][][][]432141k k k k kK ee +++=∑==[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡++++++++++++=46,636,635,642,632,641,636,535,525,515,514,523,513,532,522,515,414,413,425,315,314,323,313,322,346,236,235,225,223,242,232,222,241,246,142,141,100000000000k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k K Matlab 程序语言的编写：function Idexglobal gNode gElement gMaterial gNode=[]%gNode 同样是一个矩阵，每一行表示一个结点，第1 列是结点的x 坐标，第2 列是结点的y坐标gElement=[3 4 52 3 52 5 61 2 6 ];%gElement 是一个矩阵，每一行表示一个单元，第1 行是单元的第1 个结点号，第2 行是单元的第2个结点号。
Returnfunction k=StiffnessMatrix(ie)%计算单元刚度矩阵函数global gNode gElementk=zeros(6,6); %6x6单元刚阵E=*10^11; %材料特性u= ; %材料特性t=; %材料特性xi=gNode(gElement(ie,1),1);yi=gNode(gElement(ie,1),2);xj=gNode(gElement(ie,2),1);yj=gNode(gElement(ie,2),2);xm=gNode(gElement(ie,3),1);ym=gNode(gElement(ie,3),2); %计算节点坐标分量ai=xj*ym-xm*yj;aj=xm*yi-xi*ym;am=xi*yj-xj*yi;bi=yj-ym;bj=ym-yi;bm=yi-yj;ci=-(xj-xm);cj=-(xm-xi);cm=-(xi-xj);d=[1,xi,yi;1,xj,yj;1,xm,ym];area=det(d); %计算单元面积B=[bi 0 bj 0 bm 0 ;0 ci 0 cj 0 cm;ci bi cj bj cm bm];B=B/2/area;D=[1 u 0;u 1 0;0 0 (1-u)/2];D=D*E/(1-u^2);k=transpose(B)*D*B*t*abs(area); %计算单元刚度矩阵Returnfunction gK=AssembleStiffnessMatrix% 计算总刚阵global gElement gK iegK=zeros(12,12);for ie =1:1:4 %单元循环k=StiffnessMatrix(ie);for i=1:1:3 %节点循环for j=1:1:3 %节点循环for p=1:1:2 %自由度循环for q=1:1:2 %自由度循环m=(i-1)*2+p; %每个节点有2个自由度，i节点的第p个自由度为（i-1)*2+pn=(j-1)*2+q; %每个节点有2个自由度，i节点的第p个自由度为（i-1)*2+pM=(gElement(ie,i)-1)*2+p;N=(gElement(ie,j)-1)*2+q;gK(M,N)=gK(M,N)+k(m,n);endendendendendReturn则单元1的刚度矩阵为>> StiffnessMatrix(1)ans =+010 *0 00 00 0单元2的刚度矩阵>> StiffnessMatrix(2)ans =+010 *0 0 0 0 0 0 0 0单元3的刚度矩阵为>> StiffnessMatrix(3)ans =+010 *0 00 0 0 0单元4的刚度矩阵>> StiffnessMatrix(4)ans =+010 *0 0 0 0 0 0 0 0总刚度矩阵为ans =+011 *Columns 1 through 80 0 0 00 0 0 00 00 00 0 00 0 00 0 0 00 0 0 00 0 00 0 00 0 0 0 0 0 0 0 0 0Columns 9 through 120 00 00 0 00 0 00 0 0 0 0 02-3 在平面问题有限元分析中，（1）用到了哪些弹性力学中的基本方程答：平衡微分方程、几何方程、相容方程（形变协调方程）。
（2）力的平衡条件是如何满足的答：根据能量守恒原理，有外力所作虚功应该等于内力虚功。
也就是结构在外载荷作用下处于平衡状态则在结构上的力在任意虚功位移上所作的虚功之和等于零。
以下是用到的方程：000=+∂∂+∂∂+∂∂=+∂∂+∂∂+∂∂=+∂∂+∂∂+∂∂Z zy x Y zy x X z y x zzy zx yz y yx xzxy x στττστττσ z u x w z w y wz v y v x v y u x u zx z yz y xy x ∂∂+∂∂=∂∂=∂∂+∂∂=∂∂=∂∂+∂∂=∂∂=γεγεγε,,, （3）变形协调条件是如何满足的答：对材料进行线弹性和各向同性的假设，用弹性力学中应力-应变之间的关系得到变形协调条件。
下面是形变协调方程。
y x z y x z x z x x x z y x z y z y y z z y x z y x y x xy zxy zx yz zx y z y zx yz xy yz z y xyz xy zx xy y x ∂∂∂=⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂∂∂∂∂∂=∂∂+∂∂∂∂∂=⎪⎪⎭⎫⎝⎛∂∂+∂∂+∂∂∂∂∂∂∂=∂∂+∂∂∂∂∂=⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂∂∂∂∂∂=∂∂+∂∂εγγγγεεεγγγγεεεγγγγεε2222222222222222222,2,2,2-4 在平面三角形单元中的位移、应变、应力具有什么特征位移特征：(1)必须包含单元的刚体位移；(2)必须包含单元的常应变状态；(3)必须保证不偏惠各坐标轴；(4)必须保证单元内位移连续。
应力特征：(1)三角形单元其应力仅与单元材料和几何尺寸有关，与节点位移有关，而与单元内位置坐标无关，也即这类单元内的应力是常量。
(2)三角形单元内应力连续，但在公共边界上应力有突变，密布网格可以减少这种冲突的不合理性。
应变特征：由于简单三角形单元取线性位移模式，其应变矩阵为常数矩阵，即在这样的位移模式下，三角形单元内的应变为某一常量。
2-5 在平面三角形单元中，当尺寸逐步缩小，单元中的位移、应变、应力具有什么特征当单元尺寸逐步减小时，单元各点的应变趋于相等，这时常量应变成为主要成分，因此，位移应能反应这种常应变状态，由于应力矩阵也是常数矩阵，单元应力也是常量。
但是相邻单元一般具有不同的力，在单元的公共边上会有应力突变，随着单元尺寸的逐步减小，这种突变会急剧降低，从而不会妨碍有线单元法的简答收敛于精确解。