《有限单元法》编程作业

有限单元法考试题及答案一、单项选择题（每题2分，共10分）1. 有限元法中，单元刚度矩阵的计算是基于（）。

A. 位移法B. 势能原理C. 能量守恒定律D. 牛顿第二定律答案：B2. 在有限元分析中，以下哪项不是网格划分时需要考虑的因素？（）A. 网格数量B. 网格形状C. 材料属性D. 边界条件答案：C3. 有限元分析中，以下哪项不是结构分析的基本步骤？（）A. 离散化B. 求解C. 后处理D. 优化设计答案：D4. 在有限元分析中，以下哪种类型的单元不适用于平面应力问题？（）A. 三角形单元B. 四边形单元C. 六面体单元D. 楔形单元答案：C5. 有限元分析中，以下哪种边界条件不属于几何边界条件？（）A. 固定支座B. 压力C. 温度D. 位移答案：C二、多项选择题（每题3分，共15分）6. 有限元法中，以下哪些因素会影响单元的精度？（）A. 单元形状B. 单元数量C. 材料属性D. 网格划分答案：ABD7. 在有限元分析中，以下哪些是常见的数值积分方法？（）A. 一阶积分B. 二阶积分C. 高斯积分D. 牛顿-莱布尼茨积分答案：ABC8. 有限元分析中，以下哪些是常见的单元类型？（）A. 线性单元B. 二次单元C. 三次单元D. 非线性单元答案：ABCD9. 在有限元分析中，以下哪些是常见的后处理技术？（）A. 应力云图B. 位移云图C. 模态分析D. 热分析答案：ABC10. 有限元分析中，以下哪些是常见的非线性问题？（）A. 几何非线性B. 材料非线性C. 接触非线性D. 热应力问题答案：ABCD三、填空题（每题2分，共20分）11. 有限元法中，单元刚度矩阵的计算通常基于___________原理。

答案：势能12. 在有限元分析中，网格划分的目的是将连续的___________离散化为有限数量的单元。

答案：域13. 有限元分析中，___________是将实际问题转化为数学问题的关键步骤。

-----好资料学习有限单元法中“离散”的含义是什么？有限单元法是如何将具有无限自由度的连续介1.1质问题转变成有限自由度问题的？位移有限元法的标准化程式是怎样的？）离散的含义即将结构离散化，即用假想的线或面将连续体分割成数目有限的单元，并1（数的节在其上设定有限个节点；用这些单元组成的单元集合体代替原来的连续体，而场函点值将成为问题的基本未知量。

）给每个单元选择合适的位移函数或称位移模式来近似地表示单元内位移分布规律，即2（无限自通过插值以单元节点位移表示单元内任意点的位移。

因节点位移个数是有限的，故由度问题被转变成了有限自由度问题。

）有限元法的标准化程式：结构或区域离散，单元分析，整体分析，数值求解。

（3 ?单元刚度矩阵和整体刚度矩阵各有哪些性质？各自的物理意义是什么?两者有何区别1.3整体刚度矩阵的性单元刚度矩阵的性质：对称性、奇异性（单元刚度矩阵的行列式为零）。

个自j Kij 即单元节点位移向量中第稀疏性。

单元 Kij 物理意义质：对称性、奇异性、整体刚度 j 个自由度方向引起的节点力。

由度发生单位位移而其他位移分量为零时，在第中每一列元素的物理意义是：要迫使结构的某节点位移自由度发生单位位移，而其 K 矩阵他节点位移都保持为零的变形状态，在所有个节点上需要施加的节点荷载。

什么叫应变能？什么叫外力势能？试叙述势能变分原理和最小势能原理，并回答下述2.2问题：势能变分原理代表什么控制方程和边界条件？其中附加了哪些条件？，外力所做的功将以变形能的形式储存εσ和应变(1)在外力作用下，物体内部将产生应力起来，这种能量称为应变能。

(2)外力势能就是外力功的负值。

势能变分原理可叙述如下：在所有满足边界条件的协调位移中，那些满足静力平衡条件(3) 的位移使物体势能泛函取驻值，即势能的变分为零V=0 +δp=δ Uεδ∏此即变分方程。

对于线性弹性体，势能取最小值，即02V≥ε+δδ2∏P=δ2U 此时的势能变分原理就是著名的最小势能原理。

第三章 平面问题的有限元法作业1. 图示一个等腰三角形单元及其节点编码情况，设μ=0，单元厚度为t 。

求 1）形函数矩阵[]N ；2）应变矩阵[]B ；3）应力矩阵[]S 。

4第1题图 第2题图2. 如题图所示，结构为边长等于a 的正方形，已知其节点位移分别为：11(,)u v 、22(,)u v 、33(,)u v 、44(,)u v 。

试求A 、B 、C 三点的位移。

其中A 为正方形形心，B 为三角形形心。

3.直角边边长为l 的三角形单元，如题图所示。

试计算单元等效节点载荷列阵（单元厚度为t ，不计自重）。

第3题图 第4题图4. 如题图所示，各单元均为直角边边长等于l 的直角三角形。

试计算（1）单元等效节点载荷列阵；（2）整体等效节点载荷列阵。

已知单元厚度为t ，不计自重。

5．下列3个有限元模型网格，哪种节点编号更合理？为什么？934679121134612142(a) (b) (c)第5题图6.将图示结构画出有限元模型；标出单元号和节点号；给出位移边界条件；并计算半带宽（结构厚度为t ）。

2a(a) (b) 无限长圆筒 (c) 第6题图7. 结构如图所示，已知结构材料常数E 和 ，单元厚度为t 。

利用结构的对称性，采用一个单元，分别计算节点位移和单元应力。

第7题图答案：1. 1）形函数i x N a =， j y N a = ， 1m x y N a a=-- 2）应变矩阵[]1000101000101011011B a -⎡⎤⎢⎥=-⎢⎥--⎢⎥⎣⎦3）应力矩阵[]10001010001011111002222S a ⎡⎤⎢⎥-⎢⎥=-⎢⎥⎢⎥--⎢⎥⎣⎦2. A 点的位移为()2312A u u u =+ ， ()2312A v v v =+ B 点的位移为()24313B u u u u =++ ， ()24313B v v v v =++ C 点的位移为()1223C a u u u =+ ， ()C 1223av v v =+ 3. 单元等效节点载荷列阵为{}111100003663Tei ji jR q q q q ⎡⎤=++⎢⎥⎣⎦4. （2）整体等效节点载荷向量为{}111100006322TR qlt P qlt P P qlt qlt ⎡⎤=-⎢⎥⎣⎦7. （1） 减缩后的整体刚度方程22122122222221110222021102(1)22102x x b b ab R b ab b P v Etab a bab ab R v b a μμμμμμμμμ---⎡⎤--⎢⎥⎧⎫⎧⎫⎢⎥⎪⎪--⎪⎪⎢⎥⎪⎪-⎪⎪⎢⎥=⎨⎬⎨⎬---++⎢⎥⎪⎪⎪⎪⎢⎥⎪⎪⎪⎪⎩⎭-⎢⎥⎩⎭+⎢⎥⎣⎦ 节点位移22(1)Pb v Eatμ+=- ， 2212212b a v v bμ-+=单元应力为{}()2122201012bv E bv bv ab av μσμμ⎛⎫⎧⎫ ⎪⎪⎪-⎧⎫ ⎪⎪⎪⎪⎪=+-⎨⎬⎨⎬ ⎪-⎪⎪⎪⎪ ⎪-⎩⎭⎪⎪-⎪⎩⎭⎝⎭。

