有限元方法与MATLAB程序设计 第9章 壳的弯曲

第9章壳的弯曲§9.1平面三角形壳体单元§9.2考虑横向剪切变形影响的八结点壳体单元§9.1局部坐标系下单元刚度矩阵单元坐标系下的单元刚度矩阵平面问题刚度子矩阵平板弯曲刚度子矩阵单元坐标系下的单元刚度方程2§9.2结构坐标系下的单元刚度矩阵坐标转换矩阵3坐标转换矩阵的说明4结点位移和结点力单元坐标系结构坐标系5单元刚度矩阵结构坐标系下的单元刚度矩阵结构坐标系下的单元刚度方程某个结点周围的单元在同一个平面内，最后一个自由度的刚度就会成为零直接删除这个自由度。假设一个任意不为零的刚度值。6§9.1.4程序设计算例Nx=4;Ny=4;ne=2*Nx*Ny;nd=(Nx+1)*(Ny+1);xyz=zeros(nd,3);ndel=zeros(ne,3);for

i=0:Nxfor

j=0:Nya=i*pi/Nx/2;%结点坐标和单元信息xyz((Ny+1)*i+j+1,:)=25*[sin(a),j/Ny,cos(a)];endendEm=210e9;mu=0.0;Th=0.25;

%弹性模量，泊松比1和板厚6y523478910x①⑨3230⑧⑥⑥⑤④③②zF25x257o§9.1.4程序设计算例for

i=1:Nxfor

j=1:Nyn1=(Ny+1)*

i+j+1;n2=(Ny+1)*(i-1)+j+1;el=(Ny*(i-1)+j)*2;ndel(el,

:)=[n1,

n2,

n2-1];endenddofix=1:6*Ny+6;ndel(el-1,:)=[n2-1,n1-1,n1];dofree=setdiff(1:6*nd,dofix);F=zeros(6*nd,1);F(6*Nx*(Ny+1)+1:6:6*nd)=1e5*[0.5;ones(Ny-1,1);0.5]/zF25x25oy51234678910x⑨3230⑧⑥⑥⑤④③②①8单元刚度矩阵D=Em/(1-mu^2)*[1,mu,0;mu,1,0;0,0,(1-mu)/2];Np=[1,2,7,8,13,14];Nb=[3,4,5,9,10,11,15,16,17];M=-5:0;K=zeros(6*nd,6*nd);ke=zeros(18,18);for

el=1:neg=xyz(ndel(el,:),:)"-repmat(xyz(1,:),[3,1])";l1=sqrt(g(:,2)"*g(:,2));e1=g(:,2)/l1;e3=cross(e1,g(:,3));e3=e3/sqrt(e3"*e3);e2=cross(e3,e1);

