《有限元方法与MATLAB程序设计》-课件-4-平面问题

第4章平面问题4.1三角形常应变单元4.2矩形单元4.3有限元方法的一般讨论2基本概念刚架结构被分成若干单元，单元内假设力和位移之间的关系，即单元刚度方程。连续体结构也划分成单元，建立结点力和结点位移之间的关系3§4.1三角形常应变单元jmixyo三角单元结点坐标结点位移向量§4.1.1结点位移4§4.1.2位移模式和形函数假设位移模式jmixyo三角形单元结点位移克莱姆法则系数矩阵行列式5§4.1.2位移模式和形函数位移模式形函数6§4.1.2位移模式和形函数求三角形ijm的面积jmixyo三角单元按照第一列展开定义向量系数矩阵行列7形函数的性质jmixyo三角单元p(x,y)三角形ijm的面积灰色区域面积按照第1行展开8形函数的性质jmixyo三角单元p(x,y)9

面积坐标jmip(x,y)102.位移模式定义形函数矩阵11例4.1求形函数三角单元o结点坐标12§4.1.3应变和几何矩阵132几何矩阵几何矩阵几何子矩阵几何方程14例4.2求几何矩阵三角形单元o结点坐标几何矩阵例4.1结果比较或根据15§4.1.4应力和应力矩阵应力应变关系称为应力矩阵弹性矩阵16§4.1.5刚度矩阵和刚度方程方程回顾jmixyo三角形单元17§4.1.5刚度矩阵和刚度方程方程回顾假设结点虚位移U*，导致位移、应变和应力假设平面结构厚度为t，受到体力g

，面力q和集中力p处于平衡状态，按照虚功原理，外力和内力在虚位移上的做功相等：由虚位移的任意性18§4.1.5刚度矩阵和刚度方程等效结点力列向量单元刚度矩阵单元刚度方程实际上，由于刚度矩阵中被积函数是常数，可以直接积分出来可以写成子矩阵的形式19刚度矩阵子矩阵20例4.3求单元刚度矩阵三角形单元o例4.1-4.3的MATLAB代码21symsN[3,1]symsxyuiujumxy=[1,0;0,1;0,0];mu=0;A=det([ones(3,1),xy])/2; %单元面积式(4.8)m=[1,2,3,1,2];fori=1:3N(i)=det([ones(3,1),[x,y;xy(m(i+1:i+2),:)]])/2/A;enddisp([ui,uj,um]*N) %位移模式公式参考此处Nx=diff(N,x)';Ny=diff(N,y)';B=[kron(Nx,[1,0]);kron(Ny,[0,1]);kron(Ny,[1,0])+kron(Nx,[0,1])];disp(B)D=1/(1-mu^2)*[1,mu,0;mu,1,0;0,0,(1-mu)/2];ke=A*B'*D*B; %刚度矩阵式(4.46)disp(ke)22§4.1.6等效结点力三角形单元o1.集中力例:集中力(Fx,Fy)作用在(0.5,0.5)例4.4求结点力23结点力列阵三角形单元o2.分布面力例:均布压力24结点力列阵三角形单元o3.分布体力例:自重25解题步骤123456yxF划分有限单元；结点编号；单元编号；计算单元刚度矩阵、组装整体刚度矩阵；施加结点力；施加结点位移约束；建立刚度方程；解线性方程组，求解出结点位移；求出单元应力、主应力；求出位移约束结点的约束力。2F2F22xy2226算例对角压方板结点位移123456结点力2F2F22xy22123456yxF27算例对角压方板形成结构刚度矩阵123456jmijmiijmijmyx28形成结构刚度矩阵123456jmijmiijmijmyx29算例对角压方板形成结构刚度矩阵123456jmijmiijmijmyx2F2F22xy22单元①30算例对角压方板形成结构刚度矩阵123456jmijmiijmijmyx2F2F22xy22单元②31算例对角压方板形成结构刚度矩阵123456jmijmiijmijmyx2F2F22xy2232形成结构总刚度方程123456jmijmiijmijmyx33形成结构总刚度方程123456jmijmiijmijmyx34形成结构总刚度矩阵

单元

1234i3256j1523m2345123456jmijmiijmijmyx单元①12345612345635形成结构总刚度矩阵

单元

1234i3256j1523m2345123456jmijmiijmijmyx单元②12345612345636形成结构总刚度矩阵

单元

1234i3256j1523m2345123456jmijmiijmijmyx12345612345637形成结构总刚度矩阵

单元

1234i3256j1523m2345123456jmijmiijmijmyx12345612345638形成结构总刚度矩阵

单元

