第5讲—有限元程序设计

机械与运载工程学院湖南大学CollegeofMechanical&VehicleEngineeringHunanUniversity崔向阳有限元程序设计第5讲:有限单元法崔向阳19:46(0,0)A(0,2)B(4,0)D(4,1)C1F考虑一个平面应力问题如图所示，假设厚度h=1，材料为各项同性，杨氏模量为E=1，泊松比为ν=0，相关力和位移边界条件如图中所示，问题左端为固定约束。试用两个三角形单元分析此问题，三角形单元的网格划分如图所示。试求问题各节点位移u、v和应力σx，σy和σxy。问题描述有限单元法崔向阳19:46三角形单元程序clearallclearallfirst_timefirst_time==cputimecputime;;formatlongformatlong%%----------------------------------------------------------------------------------------%inputdataforcontrolparameters%inputdataforcontrolparameters%%----------------------------------------------------------------------------------------lengthxlengthx=4;%lengthofx=4;%lengthofx--axissideofproblemaxissideofproblemlengthy=2;%lengthofylengthy=2;%lengthofy--axissideofproblemaxissideofproblememoduleemodule=1.0;%elasticmodulus=1.0;%elasticmoduluspoissonpoisson=0.0;%Poisson'sratio=0.0;%Poisson'sratiofloadfload==--1;%thetotalload1;%thetotalload(0,0)A(0,2)B(4,0)D(4,1)C1F有限单元法崔向阳19:46三角形单元程序lx=16;%numberofelementinxlx=16;%numberofelementinx--axisaxislyly=8;%numberofelementiny=8;%numberofelementiny--axisaxisnelnel=2*lx*=2*lx*lyly;%numberofelement;%numberofelementnnelnnel=3;%numberofnodesperelement=3;%numberofnodesperelementndofndof=2;%numberof=2;%numberofdofsdofspernodepernodennodennode=(lx+1)*(ly+1);%totalnumberofnodesinsystem=(lx+1)*(ly+1);%totalnumberofnodesinsystemsdofsdof==nnodennode**ndofndof;%totalsystem;%totalsystemdofsdofsedofedof==nnelnnel**ndofndof;%degreesoffreedomperelement;%degreesoffreedomperelement%%--------------------------------------------------------------------------------------------------------------%inputdatafornodalcoordinatevalues%inputdatafornodalcoordinatevalues%%--------------------------------------------------------------------------------------------------------------x0=[];x0=[];forfori=1:lx+1i=1:lx+1forforj=1:ly+1j=1:ly+1x0=[x0;(ix0=[x0;(i--1)*1)*lengthxlengthx/lx/lx--0.5*lengthy*(1+(lx+10.5*lengthy*(1+(lx+1--i)/lx)*(1i)/lx)*(1--(j(j--1)/ly)];1)/ly)];endendendend有限单元法崔向阳19:46三角形单元程序%%--------------------------------------------------------------------------------------------------------------------------------------------%inputdatafornodalconnectivityforeachelement%inputdatafornodalconnectivityforeachelement%%nodes(i,jnodes(i,j)wherei)wherei--elementno.andjelementno.andj--connectednodesconnectednodes%%--------------------------------------------------------------------------------------------------------------------------------------------nodes=[];nodes=[];forfori=1:lxi=1:lxforforj=1:lyj=1:lynodes=[nodes;(ly+1)*(inodes=[nodes;(ly+1)*(i--1)+j(ly+1)*1)+j(ly+1)*i+ji+j(ly+1)*(i(ly+1)*(i--1)+j+1;];1)+j+1;];nodes=[nodes;(ly+1)*nodes=[nodes;(ly+1)*i+ji+j(ly+1)*i+j+1(ly+1)*(i(ly+1)*i+j+1(ly+1)*(i--1)+j+1;];1)+j+1;];endendendend有限单元法崔向阳19:46三角形单元程序%%--------------------------------------------------------------------%inputdataforboundaryconditions%inputdataforboundaryconditions%%--------------------------------------------------------------------bcdofbcdof=[];=[];bcvalbcval=[];=[];fori=1:ly+1fori=1:ly+1bcdofbcdof=[=[bcdofbcdof1+2*(i1+2*(i--1)2+2*(i1)2+2*(i--1)];1)];bcvalbcval=[=[bcvalbcval00];00];endend(0,0)A(0,2)B(4,0)D(4,1)C有限单元法崔向阳19:46三角形单元程序%%--------------------------------------------------------------------------------------------------%initializationofmatricesandvectors%initializationofmatricesandvectors%%--------------------------------------------------------------------------------------------------ff=sparse(sdof,1);ff=sparse(sdof,1);%systemforcevector%systemforcevectork=k=sparse(edof,edofsparse(edof,edof););%initializationofelementmatrix%initializationofelementmatrixkkkk==sparse(sdof,sdofsparse(sdof,sdof););%systemmatrix%systemmatrixdispdisp=sparse(sdof,1);=sparse(sdof,1);%systemdisplacementvector%systemdisplacementvectoreldispeldisp=sparse(edof,1);=sparse(edof,1);%elementdisplacementvector%elementdisplacementvectorstress=zeros(nel,3);stress=zeros(nel,3);%matrixcontainingstresscomponents%matrixcontainingstresscomponentsstrain=zeros(nel,3);strain=zeros(nel,3);%matrixcontainingstraincomponents%matrixcontainingstraincomponentsindex=sparse(edof,1);index=sparse(edof,1);%indexvector%indexvectorkinmtxkinmtx=sparse(3,edof);=sparse(3,edof);%kinematicmatrix%kinematicmatrixmatmtxmatmtx=sparse(3,3);=sparse(3,3);%constitutivematrix%constitutivematrix%%--------------------------------------------------------------------------------------------------%computematerialmatrices%computematerialmatrices%%--------------------------------------------------------------------------------------------------matmtxmatmtx==fematisofematiso(1,emodule,poisson);%constitutive(1,emodule,poisson);%constitutivematricematrice有限单元法崔向阳19:46三角形单元程序%%--------------------------------------------------------------------------------------------------%computeelementmatricesandvectorsandassemble%computeelementmatricesandvectorsandassemble%%--------------------------------------------------------------------------------------------------forforieliel=1:nel%loopforthetotalnum