2-2 图示悬臂板，属于平面应力问题，其网格图及单元、节点编号见图2-1，E=2.1×1011，u=0.28，演算其单刚阵到总刚阵的组集过程，并用MATLAB 软件计算总刚阵。

图2-1答：根据图2-1所示列出单元节点列表：i j k 1 3 5 4 2 2 5 3 3 2 6 5 4162（1）计算单元刚度阵单元1的刚度矩阵：[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=15,514,513,515,414,413,415,314,313,31k k k k k k k k k k ，[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=0000000000000000000000014,514,513,515,414,413,415,314,313,31k k k k k k k k k k ； 单元2的刚度矩阵：[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=25,523,522,525,323,322,325,223,222,22k kk k k kk k kk ，[]⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=000000000000000000000024,523,522,525,323,322,325,223,222,22k k k k k k k k k k ； 节点单元单元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=[0.0 0.01 0.5 0.01 1.0 0.01 1.0 0.0 0.5 0.00.0 0.0]%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=2.1*10^11; %材料特性u=0.28 ; %材料特性t=0.01; %材料特性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 =1.0e+010 *2.0508 0 -2.0508 0.0410 0 -0.04100 5.6966 0.0319 -5.6966 -0.0319 0-2.0508 0.0319 2.0531 -0.0729 -0.0023 0.04100.0410 -5.6966 -0.0729 5.6974 0.0319 -0.00080 -0.0319 -0.0023 0.0319 0.0023 0-0.0410 0 0.0410 -0.0008 0 0.0008单元2的刚度矩阵>> StiffnessMatrix(2)ans =1.0e+010 *2.0531 -0.0729 -2.0508 0.0319 -0.0023 0.0410-0.0729 5.6974 0.0410 -5.6966 0.0319 -0.0008 -2.0508 0.0410 2.0508 0 0 -0.04100.0319 -5.6966 0 5.6966 -0.0319 0-0.0023 0.0319 0 -0.0319 0.0023 00.0410 -0.0008 -0.0410 0 0 0.0008单元3的刚度矩阵为>> StiffnessMatrix(3)ans =1.0e+010 *0.0023 0 -0.0023 0.0319 0 -0.03190 0.0008 0.0410 -0.0008 -0.0410 0-0.0023 0.0410 2.0531 -0.0729 -2.0508 0.03190.0319 -0.0008 -0.0729 5.6974 0.0410 -5.69660 -0.0410 -2.0508 0.0410 2.0508 0-0.0319 0 0.0319 -5.6966 0 5.6966 单元4的刚度矩阵>> StiffnessMatrix(4)ans =1.0e+010 *2.0531 -0.0729 -2.0508 0.0319 -0.0023 0.0410-0.0729 5.6974 0.0410 -5.6966 0.0319 -0.0008 -2.0508 0.0410 2.0508 0 0 -0.04100.0319 -5.6966 0 5.6966 -0.0319 0-0.0023 0.0319 0 -0.0319 0.0023 00.0410 -0.0008 -0.0410 0 0 0.0008总刚度矩阵为ans =1.0e+011 *Columns 1 through 80.2053 -0.0073 -0.0002 0.0041 0 0 0 0-0.0073 0.5697 0.0032 -0.0001 0 0 0 0-0.0002 0.0032 0.4106 -0.0073 -0.0002 0.0041 0 00.0041 -0.0001 -0.0073 1.1395 0.0032 -0.0001 0 00 0 -0.0002 0.0032 0.2053 0 -0.2051 0.00410 0 0.0041 -0.0001 0 0.5697 0.0032 -0.56970 0 0 0 -0.2051 0.0032 0.2053 -0.00730 0 0 0 0.0041 -0.5697 -0.0073 0.56970 0 -0.4102 0.0073 0 -0.0073 -0.0002 0.00320 0 0.0073 -1.1393 -0.0073 0 0.0041 -0.0001-0.2051 0.0041 0 -0.0073 0 0 0 00.0032 -0.5697 -0.0073 0 0 0 0 0Columns 9 through 120 0 -0.2051 0.00320 0 0.0041 -0.5697-0.4102 0.0073 0 -0.00730.0073 -1.1393 -0.0073 00 -0.0073 0 0-0.0073 0 0 0-0.0002 0.0041 0 00.0032 -0.0001 0 00.4106 -0.0073 -0.0002 0.0041-0.0073 1.1395 0.0032 -0.0001-0.0002 0.0032 0.2053 00.0041 -0.0001 0 0.56972-3 在平面问题有限元分析中，（1）用到了哪些弹性力学中的基本方程？答：平衡微分方程、几何方程、相容方程（形变协调方程）。

（2）力的平衡条件是如何满足的？答：根据能量守恒原理，有外力所作虚功应该等于内力虚功。