西交大有限元原理及应用-大作业

1、-作者xxxx-日期xxxx西交大有限元原理及应用-大作业【精品文档】有限元原理及工程应用 大作业学 院： 机械工程学院班 级： 硕4002班小组成员： 李 追 3114001089 陈 草 3114001080【精品文档】作业题目：利用有限元方法对简支梁问题进行求解，梁的横截面为矩形，其约束情况如图1所示。已知梁的几何尺寸和物理参数如下：（1）几何尺寸：长度，截面尺寸；（2）物理参数：弹性模量GPa，泊松比，密度。图1梁及其横截面示意图要求：(1) 至少划分五个节点（四个单元）；(2) 给出单元节点信息；(3) 给出单元刚度矩阵和质量矩阵；(4) 给出总刚度矩阵和总质量矩阵；(5) 求出梁各

2、界固有频率及振型（五阶）；(6) 将所得结果与理论值进行对比，验证方法的可行性。解：由有限元知识，根据Rayleigh-Ritz法，解有限元分为四步：建立离散化、单元分析、形成总体方程、解方程，具体步骤如下：（1）建立离散化这里我们将矩形截面简支梁等分四等分，即分为六节点的五个杆单元，如图2所示：每个单元尺寸，这里只考虑杆在竖直平面的弯曲，每个节点只有y方向位移和绕z轴的旋转自由度。（2）单元分析构造一组Lagrange插值基函数，在本节点值为1，其他节点值为0。从Rayleigh-Ritz法可以看到，插值函数要p次可微，最高阶导数出现在应变能表达式中；同样，我们可以这一原则适用于基函数的选择

3、以及形状函数，否则我们将无法正确计算应变能当我们使用有限元逼近方法。梁的弯曲问题，应变能计算公式： （1-1）其中，E为弹性模量，Iz为截面惯性矩。从公式可知，位移函数必须连续，并且二阶导数平方可积。如图3，是一维杆单元模型，每个节点两个自由度，该单元含有四个自由度，即（）。本题中我们采用三次多项式插值函数： （1-2）因此，我们必须给出四个形函数（位移模式）。图3 一维杆单元模型1) 构造Hermite插值函数。选择局部坐标系(,)，其中l是单元长度，转角是挠度值的一阶导数，定义边界条件： （1-3）因此，我们给出变形的Hermite的多项式插值函数： （1-4）其中，和分别满足如下条件，对

4、应的图形如图4所示： （1-5）图4 一维Hermite插值多项式2) 基于Langrage和Hermite插值多项式，写出单元形函数 （1-6）节点位移值也可以得出 （1-7）同时，表达式（1-7）用矩阵表示为 （1-8）其中，3) 用能量表达式替代表达式中的和。动能表达式： （1-9）将（1-8）带入（1-9），得到 （1-10）从而获得质量矩阵： （1-11）带入， （1-12）应变能表达式： （1-13）刚度矩阵表达式： （1-14）带入，可以得到 （1-15）（3） 形成总体方程将每个杆单元的能量方程组装。完整梁上的的总动能和能量的和所做的总功梁上的外力作用，所有的自由度的位移矢量可

5、以给出： （1-16）将位移矢量转换为全局坐标系下的位移矢量，变换矩阵为： （1-17） （1-18） 分别以矩阵形式给出动能和质量矩阵： （1-19） （1-20）因此，总体质量矩阵为总体刚度矩阵：求固有频率。总应变能 (1-21) (1-22)利用Lagrange方程，推导简支梁自由振动方程： (1-23)这里，我们假设简支梁做简谐振动，则 (1-24)因此，特征方程为： (1-25)其中，为固有频率。（4） 解方程，有限元分析结果基于上述理论,我们获得了采用MATLAB的有限元分析程序代码。提交边界条件、材料特性和几何参数到上面的方程, 使用MATLAB代码我们得到以下结果。（1）离散简

6、支梁为5单元6节点,那么我们得到的单元质量矩阵和刚度矩阵,如下所示： （2）总质量矩阵和总刚度矩阵如下：(3) 简支梁的振动分析 表1列出了简支梁的五阶固有振动频率。从表中我们可以看出，有限元模拟分析方法和理论值在误差允许范围内是比较吻合的。计算得离散为五单元下的简支梁固有振动频率：计算值（）理论值()误差率（%）一阶固有振动频率0.0107 二阶固有振动频率三阶固有振动频率四阶固有振动频率0.2901 五阶固有振动频率离散为30单元的简支梁的振动模式如下图所示。 一阶振动图像 二阶振动图像 三阶振动图像 四阶振动图像 五阶振动图像 六阶振动图像3. 采用Matlab编写的程序代码% 利用有限