有限单元法习题答案有限单元法（Finite Element Method，简称FEM）是一种数值计算方法，用于求解工程和物理问题的数学模型。

它将复杂的连续体分割成许多简单的有限单元，通过对每个单元进行离散化，近似求解整个问题。

在实际应用中，有限单元法广泛应用于结构力学、流体力学、电磁学等领域。

在学习过程中，我们常常会遇到一些习题，下面将给出一些有限单元法习题的答案，希望对大家的学习有所帮助。

1. 有限单元法的基本原理是什么？答：有限单元法的基本原理是将连续体分割成有限个简单的单元，通过对每个单元进行离散化，建立局部方程，再通过组装得到整体方程。

通过求解整体方程，得到问题的近似解。

2. 如何选择合适的有限单元？答：选择合适的有限单元是保证计算精度的关键。

一般来说，有限单元的选择应该满足以下几个条件：简单性、合理性、适应性和可靠性。

常见的有限单元包括一维线元、二维三角形单元、二维四边形单元等。

3. 有限单元法的求解步骤是什么？答：有限单元法的求解步骤一般包括以下几个步骤：建立有限元模型、确定边界条件、选择适当的有限单元、建立单元刚度矩阵和载荷向量、组装单元刚度矩阵和载荷向量、施加边界条件、求解代数方程组、计算节点位移和应力、分析结果的准确性。

4. 有限单元法的优缺点是什么？答：有限单元法的优点包括：适用范围广、计算精度高、计算效率高、易于处理复杂边界条件等。

缺点包括：模型的精度受到有限单元的选择和网格划分的影响、计算结果的可信度需要通过验证、对计算机硬件要求较高等。

5. 有限单元法在结构力学中的应用有哪些？答：有限单元法在结构力学中的应用非常广泛，包括静力分析、动力分析、热力分析等。

例如，在静力分析中，可以通过有限单元法求解结构的受力状态；在动力分析中，可以通过有限单元法求解结构的振动特性；在热力分析中，可以通过有限单元法求解结构的温度分布等。

6. 有限单元法在流体力学中的应用有哪些？答：有限单元法在流体力学中的应用也非常广泛，包括流体流动、传热、质量传递等。

有限元方法编程【实用版1篇】目录（篇1）1.有限元方法概述2.有限元方法的编程步骤3.有限元方法的应用实例4.总结正文（篇1）一、有限元方法概述有限元方法是一种数值分析方法，广泛应用于固体力学、流体力学、热传导等领域。

它的基本思想是将待求解的连续体划分为有限个小的、简单的子区域，即单元，然后用有限个简单的方程组来代替原来的连续方程，通过求解这些方程组得到近似解。

这种方法既能降低问题的复杂度，又能保证解的精度，因此在工程界有着广泛的应用。

二、有限元方法的编程步骤1.几何建模：根据实际问题，创建待求解的几何模型。

这通常包括划分单元、计算节点坐标等步骤。

2.选择单元类型：根据问题类型和求解需求，选择合适的单元类型，如有限元、无限元、矩形单元、六面体单元等。

3.编写有限元方程：根据单元类型和几何模型，编写有限元方程。

这包括计算单元的刚度矩阵、质量矩阵、载荷矩阵等。

4.组装总方程：将所有单元的有限元方程组装成总方程，通常是一个大型的线性或非线性方程组。

5.求解方程组：使用数值方法（如有限元法、直接解法、迭代法等）求解总方程组，得到近似解。

6.后处理：对求解结果进行分析和处理，如计算应力、应变、位移等。

三、有限元方法的应用实例以一个简单的二维拉伸问题为例，假设有一个长方形板，在左右两端施加均匀拉力，求解板上各个点的应力和应变。

1.几何建模：将长方形板划分为矩形单元，计算节点坐标。

2.选择单元类型：此处采用矩形单元。

3.编写有限元方程：计算单元的刚度矩阵、质量矩阵、载荷矩阵，组装总方程。

4.求解方程组：使用有限元法求解总方程组，得到应力和应变。

5.后处理：分析应力和应变分布，验证解的正确性。

四、总结有限元方法作为一种数值分析方法，通过将连续体划分为有限个小的、简单的子区域，然后用有限个方程组来代替原来的连续方程，降低了问题的复杂度，同时保证了解的精度。

在实际应用中，有限元方法需要经历几何建模、单元选择、编写有限元方程、组装总方程、求解方程组和后处理等步骤。

有限单元法与程序编程实践与课程报告题目1.在三角形单元程序的基础上，增加矩形单元，完成混合单元（三角形和矩形）的有限元分析程序，并在此基础上添加热应力，并给出算例。

（难度系数0.8）2.在3结点三角形单元程序基础上，完成6结点三角形单元的有限元分析程序，并给出算例。

（难度系数0.85）3.写出由3结点三角形单元节点信息自动生成6结点三角形单元的单元节点信息的算法及fortran子程序，用算例验证程序的正确性。

4.在3结点三角形单元程序基础上，完成三角形单元的子结构有限元分析程序，并给出算例。

（难度系数1.0）5.在3结点三角形单元程序上，完成6结点三角形单元的子结构有限元分析程序，并给出算例。

（难度系数1.0）6.在3结点三角形单元程序基础上，完成空间四面体单元的有限元分析程序，并给出算例。

（难度系数0.8）7.在3结点三角形单元程序基础上，完成空间轴对称三角形截面环单元的有限元分析程序，（主要包括应变矩阵子程序、应力矩阵子程序、等效节点荷载子程序等）并给出算例。

（难度系数0.9）8.在3结点三角形单元程序基础上，完成平面四节点和八节点四边形等参数单元的有限元分析程序，并给出算例。

（难度系数1.0）9.在3结点三角形单元程序的基础上，完成空间轴对称四节点和八节点等参数单元的有限元分析程序，并给出算例。

（难度系数1.2）10.在现有程序的基础上，完成空间轴对称问题三节点截面环单元和四节点等参数单元组成混合单元的有限元分析程序，并给出算例。

（难度系数1.1）11.在现有程序的基础上，完成空间轴对称四节点和八节点等参数单元的有限元分析程序，并给出算例。

（难度系数1.1）12.在现有程序的基础上，完成空间二十节点等参数单元的有限元分析程序，并给出算例。

（难度系数1.2）（六面体）13.在3结点三角形单元程序基础上，薄板弯曲矩形单元的有限元分析程序，并给出算例。

（难度系数1.1）14.在3结点三角形单元程序基础上，薄板弯曲三角形单元的有限元分析程序，并给出算例。

填空题：1、利用有限单元法求解弹性力学问题时，简单来说包含结构离散化、单元分析、整体分析三个主要步骤。

2、有限单元法首先将连续体变换成为离散化结构，然后再用结构力学位移法进行求解。

其具体步骤分为单元分析和整体分析两部分。

3、每个单元的位移一般总是包含着两部分：一部分是由本单元的形变引起的，另一部分是由于其他单元发生了形变而连带引起的。

