搜索
薄板问题有限元分析例题例题：正方形薄板（如图1所示）边长2l=2m，厚度t=0.01m，材料弹性模量E=200GPa，泊松比μ=0.3。薄板四边固支，其中心受垂直集中载荷P=400N作用，试求薄板中心挠度w。【理论解()220.00560.00048922mPlwD=−=−（()32121EtDμ=−）】图1薄板结构图解：求解中采用国际单位制。完整Matlab程序%Step1划分网格Enode=[11254;22365;35698;44587];%每一行表示一个单元的信息，第1个数字表示单元编号，第2-5个数字表示该单元的结点编号。ndof=3;%薄板单元每个结点的自由度是3。Edof=caldof(Enode,ndof);%计算各单元结点所对应的自由度编号。ex=[0110;1221;1.单元划分：单元1（结点1、2、5、4）单元2（结点2、3、6、5）单元3（结点5、6、9、8）单元4（结点4、5、8、7）第1页共16页1221;0110];%每一行表示一个单元4个结点的x坐标值。ey=[0011;0011;1122;1122];%每一行表示一个单元4个结点的y坐标值。%Step2求单元刚度矩阵E=2.0e11;v=0.3;D=hooke(E,v);%E是弹性模量，v是泊松比，D是弹性矩阵。t=0.01;ep=[t];%ep=[t]表示薄板的厚度。nie=size(Enode,1);%nie表示单元总数。fori=1:nieKe(:,:,i)=platre(ex(i,:),ey(i,:),ep,D);End%Ke(:,:,i)表示第i个单元的单元刚度矩阵。2.单元刚度矩阵分析：薄板单元在局部坐标系X’O’Y’下的刚度矩阵按照以下公式计算[][][][]312eTAtKBDBdxdy=∫∫式中[]B和[]D的计算公式在附录1中。%Step3求总体刚度矩阵%利用直接刚度法自动进行刚度矩阵组集K=zeros(27);%初始化整体刚度矩阵，所有元素置为0。K=assem(Edof,K,Ke);%组集得到总体刚度矩阵。3.按单元结点编号所对应的位置进行组装可以得到整体刚度矩阵第2页共16页%Step4带入载荷及边界条件，求解总体刚度方程，得到结点位移a和支反力Q结果。f=zeros(27,1);%初始化载荷向量，所有元素置为0。f(13)=-400;%第13个自由度有外载荷-400N。bc=[10;20;30;40;50;60;70;80;90;100;110;120;160;170;180;190;200;210;220;230;240;250;260;270];%第1、2、3、4等自由度有0位移边界条件。[a]=solveq(K,f,bc);%求解刚度方程，得到结点位移结果a。4.采用降阶法求解刚度方程%Step5后处理Ed=extract(Edof,a);%提取各单元结点的位移计算结果。5.后处理板中心挠度是-0.00051705m。【理论解是-0.00048922m】第3页共16页附录1：四结点矩形薄板单元刚度矩阵薄板单元位移模式假设为：22123456322333789101112223356891011122222457891112(,)()()()(22332323xywxyxyxxyyxxyxyyxyxywxyxxyyxxyyw23)xyxxyyxyxααααααααααααθααααααααθααααααα⎧=+++++⎪++++++⎪⎪∂⎨==+++++++∂⎪⎪∂⎪=−=−−−−−−−−∂⎩yα2根据单元4个结点坐标和结点位移代入上式，可求得11~αα，将它们代回上式，整理得：{},,,,,,(,)()[]kkkxkxkykykkkijmpkijmpwxyNwNNNθθδ===++=∑∑或：{}(,)[]ewxyNδ=薄板单元形函数矩阵[][][][][]ijmNNNNN⎡⎤=⎣⎦p[](,,,)iixiyiNNNNijm⎡⎤=⎣⎦p221(1)(1)[2(1)(1)]8(1)(1)(1)(,,,)8(1)(1)(1)8iiiiiiiiixiiiiiyiiixyxxyyNxyxxyyyxyyNixyyxxxyNxxy⎧=+++−+−⎪⎪⎪⎪=−+−+⎨⎪⎪=−++⎪⎪⎩jmp薄板单元应变矩阵[][]222222222222ijmpijmpxxBBBBBNNNNNyyxyxy⎧⎫⎧⎫∂∂⎪⎪⎪⎪∂∂⎪⎪⎪⎪⎪⎪⎪⎪∂∂⎪⎪⎪⎪⎡⎤⎡==−=−⎨⎬⎨⎬⎣⎦⎣∂∂⎪⎪⎪⎪⎪⎪⎪⎪∂∂⎪⎪⎪⎪∂∂∂∂⎪⎪⎪⎪⎩⎭⎩⎭⎤⎦以[]iB为例，第4页共16页[]32332322233231626101881626110881211311312344iiiiiiiiiiiiiiiiiiiiiiiiiixxyyxyyxxyyxxyByxxyyyxyyyxyyxxxyyxyyy⎛⎞⎛⎞⎛⎞⎛⎞−++−−⎜⎟⎜⎟⎜⎟⎜⎟⎝⎠⎝⎠⎝⎠⎝⎠⎛⎞⎛⎞⎛⎞⎛⎞=−−+−+−−⎜⎟⎜⎟⎜⎟⎜⎟⎝⎠⎝⎠⎝⎠⎝⎠⎛⎞⎛⎞⎛⎞⎛⎞+−+−−−−⎜⎟⎜⎟⎜⎟⎜⎜⎟⎝⎠⎝⎠⎝⎝⎠x2231234iiiiixxxyxxx⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎛⎞−−⎢⎥⎟⎜⎢⎥⎟⎠⎝⎠⎣⎦薄板单元刚度矩阵[][][][][][][][]312iiijimipjijjjmjpeTAmimjmmmppipjpmppKKKKKKKKtKBDBdxdyKKKKKKKK⎡⎤⎡⎤⎡⎤⎣⎦⎣⎢⎥⎢⎥⎡⎤⎡⎤⎡⎤⎡⎤⎣⎦⎣⎦⎣⎦⎣⎦⎢⎥==⎢⎥⎦⎡⎤⎡⎤⎣⎦⎣⎢⎥⎢⎥⎡⎤⎡⎤⎡⎤⎡⎤⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦∫∫⎦式中，[]11122122233101001100002DDEDDDμμμμ⎡⎤⎢⎥0D⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥−⎢⎥⎢⎥−⎣⎦⎢⎥⎣⎦。令[]()11123321222331010012121100002DDtEtDDDDDμμμμ⎡⎤⎢⎥0⎡⎤⎢⎥⎢⎥⎡⎤===⎢⎥⎢⎥⎣⎦−⎢⎥⎢⎥−⎣⎦⎢⎥⎣⎦1。经积分计算得到薄板单元刚度矩阵各元素如下[]156513262114iiCCCKCCCCCC−⎡⎤⎢⎥=−⎢⎥⎢⎥−−⎣⎦、29891508018ijjiCCCKKCCCC−⎡⎤⎢⎥⎡⎤⎡⎤==⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦、[][]411121119012020immiCCCKKCCCC−⎡⎤⎢⎥==−⎢⎥⎢⎥⎣⎦、371717010016ippiCCCKKCCCC−0⎡⎤⎢⎥⎡⎤⎡⎤==−⎣⎦⎣⎦⎢⎥⎢⎥−⎣⎦、1565132162114jjCCCKCCCCCC⎡⎤⎢⎥⎡⎤=⎣⎦⎢⎥⎢⎥⎣⎦、3710717010016jmmjCCCKKCCCC⎡⎤⎢⎥⎡⎤⎡⎤==−⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦、第5页共16页411121119012020jppjCCCKKCCCC⎡⎤⎢⎥⎡⎤⎡⎤==−⎣⎦⎣⎦⎢⎥⎢⎥−⎣⎦、[]1565132162114mmCCCKCCCCCC−⎡⎤⎢⎥=−−⎢⎥⎢⎥−⎣⎦、29891508018mppmCCCKKCCCC−⎡⎤⎢⎥⎡⎤⎡⎤==−⎣⎦⎣⎦⎢⎥⎢⎥−⎣⎦、155132162114ppCCCKCCCCCC−−6⎡⎤⎢⎥⎡⎤==−⎣⎦⎢⎥⎢⎥−⎣⎦、式中，112212333311.412baCDDDababab=+++D，112212333311.4222baCDDDababab=−+−−D，112212333311.4322baCDDDababab=−−−D，112212333311.44222baCDDDababab=−−++D，221233210.252aCDDbaa=++D，111233210.262bCDDabb=++D，223320.27aCDba=+D，113320.28bCDab=+D，221233210.2922aCDDbaa=−−D，111233210.21022bCDDabb=−−D，223320.2112aCDba=−D，113320.2122bCDab=−D，22334813315abCDba=+D，11334814315baCDab=+D，22332815315abCDba=−D，11332816315baCDab=−D，22332217315abCDba=−D，第6页共16页11332218315baCDab=−D，2233219315abCDba=+D，1133220315baCDab=+D，1221CD=。第7页共16页附录2：函数原程序：主程序%----------------------------------------------------------------%PURPOSE%Analysisofaplate.%----------------------------------------------------------------%-----Topology-------------------------------------------------Enode=[11254;22365;35698;44587];ndof=3;Edof=caldof(Enode,ndof);ex=[0110;1221;1221;0110];ey=[0011;0011;1122;1122];%-----Elementstiffness----------------------------------------E=2.0e11;第8页共16页v=0.3;D=hooke(E,v);t=0.01;ep=[t];nie=size(Enode,1);fori=1:nieKe(:,:,i)=platre(ex(i,:),ey(i,:),ep,D);end%-----AssembleKeintoK---------------------------------------K=zeros(27);K=assem(Edof,K,Ke);%-----loadvectorfandboundaryconditionsbc-----------------f=zeros(27,1);f(13)=-400;bc=[10;20;30;40;50;60;70;80;90;100;110;120;160;170;180;190;200;210;220;230;240;250;260;270];%-----Solvethesystemofequationsandcomputereactions------[a]=solveq(K,f,bc);%-----Postprocess----------------------------------------------Ed=extract(Edof,a);第9页共16页%------------------------end-----------------------------------各子程序如下：（1）计算单元结点自由度编号子程序function[Edof]=caldof(Enode,ndof)%Edof=caldof(Enode,ndof)%-------------------------------------------------------------%PURPOSE%Calculatethenodedegreeoffreedom.%%INPUT:Enode:elementnode%%ndof:degreesoffreedomineachnode%%OUTPUT:Edof:topologyofthestructure%-------------------------------------------------------------Nele=size(Enode,1);Esize=size(Enode,2)-1;fori=1:NeleEdof(i,1)=i;forj=1:EsizeEdof(i,3*j-1)=3*(Enode(i,j+1)-1