1234i3256j1523m2345123456jmijmiijmijmyx12345612345639形成结构总刚度矩阵40形成结构总刚度方程41求结点位移由另外6个方程可以求解出约束力42单元应力123456jmijmiijmijmyx43其他形式的三角形单元44§4.1.7程序设计Th=1e-2; %板厚*Em=210e9; %弹性模量*Pr=0.3; %泊松比*gxy=[0,2;0,1;1,1;0,0;1,0;2,0]; %结点坐标*ndel=[3,1,2;2,5,3;5,2,4;6,3,5]; %单元信息*nd=size(gxy,1); %结点数ne=size(ndel,1); %单元数F=zeros(2*nd,1);F(2)=-1; %结点力*dofix=[1,3,7,8,10,12]; %位移约束自由度dofree=setdiff(1:2*nd,dofix); %非位移约束自由度123456jmijmiijmijmyx45结点号自由度序数结点位移向量结点i的位移在2i-1:2i位置结点号与自由度的对应关系如果单元的两个结点是i和j，它的位移在[2i-1:2i2j-1:2j]位置46应变矩阵function[StrainM,A]=PlanTriaStrain(xy) %元应变矩阵A=0.5*det([ones(3,1),xy]); %单元面积n=[1,2,3,1,2];fori=1:3

b=xy(n(i+1),2)-xy(n(i+2),2);c=xy(n(i+2),1)-xy(n(i+1),1);StrainM(:,2*i-1:2*i)=[b0;0c;cb]/A/2; %几何矩阵end47形成总刚和求解D=Em/(1-Pr^2)*[1,Pr,0;Pr,1,0;0,0,(1-Pr)/2]; %弹性矩阵K=sparse(2*nd,2*nd);forel=1:neN(2:2:6)=2*ndel(el,:);N(1:2:5)=N(2:2:6)-1; %单元自由度

[B,A]=PlanTriaStrain(gxy(ndel(el,:),:)); %几何矩阵和单元面积

K(N,N)=K(N,N)+B'*D*B*Th*A; %刚度矩阵end48形成总刚和求解U=zeros(2*nd,1);U(dofree)=K(dofree,dofree)\(F(dofree)-K(dofree,dofix)*U(dofix));fprintf('%4s%8s%10s%12s%14s\n','Node','X','Y','u','v’) %标题forj=1:nd %输出结点号，结点坐标，结点位移

fprintf('%4i%10.4f%10.4f%14.4g%14.4g\n',j,gxy(j,:),U(2*j+(-1:0))))endNodeXYuv10.00002.00000-1.561e-00920.00001.00000-6.569e-01031.00001.00001.261e-010-1.717e-01040.00000.00000051.00000.00002.677e-010062.00000.00003.192e-010049结果输出fprintf('%4s%4s%4s%4s%12s%14s%14s%14s%14s%9s\n','Elem','i','j','k','Sx','Sy','Sxy','S1','S2','An') %标题forel=1:neN(2:2:6)=2*ndel(el,:);N(1:2:5)=N(2:2:6)-1; %单元自由度

[B,A]=PlanTriaStrain(gxy(ndel(el,:),:)); %几何矩阵和单元面积

S=D*B*U(N);

c1=(S(1)+S(2))/2;c2=(S(1)-S(2))/2;c3=sqrt(c2^2+S(3)^2);fprintf('%4i%4i%4i%4i%14.4g%14.4g%14.4g%14.4g%14.4g%7.2f\n',el,ndel(el,:),S,c1+c3,c1-c3,atan2(S(3),c2)); %输出单元号，单元信息，应力和主应力end50单元信息及应力ElemijkSxSySxyS1S2An

1312-33.52-20039.19-24.76-208.825.21225317.21-30.8927.7529.88-43.5649.09352416.31-133.1016.31-133.10.0046350-36.05-11.443.324-39.38-32.4051§4.2矩形单元§4.2.1结点位移和结点力结点位移向量oaabb矩形单元52§4.2矩形单元§4.2.2形函数和位移模式母单元o1111oaabb矩形单元求解4个线性方程得到系数53§4.2矩形单元形函数位移模式§4.2.2形函数和位移模式oaabb矩形单元形函数-101-10100.20.40.60.815455§4.2.2应变和几何矩阵

