有限元课程设计

数学与统计学院信息与计算科学系《有限元课程设计》姓名：学号：专业：班级：指导老师：时间：成绩评定指导教师评语与成绩指导教师：时间：学号姓名成绩设计题目：（1）平面应力问题上的单元节点位移求法（2）双曲方程满足Dirichlet的有限元方法解答过程摘要:随着计算机技术的飞速发展，编制高效的程序求解各种偏微分方程问题已经成为可能，偏微分方程数值解法也成为科学和工程计算中的重要分支。然而，对于广大应用工作者来说，从偏微分方程模型出发，使用有限元法或有限差分法求解都要经过诸多步骤，耗费很大的工作量，才能得到数值解。进一步将结果可视化也是一个迫切需要解决的问题。目前，MATLAB中的PDEtoolbox已实现了二维偏微分方程定解问题的准确高效求解，人们只需使用图形用户界面或编写M文件，即可显示解的图形或输出解的数值。MATLAB为求解偏微分方程带来了极大方便，并在此基础上可以解决更为复杂的问题，偏微分方程在科学和工程上有着广泛的应用。有限元法是一种重要的偏微分方程数值解法。编程实现从偏微分方程到有限元求解全过程需要很好的理论基础和编程技巧,难度较高,PDEtoolbox实现了偏微分方程的有限元解法。关键词：MATLAB，PDEtoolbox，有限元问题背景：1.如下图所示的结构受集中力P作用，试按平面应力问题分析，用3节点三角形单元，求图中结构2，4节点位移，取E为常量，μ=1/6，t=1。2.详细推导出双曲方程满足Dirichlet的有限元方法解答过程，并求解正方形区域}1,1),{(yxyx上的波动方程utu22，初始条件))(arctan(cos)0(xu，))exp(cos()sin(3)0(yxdtdu，边界条件1x，Dirichlet条件0u，1y，Neumann条件0nu。详细解答第一题：第一步：对单元①，i,j,m对应节点3,2,11ijmbyy1jmibyy1mijbyy1m1m①②1m1p2341imjcxx0jimcxx1mjicxxA=1/2几何矩阵B为：00010101010000100012110111ijmijmiijjmmbbbBcccAcbcbcb弹性矩阵D为：211061036110101356150000212uEEDuuu利用公式：TrsrskBDBtA,,,rsijm求得：1111111211212113113131222211322320.72860.300010.30000.728620.51430.214310.08570.214320.30000.128610.12860.300020.51430100.2143212TTTTkBDBEkkBDBEkkBDBEkBDBEkkB3133330.51430.08570.21430.21430.72860.300010.30000.72862TTDBEkBDBE从而求得：11111121311112122231113132330.72860.30.51430.21430.30.12860.30.72860.08570.21430.12860.30.51430.21430.514300.51430.08570.08750.214300.21430.21430.21430.30.12860.514kkkkkkkEkkk30.0850.72860.30.12860.30.21430.21430.30.7286单元贡献矩阵11111121311112122231113132330000000kkkkkkkkkk注意：这里的上标代表单元的号码。第二步：对于单元②，i,j,m对应节点3,4,2.101101ijmimjjmjjimmijmjibyycxxbyycxxbyycxxA=1/2几何矩阵B为：00010100010000001012011110ijmijmiijjmmbbbBcccAcbcbcb利用公式：TrsrskBDBtA,,,rsijm求得：23333224334342223323224444222442420.51430100.214320.51430100.2143200.214310.0857020.72860.300010.30000.728620.21430.214310.02TTTTTkBDBEkkBDBEkkBDBEkBDBEkkBDBE222228570.51430.21430100.51432TkBDBE从而求得：22222232422223233342224243440.2143000.21430.21430.214300.51430.085700.08570.514300.21430.514300.514300.0857000.214300.21430.21430.21430.514300.72860.30000.08570.514300.kkkkkkkEkkk21430.30000.7286单元贡献矩阵：22222232422223233342224243440000000kkkkkkkkkk第三步：求得总刚度矩阵：1111112131121212122222323241121223132323333342224243440.72860.30000.51430.21430.30000.1286000.30000.72860.08570.21430.12860.3000000.51430.08570.72860000kkkkkkkkkkEkkkkkkkkk.51430.30000.21430.21430.21430.214300.72860.30000.21430.08570.51430.30000.12860.51430.30001.24290.30000.51430.08570.12860.30000.30000.21430.30000.94290.21430.2143000.21430.08570.51430.21430.72860.3000000.21430.51430.08570.21430.30000.7286有了总体刚度矩阵后，再形成载荷矩阵，即可得整体刚度方程，约束处理后就可求解节点位移。载荷矩阵为：11334402TpFuvuvuv位移列阵为：22440000Tuvuv形成整体平衡方程：0.72860.30000.51430.21430.30000.1286000.30000.72860.08570.21430.12860.3000000.51430.08570.728600.51430.30000.21430.21430.21430.214300.72860.30000.21430.08570.51430.30000.12860.51430.3E2240000001.24290.30000.51430.085700.12860.30000.30000.21430.30000.94290.21430.2143000.21430.08570.51430.21430.72860.3000000.21430.51430.08570.21430.30000.7286uvuv113344402/uvpuvuv处理后得到：201.00000000000001.0000000000001.68560.6210000.22180.8428000.62103.1760000.53232.643700001.0000000000001.000000000.22180.5323001.80981.0557000.84282.6437001.05573.9212uvE244000/2000000puv解以上；联立方程得22440.31051.5880.26621.3218uvpuEv第二题正方形区域}1,1),{(yxyx上的波动方程utu22，初始条件))(arctan(cos)0(xu，))exp(cos()sin(3)0(yxdtdu，边界条件1x，dirichlet条件0u，1y，Neumann条件0nu。推导：有限元法的基本步骤：[1]把微分方程定解问题转化为变分形式[2]选定单元的形状，对求解区域做剖分[3]构造基函数或单元形状函数[4]形成有限元方程[5]求解有限元方程[6]收敛性及误差分析考虑热传导方程的初边值问题(1)0,,,tGyxutu(2)Gyxyxyxu,,,0,,(3)0,0|tuG其中2222yuxuu，G是光滑边界G的平面有界域。[1]把微分方程定解问题转化为变分形式设对固定的t0tyxu,,关于(x,y)属于_2GC以GH10v乘(1)两端并积分，得(4)Gvdxdyutu0利用Green公式和边值条件(3)，得(5)Gdxdyyvyuxvxuvtu0引进双线性形式和GL2内积：(6)Gdxdyyvyuxvxuvua,则初边问题(1)-(3)的变分形式为：求GHtyxu10,,，满足(7)GHvuavtu10v,0,,(8)yxgyxuyxgyxut,0,,,,0,,10[2]对区域G作网格剖分在10H中取一n维子空间hU,0;,htUtyxuh从而Galerkin方程为(9)hhhhhUvuavtuhv,0,,(10)hhhhhhhhUvvgvyxtuvgvyxu),,(,0,,),,(,0,,10[3]构造基函数或单元形状函数在单元1,iiII考察线性差值公式及iu的系数，对每一节点ix构造山形函数(11)其他其他,0,11,,.........2,1i,0,1,11111nnnnniiiiiiiiixxxhxxxnxxxhxxxxxhxxx显然可表成的基底。任一线性相关，它们组成hh1u.,,.........UU