第8章板的彎曲§8.1四結點矩形薄板單元§8.2三結點三角形薄板單元§8.3考慮橫向剪切影響的平板彎曲單元薄板小撓度彎曲基本理論兩個平行平面和垂直于這兩個平行平面的柱面或棱柱面所圍成的物體，稱為平板，或簡稱板。上下兩個表面稱為板面；上下兩個表面之間的距離稱為板厚h；

上下兩個表面的中間平面稱為板的中面；柱面稱為側面或板邊；板厚h遠小于中面的最小尺寸的板稱為薄板，否則為厚板；作用在中面內的荷載為縱向荷載；垂直于中面的荷載為橫向荷載；薄板小撓度彎曲基本理論zdxxxdxxzzdx3物理方程4彎矩5彎矩6板的內力7板的變形應變曲率和扭率xzzuyzzv8應力內力9§8.1四結點矩形薄板單元結點位移結點力oaabb111o1§8.1.1結點力和結點位移w,

Fzyxzo10§8.1.2位移模式和形函數aabob111o111§8.1.2位移模式和形函數形函數w,

Fzyxzo幾何方程12形函數矩陣形函數性質13位移的連續性在單元邊界上3次方程的4個系數由4個參數撓度和單元邊界方向的導數連續完全確定3次方程的4個系數不能由2個參數完全確定單元邊界法線方向的導數不連續四結點矩形單元不是協調元111o11415§8.1.3應變和幾何矩陣幾何矩陣幾何方程幾何子矩陣§8.1.3應變和幾何矩陣16§8.1.3應變和幾何矩陣17§8.1.4應力和內力應力18§8.1.4應力和內力內力yxzMxMyMxyMxyQy19Qx§8.1.4剛度矩陣20§8.1.5等效結點力平板在分布集度為q的橫向力的作用下，等效結點力為，特別地，當q=q0為常數時，上式積分為2122xy2Lx2Lyo彈性模量，泊松比和%x向單元數，長度和%y向單元數，長度和%單元總數和結點總%自動形成單元信息ndel(Ny*(i-1)+j,:)=(Ny+1)*(i-1)+j+[Ny+2,1,0,Ny+1];endendF=zeros(3*nd,1);F(1)=0.25;%結點力列向量%置結點力dofix=[3*Ny+1:3*Ny+3:3*nd,3*Nx*(Ny+1)+1:3:3*nd-3];

%簡支dofix=union(dofix,[2:3*(Ny+1):3*nd,3:3:3*(Ny+1)]);Em§=281.01e.96;P程r=序0.設3;計Th=1e-2;定%義結構板厚Nx=4;Lx=4;a=Lx/Nx/2;單元半長Ny=4;Ly=4;b=Ly/Ny/2;單元半長ne=Nx*Ny;nd=(Nx+1)*(Ny+1);數ndel=zeros(ne,4);for

i=1:Nxfor

j=1:Nyy51910x7

12

17①

⑤②8

13

18③④1625234611

162115

20192423221423形成剛度矩陣和求解K=sparse(3*nd,3*nd);%結構剛度矩陣ke=PlatRectStif(Em,Th,Pr,a,b);for

el=1:neN(3:3:12)=3*ndel(el,:);N(2:3:12)=N(3:3:12)-1;N(1:3:12)=N(3:3:12)-2;K(N,N)=K(N,N)+ke;end%單元剛度矩陣U=zeros(3*nd,1); %

結點位移列向量U(dofree)=K(dofree,dofree)\F(dofree); %

求解fprintf("%12.3g\n",U(1)*

Em*Th^3/12/(1-Pr^2)/(2*Lx)^2));單元剛度矩陣function

ke=PlatRectStif(Em,Th,Pr,a,b)

%矩形板單元剛度矩陣syms

x

y

realx1=[1,-1,-1,1];y1=[1,1,-1,-1];ke=zeros(12,12);%定義符號變量%結點局部坐標%單元剛度矩陣D=Em*Th^3/12/(1-Pr^2)*[1,Pr,0;Pr,1,0;0,0,(1-Pr)/2];%彈性矩陣for

