现代设计方法43三角形三节点平面单元课件.ppt

1、三角形三节点平面单元三角形三节点平面单元结构离散化结构离散化单元分析单元分析整体分析整体分析有限元分析的基本步骤：1 结构离散化 例：图示一悬臂梁，梁的厚度为t，设泊松比m=1/3，弹性模量为E，试用三节点三角形单元进行离散。l2 m1 mAP/2P/23P/2P/2412yiijjmmx2 单元分析l单元分析的主要内容：由节点位移求内部任一点的位移，由节点位移求单元应变、应力和节点力。l单元分析的步骤：节点位移单元内部各点位移单元应变单元应力节点力(1)(2)(3)(4)单元分析(1)单元位移模式l有限元法中，通常采用位移法进行计算，即取节点的位移分量为基本未知量，单元中的位移、应变、应力等

2、物理量，都和基本未知量相关联。l 节点i的位移分量可写成 d di i=u ui i v vi i T Tl单元节点位移向量d e可写成 d d e e =d di i d dj j d dm m T T =u ui i v vi i u uj j v vi i u um m v vm m T T三角形三节点单元xymjiviuiujvjvmumuv2.1 由节点位移求单元内部任一点位移l六个位移分量需六个待定参数a1，a2，a6 l设单元内任一点位移分量是坐标(x,y)的线性函数：l (4.22)l写成矩阵的形式为：l (4.23)yxyxvyxyxu654321),(),(123456(,

3、)10 0 0(,)(,)(,)0 0 0 1u x yx yx ym x yv x yx y l设节点i，j，m坐标分别是xi，yi；xj，yj；xm，ym。把三个节点的坐标及其水平位移代入式(4.22)中得：l (a)l解得：l (b)l其中，(2)由单元节点位移d e求位移参数a 123123123iiijjjmmmuxyuxyuxy11iiijjjmmmuxyuxyAuxy21111iijjmmuyuyAuy31111iijjmmxuxuAxummjjiiyxyxyxA111123111iiijjjmmmuxyuxyuxyl对v同理可列出a4、a5、a6的方程：vi=a4+a5 xi+

4、a6 yi i,j,ml解得：a4、a5、a6l为了书写的简便与规格化，引入记号lai,bi,ci分别为行列式A中各行各个元素的代数余子式1111jjijmmjmmjijmmjimjmxyax yx yxyybyyyxcxxx mmjjiiyxyxyxA111i j m1234561()21()21()21()21()21()2iijjmmiijjmmiijjmmiijjmmiijjmmiijjmmaua ua ubub ub ucuc uc uava va vbvb vb vcvc vc v解出a1 a6结果：la=Ld e为三角形单元面积。11112211()2iijjmmjmmiijmj

5、imjixyAxyxyx yx yx yx yx yx y 00000000010002000000ijmijmijmijmijmijmaaabbbcccLaaabbbccc将a写成矩阵形式，有a=Ld el由单元节点位移d e求单元内部任一点位移d(x,y)ld(x,y)=m(x,y)a=m(x,y)Ld e00000000010001(,)00000012000000ijmijmijmijmijmijmaaabbbcccxym x yLaaaxybbbccc(,)0(,)0(,)00(,)0(,)0(,)ijmijmN x yN x yN x yN x yN x yN x y(,)Nxy)

6、(21),(ycxbayxNiiiii j mlNi,Nj,Nm是坐标的连续函数，反映单元内位移的分布状态，称为位移的形状函数，简称形函数(shape function)。矩阵N称为形函数矩阵(shape function matrix)。l形函数物理意义：Ni(x,y)，节点i单位位移，其它节点位移分量为0，单元内部产生位移分布形状mmjjiimmjjiivyxNvyxNvyxNyxvuyxNuyxNuyxNyxu),(),(),(),(),(),(),(),(,)(,)(,)(,)eu x yx yN x yv x yl在单元任一点上，三个形函数之和等于1。l形函数Ni在i点的函数值为1，

