偏微分课程课件13-有限元离散方法(二维问题).ppt

1、第八章第八章 有限元离散方法有限元离散方法 （二维问题）（二维问题）8.2 二维边值问题有限元方法二维边值问题有限元方法,(,),(,)ufx yuugx yn 1）讨论对象）讨论对象：椭圆型偏微分方程边值问题：椭圆型偏微分方程边值问题 2212()|vvVHvvdxdyxy 2）转化为变分问题）转化为变分问题.(,)()0,find uV s tD u vF vvV 其中其中(,)()D u vuvdxdyuvdsF vfvdxdygvds uuvdxdyvdsfvdxdyn ugun hVV是是(,)()D u vF vvV 找一个函数找一个函数uV,满足满足(,)()hhhhhD uvF

2、 vvV 找一个函数找一个函数 ，满足，满足3）有限元离散）有限元离散的有限维子空间的有限维子空间下面通过找基函数来构造有限维子空间下面通过找基函数来构造有限维子空间2212()|vvVHvvdxdyxy hhuV uV 定义区域离散定义区域离散1）单元编号，节点编号，节点坐标）单元编号，节点编号，节点坐标 3）边界点相关信息）边界点相关信息2）节点的局部和整体编码对应关系）节点的局部和整体编码对应关系三角形网格剖分三角形网格剖分ie uV 定义区域离散定义区域离散单元单元he 2 任一个单元的顶点不能是其他三角形的内点或任一个单元的顶点不能是其他三角形的内点或除顶点之外的边界点除顶点之外的边

3、界点1 区域被分割成有限个互不重叠的三角形单元区域被分割成有限个互不重叠的三角形单元区域三角形网格剖分要求：区域三角形网格剖分要求：3 单元最小角有下界，尽量大些单元最小角有下界，尽量大些4 任一个单元至多有两个顶点在边界上任一个单元至多有两个顶点在边界上 5 解变化剧烈处结点密些，变化平缓处结点疏些解变化剧烈处结点密些，变化平缓处结点疏些内部节点内部节点边界节点边界节点顶点顶点(节点节点)编号编号pN顶点总数顶点总数(,)ssxy坐标坐标sP(,),1ssspuu xysN 1 2 3 4 5 6 71 3 5 7 9 11 138 9 10 11 12 13 142 4 6 8 10 12

4、 14结点编号原则：（结点编号原则：（1）单元内）单元内3个编号差尽可能的小个编号差尽可能的小例如上例中，两图单元划分相同，且节点总数都等于例如上例中，两图单元划分相同，且节点总数都等于14，但两者的节点编号方式却完全不同。但两者的节点编号方式却完全不同。(a)图结点编号是按长边图结点编号是按长边进行编号，编号差进行编号，编号差d=7，若采取带宽压缩存储，若采取带宽压缩存储,则刚性矩阵存则刚性矩阵存储量储量N=14*8=112；而；而(b)图是按短边进行编号，图是按短边进行编号，d=2，N=42。显然。显然(b)的编号方式可比的编号方式可比(a)的编号方式节省的编号方式节省70个存储单元。个存

5、储单元。(a)(b)（2）先内点后边界点编号）先内点后边界点编号单元刚度矩阵单元刚度矩阵 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 单元编号单元编号hiie eN单元总数单元总数(,)i j m结点编号原则：结点编号原则：先内点后边界点，用先内点后边界点，用3个结点号逆时个结点号逆时针顺序表示针顺序表示(,)(),hhhhhD uvF vvV 定义空间：定义空间：,?a b c(,)()D u vF vvV 找一个函数找一个函数uV,满足满足找一个函数找一个函数 满足满足有限元离散有限元离散子

6、空间子空间-线性插值线性插值 ():ihhheVvCvaxbyc hhuV 定义空间：定义空间：():ihhheVvCvax by c (,)iiiuu xy iiijjjmmmaxbycuaxbycuaxbycu 记记那么那么a,b,c满足满足(,)hux yaxbyc (,)(,)(,)iijjmmNx y uNx y uNx y u (,)nx ye iPnejPmPnijmeP P P 任取一个单元任取一个单元 ,iiijjjmmmPxyPxyPxy11(,)121ijjemmxyNx yxyxy 11121iiejjmmxyxyxy iPjPmP11(,)121jmmeiixyNx