e2=e2/sqrt(e2"*e2);t=[e1,e2,e3];T=[t,zeros(3,3);zeros(3,3),t];s0=g(:,3)"*e1;d0=g(:,3)"*e2;xy=[0,0;l1,0;s0,d0];[B,A]=PlanTriaStrain(xy);ke(Np,Np)=Th*B"*D*B*A;ke(Nb,Nb)=Th^3/12*PlatTraiStif(xy,D);910单元刚度矩阵for

i=6:6:18ke(i,i)=ke(i,i)+1e0;endfor

i=1:3;

I=6*ndel(el,i)+M;for

j=1:3;J=6*ndel(el,j)+M;K(I,J)=K(I,J)+T*ke(6*i+M,6*j+M)*T";endendend输出位移U=zeros(6*nd,1);U(dofree)=K(dofree,dofree)\F(dofree);%SolveU1=pi*1e5*25^3/(4*Em*25*Th^3/12);

%解析解disp(smuintf("%14.7g,%14.7g,",U(6*nd-5),U1);§9.2考虑横向剪切变形影响的八结点壳体单元15263748o311114127856结点位移列向量11结点力列向量中面结点坐标结点中面法线任意点坐标§9.2.1单元坐标系1125263748假设结构变形导致法线转动向量定义2个与V3i垂直的向量§9.2.2位移列阵与形函数结点i处法线上任意点的位移可以用结点i的位移和相对结点i位移叠加位移模式13§9.2.3应变与几何矩阵应变几何子矩阵14应变坐标转换单元坐标系下的应变与结构坐标系下应变的关系应变坐标转换矩阵局部坐标轴方向向量15§9.2.4应力与弹性矩阵刚度矩阵结构坐标系下的弹性矩阵16局部坐标系下的弹性矩阵§9.2.5程序设计function

[gxy,ndel,dofix,F,nd,ne,Th,Em,mu]

=

Eaxm9_2Em

=

4.32e9;

mu

=

0.3;Th

=

0.25;

R

=

25;nx

=

6;

ny

=

3;

n

=

2*ny+1;nd

=

(2*nx+1)*(2*ny+1)-nx*ny;

ne

=

nx*ny;a

=

linspace(0,pi/2,2*nx+1)";x

=

sin(a)*(R+[Th,-Th]/2);y

=

linspace(0,4,2*ny+1);z

=

cos(a)*(R+[Th,-Th]/2);gxy

=

zeros(nd,3,2);for

i

=

1:2gxy(:,1,i)

=

repelem(x(:,i),[repmat([n,ny+1],1,nx),n]);gxy(:,2,i)

=

[repmat([y,y(1:2:end)],1,nx),y];gxy(:,3,i)

=

repelem(z(:,i),[repmat([n,ny+1],1,nx),n]);endzF25x25o17§9.2.5程序设计for

i=1:nx%计算单元信息for

j=1:nyn1

=

(3*ny+2)*(i-1)+2*j-1;

n2

=

(3*ny+2)*(i-1)+j+2*ny+1;n3

=

n2+ny+j+1;ndel(ny*(i-1)+j,:)

=

[n3+1,n1+2,n1,n3-1,n2+1,n1+1,n2,n3];endenddofix

=

1:5*n;F

=

zeros(5*nd,1);F(5*(nd-n+1:nd)-4)

=

2*[1,repmat([4,2],1,ny-1),4,1]/ny/3;zF25x25o1819§9.2.5程序设计function

ShelN8[xyz,ndel,dofix,F,nd,ne,Th,Em,mu]=Exam9_2;%8结点四边形壳单元。K=zeros(5*nd,5*nd);U=zeros(5*nd,1);dofree

=

setdiff(1:5*nd,dofix);D

=

ShelN8Elastic(Em,mu);for

el

=

1:neN=repelem(5*ndel(el,:),5)+repmat(-4:0,1,8);

%单元自由度%结构刚度矩阵K(N,N)

=

K(N,N)+ShelN8Stif(xyz(ndel(el,:),:,:),D,Th);endU(dofree)

=

K(dofree,dofree)\F(dofree);20§9.2.5程序设计disp("Node

X

Yfor

j

=

1:ndZ

u

v

w

alpha

beta")fmuintf("%4i%6.2f%6.2f%6.2f%10.2e%10.2e%10.2e%10.2e%10.2e\n",j,(xyz(j,:,1)+xyz(j,:,2))/2,U(5*j-4:5*enddisp("Elem

1

2

3

4

5

6

7

8Angle")Sx

Sy

Sz

Sxy

Syz

S1

S2for

el=1:ne

%计算单元应力T

=

ShelN8Rotation([0,0],xyz(ndel(el,:),:,:),Th);N

=

repelem(5*ndel(el,:),5)+repmat(-4:0,1,8);for

i

=

0:1B

=

ShelN8Strain([0,0,i],xyz(ndel(el,:),:,:),Th);%几何矩阵S=D*T*B*U(N);

%应力[Dir,S1]=eig(S([1,3;3,2]));%求特征值，计算主应力方向及大小fmuintf([repmat("%3d",1,9),repmat("%11.3e",1,8),"\n"],el,ndel(el,:),S,diag(S1),atan2d(Dir(2,1),Dir(endend21§9.2.5程序设计function

ke

=

ShelN8Stif(xyz,D,th)p

=

[-1,1]/sqrt(3);

w

=

[1,1];ke

=

zeros(40,40);for

i

=

1:2for

j

=

1:2T

=

ShelN8Rotation(p([i,j]),xyz,th);for

k

=

1:2[B,J]

=

ShelN8Strain(p([i,j,k]),xyz,th);B1

=

T*B;ke

=

ke+w(i)*w(j)*w(k)*B1"*D*B1*J;endendend单元刚度矩阵§9.2.5程序设计坐标转换矩阵function

T

=

ShelN8Rotation(p,xyz,th)V3

=

zeros(3,1);Sp

=

PlanN8ShapeFun(p);for

i

=

1:8V3

=

V3+2*(xyz(i,:,1)-xyz(i,:,2))"*Sp(i)/th;endV3

=

V3/sqrt(V3"*V3);V1

=

cross([1;0;0],V3);V1

=

V1/sqrt(V1"*V1);V2

=

cross(V3,V1);T

=[V1(1)^2,V1(2)^2,V1(3)^2,V1(1)*V1(2),V1(2)*V1(3),V1(3)*V1(1);V2(1)^2,V2(2)^2,V2(3)^2,V2(1)*V2(2),V2(2)*V2(3),V2(3)*V2(1);2*V1(1)*V2(1),2*V1(2)*V2(2),2*V1(3)*V2(3),V1(1)*V2(2)+V2(1)*V1(2),

...V1(2)*V2(3)+V2(2)*V1(3),V1(3)*V2(1)+V2(3)*V1(1);2*V2(1)*V3(1),2*V2(2)*V3(2),2*V2(3)*V3(3),V2(1)*V3(2)+V3(1)*V2(2),

...V2(2)*V3(3)+V3(2)*V2(3),V2(3)*V3(1)+V3(3)*V2(1);2*V3(1)*V1(1),2*V3(2)*V1(2),2*V3(3)*V1(3),V3(1)*V1(2)+V1(1)*V3(2),

...V3(2)*V1(3)+V1(2)*V3(3),V3(3)*V1(1)+V1(3)*V3(1)];2223function

[Sp,dSp]

=

PlanN8ShapeFun(xy)r

=

[-1,1,1,-1,0,1,0,-1];

x

=

xy(1);s

=

[-1,-1,1,1,-1,0,1,0];

y

=

xy(2);

Sp

=

zeros(8,1);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)*xi^3/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/2];endend§9.2.5程序设计形函数§9.2.5程序设计function

[B,detJ]=ShelN8Strain(p,xyz,th坐)标转换矩阵B

=

zeros(6,40);[Sp,dSp0]

=

PlanN8ShapeFun(p(1:2));J

=

ShelN8Jac