现代设计方法4-3三角形三节点平面单元.ppt

,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,4.3,平面问题的有限单元法,三角形三节点平面单元,结构离散化,单元分析,整体分析,有限元分析的基本步骤：,1,结构离散化,例：图示一悬臂梁，梁的厚度为,t,，,设泊松比,m,=1/3,，弹性模量为,E,，,试用三节点三角形单元进行离散。,l,2 m,1 m,A,P,/2,P,/2,3,P,/2,P,/2,4,1,2,y,i,i,j,j,m,m,x,2,单元分析,单元分析的主要内容：由节点位移求内部任一点的位移，由节点位移求单元应变、应力和节点力。,单元分析的步骤：,节点,位移,单元内部,各点位移,单元,应变,单元,应力,节点,力,(1),(2),(3),(4),单元分析,(1),单元位移模式,有限元法中，通常,采用,位移法,进行计算，即取节点的位移分量为基本未知量，单元中的位移、应变、应力等物理量，都和基本未知量相关联。,节点,i,的,位移分量,可写成,d,i,=,u,i,v,i,T,单元,节点位移向量,d,e,可写成,d,e,=,d,i,d,j,d,m,T,=,u,i,v,i,u,j,v,j,u,m,v,m,T,三角形三节点单元,x,y,m,j,i,v,i,u,i,u,j,v,j,v,m,u,m,u,v,2.1,由节点位移求单元内部任一点位移,六个位移分量需六个待定参数,a,1,，,a,2,，,，,a,6,设单元内任一点,位移分量是坐标,(,x,y,),的线性函数,：,(4.22),写成矩阵的形式为：,(4.23),设节点,i,，,j,，,m,坐标分别是,x,i,，,y,i,；,x,j,，,y,j,；,x,m,，,y,m,。把三个节点的坐标及其水平位移代入式,(4.22),中得：,(,a,),解得：,(,b,),其中，,(2),由单元节点位移,d,e,求位移参数,a,对,v,同理可列出,a,4,、,a,5,、,a,6,的,方程：,v,i,=,a,4,+,a,5,x,i,+,a,6,y,i,i,j,m,解得：,a,4,、,a,5,、,a,6,为了书写的简便与规格化，引入记号,a,i,b,i,c,i,分别为行列式,A,中各行各个元素的代数余子式,i j m,解出,a,1,a,6,结果：,a,=,L,d,e,为三角形单元面积。,将,a,写成矩阵形式，有,a,=,L,d,e,由单元节点位移,d,e,求单元内部任一点位移,d,(,x,y,),d,(,x,y,)=,m,(,x,y,),a,=,m,(,x,y,),L,d,e,i j m,N,i,N,j,N,m,是坐标的连续函数，反映单元内位移的分布状态，称为位移的,形状函数,，简称,形函数,(shape function),。,矩阵,N,称为形函数矩阵,(shape function matrix),。,形函数物理意义：,N,i,(,x,y,),，节点,i,单位位移，其它,节点位移分量为,0,，单元内部产生位移分布形状,在单元任一点上，三个形函数之和等于,1,。,形函数,N,i,在,i,点的函数值为,1,，在,j,点及,m,点的函数值为零。,三角形单元,i,j,m,在,ij,边上的形函数与第三个顶点的坐标无关。,形函数的性质：,例：求图示单元,和单元,的形函数矩阵,3,P,/2,P,/2,4,1,2,y,i,i,j,j,m,m,x,（,a,）,(b),(0,a,),(0,0),(,b,0),y,i,j,m,x,a,b,(,b,a,),(0,a,),y,i,j,m,x,(,b,0),分别如上图所示：,i j m,单元,如图所示。设,a,=1m,b,=2m.,(,或直接由图形可知其面积,),求系数,a,i,a,j,a,m,b,i,b,j,b,m,c,i,c,j,c,m,(0,a,),(0,0),(,b,0),y,i,j,m,x,a,b,i j m,求形函数矩阵,代入相关常数：,将,a,=1,b,=2,代入得：,(,b,a,),(0,a,),y,i,j,m,x,(,b,0),求常数,单元,如图所示,=,ab,/2,求形函数矩阵,将,a,=1,b,=2m,代入上式得：,作业：,已知三角形三节点单元坐标如图示，设单元中一点,A,的坐标,(0.5,0.2),，已知三角形三节点单元,i,节点位移,(2.0,1.0),，,j,节点位移,(2.1,1.1),，,m,节点位移,(2.15,1.05),，,1),写出单元的位移函数；,2),求,A,点的位移分量。,i(1,0),j(1,1),m(0,0),x,y,A,2.2,由节点位移求单元的应变,几何方程,i j m,i j m,简记为,e,=,B,d,e,B,可写成分块的形式：,B,=,B,i,B,j,B,m,B,称为应变矩阵,，它的元素都只与单元的几何性,质有关的常量。这种单元称为平面问题的,常应变三角形单元,。,单元应变与单元节点位移关系,(,i,j,m,),物理方程,s,=,D,e,而,e,=,B,d,e,s,=,D,B,d,e,=,S,d,e,(,求应力的表达式,),记,S,=,D,B,S,应力矩阵,:,S,=,S,i,S,j,S,m,2.3,由单元节点位移求单元的应力,1),平面应力问题：,代入,D,及,B,得,:,S,=,S,i,S,j,S,m,.,对于平面应力：,(,i,j,m,),将上式中以,E,/(1-,m,2,),代,E,，以,m,/(1-,m,),代,m,则子矩阵：,(,i,j,m,),2),对于平面应变问题：,考虑单元平衡，节点力是作用在单元上的外力，与单元应力平衡。有限元法中以,虚功方程,代替平衡方程。,节点力列阵及单元内应力列阵：,单元节点力是指单元和节点相连接的内力；,考虑节点平衡，节点力为外力，与节点外载荷平衡；,s,y,s,x,t,xy,2.4,由节点位移求单元节点力,(,求,单元刚度矩阵,),节点虚位移列阵及虚应变：,令实际受力状态在虚位移状态上做虚功，,虚功方程,：,由,e,=,B,d,e,知,e,*,=,B,d,*,e,则,e,*,T,=(,d,*,e,),T,B,T,由于,d,*,e,中的元素为常量，提至前积分号前，故：,(,对于,三角形三节点单元,，,B,和,s,为常量，单元厚度,t,也是常量；,为三角形单元面积，用,表示）,单元刚度矩阵,物理方程,s,=,D,e,简记为,F,e,=,k,e,d,e,k,e,=,B,T,D,B,t,或,=,B,T,S,t,几何方程,e,=,B,d,e,对于,平面应力问题,：,(,r,=,i,j,m,;,s,=,i,j,m,).,将上式中的,E,换成,m,换成,对于,平面应变问题,：,求单元,的单元刚度矩阵,(0,a,),(0,0),(,b,0),y,i,j,m,x,a,b,(,i,j,m,),解：（,1,）求矩阵,B,（,2,）求矩阵,S,（,3,）求矩阵,k,e,k,e,=,B,T,D,B,t,=,B,T,S,t,代入,a,=1,b=2m,得：,可算出，当,a,=,b,时单元刚度矩阵与尺寸,a,b,无关,在,单元节点力列阵,F,e,、,单元应力列阵,s,e,、,单元应变列阵,e,e,和,单元节点位移列阵,d,e,的四个列阵之间，存在五个转换关系，可得,五个转换矩阵,。,平衡条件,物理关系,几何关系,单元刚度矩阵的特性：,(1),单元刚度矩阵的物理意义：,单元刚度矩阵,k,e,表示了单元抵抗变形的能力，即表示了节点位移,d,e,与节点力,F,e,之间的关系。,k,ij,表示节点,j,发生单位位移时，其它节点位移分量均为零时，在节点,i,上产生的节点力。,(2),分块性质,:,单元刚度矩阵可以分块运算,。,(4.51),按节点进行分块，则单元刚度矩阵的分块形式可写为：,(4.52),（,3,）,对称性 单元刚度矩阵是一个对称矩阵,式,(4.51),中：,k,nl,=,k,ln,(,n,=16;,l,=16),分块形式中：,k,rs,=,k,sr,T,(,r,=,i,j,m,;,s,=,i,j,m,),（,4,）,奇异性 单元刚度矩阵是一个奇异矩阵,|,k,e,|,=0,，表明其逆矩阵不存在。,即，如果,给定,了,单元,节点位移可以得出唯一的节点力,:,F,e,=,k,e,d,e,。但只给,出节点力却无法求出确定的节点位移。因为这时的单元未考虑所受,约束，,可能存在不引起单元应力和节点力的刚性位移，这部分刚体位移由节点力是无法唯一确定的。,三节点三角形单元每行元素之和为零,例：,证明图示单元刚度矩阵,:,k,I,=,k,III,证明：由于单元刚度矩阵,k,e,=,B,T,D,B,t,可知：当两个三角形单元几何尺寸相同时，,t,值和单元面积,值均相同；当两个单元的材料性质,相同时，弹性矩阵,D,也时相同的。故,k,e,是否相同，取决于矩阵,B,是否相同。,不难验证,I,、,III,单元的上述,b,r,c,r,(,i,j,m,),值均相等。,结论：,两个单元刚度矩阵,k,e,相等的条件为：只要两单元的形状、大小，方向和单元弹性常数均相同，并且编号的方式也相同,(,如按逆时针方向编号为,i,j,m,，直角顶点编号为,m,),，则两个单元的刚度矩阵是相等的。,2.5,单元载荷的移置,.,（离散时每个单元受,载作用于节点上）,(a),原则：,将单元载荷向节点处移置，按照,虚功等效,的原则进行。对于变形体,(,包括弹性体,),，,虚功等效是指原载荷与节点载荷载在任何虚位移上做的虚功相等,。当位移模式确定后，载荷移置,(,或分解,),其结果是唯一的。,虚功等效包含了刚体体系的静力等效，当虚位移为刚体位移时，虚功等效即为静力等效，静力等效是虚功等效的特例。,离散,(b),载荷移置公式,(1),集中力,设单元,i,j,m,中任一点,M,(,x,y,),处受有集中力,P,=,P,x,P,y,T,，移置到该单元各节点处载荷列阵为,R,e,=,X,i,Y,i,X,j,Y,j,X,m,Y,m,T,(1),集中力,P,=,P,x,P,y,T,(1),集中力,假设该单元发生一微小虚位移，,M,点相应的虚位移为,f,*,，该单元各节点处相应虚位移为,d,*,，由静力等效原理，,载荷与节点等效载荷在虚位移上所作虚功相等,：,(,d,*,e,),T,R,e,=,f,*,T,P,将,f,*,=,N,d,*,e,代入上式：,有,(,d,*,e,),T,R,e,=,f,*,T,P,=,(,d,*,e,),T,N,T,P,则,R,e,=,N,T,P,R,e,=,X,i,Y,i,X,j,Y,j,X,m,Y,m,T,=,N,i,P,x,N,i,P,y,N,j,P,x,N,j,P,y,N,m,P,x,N,m,P,y,T,设单元,i,j,m,的一边受有分布的面力,P,=,P,x,P,y,将微元面积,t,d,s,上的面力合力,P,t,d,s,当作集中载荷，可得,面力的移置公式,:,(2),面力,(2),面力,(3),体力,设单元,i,j,m,受有分布体力,G,=,P,x,P,y,T,将微分体积,t,d,x,d,y,上的体力合力,G,t,d,x,d,y,当作集中载荷,d,G,同理可得：,(3),分布体力,G,=,X,Y,T,(4),三节点三角形单元上同时有体力、面力和集中力等,(a),集中力,P,=,P,x,P,y,T,(b),分布体力,G,=,X,Y,T,(c),面力,(d),单元虚位移,d,*,=,u,*,i,v,*,i,u,*,j,v,*,j,u,*,m,v,*,m,T,f,*,=,u,*,v,*,单元各节点处载荷列阵为,R,e,=,X,i,Y,i,X,j,Y,j,X,m,Y,m,T,应用,虚功等效原则,：,将,f,*,=,N,d,*,e,代入上式：,虚位移是任意的，从而矩阵,d,*,e,T,也是任意的，故：,例：设三角形单元,i,j,m,的,ij,边作用有线形分布的法向载荷，,i,和,j,两点的压力集度分别为,q,i,和,q,j,，试用公式,求其等效,节点载荷。单元厚度为,t,，节点坐标如图示。,解：,计算常数,a,i,=,x,j,y,m,x,m,y,j,=0,；,b,i,=,y,j,y,m,=,y,m,；,c,i,=,x,m,x,j,=,x,m,；,a,j,=,x,m,y,i,x,i,y,m,=,x,m,y,j,；,b,j,=,y,m,y,i,=,y,m,y,i,；,c,j,=,x,i,x,m,=,x,m,；,a,m,=,x,i,y,j,x,j,y,i,=0,；,b,m,=,y,i,y,j,=,y,i,；,c,m,=,x,j,x,i,=0.,计算形函数,由,得：,在边界,jm,和,mi,上的面力为零，故上式积分中后两项为,0,，在,ij,边上的面力分量可表示为：,代入上式中得：,计算等效节点载荷,=0,积分沿逆时针方向，有,d,s,=-d,y,所以,在,ij,边上,x,=0,，代入,N,i,N,j,N,m,中：,引入支承条件,解方程求位移,求单元应力,3,整体分析,建立整体刚度矩阵,任务：,建立整个结构的总刚度方程；引入边界条件解方程。,3.1,建 立 整 体 刚 度 矩 阵,3,P,/2,P,/2,4,1,2,y,i,i,j,j,m,m,I,II,x,单元节点受力分析,作用于节点的集中力,单元作用在节点上的力,作用在节点上等效载荷,2,由节点的平衡条件，有平衡方程：,环绕节点,2,的所有单元求和：,F,2,=,U,2,e,V,2,e,T,单元,e,作用在,2,节点上的节点力,R,2,=,X,2,Y,2,T,绕节点,2,的各单元作用在,2,节点,上的等效节点载荷及直接作用于,2,的集中力之和。,每一个单元可用节点位移表示节点力，采用分块,P89 4.46,考虑节点,2,，对应,I,单元的局部编号为,i,即,I,单元,i,节点，故,I,单元,i,节点上的节点力：,代入,得,结构的整体刚度矩阵,结构的节点位移列阵,结构的节点载荷列阵,每个节点均可得到类似的方程，每个节点可写出两个平衡方程，按节点的序号排列，用矩阵表示：,结构的节点载荷列阵,R,=,X,1,Y,1,X,2,Y,2 ,T,故,其中,l,=14,l,k,=14,整体刚度矩阵的集成规则：,(1),先求出每个单元的刚度矩阵,k,e,；,(2),将,k,e,的每个方块,k,ij,e,进行换号，换成对应的整体编号,;,(3),将换号后的子块填入整体刚度矩阵上对应的位置。,(4),若在同一位置上有几个单元的相应子块填入同一位置，则进行叠加。,图,4.24,刚度矩阵,单元号,局部号,i,j,m,i,j,m,整体号,2,4,1,4,2,3,(,书上错误,),单元刚度矩阵,k,e,的子块换号,k,22,k,24,k,21,k,44,k,42,k,43,k,42,k,44,k,41,k,24,k,22,k,23,k,12,k,14,k,11,k,34,k,32,k,33,3,P,/2,P,/2,4,1,2,y,i,i,j,j,m,m,x,i j m,2 4 1,i j m,4 2 3,i,2,j,4,m,1,i,4,j,2,m,3,解,:(,1),求各单元刚度矩阵,(P90,已求,)P90 4-50,求得单元,的刚度矩阵为：,k,I,=,k,II,(I II,单元刚度矩阵相等,),例：,4.2,。,(2),整体刚度矩阵,k,换号的单元刚度矩阵：,(3),将,、,的单元刚度矩阵填入,(P101,式,4.62),如果,I II,都有就相加便于计算机运算,总体刚度矩阵的特性：,(1),对称性,总体刚度矩阵是一个对称矩阵。因单元刚度矩阵升阶后对称性不变，由之合成的总体刚度矩阵自然是对称矩阵。,(2),奇异性,总体刚度矩阵行列式的值,|,K,|,=0,(3),稀疏性,总体刚度矩阵是一个稀疏矩阵；即矩阵中的绝大多数元素为,0,，非,0,元素只占元素总数的很小的一部分。因为只有当节点,ln,相关时,K,ln,才不是,0,，,P101,图,4.26,。,(4),带状分布规律,分布在以主对角线为中心的带状区域内。,总体刚度矩阵的特性：,3.2,节点载荷列阵,集成法求整体结构的节点载荷列阵步骤：,(1),求出每个单元的等效节点载荷；,(2),单元等效节点子向量,(,k,=,i,j,m,),换号，换成对应整体编号；,(3),换号后等效节点载荷子向量送到整体节点载荷列阵对应位置；,(4),同一位置上若有多个单元的等效节点载荷子向量，叠加；,(5),节点上有直接作用的节点载荷，按整体节点号进行叠加；,(6),若节点,k,具有水平和垂直方向的支承，支承反力为未知量，可暂设为,R,k,e,=,Q,xk,Q,yk,T,。,图,4.24,中：,R,=,Q,x,1,Q,y,1,X,2,Y,2,X,3,Y,3,Q,x,4,Q,y,4,T,=,Q,x,1,Q,y,1,0,-,P,/2,0,-,P,/2,Q,x,4,Q,y,4,T,3,P,/2,P,/2,4,1,2,y,i,i,j,j,m,m,x,3.3,引入支承条件,2,Y,2,1,3,(4.66),例：,基本未知量为节点,2,的位移，只需在方程中抽出第三、四行即可：,矩阵形式,(4.67),为了便于编程，修改后的矩阵仍保持原矩阵,K,的阶数、排列次序及矩阵对称性。式,(4.67),扩大成：,等效：,u,1,=,v,1,=,u,3,=,v,3,=0,与支承条件及前述矩阵一致,多个节点：,节点,n,的水平位移,u,n,=0,，,则,K,d,=,R,改为：,K,中的,2,n,-1,行与列中主对角元素为,1,，其他为,0,载荷向量,R,中,2,n,-1,元素置,0,若,v,n,=0,则对应,K,中,2,n,行和列,作上述修改。,仍以课本,P98,图,4.24,u,1,=,v,1,=,u,4,=,v,4,=0,3,P,/2,P,/2,4,1,2,y,i,i,j,j,m,m,x,3.4,解方程,求节点位移,K,d,=,R,将整体刚度矩阵,K,代入支承条件,图,4.24,3,P,/2,P,/2,4,1,2,y,i,i,j,j,m,m,x,设,m,=1/3,得：,得：,于是节点的位移向量为,d,=,u,1,v,1,u,2,v,2,u,3,v,3,u,4,v,4,T,=,P,/,Et,0 0 -1.50 -8.42 1.88 -8.99 0 0,T,单元,1,：当,a,=1,b,=2,m,=1/3,时，由例,4.3,3.5,求单元应力,s,e,=,S,e,d,e,同理可求单元,单元,I,t,xy,t,xy,s,x,s,y,0.281,P,/,t,1.58,P,/,t,1.580,P,/,t,0.844,P,/,t,单元,II,t,xy,t,xy,s,x,s,y,0.289,P,/,t,0.418,P,/,t,0.418,P,/,t,0.844,P,/,t,3.6,求节点力及支承反力,F,e,=,k,e,d,e,单元,的节点力,a,=1,b,=2,m,=1/3,3,P,/2,P,/2,4,1,2,y,i,i,j,j,m,m,x,单元,的节点力,同理可得,F,II,也进行验算,;,由受力图可得,验算,i,j,m,0.289,P,1.580,P,2.000,P,1.070,P,0.789,P,0.422,P,0.004,P,i,j,m,0.498,P,0.418,P,0.289,P,0.209,P,0.422,P,总结：,有限元求解弹性力学平面问题步骤如下：,(1),整理原始数据，结构离散化，对单元和节点编号。,(2),求单元刚刚度矩阵,k,e,(3),用刚度集成法，形成结构整体刚度矩阵,K,(4),求节点等效载荷，写出载荷列阵,R,(5),引入支承条件,(6),解,K,d,=,R,求出节点位移,d,(7),求单元应力,s,e,=,S,e,d,e,(8),求单元节点力,F,e,=,K,e,d,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.7559,4,线性三角形弹性力学平面问题的,matlab,程序：,(,德,),卡坦,-,韩来彬,清华大学出版社,若不让,matlab,输出运算数据，在命令行的结尾输入分号：,y=32;,z=5;,x=2*y-z;,w=3*y+4*z,w=116,matlab,区分大小写,x=1,x=1,X=2,X=2,x,x=1,使用,help,命令可以获得所有,matlab,命令的详细用法,help inv,下面的例子显示如何输入矩阵并实现简单的矩阵运算。,x=1 2 3;4 5 6;7 8,1 2 3,4 5 6,7 8 9,x=,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.5926,A=2 -1 3 0,1 5 -2 4,2 0 3 -2,1 2 3 4,采用如下方法也可求,x=inv(A)*b,x=1.9259,-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.,再用,plot(x,y),绘图。,x=1 2 3 4 5 6 7 8 9 10,y=x,.,2,plot(x,y),二、线性三角元,1.,基本方程,.,Linear triangular element,线性形函数,.,系数：弹性模量,E,泊松比,m,厚度,t,单元刚度阵：,k,=,B,T,D,B,t,2,=,x,i,(,y,j,-,y,m,)+,x,j,(,y,m,-,y,i,)+,x,m,(,y,i,-,y,j,),a,i,b,i,c,i,(P79 4.26).,平面应力,平面应变,显然线性三角形元有,6,个自由度,每个节点,2,个自由度。,对一个有,n,个节点的结构，整体,K,为,2n*2n,K,d,=,F,边界条件手动赋值,采用高斯消去求解,s,=S,d,求得单元应力矩,2,，用到的,Matlab,函数,线性三角形元用到的,5,个,matlab,函数分别为：,(1)LinearTriangleElementArear(xi,，,yi,，,xj,，,yj,，,xm,ym),该函数根据给出的,i,j,m,的节点坐标返回单元的面积。,j,m,节点坐标已知的线性三角形单元的刚度矩阵。,p=1,表明函数用于平面应力情况。,P=2,表明用于平面应变情况。返回,6*6,的单元刚度矩阵,k,。,(3)LinearTriangleElementAssemble(K,k,i,j,m),将连接,i,j,m,线性三角形元的单元刚度矩阵,k,集成到整体刚度矩阵,K,。返回,2n*2n,的整体刚度矩阵,K,。,(4)LinearTriangleElementStresses(E,NU,t,x,i,y,i,x,j,y,j,x,m,y,m,p,u),u,为单位位移矢量,返回单元应力,(5)LinearTriangleElementPStresses(sigma),计算单元主应力，返回,3*1.sigma1,sigma2,theta,T,例：求解如图所示的受均布载荷作用的薄平板结构，,将平板离散比化为两个线性三角形元，如右图示,:,E=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,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,）引入边界条件,边界条件如下：,u,1,x,=,v,1,y,=,u,4,x,=,v,4,y,=0,F,2,x,=9.375,F,2,y,=0,F,3,x,=9.375,F,3,y,=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,源代码,