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

第4章平面问题三角形常应变单元矩形单元有限元方法的一般讨论基本概念刚架结构被分成若干单元，单元内假设力和位移之间的关系，即单元刚度方程。连续体结构也划分成单元，建立结点力和结点位移之间的关系2§4.1三角形常应变单元jmixyo三角单元结点坐标结点位移向量§4.1.1结点位移3§4.1.2位移模式和形函数假设位移模式jmixyo三角形单元系数矩阵行列式结点位移4克莱姆法则§4.1.2位移模式和形函数位移模式形函数5§4.1.2位移模式和形函数jmixyo三角单元求三角形ijm的面积定义向量按照第一列展开系数矩阵行列6形函数的性质jxyoimp(x,y)三角单元三角形ijm的面积灰色区域面积按照第1行展开7形函数的性质ixyomp(x,y)j三角单元8面积坐标jmip(x,y)92.位移模式定义形函数矩阵10例4.1求形函数o三角单元结点坐标11§4.1.3应变和几何矩阵122几何矩阵几何矩阵13几何子矩阵几何方程例4.2求几何矩阵三角形单元o结点坐标14几何矩阵例4.1结果比较或根据§4.1.4应力和应力矩阵应力应变关系称为应力矩阵弹性矩阵15§4.1.5刚度矩阵和刚度方程方程回顾jmixyo三角形单元16§4.1.5刚度矩阵和刚度方程方程回顾假设结点虚位移U*，导致位移、应变和应力假设平面结构厚度为t，受到体力g

，面力q和集中力p处于平衡状态，按照虚功原理，外力和内力在虚位移上的做功相等：由虚位移的任意性17§4.1.5刚度矩阵和刚度方程等效结点力列向量单元刚度矩阵实际上，由于刚度矩阵中被积函数是常数，可以直接积分出来单元刚度方程可以写成子矩阵的形式18刚度矩阵子矩阵19例4.3求单元刚度矩阵o三角形单元20例4.1-4.3的MATLAB代码syms

N

[3,1]syms

x

y

ui

uj

umxy

=

[1,0;

0,1;

0,0];mu

=

0;A

=

det([ones(3,1),xy])/2;m

=

[1,2,3,1,2];for

i

=

1:3%单元面积式(4.8)N(i)

=

det([ones(3,1),[x,y;xy(m(i+1:i+2),:)]])/2/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)21§4.1.6等效结点力三角形单元o例4.4

求结点力

1.

集中力例:集中力(Fx,Fy)作用在(0.5,0.5)22结点力列阵三角形单元o2.分布面力例:均布压力23结点力列阵三角形单元o3.分布体力例:自重24解题步骤23456yxF划分有限单元；结点编号；单元编号；计算单元刚度矩阵、组装整体刚度矩阵；施加结点力；施加结点位移约束；建立刚度方程；解线性方程组，求解出结点位移；求出单元应力、主应力；求出位移约束结点的约束力。122xy2F2F2225算例对角压方板结点位移结点力1234562y2F2F22

x2123456yxF26算例对角压方板形成结构刚度矩阵5m4i6ji

jm2

m

i

3ji

my1

jx27形成结构刚度矩阵5m4i6i

jm2

m

i

3ji

mjy1

jx28算例对角压方板形成结构刚度矩阵m4ii6jm5j2

m

i

3ji

my1

jx2y2F2F22

x2单元①29算例对角压方板形成结构刚度矩阵m4ii6jm5j2

m

i

3ji

my1

jx2y2F2F22

x2单元②30算例对角压方板形成结构刚度矩阵m4ii6jm5j2

m

i

3ji

my1

jx2y2F2F2312

x2形成结构总刚度方程5m4i6i

jmyjm

i

3ji

mjx32形成结构总刚度方程5m4i6i

jm2

m

i

3ji

mjy1

jx33形成结构总刚度矩阵单元1234i3256j1523m23455m4i6i

jm2

m

i

3ji

mjy1

jx单元①45612345634123形成结构总刚度矩阵单元1234i3256j1523m23455m4i6i

jm2

m

i

3ji

mjy1

jx单元②45612345635123形成结构总刚度矩阵单元1234i3256j1523m23455m4i6i

jm2

m

i

3ji

mjy1

jx12345612345636形成结构总刚度矩阵单元1234i3256j1523m23455m4i6i

jm2

m

i

3ji

mjy1

jx12345612345637形成结构总刚度矩阵单元1234i3256j1523m23455m4i6i

jm2

m

i

3ji

mjy1

jx12345612345638形成结构总刚度矩阵39形成结构总刚度方程40求结点位移由另外6个方程可以求解出约束力41单元应力5m4i6i