4、每个单元的应变一般总是包含着两部分：一部分是与该单元中各点的位置坐标有关的，是各点不相同的，即所谓变量应变；另一部分是与位置坐标无关的，是各点相同的，即所谓常量应变。

5、为了能从有限单元法得出正确的解答，位移模式必须能反映单元的刚体位移和常量应变，还应当尽可能反映相邻单元的位移连续性。

6、为了使得单元内部的位移保持连续，必须把位移模式取为坐标的单值连续函数，为了使得相邻单元的位移保持连续，就不仅要使它们在公共结点处具有相同的位移时，也能在整个公共边界上具有相同的位移。

7、在有限单元法中，单元的形函数N i在i 结点N i= 1 ；在其他结点N i= 0 及∑N i= 1 。

8、为了提高有限单元法分析的精度，一般可以采用两种方法：一是将单元的尺寸减小，以便较好地反映位移和应力变化情况；二是采用包含更高次项的位移模式，使位移和应力的精度提高。

9、在有限单元法中，结点力是指结点对单元的作用力。

（√）10、在平面三结点三角形单元的公共边界上应变和应力均有突变。

（√）11、形函数N i(xi,yi)= ＿＿(i=j)N i(xi,yi)= ＿＿(i≠j)简答题：1、有限元分析的基本思路答：首先，将物体或求解域离散为有限个互不重叠仅通过节点互相连接的子域（即单元），原始边界条件也被转化为节点上的边界条件，此过程称为离散化。

其次，在单元内，选择简单近似函数来分片逼近未知的求解函数，即分片近似。

具体做法是在单元上选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，这是有限元法的创意和精华所在。

有限单元法与程序Ansys 上机实践题目1. 图1示为一带圆孔的单位厚度的正方形平板，在x 方向作用均布压力0.25Mpa ，试用矩形单元和六节点三角形单元对平板进行有限元分析，绘出变形及应力云图，并与弹性力学理论解对比，说明有限元数值解的特性。

图12. 试生成图2所示双孔板拉伸问题恰当的有限元网格（L =1000mm ），计算孔边最大应力σmax ，绘出变形及应力云图。

并通过网格细分说明结果的精度。

图23. 两端简支、狭长矩形截面直梁，顶面受到集度为q 的均布荷载作用，求中间截面的应力分布及梁中性线上的剪应力分布曲线和挠度曲线。

梁长为80cm ，高度10cm ，厚度1cm ，材料E=2×10^5 MPa ，γ=0.3，q=10MPa ，计算要求：（1）绘出梁中间截面上的应力分布图，并标明最大值；（2）绘出梁中性线的剪应力的分布曲线和挠度曲线，并标明最大值；（3）比较不同单元形式和网格密度对计算结果的影响，分析位移解的收敛性和收敛速度及位移解的下限性；（4）与材料力学及弹性力学解答比较，并分析造成差别的原因。

4. 如图4示带方孔（边长为80mm ）的悬臂梁，其上受部分均布载荷（p=10Kn/m ）作用，试采用一种平面单元，对图示两种结构进行有限元分析，并就方孔的布置进行分析比较，如将方孔设计为圆孔，结果有何变化？（板厚为1mm ，材料为钢）σ0σ图45. 如图5所示八角形管道，取图中八角形管道的22.5°区域建立有限元模型，求出这部分变形后的形状，并给出平面内的最大剪应力分布。

图56. 如图6所示的扭矩臂是一个汽车零件，它被固定在左端螺栓孔处。

求最大的Von Mises应力的位置和大小。

图67. 如图7所示一铜制槽环，在受到两个大小都为P的力作用下，宽度为3mm的右缺口恰好闭合，求荷载P的大小和变形后的形状。

图78.对图8示具有钢衬套的混凝土压力容器进行应力分析.取混凝土的E=30GPa，μ=0.15. 钢衬套的E=205GPa，μ=0.25.钢衬套的厚度为50mm，压力P=700KPa。

西北工业大学22春“公共课”《有限元及程序设计》作业考核题库高频考点版（参考答案）一.综合考核(共50题)1.离散化过程应遵循()的原则。

A.位移近似B.几何近似C.物理近似D.形状近似参考答案：BC2.三角形单元三个节点编码应按()编排。

A.顺时针B.逆时针C.线性D.非线性参考答案：B3.剪力墙墙体受自重作用属于平面应变问题，天然气管道属于平面应力问题。

()A.错误B.正确参考答案：A4.三结点三角形单元是常应力单元，但不是常应变单元。

()A.错误B.正确参考答案：A总体刚度矩阵具有()性质。

A.对称性B.稀疏性C.带状分布性D.奇异性参考答案：ABCD6.下列不属于提高单元精度的方法是()。

A.增加单元结点数目B.在单元内增设结点C.减少单元结点数目D.设等参元参考答案：C7.弹性力学问题的求解方法有()。

A.按应变求解B.按应力求解C.按体力求解D.按位移求解参考答案：BD8.将单元内所受到的荷载(包括体力、面力和集中力)通过静力等效的原则移置到结点上的过程称等效结点荷载。

()A.错误B.正确参考答案：B9.弹性力学的基本假设有()。

A.假设物体是连续的D.假设物体内无初应力E.假设物体的变形是很小的参考答案：ABCDE10.解平面应力和平面应变问题采用的应力矩阵相同。

()A.错误B.正确参考答案：B11.有限元的单元局部编码，应按逆时针规律编排。

()A.错误B.正确参考答案：B12.薄板小挠度弯曲理论的基本假定是()。

A.直法线假定B.法向位移假定C.中面位移假定D.板内无挤压假定参考答案：ACD13.应力函数线性项可解决矩形板拉伸问题。

()A.错误B.正确参考答案：AA.错误B.正确参考答案：B15.应变分量中正应变以伸长为正，剪应变以夹角小为正。

()A.错误B.正确参考答案：B16.矩形薄板单元是完全协调单元。

()A.错误B.正确参考答案：A17.下列属于高精度空间单元的有()。

A.10结点30自由度四面体单元B.20结点60自由度六面体单元C.6结点三角形单元D.4结点48自由度四面体单元参考答案：ABD18.广义结点力是垂直于x轴和y轴的弯矩和单位长度上的扭矩。

有限单元法及程序设计有限单元法（Finite Element Method，FEM）是一种用于数值分析和计算的方法，广泛应用于工程和科学领域。

它通过将连续问题离散化成有限个小单元，并在每个小单元上建立数学模型来近似求解问题。

本文将介绍有限单元法的基本原理、步骤以及程序设计方面的注意事项。

一、有限单元法基本原理有限单元法的基本原理是将连续的物理区域划分为有限个离散的小单元，每个小单元内的场量近似表示为一些插值函数的线性组合。

通过对这些小单元进行逐个求解，最终得到整个问题的近似解。

有限单元法的核心思想是利用局部性原则，将整个问题分解成多个小问题。

每个小问题只涉及到相邻的单元，在确定了边界条件和材料特性后，可以进行独立的求解。

最后通过组合各个小问题的解，得到整个问题的解。

二、有限单元法步骤有限单元法的求解过程主要包括几个基本步骤，具体如下：1. 离散化：将连续的物理区域划分为有限的小单元。

常用的小单元形状包括三角形、四边形、六边形等。

2. 建立数学模型：在每个小单元上建立数学模型，通常使用插值函数来近似表示物理量。

