有限单元法实验
有限元实验报告
有限元实验报告一、实验目的本实验旨在通过有限元方法对一个复杂的工程问题进行数值模拟和分析，从而验证理论模型的正确性，优化设计方案，提高设计效率。
二、实验原理有限元方法是一种广泛应用于工程领域中的数值分析方法。
它通过将连续的求解域离散化为由有限个单元组成的集合，从而将复杂的偏微分方程转化为一系列线性方程组进行求解。
本实验将采用有限元方法对一个具体的工程问题进行数值模拟和分析。
三、实验步骤1、问题建模：首先对实际问题进行抽象和简化，建立合适的数学模型。
本实验将以一个简化的桥梁结构为例，分析其在承受载荷下的应力分布和变形情况。
2、划分网格：将连续的求解域离散化为由有限个单元组成的集合。
本实验将采用三维四面体单元对桥梁结构进行划分，以获得更精确的数值解。
3、施加载荷：根据实际工况，对模型施加相应的载荷，包括重力、风载、地震等。
本实验将模拟桥梁在车辆载荷作用下的应力分布和变形情况。
4、求解方程：利用有限元方法，将偏微分方程转化为线性方程组进行求解。
本实验将采用商业软件ANSYS进行有限元分析。
5、结果后处理：对求解结果进行可视化处理和分析。
本实验将采用ANSYS的图形界面展示应力分布和变形情况，并进行相应的数据处理和分析。
四、实验结果及分析1、应力分布：通过有限元分析，我们得到了桥梁在不同工况下的应力分布情况。
如图1所示，桥梁的最大应力出现在支撑部位，这与理论模型预测的结果相符。
同时，通过对比不同工况下的应力分布情况，我们可以发现，随着载荷的增加，最大应力值逐渐增大。
2、变形情况：有限元分析还给出了桥梁在不同工况下的变形情况。
如图2所示，桥梁的最大变形发生在桥面中央部位。
与理论模型相比，有限元分析的结果更为精确，因为在实际工程中，结构的应力分布和变形情况往往受到多种因素的影响，如材料属性、边界条件等。
通过对比不同工况下的变形情况，我们可以发现，随着载荷的增加，最大变形量逐渐增大。
3、结果分析：通过有限元分析，我们验证了理论模型的正确性，得到了更精确的应力分布和变形情况。
《弹性力学问题的有限单元法》
《弹性力学问题的有限单元法》弹性力学问题的有限单元法（FiniteElementMethod，简称FEM）是一种经典的多学科跨领域的计算方法，它用于估算连续体结构中非线性材料力学性能，如强度、刚度和破坏。
有限单元法已成为工程和材料科学中最重要的数值计算方法，可用于解决各种复杂多学科优化和设计问题。
有限单元法的基本思想是把复杂的连续体结构划分成许多小的、较容易处理的有限元素，而不是像一般的解析方法那样求取整体的解析解。
基于有限元素重要的性质，即小元素经过一系列的连接后就可以构成整个结构的模型，有限单元法的本质是数值分析，也就是根据模型的物理知识，选择有效的数值化方法，用数值计算的方法求解所要求的结果，从而使这些数值计算结果符合实际结构物理知识。
有限单元法是一种有效计算弹性力学问题的方法，它可以用来求解任意形状的结构问题，无论是有边界条件还是无边界条件，无论是线性或者非线性的形状变化，有限单元法都能够有效地应用。
其优势在于以节省计算时间和消耗的成本，在特殊的材料条件下，它可以比较快速地获得弹性力学问题的有效精确解。
其精度依赖于计算模型元素的类型、形状和几何尺寸等，因此通常需要调节元素的类型、形状和尺寸，以满足计算需要。
在计算机技术的发展下，有限单元法的计算能力越来越强大，可以对更多的复杂问题进行分析，可以更有效地解决工程设计中的实际问题。
由于计算机可以模拟各种变形和应力的变化，因此有限单元法可以为工程设计和材料研究提供更可靠的结果。
有限单元法在工程应用中的实际作用是显而易见的。
它不仅可以用来计算弹性结构中的材料力学特性，还可以分析复杂结构的动态响应。
此外，有限单元法还可以用来计算弹性结构中的表面张力、刚度，以及各种材料的裂缝扩展。
通过有限单元法的应用，可以获得有效的数值结果，从而提高设计效果和工程安全性。
因此，有限单元法对于材料科学和工程设计都具有重要价值，今后还将发挥更多的功能。
有限单元法是多学科跨学科的计算方法，它可以用来有效地分析复杂形状结构的力学特性，计算出精确的结果，从而提高工程设计的效果和安全性。
有限单元法第二版课程设计 (2)
有限单元法第二版课程设计一、设计背景有限单元法是一种常用的分析方法，广泛应用于工程学和自然科学领域。
为了进一步提高学生的有限单元法水平，本次课程设计旨在设计一个较为完整的有限单元方法分析项目。
二、设计目标通过本次课程设计，旨在让学生深入了解有限单元方法的原理和实现过程，提高学生的分析和解决实际问题的能力。
三、设计内容本次课程设计的主要内容包括以下三个部分：1.有限单元法的基础知识学习是本次课程设计的首要任务。
学生应该充分掌握有限单元法的基本原理、有限单元法的应用领域、有限单元法的基本步骤、有限单元法的精度等相关知识，为后续的分析工作奠定基础。
2.本次课程设计的重点是学生自行选择一个实际问题，并建立相应的有限单元模型，进行静态、动态或热力学分析。
学生应该根据具体情况选择不同的求解方法，如使用有限元软件求解或自编程求解。
3.在完成有限元分析后，学生应该对结果进行分析和讨论。
包括模型的合理性、分析结果的精度和可靠性等等，对分析结果进行进一步的解释和讨论。
四、设计要求1.本次课程设计应该由每个学生独立完成，不得相互抄袭和抄袭现成的模型。
2.学生自行选择并设计仿真模型，可以是自行查找的数据或者自己设计的模型。
3.分析结果应该以文本的方式进行输出，要求输出结果应该包括模型的详细说明、分析结果和分析结论等内容。
4.报告应该能够详细说明分析流程，从建模、求解到结果的呈现，必须清晰且易于理解。
5.学生应该按照教师要求的时间和形式，将完成的报告提交给教师评分。
五、总结有限单元法是一种重要的计算方法，对于提高学生的工程实践能力和实际应用技能有着重要的作用。
通过本次课程设计的学习，有助于学生深入理解有限单元法，将学校理论与实际问题相结合，为将来的工作打下坚实的基础。
有限单元法
ANSYS软件是融结构、流体、电场、磁场、声场分析于一体 软件是融结构、流体、电场、磁场、 软件是融结构 的大型通用有限元分析软件。由世界上最大的有限元分析软 的大型通用有限元分析软件。 件公司之一的美国ANSYS开发，它能与多数 开发， 软件接口， 件公司之一的美国 开发 它能与多数CAD软件接口， 软件接口 实现数据的共享和交换，是现代产品设计中的高级CAE工具 实现数据的共享和交换，是现代产品设计中的高级 工具 之一。 之一。
三、应用实例
运用有限元分析软件ANSYS对齿轮做了参数设计与网格划分、模态分析， 求出了齿轮的前5阶模态，得到相应的固有频率和振型图。利用振型图，直 观地分析了圆柱齿轮的动态特性并发现薄弱环节，为齿轮的动态性能实验、 设计和维护提供理论依据，同时也可为结构系统的动态响应计算和分析打下 基础。
谢谢！
一、有限单元法的分析过程 直接法的分析过程： 直接法的分析过程： 结构的离散化： ⑴结构的离散化：用有限个方位不同但几何性质及物理性质均相似
的单元组成的集合体代替原来的连续体和结构。
⑵单元分析： 单元分析： ①选取位移模式：{ f } = [N ]{δ }e ②利用几何方程： {ε } = [B ]{δ } ③利用物理方程：{σ } = [S ]{δ } ④利用平衡条件： {R} = [K ]{δ }
二、介绍一款有限元软件（ANSYS
模型变形前后的对比图
二、介绍一款有限元软件（ANSYS） 介绍一款有限元软件（ANSYS）
模型的应形图
模型的应力图
三、应用实例
齿轮传动因其诸多优点，成为机械传动中最主要的一类传动，齿轮副在工 作时，在内部和外部激励下将发生机械振动 。对于一个齿轮，它具有固定 的振动特性，一般包括固有频率和振型，它是系统的动态特性之一 。在这 里，我们对一个给定的齿轮进行模态分析，得到其固有频率和振型。 问题描述：齿轮的齿数Z=24，模数m=3.5，齿宽b=15mm。通过ANSYS软件 求其模态分析。
有限单元法
有限单元法人们常说：“教学有法，教无定法。
”的确，要提高语文课堂教学的质量，要提高学生的语文素养，教师不能一味地把知识灌输给学生，而应该为学生营造轻松、自主、开放的课堂氛围，从而提高学生的学习兴趣。
如何将课堂活动落到实处？语文老师们苦思冥想，找出了许多种教学方法，但这些教学方法都存在一个共同的问题：一节课下来，学生的知识似乎没有增加多少，他们好像只懂得了听讲，对知识点不求甚解，效果可见不佳。
那么怎样才能让学生在有限的时间内既扎实基础又培养能力呢？有限单元法可以助你一臂之力。
这就是有限单元法。
在上《夏天里的成长》这篇课文时，我把全班分成了三组，每一组负责查阅《大自然的语言》《夏天里的成长》和《童年的水墨画》三篇课文。
每个小组安排一名组员负责摘抄三篇课文中具有代表性的段落，并把它们进行分类整理，写出自己的感受。
这一环节引导学生在课外对课文进行深入地了解，发挥了课本学习的延伸作用。
《夏天里的成长》一课中，安排了三次关于“蝉鸣”的交流讨论，我告诉学生“不同的季节会听到不同的蝉声，我们所熟悉的蝉声就来自这个春天……”“请大家拿出各自的工具书，通过字典或百度来了解一下‘蝉’这个字的含义。
”“‘鸣’的古意是什么？”通过交流与探讨，同学们纷纷表示会收集“鸣”的资料，丰富自己的知识。
整个过程轻松愉快，活跃了课堂气氛，培养了学生读书的好习惯。
除了这些，我还用了有限单元法设计了“一石激起千层浪”这一环节，精心创设教学情境，使学生置身于具体的情境之中，受到熏陶，得到启迪。
在交流讨论时，有同学提出“有的蝉是好几年才叫一次的，一辈子就叫一回，也有的蝉在一年中的不同时候都叫……那么蝉为什么叫的次数不同呢？”面对这样的问题，我们没有急于给出答案，而是鼓励学生继续查阅资料，多思考，相信他们肯定会带着这个问题走进下一课。
这一环节的设计巧妙利用了网络资源，拓宽了学生的视野，开阔了学生的思路，学生仿佛一下子解开了心中的疑惑，收获良多。
有限单元法大作业
有限单元法大作业问题：一薄板，尺寸为10x10，厚为0.1，中心受集中力为400，四个角点为简支，求薄板的位移场。
分析：根据题目要求，可以选用四节点板单元。
对题目模型分网如下：图1：分网图单元大小为2.5×2.5，由于单元仅在中心处受集中载荷，结构是对称的，因此在计算时只需对其四分之一进行计算，然后再将计算结果扩充到整个板就可以了。
源程序如下：%用四节点平面板单元求解平板中心受集中载荷问题%------------------------------------% 输入控制参数%------------------------------------clear %清除workplace残留数据nel=4; % 单元数nnel=4; % 每个单元的节点数ndof=3; % 每个节点的自由度数nnode=9; % 系统总节点数sdof=nnode*ndof; % 系统总自由度数edof=nnel*ndof; % 每个单元的自由度数emodule=206e9; % 杨氏弹性模量poisson=0.3; % 泊松比t=0.1; % 薄板厚度nglxb=2; nglyb=2; % 弯曲对应的2x2高斯拉格朗日积分nglb=nglxb*nglyb; % 弯曲对应的每个单元的高斯积分点数nglxs=1; nglys=1; % 剪切对应的1x1高斯拉格朗日积分ngls=nglxs*nglys; % 剪切对应的每个单元的高斯积分点数%---------------------------------------------% 输入节点坐标值% gcoord(i,j) i：节点号 j：x,y值%---------------------------------------------gcoord=[0.0 0.0; 2.5 0.0; 5.0 0.0;0.0 2.5; 2.5 2.5; 5.0 2.5;0.0 5.0; 2.5 5.0; 5.0 5.0];%---------------------------------------------------------% 每个单元对应的节点号（逆时针排列）% nodes(i,j) i:节点号 j:对应的单元号%---------------------------------------------------------nodes=[1 2 5 4; 2 3 6 5; 4 5 8 7; 5 6 9 8];%-------------------------------------% 输入边界条件%-------------------------------------bcdof=[1 2 3 4 6 7 9 11 12 16 20 21 23 25 26]; % 约束的自由度bcval=zeros(1,15); % 对应的值%----------------------------------------------% 初始化矩阵和矢量%----------------------------------------------ff=zeros(sdof,1); % 载荷矢量kk=zeros(sdof,sdof); % 系统刚度矩阵disp=zeros(sdof,1); % 系统位移矢量index=zeros(edof,1); % 每个单元所包含的自由度kinmtpb=zeros(3,edof); % 弯曲几何函数矩阵matmtpb=zeros(3,3); % 弯曲材料系数矩阵kinmtps=zeros(2,edof); % 剪切几何函数矩阵matmtps=zeros(2,2); % 剪切材料系数矩阵%----------------------------% 载荷矢量%----------------------------ff(27)=100; % 结点9所受的集中载荷%-----------------------------------------------------------------% 单元刚度矩阵计算及其组合%-----------------------------------------------------------------%% 弯曲相关计算%[pointb,weightb]=swp2(nglxb,nglyb); % 积分点和权系数matmtpb=sbm(emodule,poisson)*t^3/12; %弯曲材料系数%% 剪切相关计算%[points,weights]=swp2(nglxs,nglys); % 积分点和权系数shearm=0.5*emodule/(1.0+poisson); % 剪切模量shcof=5/6; % 剪切修正因数matmtps=shearm*shcof*t*[1 0; 0 1]; % 剪切材料系数矩阵for iel=1:nel % 对所有单元数的循环for i=1:nnelnd(i)=nodes(iel,i); % 当前单元对应的节点xcoord(i)=gcoord(nd(i),1); % 节点对应的x坐标值ycoord(i)=gcoord(nd(i),2); % 节点对应的y坐标值endk=zeros(edof,edof); % 初始化单元刚度矩阵kb=zeros(edof,edof); % 初始化弯曲刚度矩阵ks=zeros(edof,edof); % 初始化剪切刚度矩阵%------------------------------------------------------% 弯曲相关计算%------------------------------------------------------for intx=1:nglxbx=pointb(intx,1); % x轴积分点坐标wtx=weightb(intx,1); % 权系数for inty=1:nglyby=pointb(inty,2); % y轴积分点坐标wty=weightb(inty,2) ; % 权系数[shape,dhdr,dhds]=ssf(x,y); % 计算形函数和对其相应的求导jacob2=sjc(nnel,dhdr,dhds,xcoord,ycoord); % 计算雅可比行列式detjacob=det(jacob2); % 计算雅可比行列式的值invjacob=inv(jacob2); % 求雅可比行列式的逆[dhdx,dhdy]=sxy(nnel,dhdr,dhds,invjacob); % 计算ddhdr,dhds在迪卡尔坐标下的值kinmtpb=sbB(nnel,dhdx,dhdy); % 计算弯曲几何函数矩阵%--------------------------------------------% 计算弯曲刚度矩阵%--------------------------------------------kb=kb+kinmtpb'*matmtpb*kinmtpb*wtx*wty*detjacob;endend % 结束弯曲刚度矩阵的计算%------------------------------------------------------% 剪切相关计算%------------------------------------------------------for intx=1:nglxsx=points(intx,1); % x轴积分点坐标wtx=weights(intx,1); % 权系数for inty=1:nglysy=points(inty,2); % y轴积分点坐标wty=weights(inty,2) ; % 权系数[shape,dhdr,dhds]=ssf(x,y); % 计算形函数和对其相应的教学求导jacob2=sjc(nnel,dhdr,dhds,xcoord,ycoord); % 计算雅可比行列式detjacob=det(jacob2); % 计算雅可比行列式的值invjacob=inv(jacob2); % 求雅可比行列式的逆[dhdx,dhdy]=sxy(nnel,dhdr,dhds,invjacob); % 计算dhdr,dhds在迪卡尔坐标下的值kinmtps=ssB(nnel,dhdx,dhdy,shape); % 计算剪切几何函数矩阵%----------------------------------------% 计算剪切刚度矩阵%----------------------------------------ks=ks+kinmtps'*matmtps*kinmtps*wtx*wty*detjacob;endend % 结束剪切刚度矩阵的计算%--------------------------------% 计算单元刚度矩阵%--------------------------------k=kb+ks;index=etsd(nd,nnel,ndof);% 单元对应的系统自由度号kk=ask(kk,k,index); % 合成系统刚度矩阵end%-----------------------------% 加载边界条件%-----------------------------[kk,ff]=dbc(kk,ff,bcdof,bcval);%----------------------------% 求解%----------------------------disp=kk\ff;num=1:1:sdof;nodedisp=[num' disp] % 输出节点位移%----------------------------% 后处理%----------------------------result=zeros(25,3); %初始化displace(75)=0; %将输出的结果从5x5的四分之一板扩充到10x10的全板a=re(1,0,disp);a(75)=0;displace=displace+a;a=re(10,6,disp);a(75)=0;displace=displace+a;a=re(19,12,disp);a(75)=0;displace=displace+a;for i=1:15;displace(45+i)=displace(15+i);displace(60+i)=displace(i);end[result]=dtxy(displace); %将节点位移以节点顺序输出，i->节点号，j->对应位移exgcoord=excoord; %将节点坐标从5x5的四分之一板扩充到10x10的全板[finresult]=agr(exgcoord,result); %将全板的节点坐标和节点位移对应起来，i->节点号，j->对应坐标和位移for i=1:25;finresultin(1,i)=sqrt(finresult(3,i)^2+finresult(4,i)^2+finresult(5,i)^2); %求节点位移endZ=arrayfin(finresultin); %排列节点位移[X,Y]=meshgrid(0:2.5:10,0:2.5:10);surf(X,Y,Z); %画板变形图以下为子程序：function [point2,weight2]=swp2(nglx,ngly)%-------------------------------------------------------------------% 目的:% 求二维高斯积分的积分点和权系数%% 变量：% nglx - x轴高斯积分点数% ngly - y轴高斯积分点数% point2 - 高斯积分点坐标% weight2 - 权系数%-------------------------------------------------------------------% 确定x,y轴最大的积分点数if nglx > nglyngl=nglx;elsengl=ngly;end% 初始化point2=zeros(ngl,2);weight2=zeros(ngl,2);% 求出相应的积分点坐标和权系数[pointx,weightx]=swp1(nglx);[pointy,weighty]=swp1(ngly);% 二维积分for intx=1:nglxpoint2(intx,1)=pointx(intx);weight2(intx,1)=weightx(intx);endfor inty=1:nglypoint2(inty,2)=pointy(inty);weight2(inty,2)=weighty(inty);endfunction [matmtrx]=sbm(elastic,poisson)%------------------------------------------------------------------------% 目的:% 弯曲材料系数%% 变量:% elastic - 弹性模量% poisson - 泊松比%------------------------------------------------------------------------ matmtrx= elastic/(1-poisson*poisson)* ...[1 poisson 0; ...poisson 1 0; ...0 0 (1-poisson)/2];function [shapeq4,dhdrq4,dhdsq4]=ssf(rvalue,svalue)%------------------------------------------------------------------------ % 目的：% 在自然坐标下计算形函数和对其相应的求导%% 变量：% shapeq4 - 四节点单元形函数% dhdrq4 - 形函数对r求导% dhdsq4 - 形函数对s求导% rvalue - 对应点的r坐标值% svalue - 对应点的s坐标值%% 说明：% 第一个点自然坐标(-1,-1), 第二个点自然坐标 (1,-1) % 第三个点自然坐标(1,1), 第四个点自然坐标(-1,1)%------------------------------------------------------------------------ % 形函数shapeq4(1)=0.25*(1-rvalue)*(1-svalue);shapeq4(2)=0.25*(1+rvalue)*(1-svalue);shapeq4(3)=0.25*(1+rvalue)*(1+svalue);shapeq4(4)=0.25*(1-rvalue)*(1+svalue);% 相应的求导dhdrq4(1)=-0.25*(1-svalue);dhdrq4(2)=0.25*(1-svalue);dhdrq4(3)=0.25*(1+svalue);dhdrq4(4)=-0.25*(1+svalue);dhdsq4(1)=-0.25*(1-rvalue);dhdsq4(2)=-0.25*(1+rvalue);dhdsq4(3)=0.25*(1+rvalue);dhdsq4(4)=0.25*(1-rvalue);function [jacob2]=sjc(nnel,dhdr,dhds,xcoord,ycoord)%------------------------------------------------------------------------ % 目的：% 求二维雅可比行列式% 变量：% jacob2 - 二维雅可比行列式% nnel - 每个单元节点数% dhdr - 自然坐标r对形函数的求导% dhds - 自然坐标s对形函数的求导% xcoord - 节点的x坐标值% ycoord - 节点的y坐标值%------------------------------------------------------------------------ jacob2=zeros(2,2);for i=1:nneljacob2(1,1)=jacob2(1,1)+dhdr(i)*xcoord(i);jacob2(1,2)=jacob2(1,2)+dhdr(i)*ycoord(i);jacob2(2,1)=jacob2(2,1)+dhds(i)*xcoord(i);jacob2(2,2)=jacob2(2,2)+dhds(i)*ycoord(i);endfunction [dhdx,dhdy]=sxy2(nnel,dhdr,dhds,invjacob)%------------------------------------------------------------------------ % 目的：% 求dhdr,dhds在迪卡尔坐标下的值%% 变量：% dhdx - 形函数对迪卡尔坐标x的求导% dhdy - 形函数对迪卡尔坐标y的求导% nnel - 每个单元节点数% dhdr - 形函数对自然坐标r的求导% dhds - 形函数对自然坐标s的求导% invjacob - 求雅可比行列式的逆%------------------------------------------------------------------------for i=1:nneldhdx(i)=invjacob(1,1)*dhdr(i)+invjacob(1,2)*dhds(i);dhdy(i)=invjacob(2,1)*dhdr(i)+invjacob(2,2)*dhds(i);endfunction [kinmtpb]=sbB(nnel,dhdx,dhdy)%--------------------------------------------------------------------------% 目的:% 计算弯曲几何函数矩阵%%% 变量:% nnel - 每个单元节点数% dhdx - 形函数对迪卡尔坐标x的求导% dhdy - 形函数对迪卡尔坐标y的求导%--------------------------------------------------------------------------for i=1:nneli1=(i-1)*3+1;i2=i1+1;i3=i2+1;kinmtpb(1,i1)=dhdx(i);kinmtpb(2,i2)=dhdy(i);kinmtpb(3,i1)=dhdy(i);kinmtpb(3,i2)=dhdx(i);kinmtpb(3,i3)=0;endfunction [kinmtps]=ssB(nnel,dhdx,dhdy,shape)%------------------------------------------------------------------------ % 目的：% 计算剪切几何函数矩阵%% 变量：% nnel - 每个单元节点数% dhdx - 形函数对迪卡尔坐标x的求导% dhdy - 形函数对迪卡尔坐标y的求导% shape - 形函数%------------------------------------------------------------------------for i=1:nneli1=(i-1)*3+1;i2=i1+1;i3=i2+1;kinmtps(1,i1)=-shape(i);kinmtps(1,i3)=dhdx(i);kinmtps(2,i2)=-shape(i);kinmtps(2,i3)=dhdy(i);endfunction [index]=etsd(nd,nnel,ndof)%----------------------------------------------------------% 目的:% 单元对应的系统自由度号%% 变量：% index - 和单元"iel"对应的系统自由度数% iel - 要确定系统自由度的单元% nnel - 每个单元的节点数% ndof - 每个节点的自由度%-----------------------------------------------------------edof = nnel*ndof;k=0;for i=1:nnelstart = (nd(i)-1)*ndof;for j=1:ndofk=k+1;index(k)=start+j;endendfunction [kk]=ask(kk,k,index)%----------------------------------------------------------% 目的:% 合成系统刚度矩阵%% 变量:% kk - 系统刚度矩阵% k - 单元刚度矩阵% index - d.o.f. vector associated with an element %-----------------------------------------------------------edof = length(index);for i=1:edofii=index(i);for j=1:edofjj=index(j);kk(ii,jj)=kk(ii,jj)+k(i,j);endendfunction [kk,ff]=dbc2(kk,ff,bcdof,bcval)%----------------------------------------------------------% Purpose:% 为方程[kk]{x}={ff}确定边界条件%% 变量:% kk - 系统刚度矩阵% ff - 载荷% bcdof - 边界节点对应的自由度% bcval - 边界节点对应的自由度的值%-----------------------------------------------------------n=length(bcdof);sdof=size(kk);for i=1:nc=bcdof(i);for j=1:sdofkk(c,j)=0;endkk(c,c)=1;ff(c)=bcval(i);endfunction [u]=re(a,b,disp)%------------------------------------------------------------------------% 目的：% 从5x5的四分之一板扩充输出结果到10x10的全板%% 变量：% disp - 扩充前的节点位移矩阵% u - 扩充后的节点位移矩阵%------------------------------------------------------------------------for i=a:1:a+8;u(1,i+b)=disp(i);endfor i=a+9:1:a+11;u(1,i+b)=disp(i-6);endfor i=a+12:1:a+14;u(1,i+b)=disp(i-12);endfunction result=dtxy(displace)%------------------------------------------------------------------------% 目的：% 将节点位移以节点顺序输出，i->节点号，j->对应位移%% 变量：% result - 按节点顺序输出后的矩阵% displace - 按自由度顺序输出的矩阵%------------------------------------------------------------------------for i=1:25;if i==1;result(1,1)=displace(i);elseresult(i,1)=displace(i*3-2);endif i==1;result(1,2)=displace(2);elseresult(i,2)=displace(i*3-1);endif i==1;result(1,3)=displace(3);elseresult(i,3)=displace(3*i);endendfunction exgcoord=excoord%------------------------------------------------------------------------ % 目的：% 将节点坐标从5x5的四分之一板扩充到10x10的全板%% 变量：% exgcoord - 扩充后的节点坐标矩阵%------------------------------------------------------------------------x=0:2.5:10;y=0:2.5:10;exgcoord=zeros(25,2);for i=1:5;exgcoord(i,1)=x(i);exgcoord(i,2)=y(1);endfor i=6:10exgcoord(i,1)=x(i-5);exgcoord(i,2)=y(2);endfor i=11:15exgcoord(i,1)=x(i-10);exgcoord(i,2)=y(3);endfor i=16:20exgcoord(i,1)=x(i-15);exgcoord(i,2)=y(4);endfor i=21:25exgcoord(i,1)=x(i-20);exgcoord(i,2)=y(5);endfunction finresult=agr(exgcoord,result)%------------------------------------------------------------------------% 目的：% 将全板的节点坐标和节点位移对应起来，i->节点号，j->对应坐标和位移%% 变量：% finresult - 节点坐标和对应位移对应起来的矩阵% result - 节点位移矩阵% exgcoord - 节点坐标矩阵%------------------------------------------------------------------------finresult=zeros(5,25);for i=1:5;if i<=2;finresult(i,:)=exgcoord(:,i);elsefinresult(i,:)=result(:,i-2);endend%------------------------------------------------------------------------% 目的：% 排列节点位移%% 变量：% finresultin - 节点坐标和对应位移对应起来的矩阵% Z - 按节点顺序排列的节点位移矩阵%------------------------------------------------------------------------for i=1:25;if i<=5;Z(1,i)=finresultin(i);elseif i>5&i<=10;Z(2,i-5)=finresultin(i);elseif i>10&i<=15;Z(3,i-10)=finresultin(i);elseif i>15&i<=20;Z(4,i-15)=finresultin(i);elseif i>20&i<=25;Z(5,i-20)=finresultin(i);endend计算结果输出如下。
有限单元法
弹塑性有限单元法汽车车体冲压件通常使用弹塑性材料，在冲压过程中，这种材料的变形成形过程非常复杂，一般用刚塑性FEM与弹塑性FEM二种方法来评价整个冲压成形过程。
在刚塑性FEM中，忽略弹性变形，仅将塑性应变作为计算指标。
因此，在冲压成形过程中，当材料放置到模具上因自重产生的弯曲挠度，从模具中取出冲压件厚产生的弹性恢复等材料变形不能进行计算。
因此，有人提出了根据刚塑性FEM的计算结果，再用弹性FEM计算其卸载过程，但是，刚塑性FEM很难正确地预测在冲压过程中产生的缺陷。
弹塑性FEM可以再空间上时间上交替考虑弹性变形与塑性变形，从理论上讲可以正确地描述整个冲压过程，所以弹塑性FEM可以说是评价冲压过程的最好解析方法。
现有的弹塑性FEM，根据其时间积分方法的不同，可分为“静态显函数法”“静态隐函数法”和“动态显函数法”。
讲加速度项加入平衡方程式求解的称为动态，反之，平衡方程式中不包含加速度项的解法称为静态。
隐函数与显函数是常微分方程数值计算方法中的数学用语。
显函数求方程的解不需要反复计算，而隐函数常微分方程求解时需要迭代多次逼近其解。
显函数解法要求增分补偿不能取得太大，解析冲压成形过程需要较多的计算解析次数。
隐函数解法通常可以保证应力平衡方程式成立，因而增分步长可以取得较大些，以减少解析计算次数。
各种弹塑性FEM的优缺点如下：动态显函数法：该方法求解各节点的独立性运动方程以获得节点变形，因而不需要组成刚度矩阵，即使单元划分得再细，节点再多也占用的计算内存较少，并且每一模拟步骤的计算速度也比其他方法快，因此可以计算对象的单元分割得很细。
但是这种方法是用动态的冲击求解变形问题，时间增量需控制在10-6秒以下，要模拟一秒钟的冲压过程，就需要计算10^6次，实际上这种计算方法十分耗时，为减少运算时间，常常将物理意义不十分清楚的衰减项加入到运动过程，人为地将质量附以加权常数以减少模拟计算次数。
此外，即使在方程式中加入了衰减项，应力值还是会发生振动，增加了弹性恢复计算的难度。
有限单元法基本步骤示例
有限单元法在电磁 场分析中用于求解 电磁场方程
通过对电磁场进行 离散化处理，将连 续的电磁场转换为 离散的有限单元
通过对有限单元进 行数学建模和求解 ，得到电磁场的分 布情况和相关参数
有限单元法在电磁 场分析中具有广泛 的应用，如电磁场 仿真、天线设计、 电磁兼容性分析等
06 总结与展望
总结有限单元法的优势和不足
关系等。
最小势能原理
定义：最小势能原理是指在物理系统中，系统的总势能总是趋向于最小 值 应用：有限单元法中，通过最小化总势能来求解物理问题
优势：能够考虑系统的约束条件，得到精确解
局限性：对于非线性问题，可能会出现求解困难的情况
虚功原理
定义：虚功原理是有限单元法的基本原则之一，它指出在结构分析中，如果结 构受到的载荷是虚的（即不真实的），则结构的位移和应力也将是虚的。
优势：适用于复杂形状 和边界条件的离散化， 能够解决各种工程问题。
不足：计算量大，需要 高性能计算机支持；对 初学者来说，掌握难度 较大。
展望有限单元法未来的发展方向和应用前景
研究方向：随着科技的不断进步，有限单元法在理论和应用方面将会有更多的突破和创新，例 如开发更加高效、精确的算法和模型，以解决更加复杂的问题。
进。
05 有限单元法的应用实例
结构分析中的应用
有限单元法在结构分析中用于建立离散化模型，将连续的结构离散为有限 个单元，以便进行数值计算和分析。
有限单元法在结构分析中可以模拟各种复杂的结构和边界条件，例如桥梁、 高层建筑和核反应堆等。
有限单元法在结构分析中可以用于评估结构的强度、刚度和稳定性等性能 指标，为结构设计提供依据。
应用领域：随着工业和科技的不断发展，有限单元法将会被应用到更多的领域中，例如航空航 天、汽车、建筑、生物医学等，解决各种复杂的问题和挑战。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验课程名称
弹性有限元法及应用
开课实验室
DS1442 材料学院 年级 08 号 专业班 材控 1 班
学
院
学 生 姓 名
杨 阳 2009 2010
学
20084062 2 学期
开 课 时 间
至
ቤተ መጻሕፍቲ ባይዱ
学年第
总 成 绩 教师签名
材料学院制 材料学院制
《有限单元法上机实验》实验报告 有限单元法上机实验》
开课实验室： DS1442 开课实验室： 学院 课程 有限单元法实验 名称 教 师 评 教师签名： 语 年 一、实验目的 1、掌握有限单元法分析原理；知道如何进行有限元离散、单元分析和整体分析，并得出 月 日 名 称 材料学院 年级、专业、班 实验项目 I-deas 的熟悉与初步使用 指导教师 徐戊娇 08 级材控 1 班 姓名 2010 年 杨 阳 成绩 5 月 25 日
五、实验过程原始记录(数据、图表、计算等) 实验过程原始记录(数据、图表、计算等) 模拟载荷实验过程的一些截图： 模拟载荷实验过程的一些截图： 实验过程的一些截图
模型构建的一些截图：
六、实验结果及分析 悬臂梁模拟载荷，有限元分析结果： 悬臂梁模拟载荷，有限元分析结果：
模型建立，实例构建结果： 模型建立，实例构建结果：
三、使用仪器、材料 使用仪器、 EDS 下 I-DEAS10 软件
四、实验步骤 （一）对 I-DEAS 软件进行初步介绍: 1.介绍软件界面的不同分区，包含图形区、执行命令输出信息区、提示区和功能按钮区。 2.介绍软件的不同模块，主要介绍设计和模拟分析模块可以实现的基本功能。 （二）学习 CAD 的基本知识和 3D 建模的基本步骤，为 CAE 分析做好模型准备 （三）有限元分析的一般格式： 1.对结构进行离散化，选择恰当的单元和适当的单元尺寸划分结构网络； 2.定义性函数 N，进行通过单元节点量描述单元域内连续的变量； 3.行程单元矩阵，建立单元与外界环境的平衡关系； 4.集成，将覆盖结构全域的所有单元的刚度矩阵和节点外力对平衡的贡献集成，建立整体结构 的平衡方程； 5.求解平衡方程，求解已知位移的节点反作用力和未知的节点位移； 6.回带，根据计算出的节点变量，带入第二步的表达式中，获得单元应变、应力或热流等量。 (四)用 I-DEAS 软件进行构造有限元分析的步骤： 1.运用软件的设计模块构建所分析结构的 3D 实体模型，为结构 FEM 做准备； 2.前处理部分： *划分网络，确定单元类型和单元的结构尺寸； *构建边界条件集，包含载荷集、约束集、接触集和温度分布等结构； *确定求解类型，是线性静态问题、非线性静态或动态问题等。 3.模型求解； 4.后处理部分：显示求解结果，分析结构的位移、应力和应变分布云图。
相应的分析结果。
2、学会用 I-DEAS 建模并进行有限元的分析。通过上机练习，让学生学习 3D 实体建模的
基本知识。了解有限元法作为一种数值计算方法解决工程实际问题的全过程。
二、实验原理 用 I-deas 软件进行结构偶有限元分析的步骤 前处理程序：前处理程序的功能是根据一个实际问题的物理模型，用尽可能接近自然语言的方式，向 计算机输入尽可能少的定义有限元模型和控制分析过程的数据， 自动生成有限元分析主体程序所需要 的全部输入文件。 几何造型：几何造型是指在选定的坐标（常用的是直角坐标系、圆柱坐标系和球坐标系）内通过几何 元素（点、线、面、体）的生成和对他们的编辑，生成有限元分析对象的几何构造和图形。 网格生成： 1. 转换生成法：这是直接将几何元素的点、线、面和体直接转化成有限元的结点、线单元、面单元 和体单元。 2. 自动生成法：这是在商业软件中常被采用的方法。他在任意形状的平面或空间裁剪曲面上可生成 三角形或四边形单元;对单元几何实体或由曲面围成的封闭空间，可生成四面体或六面体实体单 元。 后处理程序： 后处理程序的功能是对用户在前处理程序中指明需要输出的计算结果进一步处理和图形 显示。主要计算结果中通常需要处理的物理量是位移（矢量或各个分量） 、应力（等效应力或各个分 量） 和温度等。 为了清晰地观察变形图和未变形图对比的显示结果， 变形图可以由用户给出放大倍数。