有限单元法大作业

有限单元法大作业问题：一薄板，尺寸为10x10,厚为0.1,中心受集中力为400,四个角点为简支，求薄板的位移场。分析：根据题目要求，可以选用四节点板单元。对题目模型分网如下：图1：分网图单元大小为2.5X2.5,由于单元仅在中心处受集中载荷，结构是对称的，因此在计算时只需对其四分之一进行计算，然后再将计算结果扩充到整个板就可以了。源程序如下：%用四节点平面板单元求解平板中心受集中载荷问题%输入控制参数%clear%清除workplace残留数据nel=4;%单元数nnel=4;%每个单元的节点数ndof=3;%每个节点的自由度数nnode=9;%系统总节点数sdof=nnode*ndof;%系统总自由度数edof=nnel*ndof;%每个单元的自由度数emodule=206e9;%杨氏弹性模量poisson=0.3;%泊松比t=0.1;%薄板厚度nglxb=2;nglyb=2;%弯曲对应的2x2高斯拉格朗日积分nglb=nglxb*nglyb;%弯曲对应的每个单元的高斯积分点数nglxs=l;nglys=l;%剪切对应的lxl高斯拉格朗日积分ngls=nglxs*nglys;%剪切对应的每个单元的高斯积分点数%输入节点坐标值%gcoord(i,j)i：节点号j：x,y值%gcoord=0.00.0;2.50.0;5.00.0;0.02.5;2.52.5;5.02.5;0.05.0;2.55.0;5.05.0;%每个单元对应的节点号(逆时针排列)%nodes(i,j)i:节点号j:对应的单元号%nodes=l254;2365;4587;5698;%输入边界条件%bcdof=l234679lll2l6202l232526;%约束的自由度bcval=zeros(l,l5);%对应的值%初始化矩阵和矢量%ff=zeros(sdof,l);%kk=zeros(sdof,sdof);%disp=zeros(sdof,l);%index=zeros(edof,l);%kinmtpb=zeros(3,edof);%matmtpb=zeros(3,3);%kinmtps=zeros(2,edof);%matmtps=zeros(2,2);%载荷矢量系统刚度矩阵系统位移矢量每个单元所包含的自由度弯曲几何函数矩阵弯曲材料系数矩阵剪切几何函数矩阵剪切材料系数矩阵%载荷矢量%单元刚度矩阵计算及其组合%弯曲相关计算pointb,weightb=swp2(nglxb,nglyb);%积分点和权系数matmtpb=sbm(emodule,poisson)*t3/12;%弯曲材料系数%积分点和权系数%剪切模量%剪切修正因数%剪切材料系数矩阵%剪切相关计算points,weights=swp2(nglxs,nglys);shearm=0.5*emodule/(1.0+poisson);shcof=5/6;matmtps=shearm*shcof*t*10;01;foriel=1:nel%对所有单元数的循环fori=1:nnelnd(i)=nodes(iel,i);%当前单元对应的节点初始化单元刚度矩阵初始化弯曲刚度矩阵初始化剪切刚度矩阵xcoord(i)=gcoord(nd(i),1);%节点对应的x坐标值ycoord(i)=gcoord(nd(i),2);%节点对应的y坐标值endk=zeros(edof,edof);%kb=zeros(edof,edof);%ks=zeros(edof,edof);%弯曲相关计算%forintx=1:nglxbx=pointb(intx,1);wtx=weightb(intx,1);forinty=1:nglyby=pointb(inty,2);wty=weightb(inty,2);%x轴积分点坐标%权系数%y轴积分点坐标%权系数shape,dhdr,dhds=ssf(x,y);%计算形函数和对其相应的求导jacob2=sjc(nnel,dhdr,dhds,xcoord,ycoord);%计算雅可比行列式detjacob=det(jacob2);invjacob=inv(jacob2);%计算雅可比行列式的值%求雅可比行列式的逆dhdx,dhdy=sxy(nnel,dhdr,dhds,invjacob);%计算ddhdr,dhds在迪卡尔坐标下的值kinmtpb=sbB(nnel,dhdx,dhdy);%计算弯曲几何函数矩阵%计算弯曲刚度矩阵%kb=kb+kinmtpb*matmtpb*kinmtpb*wtx*wty*detjacob;endend%结束弯曲刚度矩阵的计算%剪切相关计算%forintx=1:nglxsx=points(intx,1);wtx=weights(intx,1);forinty=1:nglysy=points(inty,2);wty=weights(inty,2);%x轴积分点坐标%权系数%y轴积分点坐标%权系数shape,dhdr,dhds=ssf(x,y);%计算形函数和对其相应的教学求导jacob2=sjc(nnel,dhdr,dhds,xcoord,ycoord);%计算雅可比行列式detjacob=det(jacob2);%计算雅可比行列式的值invjacob=inv(jacob2);%求雅可比行列式的逆dhdx,dhdy=sxy(nnel,dhdr,dhds,invjacob);%计算dhdr,dhds在迪卡尔坐标下的值kinmtps=ssB(nnel,dhdx,dhdy,shape);%计算剪切几何函数矩阵%计算剪切刚度矩阵%ks=ks+kinmtps*matmtps*kinmtps*wtx*wty*detjacob;endend%结束剪切刚度矩阵的计算%计算单元刚度矩阵%k=kb+ks;index=etsd(nd,nnel,ndof);%单元对应的系统自由度号kk=ask(kk,k,index);%合成系统刚度矩阵end%加载边界条件%kk,ff=dbc(kk,ff,bcdof,bcval);%求解%disp=kkff;%输出节点位移%初始化%将输出的结果从5x5的四分之一板扩充到10x10num=1:1:sdof;nodedisp=numdisp%后处理%result=zeros(25,3);displace(75)=0;的全板a=re(1,0,disp);a(75)=0;displace=displace+a;a=re(10,6,disp);a(75)=0;displace=displace+a;a=re(19,12,disp);a(75)=0;displace=displace+a;fori=1:15;displace(45+i)=displace(15+i);displace(60+i)=displace(i);endresult=dtxy(displace);%将节点位移以节点顺序输出，i-节点号，j-对应位移exgcoord=excoord;%将节点坐标从5x5的四分之一板扩充到10x10的全板finresult=agr(exgcoord,result);%将全板的节点坐标和节点位移对应起来，i-节点号，j-对应坐标和位移fori=1:25;finresultin(1,i)=sqrt(finresult(3,i)A2+finresult(4,i)A2+finresult(5,i)A2);%求节点位移endZ=arrayfin(finresultin);%排列节点位移X,Y=meshgrid(0:2.5:10,0:2.5:10);surf(X,Y,Z);%画板变形图以下为子程序：functionpoint2,weight2=swp2(nglx,ngly)%目的:%求二维高斯积分的积分点和权系数%变量：%nglx-x轴高斯积分点数%ngly-y轴高斯积分点数%point2-高斯积分点坐标%weight2-权系数%确定x,y轴最大的积分点数ifnglxnglyngl=nglx;elsengl=ngly;end%初始化point2=zeros(ngl,2);weight2=zeros(ngl,2);%求出相应的积分点坐标和权系数pointx,weightx=swp1(nglx);pointy,weighty=swp1(ngly);%二维积分forintx=1:nglxpoint2(intx,1)=pointx(intx);weight2(intx,1)=weightx(intx);endforinty=1:nglypoint2(inty,2)=pointy(inty);weight2(inty,2)=weighty(inty);endfunctionmatmtrx=sbm(elastic,poisson)%目的:%弯曲材料系数%变量:%elastic-弹性模量%poisson-泊松比%matmtrx=elastic/(1-poisson*poisson)*.1poisson0;.poisson10;.00(1-poisson)/2;functionshapeq4,dhdrq4,dhdsq4=ssf(rvalue,svalue)%目的：%在自然坐标下计算形函数和对其相应的求导%变量：%shapeq4-四节点单元形函数%dhdrq4-形函数对r求导%dhdsq4-形函数对s求导%rvalue-对应点的r坐标值%svalue-对应点的s坐标值%说明：%第一个点自然坐标(-1,-1),第二个点自然坐标(1,-1)%第三个点自然坐标(1,1),第四个点自然坐标(-1,1)%形函数shapeq4(1)=0.25*(1-rvalue)*(1-svalue);shapeq4(2)=0.25*(1+rvalue)*(1-svalue);shapeq4(3)=0.25*(1+rvalue)*(1+svalue);shapeq4(4)=0.25*(1-rvalue)*(1+svalue);%相应的求导dhdrq4(1)=-0.25*(1-svalue);dhdrq4(2)=0.25*(1-svalue);dhdrq4(3)=0.25*(1+svalue);dhdrq4(4)=-0.25*(1+svalue);dhdsq4(1)=-0.25*(1-rvalue);dhdsq4(2)=-0.25*(1+rvalue);dhdsq4(3)=0.25*(1+rvalue);dhdsq4(4)=0.25*(1-rvalue);functionjacob2=sjc(nnel,dhdr,dhds,xcoord,ycoord)%目的：%求二维雅可比行列式%变量：%jacob2-二维雅可比行列式%nnel-每个单元节点数%dhdr-自然坐标r对形函数的求导%dhds-自然坐标s对形函数的求导%xcoord-节点的x坐标值%ycoord-节点的y坐标值%-jacob2=zeros(2,2);fori=1:nneljacob2(1,1)=jacob2(1,1)+dhdr(i)*xcoord(i);jacob2(1,2)=jacob2(1,2)+dhdr(i)*ycoord(i);jacob2(2,1)=jacob2(2,1)+dhds(i)*xcoord(i);jacob2(2,2)=jacob2(2,2)+dhds(i)*ycoord(i);endfunctiondhdx,dhdy=sxy2(nnel,dhdr,dhds,invjacob)%目的：%求dhdr,dhds在迪卡尔坐标下的值%变量：%dhdx-形函数对迪卡尔坐标x的求导%dhdy-形函数对迪卡尔坐标y的求导%nnel-每个单元节点数%dhdr-形函数对自然坐标r的求导%dhds-形函数对自然坐标s的求导%invjacob-求雅可比行列式的逆%-fori=1:nneldhdx(i)=invjacob(1,1)*dhdr(i)+invjacob(1,2)*dhds(i);dhdy(i)=invjacob(2,1)*dhdr(i)+invjacob(2,2)*dhds(i);endfunctionkinmtpb=sbB(nnel,dhdx,dhdy)%目的:%计算弯曲几何函数矩阵%变量:%nnel-每个单元节点数%dhdx-形函数对迪卡尔坐标X的求导%dhdy-形函数对迪卡尔坐标y的求导%-fori=1:nneli1=(i-1)*3+1;i2=i1+1;i3=i2+1;kinmtpb(1,i1)=dhdx(i);kinmtpb(2,i2)=dhdy(i);kinmtpb(3,i1)=dhdy(i);kinmtpb(3,i2)=dhdx(i);kinmtpb(3,i3)=0;endfunctionkinmtps=ssB(nnel,dhdx,dhdy,shape)%目的：%计算剪切几何函数矩阵%变量：%nnel-每个单元节点数%dhdx-形函数对迪卡尔坐标X的求导%dhdy-形函数对迪卡尔坐标y的求导%shape-形函数%-fori=1:nneli1=(i-1)*3+1;i2=i1+1;i3=i2+1;kinmtps(1,i1)=-shape(i);kinmtps(1,i3)=dhdx(i);kinmtps(2,i2)=-shape(i);kinmtps(2,i3)=dhdy(i);endfunctionindex=etsd(nd,nnel,ndof)%目的:%单元对应的系统自由度号%变量：%index-和单元”iel”对应的系统自由度数%iel-要确定系统自由度的单元%nnel-每个单元的节点数%ndof-每个节点的自由度%edof=nnel*ndof;k=0;fori=1:nnelstart=(nd(i)-1)*ndof;forj=1:ndofk=k+1;index(k)=start+j;endendfunctionkk=ask(kk,k,index)%目的:%合成系统刚度矩阵%变量:%kk-系统刚度矩阵%k-单元刚度矩阵%index-d.o.f.vectorassociatedwithanelement%edof=length(index);fori=1:edofii=index(i);forj=1:edofjj=index(j);kk(ii,jj)=kk(ii,jj)+k(i,j);endendfunctionkk,ff=dbc2(kk,ff,bcdof,bcval)%Purpose:%为方程kkx=ff确定边界条件%变量:%kk-系统刚度矩阵%ff-载荷%bcdof-边界节点对应的自由度%bcval-边界节点对应的自由度的值n=length(bcdof);sdof=size(kk);fori=1:nc=bcdof(i);forj=1:sdofkk(c,j)=0;endkk(c,c)=1;ff(c)=bcval(i);endfunctionu=re(a,b,disp)%目的：%从5x5的四分之一板扩充输出结果到10x10的全板%变量：%disp-扩充前的节点位移矩阵%u-扩充后的节点位移矩阵%-fori=a:1:a+8;u(1,i+b)=disp(i);endfori=a+9:1:a+11;u(1,i+b)=disp(i-6);endfori=a+12:1:a+14;u(1,i+b)=disp(i-12);endfunctionresult=dtxy(displace)%目的：%将节点位移以节点顺序输出，i-节点号，j-对应位移%变量：%result-按节点顺序输出后的矩阵%displace-按自由度顺序输出的矩阵%fori=1:25;ifi=1;result(1,1)=displace(i);elseresult(i,1)=displace(i*3-2);endifi=1;result(1,2)=displace(2);elseresult(i,2)=displace(i*3-1);endifi=1;result(1,3)=displace(3);elseresult(i,3)=displace(3*i);endendfunctionexgcoord=excoord%目的：%将节点坐标从5x5的四分之一板扩充到10x10的全板%变量：%exgcoord-扩充后的节点坐标矩阵%x=0:2.5:10;y=0:2.5:10;exgcoord=zeros(25,2);fori=1:5;exgcoord(i,1)=x(i);exgcoord(i,2)=y(1);endfori=6:10exgcoord(i,1)=x(i-5);exgcoord(i,2)=y(2);endfori=11:15exgcoord(i,1)=x(i-10);exgcoord(i,2)=y(3);endfori=16:20exgcoord(i,1)=x(i-15);exgcoord(i,2)=y(4);endfori=21:25exgcoord(i,1)=x(i-20);exgcoord(i,2)=y(5);endfunctionfinresult=agr(exgcoord,result)%目的：%将全板的节点坐标和节点位移对应起来，i-节点号，j-对应坐标和位移%变量：%finresult-节点坐标和对应位移对应起来的矩阵%result-节点位移矩阵%exgcoord-节点坐标矩阵%finresult=zeros(5,25);fori=1:5;ifi=2;finresult(i,:)=exgcoord(:,i);elsefinresult(i,:)=result(:,i-2);endend%目的：%排列节点位移%变量：%finresultin-节点坐标和对应位移对应起来的矩阵%Z-按节点顺序排列的节点位移矩阵%fori=1:25;ifi5&i10&i15&i20&i=25;Z(5,i-20)=finresultin(i);endend计算结果输出如下。1至9节点的位移为：节点号自由度号1位移值1.379E-17128.807E-1831.639E-1743.668E-18252.072E-0567.656E-187-2.392E-17382.569E-059-1.954E-17102.072E-05411-1.448E-1812-4.538E-18131.559E-055141.559E-05154.542E-05161.250E-206172.620E-05186.487E-05192.569E-057200.000E+00210.000E+00222.620E-058230.000E+00246.487E-05250.000E+009260.000E+00279.778E-05由表可知，节点9，也就是薄板中心集中载荷加载处的位移最大。将输出结果扩充到全板，并作图如下：-4kmo