插值函数的选择对求解结果的准确性和效率有重要影响。

3. 形成总体方程：根据物理规律和边界条件，利用适当的数学方法推导出总体方程。

常见的总体方程包括稳定性方程、运动方程等。

4. 矩阵装配：将每个小单元的局部方程装配成整个系统的总体方程。

这一步骤常常需要对单元进行编号和排序，以便正确地装配矩阵。

5. 边界条件处理：根据实际问题的边界条件，对总体方程进行修正。

边界条件的处理通常包括施加约束和设定边界值。

6. 求解方程：通过数值方法，如有限差分法或有限元法，求解总体方程。

常用的求解方法包括直接法和迭代法。

7. 后处理：对求解结果进行计算和分析，以获得实际问题的有用信息。

后处理包括输出位移、应力、应变等字段，以及进行可视化展示。

三、程序设计注意事项在进行有限单元法的程序设计时，需要充分考虑以下几个方面的注意事项：1. 算法选择：根据问题的特点和求解需求，选择合适的有限单元类型、插值函数和数值解法。

航空航天结构中的有限元方法编程大作业姓名熊蕾班级110514学号11051121目录一、使用程序介绍 (1)二、数据文件生成程序 (1)三、计算程序 (2)四、例题求解 (4)例1: (4)例2： (5)例3： (6)五、程序不足之处与改进方法 (7)1、数据文件的可读性 (7)2、数据长度混乱 (7)3、输入量很难从图形得出 (8)4、算法的性质 (8)六、参考文献 (8)一、使用程序介绍由于有限元刚度方程求解过程是一个矩阵运算的过程，使用matlab编写程序会比C语言更加方便。

MATLAB是近年来得到迅速发展的数学软件，它将高性能的数值计算和可视化集成在一起，并提供了大量内置函数，被广泛地应用于科学计算、控制系统、信息处理等领域的分析、仿真和实际工作。

二、数据文件生成程序为了方便计算程序读取数据，我编写了以下程序来生成txt数据文件，这样所得的数据文件中数据储存的格式正是计算时需要的格式。

