5、利用线性方程组的数值解法，对结构的有限元方程进行求解，得到所有各结点的位 移向量。最后根据需要求解单元应力。
3.16一长方形薄板如图所示。其两端受均匀拉伸P。板长12cm,宽4cm，厚1cm。材料
E?2.0?105MPa，泊松比??0.3。均匀拉力p?5MPa。使用有限元法求解板的内应力，并和
精确解比较（提示：可利用结构对称性，并用2个三角形单元对结构进行离散）。
解：
解：结点编号 1 2 3 4 单元号 X坐标 0 12 0 12 相邻结点 Y坐标 0 0 4 4
1x1y12??1x2y2?0.0024m2平面三角形单元的面积均为1
x3 y3
?2.19780.65930?D???0.65932.19780??1011应力矩阵为：??000.7692???
单元1的应变距阵为：
??16.6667016.6667000?B(1)???0?2500025????25?16.6667016.6667250???
单元1的单元刚度矩阵为：
1 2 1 3 2 2 3 4
K(1)?1.3095?0.7143???0.7326????0.3846??0.5729????0.32970.71431.9048?0.3297?0.2584?0.3846?1.6484?0.7326?0.3846?0.5769?0.3297??0.3297?0.2564?0.3846?1.6484??0.7326000.3297?9??1000.25640.38460??00.38460.57690?0.3297001.6484??
单元2的应变距阵为：
B(2)016.66670?16.66670?0????0?2502500???2516.66670?16.6667???250?
单元2的单元刚度矩阵为：
K(2)0?0.5769?01.6484???0.5769?0.3297????0.3846?1.6484?00.3297?0??0.38460.3846???0.3297?1.64840.32970?1.30950.7143?0.7326?0.3846?9??100.71431.9048?0.3297?0.2564???0.7326?0.32970.73260??0.3846?0.256400.2564??
?0.5769?0.38460总刚度矩阵为：
?1.3095?0.7143???0.7326??0.3846K????0.5769???0.3297?0?0??位移分量为：载荷列阵为：因为
0.71431.9048?0.3297?0.2564?0.3846?1.648400?0.7326?0.32971.3095000.7143?0.5769?0.3846?0.3846?0.256401.90480.71430?0.3297?1.6484?0.5769?0.384600.71431.30950?0.7326?0.3297?0.3297?1.64840.7143001.9048?0.3846?0.256400?0.5769?0.3297?0.7326?0.38461.30950.71430??0??0.3846???1.6484??109?0.3297???0.2564?0.7143??1.9048???8?1??0,v1,u2,v2,0,v3,u4,v4?
R8?1??Fx1,0,1000,0,Fx3,0,1000,0?
R8?1?K8?8??8?1
?5????0,?0.0649,0.15,?0.0369,0,?0.0525,0.15,?0.0726?10m 8?1可以得
单元1的单元应力：单元2的单元应力：
MPa?x?5.7000MPa ?y?2.3332MPa ?xy?0.3245 MPa?x?4.9505MPa ?y??0.1648MPa ?xy??0.2583
长方形薄板内应力的精确解为：拉应力5MPa，用有限元法求解出的结果与精确解大致相等。
3.17 验证三角形单元的位移差值函数满足
Ni?xj,yj???ij及
Ni?Nj?Nm?11x12A?1x2。
y1y2y3，
解：平面三角形形函数为：
Ni?12A?ai?bix?ciy?，其中，
1x3a1,b1,c1;a2,b2,c2;a3,b3,c3,分别是行列式2A中的第一行，第二行和第三行各元素的代数余子式。行列
式中，任一行的元素与其相应的代数余子式的乘积之和等于行列式的值，而任一行的元素与其它行对应元素的代数余子式乘积之和为零，故有： 当
N1?x1,y1???12???a1?b1x1?c1y1??1，同时有，
??N2?x2,y2???12???a2?b2x2?c2y2??0???N3?x3,y3???12???a3?b3x3?c3y3??0
同理也有：
N2?x1,y1??0,N2?x2,y2??0,N2?x3,y3??0N3?x1,y1??0,N3?x2,y2??0,N3?x3,y3??0,
。
?N1?x,y??N2?x,y??N3?x,y??1，即
Ni?Nj?Nm?13.18 推导如图所示的9节点矩形单元的形函数。
解:三维杆单元的形状函数，
?1??3??1??2?2??3?N1x???x1?x2,N?,N?5x2x??x1?x3??x2?x1??x2?x3??x3?x1??x3?x2?x?xx?xx?xx?xx?xx?x ①
在局部坐标系中令节点1,5,2所对应的的形函数：
?N1x?????0??a/2???0?a????x??a/2???x?a?2x?a??x?a?a2x1?0,x2?a/2,x3?a带入①式得到节点1,5,2仅在x方向上
②
同理可得：
N5x??x?0??x?a????a/2??0?????a/2??a???x?x?a??a2/4
由
y1,y2,y3，即节点2,6,3，可得到沿着全局坐标系y轴的形状函数(通过变量轮换)，节点1的形函数
?2x?a??x?a??2y?b??y?b?a2b2即x，y方向的乘积：
N1?N1xN1y?由此可得：
?9??7?N5y???y5?y9???y5?y7?y?yy?y2y?b??y?b?b2
a2b2?N5?N5xN5y??同理可整理得：
4x?2x?a??2b?y??y?b?
N2?N2xN2y?N3?ab22x?2x?a??2y?b??y?b?a2b2
?2x?a??x?a??2y?b?ab22x?2x?a?y?2y?b?，
N4?y，
N6??，
4xy?2x?a??y?b?a2b216xy?x?a??y?b?a2b2， ，
N7??4xy?x?a??2y?b?ab22，
N8??4y?2x?a??x?a??y?b?ab22N9?3.19 如图所示为一个桁架单元，端点力为[U1，U2]，端点位移为[u1，u2]，设内部任一点的轴向位
移u是坐标x的线性函数：
u?a1?a2x
推导其形函数矩阵N。
?a?u?x???1x??1?u?a1?a2x，写成向量形式为?a2?，设两个节解：轴向位移u是坐标x的线性函数，
点的坐标为
xi,xj，代入向量形式的位移函数解
a1,a2得：
?a1??1xi??ui???????u?1xaj??j??2??
?1xi??ui?u?x???1x????u?1xj??j??则由位移函数可得形函数为：
?1?1
4.1 答：轴对称三角形环单元不是常应变单元，如果弹性体的几何形状、约束条件及载荷都对称于
某一轴，则所有的位移应变及应力也是对称于此轴，这样问题称为轴对称。轴对称三角形环单
?1xi?N??1x????1xj???111xi1xj??xj?xx?xi?????Ni?xj?xNj????x?x?ji?x?xi??xj?xi??
元与平面常应变单元是不同的，轴对称三角形环单元的应变不是常数矩阵，其应变矩阵B=[Bi
bi0acz1fi0Bj Bm],其中Bi=[]，fi?i?bi?i（i,j,m）。应变分量?r，?z，?rz都是常量，
rr2?0cicibi但环向应变??不是常量，它与fi，fj，fm中的r和z有关。
4.2 答：轴对称问题中，刚度自由度：环向位移，径向位移，轴向位移。以三角环单元平均半径、
平均高度进行计算的单元刚度矩阵，配合以精确积分所得的等效结点载荷矩阵，计算的结果还是不错的！
4.3 轴对称问题的两个单元a和b，设材料的弹性模量为E，泊松比为μ = 0.15，试手算这两个单元的
刚度矩阵。
解：对于a单元，由题可知：
ria?lzia?0rja?0zja?lrma?0zma?0zia1l11ra?(ria?rja?rma)?l3311za?(zia?zja?zma)?l33单元a的截面面积为
?a?11rja21rma1ria11zja?10l?l222zma1000aia?rjarmazjazma?rjazma?rmazja?0aja?rmariazmazia?rmazia?riazma?0bia??1zja1zma?zja?zma?lbja??1zma1zia?zma?zia?0cia?1rja1rmariarja1ria1rja??(rja?rma)?02cja?1rma1ria1zia1zja??(rma?ria)?lama?ziazja?riazja?rjazia?lbma???zia?zja??lcma???(ria?rja)??l