i=1:4for

j=1:4xi=x1(i);x0=x*xi;yi=y1(i);y0=y*yi;Ni=(1+x0)*(1+y0)*(2+x0+y0-x*x-y*y)/8;Nix=-b*yi*(1+x0)*(1+y0)*(1-y*y)/8;Niy=

a*xi*(1+x0)*(1+y0)*(1-x*x)/8;N=[Ni,Nix,Niy];Nixx=diff(diff(N,x),x)/a^2;Nixy=diff(diff(N,x),y)/a/b;Niyy=diff(diff(N,y),y)/b^2;Bi=-[Nixx;Niyy;2*Nixy];

%結合矩陣24單元剛度矩陣xj=x1(j);x0=x*xj;yj=y1(j);y0=y*yj;Nj=(1+x0)*(1+y0)*(2+x0+y0-x*x-y*y)/8;Njx=-b*yj*(1+x0)*(1+y0)*(1-y*y)/8;Njy=

a*xj*(1+x0)*(1+y0)*(1-x*x)/8;N=[Nj

Njx

Njy];

%Njxx=diff(diff(N,x),x)/a^2;Njxy=diff(diff(N,x),y)/a/b;Njyy=diff(diff(N,y),y)/b^2;Bj=-[Njxx;Njyy;2*Njxy];

%ke(3*i-3+(1:3),3*j-3+(1:3))=a*b*double(int(int(Bi"*D*Bj,x,-1,1),y,-1,1));endend25§8.2三角形薄板單元§8.2.1位移模式與形函數231xyo當對單元的兩個邊界分別平行于x和y軸時，確定系數的矩陣奇異。為此，引入面積坐標p(x,y)26§8.2三角形薄板單元§8.2.1位移模式與形函數231p(x,y)xyo一次二次三次在3個結點為零27位移模式28位移模式29形函數30導數關系31導數關系32導數關系33§8.2.1位移模式與形函數34§8.2.1位移模式與形函數35§8.2.1位移模式與形函數指標輪換得到結點2和3的形函數36位移的連續性31p(x,y)2xyo3次方程的4個系數由4個參數撓度和單元邊界方向的導數連續完全確定沿法向的導數是2次方程，有3個系數，不能由2個參數完全確定單元邊界法線方向的導數不連續三結點三角形單元不是協調元37§8.2.3應變幾何矩陣38§8.2.3應變幾何矩陣39§8.2.4單元剛度矩陣40§8.2.5等效結點力平板在分布集度為q的橫向力的作用下，等效結點力為，特別地，當q=q0為常數時，上式積分之一為41§8.2.7邊界條件典型邊界條件簡支邊固定邊自由邊yxo分布力偶42§8.2.6程序設計function

PlatTria%三角形薄板單元Em=210e9;mu=0.3;Th=2e-3;

%彈性模量，泊松比和板厚Lx=4;Nx=4;Ly=4;Ny=4;%x軸方向長度，單元數%y軸方向長度，單元數%單元總數，結點總數ne=2*Nx*Ny;

nd=(Nx+1)*(Ny+1);gxy=zeros(nd,2);for

i=0:Nxfor

j=0:Nygxy((Ny+1)*i+j+1,:)=[i*Lx/Nx,j*Ly/Ny];endend%置結點坐標xy15910①②③④⑧⑨32236

11

16

21715

20

251917242322128

13

1814⑤4

⑥xy2Lx2Lyo4344§8.2.6程序設計K=sparse(3*nd,3*nd);D0=Em*Th^3/12/(1-mu^2);%板剛度D=D0*[1,mu,0;mu,1,0;0,0,(1-mu)/2];for

el=1:neN=kron(3*ndel(el,:),[1,1,1])-kron([1,1,1],[2,1,0]);ke=PlatTraiStif(gxy(ndel(el,:),:),D);K(N,N)=K(N,N)+ke;EndU=zeros(3*nd,1);U(dofree)=K(dofree,dofree)\F(dofree);w0

=

0.0116;w1

=

U(1)/(4*Lx^2/(Em*Th^3/12/(1-mu^2)));fprintf("中心撓度誤差：(%6.3f/%6.3f-1)*100=%5.2f%%\n",w1,w0,(1-w1/w0)*1e2);§8.2.6程序設計function

ke=PlatTraiStif(xy,D)b=zeros(3,1);

c=zeros(3,1);syms

L1

L2

realA=det([xy,[1;1;1]])/2;for

i1=1:3i2=mod(i1,3)+1;i3=mod(i2,3)+1;b(i1)=xy(i2,2)-xy(i3,2);c(i1)=xy(i3,1)-xy(i2,1);endL3=1-L1-L2;

L=[L1,L2,L3];%三角形薄板單元剛度矩陣T=[b(1)*b(1),

b(2)*b(2),

2*b(1)*b(2);c(1)*c(1),

c(2)*c(2),

2*c(1)*c(2);2*b(1)*c(1),

2*b(2)*c(2),

2*b(1)*c(2)+

2*c(1)*b(2)]/4/A^2;45§8.2.6程序設計for

i1=1:3i2=mod(i1,3)+1;i3=mod(i2,3)+1;Sp=[L(i1)*(1+L(i1)*L(i2)+L(i1)*L(i3)-L(i2)^2-L(i3)^2);L(i1)^2*(b(i2)*L(i3)-b(i3)*L(i2))+(b(i2)-b(i3))*L(i1)*L(i2)*L(i3)/2;L(i1)^2*(c(i2)*L(i3)-c(i3)*L(i2))+(c(i2)-c(i3))*L(i1)*L(i2)*L(i3)/2]";dNdL1=diff(Sp,L1);dNdL2=diff(Sp,L2);B(:,3*i1-2:3*i1)=T*[diff(dNdL1,L1);diff(dNdL2,L2);2*diff(dNdL1,L2)];endke=2*A*double(int(int(B"*D*B,L2,0,1-L1),L1,0,1));end46§8.3考慮橫向剪切影響的平板彎曲單元§8.3.1結點力和結點位移1234567847Hencky板理論48§8.3.2位移模式1234567849§8.3.3應變與幾何矩陣(1)應變50§8.3.3應變與幾何矩陣(1)應變(2)幾何矩陣51§8.3.4應力、內力(1)應力52§8.3.4應力、內力(2)內力53§8.3.5剛度矩陣54§8.3.7程序設計y12923

