EulerBernoullibeamtheroy(歐拉伯努利梁理论与有限元计算).docVIP

Euler-Bernoulli beam 一、理论部分 Euler-Bernoulli beam 假设 由(1.1)式可得 虚位移原理 其中 令，为单元长度，则上式成为 单元节点位移取为 令 其中形函数 分别对式(1.7)、(1.8)求一阶和两阶导数得 由，可得 将(1.6)、(1.11)式代入(1.4)式，可得 刚度矩阵 二、算例 本文所举算例为悬臂梁端部受竖向载荷，挠度理论曲线为 具体数据如下： ，，，，， 划分10个单元，11个结点。 由上图可以看出，数值模拟结构跟理论解吻合非常好。 源代码(Matlab) % program_main initial_statistic; le=pretreatment(ND,NE,lamda,x,y); [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]; function 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); end end function [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,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; End function 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),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*l2)+Ke(4:6,4:6); end end function 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); end end function [x,Hw]=Ngauss(N) % 给定高斯点个数返回高斯点出坐标和权系数值； s