下载文档原格式
河北工程大学《有限元及程序设计》课程设计说明书课程设计题目：平面钢架有限元程序功能修改设计副标题：1增加非节点荷载为集中力偶、纵向均布力以及线性分布力时固端反力的计算功能2 结果以文件形式输出指导教师：班级：姓名：学号：摘要有限元法是现代工程数值分析中应用广泛的一种方法，本文根据线性有限元理论对受有五种不同荷载的三杆钢架进行静力分析，将结构离散为三个单元四个节点。
分别建立每个单元的单元刚度矩阵和节点荷载列阵，根据刚度集成法建立了结构的总刚度矩阵和节点荷载列阵，得出结构的平衡方程，并用对角线元素置一法引入边界条件，用高斯消元法求解平衡方程。
最后编写C程序求解此问题，并通过与手算的比较验证了程序的准确性，通过增加一个杆件说明程序的通用性键词：有限元法杆平面刚架刚度矩阵对角线元素置一法 C语言A b s t r a c tFinite Element Method (FEM) makes an extensive use in the numerical analysis of modern construction. In this paper, we study Static state of the Plane frame which is loaded by five different kinds of loads based on linear finite element theory loads; the structure is divided into three units and four nodes. The stiffness matrix array node load of each unit is set up, According to Stiffness integration method we establish the total stiffness matrix and load node array of the structure which is aimed drawing the balance equation of the structure, and then introduce the boundary conditions by buy-one-diagonal elements .we can use Gaussian elimination method for solving equilibrium equations. Finally, we program procedures for the preparation of C to solve this problem, and through comparison with the hand count to verify the accuracy of the procedure, by adding a bar to descript generic property of the procedureKey words: Finite Element Method (FEM) Planar rigid frame Buy-one-diagonal elements Stiffness matrices C language目录设计题目说明-------------------------------------------------4 1.用有限元法进行手算----------------------------------------41.1化分单元，选取坐标系----------------------------------------------4 1.2求局部坐标系下各单元的单元刚度矩阵--------------------------------5 1.3求整体坐标系下各单元的单元刚度矩阵--------------------------------5 1.4求整体刚度矩阵----------------------------------------------------6 1.5求非节点荷载引起的等效节点荷载及节点荷载列阵----------------------6 1.6列整刚方程，求节点位移--------------------------------------------9 1.7求单元内力----- --------------------------------------------------101.7.1转换位移列阵- ---------------------------------------------------101.7.2求内力- --------------------------------------------------------111.8画内力图并列表--------------------------------------------132.程序设计与上机调试结果-------------------------------------142.1说明与结果- ------------------------------------------------------14 2.2程序设计中一些问题的描述- ----------------------------------------142.2.1数字描述--------------------------------------------------------142.2.2程序总框图------------------------------------------------------16主要结论-----------------------------------------------------16 设计心得体会---------------------------------------------------------16参考文献-----------------------------------------------------17附录C程序源代码及修改注释- --------------------------------17求图示平面刚架节点位移及各杆的内力错误!具体参数：面积弹性模量E=27/101.2mKN ⨯惯性矩I=421016667.4m-⨯1用有限元法进行手算1.1划分单元，标出单元号码及节点号码；选取整体坐标系O x y ，局部坐标系Oxy ，并标上单元的局部节点码i(1),j(2)，见下页图。
结构有限元分析程序设计绪论§0.1 开设“有限元程序设计”课程的意义和目的§0.2 课程特点§0.3 课程安排§0.4 课程要求§0.5 基本方法复习$0.1 意义和目的1.有限元数值分析技术本身要求工程设计研究人员掌握1). 有限元数值分析技术的完善标志着现代计算力学的真正成熟和实用化,已在各种力学中得到了广泛的应用。
比如：，已杨为工程结构分析中最得以收敛的技术手段，现代功用大致有：a). 现代结构论证。
对结构设计从内力，位移等方面进行优劣评定，从而进行结构优化设计。
b)可取代部份实验，局部实验+有限元分析，是现代工程设计研究方法的一大特点。
c)结构的各种功能分析（疲劳断裂，可靠性分析等）都以有限元分析工具作为核心的计算工具。
2). 有限元数值分析本身包括着理论+技术实现（本身功用所绝定的）有限元数值分析本身包括着泛函理论+分片插值函数+程序设计2. 有限元分析的技术实现（近十佘年的事）更依赖于计算机程序设计有限元分析的技术取得的巨大的成就，从某种意义上说，得益于计算机硬件技术的发展和程序设计技术的发展，这两者的依赖性在当代表现得更加突出。
（如可视化技术）3.从学习的角度，不仅要学习理论，而且要从程序设计设计角度对这些理论的技术实现有一个深入的了解，应当致力于掌握这些技术实现能力，从而开发它，发展它。
（理论本身还有待于进一步完美相应的程序设计必须去开发）4.程序设计不仅是实现有限元数值分析的工具和桥梁，而且在以下诸方面也有意义：1). 精通基本概念，深化理论认识；2). 锻炼实际工程分析，实际动手的能力；3). 获得以后工作中必备的工具。
（作业+老师给元素库）目的：通过讲述有限元程序设计的技术与技巧，便能达到自编自读的能力。
§0.2 课程特点总描述：理论+算法+数据结构（程序设计的意义）理论：有限元算法，构造，步骤，解的等外性，收敛性，稳定性，误差分析算法；指求解过程的技术方法，含两方面的含义；a. 有限元数值分析算法，b, 与数据结构有关的算法（总刚稀疏存贮，提取，节点优化编号等）数据结构：指各向量矩阵存贮管理与实现，辅助管理结构（指针，数据记录等）具体特点：理论性强：能量泛函理论+有限元构造算法+数据结构构造算法内容繁杂：理论方法+技术方法+技术技巧技巧性强：排序，管理结构（指针生成，整型运算等）§0.3 课程安排①. 单元刚度矩阵及元素设计（单元刚阵算法，杆梁平面分析，板弯非协调元等）②. 总刚的形式及程序设计（单刚提前准备，技术复杂）③. l边界条件及程序设计（等效荷载计算，位移边界条件置入，多工况的对称性）④. 总刚线性方程组求解（LDL T分解，分块算法，子结构算法，波前法）⑤．单元应力计算+应力处理与改善。
有限元程序设计课程设计一、课程目标知识目标：1. 掌握有限元分析的基本原理，理解有限元方法在工程问题中的应用。
2. 学会使用至少一种有限元分析软件，并能正确进行前处理、计算及后处理操作。
3. 掌握编写有限元程序的基本步骤，理解数据结构、算法在有限元程序设计中的作用。
技能目标：1. 能够运用所学知识解决简单的工程问题，通过有限元方法进行力学分析。
2. 具备独立操作有限元软件的能力，完成模型建立、计算及结果分析的完整流程。
3. 能够根据实际问题需求，编写简单的有限元程序，提高编程实践能力。
情感态度价值观目标：1. 培养学生对工程问题的探究精神，激发学生主动学习的兴趣。
2. 增强学生的团队合作意识，培养沟通协调能力，提高解决实际问题的能力。
3. 使学生认识到有限元技术在工程领域的重要价值，树立正确的科技观。
课程性质：本课程为专业选修课，旨在让学生掌握有限元程序设计的基本方法，提高解决工程问题的能力。
学生特点：学生具备一定的编程基础，对有限元分析有初步了解，但实践能力较弱。
教学要求：注重理论与实践相结合，强调学生动手实践，培养解决实际问题的能力。
通过本课程的学习，使学生能够将所学知识应用于工程实践，提高综合素养。
二、教学内容1. 有限元分析基本原理：包括有限元离散化方法、变分原理、刚度矩阵和质量矩阵的构建等。
教材章节：第一章 有限元分析概述，第二章 有限元离散化方法。
2. 有限元软件操作：介绍主流有限元软件的功能、操作流程，以ANSYS为例进行实践教学。
教材章节：第三章 有限元软件及其应用。
3. 有限元程序设计：讲解有限元程序设计的基本步骤、数据结构、算法实现等。
教材章节：第四章 有限元程序设计基础，第五章 数据结构及算法。
4. 实践案例：选取具有代表性的工程问题，指导学生运用有限元软件和编程技能解决问题。
教材章节：第六章 实践案例。
5. 课程项目：分组进行项目实践，要求学生完成项目报告和成果展示。
教材章节：第七章 课程项目与实践。
有限元单元法程序设计有限元单元法程序设计是一种广泛应用于工程领域的数值计算方法，它能够模拟复杂结构的受力情况并计算出相应的应力、变形等物理量。
本文将从有限元单元法的基本原理、程序设计流程、关键步骤等方面入手，为您详细介绍有限元单元法程序设计的相关内容。
一、有限元单元法基本原理有限元单元法是一种工程结构分析的数值计算方法，它基于弹性力学原理，将结构划分为有限个小单元（有限元）进行离散化处理，通过对各个单元的力学行为进行分析来描述整个结构的受力情况。
有限元单元法的基本原理可以总结为以下几个步骤：1. 将结构离散化为有限个小单元，每个单元内的应力、变形等物理量满足弹性力学理论。
2. 建立每个单元的位移与节点力之间的关系，通常采用单元刚度矩阵来描述。
3. 根据整个结构的连接条件和边界条件，组装各个单元的刚度矩阵，形成整个结构的刚度矩阵。
4. 应用外载荷和边界条件，求解整个结构的位移场，并由此计算出应力、变形等物理量。
二、有限元单元法程序设计流程有限元单元法程序设计通常包括以下几个关键步骤，我们将逐步介绍其设计流程：1. 确定结构的几何形状和材料性质，将结构进行离散化处理，确定有限元的类型和数量。
2. 建立单元刚度矩阵的表达式，通常采用弹性力学理论和数值积分方法来进行推导和计算。
3. 将各个单元的刚度矩阵组装成整个结构的刚度矩阵，考虑节点之间的连接关系以及边界条件的处理。
4. 应用外载荷和边界条件，求解整个结构的位移场，并计算出节点处的应力、变形等物理量。
5. 对程序进行稳定性和准确性的验证，包括收敛性分析、误差估计等。
6. 编写相应的有限元单元法程序，实现结构的建模、求解和结果输出等功能。
三、有限元单元法程序设计的关键步骤在有限元单元法程序设计中，有几个关键的步骤需要特别重视：1. 单元选择和刚度矩阵的建立：选择适合结构特点的有限元类型，建立单元的刚度矩阵表达式，考虑单元的形函数、应变-位移关系等。
2. 结构刚度矩阵的组装：将各个单元的刚度矩阵通过节点的连接关系组装成整个结构的刚度矩阵，考虑节点自由度的排序和边界条件的处理。
附录有限元教学程序及使用说明§A-1平面三角形３结点有限元程序1、程序名：FEM3.FOR，FEM3.EXE2、程序功能本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。
主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结点的支座反力。
程序采用Visual Fortran 编制而成，输入数据全部采用自由格式。
3、程序流程及框图图A-1 程序流程图277278图A-2 程序框图其中，各子程序的功能如下：INPUT ——输入结点坐标、单元信息和材料参数； MR ——形成结点自由度序号矩阵；FORMMA ——形成指标矩阵MA （N ）并调用其他功能子程序，相当于主控程序； DIV ——取出单元的3个结点号码和该单元的材料号并计算单元的b i ，c i 等； MGK ——形成整体劲度矩阵并按一维存放在SK （NH ）中； LOAD ——形成整体结点荷载列阵F ； OUTPUT ——输出结点位移或结点荷载；TREAT ——由于有非零已知位移，对K 和F 进行处理； DECOMP ——整体劲度矩阵的分解运算； FOBA ——前代、回代求出未知结点位移 ； ERFAC ——计算约束结点的支座反力； KRS ——计算单元劲度矩阵中的子块K rs 。
4、程序使用说明当程序开始运行时，按屏幕提示，键入数据文件的名字。
在运行程序之前，必须根据程序中输入要求建立一个存放原始数据的文件，这个文件的名字由少于8个字符或数字组成。
数据文件包括如下内容： ⑴ 总控信息，共一条，9个数据NP ，NE ，NM ，NR ，NI ，NL ，NG ，ND ，NC NP ——结点总数； NE ——单元总数；NM ——材料类型总数；279NR ——约束结点总数；NI ——问题类型标识，0为平面应力问题，1为平面应变问题； NL ——受荷载作用的结点的数目；NG ——考虑自重作用为1，不计自重为0； ND ——非零已知位移结点的数目；NC ——要计算支座约束反力的结点数目。
第十章有限元程序设计有限元方法作为一门系统的技术，仅学会了它的基本理论是远远不够的，只有形成完整的计算程序，问题才最终得到了解决。
完成这样的有限元程序设计是一项工作量很大的工程。
本章就是要结合简单的有限元教学程序FEMED，简要介绍有限元程序设计技术。
FEMED是专为有限元程序设计教学编制的程序，它不包含复杂的前后处理功能，可进行平面问题及平面桁架的线弹性静力分析，在程序结构上与大型程序类似，具有计算单元的任意扩充功能，在方程的组集和求解上也采用了较为流行的变带宽存储方式。
有限元程序大致可分为两类，第一类是专用程序，主要用于研究或教学，一般这类程序规模较小，前后处理功能较弱。
用于研究的程序能够解一些特殊的问题，满足研究工作的需要。
而教学程序则是为了学生了解有限元的主要结构和设计方法设计的，程序比较简单，FEMED就属于这类程序。
第二类是大型通用程序，是大型结构分析的得力工具，目前国际上流行的大约有2000多种。
常用的有NASTRAN、MARC、ANSYS、ADINA和ABAQUS等。
这类程序一般前后处理功能比较强，有友好的界面，能进行大型计算，但往往无法完成具有特殊要求的计算。
通过本章的学习，使读者初步掌握有限元编程的基本方法，具有开发特殊功能的专用程序或为通用程序开发具有特殊功能的计算模块的能力。
§10.1有限元程序的基本结构有限元程序一般包括三项基本内容：前处理、结构分析和后处理。
早期有限元分析软件的研究重点在于推导新的高效率求解方法和高精度的单元，随着数值分析方法的逐步完善，尤其是计算机内存和运算速度的飞速发展，整个计算系统用于求解运算的时间越来越少，加之求解问题的日益大型化和复杂化，使得数据准备和运算结果的表现问题日益突出。
因此目前几乎所有的商业化有限元程序系统都有功能很强的前后处理模块，这直接关系到分析软件的可推广性。
它是商用有限元软件不可或缺的部分，但它不是有限元的中心部分，在本书中不作详细介绍。
10.1.1 前处理在结构分析之前，必须完成的工作就是前处理的任务，它主要包括以下内容：结构造型、选取单元类型，、网格剖分，从而形成节点数、节点编码、节点坐标、单元数、单元节点编码等。
另外根据所解题目的不同，还需读入不同的材料参数、边界约束条件及载荷工况等。
有限元前处理程序的功能就是为用户提供一种工具，使其能尽可能方便地完成上述工作。
而这一问题解决的好坏直接关系到程序使用者需付出的工作量，而网格剖分质量还直接影响计算精度。
很多程序都建立了对用户非常友好的图形用户界面，使用户能以可视图形方式直观快速地进行结构造型和网格自动划分，生成有限元分析所需数据。
在上述工作中工作量最大的就是结构造型和网格剖分部分，目前越来越倾向于由专用软件完成上述工作，所以当今所有的商业化有限元软件开发商都开发了和著名的CAD软件的接口，从而极大地提高了设计水平和效率。
今天，工程师可以在集成的CAD和FEA软件环境中快捷地解决一个在以前无法应付的复杂工程分析问题。
在FEMED程序中，只包含了很简单的节点、网格及给定位移的自动生成。
所有前处理功能是由INPUT及相关子程序完成的。
10.1.2 结构分析结构分析部分是有限元程序的主体部分，主要需计算插值函数矩阵[N]，插值函数的导数矩阵[B]，进而进行数值积分得到单元刚度矩阵，组集总刚度矩阵及总载荷向量，解线性代数方程组，得到节点位移。
对于动力问题则需计算质量矩阵，求解特征值问题。
而对于动力响应问题及非线性问题，则需进行多增量步的计算。
在得到节点位移的基础上还需进行应力的计算。
在本章中重点介绍的就是结构分析中的静力分析部分，在FEMED程序中主要包含的也就是这部分。
通过这部分的学习，希望读者能掌握有限元程序的基本编程技巧。
10.1.3 后处理在大型分析软件中，程序的后处理功能也是非常重要的，在有限元结构分析完成后，通过强有力的后处理图形功能，可以给出各应力分量、等效应力等的分布云图，结构的变形图，振型图和实图10.1时动画等，使得使用者对各物理量在整个结构中的变化情况有一个全面的认识，也可检查网格剖分和所加载荷及约束是否有误。
否则面对输出成千上万个计算结果数据，往往使人如坠雾中，不得要领，需花费大量的时间进行数据的整理和解释。
由于是简单的教学程序，在FEMED程序中不包含后处理功能，仅仅是将计算结果输出。
图10.1给出了一般求解静力学问题的程序框图，FEMED程序的结构基本如图所示，但没有包含带宽优化部分。
§10.2 数组的半动态分配在有限元程序中需开辟许多不同功能的数组，而这些数组的大小则与所解题目的大小及类型有关，如节点坐标数组、节点位移数组等与节点数相关，而单元节点编号数组等则与单元数相关。
此外，上述数组的大小还和自由度维数等相关。
对于一个有限元程序而言，应具有求解各类问题及各种规模问题的机动性，因此求解不同的问题，应定义不同大小的数组。
FORTRAN 语言规定数组的大小必须预先给定，在子程序中可以定义大小可变化的可调数组，但其内存分配必须由该子程序的上一级模块确定。
最简单的办法就是每个数组都定义得足够大，使得要求解的问题都能被满足，但这显然会带来内存的巨大浪费，因此，一般采用数组半动态分配的方法来解决上述问题，FEMED 也采用了这一方法。
在主控程序PCONTR 中开辟了两个大数组A(*)和M(*)，分别存储实型量和整型量，将各数组中的数顺序存放在大数组中，在读入控制信息NJ （节点数）、NE （单元数）、NUMMAT （材料数）、NDF （节点自由度数）、NDM （空间维数）、NEN （单元节点数）等之后，就可进行存储单元的分配。
在图10.2中以大数组A 的存储为例，先存放节点坐标，占用NJ ×NDM 个存贮单元，再存放材料参数，占用NUMMA T ×20个存贮单元……图10.2在上述过程中，最重要的是计算出各数组存放的起始位置又称作地址，如NXC 为XC 数组在A 中的地址，ND 为D 数组在A 中的地址……，这项工作由POINT 子程序完成。
需要分别计算出大数组A 和M 中所存数组的所有地址，最后计算出A 和M 数组需要的总长度，与预设的数组长度作比较，若没有超界计算继续进行，若超界则计算停止，显示错误信息。
需要指出的是，在上述数组的半动态分配时，将一系列存储容量不定的一维、二维数组均用地址对应关系分配到事先容量固定的大数组A 和M 中去，对于大数组A 和M 来说存储量是不变的，但各数组在大数组中的存储分配却是动态的。
我们把这种关系称为数组的半动态分配，这一过程是通过哑实结合完成的。
这些数组称为可调数组，即其数组长度是可调的。
在子程序的数组说明中，一维数组写成下列形式是等价的，对数组长度没有限制：F(*)，F(1)，F(NJ)，对于二维数组，则只可一个方向自由变动。
由于FORTRAN 规定二维数组先变行后变列，因此行的长度必须是确定的。
，在子程序的数组说明中，下列形式是等价的：XC(2,*)，XC(2,1)，XC(2,NJ)，XC(NDM,*)，(NDM=2)注意，只有对上一模块中已定义过的数组才可作如上的数组说明，否则，则必须准确确定数组的长度。
§10.3 构造单元刚度矩阵及总刚度矩阵的组集构造单元刚度矩阵是由单元模块完成的，不同的单元调用不同的模块，在FEMED 中提供了平面桁架单元（TRUSS ）和平面单元（PLANE ）两个单元计算模块。
读者可根据需要随意开发单元类型，通过ELEMLIB 子程序接入程序中。
针对不同的单元，在不同的单元模块中构造不同的[B]、[D]矩阵，单元刚度矩阵可表示为：VdV =⎰e e T K B DB（10.1）对于常应变单元，上述被积函数为常数，直接用被积函数值乘积分域的大小即可。
而对于等参元等非常应变单元，则需对上式做GAUSS 积分，最终得到一个二维的单元刚度矩阵。
将单元刚度矩阵正确组集为总刚度矩阵，这一工作是由ADDSTF 子程序完成的。
下面将举例说明总刚度矩阵的组集方法。
例1：平面桁架结构如图示，共有5个节点，6个单元。
每个单元形成一个4×4的单元刚度矩阵，总刚度阵为10×10。
图10.3以1号单元为例，其左节点为2，右节点为4，其4×4的单元刚度矩阵的各个参数应迭加在(10.2)式总刚度矩阵中的图示位置上。
需要指出的是，总刚度矩阵的一个位置上往往有多个单元刚度矩阵做贡献，如与3号节点对应的(5,5)，(6,5)，(6,6)位置上就有2、3、4和6号单元都作了贡献。
全部单元组集完毕后，总刚度矩阵如式（10.2），形成一对称带状阵，即仅对角元附近的元素不为零，而远离对角线的元素皆为零。
解题越大带状性质越明显。
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯=00000000444342410033323100222111][称对K（10.2）若储存的是总刚度矩阵的上三角元素，则带宽为每行从对角元，到最后一个非零元素的个数；若储存的是总刚度矩阵的下三角元素，则从第一个非零元素到对角元元素的个数称为这行的带宽。
每行的带宽各不相同，且与节点编号的方式有关，在大型商业软件中，一般都包含带宽优化模块，就是为了减少总体带宽。
§10.4 总刚度矩阵的变带宽存储总刚度阵占用的空间是有限元占用空间的绝大部分，一个软件的解题规模往往是受此限制。
所以，一直以来人们在这方面想了很多的办法。
很显然，将总刚度阵作为一个二维数组，完全储存起来是最直接和最自然的方法，但这也是占用空间最大、最不合理的方法，一般很少采用。
考虑到总刚度阵具有带状稀疏对称的特征，使用半带宽存储无疑是一个最自然的改进，它利用一个二维数组，只存储每行从对角元到这一行的最后一个非零元，如图10.4所示⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯⨯0 000000000000000000000000000000000000称对图10.4由10.4的示意图知，原10×10的矩阵现在可用10×4的矩阵存储，存储空间减少了很多，但若其中有个别行的带宽特别大，则这种方法依然有存贮空间的大量浪费。
为解决上述问题，一般多采用变带宽存储又称一维压缩存储的方法，这样刚度阵中大量的零既不需要存储也不参与计算，大大减少了存储空间和计算量，是目前使用较为广泛的方法。
在有限元的发展历史中波前法也是一种不可不提的计算方法，当计算机的内存非常有限时，即使采用变带宽存储也无法求解较大规模的问题，波前法曾经起到了不可忽视的作用。
