有限元分析程序设计

结构有限元分析程序设计绪 论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边界条件及程序设计（等效荷载计算，位移边界条件置入，多工况的对称性）. 总刚线性方程组求解（LDLT分解，分块算法，子结构算法，波前法） 单元应力计算+应力处理与改善。. 数据处理（数据分类，压缩存贮，节点排序方法）, 变带宽存贮的优化理论，图论的理论，有限元的图结构，存贮管理复核）. 有限元议程全稀疏管理与求解策略。说明： 仅线性部份，复材，接触，弹塑性等不包括，基本部份。实践性作业安排：1. 作业： 总的结构管理程序+子功能模块的编程，一个题的计算实践2. 送有限元元素库。0.4 课程要求1. 先行要求2. 作业要求（计算机编程得出正确结果）3. 课程校核要求0.5 基本方法复习0.5.1 结构化程序设计方法0.5.2 有限元分析方法回顾0.5.3 Fortran语言回顾0.5.4 结构化程序设计1. 基本结构：构成一个问题从输入到求解输出的基本程序形式：ImputProcessOutput(输入) （处理） （输出）三种基本形式：a). 树形结构（顺序执行结构） Imputp1 . pn Output (多用于程序结构 call1,call2.)b). t选择结构（条件执行结构） 内部算法实现: IF.THEN; SWHICH,CASEc).循环结构(反复执行结构)特点: 结构特征简单明了,易读,易调试. 尽量少用GO TO 语句2. 整体结构(算法语言系统结构)积木式(Fortran): 每个设定的功能分析团体的一个模块,每个模块又称作整体结构的素材,主模块象积木一样堆积.语体不联系,但有通讯方法沟通模块间信息,各种模块有各自的特征语体,main progame,subroutine sub. 嵌套式(Pascall): 主模块与子模块相互嵌套,各模块的特征,语体相同 procedure main procefure ssub1 . procedure sub2 . End sub1 end sub2 . End main 函数式(c语言): 主要特点是功能模块作为库函数调用,需用时在库内调用,每一个函数有表征语句,这种语体接近自然思维,而且对系统资源的调配应用更完善. 面向对象的程序设计: 实施过程的可视化+控制性 3. 结构化程序设计方法a). TopDown(自上面下),系统性强,选择性强.b). Critical Component First (关键部份优先),先抓主要矛盾,分清重,缓,急.c). 独立调试,总体联调,(软件设计的社会化作业).4. 程序设计要点 a). 自觉有意识地设计一个良好的程序结构,做到:易读,易懂,易管理,易修改,易发展. b). 做到逻辑清晰,说明完整. c). 要有工艺设计概念有框图,有步骤.5). 结构化程序设计原则 a). 尽可能通用性好（适应各种规模的复题，？的扩大依据程序设计指标而定） b).整体精炼，清晰；避免GOTO。c).省机时，省存贮，计算精度高，（算法上下功夫，要理论分析加技巧）d).输入数据少，格式简单。e).输出结果简明，忌讳打印过多（与具体调试过程不一样）。f).易读易维护，易发展。0.5.2.有限无方法求解过程回顾一. 力学模型的分级管理有限无程序对力学模型的数据按一级：结构级（有点广义，不仅指具体结构，也指模型题目的规模） 二级：单元级三级：节点级2.基本关系3. 描述参数 A). 节点描述参数 (1). 节点位置(总体坐标系下的坐标). (2). 节点局部坐标(按节点的约束方向制定的特殊坐标系x，y，z,v如斜支撑) (3). 节点的性质(自由,固定,指定位移,从属其它节点). (4). 节点力:(Fx ,Fy , Fz ,Mx ,My ,Mz ) (5). 节点位移:(u , v , w, x ,y ,z )B. 单元描述参数 (1). 材料特性参数不清 E, G, D (2). 节点的几何刚度参数(即面积A,板厚H,梁抗弯模量I) (3). 单元的局部坐标. (用于应力分析等,如图形曲面) (4). 单元的节点编号 (5). 单元的几何矩阵营 (节点变形与应力关系矩阵) (6). 单元刚度矩阵 K (7). 单元的应力,应变向量,(有限元分析多用向量,而不用矩阵(张量)结构描述参数 单元总数,节点总数,单元娄型总数,结构材料种娄数,节点自由度数(控制题目规模)二. 基本公式系统1. 单元刚度计算公式 2. 单元刚阵组合 K=ATKA 3. 单元节点荷载计算4. 节点荷载组装: 5. 位移约束关系:6. 总刚方程解:7. 应变计算:8. 应力计算:9. 支撑反力计算:输入边界条件(对称条件)形成各荷载工况的节点荷载阵 总刚分解 回代求出位移及输出 计算应变、应力形成单元刚阵单刚向总刚投放坐标变换输入原始参数计算总刚规模形成总刚方程向总节点荷载阵投放形成单元荷载阵调整几何、弹性矩阵调整单元位移列阵三. 有限元分析的模块组织.四 结构分析的原始输入数据 1题目规模 节点数目：NNP 单元数目： NE 2节点数据 单元人坐标：XE（NNP，3） 3单元数据 单元节点编号 ： ME（NE，3）、ME（NE，2） 材料特性： E、N 单元几何参数： I、RI（惯矩） 4荷载数据 外荷载作用点，坐标及大小： PA（NNP，1）0. 5.3 Fortran语言回顾1子模块（子程序） subroutinea. 特点：独立性强，只要输入输出接口，象一个黑匣子，与外界无关。b. 作用：完成一个独立的功能（求应力，矩阵分解，投放等）c. 格式： subroutine function(ip1,ip2,rp1,rp2,io1,io2,ro1,ro2).(其中ip1,ip2,rp1,rp2,是输入形参，io1,io2.rp1,rp2是输出形参)2数据传递形式 1）. COMMON 公共块语句传递，（公共块的内容不能作为形参）a. 公共块分为无名公共块和有名公共块b. 公共块的参数不能作为子程序的参数出现，c. 公共块名一致，其内容在不同公共块中可以标志符不同（但其长度应一致）d. 通用原始数据放入公共块（作为实参错误率大）e. 尽可能不放数值，安息组一般可作成可调长度 f. 格式 Common/comm/. Subroutine fun() Common/comm/.2).形参实参对应a. 实参不能开辟存贮单元，子程序内定义语句中的形参数组由主程序定义，在子程序中仅形式定义（即仅说明是数组，因而大小无所谓）b. 格式：Dimension RP(1000),RO(1000) . Call sub1(RP,RO) END DIMENTION IBANK SUBROUTINE SUB1(RP,RO,NE)DIMENTION RP(1),RO(1),SP(50)DIMENTION RP(NE,1),RO(NE,1) (形参的动态定义，实参不能)3）. 数组长度自动调整方法。 PROGRAM MAIN INPLICIT REAL*8 (AH,O-Z) CHARACTAR*20 TR COMMON/COMM/. DIMENTION IBANK( ),RBANK( ),IP1( ),IP2( ) IP1(1)= IP(N)=. IP2(1)= IP2(N)=. CALL SUB1(IBANK(IP1(1),IBANK(IP1(N),RBANK(IP2(1),.) . END SUBROUTINE SUB1(II1,IO2,.RI1,.RO1.NE) DIMENTION II1(1) ,IO2(NE,1), RI1(1), RO1(1)*.公共块分无名与有名公共块：公共块参数不能子程序形参中出现。公共块一致，其内容在不同程序中可以标识符不同。（但长度要一致）*.公共块使用原则：1) 通用原始数据代入公共块，作为实参出错率大2) 尽可能不放数组，数组要作成可调的（自动可调）。数组长度自动调整方法。1.单元刚阵几元素设计主讲内容：工程常用元素的单元刚阵计算与编程。1.1杆、梁的单元刚阵1.2平面问题的三角形、矩形单元刚阵1.3板弯问题的三角形、矩形单元刚阵1.4非、拟协调之介绍及分片检验方法1.0本讲内容序1.0.1单元刚阵在总刚中的地位和作用由有限元分析过程可以知道，最主要的两步骤：1) 进行单元的特性分析，建立单元级的刚度方程（平衡方程），得到单元的刚度矩阵；注：中各元素的位置与和要有严格的对应关系2) 进行结构整体分析，集合所有单元的刚度方程建立全结构的刚度方程，设得到全结构的刚度矩阵。即A为布尔矩阵总刚由单元刚阵按节点编号顺序投放而成，故单元刚阵是整个有限元方法的第一步；是总装刚阵的原材料。解释：单元节点平衡；结构总刚对号和迭加的原理1.0.2单元刚阵形式1. 基本公式：（1.1）成为几何矩阵，作用是把应变与位移联系起来反映了单元的形变特性与性质。固而是坐标的函数，规模与节点数目（自由度数目）有关，不同类型的单元有不同的。（依据设定的位移形函形式与节点位移的性质）。对常应元是常数。为弹性矩阵，只取决于材料特性E，与坐标无关，对于平面应力问题： （1.2）对于平面应变问题，将上式E;代入上式即得。2. 形成方法1) 直接赋值。有了积分出来的 各元素显式，即可逐项（元素）直接赋值。在赋值时，有数字变化规律的应按特点赋值，以缩短程序长度，提高执行效率。2) 矩阵乘法。用矩阵乘法进行赋值。1.1杆、梁的单元刚阵1.1.1杆元xy对于各种变截面等轴力杆其单元刚阵显示：uj(Pj)y（1.3）为等效截面面积。ui(Pi)x对于不同的截面变化情况可由材料力学方法计算得到。*对于等截面杆，可设定形状函数，利用公式（1.1）计算，即由：,(1.4) (1.5)*对于变截面杆，采用假设形状函数的方法不妥（误差大），因此，常用材料力学推倒方法：即公式： (等轴力)（1.6）相对于1点：同理，相对于2点：于是得：（1.6）令：故得：（1.3）式思考题：变截不均的弯轴力杆的刚阵形式如何？对于满足杆元截面积随坐标单调变化的函数，Aeq计算见P2830*形成方法可采用直接赋值语句完成见P30FAKE11.1.2梁元应用公式（1.1）进行直接积分，得到显式，然后用上赋值语句编程。1. 梁元形 及刚阵显式工程上使用的梁元一般为复合受力状态，包括轴向、扭转及各平面内的y弯曲，因此对于一般的三维空间中的梁元节点位移包括：线性位移、及转角位移、；节点力 位包括：三个方向的力、和三个方向的弯矩。yxzxz故梁元每个节点有6个自由度，即：1) 对于轴向拉伸与扭转可采用线性形函描述即：(1.7)应变定义（1.8）截面极惯矩（1.9）注意：这里的体积分、面积分为常数包括在弹性阵内(1.10)2) 在xy平面内的弯曲，不计剪切，节点位移为；节点力为形 的选取是依据弯曲理论（材料力学）的基本假设的两直法线公设所确定即：（1.11）进行积分，可得三次挠度曲线（1.12）代入节点位移系统可得到形状函数（1.13）（1.14）*应变定义：梁的曲率定义为广义定义应变，梁的弯矩定义为广义应力，即：（1.15）按梁的挠度方程，可有广义的应力应变关系(1.16)其中代入（1.13式），得：（1.17）（1.18）由式（1.1），于是：（1.19）（1.20）3) 在yz平面内的弯曲，节点位移与力系统为：形函推导过程类似。刚阵显式：（1.21）4) 综合上述三种单独刚度，可以得空间梁元素在局部坐标系下复合状态的刚阵显示见P4243间，即书公式（224），（注意节点位移、节点力的排序）。2. 形成梁元刚阵的子程序设计方法：直接赋值。见P4344FAKE41.2平面问题的三角元、矩形单元刚阵1.2.1平面应力三角元。（常应变元）1. 公式。(1.22)Qy(1.23)m其中：在板平面内TjR转换123iPx2. 形成B阵的轮换循环赋值方法。I=1,3I1=I+1-3*(I/3)I2=I+2-3*(I+1)/3)KI=ME(K,I1)KJ=ME(K,I2)B(1,2*I-1)=XE(KI,2)-XE(KJ,2)B的第1行1,3,5之元素B(2,2*I)=XE(KJ,1)-XE(KI,1)B的第2行2,4,6之元素B(3,2*I-1)=B(2,2*I)PP34B(3,2*I)=B(1,2*I-1)未完当平面元处于空间位置时，（即节点坐标为三维数据）。平面内的几何参数需要计算，计算公式：y1） PQ长度：mQT2） PQ的方向余弦：RjiPx3） PT长度（按线段投影公式）4） TR长度：5） QT长度：6） TR的方向余弦7） 面积8） 刚度阵中的参数9） 在中各长度的取值记号规则：利用上面的子程序，可直接赋值形成局部坐标下的阵。3. 局部系下刚阵计算与编程。4（1.24）31mP 其中C为板厚利用矩阵乘法可直接运算得到间见P351.2.2平面应力四节点单元刚度i1 单元平面内局部下形状函数。2j节点位移和节点力系统：形状函数：（1.25）（1.26）（1.27）显见得不再是常数阵，说明应变在平面内线性变化。2.矩阵平面应力刚阵形成公式：由于不再是 常数阵，只能手工积分每一元素，设得到刚阵显式，然后用赋值方法编程。见公式PP38-39（2-19）刚度形成子程序FAKE3见PP39-40思考题：变轴力、变截面杆单元刚阵的形成方法。（用柔度方法）作业：（1）阅读有限元PP27-60内容（2）有一空间梁元，截面尺寸：长度1m. 编程计算单元刚度矩阵。1.3.1板弯问题的三角元、矩形单元刚阵作为作业、阅读理解，不作掌握要求 1.4 协调、非协调、广义协调及分片检验1、4、0 引 以有限元数值分析的技术实现为目的本门课程，不仅要求学生能够进行实际的工程运算；另一方面也需要对解的收敛及精确性有所了解，是能从细节计算到理论性质都有所把握，这样，才能做到全面深入有助于对解结果得理论分析，此为基本之目的。1、4、1 协调、非协调介绍 位移法有限元以Ritz的结构最小有限元为基础，该原理在数学上是一个泛函极值（变分）问题，系统势能可以表为以下数学形式：面力势能内力势能体力势能=1/2 - - (1) 0 。表述为：在所有满足内部连续性和运动学边界条件的位移中满足平衡方程的位移使系统势能取驻值。如果驻值是极小点的，则平分行是稳定阶。又：对于精确于问题的位移函数，系统势能的变分可求得关于问题应满足的所有微分方程：平衡方程边界条件（几何关系及物理方程是自然满足的）遗憾的是精确位移难得寻找，故一般采用泛函的极小化序列逼近方法。类似于傅立叶级数逼近函数那样，把无穷维空间用有限空间去逼近。在有限元当中，当元素尺寸趋近于0时（即节点数目或节点自由度数趋于时），最后的解答若能无限逼近准确解，那么这样的位移函数（或形状函数）就称为收敛的，因此从收敛性及算收敛速度方面提出几点对形状函数的要求：、函数本身及其导数应在元素上连续，并含有常数部分；、元素之间的位移协调，不仅节点处的位移应当协调，沿整个内边界上的位移也应当协调（或称相容 ）。 、多项式的项数越多越好，因用高次比低次多项式收敛快。、含有刚体位移（平动包含常数项，转动包含线性项）。协调之： 即满足、条件的形状函数的元素，当然能满足3) 4)条件协调 元的收敛率就更高。 协调元的性质：1) 能够以单调趋势逼近于正确解。如曲线.2) 势能总是大于最小状态，故解得上界。3) 近似刚度k偏大，即元素偏“硬”。4) 近似的位移偏小，即求得位移的下界。能够以单调趋势逼近于正确解。如曲线.势能总是大于最小状态，故解得上界。给点数近似刚度k偏大，即元素偏“硬”。近似的位移偏小，即求得位移的下界非协调元：在弹性力学中，如板弯曲，相邻元素不仅要求位移本身连续，而且要求位移的导数连续（板弯边界上的相容性）。而在工程上能够保证导数相容的形变往往难以找到，以致工程上只能采用违反相容原则的一些形状函数，由违反相容原则的形函所构成的元素称为非协调元。非协调元性质：、不能以最小位能原理作为它的理论基础。、解的趋势可能收敛，可能不收敛，（取决于网格划分）。对于收敛的趋势也未必满足单调性。可能收敛曲线如图中。QQmjik例： 三角板元： 节点位移系统： 位移形函：W=CL+CLL+CLL+CL+CLL+CLL+CL+ CLL+CLL+CLLL L =( a+ bx + dy )/2a= xy - xy ; b= y- y ; d= - x+ x 。代入节点位移参数：可得：w= N N 变换一下写法：w= H HH H H H H HH H H = L(3 2L) 7LLL H = L ( dL dL ) + ( d d )LLL H = L ( bL bL ) + ( b b) LLL H = 27 LLLi +2=k i j k 下标按循环计算 上述元素只能在结构上做到位移导数连续，在边界上其他点处，位移的法向导数并不连续，这因为：由于法向导数是一个完全的二次多项式，在元素的每条边上，其变化规律位一条二次抛物线，需要三个点上法向导数的相等条件才能维一确定，故相邻两条曲线一般不全重合。故所举三角板弯元为非协调元。例 书P的矩形元，由于坐标的交叉双乘积（不完备），可发现不该是w 或其导数， 都是连续的，这样只要节点的这些参数相同，边界上的这些是没有问题的，但展开N的项，可以发现xy项，或者说缺少了代表热率变形的一项，因此，作为形状函数，是不能保证向正确的解答收敛，因而是非协调元。改进方案之一，是在节点处增加节点参数，并采用完全的埃尔米特三次多项式。1、4、2 非协调元的排先检检验协调元虽可以保证总位能从上往下地正确结果单调收敛，但往往过于复杂，使用麻烦。在工程上往往使用形式简单的非协调元，自然，最小位能原理对此不再适用，那么在什么条件下，这类元素才能导致向正确解收敛呢？ Irons 提出了一个称作“拼片实验”（patch test） 的检验方法。实践表明，这种检验方法是有效的，但“拼片实验”的理论证明尚不清楚。拼片试验内容为“假设由若干元素拼成的一个任意拼片处于等应力状态，这时，其位移函数w( x , y )一般可用一m阶完全多项式函数 P( x , y ) 表示，（入在薄板问题中，m=2） ，而且，在这一拼片的边界上，也设置了符合等应变状态的位移边界条件。然后，将需要检验的某种元素按此条件进行计算。如果最后得到的有限元法解答能和P( x , y ) 一致，那么，称这种元素能够通过拼片试验，而通过拼片试验的元素将给出收敛的结果。注：P( x , y ) 至少应能代表各种等应变状态。如： 拼片： 9节点参数三角变元 *拼片test实 际上成为非协调元的收敛准则（1） 注释： 、两种拼片均有满足各种等应变状态的边界条件。、用非协调元解上述问题。、对比，不一致。、一致，不一致。、说明通过，的网格划分不通过，但误差不大，尚可用。、对于更不规则的网格，其误差可能更大，大到不宜采用。 1、4、3 广义协调元简介1、 概念： 广义协调元利用对变分原理的改造，使边界位移在位能条件降低，从而获得比非协调元有更好性质的单元，即对于不规则网格能通过分片检验；提高计算精度。2、 原理： 用能量和加权线量法联合改造。、势能原理： （b）设单元边界真实变量为 ( x , y ),则在单元边界上每一点若均能满足，= 0 ( c )则单元协调，位移跨越单元连续，若不能保证（c）式成立，则为非协调元，其能量泛函需要做如下修正：) d s ( d )( d ) 式增加了单元不协调位移的能量贡献项。为Lagrange 乘子，其物理意义为单元边界的函数。( d )式可通过带约束的广义变分原理求解。如果用加权残量法，放松单元向位移协调条件，即（）式中T为权函数，以式（c）为约束条件建立单元位移形函的矩阵，可使单元向位移协调条件在某种积分意义下得到满足。、（e）式实际上是一个统一的形式，如取T权函数为函数，即得加权残方程为：0、全Tn为边界应力，则权残方程为：表示常量应力，等价于Irons的分片检验，是非协调远的收敛准则。、具体广义协调元的构造方法可参见：龙驭球著“新型有限元之引论”1.5曲线坐标系及等参元概念1、 5、0 引在工程有限元分析中，有几方面的原因需要需要采用高精度单元：、提高计算精度：、适应不规则的几何边界条件；、用较粗的单元数据准备，获得较高的计算效益；为获得高计算精度采用通常意义下的单元，本身具有某种缺陷：、增加节点参数提高位移参数的阶次，使计算复杂，并增加求解规模、只适用于规则外形的结构，复杂边界难以准确描述。因此，目前多采用等参元在二维问题中得到了极大的成功，以后，在三维问题和板壳计算中，也取得了很好的效果。举例说明对受拉带孔板应力集中的计算结果（元与等参元对比）单元类型 m单元的相对边长应力集中系数K 3.02误差 6810 0.066 0.16 0.008 2.567 2.806 2.860157.15.36等参元 6 0.0662.8695.06m等分0w等参元的优越性：w、精度高（有线性，也有二次等）；、匹配各种复杂的直线及曲线边界；、构造单元方法统一，推广其它类型应用不困难。4.刚度计算、公式；K= = =、Guass消元法消去中间节点，得聚缩刚阵； 、程序系统不再介绍。2总体刚阵的形成及程序设计从有限元的计算过程可知，形成单刚、集合总刚、置边界条件、求等效载荷列阵后，才能得到位移，应变及应力，课程正是按此工作进行，最终达到我们的目的。 2.1总刚的特征及其组装时所需要效能的问题1 总刚特性及投放过程总体结构刚阵的形成是各单元刚阵按节点序排列的集合矩阵（对号入座）。、对称性：依据 Betti互易原理（讲：在结构体的两点1、2作用单位位移位移，（）先1点作用,再2点作用 ( 2 ) 先2点作用,再1点作用 、带状性：（刚阵元素较集中在对角线的位置）一个节点自由度对应总刚的一行（一个节点不可能和结构的所有节点相连），这一行总是和相关的节点刚阵元素相关，这也是总刚的组织进程反映出的特性。、稀疏性：不相关的单元（指节点没有任何联系的单元）没有共同联系的刚性元素，故在对号累加中始终为0。1例： 节点（ 1 ，4）分属两个单元，32 故刚阵中无K系数，同理 （1，6）、（2，6）、（3，4）、（4，6） 4654 )投放过程：*总刚按节点排序（不与单元节点序相关） 考虑的问题：总的按节点排序： 2、投放。2、 考虑的问题： 坐标转换 ；主从节点转换；有效方式及其效率，依据总刚的特点，寻找节约的有序方法（由于对称可仅有一半）存贮效率措施使用存贮单元的多少，用总刚规模替标度量。总刚规模：存贮总刚单元的总数及其存贮的形式。*存贮方式 1、三角存放 （一维存） 总刚规模：KK= 一维存时，要有一个投放算法，使用时，要有一个调用算法。 n m方式：可一维存，可二维存（不好）切齐即可2维存*存贮方式2、 等带宽存贮： KK = m x n ；m 为最大带宽。 按节点编号的最大差计算；由于带状性质 mn （可一维，也可二维） *变带宽存贮： 可将每行第1个非0元素以外的0元素全部排除在存贮单元之外，需要专门的程序计算 KK. 只能一维存，还需要一个每行第1非0元素的管理。 *全希疏存放： 仅存贮各单元的非0元素及分解出的非0元素。需要记每一行的首尾特征，行元素的列位置等，则要有复杂的管理机构。管理机构也消耗单元，效率以远小于其它存方式消耗单元总数为恰当。求解效率： 一般不计编程的效率（程序编的优劣不计），只计存贮效率和求解效率。求解效率当然以机时计算 ，效率成本（即存贮求解的消耗成本）。求解效率与存贮效率要综合评价，但存贮效率一般放在第一位考虑。设：执行的时间为T，存贮量为S，则一个存贮与求解的成本估价是T, S的复杂函数，即： C( S, T ) = T * P( S ) (P 是大于1的多项式)。重要的存贮效率，故上述函数反映了节省存贮资源的重要性，因为减少存贮要求后，也许不再需要动用外存，从而节省了数据传输的费用。T实际上也是S的函数，存入更多的元素，也就需要更多的计算，因此，该意义上讲也要省存贮。3、 影响总刚规模的因素： 除存贮方式的因素之外，另一个影响总刚规模的因素是节点编号的顺序（按带状存贮时，也按全稀疏存贮（编号仅是一个因素），可能影响非0元的分解时的增加策略） 节点编号对带宽的影响：变带宽KK=17总纲形式321 带宽的主要影响因素是节点单元间的编号，差越大，则带宽宽，若差越小，各非0元集中在对角线附近，便会使带宽减少，这一点构成了我们过行节点编号的基本原则，举例说明：最大带宽m = 6KK=36 1 2 3 4 5 6 456m = 4KK=24171 2 3 4 5 6 3216541 2 3 4 5 6 53115m = 3KK=18246由此可见，取定了简化模型之后，对节点如何编号是一个很重要的问题，可直接越小到存贮规模的大小和机时的多少，但如何把节点的编号的顺序便的较优，是一个很复杂的问题，这将在以后专门研究这里给出一个编号的原则（工程准则），即是“优先沿结构较宽的方向按坐标进行编号。 节点单元编号的自动形成所谓节点单元编号就是给XE,ME赋值，对于一般结构的分析问题，都有很多的节点单元，若仅用人力在图纸上划分单元和计算节点坐标，则是很麻烦的，也是易于出错的，因此，这部分工作尽可能地去让计算机完成。自动编序（又称自动划分单元），当然也需要把结构的简化模型取定，然后给某个编序原则，（当然应尽量靠近编排优序的原则，当然也可以依靠算法去优化），这一原则也需要设计者事先在草图上示意性标处，以便检验这原则的实施正确性。（也可以用图显方式对照草图检验），最后即可编制程序。举一例：yx自动说明自动编序实施方法：1591317(13)(19)(15)(21)(14)(20)26101418JA(JB=3)(16)(22)37111519(11)(17)(23)(12)(18)(24)48121620JB (分4等)SUBROU XEME( JA , JB , A , B , XE , ME , NP) NP=(JA + 1) * (JB + 1)DIMENDION XE( NP , 1 ) ME( NP , 3 )D 1 N = 1 NPXE( N , 1 ) = B * ( JB - ( N -1 )/( JA + 1 )XE( N ,2 ) = -A * ( JA + 1 ) * (FLOAT ( N - 1 )/( JA + 1) - ( N -1 )/( JA + 1 )D 2 I = 1 JAD 2 J = 1 JBKI = 2 * ( I + ( J - 1 ) * JA ) 产生偶数单元编号KI = 2 、8、14、20第一行ME ( KI , 1 ) = I + ( J - 1 ) * ( JA + 1 ) + 1ME ( KI , 2 ) = ME( KI , 1 ) - 1ME( KI , 3 ) = ME( KI , 1 ) + JA + 1 ME( KI - 1 , 1 ) = ME ( KI , 1 )ME( KI - 1 , 2 ) = ME ( KI , 3 ) - 1 ME( KI - 1 , 3 ) = ME ( KI , 3 ) RETURN END 2.2 总刚压缩存贮的实现 由上述分析可知，变带宽存贮 KK 最小，显然最省机时，但由于带宽是变的，故只能用一维数组依次存放总刚阵中各行的第1 个非0元素到对角线的各元素，令该数组称为AKM(kk). 要使用AKM(kk)，还必须另外一个整形数组来记忆原总刚阵中K中每一行（对应一个自由度） 的非0 元素到对角线的个数，该数组计为KDKM(NID),NID 为总自由度数。 在实用中，我们并不是真正去记忆每行的非0 元素数，而是记忆每行非0元素的累加数，也即每行最后一个元素在一维数组中的位置，这样给使用带来了很大的方便。所以在机器程序中记忆位置的编号，经常采用后一种方法。按这种方法，例2的KDKM（I）的元素应为1，3，5，9，13，17（KK）。而第3个例题则应为1，3，6，9，12，15（KK）。显然这时KDKM( NID ) = KK (最后一个数)。1 计算原理： P， 点。解释：行快是以每一个节点所对应的自由度数为行数，列号从1到最后一列，2维是 Row = 2 ; 3 维是 Row = 3 .A、 只有与节点之构成的相关点（与节点相关的那些单元的编号），才在总刚 阵中的第i 行块上提供非0子块。B、 第i行块上非0子块的分布情况与节点 i的相关节点编号的差值有关。带宽对于下 阵，则是相关节点最小编号列对角线子块。对上阵存贮，则是相关节点中最大编号列对角线子块（从右到左）。2、框图：书P 。2. 出口刚阵的形成计算由式再由Gauss数值积分，可知几个分块刚阵的计算如下： Kd已在前节介绍过了。于是取代计算为：计算出后，即可按式求解出口刚阵上述是理论求解过程，要求的逆及矩阵乘，矩阵求解是费时的，因此，一般真实计算时，并不沿循理论理论求解过程而是采用高斯消去方法，下面讲述此过程。由式可知，为单元内节点的等效载荷列阵（实际并不指出是那个具体节点，相当于不增加节点。），因形成时仅对单元边界节点进行，故，所以式可写为： 利用矩阵的初等变换，消去阵，使在的第1式中得到： 由于与为对称转置关系，所以消去过程查找元素仅在中去取，十分方便，消去过程如下： 从中的最后一行开始，消