几何矩阵应变56§4.2.3应力和应力矩阵应力矩阵应力57§4.2.4单元刚度矩阵58§4.2.4单元刚度矩阵矩形单元MATLAB代码59symsxy[3,2]realsymsxyymAEmutrealm=[1,2,3,1,2];fori=1:3N0(i)=simplify(det([ones(3,1),[x,y;xy(m(i+1:i+2),:)]])/2/A); %形函数(4.14)enddisp(N0)Nm=kron(N0,eye(2));disp(Nm) %形函数矩阵(4.17)Nx=diff(N0,x);Ny=diff(N0,y);B0=[kron(Nx,[1,0]);kron(Ny,[0,1]);kron(Ny,[1,0])+kron(Nx,[0,1])];disp(B0) %几何矩阵(4.31)D=E/(1-mu^2)*[1,mu,0;mu,1,0;0,0,(1-mu)/2];k0=A*t*B0'*D*B0;disp(k0) %刚度矩阵60§4.2.5程序设计单元划分及结点编号61算例结点编号（绿字）

单元编号（蓝字）xy10MPa10MPa10MPa10MPaxy123456①②③④⑤716212223242581514131362结点力xy10MPa10MPa10MPa10MPaMPa63§4.2.5程序设计Em=210e9;Pr=0.3;Th=1e-2; %弹性模量，泊松比，板厚度Nx=4;Lx=4.5;a=Lx/Nx/2; %水平方向单元数量，长度和单元半长Ny=4;Ly=3;b=Ly/Ny/2; %竖直方向单元数量，长度和单元半长ne=Nx*Ny; %单元总数nd=(Nx+1)*(Ny+1); %结点总数ndel=zeros(ne,4); fori=1:Nx forj=1:Nyndel(Ny*(i-1)+j,:)=(Ny+1)*(i-1)+j+[Ny+2,1,0,Ny+1]; %单元信息

endendF=zeros(2*nd,1); F(2*(Ny+1)*Nx+1:2:2*nd)=[0,6,12,18,11]*3e3/94; %结点力列向量dofix=[1:2:9,2,11:10:41];dofree=setdiff(1:2*nd,dofix);输入结构参数64形成结点力和刚度矩阵x=[1,-1,-1,1];y=[1,1,-1,-1]; %结点局部坐标K=sparse(2*nd,2*nd); forel=1:nefori=1:4 forj=1:4c1=b/a*x(i)*x(j)*(1+y(i)*y(j)/3);c3=x(i)*y(j);

c2=a/b*y(i)*y(j)*(1+x(i)*x(j)/3);c4=x(j)*y(i);c0=(1-Pr)/2;ke(2*i+(-1:0),2*j+(-1:0))=[c1+c0*c2,Pr*c3+c0*c4;Pr*c4+c0*c3,c2+c0*c1];

endendN(2:2:8)=2*ndel(el,:);N(1:2:7)=N(2:2:8)-1; %单元自由度

K(N,N)=K(N,N)+Em*Th/4/(1-Pr^2)*ke; %安装总体刚度矩阵end65求解与结果输出U=zeros(2*nd,1); %结点位移列向量U(dofix)=0; %约束位移（如果不是固定填实际具体值）U(dofree)=K(dofree,dofree)\(F(dofree)-K(dofree,dofix)*U(dofix));fprintf('%4s%8s%10s%12s%14s\n','Node','u','v')%标题fori=1:nd %输出结点号，结点位移

fprintf('%4i%12.4g%12.4g\n',i,U(2*i+[-1;0]))end66求解与结果输出fprintf('%4s%4s%4s%4s%4s%12s%14s%14s%14s%14s%9s\n','Elem','i','j','l','m','Sx','Sy’,'Sxy’,'S1','S2','An') %输出标题D=Em/(1-Pr^2)*[1,Pr,0;Pr,1,0;0,0,(1-Pr)/2];%弹性矩阵forel=1:nefori=1:4B(:,2*i-1:2*i)=[x(i)/a,0;0,y(i)/b;y(i)/b,x(i)/a]/4; %几何矩阵

endN(2:2:8)=2*ndel(el,:);N(1:2:7)=N(2:2:8)-1; %单元自由度

S=D*B*U(N); %应力列向量

c1=0.5*(S(1)+S(2));c2=0.5*(S(1)-S(2));c3=sqrt(c2^2+S(3)^2);fprintf('%4i%4i%4i%4i%4i%14.4g%14.4g%14.4g%14.4g%14.4g%7.2f\n‘,el,ndel(el,:),S,c1+c3,c1-c3,atan2(S(3),c2));end %输出单元号，单元信息，应力列向量，主应力67

结点位移输出Nodeuv...210-1.697e-006225.701e-007-1.711e-006231.138e-006-1.751e-006241.703e-006-1.818e-006252.247e-006-1.914e-006ElemijlmSxSySxyS1S2An172161.258e+004-15.966.11.258e+004-15.980.15283273.774e+004-29.62-0.096123.774e+004-29.62-0.00394386.288e+004-19.46-10.366.288e+004-19.46-0.024105498.803e+004-4.462-5.6418.803e+00