毕业论文-基于有限单元法的苎麻茎秆力学性能分析.doc

武汉科技学院毕业设计（论文）基于有限单元法的苎麻茎秆力学性能分析院系名称 机电工程学院 指导老师 专业班级 机设 054 指导教师职称 副教授 学生姓名 院系毕业设计（论文）工作领导小组组长签字 摘 要苎麻茎秆主要由韧皮纤维、木质部组成。本文重点研究了茎秆、韧皮纤维、木质部三者的力学特性及外载荷作用下的应力应变规律。假定苎麻茎秆是由韧皮纤维、木质部两种异质材料复合而成的复合材料。在获取了平均几何尺寸、平均工程常数的情况下，将苎麻茎秆、韧皮纤维、木质部抽象为空心圆柱体的力学模型。在Comsol软件环境下通过有限单元法建立力学模型，并对苎麻茎秆进行了径向均布静载荷、切向均布静载荷和轴向均布静载荷等不同边界作用下的模拟仿真分析。主要得出了苎麻茎秆在径向、切向和轴向上施加60N500N不同压力下的位移变化、等效应力变化和接触应力变化，找出其应力应变的规律，为进一步求出最佳切断速度提供参考和理论依据，推动苎麻机和机剥麻的开发设计。关键词：:苎麻茎秆; 韧皮纤维; 有限单元法建模; 等效应力应变分析 ABSTRACRamie stalk is composed of bast fibre and xylem mainly. Research on their mechanical properties and their change of stress and strain under different loads are emphasized. Supposing that ramie stalk is made up of two kinds of different material (bast fibre and xylem). On the basis of getting the average size and the average elastic constants. Its assumed that ramie stalk, bast fibre and xylem are the geometrical model of hollow cylinder. Establish the mechanical model through finite element method in the environment of the Comsol Software, the stress and strain distributions of ramie stalk under the transversely even-distributed load, tangential even-distributed load and axial even-distributed load are analyzed through Comsol Software. The major research is that the law of the changeable stress distribution of xylem and bast fibre is analyzed under the loads between 60N to 500N; find out the law of the movement and the stress, In order to further gets the speed of the best cut off to provide the best reference and theoretical basis, promote and drive the design of machine.Key words: ramie stalk; bast fibre ; mechanical model; stress and strain distribution目 录1 绪论11.1 研究目的及意义11.2 植物茎秆力学性能研究现状21.3 研究的内容和目标22 苎麻茎秆力学模型建立42.1苎麻茎秆力学模型建立的有关假定42.2 苎麻茎秆弹性力学分析52.3有限单元数值算法概述92.4本章小结113 苎麻茎秆有限元模拟分析123.1 有限元分析的一般过程123.2线性有限元刚度方程133.3 非线性有限元求解163.4 几何非线性193.5 边界非线性接触问题213.5.1 接触问题的提法233.5.2 接触问题的描述233.5.3 接触界面条件243.5.4 接触约束算法263.5.5 摩擦问题的处理273.6本章小结284 苎麻茎秆静载荷模拟分析304.1 Multiphysics Comsol简介304.2 Multiphysics Comsol 力学分析过程304.2.1.苎麻茎秆结构模型304.2.2 建立力学模型314.2.3 网格的划分及实体生成314.2.4 设置求解域-材料的属性324.2.5 设置边界条件334.2.6 设置其他相关量并求解344.2.7 查看分析结果344.3 本章小结425 结束语43参考文献44致谢47武汉科技学院2009届毕业设计论文1 绪论1.1 研究目的及意义苎麻是我国特有的传统特色的经济作物，其纤维是纺织工业的优质原料。苎麻单纤维较粗而长，是天然纤维中强度最强的纤维，能加工成高支数的纯麻细薄织物，且具有凉爽、散热快和天然的抗菌特性，广泛受到国际消费市场的青睐。我国85%以上的苎麻纺织品用于出口。主要出口欧洲市场、东南亚市场、中东市场、北美市场。尤其是苎麻纱线及其纺织品和棉麻混纺织品，无论是数量上还是价格，在国际市场上均处于垄断地位，没有一个国家能与中国竞争。据海关统计，2005年l 6月份麻类纤维纺织厂物制品进出口总额5.73亿美元。同比增长了12.11%，实现了贸易顺差1.05亿美元。同时，由于欧美对棉花和化纤产品设限种类较多，为适应市场，减少贸易摩擦，棉花和化纤企业将进一步调整产品结构，苎麻纤维作为一种天然生态纤维，在一定程度上会替代部分棉花和化纤产品，必将会对麻类产品产生巨大的需求。近年来，全国苎麻产业呈现良好的发展势头。从全国苎麻市场价格上看，2004年原麻均价由2002年的5000元/吨左右上涨到9000元/吨左右，纯麻纱价格由25000元/吨左右上涨到35000元/吨左右。不到两年时间，苎麻产品市场价格每吨已上涨了近万元，充分显示了我国苎麻纺织业的市场前景。从全国种植面上，2004年全国苎麻种植面积达到10万h耐，总产量23.3万吨，分别较上年增长15%20%。从湖北省苎麻产销形势上看，全省原麻总产量4万多吨，全省苎麻加工企业30多家，年加工能力13万吨以上，现在的苎麻产量还不足需求量的三分之一。同时，随着我国人民生活水平的提高，崇尚自然、绿色消费已成为时尚，人们的衣着向高档化、时尚化、个性化、舒适化方向发展，麻纺服装、室内装饰、卫生保健用品等也将迎来巨大的消费市场。为了获得芝麻的韧皮，需要从麻秆上将皮层(包括麻皮和韧皮)剥下，然后再将麻皮刮去。前者称为剥皮，后者称为刮青。经剥皮和刮青后取得的韧皮经晒干后称为原麻，即是苎麻纺织厂的原料。苎麻一般每年收获3次，头麻的产量高，纤维质量好。剥麻是苎麻生产中的关键环节，也是苎麻产业化发展中的制约因素，剥麻季节时间要求紧，技术要求高，但至今仍普遍采用传统的人工剥麻。这种方式一是用工多，工效低;二是劳动强度大;三是与水稻作物争劳力。471.2 植物茎秆力学性能研究现状植物秸秆是自然界中常见的植物茎干结构。茎的形式分为实心和管状两种，实心茎通常的形式是在轻巧的泡沫型的中心周围围绕着连续的管状外皮，管状的茎是由不同直径的管状圆柱所构成。且芯对该植物秸秆的强度及稳定性具有不可忽略的影响。某些有芯植物秸秆还具有规律分布的节结构，用以增加其自身的抗倒伏、抗冲击等性能。这些植物秸秆结构可以抽象为具有各向异性特征的复合材料，外层包裹具有各向异性特征的柱状内芯及节，且各部分材料特性差异很大。对于植物材料的简单力学参数已有很多相关报道，试验方法多种多样，涉及的植物种类亦是非常繁多。试验方法根据研究目的的不同主要有静力学方法和动力学方法。力变形关系：研究农产品在自然状态下的力变形特性可得到在工程分析和设计中有参考价值的数据。由于原料在运输和加工时一般以自然状态承受各种机械作用，所以开展完整形态下物料力学特性的研究有着重要的理论价值和现实意义。国内有人曾研究了完整形态的苹果在静载作用下的变形规律。应力应变关系：测定物料的应力应变关系曲线时，需将物料制成试件，采用传统的工程材料试验方法，依照应力应变关系提供试验数据。这样类型的试验有同轴拉伸、同轴压缩、剪切以及弯曲等。同轴拉伸的试验结果有：极限抗拉强度的变化。该变化大部分是茎秆含水量和干物质容积密度综合影响的结果，与含水量成反比，与干物质密度成正比。同轴压缩试验表明：茎秆试件的破坏类型取决于植株试件的直线性以及长度和直径比，一般试件长度取为1.97.6cm。当长度和直径之比大时，大部分茎秆是弯曲破坏。1.3 研究的内容和目标本研究结合计算机数字模拟技术。通过对苎麻茎秆、韧皮纤维、木质部进行了拉伸、压缩、弯折试验，分析了三者在不同方向上力学性能的差异，在此基础上建立苎麻茎秆的本构关系和参数的选择确定，再用有限元软件进行模拟分析。具体研究内容:(1)建立苎麻茎秆的本构关系在参考木材、甘蔗物理模型及本构关系建立的基础上，将苎麻茎秆抽象为具有韧皮纤维、木质部结构的空心圆柱体线弹性力学模型，且每一组成部分的材料为具有横观各向同性的性质，确立韧皮纤维、木质部为横观各向同性的本构关系。同时运用复合材料宏观力学分析方法，不考虑韧皮纤维、木质部之间的具体差别，将苎麻茎秆抽象为材料均匀的线弹性横观各向同性的本构关系。从理论上推导了适用于苎麻茎秆、韧皮纤维、木质部以纵向对称轴的横观各向同性材料的柔度矩阵，并得到了各工程弹性常数之间的关系。(2)建立苎麻茎秆的计算模型并进行有限元分析根据已有苎麻茎秆、木质部、韧皮纤维的平均几何尺寸建立几何模型，在Comsol软件中，选择相关的单元和施加约束，对苎麻茎秆两种不同材料(韧皮纤维、木质部)在不同载荷下和不同方向下进行静态分析，分析其在不同载荷下和不同方向下应力应变及位移的变化分布规律。2 苎麻茎秆力学模型建立2.1苎麻茎秆力学模型建立的有关假定(1)几何结构形状的假设苎麻茎秆个体之间及不同部位几何形状存在较大差异，因此，必须对其抽象简化。假定茎秆截面为某一尺寸的圆形截面，壁厚尺寸均匀。整杆为空心直杆，由两种不同的材质构成。(2)材料性能及边界条件假设考虑宏观特性，忽略细观结构的不连续性及其细观缺陷；韧皮与木质部在界面处应变连续，不发生相对位移；各组分材料是均质、线弹性，横观各向同性性质，考虑韧皮与木质部材料粘性特性特征；苎麻整秆宏观均质线弹性，横观各向同性；整秆及组分材料均处于连续变形状态。材料结构对称，即将所研究材料视为结构通直。材料各部分纤维均沿轴向整齐排列。苎麻茎秆边界条件为受横向均布载荷的作用。 图2-1（苎麻茎秆剖面图） 图2-2(苎麻茎秆三维视图)由以上假设，苎麻茎秆计算结构简化为如图2-1和图2-2所示几何物理结构。其中图2-3给出了便于计算的载荷示意图。其载荷为横向均布载荷，载荷密度为P(N/m)。图2-3(苎麻茎秆载荷示意图)2.2 苎麻茎秆弹性力学分析由上面的假设可知，苎麻茎秆的本构关系为横观各向同性，弹性力学的平衡方程、几何方程同样能运用。由此可推导出由位移表示的解析解。由苎麻茎秆横观各向同性材料的本构关系需要考虑如下几个参数：E1 ：茎秆轴向弹性模量；E2 ：茎秆横向弹性模量；12=13：对应茎秆轴向拉伸时的泊松比；23 ：茎秆横向同性面内泊松比；G12=G13 ：横向扭转抗剪模量。为方便推导，有1代替12=13，2代替23。有弹性理论，在不考虑惯性力和体积力的情况下，平衡方程张量表示为，其中i，j=1，2，3。在柱坐标系中形式如下： (2-1)同样，柱坐标系下的几何方程为： (2-2) (2-3)由坐标转换关系，可得由横观各向同性的本构关系的工程常数表达的物理方程为： (2-4) (2-5)将本构方程改写为用应变表达应力的形式如下： (2-6)再将几何方程代入上式得由位移表达的应力如下： (2-7)将位移表达式应力方程代入平衡方程即是所求的用位移表达的平衡方程，其形式如下：(2-7)上式即为由位移函数表示的基本方程组有推导过程中并未考虑材料种类和材料相互作用部位，因此可普遍应用。边界条件：如物体表面上的外力，已知，则边界条件为： (2-8)如物体表面上的位移，为已知，则边界条件为u=u1 v=u2 w=u3。也有混合边界问题，部分表面上已知外力，另一部分表面上已知位移。初始条件：当解动力学问题时，要有下列初始条件：当t=0时间时， (2-9)上式中F1,F2,F3,1,2,3为已知函数。由前面的假设苎麻茎秆为线弹性体，只是材料的本构关系不同，上面推导的偏微分方程组同样适用，给定苎麻茎秆的边界条件即可求出苎麻茎秆在外力作用下的应力、应变及其变化规律。苎麻茎秆抽象为木质部、韧皮纤维两种材料组成，其边界条件包括边界条件和木质部与韧皮纤维相互接触的边界条件。接触条件中假设苎麻茎秆韧皮纤维与木质部在接触面上保持完全接触“即互脱离、也不互相滑动”。这样，在接触面上就有应力和位移两方面的接触条件。应力方面的接触条件是：两弹性体在接触面上的正应力相等，切向位移也相等。对平面问题来说，在通常的边界面上，有两个边界条件。现在看到，在接触面上，有四个接触条件，条件并没有增多或减少，因为接触面是两个弹性体的同样形状的边界。在苎麻茎秆皮、木质部结合处连续条件表达为： (2-10)图2.3中所示边界条件为： (2-11)这样，便得到了所述问题的完整的微分解法。2.3有限单元数值算法概述从弹性力学观点看，在外界因素影响下，物体内产生的应力应变和位移(基本未知量)是某个坐标的连续函数、为无限自由度问题。从物体中任取一点微分单元体加以研究，从而得到平衡微分方程、几何方程和物理方程。这些方程分别体现出了应力与体力、位移与应变、应力与应变之间的相互关系。在给定边界条件下，弹性力学问题也称为偏微分方程的边值问题。可按位移或应力为未知量这两条途径来求解弹性力学问题。从科学的结果和形式上的分析、弹性力学存在解析解和数值解。已经发展的偏微分方程数值分析方法可以分为两大类。一类以有限差分法为代表，其特点是直接求解基本方程和相应定解条件的近似解。一个问题的有限差分法的求解步骤归纳为:首先将求解域划分为网格，然后在网格的结点上用差分方程来近似微分方程。当采用较密的网格，即较多的结点时，近似解的精度可以得到改进。借助于有限差分法，能够求解相当复杂的问题，特别是求解方程建立于固结在空间的坐标系(欧拉(Euler)坐标系)的流体力学问题，有限差分法有自身的优势因此在流体力学领域内，至今仍占支配地位。但是对于固体结构问题，由于方程通常建立于固结在物体上的坐标系(拉格朗日(Lagrange)坐标系)和形状复杂，则采用另一类数值分析方法有限元法则更为适合。有限元法从方法的建立途径方面考虑，它区别于有限差分法，即不是直接从问题的微分方程和相应的定解条件出发，而是从与其等效的积分形式出发。等效积分的一般形式是加权余量法，它适用于普遍的方程形式。利用加权余量法的原理，可以建立多种近似解法。如果原问题的方程具有某些特定的性质，则它的等效积分形式的伽辽金法可以归结为某个泛函的变分。相应的近似解法实际上是求解泛函的驻值问题。里兹法就是属于这一类近似解法。有限元法区别于传统的加权余量法和求解泛函驻值的变分法，该法不是在整个求解域上假设近似函数，而是在各个单元上分片假设近似函数。这样就克服了在全域上假设近似函数所遇到的困难，是近代工程数值分析领域的重大突破。有限元法的特性为:(1)对于复杂几何构形的适应性。由于单元在空间可以是一维、二维或三维的，而且每一种单元可以有不同的的形状。这样一来，工程实际中遇到的非常复杂的结构或构造都有可能离散为由单元组合体表示的有限元模型。(2)对于各种物理问题的可应用性。由于用单元内近似函数分片地表示全求解域的未知场函数，并未限制场院函数所满足的方程形式，也未限制各个单元所对应的方程必须是相同的形式，所以尽管有限元法开始是对线弹性的应力分析问题提出的很快就发展到弹塑性问题、粘弹性问题、动力问题、屈曲问题等。(3)建立于严格理论基础上的可靠性。因为用于建立有限元方程的变分原理或加权余量法在数学上己证明是微分方程和边界条件的等效积分形式。只要原问题的数学模型是正确的，同时用来求解有限元方程的算法是稳定、可靠的，则随着单元数目的增加，即单元尺寸的缩小，或者随着单元自由度数目的增加及插值函数阶次的提高，有限元解的近似程度将不断地被改进。如果单元是满足收敛准则的，则近似解最后收敛于原数学模型的精确解。(4)适合计算机实现的高效性。由于有限元分析的各个步骤可以表达成规范化的矩阵形式，最后导致求解方程可以统一为标准的矩阵代数问题，特别适合计算机的编程和执行。随着计算机软件技术的高速度发展，以及新的数值计算方法的不断出现，大型复杂问题的有限无分析已面成为工程技术领域的常规工作。2.4本章小结(l) 苎麻茎秆的几何尺寸差异悬殊，为了研究分析，将苎麻茎秆、木质部、韧皮纤维分别抽象为平均尺寸的空心圆柱体，并视苎麻茎秆是由木质部、韧皮纤维两种不同材料属性构成的复合材料。(2) 按照复合材料力学理论假定茎秆、木质部、韧皮纤维是均质线弹性，并给出了静态载荷的力学模型，受横向均布载荷的作用。(3) 介绍了物理方程为横观各向同性本构关系，在直角坐标系和圆柱坐标系下所求的通解，并根据力学模型以及苎麻茎秆木质部与韧皮纤维粘附接触，给定了苎麻茎秆受均布横向载荷的边界条件和接触边界条件。3 苎麻茎秆有限元模拟分析3.1 有限元分析的一般过程(1) 有限元法的基本思想有限元分析(FEA，Finite Element Analysis)方法，简称有限元法，是在计算机技术和数值分析方法支持下发展起来的，为解决复杂的工程分析计算问题提供了有效的途径。有限元法的基本思想是将问题的求解域划分为一系列单元，单元之间仅靠节点连接。单元内部点的待求量可由单元节点量通过选定的函数关系插值求得。由于单元形状简单，易于由平衡关系或能量关系建立节点量之间的方程式，然后将各个单元方程“组集”在一起而形成总体代数方程组，计入边界条件后即可对方程组求解。单元划分越细，计算结果就越精确。有限元分析是利用数学近似的方法对真实物理系统(几何和载荷工况)进行模拟，再用简单而又相互作用的元素(即单元)就可以用有限数量的未知量去逼近无限未知量的真实系统。有限元法是在20世纪50年代中期作为结构分析的矩阵法的推广应用到固体力学中的。结构分析的矩阵法是分析含有大量构件的结构系统的分析方法。这些结构系统的构件是简单的拉压杆件以及受弯曲或者扭转的梁，它们在有限个结点上的位移和内力的关系，列出方程组，可以以节点位移或节点内力作为未知数，有时也可以以节点位移和内力混合作为未知数。根据所用的未知数的不同，矩阵分析法可分为位移法、力法或混合法。对离散的结构系统写出方程对于具备结构力学知识的人是熟悉的，采用矩阵符号推导和掌握这些方程也最方便，而对于大型代数方程的求解则可以交给计算机去完成。有限元法也是在电子计算机出现后才真正用于工程分析之中的。由于有了计算机这个强有力的工具，现在有限元法的应用已由弹性力学平面问题扩展到空间问题、板壳问题，由静力平衡问题扩展到稳定性问题、动力问题和波动问题，分析的对象从弹性材料扩展到塑性、粘弹性、粘塑性和复合材料，从固体力学扩展到流体力学、传热学、电磁学等领域。(2) 有限元法的基本要素节点(Node)：:节点是构成有限元系统的基本对象，也就是整个工程系统中的基本点。它包含了坐标位置以及具有物理意义的自由度信息。单元(Element)：单元是由节点与节点相连而成，是构成有限元系统的基础。一个有限元系统必须有至少一个以上的单元。单元和单元之间由各节点相互连接。在具有不同特性的材料和不同的具体结构当中，可选用不同种类的单元，单元中包含物理对象的各种特性。传统的结构矩阵分析中，结构构件的节点力和节点位移之间的关系是精确导出的，而在有限元法中，是根据单元内近似的位移函数导出。单元类型不同，位移函数也不一样。通用有限元分析软件中都要提供多种可供选择的单元，例如梁单元、平面单元、体积单元等。在工程分析时，选择适当的单元可以大大提高计算的效率和精度。自由度(DOF，Degree of Freedom)包括系统自由度和节点自由度。整个系统的自由度，在分析中需要进行适当的约束，系统中每个节点都有各自的节点坐标系和对应的节点自由度，对于不同有单元上的节点，具有不同的自由度。(3) 有限元法解题基本步骤利用有限元法进行结构分析，可采用以下解题步骤：单元剖分和插值函数的确定。根据构件的几何特性、载荷情况及所要求的变形点，建立由各种单元所组成的计算模型。也就是对整个结构进行离散化，再按单元的性质和精度要求，写出表示单元内任意点的位移函数，利用节点位移表示单元体内任意点位移的插值函数。单元特性分析。根据位移插值函数，由应变与位移关系求得单元刚度矩阵。单元组集。把各单元按节点组集成与原结构相似的整体结构，得到整体结构的节点力与节点位移的关系，即整体结构平衡方程组。求解有限元方程。引入支承条件，根据不同的计算方法求解有限元方程，得出各节点的位移。计算单元内的应力和应变。根据结点位移计算相应的节点应力以及单元内部任意点的应力和应变。3.2线性有限元刚度方程从单元到结构考虑到能量泛函极值时，就可得到刚度方程，这里，仅以平面三角元为例来说明建立刚度方程的一般过程。任意平面三角元，如图3-1所示，任一点的位移和结点力分别为图3-1(平面三角单元) (3-1) (3-2)单元结点位移和结点力分别为 (3-3) (3-4)单元能量泛函为 (3-5)式中，称为弹性体内任一点的应变列阵；称为弹性体内任一点的应力列阵。若单元位移函数为则 (3-6)式中，单元的形函数矩阵。该式为位移函数与单元结点位移之间的关系式。 (3-7)式中，单元应变矩阵。该式为应变与单元结点位移之间的关系式。 (3-8)式中，应力矩阵；弹性矩阵。该式为应力与单元结点位移关系。将式(3-7)和(3-8)代入(3-5)式，则有 (3-9)根据势能原理，将上式对取变分，可得 (3-10)由于的任意性，故有 (3-11)从而可知单元结点位移与结点力之间的关系为 (3-12)这即是单元的平衡方程，且(3-13)称为单元刚阵。将单元能量泛函式对结构的所有离散单元予以叠加，则可得到整体能量泛函：将式(3-13)代入上式，可得 (3-14)按照单元刚度阵组装成整体结构刚度阵的”对号入座”规律以及结点力与结点载荷在结点上应力平衡的条件，可以得到 (3-15) (3-16)式中， (S为结点总数)称为整体结构的结点位移； （n=2s）称为整体结构的已知结点载荷。将式（3-15）和（3-16）代入（3-14）式，可得 (3-17)式中，式即为整体结构能量泛函的工程使用显式。根据势能原理，将式对变分，可得式即是整体结构的刚度方程。整体刚阵具有三个特性，即对称性、正定性和稀疏性。虽然刚度方程式是针对平面三角元提出来的，但其思路、做法与形式具有普遍意义，无论对哪一类型单元都可按照上述过程推导出刚度方程，且形式完全相同。3.3 非线性有限元求解有限元是随着计算机的发展而出现的一种基于变分原理来解偏微分方程的数值计算方法。通过把连续体分成有限个节点连接的单元，来解决应力与应变的分布。由于可独立分析每一个单元体，从而可解决初始条件分布不均的问题，适用于以大应变和非线性的材料行为为特征的塑性加工问题。有限元可分为两大类：线弹性有限元法和非线性有限元法。其中线弹性有限元法是非线性有限元法的基础，二者不但在分析方法和研究步骤上有相似之处，而且后者常常要引用前者的某些结果。线弹性有限元以理想的弹性体作为研究对象，所考虑到变形建立在小变形的基础上。这类问题材料的应力应变呈线性关系，满足广义胡克定律，应变与位移也是线性关系。线弹性有限元问题是求解线性方程组的问题，只需要较少的计算时间。非线性有限元法所求解的非线性问题主要有：材料非线性：即材料的应力与应变是非线性关系，在工程中普遍的材料非线性问题的非线弹性、弹塑性、粘塑性和蠕变等；几何非线性：是由于位移之间存在非线性关系引起的，这类问题包括大位移大应变问题和大位移小应变问题；非线性边界(接触非线性)：由于接触和摩擦的作用，接触边界属于高度非线性边界。非线性问题求解主要表现在：(1)非线性问题的方程式非线性的，一般需要迭代求解；(2)非线性问题不能采用叠加原理；(3)非线性问题不总有一致的解，尽管几何、材料及接触等均正确，有时也甚至没有解。这三方面的因素使非线性问题的求解过程更加复杂、计算时间更长。因此非线性有限元程序不仅需要做复杂的算式和有效的数据管理，而且必须包含合理的逻辑来指导求解过程。非线性的求解流程见图3-2。工程中所有的问题都是非线性的，为适应工程问题的需要，在解决某些具体问题时，往往忽略一些次要因素，将它们近似地作为线性问题处理，这也是完全合理的。但并不是一切问题都可以简化为线性问题。非线性问题的分类：线性弹性力学基本方程的判别：(1) 表征材料应力应变关系的本构方程式线性的。(2) 描述应变和位移之间关系的几何方程是线性的。(3) 建立于变形前状态的有限元平衡方程是线性的。(4) 结构的边界条件是线性的。工程实际问题中，上述4条往往不能同时满足，通常习惯上把不满足上述条件(1)的称为材料非线性；条件(2)、(3)不能满足时称为几何非线性；(4)不满足时称为边界非线性。对于苎麻茎秆的有限元分析，既存在几何非线性，也存在边界非线性问题；接下来针对这两个问题进行阐述。读取输入数据、空间分配、数据检查计算等效节点载荷矢量组集矩阵矩阵求解应力计算收敛结果输出网格自适应下一个增量步停止施加载荷增量是否迭代循环图3-2(非线性求解流程图)3.4 几何非线性几何非线性是有结构变形的大位移所造成。我们在讨论线性弹性力学问题时，均隐含一个假设：结构在外载作用下产生的位移及应变都是很小的，在建立结构或微元体的平衡条件时，可以不考虑物体位置和形态的变化，也就是用变形前的状态建立平衡条件。同时还认为应变与变形之间存在线性关系。因此，在线性有限元计算中，仍假定结构加载过程中单元的几何形态之间存在形态基本不变。这实质上是一种线性近似，它的近似包含了两个方面：一是应变与位移之间作了线性化处理，而忽略高阶应变的小量，即=B，其中B为线性应变矩阵；二是把平衡方程的坐标系建立在平衡前初始坐标上，即将结构变形后平衡状态用变形前初始结构平衡状态不作任何修正地加以描述，这就是所谓小变形假设的近似处理。几何非线性关注的问题就是将上述两个近似处理恢复事物的本来面目，对上述第一个近似线性化处理撤消，即将B矩阵由线性矩阵转变为包含高阶微量的非线性矩阵B ,这就是习惯上所谓几何小变形非线性问题，如板壳横向弯曲的挠度问题，它可以表述为结构在加载过程中不能忽略小应变的有限转动的弹性力学问题。另一类几何非线性问题是指有限变形问题，如橡胶等高分子结构材料，它们即使在弹性状态下，也可能产生很大的变形和位移，它们的变形过程已经不可能直接用初始状态加以描述，且平衡状态的几何位置还是未知的。同时解决这了问题的另一个难度还在于它们的应力、应变定义和度量准则与线性问题有所变化，因此必须给出新的定义。由此给几何大变形非线性问题的方程建立和求解带来复杂性。之前所讨论的问题都是基于小变形的假设。它包含了两方面的内容：一是假定物体发生的位移远小于物体自身的几何尺度，在此前提下，建立结构或微元体的平衡条件时可以不考虑物体的位置和形状的变化，因此分析中不必区分变形前和变形后的位形；二是假设在加载和变形过程中的应用可用一阶微量的线性应变进行度量，即应变与位移成一阶线性关系。通常将满足小变形假设的即小变形、小转动，且应变不能说转动的高阶微量的问题称为几何线性问题。工程中碰到很多不符合小变形假设的几何非线性问题，它们主要包括小应变有限转动等有限变形问题以及小应变、小转动，但是应变与转动相比是高阶微量的问题。我们顺次称它们小变形几何非线性问题、有限变形几何非线性问题和结构稳定性的许多初始屈曲问题。小变形几何非线性有限元方程的建立和求解之前介绍过的建立有限元单元平衡方程时采用应变能泛函求极值，单元体的应变能等于外力所做的功与单元变形所产生的应变能之和。如果用F表示内力和外力矢量的总和，那么由位移增量表示的应变能泛函可以写成如下表达式： (3-18)式中F代表所有载荷列阵；d为位移列阵；d为应变列阵，为应力列阵。如果用应变的矢量形式写成位移增量和应变增量的关系，有 (3-19)将上式代入式（3-18），并对应变能泛函求极值，得到非线性问题的一般平衡方程式 (3-20)上式中积分运算时应用通常的方法，由各个元素的积分对于节点平衡所做的奉献综合而成。式(9-2)的推导，仅应用了最小势能原理，没有涉及材料、几何时线性还是非线性问题，因此不论位移是大还是小，都应完全适用。在线性有限元中，由于应变与位移之间及应力与应变之间的关系均是线性的，即=B =D将上述关系代入式(3-20)，因此可以得到 (3-21)现在考察小变形几何非线性问题，应变和位移的关系是非线性的，因此矩阵B是的函数，为简便起见，可以写成 (3-22)式中B0作为线性应变分析的矩阵项；BL是大位移应变矩阵项，并取决于，它是由非线性变形引起的。一般地说，BL是位移列阵的线性函数。而应力应变关系还仍然认为是线性弹性关系，则有 (3-23)式中D为材料弹性矩阵；0为初应变列阵；0为初应力列阵。很显然，式(3-22)的解可以通过迭代方法求得。如果采用牛顿拉斐逊(NewtonRaphson)方法，必须寻求d和d之间的关系。通过式(3-20)，取得微分，于是有 (3-24)如不考虑初应变和初应力的影响，得到 (3-25)由式(3-22)可知 (3-26)所以式(3-24)可以改写为 (3-27)式中，且K0表示通常的小位移动线性刚度矩阵，即，表示大位移所引起的影响矩阵，一般称大位移矩阵或初位移矩阵，而式(3-24)右端第一项可表述为，它是应力水平的对称矩阵，称它为初应力矩阵或几何刚度矩阵。这样，完成了小应变几何非线性问题的非线性方程组的建立。3.5 边界非线性接触问题边界非线性包括两个结构物的接触边界随加载和变形而改变引起的接触非线性(它又包含有摩擦接触和无摩擦接触)，也包括非线性弹性地基的非线性边界条件和可动边界问题等。两个物体相互接触后，随着两个物体间接触合力的变化，他们之间的接触面大小、接触处的应力均会发生变化。这些变化不仅与接触面力的大小有关，而且与两个物体的各自材料性质有关。即使材料性质是线性弹性的，接触问题仍然表现出强非线性性质。如果材料性质是非线性的，接触非线性性质表现更为强烈与复杂。非线性弹性地基的边界和可动边界非线性是十分明显的。此外，碰撞问题是一类与边界质点速度有关的边界非线性问题。边界非线性中有相当一部分问题往往不再遵循最小势能原理，而呈现耗散特性。例如考察摩擦边界和碰撞问题。此外上述问题相互耦合的非线性问题更为复杂。在非线性问题求解中，它们有几个共同的特点。对于一般非线性方程或方程组，到目前为止，尚未找到一种理论上能精确求解的方法，现在均采用近似解法，其中数值解法是近似解法之一，也是采用最多，应用最广的一种。数值近似解法具有以下特征：(1) 非线性问题的解不一定是唯一的。(2) 解的收敛性事前不一定能得到保证，还可能出现不稳定状态，如振荡现 象，甚至于发散。(3) 非线性问题的求解过程比线性问题更为复杂和困难，结果的处理也更为复杂。接触问题在工程中处处可见。例如齿轮啮合、滚珠轴承、螺栓连接，以及在金属冲压成型过程中薄板与模具之间的接触、摩擦以及滑动过程均是接触非线性问题。早在18世纪的经典弹性力学问题中已开始研究最简单的赫兹(Hertz)接触和库伦(Coulomb)摩擦。但由于其高度非线性性质，使得研究一直难以深入。随着计算机技术的发展，带动了计算数学和有限元技术的发展，使得这类问题具备了实际求解的可能。但是，与材料非线性和几何非线性相比较，有限元在求解接触问题和摩擦问题方面还是相对滞后的。弹性接触问题的描述和求解仍属有待于进一步研究的领域。到目前为止，还没有一致公认的具体解法。求解接触问题的研究过程可分为三个相互重叠的阶段。在第一阶段中只涉及刚体，用牛顿定律和欧拉定律建立运动微分方程，摩擦用库伦定律描述。碰撞则用牛顿恢复系数来处理。第二阶段是以赫兹的开创性论文为标志。首先讨论弹性体的静态接触问题，描述如何计算接触区域和接触力。第三阶段是以数值计算为本，借助于计算机以及接触过程的近似处理，人们能对复杂的集合形体和对复杂的材料本构关系其接触过程进行数值模拟，主要是使用有限元。3.5.1 接触问题的提法弹性力学的基本方程包括平衡方程、几何方程、本构方程以及边值或初值条件。材料非线性和几何非线性分别由本构方程的非线性性质而引起；而接触非线性问题，它是有边界条件的非线性性质引起。这主要表现在两个方面：一是接触表面的改变，即有、自由表面的一部分转变为接触边界，或者反之，有接触边界放松接触面而蜕变成为自由边界；二是接触面的变形、摩擦和滑移，可能表现出强烈的非线性性质，同时滑动时的摩擦力也会表现出非线性性质。物体相互接触，随着载荷的增加或减少，在接触面会出现弹塑性变形及塑性卸载情况。凡此种种，不但使问题成为高度非线性，而且使其过程成为不可逆。从而导致接触问题与摩擦问题耦合的高度非线性，这使其求解过程更为困难和复杂。对于这类问题，若要求得解析结果是十分困难，甚至是不可能的，只能采用数值模拟方法。其中非线性有限元方法是目前可以给出这类问题近似解的最有效途径之一。3.5.2 接触问题的描述接触可分为法向接触和切向接触。在切向接触中摩擦十分重要且又异常复杂，它又分为静摩擦和动摩擦。由于非穿透条件和压应力条件，接触问题中含有不等式约束。处理弹性接触问题最根本的方法是罚函数法和拉格朗日乘子法。接触问题属非线性问题，通常只能用数值方法，例如牛顿-苪弗逊迭代求解。在非线性有限元中，为了确定切线刚度矩阵和接触余量，描述接触问题的非线性方程组需要线性化。在接触过程的计算中，根本性的问题是如何满足非穿透性条件。如果不考虑物体自身的相互接触，两个物体间的相互接触可用下图形象化表示。图3-3(两个物体间的相互接触)物体i(i=1,2)分别具有内部区域i和边界i。在体力和面力的作用下，物体在空间中运动。在任一时刻每个物体的边界由三部分组成，即 (3-28)其中u表示已给定位移的边界，q表示已给定应力的边界，c表示接触边界。非穿透性条件是指在任一时刻任一位移物体的内点不可能同时属于另一个物体，即区域1和2的交集为空集：12= (3-29)这一非穿透性条件虽然在物理上很好理解，但在数学上却难以处理，对于边界条件则有 (3-30)由于物体的位移和应变影响着接触区域的位置和大小，即使对几何线性问题和线弹性材料，接触条件也导致非线性问题。在求解非线性的接触问题中，通常需要使用迭代法进行迭代运算，直至接触区域和物体的位移场收敛为止。非穿透性条件亦可借助于一个定义在边界上的距离函数gN来描述，即gN0。此外在接触区域上还有单向应力条件，法向应力必须为压应力，即pN0。综合非穿透性条件和压应力条件，则得描述连续体法向接触的互补问题 (3-31)上式中还表示只有接触的区域，即gN0，才能传递作用力。对于带有摩擦的接触问题也有类似的互补条件。接触区域不能承受拉应力，因而导致单面约束，以使对给定接触区域的接触问题也成为非线性问题。3.5.3 接触界面条件对于接触问题，除了其场变量需满足固体力学基本方程、给定的边界条件以及动力问题的初始条件外，还要满足接触面上的接触条件，主要为不可侵切条件和摩擦条件。对于接触或将要接触的两个物体，其界面接触状态可分为分离、粘结接触和滑动接触三种。对于这三种情况，接触界面的位移和力条件是各不相同的，正式由于实际的接触状态在这三种情况中的转化，导致了接触问题的高度非线性特点。按照Hughes给出的方法建立如图所示的局部坐标系1，2，3，其单位基矢量分别为e1、e2和e3=n，其中n为物体A在接触点处的表面的单位外法线矢量。(1) 分离状态：位移条件 (3-32)面力条件 (3-33) (3-34)(2)粘结接触状态： (3-35) (i=1,2) (3-36)面力条件 (3-37) (3-38) (3-39)(3)滑动接触状态：位移条件面力条件 (库仑摩擦理论)以上各式中，u为位移矢量，pi、gi分别为三个局部坐标方向的接触面力和间隙量，上标A、B分别表示物体A和B，d30为接触点对初始间距，为摩擦因数。接触类型判断准则表4-1类型准 则tntn+1分离分离g30接触g30粘结粘结分离滑动滑动粘结分离滑动3.5.4 接触约束算法在将物体离散为有限元之后，接触点对待间隙量可以用相应节点位移和初始间隙来表示， (3-40)对于一般的点面接触情况，利用单元位移模式也可得到上述形式，这在很多文献中有叙述。对乘子也在单元一级作离散和近似，在考虑了整体和局部坐标系的转换关系后，可得 (3-41)其中，为相关节点处接触内力列阵。这样，离散化后弱形式中的接触约束项的一阶变分可写为 (3-42)其中， (3-43)这样，对于有限元离散系统可以得到相应于各种接触约束算法的代数方程组Lagrange乘子法 (3-44)罚方法： (3-45)增广Lagrange乘子法 (3-46)摄动Lagrange乘子法 (3-47)在实际的接触问题的分析过程中，往往需采用增量迭代法进行求解，在每一迭代步中需检测接触对的接触状态，并将相应的界面条件引入到系统方程组中。对于动接触问题，还需考虑惯性力和阻尼力，弱形式须以修正的Hamilton原理表示，由于前面所给出的界面条件也适用于动接触问题，故在有限元列式中，需加入相应的惯性力和阻尼力项(3-48)。3.5.5 摩擦问题的处理通常采用的摩擦模型主要有库仑摩擦模型和剪切摩擦模型。(1)库仑摩擦 库仑摩擦模型为： (3-49) 节点合力的形式为： (3-50) (3-51)