7、元方法求解简支梁的振动问题 %clc ; clear all;syms x l rho b t EA = b*t;I = b*t3/12;n = input ('Please input the number of discrete elements n = '); % 输入离散化单元的数量n% 定义形函数N1 = 1-3*(x/l)2+2*(x/l)3; N2 = (x/l-2*(x/l)2+(x/l)3)*l;N3 = 3*(x/l)2-2*(x/l)3;N4 = (x/l)3-(x/l)2)*l;% %求解单元质量矩阵、刚度矩阵以及总质量矩阵和刚度矩阵N = N1,N2,

8、N3,N4;Me0 = int(N'*N,x,0,l); Me = rho*A*Me0Ke0 = int(diff(N.',x,2)*diff(N,x,2),x,0,l) Ke = (E*I)*Ke0;M0 = zeros(2*(n+1),2*(n+1);K0 = zeros(2*(n+1),2*(n+1);for i=1:1:n ae = zeros(4,2*(n+1); for j=1:1:4 ae(j,2*i+j-2) = 1; % 定义坐标变换矩阵 end M0 = M0+ae.'* Me * ae; K0 = K0+ae.'* Ke *ae; endd

9、isp('The element mass matrix Me = '); disp(Me); disp('The element stiffness matrix Ke = ');disp(Ke); disp('The total mass matrix M = ');disp(M0); disp('The total stiffness matrix K = ');disp(K0); %Me1 = matlabFunction(Me);Ke1 = matlabFunction(Ke);M1 = matlabFunction(M

10、0); % 将质量符号矩阵转换为代数矩阵K1 = matlabFunction(K0); % 将刚度符号矩阵转换为代数矩阵% 输入的几何参数和物理常数L = 0.4; % 梁的长度, mb = 0.02; % 梁的宽度, mt = 0.002; %梁的厚度, m E = 0.7*1011; % 梁的弹性模量, GParho = 2700; %梁的密度, kg/m3l = L/n;MeNumeric = Me1(b,l,t,rho);KeNumeric = Ke1(E,b,l,t);MNumeric = M1(b,l,t,rho);KNumeric = K1(E,b,l,t);disp('

11、;The numerical element mass matrix Me= ');disp(MeNumeric); % 输出单元质量矩阵disp('The numerical element stiffness matrix Ke = ');disp(KeNumeric); % 输出单元刚度矩阵disp('The numerical total mass matrix M = ');disp(MNumeric); %输出总质量矩阵disp('The numerical total stiffness matrix K = ');disp

12、(KNumeric); %输出总刚度矩阵% 简支梁的振动求解%对于简支梁考虑约束条件： ;%采用消元法：消去整体质量矩阵的第1行，第1列；倒数第二行和倒数第二列% 消去整体刚度矩阵的第1行，第1列；倒数第二行和倒数第二列Mssb = MNumeric(2:2*n,2*(n+1),2:2*n,2*(n+1);Kssb = KNumeric(2:2*n,2*(n+1),2:2*n,2*(n+1);Xssb,lamdassb=eig(Kssb,Mssb);omegassb=sort(sqrt(diag(lamdassb);XXssb=zeros(n,n);XXssb(1,:)=0;XXssb(n,:

13、)=0;for i=2:n-1 for j=1:n XXssb(i,j)=Xssb(2*(i-1),j); endend% 根据振动力学方程得到简支梁固有频率betalssb = ;omega_realssb = ;errorssb = ;for i=1:10 betalssb(i) = pi*i;omega_realssb(i) = (betalssb(i)/L)2*sqrt(E*(b*t3/12)/(rho*b*t);enddisp(' True output = ');disp(omega_realssb); % 输出理论固有频率disp(' Calculated

14、 Output = ');disp(omegassb); % 输出计算固有频率值for i = 1:10errorssb (i) = (abs(omega_realssb(i)-omegassb(i)/omega_realssb(i)*100; % 输出误差率enddisp(' Output percentage error = ');disp(errorssb); % % 画出前五阶振型r = 1:n;for i = 1:n figure('Color',1 1 1); plot(r,XXssb(:,i); xlabel('itxrm/itl','FontName','Times New Roman','FontSize',20); ylabel('urm(itxrm)','FontName','Times New Roman','FontSize',20); set(gca,'Linewidt