7、在j点及m点的函数值为零。l三角形单元i,j,m在ij边上的形函数与第三个顶点的坐标无关。形函数的性质：例：求图示单元和单元的形函数矩阵3P/2P/2412yiijjmmx （a）(b)(0,a)(0,0)(b,0)yijmxab(b,a)(0,a)yijmx(b,0)分别如上图所示：)(21),(ycxbayxNiiiii j mabyxyxyxyxyxyxjmmjmiimijji21)(21 单元如图所示。设a=1m,b=2m.l(或直接由图形可知其面积)求系数ai,aj,am,bi,bj,bm,ci,cj,cm0000ijmmjjmiimmijjiijmjmjmijimjjimmjiax

8、 yx yax yx yax yx yabbyyabyybyyacxxcxxbcxxb (0,a)(0,0)(b,0)yijmxab)(21),(ycxbayxNiiiii j mbxycxbayxNiiii)(21),(ayycxbayxNjjjj)(21),(aybxycxbayxNmmmm1)(21),(求形函数矩阵l代入相关常数：00000000100001ijmijmNNNNNNNxyxybabaxxyaaybb yxyxyxyxN210020021002将a=1,b=2代入得：bxxcbxxcxxcayybyybayybabyxyxaabyxyxaabyxyxaijmmijjmij

9、imimjmjiijjimmiimjjmmii00(b,a)(0,a)yijmx(b,0)求常数l单元如图所示=ab/2bxycxbayxNiiii1)(21),(ayycxbayxNjjjj1)(21),(aybxycxbayxNmmmm1)(21),(求形函数矩阵)1(010100)1(0101aybxaybxaybxaybxNyxyxyxyxN2101021002101021将a=1,b=2m代入上式得：作业：l已知三角形三节点单元坐标如图示，设单元中一点A的坐标(0.5,0.2)，已知三角形三节点单元i节点位移(2.0,1.0)，j节点位移(2.1,1.1)，m节点位移(2.15,1.

10、05)，1)写出单元的位移函数；2)求A点的位移分量。i(1,0)j(1,1)m(0,0)xyA2.2 由节点位移求单元的应变 mmjjiimmjjiivyxNvyxNvyxNyxvuyxNuyxNuyxNyxu),(),(),(),(),(),(),(),(xvyuyvxuxyyxmmjjiimmjjiimmjjiimmjjiivxyxNvxyxNvxyxNuyyxNuyyxNuyyxNvyyxNvyyxNvyyxNuxyxNuxyxNuxyxN),(),(),(),(),(),(),(),(),(),(),(),(几何方程)(21),(ycxbayxNiiiii j miiiicyNbxN

11、21,21i j m00010002iixijmjyijmjxyiijjmmmmuvbbbucccvcbcbcbuv(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)jimijmxjimyijmxyjimijmjimijmNx yN x yNx yuuuxxxNx yN x yNx yvvvyyyNx yN x yNx yuuuyyyNx yN x yNx yvvvxxx简记为e=Bd eB可写成分块的形式：B=Bi Bj Bm B称为应变矩阵，它的元素都只与单元的几何性质有关的常量。这种单元称为平面问题的常应变三角形单元。单元应变与单元节点位移关系iiiiibccbB0

12、021(i,j,m)l物理方程 s=De 而 e=Bde ls =DBd e =S d e (求应力的表达式)l记 S=DB S应力矩阵:S=Si Sj Sm2.3 由单元节点位移求单元的应力iiiiiiibccbcbES2121)1(221)平面应力问题：l代入D及Bl得:S=Si Sj Sm.l对于平面应力：l (i,j,m)l将上式中以 E/(1-m2)代 E，以 m/(1-m)代 m 则子矩阵：l (i,j,m)2)对于平面应变问题：iiiiiiibccbbbES)1(221)1(22111)21)(1(2)1(l考虑单元平衡，节点力是作用在单元上的外力，与单元应力平衡。有限元法中以虚

13、功方程代替平衡方程。l节点力列阵及单元内应力列阵：单元节点力是指单元和节点相连接的内力；考虑节点平衡，节点力为外力，与节点外载荷平衡；sysxtxymmjjiieVUVUVUFxyyxtsss2.4 由节点位移求单元节点力(求单元刚度矩阵)节点虚位移列阵及虚应变：*xyyxt*mu*iu*ju*mv*iv*jv*mmjjiievuvuvu*xyyxl令实际受力状态在虚位移状态上做虚功，虚功方程：l由 e=Bd e 知 e*=Bd*el则 e*T=(d*e)TBTl由于d*e中的元素为常量，提至前积分号前，故：l(对于三角形三节点单元，B和s 为常量，单元厚度t也是常量；Ayxdd为三角形单元面

14、积，用表示）yxtFATeTedd)(*syxtBFATTeeTedd)()(*syxtBFATeddsmmmjmijmjjjiimijiiekkkkkkkkkk单元刚度矩阵tBFTes eTFBDt物理方程物理方程 s s =D De e 简记为 Fe=ked eATeyxtBDBkddke=BTDBt 或 =BTSt TeBD Bt几何方程几何方程 e e =B Bd d e esrsrsrsrsrsrsrsrsTrrsbbcccbbcbccbccbbEttBDBk21212121)1(42对于平面应力问题：(r=i,j,m;s=i,j,m).21E1将上式中的E 换成换成,srsrsrs

15、rsrsrsrsrrsbbcccbbcbccbccbbtEk)1(221)1(2211)1(2211)1(221)21)(1(4)1(对于平面应变问题：求例4.2（p84）单元的单元刚度矩阵baabaabbbcbcbccccbbbBmmjjiimjimji11011010100001000100000021(0,a)(0,0)(b,0)yijmxabiiiiibccbB0021(i,j,m)ijmmjijmimjax yx ybyycxxl解：（1）求矩阵B2100010112EDbaababbbababEBDS21210212101100100112（2）求矩阵S21001001100001

16、1 001210011110011222011210ebbbabaabtEkbbbaabaabbaaba（3）求矩阵kelke=BTDBt=BTSt2222222222222222211001110022211100222112(1)00111111()()2222111111()()22221212bbabbabababEtabaaabaabababaabbaabababbaabbabbaababab)811()412(141812)2141(22141411002214108104112称对Etkl代入a=1,b=2m 得：l可算出，当a=b时单元刚度矩阵与尺寸a,b无关在单元节点力列阵F

17、e、单元应力列阵s e、单元应变列阵e e和单元节点位移列阵d e的四个列阵之间，存在五个转换关系，可得五个转换矩阵。eF16e13se13e16tBT3663BDStBDBkTe6633D63B平衡条件物理关系几何关系单元刚度矩阵的特性：(1)单元刚度矩阵的物理意义：单元刚度矩阵ke表示了单元抵抗变形的能力，即表示了节点位移d e与节点力Fe之间的关系。kij表示节点j发生单位位移时，其它节点位移分量均为零时，在节点i上产生的节点力。(2)分块性质:单元刚度矩阵可以分块运算。eemmjjiimmjjiiekvuvuvukkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkV

18、UVUVUF666564636261565554535251464544434241363534333231262524232221161514131211eemjimmmjmijmjjjiimijiimjiekkkkkkkkkkFFFF(4.51)按节点进行分块，则单元刚度矩阵的分块形式可写为：(4.52)（3）对称性 单元刚度矩阵是一个对称矩阵 式(4.51)中：knl=kln (n=16;l=16)分块形式中：krs=ksrT (r=i,j,m;s=i,j,m)（4）奇异性 单元刚度矩阵是一个奇异矩阵|ke|=0，表明其逆矩阵不存在。即，如果给定了单元节点位移可以得出唯一的节点力:Fe=

19、ked e。反之，如果给出节点力却无法求出确定的节点位移。因为这时的单元未考虑所受约束时，可能存在不引起单元应力和节点力的刚性位移，这部分刚体位移由节点力是无法唯一确定的。三节点三角形单元每行元素之和为零三节点三角形单元每行元素之和为零例：证明图示单元刚度矩阵:kI=kIIIl证明：由于单元刚度矩阵l ke=BTDBtl可知：当两个三角形单元几何尺寸相同时，t值和单元面积值均相同；当两个单元的材料性质相同时，弹性矩阵D也时相同的。故ke是否相同，取决于矩阵B是否相同。l不难验证,I、III单元的上述br,cr(i,j,m)值均相等。l结论：l两个单元刚度矩阵ke相等的条件为：只要两单元的形状、

20、大小，方向和单元弹性常数均相同，并且编号的方式也相同(如按逆时针方向编号为 i,j,m，直角顶点编号为m)，则两个单元的刚度矩阵时相等的。mmjjiimjimjibcbcbccccbbbB000000212.5 单元载荷的移置.（离散时每个单元受载作用于节点上）(a)原则：将单元载荷向节点处移置，按照虚功等效的原则进行。对于变形体(包括弹性体)，虚功等效是指原载荷与节点载荷载在任何虚位移上做的虚功相等。当位移模式确定后，载荷移置(或分解)其结果是唯一的。虚功等效包含了刚体体系的静力等效，当虚位移为刚体位移时，虚功等效即为静力等效，静力等效是虚功等效的特例。离散(b)载荷移置公式(1)集中力 设

21、单元i,j,m中任一点M(x,y)处受有集中力P=Px PyT，移置到该单元各节点处载荷列阵为Re=Xi Yi Xj Yj Xm YmT (1)集中力P=Px PyT(1)集中力l假设该单元发生一微小虚位移，M点相应的虚位移为 f*，该单元各节点处相应虚位移为d*，由静力等效原理，载荷与节点等效载荷在虚位移上所作虚功相等：l (d*e)TRe=f*TPl 将 f*=N d*e 代入上式：l有(d*e)TRe=f*TP=(d*e)TNTPl则 Re=NTPlRe=Xi Yi Xj Yj Xm YmT =NiPx NiPy NjPx NjPy NmPx NmPyTl设单元i,j,m的一边受有分布的

22、面力P=Px Py将微元面积tds上的面力合力Ptds当作集中载荷，可得面力的移置公式:(2)面力sTestPNRd(2)面力 yxPPP(3)体力l设单元i,j,m受有分布体力G=Px PyT将微分体积tdxdy上的体力合力Gtdxdy当作集中载荷dG同理可得：d dTeRNG txy(3)分布体力G=X Y T(4)三节点三角形单元上同时有体力、面力和集中力等(a)集中力P=Px PyT(b)分布体力G=X Y T(c)面力 yxPPP(d)单元虚位移d*=u*i v*i u*j v*j u*m v*mT f*=u*v*l单元各节点处载荷列阵为 Re=Xi Yi Xj Yj Xm YmTl

23、应用虚功等效原则：l将 f*=Nd*e 代入上式：l虚位移是任意的，从而矩阵d*eT 也是任意的，故：stPfyxtGfPfRTATeTeTeddd)()(*stPNyxtGNPNRTTeATTTeTTeeTed)(dd)()()(*stPNyxtGNPNRTATTTedddl例：设三角形单元i,j,m的ij边作用有线形分布的法向载荷，i和j两点的压力集度分别为qi和qj，试用公式l求其等效节点载荷。单元厚度为t，节点坐标如图示。sTedsPNRiiqjq),0(jy),(mmyxm)0,0(jxy11(,)()()22iiiimmN x yab xc yy xx y)(21)(21),(yx

24、xyyyxycxbayxNmimimjjjjxyycxbayxNimmmm21)(21),(mjimjiNNNNNNyxN000000),(l解：计算常数ai=xj ymxm yj=0；bi=yj ym=ym；ci=xmxj=xm；aj=xm yi xi ym=xm yj；bj=ym yi=ym yi；cj=xi xm=xm；am=xi yj xj yi=0；bm=yi yj=yi；cm=xj xi=0.l计算形函数l由l得：l在边界jm和mi上的面力为零，故上式积分中后两项为0，在ij边上的面力分量可表示为：l代入上式中得：stPNstPNstPNRmijmijLTLTLTeddd stPN

25、RTed0_jjiiqqqyyYXP2 22 22 2 0iijjeijLijmN IyqqqyRN ItdsN I计算等效节点载荷 000)(0000000000dyXNdyXNdyXNtytdYXNNNNNNRiiiiymyjyiymmjjiie=0 TjijiieqqqqtyR 000)2(0)2(6积分沿逆时针方向，有ds=-dy 所以l在ij边上 x=0，代入 Ni Nj Nm中：引入支承条件解方程求位移求单元应力3 整体分析建立整体刚度矩阵l任务：建立整个结构的总刚度方程；引入边界条件解方程。3.1 建 立 整 体 刚 度 矩 阵3 3P/2P/24 41 12 2yiijjmmI

26、 IIIIIx单元节点受力分析作用于节点的集中力单元作用在节点上的力作用在节点上等效载荷2 2X2Y2Y2XxP22V2VyP22U2Ul由节点的平衡条件，有平衡方程：l 环绕节点2的所有单元求和：lF2=U2e V2eT 单元e作用在2节点上的节点力lR2=X2 Y2T 绕节点2的各单元作用在2节点上的等效节点载荷及直接作用于2的集中力之和。222222220,0,eexxeeeeyyeeFUXPXFVYPY eRF22e mjimmmjmijmjjjiimijiimjikkkkkkkkkFFF每一个单元可用节点位移表示节点力，采用分块P89 4.46l考虑节点2，对应I单元的局部编号为i即

27、I单元i节点，故I单元i节点上的节点力：l代入l得,iiiiijjimminnn i j mFkkkk 22RFe2,2Rkemjinnn RK结构的整体刚度矩阵结构的节点位移列阵结构的节点载荷列阵l每个节点均可得到类似的方程，每个节点可写出两个平衡方程，按节点的序号排列，用矩阵表示：l结构的节点载荷列阵 R=X1 Y1 X2 Y2 Tl故l其中l l=14l l,k=14 4321444342413433323124232221141312114321KKKKKKKKKKKKKKKKRRRR lllYXR lllvuklklklkllkKKKKK2,212,22,1212,12整体刚度矩阵的

28、集成规则：(1)先求出每个单元的刚度矩阵ke；(2)将ke的每个方块kije进行换号，换成对应的整体编号;(3)将换号后的子块填入整体刚度矩阵上对应的位置。(4)若在同一位置上有几个单元的相应子块填入同一位置，则进行叠加。图4.24 刚度矩阵 单元号局部号ijmijm整体号241423(书上错误)单元刚度矩阵ke的子块换号k22k24k21k44k42k43k42k44k41k24k22k23k12k14k11k34k32k333P/2P/2412yiijjmmx111412414442212422kkkkkkkkkk333234232224434244kkkkkkkkkk1112142122

29、222324243233344142424344440 0kkkkkkkkkkkkkkkkkkkki j m2 4 1i j m4 2 3i 2j 4m 1i 4j 2m 3l解:(1)求各单元刚度矩阵(P90已求)P90 4-50l求得单元的刚度矩阵为：lkI=kII (I II单元刚度矩阵相等)例：4.2。8941141812414232214141121002412102141081410418102412004112111412414442212422EtkkkkkkkkkkI333234232224434244kkkkkkkkkk111412414442212422kkkkkkkkk

30、k333234232224434244kkkkkkkkkk44444342424134333224242322222114121100kkkkkkkkkkkkkkkkkkk(2)整体刚度矩阵kl换号的单元刚度矩阵：(3)将、的单元刚度矩阵填入(P101 式4.62)111213142122238 8313233414224343444kkkkkkkkkkkkkkkkkl如果I II都有就相加便于计算机运算24104104k 341144128k 总体刚度矩阵的特性：l(1)对称性 总体刚度矩阵是一个对称矩阵。因单元刚度矩阵升阶后对称性不变，由之合成的总体刚度矩阵自然是对称矩阵。l(2)奇异性

31、总体刚度矩阵行列式的值|K|=0l(3)稀疏性 总体刚度矩阵是一个稀疏矩阵；即矩阵中的绝大多数元素为0，非0元素只占元素总数的很小的一部分。因为只有当节点ln相关时Kln才不是0，P101图4.26。l(4)带状分布规律 分布在以主对角线为中心的带状区域内。总体刚度矩阵的特性：3.2 节点载荷列阵 集成法求整体结构的节点载荷列阵步骤：(1)求出每个单元的等效节点载荷；(2)单元等效节点子向量(k=i,j,m)换号，换成对应整体编号；(3)换号后等效节点载荷子向量送到整体节点载荷列阵对应位置；(4)同一位置上若有多个单元的等效节点载荷子向量，叠加；(5)节点上有直接作用的节点载荷，按整体节点号进

32、行叠加；(6)若节点k具有水平和垂直方向的支承，支承反力为未知量，可暂设为Rke=Qxk QykT。图4.24中：lR=Qx1 Qy1 X2 Y2 X3 Y3 Qx4 Qy4 Tl =Qx1 Qy1 0 -P/2 0 -P/2 Qx4 Qy4 T3P/2P/2412yiijjmmx 4321444342413433323124232221141312114321KKKKKKKKKKKKKKKKRRRR3.3 引入支承条件2Y21333211226665646362615655545352514645444342413635343332312625242322211615141312110000

33、0yxyxQQYQQvukkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk(4.66)例：基本未知量为节点2的位移，只需在方程中抽出第三、四行即可：矩阵形式 (4.67)22442432342330Yvkukvkuk222444334330Yvukkkk0000010000001000000000000000010000001233221144433433Yvuvuvukkkk为了便于编程，修改后的矩阵仍保持原矩阵K的阶数、排列次序及矩阵对称性。式(4.67)扩大成：l等效：u1=v1=u3=v3=0 与支承条件及前述矩阵一致l多个节点：节点n的水平位移un=0，则Kd=

34、R改为：K中的2n-1行与列中主对角元素为1，其他为0 载荷向量R中2n-1元素置0l若vn=0 则对应K中2n行和列作上述修改。002020001000000001000000000000000000000000000010000000014433221166656463565554534645444336353433PPvuvuvuvukkkkkkkkkkkkkkkkl仍以课本P98 图4.24lu1=v1=u4=v4=03P/2P/2412yiijjmmx3.4 解方程 求节点位移 Kd=R222333 2 110042490182213 2110224411912448upvEtuvp

35、 l 将整体刚度矩阵K代入支承条件l 图4.243P/2P/2412yiijjmmx202013412272413073233322PPvuvuEt称对设m=1/3得：得：99.888.142.850.13322EtPvuvu于是节点的位移向量为d=u1 v1 u2 v2 u3 v3 u4 v4T =P/Et0 0 -1.50 -8.42 1.88 -8.99 0 0T单元1：当a=1,b=2,m=1/3时，由例4.33.5 求单元应力s e=Sed e613103161016110061312131002189IES580.1281.0844.0000042.850.161310316101

36、6110061312131002189tpEEtPs同理可求单元418.0289.0844.023.254.150.4163tPps单元Itxytxysxsy0.281P/t1.58P/t1.580P/t0.844P/t单元IItxytxysxsy0.289P/t0.418P/t 0.418P/t0.844P/t121331161121613112761316141161100616131031610121610611210614161004189Etk3.6求节点力及支承反力 Fe=ked e单元的节点力 a=1,b=2,m=1/33P/2P/2412yiijjmmx单元的节点力 PPPPP

37、PF070.1000.2281.0580.1789.0422.000yxFFPRPRyx07.100.211PRPRyx071.0000.244同理可得FII 也进行验算;由受力图可得 验算ijm0.289P1.580P2.000P1.070P0.789P0.422P0.004Pijm0.498P0.418P0.289P0.209P0.422P总结：l有限元求解弹性力学平面问题步骤如下：(1)整理原始数据，结构离散化，对单元和节点编号。(2)求单元刚刚度矩阵ke(3)用刚度集成法，形成结构整体刚度矩阵K(4)求节点等效载荷，写出载荷列阵 R(5)引入支承条件(6)解K d=R 求出节点位移d(

38、7)求单元应力s e=S ed e(8)求单元节点力 Fe=Ked e(9)整理结果，作节点位移图及应力图。一、matlab 基础在 matlab 的提示符 符“”下输入命令.3*4+5 ans=17 cos(30*pi/180)ans=0.8660 x=4 x=4 2/sqrt(3+x)ans=0.75594 线性三角形弹性力学平面问题的matlab程序：(德)卡坦-韩来彬清华大学出版社若不让 matlab 输出运算数据，在命令行的结尾输入分号：y=32;z=5;x=2*y-z;w=3*y+4*z w=116 matlab区分大小写x=1 x=1X=2 X=2x x=1 使用help命令可以

39、获得所有matlab命令的详细用法 help inv下面的例子显示如何输入矩阵并实现简单的矩阵运算。x=1 2 3;4 5 6;7 8 1 2 34 5 67 8 9x=y=2;0;-3 y=2 0 -3 w=x*y w=-7 -10 -13求解下面的联立方程组采用高斯消去法解方程组:A=2 -1 3 0;1 5 -2 4;2 0 3 -2;1 2 3 4 b=3;1;-2;2 b=3 1 -2 2 x=Ab.x=1.9259 -1.8148 -0.8889 1.5926A=2 -1 3 0 1 5 -2 4 2 0 3 -2 1 2 3 4采用如下方法也可求 x=inv(A)*b x=1.9

40、259 -1.8148 -8.8889 1.5926时间长，矩阵大尤其长 D=1 2 3 4 5;2 4 6 8 9;2 4 6 2 4;1 7 2 3-2;9 0 2 3 1 D=1 2 3 4 5 2 4 6 8 9 2 4 6 2 4 1 7 2 3 -2 9 0 2 3 1可以以矩阵中提取 24 行，35列作子矩时：E=D(2:4,3:5)E=6 8 9 6 2 4 2 3 -2 提取D的第3列：F=D(1:5,3)F=3 6 6 2 2 提取D的第2行作为子矩阵 G=D(2,1:5)G=2 4 6 8 9提取D中的4行 3列 H=D(4,3)H=2 绘制y=f(x)定义 x.y.再用

41、plot(x,y)绘图。x=1 2 3 4 5 6 7 8 9 10 y=x.2 plot(x,y)12345678102030405060yx)(iiyxi)(jjyxj)(mmyxm平面应变弹性问题平面应力弹性问题 二、线性三角元1.基本方程.Linear triangular element 线性形函数.系数：弹性模量E 泊松比m 厚度t单元刚度阵：k=BTDBt2=xi(yj-ym)+xj(ym-yi)+xm(yi-yj)mmjjiimjimjibcbcbccccbbbB00000021ai bi ci (P79 4.26).平面应力 2100010112ED21EE1 平面应变 显然

42、线性三角形元有6个自由度每个节点2个自由度。对一个有n个节点的结构，整体K为2n*2n Kd=F 边界条件手动赋值 采用高斯消去求解 s=Sd 求得单元应力矩2，用到的Matlab函数线性三角形元用到的5个matlab函数分别为：(1)LinearTriangleElementArear(xi，yi，xj，yj，xm,ym)该函数根据给出的i,j,m的节点坐标返回单元的面积。http:/ K。(4)LinearTriangleElementStresses(E,NU,t,xi,yi,xj,yjxm,ym,p,u)u为单位位移矢量,返回单元应力(5)LinearTriangleElementPS

43、tresses(sigma)计算单元主应力，返回3*1.sigma1,sigma2,thetaT 例：求解如图所示的受均布载荷作用的薄平板结构，将平板离散比化为两个线性三角形元，如右图示:0.25mw0.5m12439.375KN9.375KNE=210MPa m=0.3 t=0.025 w=3000 KN/m2根据有限元求解步骤，采用线性三角形单元法求解：（1）离散比单元 节点i 节点j 节点m 1 3 4 1 2 3(2)写出单元刚度矩调用LinearTriangleElementStiffness(E,NU,t,0,0,0.5,0.25,0,0.25,1)求出k1 k2 E=210e6

44、NU=0.3 t=0.025 k1=LinearTriangleElementStiffness(E,NU,t,0,0,0.5,0.25,0,0.25,1)k2=LinearTriangleElementStiffness(E,NU,t,0,0,0.5,0,0.5,0.25,1)（3）集成整体刚度矩阵 K=zeros(8,8)K=LinearTriangleAssemble(K,k1,1,3,4)K=LinearTriangleAssemble(K,k2,1,2,3)（4）引入边界条件 yxyxyxyxyxyxyxyxFFFFFFFFvuvuvuvuk4433221144332211 边界条件

45、如下：u1x=v1y=u4x=v4y=0 F2x=9.375,F2y=0,F3x=9.375,F3y=0代入上式，并采用消去。（5）k=K(3:6,3:6)f=9.375;0;9.375;0 u=kf（6）后处理：U=0;0;u;0;0 F=K*U u1=U(1);U(2);U(5);U(6);U(7);U(8)u2=U(1);U(2);U(3);U(4);U(5);U(6)单元应力 sigmal1=LinearTriangleElementStresses(E,NU,t,0,0,0.5,0.25,0,0.25,1,u1)sigmal2=LinearTriangleElementStresses(E,NU,t,0,0,0.5,0,0.5,0.25,1,u2)S1=LinearTriangleElementPStresses(sigmal1)S2=LinearTriangleElementPStresses(sigmal2)附：matlab源代码