EulerBernoulli beam theroy欧拉伯努利梁理论与有限元计算

1、Euler-Bernoulli beam一、理论部分Euler-Bernoulli beam 假设由(1.1)式可得虚位移原理其中令，为单元长度，则上式成为单元节点位移取为令其中形函数分别对式(1.7)、(1.8)求一阶和两阶导数得由，可得将(1.6)、(1.11)式代入(1.4)式，可得刚度矩阵二、算例本文所举算例为悬臂梁端部受竖向载荷，挠度理论曲线为具体数据如下：，划分10个单元，11个结点。由上图可以看出，数值模拟结构跟理论解吻合非常好。源代码(Matlab)% program_maininitial_statistic;le=pretreatment(ND,NE,lamda,x,y);

2、Dis,K=Euler_Bernoulli_beam(ND,NE,lamda,EA,EI,le,P,Node_c,Dis_c);post(Dis,x,y,ND,NE);% initial statistic;ND=11;NE=10;EA=2.1e8*ones(NE,1);EI=4.375e4*ones(NE,1);x=(0:0.1:1)'y=zeros(11,1);lamda=1,2;2,3;3,4;4,5;5,6;6,7;7,8;8,9;9,10;10,11;P=zeros(3*ND,1);P(32)=-1.0e4;Node_c=1,2,3'Dis_c=0,0,0'f

3、unction le=pretreatment(ND,NE,lamda,x,y)% the program of pretreatment;for i=1:NE le(i)=sqrt(x(lamda(i,1)-x(lamda(i,2)2+(y(lamda(i,1)-y(lamda(i,2)2);endendfunction Dis,K=Euler_Bernoulli_beam(ND,NE,lamda,EA,EI,le,P,Node_c,Dis_c)% solve the Euler Bernoulli beam problem by the finite element method;% ND

4、,NE,lamda the information of Node,Element;% EA,EI,le the properties of material;K=Stiffness(ND,NE,lamda,EA,EI,le);K,P=Import_c(K,P,Node_c,Dis_c);Dis=inv(K)*P;Endfunction K=Stiffness(ND,NE,lamda,EA,EI,le)% 集成总体刚度矩阵；% ND节点数，NE单元数，lamda单元节点定位向量；K=zeros(3*ND);Ke=zeros(6);for ie=1:NE Ke=Estiffness(EA(ie)

5、,EI(ie),le(ie),2); l1=lamda(ie,1);l2=lamda(ie,2); K(3*l1-2:3*l1,3*l1-2:3*l1)=K(3*l1-2:3*l1,3*l1-2:3*l1)+Ke(1:3,1:3); K(3*l1-2:3*l1,3*l2-2:3*l2)=K(3*l1-2:3*l1,3*l2-2:3*l2)+Ke(1:3,4:6); K(3*l2-2:3*l2,3*l1-2:3*l1)=K(3*l2-2:3*l2,3*l1-2:3*l1)+Ke(4:6,1:3); K(3*l2-2:3*l2,3*l2-2:3*l2)=K(3*l2-2:3*l2,3*l2-2:3

6、*l2)+Ke(4:6,4:6);endendfunction Ke=Estiffness(EA,EI,le,N)% 求单元刚度矩阵；% N为高斯点个数；x,Hw=Ngauss(N);Ke=zeros(6);for i=1:N Bu=-0.5,0,0,0.5,0,0' Bv=0,1.5*x(i),(0.75*x(i)-0.25)*le,0,-1.5*x(i),(0.75*x(i)+0.25)*le' Ke=Ke+(Bu*Bu'*2/le*EA+Bv*Bv'*(2/le)3*EI)*Hw(i);endendfunction x,Hw=Ngauss(N)% 给定高斯

7、点个数返回高斯点出坐标和权系数值；switch N case 1 x=0;Hw=2; case 2 x=-0.577350269189626,0.577350269189626'Hw=1,1' case 3 x=-0.774596669241483,0,0.774596669241483' Hw=0.555555555555556,0.888888888888889,0.555555555555556' case 4endendfunction K,P=Import_c(K,P,Node_c,Dis_c)% Import boundary condition;a

8、=1.0e10;for i=1:length(Node_c) %K(Node_c(i),Node_c(i)=a*K(Node_c(i),Node_c(i); %P(Node_c(i)=K(Node_c(i),Node_c(i)*Dis_c(i); K(Node_c(i),:)=0;K(:,Node_c(i)=0;K(Node_c(i),Node_c(i)=1; P(Node_c(i)=0;endendfunction post(Dis,x,y,ND,NE)% the program of post treatment;% x,y the coordinate of Nodes;x=x+Dis(1:3:(3*ND);y=y+Dis(2:3:(3*ND);figure,plot(x,y,'-o');% 本程序以悬臂梁端部承受集中力为例，解析解为% y=-0.5*Pl/EI*x2+1/6*P/EI*x3;P=1.0e4;l=1;EI=4.375e4;x_s=0:0.1:1;y_s=-0.5*P*l/EI*x_s.2+1/6*P/EI*x_s.3;hold on,plot(x_s,y_s,'r')xlabel('x(m)')ylabel('the defl