format longnn = input('请输入节点总数\n')%节点数nm = input('请输入单元数\n')%单元数L = input('请输入各杆件的长度，用分号隔开，单位为:m\n')%各单元长度（单位：米）E = input('请输入弹性模量，单位为:GPa\n')%弹性模量A = input('请输入杆件截面积，单位为:m^2\n')%截面积Q = input('请输入各杆倾角，用分号隔开，单位为:度\n')%各杆倾角JM = input('请输入各杆对应的节点(同一杆的节点序号用空格隔开，不同杆之间用分号隔开）\n')%各杆对应的节点(同一杆的节点序号用空格隔开，不同杆之间用分号隔开）nf=4%外载荷数量F=input('请输入外载荷，单位为:kN\n')%外载荷矩阵bc=input('请输入边界条件，被约束自由度序号的元素置1，未被约束置0')%边界条件，有位移约束的置1data=zeros(11,2*nn)data(1,1)=nndata(2,1)=nmfor i=1:nmdata(3,i)=L(i,1)enddata(4,1)=Edata(5,1)=Afor i=1:nmdata(6,i)=Q(i,1)endfor i=1:nmdata(7,i)=JM(i,1)endfor i=1:nmdata(8,i)=JM(i,2)enddata(9,1)=nffor i=1:2*nndata(10,i)=F(1,i)endfor i=1:2*nndata(11,i)=bc(1,i)endfilename=input('请输入保存的文件名\n','s')fid=fopen(filename,'w');fprintf(fid,'%f %f %f %f %f %f %f %f %f %f %f %f\n',data)fclose(fid)运行这个.m文件，在指引下输入各项数据，最后能在工作目录下生成文件名可以自定义的txt数据文件。

C 《有限单元法基本原理和数值方法》一书的源程序C********************************************************************* C *C PL ---- PROGRAM OF PLANE PROBLEM 96.1 *C *C********************************************************************* CC-------- 输入数据顺序--------C 1.NG 1整型C NG 结构结点总数C NG=0 则停止运行C 2.NE,MC,NX,NB,ND,EO,VO,T 5整型3实型C NE 结构单元总数C MC 计算控制类型参数C MC=0 平面应力C =1 平面应变C NX 作用载荷组数C NB 给定位移个数C ND 结构刚度矩阵的半带宽C EO 弹性模量C VO 泊松比C T 单元(结构)的厚度C 3.NWA NWE,NWK,NWP,NWD 5整型C 输出控制参数C =1 输出C =0 不输出C NWA 单元参数的输出控制参数C NWE 单元刚度矩阵的输出控制参数C NWK 结构刚度矩阵的输出控制参数C NWP 载荷向量的输出控制参数C NWD 结点位移的输出控制参数C 4.IJM(3,NE) 单元结点编码数组 3×NE整型C IJM(1,I); IJM(2,I); IJM(3,I)C 第I个三角形单元的结点编号,按结点编号顺序填写C 5.XY(2,NG) 结构结点坐标数组 2×NG实型C 6.MB(2,N,ZB(N 2×NB整型,NB实型C MB(1,I)---第I个给定位移所在的结点号C MB(2,I)=1--给定X方向位移C =0--给定Y方向位移C ZB(I)----给定位移值(以坐标正向为正)C 7.NF,NP 2整型C NF-----作用在结点上的集中载荷(坐标方向)的个数C NP-----作用均布侧压的单元边数C 若 NF>0 则填写C 8.MF(2,NF),ZF(NF) 2×NF整型,NF实型C MF(1,I)---第I个集中载荷所在的结点号C MF(2,I)=1--给定X方向集中力C =0--给定Y方向集中力C ZF(I)-----作用的集中力值C 若 NP>0 则填写C 9.MP(2,NP),ZP (NP) 2×NP整型,NP实型C MP(1,I)----第I个载荷作用边的起始结点号C MP(2,I)----第I个载荷作用边的起始结点号C ZP(I)------第I个均布载荷值CC 若 NX>1 重复 7.-9. (NX-1) 次CC 最后 NG=0 表示数据结束CC-------- 输出数据顺序--------C 1.IJM(3,NE) 单元结点编码数组 3×NE整型C IJM(1,I); IJM(2,I); IJM(3,I)C 第I个三角形单元的结点编号,按结点编号顺序填写C 2.XY(2,NG) 结构结点坐标数组 2×NG实型CC 若NWA=1,则输出C 3.I,B(7) 单元参数 NE行,1×NE整型,7×NE实型C 每行结构为:'NE='+单元号+Bi+Bj+Bm+Ci+Cj+Cm+A CC 若NWE=1,则输出C 4.IO,EK(6×6)单元刚度阵 NE行,1×NE整型,6×6×NE实型C 每行结构为:'NE='+单元号+EK(单元刚度阵)CC 若NWK=1,则输出C 5.SK(NT,ND) 结构刚度矩阵 NT×ND=2NG×ND实型CC 若NWD=1,则输出C 6.I,B 结点位移数据 NG行,1×NG整型,2×NG实型C 每行结构为:单元号+U+VCC 7.S1,S2,S3,X1,X2,CTAC 单元应力数据 6×NE实型C 分别代表σx,σy,τxy,σ1,σ2和主应力方向CC 若 NX>1 重复 6.-7. (NX-1) 次CC--------可调数组分配--------CC 实型数组 C(100000) 整型数组 IA(100000)C C(1) XY(2,NG) IA(1) IJM(3,NE)C C(N1) ZB(N IA(M1) MB(2,MC C(N2) BCA(7,NE) IA(M2) MF(2,NC C(N3) SK(NT,ND) IA(M3) MP(2,NP)C C(N4) F(NT) IA(MEND) 下限C C(N5) ZF(NF)C C(N6) ZP(NP)C C(NEND) 下限CC--------程序停止代码--------C 0 正常停止C 111 数组C越界C 222 数组C/IA越界C 333 单元面积非正C 444 结构刚度矩阵主元非正C---------------------------------------------------------------------- CC 主程序CDIMENSION C(500000),IA(50000),EK(36)CHARACTER*12 IN,OUTC IN和OUT为输入文件和输出文件的文件名WRITE(*,*)' 'WRITE(*,*)' PLEASE INPUT THE INPUT-FILE NAME (A<12)'WRITE(*,*)' 'READ(*,5) INC 输入输入文件的文件名WRITE(*,*)' 'WRITE(*,*)' PLEASE INPUT THE OUTPUT-FILE NAME (A<12)'WRITE(*,*)' 'READ(*,5) OUTC 输入输出文件的文件名5FORMAT(A12)OPEN(5,FILE=IN, STATUS='OLD')OPEN(6,FILE=OUT,STATUS='UNKNOWN')C 打开对应的输入和输出文件10READ(5,*) NGIF(NG.EQ.0) STOPC 输入结构结点数;如果结点数为0则停止运行READ(5,*)NE,MC,NX,NB,ND,EO,VO,TC 按顺序输入结构单元数,问题类型参数,载荷组数,给定位移个数C 结构刚度阵的半带宽,弹性模量,泊松比和结构厚度READ(5,*)NWA,NWE,NWK,NWP,NWDC 按顺序输入各输出控制参数NT=2*NGC 确定总刚度矩阵阶数NTCC 计算变界数组的下限CM1=3*NE+1M2=M1+2*NBN1=2*NG+1N2=N1+NBN3=7*NE+N2N4=N3+NT*NDN5=N4+NTC 得到各变界数组在一维大数组中的起始元素编号CC 检验实型数组C的下限CNEND=N5IF(NEND.LE.500000) GOTO 35WRITE(*,*)'*** EXCEED THE LIMIT OF ARRAY C(IN THE MIDDLE)!! ***'WRITE(*,30) NEND30FORMAT(/,'******** NEND=',I6,1X,'>80000 ********')STOP111C 若C下限超出500000,则给出错误信息并停止运行CC 数据输入C35CALL INPUT(NE,NG,NB,IA(1),C(1),IA(M1),C(N1))C 调用INPUT子程输入数据WRITE(*,40)40FORMAT(/10X,'##### INPUT PASSED #####')C 显示提示信息IF(MC.EQ.0) GOTO 45C 检验是否是平面应力问题CC 平面应变问题CE=EO/(1.0-VO*VO)V=VO/(1.0-VO)C 平面应变问题时,先进行弹性常数替换GOTO 50CC 平面应力问题C45 E=EOV=VO50 NX1=NXA1=E/(1.0-V*V)/4.0A2=0.5*(1.0-V)C 初始化NX1,A1和A2 / NX1为剩余载荷的组数CC 计算单元参数CCALL ABC(NE,NG,NWA,IA(1),C(1),C(N2))WRITE(*,55)55FORMAT(/10X,'##### ABC PASSED #####')C 调用ABC子程计算单元参数并显示提示信息CC 集成结构刚度矩阵KCDO60 I=N3,N4C(I)=0.060CONTINUEC 初始化结构刚度矩阵SKDO65 K=1,NEC 遍历结构的所有单元IO=KCALL KE(IO,NE,NWE,T,A1,A2,V,EK(1),C(N2))CALL SUMK(IO,NE,ND,NT,IA(1),C(N3),EK(1))C 调用KE子程计算出单元刚度阵并调用SUMK子程将其集成到结构刚度阵中65CONTINUEWRITE(*,70)70FORMAT(/10X,'##### SUMK PASSED #####')C 显示提示信息CALL CHECK(NT,ND,NWK,C(N3))C 调用CHECK子程检验结构刚度阵中的主元是否非正WRITE(*,75)75FORMAT(/10X,'##### CHECK PASSED #####')C 显示提示信息80READ(5,*) NF,NPC 输入集中载荷个数NF和均布载荷个数NPCC 再次计算变界数组的下限C 并检验实型数组C和整型数组IA的下限CM3=M2+2*NFN6=N5+NFNEND=N6+NP-1MEND=M3+2*NP-1C 计算C和IA的下限NM=0IF(NEND.LE.500000) GOTO 85WRITE(*,*)'*** EXCEED THE LIMIT OF ARRAY C (AT THE END)!! ***' WRITE(*,30) NENDNM=185IF(MEND.LE.50000) GOTO 95WRITE(*,*)'*** EXCEED THE LIMIT OF ARRAY IA (AT THE END)!! ***' WRITE(*,90) MEND90FORMAT(/,'******** MEND=',I6,1X,'>500 ********')STOP22295IF(NM.EQ.1) STOP222C 检验两个数组的下限,若下限超出则给出错误信息并停止运行CC 集成结构结点载荷列阵PCDO100 I=N4,N5C(I)=0.0100CONTINUEC 初始化结构结点载荷列阵FIF(NF.GT.0) CALL PF(NF,NP,NT,NWP,C(N4),IA(M2),C(N5))WRITE(*,105)105FORMAT(/10X,'##### PF PASSED #####')C 若集中载荷个数>0,调用PF子程集成各结点集中力并显示提示信息IF(NP.GT.0) CALL PP(NP,NT,NG,NWP,C(1),C(N4),IA(M3),C(N6))WRITE(*,110)110FORMAT(/10X,'##### PP PASSED #####')C 若均布载荷个数>0,调用PP子程集成各均布载荷并显示提示信息CC 引入给定位移CCALL DBC(NT,ND,NB,NX,NX1,C(N3),C(N4),IA(M1),C(N1))WRITE(*,115)115FORMAT(/10X,'##### DBC PASSED #####')C 调用DBC子程引入给定位移消除系数矩阵的奇异性,并显示提示信息CC 求解线性方程组KA=PCCALL GAUSS(NT,ND,NWD,NX,NX1,C(N3),C(N4))WRITE(*,120)120FORMAT(/10X,'##### GAUSS PASSED #####')C 调用GAUSS子程求解线性方程组并显示提示信息CC 计算单元应力CCALL STRESS(NE,NT,A1,A2,V,IA(1),C(N2),C(N4))WRITE(*,125)125FORMAT(/10X,'##### STRESS PASSED #####')C 调用STRESS子程输出各单元应力并显示提示信息NX1=NX1-1IF(NX1.GT.0) GOTO 80C 剩余载荷组自减1;若还有载荷剩余则继续计算GOTO 10C 重新输入ENDC********************************************************************* C 1 *C 子过程名称: INPUT *C 子过程功能: 按顺序输入并输出计算所需数据 *C *C********************************************************************* SUBROUTINE INPUT(NE,NG,NB,IJM,XY,MB,ZC 形参说明C 输入:C NE 整型,结构单元总数C NG 整型,结构结点总数C NB 整型,给定位移的个数C IJM(3,NE) 整型,单元结点编码数组C XY(2,NG) 实型,结构结点坐标数组C MB(2,N 整型,位移约束信息数组C ZB(N 实型,位移约束数值数组DIMENSION IJM(3,NE),XY(2,NG),MB(2,N,ZB(NCREAD(5,*) ((IJM(I,L),I=1,3),L=1,NE)C 输入单元结点编码数组READ(5,*) ((XY(I,J),I=1,2),J=1,NG)C 输入结构结点坐标数组READ(5,*)((MB(I,L),I=1,2),L=1,N,(ZB(L),L=1,NC 输入位移约束信息和数值数组WRITE(6,20)((IJM(M,I),M=1,3),I=1,NE)20FORMAT(1X,4(3I4,3X),3I4)C 输出单元结点编码数组WRITE(6,40)((XY(M,I),M=1,2),I=1,NG)40FORMAT(1X,6E12.5)C 输出结构结点坐标数组RETURNENDC********************************************************************* C 2 *C 子过程名称:ABC *C 子过程功能:根据各单元的结点坐标, *C 计算并输出所有单元的各参数 *C *C********************************************************************* SUBROUTINE ABC(NE,NG,NWA,IJM,XY,BCA)C 形参说明C 输入:C NE 整型,结构单元总数C NG 整型,结构结点总数C NWA 整型,单元参数输出控制,0－不输出,1－输出C IJM(3,NE) 整型,单元结点编码数组C XY(2,NG) 实型,结构结点坐标数组C 输出:C BCA(7,NE) 实型,结构单元参数数组C 变量说明C X(2,5) 实型,当前计算单元的结点坐标数组C B(7) 实型,当前计算单元的单元参数数组DIMENSION IJM(3,NE),XY(2,NG),BCA(7,NE),X(2,5),B(7)CIF(NWA.EQ.1) WRITE(6,5)5FORMAT(/10X,'PARAMETERS OF ELEMENTS BCA(7,NE)'/)C 若控制打开,则输出单元参数提示信息DO80 I=1,NEC 遍历所有单元DO10 K=1,3C 遍历单元内的3个结点K1=IJM(K,I)C 取结点在结构中对应的结点号DO10 J=1,2X(J,K)=XY(J,K1)C 得到当前结点的坐标值10CONTINUEDO20 J=1,2X(J,4)=X(J,1)X(J,5)=X(J,2)20CONTINUEC 为编程方便,每单元多存两个结点坐标DO30 K=1,3B(K)=X(2,K+1)-X(2,K+2)C 计算BmB(K+3)=X(1,K+2)-X(1,K+1)C 计算Cm30CONTINUEB(7)=(B(1)*B(5)-B(4)*B(2))*0.5C 计算单元面积AIF(NWA.GT.0) WRITE(6,40)I,B40FORMAT(1X,'NE=',I3,/3X,7E10.4)C 若输出控制打开,则输出单元号和对应的参数IF(B(7).LE.0.0) GOTO 60C 若当前单元面积为负,则出错DO50 J=1,7BCA(J,I)=B(J)50CONTINUEC 将当前单元参数数组中的值传送给输出数组GOTO 8060WRITE(6,70)I,(IJM(J,I),J=1,3)70FORMAT(/5X,'ELEMENT',I5,5X,'AREA IS NONPOSITIVE',5X,'IJM',3I5) STOP333C 显示出错信息并停止运行80CONTINUERETURNENDC********************************************************************* C 3 *C 子过程名称:KE *C 子过程功能:根据结构单元参数数组,计算出 *C 指定单元的单元刚度矩阵 *C *C********************************************************************* SUBROUTINE KE(IO,NE,NWE,T,A1,A2,V,EK,BCA)C 形参说明C 输入:C IO 整型,计算单元编号C NE 整型,结构单元总数C NWE 整型,刚度矩阵输出控制,0－不输出,1－输出C IJM(3,NE) 整型,单元结点编码数组C XY(2,NG) 实型,结构结点坐标数组C T 实型,单元厚度C A1 实型,材料系数,A1=E/(4*(1-V**2))C A2 实型,材料系数,A2=(1-V)/2C V 实型,泊松比C BCA(7,NE) 实型,结构单元参数数组C 输出:C EK(6,6) 实型,单元刚度矩阵DIMENSION B(7),BCA(7,NE),EK(6,6)DO10 I=1,7B(I)=BCA(I,IO)10CONTINUEC 得到计算单元的参数A=A1/B(7)*TDO20 I=1,3DO20 J=I,3C 将单元刚度矩阵分块成3×3个子矩阵I1=2*IJ1=2*JEK(I1-1,J1-1)=A*(B(I)*B(J)+A2*B(I+3)*B(J+3))EK(I1-1,J1)=A*(V*B(I)*B(J+3)+A2*B(I+3)*B(J))EK(I1,J1-1)=A*(V*B(I+3)*B(J)+A2*B(I)*B(J+3))EK(I1,J1)=A*(B(I+3)*B(J+3)+A2*B(I)*B(J))C 计算每个子矩阵各元素的值20CONTINUEDO30 I=3,6DO30 J=1,IEK(I,J)=EK(J,I)30CONTINUEC 根据对称性得到左下角矩阵的值IF(NWE.EQ.0) GOTO 60C 若输出控制关闭,则直接结束子过程WRITE(6,40) IO40FORMAT(/1X,'EK NE=',I5)WRITE(6,50)EK50FORMAT(1X,6E11.4)C 输出单元号和对应的刚度矩阵60RETURNENDC********************************************************************* C 4 *C 子过程名称:SUMK *C 子过程功能:将指定单元的单元刚度矩阵集成到结构刚度 *C 矩阵中,结构刚度矩阵以二维等带宽方式存储 *C *C********************************************************************* SUBROUTINE SUMK(IO,NE,ND,NT,IJM,SK,EK)C 形参说明C 输入:C IO 整型,计算单元编号C NE 整型,结构单元总数C ND 整型,结构刚度矩阵的半带宽C NT 整型,结构刚度矩阵的阶数C IJM(3,NE) 实型,单元结点编码数组C EK(6,6) 实型,单元刚度矩阵C 输出:C SK(NT,ND) 实型,结构刚度矩阵DIMENSION IJ(3),SK(NT,ND),IJM(3,NE),EK(6,6)CDO10 I=1,3IJ(I)=IJM(I,IO)10CONTINUEC 取出计算单元的结点号DO20 I=1,3DO20 J=1,3C 遍历单元刚度矩阵的3×3个子矩阵IF(IJ(I).GT.IJ(J)) GOTO 20C 如果是下三角元素,则不储存M=2*IJ(I)-1N=2*(IJ(J)-IJ(I))+1C 得到在结构刚度矩阵中等带宽储存的行列码MO=2*I-1NO=2*J-1C 得到对应在单元刚度矩阵中的行列码SK(M,N)=SK(M,N)+EK(MO,NO)SK(M,N+1)=SK(M,N+1)+EK(MO,NO+1)SK(M+1,N)=SK(M+1,N)+EK(MO+1,NO+1)IF(IJ(I).EQ.IJ(J)) GOTO 20C 不存储主子块的下三角元素SK(M+1,N-1)=SK(M+1,N-1)+EK(MO+1,NO)C 将单元刚度矩阵的元素叠加到结构刚度矩阵中20CONTINUERETURNENDC********************************************************************* C 5 *C 子过程名称:CHECK *C 子过程功能:检验结构刚度矩阵的主元并输出结构刚度矩阵, *C 如果主元非正,则输出错误信息并停止运行 *C *C********************************************************************* SUBROUTINE CHECK(NT,ND,NWK,SK)C 形参说明C 输入:C NT 整型,结构刚度矩阵阶数C ND 整型,结构刚度矩阵半带宽C NWK 整型,结构刚度矩阵输出控制,0－不输出,1－输出C SK(NT,ND) 实型,结构刚度矩阵C 变量说明C M 整型,错误主元个数DIMENSION SK(NT,ND)CIF(NWK.EQ.0) GOTO 30WRITE(6,20) ((SK(I,J),I=1,NT),J=1,ND)20FORMAT(1X,5E13.6)C 若输出控制打开,则输出结构刚度矩阵30 M=0C 置错误个数为0DO50 I=1,NTC 遍历结构刚度矩阵各主元IF(SK(I,1).GT.1E-10) GOTO 50WRITE(6,40)I,SK(I,1)40FORMAT(/10X,'MAIN ELEMENT IS NONPOSITIVE NT=',I4,5X,E12.6)M=M+1C 若主元非正,则输出错误信息,并使错误个数+150CONTINUEIF(M.GT.0) GOTO 60GOTO 7060STOP444C 若错误个数>0则停止运行70RETURNENDC********************************************************************* C 6 *C 子过程名称F *C 子过程功能:将结点集中力装入等效载荷列阵, *C 同时输出等效载荷列阵 *C *C********************************************************************* SUBROUTINE PF(NF,NP,NT,NWP,F,MF,ZF)C 形参说明C 输入:C NF 整型,坐标方向上集中载荷的个数C NP 整型,作用均布侧压的边数C NT 整型,等效载荷列阵元素个数C NWP 整型,载荷输出控制,0－不输出,1－输出C MF(2,NF) 整型,作用于结点上集中载荷的信息数组C MF(1,I):第I个载荷作用的结点号C MF(2,I):第I个载荷作用的方向,0－Y向,1－X向C ZF(NF) 实型,作用于结点上集中载荷的值C 输出:C F(NT) 实型,等效载荷列阵DIMENSION MF(2,NF),ZF(NF),F(NT)CREAD(5,*) ((MF(I,L),I=1,2),L=1,NF),(ZF(L),L=1,NF)C 输入集中载荷信息和数值数组DO40 I=1,NFC 遍历所有集中载荷N=2*MF(1,I)-MF(2,I)C 得到载荷在结构等效载荷列阵中的对应顺序F(N)=F(N)+ZF(I)C 将集中载荷叠加到结构等效载荷列阵中40CONTINUERETURNENDC********************************************************************* C 7 *C 子过程名称P *C 子过程功能:计算均布侧压的等效结点载荷并装入等效载荷列阵, *C 同时输出等效载荷列阵 *C *C********************************************************************* SUBROUTINE PP(NP,NT,NG,NWP,XY,F,MP,ZP)C 形参说明C 输入:C NP 整型,作用均布侧压的边数C NT 整型,等效载荷列阵元素个数C NG 整型,结构结点总数C NWP 整型,载荷输出控制,0－不输出,1－输出C MP(2,NF) 整型,作用于单元边上均布载荷的信息数组C MF(1,I):第I个载荷作用边的起始结点号C MF(2,I):第I个载荷作用边的终止结点号C ZP(NF) 实型,作用于单元边上均布载荷的值C 输出:C F(NT) 实型,等效载荷列阵DIMENSION MP(2,NP),ZP(NP),XY(2,NG),F(NT)CREAD(5,*) ((MP(I,L),I=1,2),L=1,NP),(ZP(L),L=1,NP)C 输入均布载荷信息和数值数组DO40 I=1,NPC 遍历所有均布载荷N1=MP(1,I)C 得到载荷的起始结点号N2=MP(2,I)C 得到载荷的终止结点号PX=XY(2,N1)-XY(2,N2)PY=XY(1,N2)-XY(1,N1)PX=.5*ZP(I)*PXPY=.5*ZP(I)*PYC 得到等效结点载荷 PX=qt(Yi-Yj)/2,PY=qt(Xi-Xj)/2F(2*N1-1)=F(2*N1-1)+PXF(2*N1)=F(2*N1)+PYF(2*N2-1)=F(2*N2-1)+PXF(2*N2)=F(2*N2)+PYC 将均布载荷的等效结点载荷叠加到结构等效载荷列阵中40CONTINUERETURNENDC********************************************************************* C 8 *C 子过程名称BC *C 子过程功能:引入给定位移的边界条件,消除系数矩阵的奇异性 *C *C********************************************************************* SUBROUTINE DBC(NT,ND,NB,NX,NX1,A,B,MB,ZC 形参说明C 输入:C NT 整型,系数矩阵阶数C ND 整型,系数矩阵的半带宽C NB 整型,给定位移的个数C NX 整型,载荷的总组数C NX1 整型,载荷的剩余组数C A(NT,ND) 实型,系数矩阵 (兼输出)C B(NT) 实型,等效结点载荷列阵 (兼输出)C MB(2,N 整型,给定位移的信息数组C MB(1,I):第I个给定位移的结点号C MB(2,I):第I个给定位移的方向,0－Y向,1－X向C ZP(N 实型,给定位移的值DIMENSION MB(2,N,ZB(N,A(NT,ND),B(NT)CDO60 I=1,NBC 遍历所有给定位移N=2*MB(1,I)-MB(2,I)C 取要修改的方程的序数Z=ZB(I)C 取对应位移的值IF(ABS(Z).LT.1E-10) GOTO 20C 若位移为0,用对角元素改1法,否则用对角元素乘大数法IF(NX.NE.NX1) GOTO 10C 若不是第1组载荷,则只在B中引入给定位移A(N,1)=A(N,1)*1E+1510 B(N)=A(N,1)*ZC 对角元素乘大数,并修改对应的等效结点载荷GOTO 60C 以下为对角元素改1法20IF(NX.NE.NX1) GOTO 50C 若不是第1组载荷,则只在B中引入给定位移A(N,1)=1.0C 将对角元素置1DO30 J=2,NDA(N,J)=0.030CONTINUEC 将对应行元素置0DO40 K=2,NDIF(N.LT.K) GOTO 50M=N-K+1A(M,K)=0.040CONTINUEC 将对应列元素置050 B(N)=0.0C 等效结点载荷置060CONTINUERETURNENDC********************************************************************* C 9 *C 子过程名称:GAUSS *C 子过程功能:使用高斯消元法求解线性方程组 *C *C********************************************************************* SUBROUTINE GAUSS(NT,ND,NWD,NX,NX1,A,C 形参说明C 输入:C NT 整型,结构刚度矩阵阶数C ND 整型,结构刚度矩阵半带宽C NWD 整型,载荷向量输出控制,0－不输出,1－输出C NX 整型,载荷的总组数C NX1 整型,载荷的剩余组数C A(NT,ND) 实型,系数矩阵C B(NT) 实型,等效结点载荷列阵 (兼输出)CDIMENSION A(NT,ND),B(NT)CCC 消去过程CN=NT-1IF(NX.EQ.NX1) GO TO 10ASSIGN 50 TO MGO TO 2010 ASSIGN 30 TO M20DO60 K=1,NM1=ND-1IF((M1).GT.(NT-K)) M1=NT-KDO60 L=1,M1C=A(K,L+1)/A(K,1)IF(ABS(C).LT.1E-18) GO TO 60GO TO M,(30,50)30 M2=ND-LDO40 J=1,M2A(K+L,J)=A(K+L,J)-C*A(K,J+L)40CONTINUE50 B(K+L)=B(K+L)-C*B(K)60CONTINUECC 回代过程CB(NT)=B(NT)/A(NT,1)DO80 K=1,NI=NT-KM1=NDC=B(I)IF((K+1).LT.(ND)) M1=K+1DO70 J=2,M1L=I+J-1C=C-A(I,J)*B(L)70CONTINUEB(I)=C/A(I,1)80CONTINUECC 输出求解结果CIF(NWD.EQ.0)GOTO 150N=NT/2N11=N/2IF(FLOAT(N11)+.3-FLOAT(N)/2.0.GT.1E-7)N=N-1DO140 I=1,N,2NT2=NT/4NT3=NT/2FL1=FLOAT(NT2)+.3-FLOAT(NT3)/2IF((FL1.LT.0).AND.(I.EQ.N)) GOTO 120J=2*I-1K=2*II1=I+1J1=J+2K1=K+2WRITE(6,110) I,B(J),B(K),I1,B(J1),B(K1) 110FORMAT(1X,2(I3,3X,E11.5,2X,E11.5,5X))GOTO 140120 J=2*I-1K=2*IWRITE(6,130)I,B(J),B(K)130FORMAT(1X,I3,3X,E11.5,2X,E11.5)140CONTINUE150RETURNENDC********************************************************************* C 10 *C 子过程名称:STRESS *C 子过程功能:根据结构各结点的位移,输出各单元的 *C 单元应力,主应力和应力主方向 *C *C********************************************************************* SUBROUTINE STRESS(NE,NT,A1,A2,V,IJM,BCA,F)C 形参说明C 输入:C NE 整型,结构单元总数C NT 整型,结点位移列阵元素个数C A1 实型,材料系数,A1=E/(4*(1-V**2))C A2 实型,材料系数,A2=(1-V)/2C V 实型,泊松比C IJM(3,NE) 整型,单元结点编码数组C BCA(7,NE) 实型,结构单元参数数组C F(NT) 实型,结构结点位移列阵C 变量说明C B(7) 实型,当前计算单元的单元参数数组C R(6) 实型,当前计算单元的结点位移数组C A 实型,A=E/(2*(1-V**2)*Ae)C S1,S2,S3 实型,当前计算单元的应力σx,σy,τxyC X1,X2 实型,当前计算单元的主应力σ1,σ2C CTA 实型,当前计算单元的主应力方向DIMENSION IJM(3,NE),BCA(7,NE),F(NT),B(7),R(6)CWRITE(6,5)5FORMAT(/,10X,' ',/)DO60 I=1,NEC 遍历所有单元S1=0.S2=0.S3=0.C 当前单元的应力值初始化DO20 J=1,7B(J)=BCA(J,I)C 取当前单元的单元参数20CONTINUEA=2*A1/B(7)C 得到E/(2*(1-V**2)*Ae)DO30 J=1,3C 遍历单元内的3个结点N=IJM(J,I)*2R(2*J-1)=F(N-1)R(2*J)=F(N)C 取对应结点的结点位移30CONTINUEDO40 J=1,3K=2*JS1=S1+A*(B(J)*R(K-1)+V*B(J+3)*R(K))S2=S2+A*(V*B(J)*R(K-1)+B(J+3)*R(K))S3=S3+A*A2*(B(J+3)*R(K-1)+B(J)*R(K))40CONTINUEC 计算当前单元的应力C σx=A*Σ(Bi*Ui+V*Ci*Vi),σx=A*Σ(V*Bi*Ui+Ci*Vi)C τxy=A*Σ(A2*Ci*Ui+A2*Bi*Vi)P=.5*(S1+S2)Q=.5*(S1-S2)X1=P+SQRT(Q*Q+S3*S3)X2=2*P-X1C 计算当前单元的主应力CTA=0.IF (S3.GT.0) CTA=ATAN((X1-S1)/S3)C 计算当前单元的主应力方向WRITE(6,50) S1,S2,S3,X1,X2,CTA50FORMAT(1X,E10.4,2X,E10.4,2X,E10.4,2X,E10.4,2X,E10.4,2X,F8.4) C 按顺序输出当前单元的σx,σy,τxy,σ1,σ2和主应力方向60CONTINUERETURNEND。

有限单元法大作业问题：一薄板，尺寸为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计算结果输出如下。

有限元作业1.输出1-100之间能被3整除的数,并求这些数的累计之和。

解答：在workplace中输入如下代码：READ*,NN=100SUM=0.0DO I=1,NIF (MOD(I,3)==0) THENSUM=SUM+IWRITE(6,*)IELSESUM=SUM+0END IFEND DOWRITE(6,*) SUMEND按Ctrl+F5键，输入任意一个N，因为N为给定的100，无论输入任意数据都行，在界面中输入数字1，得出结果为：2.某大学按学分收费,如果不超过12学分,应交纳学费4000元,如果超过12学分,则每超过1个学分,加收500元学费,输出学分和应缴纳的学费.解答：在workplace中输入如下代码：READ*,NSUM=4000DO I=1,NIF(I>=12)THENSUM=4000+(I-12)*500ENDIFWRITE(6,*)I,SUMENDDOEND按Ctrl+F5键，输入任意一个N ，在界面中输入数字20，得出结果为：3.用迭代法求解方程：038223=++-x x x 解答：在workplace 中输入如下代码：READ *,X0 N=0 X0=3 DO I=1,NF1=X*X*X-2*X*X+8*X+3 F2=3*X*X-4*X+8 X=X0-F1/F2 N=N+1IF (ABS (X-X0)>=0.00001)THEN X=X ENDIF ENDDO WRITE (6,*)X,F1 END按Ctrl+F5键，输入任意一个N ，在界面中输入数字0.02，得出结果为：4. 编写任意大小的两个矩阵相乘的程序，并算一具体实例5. 哥德巴赫提出，一个不小于6的偶数必定能表示为两个素数之和，请将6-100之间的全部偶数表示为两个素数之和6. 已知2)sinh(xx e e x --=利用级数展开式求出x=1~5时的sinh(x)7.已知矩阵，求出该矩阵对角线上元素的和⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=359162173895079684835296153********2A 8. 利用Fortran 语言编译一程序打印乘法口诀表 9. 设计一个程序计算：S=9!+ 7!+ 5!+3!+1! 10. 编译一程序求逆矩阵并用一实例验证 11. 不限方法求积分⎰10dxe x12. 已知杆的E I =6，画连续梁的内力图kN.m313.画刚架的内力图14.已知各杆的弹性模量和横截面面积均相同，为E=29.5´104N/mm2, A=100mm2,求各杆的内力。