7、yxyxy 11(,)121miiejjxyNx yxyxy (,)iNx ye上线性插值的基函数上线性插值的基函数单元单元e的面积的面积(,)x yxyiNiPjPmP10,(,)1,illlj mNxyli 12iiiea xb yc 11,11jjiimmjjimmyxabyxxycxy 11(,)121ijjemmxyNx yxyxy 0,(,)1,illlj mNxyli 即，子空间即，子空间Vh定义为定义为(,)(,)(,)hiijjmmuNx y uNx y uNx y u ():(,)(,)(,)ihiijjmmeVu CuN x y uN x y uN x y u ,ijmN

8、NNN ,Teijmuu u u eNu nnnnnnnnhheeeeeh hhheeeuv dxdyu v dsfv dxdygv ds 4）化为代数方程组）化为代数方程组(4.1)单元分析单元分析(4.2)总体合成总体合成(,)(),hhhhhD uvF vvV 找一个函数找一个函数hhhhhhuv dxdyu v dsfv dxdygv ds hhuV 满足满足nijmeP P P 任取一个单元任取一个单元单元分析单元分析(,)(,)(,)hiijjmmuNx y uNx y uNx y u neNv neNu 记函数记函数 和和 在结点在结点的取值分别为的取值分别为sP,ssu vsi

9、 j m(,)(,)(,)hiijjmmvNx y vNx y vNx y v ,iiijjjmmmPxyPxyPxy ,niijmejmuNNNNuuuhuhvhhhuxuuy (,)(,)(,)hiijjmmuNx y uNx y uNx y u (,)(,)(,)(,)(,)(,)jimijmjimijmN x yN x yNx yuuuxxxN x yN x yNx yuuuyyy hhhvxvvy hhhhuvuvxxyy hhuv nneeBuBv ,njimiejjimmNNNuxxxBuuNNNuyyy nhevNv nnheeuNu 在在上上()()()()nnnnnnnnn

10、nnnnnTeeeeTeeeeTTeeeeeeBvBu dxdyNvNu dsNvfdxdyNvgds nnnnnnnnhhhheeeehheeeeuv dxdyu v dsfv dxdygv ds()nnhhTeeu vNvNu nnheheuBuvBv()nnhhTeeuvBvBu nnnnnnnnnnnnnnTTeeeeTTeeeeTTTTeeeeeevB Bdxdy uvN Nds uvN fdxdyvN gds eKeKeFeF nnnnnnnnTTeeeeeeeeeKKFFvuv neK 单元刚度矩阵neF 单元荷载面元刚度矩阵面元刚度矩阵线元刚度矩阵线元刚度矩阵面元荷载向量面元荷

11、载向量线元荷载向量线元荷载向量14nnnTeeiiijmjjeijmmmeeeiiijimeeejijjjmeeemimjmmKB Bdxdyabaaaabbbbabkkkkkkkkk 1,.4estststeka ab bs ti j m 其中面元刚度矩阵的计算面元刚度矩阵的计算 1,0.ijijijiP PjP PmP PtNltNlN nijmeP P P ijnePP 设设,:0,:,ijijPPlPtPtl (,)1,(,)0,iiiijjNxyNxy 0nnnnnnnnnnlTeeeeiiijimeeejijjjmeeemimjmmKN Ndtkkkkkkkkk 200201,1,

12、0.nnnnnnnnnleiileeijjilejjeeeeemimjmmjmimtkdtlttkkdtlltkdtlkkkkk 其其中中线元刚度矩阵的计算线元刚度矩阵的计算 00,.1,0nnnnnnnnnnnneieeTejiseeemleileTejemFFNfdxdyFFN fdxdy si j mFtgdtlFtFN gdsFgdtlF 其中面元荷载的计算面元荷载的计算线元荷载的计算线元荷载的计算 312123123123!22!eeNNNdxdy 可可以以利利用用的的公公式式总体合成总体合成 ,nnTTeijmeijmuu u uvu u u nijmeP P P 记函数记函数 和