jm2

m

i

3ji

mjy1

jx42其他形式的三角形单元4344§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);自由度%板厚*%弹性模量*%泊松比*%结点坐标*%单元信息*%结点数%单元数%结点力*%位移约束自%非位移约束5m4i6i

jm2

m

i

3ji

mjy1

jx结点号45自由度序数结点位移向量结点i的位移在2i-1:2i位置结点号与自由度的对应关系如果单元的两个结点是i和j，它的位移在[2i-1:2i

2j-1:2j]位置应变矩阵function

[StrainM,A]=PlanTriaStrain(xy)元应变矩阵A=0.5*det([ones(3,1),xy]);单元面积n=[1,2,3,1,2];for

i=1:3b=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)=[b

0;0

c;c

b]/A/2;几何矩阵end%%%46形成总刚和求解D=Em/(1-Pr^2)*[1,Pr,0;Pr,1,0;0,0,(1-Pr)/2];%弹性矩阵K=sparse(2*nd,2*nd);for

el=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;%刚度矩阵end47形成总刚和求解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’)标题for

j=1:nd输出结点号，结点坐标，结点位移%fprintf("%4i%10.4f%10.4f%14.4g%14.4g

\n

",j,gxy(j,:),U(2*j+(-1:0))))u48enNdode123456X0.00000.00001.00000.00001.00002.0000Y2.00001.00001.00000.00000.00000.0000001.261e-01002.677e-0103.192e-010v-1.561e-009-6.569e-010-1.717e-010000结果输出fprintf("%4s%4s%4s%4s%12s%14s%14s%14s%14s%9s\n","Elem","i","j","k","Sx","Sy","S%S2","An")标题for

el=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));输出单元号，单元信息，应力和主应力end4950单元信息及应力Elemi

j

k

SxSySxyS1S2

An1

312-33.52-20039.19-24.76-208.825.212

25317.21-30.8927.7529.88-43.5649.093

52416.31-133.1016.31-133.10.004

6350-36.05-11.443.324-39.38-32.40§4.2矩形单元§4.2.1结点位移和结点力结点位移向量oaabb矩形单元51§4.2矩形单元§4.2.2形函数和位移模式1o11

1母单元aabob矩形单元求解4个线性方程得到系数52§4.2矩形单元位移模式§4.2.2形函数和位移模式形函数oaabb矩形单元53形函数54几何矩阵§4.2.2应变和几何矩阵应变55§4.2.3应力和应力矩阵应力应力矩阵56§4.2.4单元刚度矩阵57§4.2.4单元刚度矩阵5859矩形单元MATLAB代码syms

xy[3,2]

real

syms

x

y

ym

A

E

mu

t

realm

=

[1,2,3,1,2];for

i

=

1:3N0(i)

=

simplify(det([ones(3,1),[x,y;xy(m(i+1:i+2),:)]])/2/A);enddisp(N0)Nm

=

kron(N0,eye(2));%形函数(4.14disp(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)

%刚度矩阵§4.2.5程序设计单元划分及结点编号60算例结点编号（绿字）单元编号（蓝字）xy10MPa10MPa10MPa10MPaxy564③3②2①1④⑤716212223242581514131361结点力xy10MPa10MPa10MPa10MPaMPa62Em=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);for

i=1:Nxfor

j=1:Nyndel(Ny*(i-1)+j,:)=(Ny+1)*(i-1)+j+[Ny+2,1,0,Ny+1];%单元信息endend63§4.2.5程序设计输入结构参数64形成结点力和刚度矩阵x=[1,-1,-1,1];

y=[1,1,-1,-1];结点局部坐标K=sparse(2*nd,2*nd);for

el=1:nefor

i=1:4for

j=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;

%%%65求解与结果输出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")%标题%for

i=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""Sxy’,

"S1","S2","An")输出标题D=Em/(1-Pr^2)*[1,Pr,0;Pr,1,0;0,0,(1-Pr)/2];%弹性矩阵for

el=1:nefor

i=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);67结点位移输出Node...uv210-1.697e-006225.701e-007-1.711e-006231.138e-006-1.751e-006241.703e-006-1.818e-006252.247e-006-1.914e-006Elem

i

jl

mSxSySxyS1S2An1

72161.258e+004-15.966.11.258e+004-15.980.152

83273.774e+004-29.62-0.096123.774e+004-29.62-0.003

94386.288e+004-19.46-10.366.288e+004-19.46-0.024105498.803e+004-4.462-5.6418.803e+004-4.463-0.01...§4.3有限元方法的一般讨论§4.3.1位移模式协调性：位移在单元内以及单