有限元平面钢架.doc

平面钢架有限元作业2017-7-51401S205 付小龙1401s207 胡新帅1401s208 贾伟波目录第一章 问题重述3一、题目内容：3二、题目要求：3第二章 原理说明4一、平面刚架的解题思路（步骤）：4二、单元刚度矩阵的建立4三、总体刚度矩阵9四、位移的求解9第三章 手工算法求解11一、单元离散化11二、单元分析11三、单元组装13四、边界条件引入及组装总体方程14五、求解整体刚度方程，计算节点2的位移和转角14六、求节点1、3支撑反力15七、设定数据，求解结果15八、轴力图、弯矩图、剪力图的绘制16第四章 matlab编程计算18一、流程图18二、输入数据18三、计算单元刚度矩阵18四、建立总体刚度矩阵：19五、计算未约束点位移19六、计算支反力19七、输出数据19八、编程19第五章 ANSYS分析部分20一、参数预处理设置20二、建立数学模型21三、添加载荷和约束24四、分析求解27五、绘制图像29第六章 结果比较32第七章 心得体会33胡新帅33付小龙33贾伟波33附录35一、matlab原程序35二、matlab运行结果36第一章 问题重述一、题目内容：图示平面钢架结构图1.1 题目内容二、题目要求：（1）采用平面梁单元进行有限元法手工求解，要求写出完整的求解步骤，包括： a）离散化：单元编号、节点编号；b）单元分析：单元刚度矩阵，单元节点等效载荷向量；c）单元组装：总体刚度矩阵，总体位移向量，总体节点等效载荷；d）边界条件的引入及总体刚度方程的求解；e）B点的位移，A、C处支撑反力，并绘制该结构的弯矩图、剪力图和轴力图。（2）编制通用平面钢架分析有限元Matlab程序，并计算盖提，与手工结果进行比较；（3）利用Ansys求解，表格列出B点的位移，A、C处支反力，绘制弯矩图、剪力图和轴力图，并与手算和Matlab程序计算结果比较。（4）攥写报告，利用A4纸打印；（5）心得体会，并简要说明各成员主要负责完成的工作。第二章 原理说明一、平面刚架的解题思路（步骤）：（1）结构离散化：将每个杆件作为一个单元，杆件与杆件之间的联结点为节点，并进行编码。（2）单元刚度矩阵的建立：具体内容有：建立局部坐标系、单元刚度矩阵的导出、坐标转换、单元刚度矩阵的坐标变换、单元刚度矩阵的特性。（3）求总体刚度矩阵：总体刚度矩阵的分块形式、总体刚度矩阵的集成、总体刚度矩阵的特性。（4）位移的求解：支撑条件的引入、非节点载荷的处理。（5）求应变、应力二、单元刚度矩阵的建立（1）局部坐标系与整体坐标系：由于各杆件的取向不同，所以对各杆单元必须先建立各自的局部坐标系。以节点1(e)为原点，由1(e)到2(e)的方向x轴正向的右手坐标系就叫做该单元的局部坐标系，用英文小写字母x和y表示。如图2-1所示。图2-1 一维单元（钢架杆单元）为便于进行整体分析，要求整个结构的所有杆件都采用同一个坐标系，如图所示。该坐标系就称为整体坐标系，用英文大写字母X和Y表示：在局部坐标系中，杆单元每个节点都有三个位移分量，如图(a)所示。每个节点都有三个力分量，如图（b）所示图中 轴向位移、轴向力 横向位移、横向力 角位移、弯矩图中所示各力的方向为正。每个单元在两端共有六个节点位移和六个节点力分量（2）单元刚度矩阵的导出：对于刚架杆单元，由于位移和力向量都不再是二维而是六维，所以单元刚度矩阵也不再是22阶而是66阶，其基本形式为：求单元刚度矩阵的核心是求各单元刚度系数。下面求单元刚度矩阵。的第一列对应于只有节点1(e)在轴向方向产生单位位移时的情况，如下图所示。这时除了外，两节点的其余位移分量都为零。要产生这个变形并能保持该杆单元的平衡，只需要在两个节点处施加两个轴向力他们分别为刚度系数。很明显，第一列的其它元素均为零。的第二列对应于只有节点1(e)在横向方向产生单位位移时的情况，如下图所示。这时除了外，两节点的其余位移分量都为零。它相当于自由端绕度为1、转角为0的一根悬臂梁。要产生这个变形并能保持该梁的平衡，在节点1(e)处要施加横向力在节点2(e)处要施加横向力和弯矩，他们分别为刚度系数、。同样。的第三列对应于只有节点1(e)在产生单位角移时的情况，如下图所示。这时除了外，两节点的其余位移分量都为零。它相当于一根悬臂梁，自由端有铰链支撑，在节点1(e)处产生单位转角。要产生这个变形并能保持该梁的平衡，在节点1(e)处要施加横向力：在节点2(e)处要施加横向力和弯矩： 他们分别为刚度系数、。同样。中的后三列元素也可以用同样的方法导出。局部坐标系中杆单元(e)的节点力和节点位移之间的关系式为：单元刚度矩阵是一个66阶的矩阵（3）坐标变换：前面推导时采用的是局部坐标系，在研究刚架或桁架结构时，由于各杆的取向不同，存在若干个局部坐标系，这就给整体分析带来了麻烦。因此必须通过坐标变幻，将在局部坐标系中获得的单元刚度矩阵转换成整体坐标系中的单元刚度矩阵。假设一个杆单元两个节点的局部码1(e)和2(e)所对应的总码分别为I和j，则可写出一个坐标转换的矩阵方程（过程略）：式中单元（e）的坐标转换矩阵，它是一个正交矩阵，即其逆阵等于其转置阵 所以上式又可写为 局部坐标系x、y和整体坐标系X、Y两者之间相差的角度。规定由X轴到x轴以逆时针方向为正。单元（e）在局部坐标系中按节点局部码标记的节点位移向量单元（e）在整体坐标系中按节点局部码标记的节点位移向量对于单元节点力向量也存在上述关系，即式中分别为与相对应的节点力向量（4）单元刚度矩阵的坐标变换整体坐标系中单元节点力和节点位移之间的关系式式中整体坐标系中单元(e)的刚度矩阵。因为所以。三、总体刚度矩阵前面首先推导出了局部坐标系中的单元刚度矩阵，然后通过坐标转换获得了整体坐标系中的单元刚度矩阵。在得到所有单元在整体坐标系中的单元刚度矩阵的基础上，可以建立总体刚度矩阵。（1） 总体刚度矩阵的集成由整体坐标系中的单元刚度矩阵的子矩阵集成总体刚度矩阵的步骤如下：a） 对一个有n个节点的结构，将总体刚度矩阵K划分为nn各子区间，然后按节点总码的顺序进行编号。b） 将整体坐标系中单元刚度矩阵的各子矩阵根据其下标的两个总码对号入座，写在总体刚度矩阵相应的子区间。c） 同一子区间内的子矩阵相加，成为总体刚度矩阵中的相应的子矩阵。（2） 总体刚度矩阵的特性a) 对称性：因为由此特性，在计算机中只需存储其上三角部分.b) 奇异性：物理意义仍为在无约束的情况下，整个结构可做刚体运动。c) 稀疏性：K中有许多零子矩阵，而且在非零子矩阵中还有大量的零元素，这种矩阵称为稀疏矩阵。大型结构的总体刚度矩阵一般都是稀疏矩阵。d) 分块性：这个性质已经利用过，在此不再叙述。四、位移的求解从前面例子中看出获得总体刚度矩阵后，即使已知节点力向量，仍然不能求出节点位移向量，因为K是一个奇异矩阵。为此，必须引入支撑条件。另一方面，与前面例子不同的是对刚架结构的杆单元，载荷不仅可以施加在节点上，还可能施加在单元上，这种载荷称为非节点载荷。对非节点载荷，必须将他们移植到节点上才能获得已知的节点力向量。所以，在此介绍支撑条件的引入以及非节点载荷的处理。（1） 支撑条件的引入刚架结构的特点，从支撑条件的角度可以分为两类：一类是在支撑处的位移向量为零；另一类是在无支撑处，其位移向量未知。例如一个有三个节点的钢架结构，节点1和节点2是在无支撑处，其节点位移向量分别为，节点3是在支撑处，其位移向量。其总体平衡方程可写为：可简化成：求出后代入可求得节点3处的支撑反力。（2） 非节点载荷的处理在刚架结构以及其他较复杂的结构上，他们所受的载荷可以直接作用在节点上，又可以不直接作用在节点上而作用于单元节点间的其他位置上。后一种情况下的载荷称为非节点载荷。有限元分析时，总体刚度方程中所用到的力向量 是节点力向量。因此在进行整体分析前应当进行载荷的移植，将作用于单元上的力移植到节点上。移植时按静力等效的原则进行。处理非节点载荷一般可直接在整体坐标系内进行，其过程为：a) 将各杆单元看成一根两端固定的梁，分别求出两个固定端的约束反力。其结果可直接利用材料力学的公式求得。b) 将各固定端的约束反力变号，按节点进行集成，获得各节点的等效载荷。第三章 手工算法求解一、单元离散化将平面梁离散为两个单元，单元编号分别为和，节点号分别为1、2、3；如图3-1所示：图3-1单元离散化示意图二、单元分析首先建立整体坐标系与局部坐标系如图3-1所示；1、求单元刚度矩阵单元局部坐标系的单元刚度矩阵为：由于单元局部坐标系与整体坐标系的夹角为：，则单元的局部坐标变换矩阵为：因此在总体坐标系下的单元的刚度矩阵为：单元局部坐标系的单元刚度矩阵为：由于单元局部坐标系与整体坐标系的夹角为,则。2、求单元节点等效载荷向量将P等效于单元两侧节点1,2上：将均布载荷等效在单元两侧的节点2，3上：与作用在节点上的力叠加为整体坐标系下的节点载荷：三、单元组装将两个整体坐标系下的单元刚度矩阵组装成的整体刚度矩阵如下：四、边界条件引入及组装总体方程由于节点1、3为固定约束，所以节点1和3的x、y方向的位移以及转角均为0，节点2无位移约束，不存在支反力，所以力约束即为外力约束。 五、求解整体刚度方程，计算节点2的位移和转角提取节点2位移的相关要素：求得：六、求节点1、3支撑反力根据总体方程，提取求解节点1支承反力所需方程：根据总体方程，提取求解节点3支承反力所需方程：七、设定数据，求解结果设定各个数据：杨氏模量：泊松比：力：截面面积：惯性矩：将以上数据代入各式：节点2的位移和转角：节点1支撑反力：节点3支撑反力：八、轴力图、弯矩图、剪力图的绘制轴力图：图3-2 轴力图剪力图：图3-3 剪力图弯矩图图3-4 弯矩图第四章 matlab编程计算一、流程图图 4.1 总体流程图二、输入数据输入每个单元的杨氏模量，惯性矩，单元长度，单元截面积以及单元的旋转角度并采用矩阵的形式进行输入。三、计算单元刚度矩阵图3.2 单元刚度矩阵生成流程图考虑到每个单元的刚度矩阵与坐标变换的矩阵形式相同，只是数据不同，故采取建立模板，利用subs()，函数来带入不同单元的值，生成一系列单元刚度矩阵，并用一个三维数组存储这些矩阵。四、建立总体刚度矩阵：考虑到每个单元刚度矩阵都是66的形式，表述了2个节点间的相互关系；故建立元胞数组，并使元胞数组的阶数与节点个数相同，利用元胞数组存储节点间关系。首先建立与节点个数相同阶数的空元胞数组，之后检索每个单元刚度矩阵对应的2个节点间的关系，将其分离成4个33的矩阵，按节点与单元对应关系，存储到元胞数组中。最后将元胞数组展开形成的大矩阵即为总体刚度矩阵。五、计算未约束点位移利用总体位移与外力间的关系，采用矩阵求解，求取非约束点的位移。并针对结果进行对应处理，使结果与作用点、作用形式对应。六、计算支反力利用约束点位移皆零的特点，简化总体刚度矩阵，同时由于部分节点的部分方向上为内力而非支反力，再度简化总体刚度矩阵。利用两次简化后的刚度矩阵与计算出的位移结果相乘，求得不计直接作用在节点约束方向上时的支反力，将结果加上由于直接作用在节点约束方向上时产生的支反力，即为最后的支反力结果。七、输出数据将计算所得的未约束点位移与支反力，采用与输入方式相似的方式进行处理并进行输出。见附录二。八、编程见附录一 第五章 ANSYS分析部分一、参数预处理设置选择单元类型ANSYS Main Menu: Preprocessor Element TypeAdd/Edit/Delete Add beam：2D elastic 3OK (返回到Element Types 窗口) Close2、定义材料参数：ANSYS Main Menu: Preprocessor Material Props Material ModelsStructural Linear ElasticIsotropic: EX:3e10 (弹性模量)，PRXY:0.3(泊松比)OK3、定义单元截面积和惯性矩：ANSYS Main Menu: Preprocessor Real constant Add Type beam 3 Ok Cross-sectional area AREA:0.05(横截面积) Area moment of inteia IZZ:1(惯性矩) OK二、建立数学模型建立三点坐标：ANSYS Main Menu: Preprocessor Modeling CreatKeypoint In Active CSNode number： 1 X:0,Y:0,Z:0 Apply Node number ： 2 X:0,Y:1,Z: Apply Node number ： 3 X:2,Y:1,Z:0OK 三点确立后的模型：构造直线模型;ANSYS Main Menu: Preprocessor Modeling CreatLine linesstraight line 依次连接各点Ok设置所有线上划分单元的个数ANSYS Main Menu: Preprocessor Meshing Size Cntrls ManualSize Lines All Lines No.of element divisions ,10OK划分后：划分单元ANSYS Main Menu: Preprocessor Meshing Mesh Lines 单击 pick all OK注：此时从PointCtrls Number.中可以选择显示单元编号或者节点编号三、添加载荷和约束为1号关键点和3号关键点添加约束ANSYS Main Menu: Preprocessor SolutionDefine loads Apply Structural Displacement On nodes 选择1关键点ALL DOFApplyOn nodes 选择3关键点ALL DOFOK2.添加顶部均布载荷：ANSYS Main Menu: Preprocessor SolutionDefine loads Apply Structural Pressure On beams 选择顶部所有的单元VALI pressure value node I :1000 VALJ pressure value node J :1000 OK添加力矩和力：ANSYS Main Menu: Preprocessor SolutionDefine loads Apply Structural Force/ MonmentOn nodes 选择2节点Apply LAB MZVALUE 100 . (输入力矩)On nodes 选择8节点Apply LAB FX VALUE 1000 （输入力）所有载荷约束添加完成后：四、分析求解ANSYS Main Menu: Solution Solve Current LS OK Should the Solve Command be Executed? Y Close (Solution is done! ) 关闭窗口求解位移ANSYS Main Menu: General Postproc List resultNodal solutionDOF solutionX-component of displacement Apply Y-component of displacementOK2、支反力：ANSYS Main Menu: General Postproc List resultReaction SoluAll itemsOK五、绘制图像1、设置参数ANSYS Main Menu: General Postproc Element Table Difine Table Add 在user label for item 中输入FX-I, 在Results data item中选择By sequence num并输入 smisc,1Apply在user label for item 中输入FX-J, 在Results data item中选择By sequence num,并输入 smisc,7Apply在user label for item 中输入FY-I, 在Results data item中选择By sequence num,并输入 smisc,2Apply在user label for item 中输入FY-J, 在Results data item中选择By sequence num,并输入 smisc,8Apply在user label for item 中输入MZ-I, 在Results data item中选择By sequence num,并输入 smisc,6Apply在user label for item 中输入MZ-J, 在Results data item中选择By sequence num,并输入 smisc,12OK2、图像输出a）轴力图：ANSYS Main Menu: General Postproc Plot resul tContour plot Line Elem Res 选择 FX_I ，FX_J Applyb）剪力图：ANSYS Main Menu: General Postproc Plot resul tContour plot Line Elem Res 选择 FY_I FY_J Apply图4.12 剪力图c）弯矩图：ANSYS Main Menu: General Postproc Plot resul tContour plot Line Elem Res 选择 MZ_I ， MZ_J OK第六章 结果比较结果手算MATLABANSYS-0.4142210-8-0.4142210-8-0.4142210-8-0.3302310-7-0.3302310-7-0.3302310-7-1033.1-1033.1-1033.149.53549.53549.53564.50164.50164.5013.10673.10673.10671950.51950.51950.5-1462.3-1462.3-1462.3通过对比知道，三种方式的结果完全一样突出了结果的正确性。第七章 心得体会胡新帅主要负责 第二章原理说明 及 第四章matlab编程计算 的编撰工作。“以前会有限元的人只要在火车站喊一声，就会有人抢走你。”当看到这篇文章时，我特别好奇：有限元是什么？当时为什么这么奇缺？经过一学期的熏陶与学习，发现有限元是解决一类问题的通解，相对传统学科它的优势就在于对问题的限定条件更少。有限元法（Finite Element Method）是基于近代计算机的快速发展而发展起来的一种近似数值方法, 用来解决力学，数学中的带有特定边界条件的偏微分方程问题（PDE）。而这些偏微分方程是工程实践中常见的固体力学和流体力学问题的基础。有限元和计算机发展共同构成了现代计算力学的基础。所以有限元是在复杂区域（像汽车、船体结构、输油管道）上解偏微分方程的一个很好的选择，也就造成了当时社会上对有限元人才的渴求与此同时，我的matlab编程能力也得到了一定的提升。不仅熟悉了matlab的基本操作，而且学会了各种矩阵的操作，更加懂得了matlab的在具体应用中的重要性。付小龙听说有限元分析这个课程是李建宇老师在给我们上材料力学的时候提到的，李老师一直再强调有限元分析这个重要性。到了这学期选了有限元分析这门课后，才真正领会了有限元的神奇之处。在上课的时候听这个课，确实感觉很枯燥，而且有点抽象，因为涉及到矩阵，就有些抽象，听下来一遍也感觉没掌握多少，通过最后的这个大作业才又搞懂了一些，算一遍确实是有好处的，弄清楚了有限元的基本原理，用有限元分析解决了大作业这个问题，感觉这个方法用起来还是挺方便的，怪不得外界需求量这么大。贾伟波刚开始接触有限元这门课的时候，感觉接受难度还是挺大的，毕竟没有理论力学和材料力学的学习基础，虽然原理能听懂，但是概念大部分都不理解。随着课程不断深入，开始对各种矩阵运算有浅显的了解。到现在，粗略学会杆与梁的有限元分析，对板的研究也略解皮毛。说到收获，感觉最大的还是学会了ansys这款软件。相比于复杂的手算，冗长的编程，ansys不仅操作更简便，分析更全面，而且可靠性更高，出错率更低，运算更方便。运用ansys的最大难点是如何将实际问题转换成模型，在本次作业中也遇到了由于题示所给参数不全而煞是苦恼的情况，不过好在有机械系同学的帮忙，解决了参数间的关系问题。 说实话，从虚位移部分开始就基本又听不懂了。虽然高中时接触过虚功原理，但并没有做过系统，理论的分析，加上缺少大学力学基础知识，公示概念依旧晦涩难懂。这应该是学习有限元最大的缺憾之一。学完之后，印象最深的是老师每次上课前都重复的那段话、那张图：有限元基本的原理和思想。学会了节点，单元，整合，略微了解了超级复杂的弹性力学分析，总之，有限元这节课受益匪浅。附录一、matlab原程序clcclearclose allsyms fai e l a i pi p%本题中数据jd_num=3;jdzh=44*p/125 0 -12*p*l/125 , 81*p/125 , -p , -67*p*l/750 , 0 , -p , p*l/3; %整体坐标系下的节点载荷jD=1,2; 2,3; %每个单元在整体坐标下节点顺序k=e*a/l 0 0 -e*a/l 0 0; 0 12*e*i/l3 6*e*i/l2 0 -12*e*i/l3 6*e*i/l2; 0 6*e*i/l2 4*e*i/l 0 -6*e*i/l2 2*e*i/l; -e*a/l 0 0 e*a/l 0 0; 0 -12