13、和 在结点在结点的取值分别为的取值分别为,1,sspu vsN sP 1212,ppTTNNuu uuvv vv nnnnnnnnnTeeeeeeeeevKKuFFv其中其中huhvnniejemvvvC vv300100000100010nnepeCNC 记记为为的的矩矩阵阵)(95m101()i102()jiPjPmP10110295则则()nnnnnnTTTeeeeeeC vK C uv C K C u ,nnnnnnnnnnneeeiiijimieeeTjijjjmeeeijmjeeemmimjmmkkkuvK uv v vkkkuukkk 0000000000000000nnnnnn

14、nnneeeiiijimieeejijmjijjjmeeemmimjmmkkkuuvvvkkkukkk 000000000000000000000nnnneeiiijieejijmjijjmkkuuvvvkku nijmeP P P ijnePP 设设 0()()nnnnnnlTTTehheeeeu v dsNvNudsv C K C u ,()nnnnnnnnnnnTeTTeeTeeeijmeeTeeeTTTeeijmv FvF F FC v Fv C FvFFF =,00nnnnnnnnnTTeeeeTTTTTeeeijeeijv FvFFv C FvFF nijePP nnnnnnnnn

15、TTeeeeeeeeevKKuvFF nnnnnneeeTTeeevuKKFFv TTvuKFv nnTeeFN fdxdy nnTeeFN gds 若若讨论对象讨论对象：椭圆型偏微分方程边值问题：椭圆型偏微分方程边值问题 10222()|,0VHvvvvdxdyvxy 转化为变分问题转化为变分问题 5）约束边界处理）约束边界处理.(,)()0,find uV s tD u vF vvV 其中其中(,),()D u vuvdxdy F vfvdxdy 0():,1,2,0,iihhhhePVvCvaxvilby c 0(,)(),hhhhhD uvF vvV 找一个函数找一个函数 满足满足0h

16、huV 有限元化有限元化试探函试探函数空间数空间边界节边界节点编号点编号 111122122212222210,0,0,0,0,0,plNIIIIIIIIIITIIIIIITIIIIIIIIIIIIvvvuFvuFvuFKKuFvKKuFvK uK uFK uFK u 若不在开头，总刚矩阵划去边界节点相应的行与列若不在开头，总刚矩阵划去边界节点相应的行与列TTv KuFv n有限元法内容总结有限元法内容总结1 1 明确微分方程明确微分方程2 2 转化为变分形式；转化为变分形式；3 3 有限元离散；有限元离散；区域剖分；确定单元基函数区域剖分；确定单元基函数；4 4 转化为代数方程转化为代数方程

17、单元分析；总体合成；边界条件的处理单元分析；总体合成；边界条件的处理；5 5 解代数方程。解代数方程。总结：实际计算步骤总结：实际计算步骤1）单元剖分）单元剖分nijmeP P P 节点编号和坐标节点编号和坐标 单元编号与对应节点单元编号与对应节点(i,j,m)边界条件节点编号边界条件节点编号2）有限元离散）有限元离散基函数基函数单元号和单元顶点单元号和单元顶点i处参数处参数(ai,bi,ci)计算单元刚度矩阵、单元荷载向量计算单元刚度矩阵、单元荷载向量 合成总刚度矩阵与总荷载向量，形成代数方程组合成总刚度矩阵与总荷载向量，形成代数方程组3）根据边界调整）根据边界调整K和和F，即，即“掐头或去

18、尾掐头或去尾”4）求解上代数方程组，得数值解）求解上代数方程组，得数值解202020,(,)(0,2)50,1000,0yyxxux yuuuuxx 例题：例题：分析分析解：解：平面区域上无热平面区域上无热源的定常温度场，源的定常温度场，两边给出温度值，两边给出温度值，另两边给出绝热另两边给出绝热条件条件采用三角形剖分，采用三角形剖分，共共16个单元，个单元，15个节点个节点11()()()()nnnnnnnnnnnnnnTeeeTeeeeTeeeTeeeeBvBu dxdyNvfdxdNvNu dsNvgdsy 其中其中1(,)|02,02x yxy 202020,(,)(0,2)50,10

19、00,0yyxxux yuuuuxx 1）节点坐标与编号）节点坐标与编号 1(,)2,sssseNx ya xb ycsi j k 2）基函数与相应参数）基函数与相应参数单元单元信息信息结点号结点号(i,j,k)参数值参数值(4,1,5)(-0.5,1)(0,-1)(0.5,0)(2,5,1)(0.5,-1)(0,1)(-0.5,0)(5,2,6)11,11jjjjiiimmmmyxxyabcyxxy 通过计算可知，奇数单元上参数值相同通过计算可知，奇数单元上参数值相同,偶数单元上偶数单元上参数值相同参数值相同单元单元信息信息结点号结点号(i,j,k)(4,1,5)(-0.5,1)(0,-1)

20、(0.5,0)(2,5,1)(0.5,-1)(0,1)(-0.5,0)(5,2,6),()()()iijjkka ba ba b1541541415415(,)(,)(,)(,)(,)(,)0.500.51 2 2110eeNx yNx yNx yxxxBNx yNx yNx yyyyaaabbb 1212()0.51iiejjmmxyxyarea exy 1541440101edxdy 54114404101 4 15 4 1 5单元单元1刚度矩阵刚度矩阵154114404101eK 奇数单元刚度矩阵奇数单元刚度矩阵奇数单元上的刚度矩阵等于奇数单元上的刚度矩阵等于 y 13 14 15 (1

21、3)(14)(15)(16)(1)(2)(3(1)(2)(3）(4)(4)1 2 3 0 1 2 x354114404101 eK 1eK1315,eeeKKK5 26 5 2 62222TeeeeKB B dxdy 54114404101 2 51 2 5 12251251 0.500.512 2110eeaaaBbbb 单元单元2的刚度矩阵的刚度矩阵254114404101eK 3 62 3 6 2偶数单元刚度矩阵偶数单元刚度矩阵偶数单元上的刚度矩阵等于偶数单元上的刚度矩阵等于 y 13 14 15 (13)(14)(15)(16)(1)(2)(3(1)(2)(3）(4)(4)1 2 3

22、0 1 2 x454114404101eK 2eK2416,eeeKKK12345 1 2 3 4 00000000110154 01400005 00004 4 扩充单元刚度矩阵扩充单元刚度矩阵254114404101eK 2 51 2 5 1叠加得总体刚度矩阵叠加得总体刚度矩阵 y 13 14 15 (13)(14)(15)(16)(1)(2)(3(1)(2)(3）(4)(4)1 2 3 0 1 2 x半带宽为半带宽为4带宽为带宽为7对称对称3）边界处理）边界处理0250,100yyuu y 13 14 15 (13)(14)(15)(16)(1)(2)(3(1)(2)(3）(4)(4)1

23、 2 3 0 1 2 x12350uuu 131415100uuu 123131415,TTIIIIIIuuuuuuuuu 0Tv Ku 1112132122233132330,v,00ITIIIIIIIKKKuKKKuKKKu 2122230TIIIIIIIIvKuKuKu 222123IIIIIIKuKuKu 21123231314154,8,4,0000,4,8,4TITIIIKuuuuKuuuu 1112132122233132330,v,00ITIIIIIIIKKKuKKKuKKKu其中其中总体刚度矩阵总体刚度矩阵4567891011121020450220208100021000450400102040148022020804402100040400100010080020220040210100uuuuuuuuu 得代数方程组得代数方程组利用软件求解即可利用软件求解即可课堂练习课堂练习1、请给出有限元法程序设计流程。、请给出有限元法程序设计流程。2、