28x④

22⑧

3651

56

6513

21③

2016

64435678145737

4240

474641

494850525311

17

25

31

39

45①16

⑤30

4410

15

24

29

38

435412

19

26

33②

18

325527

3534636261605958xyLLo5556§8.3.7程序設計function

[gxy,ndel,dofix,F,nd,ne,Th,Em,Pr,a,b,L]=PlatN8ExamEm=210e9;Pr=0.3;Th=1.6;a=1;b=1;L=16;ndel=[17,

3,

1,15,11,

2,10,16;19,

5,

3,17,12,

4,11,18;…

…65,51,49,63,56,50,55,64];gxy(:,1)=[zeros(9,1);ones(5,1);2*ones(9,1);3*ones(5,1);4*ones(9,1);5*ones(5,1;7*ones(5,1);8*ones(9,1)];gxy(:,2)=[0:8,0:2:8,0:8,0:2:8,0:8,0:2:8,0:8,0:2:8,0:8]";ne=size(ndel,1);nd=max(max(ndel));left=3*(1:9);down=3*[1:14:nd,10:14:nd];top=3*[9:14:nd,14:14:nd];right=3*(57:nd);dofix=[left,down-1,top-2,top-1,top,right-2,right-1,right];F=zeros(3*nd,1);F(1)=1/4;§8.3.7程序設計function

PlatRectN8

%8結點四邊形板單元[ndel,dofix,F,nd,ne,th,Em,pr,a,b,L]=PlatN8Exam;K=zeros(3*nd,3*nd);U=zeros(3*nd,1);dofree=setdiff(1:3*nd,dofix);

%無約束自由度D1=Em*th^3/12/(1-pr*pr)*[1,pr,0;pr,1,0;0,0,0.5-pr/2];%%單元剛度矩D2=1.0*Em*th/2/(1+pr)*eye(2);ke=PlatRectN8Stif(D1,D2,a,b);陣for

el=1:neN=kron(3*ndel(el,:),ones(1,3))-repmat([2,1,0],[1,8]);%組裝結構剛K(N,N)=K(N,N)+ke;度矩陣endU(dofree)=K(dofree,dofree)\F(dofree);5758輸出位移disp("Node

w

Rx

Ry")

%輸出標題for

j=1:ndfprintf("%4i%11.2e%11.2e%11.2e%11.2e%11.2e\n",j,U(3*j-2:3*j)));w0=5.6e-3*16^2/(Em*th^3/12/(1-pr*pr));w1=U(1);fprintf("Error=(%10.3g/%10.3g-1)*100=%5.2f%%\n’,

w1,w0,(1-w1/w0)*1e2));59計算內力disp("Node

Mx

My

Mxy

QyMoment=zeros(3,ne);Shear=zeros(2,ne);Qx")%彎矩和剪力for

el=1:ne

%計算單元應力Ue=U(kron(3*ndel(el,:),ones(1,3))-repmat([2,1,0],[1,8]));%幾何矩陣%彎矩%剪力[B1,B2]=PlatN8Strain([0,0],a,b);Moment(:,el)=D1*B1*Ue;Shear

(:,el)=D2*B2*Ue;endfor

el=1:nefprintf("%3i%11.2e%11.2e%11.2e%11.2e%11.2e%11.2e%11.2e\n",el,Moment(:,el),Shear(:,el));end幾何矩陣function

[B1,B2,detJ]=PlatN8Strain(xy,a,b)B1=zeros(3,24);B2=zeros(2,24);%形函數及其導數%雅可比矩陣及其行列式值%形函數導數[Sp,dSp]=PlanN8ShapeFun(xy,a,b);J=dSp*xy;detJ=det(J);dx=J\dSp;for

i=1:8B1(:,3*i-2:3*i)=[0,0,dx(i,1);0,-dx(i,2),0;0,-dx(i,1),dx(i,2)];B2(:,3*i-2:3*i)=[dx(i,2),-Sp(i),0;dx(i,1),0,Sp(i)];end60剛度矩陣%高斯積分點和積分權重%高斯積分function

ke=PlatRectN8Stif(D1,D2,a,b)p=sqrt(0.6)*[-1,0,1];w=[5,8,5]/9;ke=zeros(24,24);for

r=1:3for

s=1:3[B1,B2,J]=PlatN8Strain(p([r,s]),a,b);

%幾何矩陣ke=ke+a*b*w(r)*w(s)*(B1"*D1*B1+B2"*D2*B2)*J;%單元剛度矩陣endend61形函數function

[Sp,dSp]=PlanN8ShapeFun(xy,a,b)r=[1,-1,-1,1,

0,-1,0,1];

s=[1,1,-1,-1,

1,0,-1,0];Sp=zeros(8,1);

x=xy(1);

y=xy(2);for

i=1:8xi=r(i);

yi=s(i);

x0=xi*x;

y0=yi*y;Sp(i)=(1+x0)*(1+y0)*(x0+y0-1)*xi^2*yi^2/4...+(1-x^2)*(1+y0)*(1-xi^2)*yi^2/2+(1-y^2)*(1+x0)*(1-yi^2)*xi^2/2;endif

nargout>1dSp=zeros(2,8);for

i=1:8xi=r(i);

yi=s(i);

x0=xi*x;

y0=yi*y;dSp(:,i)=[(2*x0+y0)*(1+y0)*xi^3*yi^2/4-x*(1+y0)*(1-xi^2)*yi^2+(1-y^2)*(1-yi^2)(2*y0+x0)*(1+x0)*yi^3*xi^2/4-y*(1+x0)*(1-yi^2)*xi^2+(1-x^2)*(1-xi^2)*yi^3/endend6263輸出計算結果Node

MxMy

Mxy

Qy

Qx1.37e-01

-2.59e-02

-5.57e-02

-5.57e-023.34e-02

-1.66e-02

-