ANSYS结构有限元高级分析方法与范例应用（第三版） - 万水书苑-出版资源网

ANSYS 结构分析的理论背景及功能概述

本章导读 有限元分析软件不同于一般的应用软件，它对使用者的知识水平提出了较高要求。有限 元软件的使用者应当具备必要的理论基础，这样才能对计算过程中的各种参数作出正确的设 置，并对程序的计算结果作出正确的判断和评价。 ANSYS 结构分析软件功能全面，涵盖线性静力分析、特征值屈曲分析、模态分析、谐响 应分析、瞬态分析、响应谱分析、随机振动分析、子结构分析、子模型分析、热传导计算等分 析类型，可以模拟材料非线性、几何非线性、接触分析等各类非线性力学行为，还提供了复合 材料分析、断裂力学分析、热—结构耦合分析、声场以及声—结构耦合分析、优化设计、概率 设计以及并行求解计算等技术，形成了完善的高级结构分析解决方案。本章将围绕 ANSYS 计 算程序应用的各个方面介绍与之相关的理论背景知识。 本章包括如下主题：  有限单元法的基本思想与解题过程  ANSYS Mechanical 的分析能力及一般分析流程  ANSYS 结构静、动力分析的理论基础  ANSYS 热传导分析的理论基础  ANSYS 流固耦合分析的理论基础  ANSYS 结构非线性分析的基本概念和算法

1.1 有限单元法的基本思想与解题过程

1.1.1 起源于杆系结构分析的有限单元法

ANSYS 软件的理论基础是有限单元法（Finite Element Method）。有限单元法是结构工程 师和应用数学研究人员共同智慧的结晶，其基本概念起源于杆系结构的矩阵分析法。1956 年，

Turner 和 Clough 等人首次将刚架分析的矩阵位移法应用于飞机结构的分析中。之后的 1960 年，Clough 将这种处理问题的思路推广到求解弹性力学的平面应力问题，给出了三角形单元 求解弹性平面应力问题的正确解答，并且首次提出“有限单元法”的名称。之后，应用数学家 和力学家们则通过研究找到有限元方法的数学基础——变分原理，进而将这一方法推广应用于 求解各种数学物理问题，如热传导、流体力学、电磁场以及各种耦合场问题。有限单元法目前 已经发展成为一种多物理场通用的数值计算方法。 为了介绍有限单元法的基本概念，首先来回顾一下杆系结构的矩阵分析方法。杆系结构 的矩阵分析方法以杆端位移作为基本未知量，因此又被称为矩阵位移法。该方法的计算过程可 以概括为以下几个步骤： （1）结构离散化。 将整个结构离散为若干个杆件单元，通常情况下一个结构杆件被作为一个单元，也可根 据杆件截面变化以及支座和加载情况等将一个构件划分为多个单元。 （2）单元特性分析。 单元分析是有限元分析的重要一环，也是结构整体分析的基础。单元分析的任务是得到 单元节点力与节点位移之间的关系方程——单元刚度方程，此方程的系数矩阵被称为单元刚度 矩阵，单元刚度矩阵中的元素称为刚度系数，刚度系数的数值等于发生单位位移所引起的力。 对弹性结构而言，刚度系数的实质为弹簧的刚度系数，即单位变形量引起的弹性回复力。为说 明刚度系数的意义，下面给出几个简单的例子。 图 1-1 为只承受轴向力的等截面直杆，长度为 l，截面积为 A，弹性模量为 E。直杆的右 端发生单位微小轴向位移时，由材料力学可知所需施加的力为 F=EA/l，因此该直杆的轴向抗 拉（压）刚度系数为 EA/l，其实质就是一个刚度系数为 EA/l 的轴向弹簧。 图 1-1 杆的刚度系数 图 1-2 为一等截面的直梁，跨度为 l，材料弹性模量为 E，截面惯性矩为 I，梁左端固定， 右端受到竖向支撑的约束。现在如果强迫梁的右端发生单位微小支座位移时，竖向支杆中的约 束反力为_{F3EI l/} 3_{，因此此梁的右端抗侧向变形刚度为3EI l 。/} 3 图 1-2 梁端的抗侧移刚度

图 1-3 中，等截面梁的左端发生单位微小的转角时，需要施加的力矩为M 3EI l/ ，因此 该梁左端的平面内转动刚度为 3EI/l。 图 1-3 梁端的抗转动刚度 在单元分析时，单元上作用的各种分布载荷需要按照静力等效原则移置到杆件的节点上， 形成等效节点载荷。在矩阵位移分析中，由于各杆件放置角度不同，还需将单元刚度方程（即 单元刚度矩阵和单元等效节点力）转换至结构整体坐标系中。 （3）结构分析。 结构分析则是通过结构整体的平衡关系和变形协调条件将单元刚度矩阵按照节点位移编 号和对号入座的方式组合形成结构的总体刚度矩阵，同时各杆件上作用载荷形成等效节点载荷 并转换至整体坐标系下，形成结构的节点载荷向量与节点位移向量之间的关系。 （4）引入边界条件。 未引入边界约束条件的总体刚度矩阵是一个奇异矩阵，通常采用划行划列降阶法、改 0 置 1 法、对角元素充大数法等方法引入边界条件。如结构的约束足够，不论是静定还是超静定 系统，引入边界条件后的总体刚度方程均可解。 （5）求解线性方程组得到节点位移。 求解结构总体刚度方程组，即可得到杆端位移（节点位移）。 （6）计算杆端力。 根据节点位移计算等效节点载荷引起的杆端力，叠加固端约束杆端力，得到实际结构各 杆件的杆端力。 上述过程中最为核心的两个环节是单元特性分析和结构分析。下面以图 1-4 所示的简单桁 架结构的矩阵分析为例，向读者介绍相关的一系列基本概念。 图 1-4 桁架结构有限元分析的原理

如图 1-4 所示是由两根杆件①和②组成的平面桁架，杆件的截面积均为 A，材料的弹性模 量都是 E，长度分别为 L1和 L2。图中，Xi和 Yi分别为桁架结构在第 i 个节点处所受到的水平 以及竖向外力，ui和 vi分别为结构中节点 i 的水平和竖向节点位移，u 和_ij v 分别为单元_ij j 杆端 节点 i 的水平和竖向杆端位移（上标 j 表示单元号，节点编号采用结构中的节点编号）， j i X 和 j i Y 分别为单元j 的杆端节点 i 处的水平和竖向杆端力。 由于杆件在节点处均为铰支，因此杆端无力矩作用，每个节点的力和位移均有两个分量， 或者说，一个节点具有 2 个自由度。我们取以上与坐标方向平行的正交分量分析杆端节点力和 节点位移之间的关系。 对于杆件①，设其轴线方向（由节点 1 指向节点 2）与 x 轴正向成角，则其杆端力与杆 端位移之间的关系可以表示为： 1 2 2 1 1 1 1 2 2 1 1 1 1 2 2 1 1 2 2 1 2 2 1 2 2

cos cos sin cos cos sin cos sin sin cos sin sin

cos cos sin cos cos sin cos sin sin cos sin sin

X u Y _EA v L X u Y v                                        _{ }     _{ }        _{ }                _  _      其中，杆端力与杆端位移之间的系数矩阵称为单元刚度矩阵，我们用 kij来表示刚度矩阵 中位于第 i 行 j 列的元素。下面，我们以位于第 1 列的各元素 ki1（i=1，2，3，4）为例说明刚 度矩阵元素的计算方法和力学意义。 首先，令单元①节点力与节点位移关系式中的 1 1 1 u  ，v₁1v1₂ u₂1  ，则得到：0 X₁1k₁₁， 1 1 21 Y k ， 1 2 31 X k ， 1 2 41 Y k ，这表明：当单元①的节点 1 发生单位微小水平位移（ 1 1 u =1）且 同时约束单元①的其他所有节点位移（ 1 1 1 1 2 2 0 v v u  ）时，单元①在各节点处受到杆端力， 这时的杆端力组成一个平衡力系，它们在数值上（以与坐标轴正向一致为正值）就等于刚度矩 阵中的元素，表征单元①抵抗节点 1 水平位移 1 1 u 的刚度。这些力的数值很容易由材料力学求 得，如图 1-5 所示。 图 1-5 刚度系数的意义图示 当 1 1 1 u  ，其余节点位移分量v₁1v₂1 u₂1 时，杆单元①的缩短量为0  l cos ，于是杆

件受到的轴向压力为EAcos/L₁，这就是杆件①在节点 1 处受到的节点力，其在 x 方向和 y 方向的分量分别为： 2 11 cos / 1 k EA  L ，k₂₁EAcos sin /L₁ 杆件①在节点 2 处受到的力，其大小等于杆件的轴力，且方向与节点 1 处的节点力相反， 其在 x 方向和 y 方向的分量分别为： 2 31 cos / 1 k  EA  L，k₄₁ EAcos sin /L₁ 继续对各位移分量做类似的分析，即可得到刚度矩阵的全部元素，完成桁架单元①特性 的分析。实际上，通过这一操作，我们已经得到了平面桁架杆单元的一般刚度特性，对于桁架 单元②，其相当于倾角为 90 度的特例，我们可直接得到单元刚度方程为： 2 2 2 2 2 2 2 2 2 2 2 3 3 2 2 3 3 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 1 X u Y _EA v L X u Y v                                         _    _{ }      单元分析结束后，就需要进行结构分析，将单元的刚度特性集成得到结构的刚度方程， 根据几何条件——相邻单元公共节点位移的协调关系，即： 1 1 1 2 1 2 2 2 1 1 1 1 2 2 2 2 2 2 3 3 3 3 u u，v v；u u u；v v v；u u，v v 又根据节点处力的平衡条件，作用于某节点上的外力应等于包含该节点的各单元的杆端 节点力的合力，即： 1 1 1 2 1 2 2 2 1 1 1 1 2 2 2 2 2 2 3 3 0 3 3 X X，Y Y；X X X ，Y Y Y ；X X  ，Y Y 以上 3 个节点处沿坐标轴方向的力的平衡条件就是结构的节点力与节点位移的关系，代 入单元杆端力与杆端位移的关系（两个单元的刚度方程），即得到总体刚度矩阵的形式。 显然，总体刚度矩阵是一个 6×6 的矩阵。由于节点 2 同时连在杆件①和②上，因此两根 杆件将共同抵抗公共节点 2 的变位，这在结构总体刚度矩阵中表现为两个单元对应刚度系数的 叠加，即： 1 1 1 11 12 1 1 2 2 2 2 22 22 23 2 2 2 3 32 33 3               _{ }                      _{ }  1 1 F K K 0 Δ F K K K K Δ F 0 K K Δ 其中F_i{X Y_i, }_i T，Δ_i { , }u v_{i i} T，K 表示单元m_ij m 的节点 i 和节点 j 之间的刚度关系（单 元刚度矩阵的子块），即节点 i 的单位位移引起节点 j 处的杆端力。 这样我们就完成了结构分析，得到总体刚度方程，其系数矩阵就是总体刚度矩阵。在引 入边界约束条件之前，总体刚度矩阵是奇异的。通过前述方法在总体刚度方程组中引入约束条 件，求解线性方程组即可得到节点位移，进而得到杆端力等结果。具体过程在相关的教程中都 有介绍，此处不再展开叙述。 回顾上述杆系结构矩阵分析过程，其基本思路可以概括为：每一个单元的力学特性被看 作是构成整体结构的砖瓦，总体结构作为一系列单元组成的集合体，通过单元特性的组合装配

即可提供总体结构的力学特性。尽管本节是以平面桁架结构为例介绍杆系结构矩阵分析的基本 实现过程，但其中的一些基本概念同样适用于其他杆系结构系统或其他各种弹性结构系统的有 限元分析。理解这些概念对读者正确应用 ANSYS 软件进行结构建模和分析是十分必要的。 1.1.2 有限单元法处理问题的一般过程及算法特点 在工程计算中，很多问题的解决都与求解微分方程有关。比如材料力学中梁的弯曲问题 可以归结为求解以挠度为未知函数的常微分方程的问题，弹性体的变形问题可以归结为求解以 位移为未知场变量的偏微分方程的问题，固体的热传导问题可以归结为求解温度为未知场变量 的偏微分方程（热传导方程）的问题。由于各种因素，这些问题中可以通过解析方法来求解的 仅仅是其中很少的一部分，大量实际工程问题只能借助数值方法来求解。另一方面，由于连续 的求解域包含无限多个微元体，这与数值求解的实现过程发生矛盾。有限元方法通过将连续的 求解域离散化的方法来解决这一矛盾，该方法的核心思想可以概括为：首先将系统离散为单元， 再将单元特性组合得到系统的特性。通过求解域的离散，待分析的系统成为具有有限个自由度 的元素（或单元）的集合体，数学物理微分方程问题便可简化为适合于数值求解的结构型问题。 有限元方法求解一般物理问题的基本实现过程可以概括为以下几个环节： （1）创建求解域的几何模型。 创建几何模型通常是花费时间较多的环节，该环节的目标是创建一个与求解域形状相一 致的几何模型。可以通过 3D 的 CAD 设计软件创建几何模型之后导入有限元分析软件中，也 可以直接使用 ANSYS 等有限元软件的前处理（建模）程序，如 DM、SCDM 等来创建。创建 几何模型时，可以根据问题特点或划分网格的需要对实际结构进行各种合理的简化，如去除表 面突起、去除小孔面或倒角面等。 （2）几何模型离散化（Meshing）。 几何模型离散化是形成有限元分析模型最为关键的一个环节。在这个环节中，将连续的 求解域离散化为网格状的分块区域的集合体，因此离散化过程又被形象地称为划分网格（即 Meshing）。网格化后每一个分块所包围的区域称为一个单元（Element），整个离散的结构模型 由有限个单元组成，被称为有限元模型（Finite Element Model）。相邻的单元彼此之间只是在 一些特定的点处连接在一起，这些点称为节点（Node）。每个节点处的未知的场变量被称为自 由度（Degree of Freedom），是有限元方程求解的基本未知量。不同学科类型的单元，其节点 具有不同的自由度。对于热分析，节点自由度为温度。对于结构分析，自由度即节点所具有的 各种位移自由度，如平移自由度或转动自由度等。对其他物理场的分析也都有相应的场变量作 为节点自由度。

ANSYS 软件提供了功能强大的网格划分程序 ANSYS Mesh，基于 ANSYS Mesh 可以对各 种线、面、体等几何对象进行网格划分。对于一维问题中线段划分为多个小段的过程也习惯地 称为 Meshing。尽管不像二维或三维问题那样形成网格状的单元组合体，但这些划分后的线段 仍被称为网格。在 ANSYS 或其他程序中，Meshing 实际上代表的是一个过程，即借助于几何 模型形成有限单元组合体的过程。 单元是组成有限元分析模型的基本元件。各种工程结构中的杆、梁、索、板、壳、膜以及 三维块体等都可离散化为有限数量单元的组合体，即有限元分析模型。组成结构模型的单元可 以各式各样，但是各种单元又有其作为结构分析单元的共性。随着有限元方法的不断发展，结

构分析的单元类型日渐丰富，每一类的单元都有很多不同的算法，可用于模拟处于各种不同形 态以及不同工作状态下的结构构件。表 1-1 列出了与结构分析相关的常见单元类型及应用场合。 表 1-1 结构分析单元的类型 单元类型 单元形状 单元应用场合 梁、杆、索 直线段、曲线段 用于模拟二维尺度远远小于第三维的结构构件，如 梁、桁架杆、拉索、一般弹簧、阻尼器等 板、壳、模 三角形、四边形 可以直边或曲边 用于模拟二维尺度远远大于第三维的结构构件，模 只承受面内的力，板只承受面外的力，壳则受到面 内外载荷的共同作用 面单元 三角形、四边形 可以直边或曲边 用于模拟平面应力、平面应变、轴对称断面等可简 化为二维问题的弹性体变形问题 三维体单元 四面体、六面体 单元的各边可以是直线也可以是曲线 用于模拟一般的不能简化的三维弹性体变形问题 连接单元 线状的单元 用于模拟两点之间的特殊连接，如各种非线性的弹 簧、阻尼器等 注：表中的梁、板类型的单元常被称为结构单元，而平面单元、体单元则被称为实体单元。 （3）单元分析。 按照分块近似的思路，选择一个相对简单的函数来近似地描述 Mesh 后每一个单元内部的 未知场变量的分布规律。单元内部场变量的插值可以选用线性插值函数，也可以通过二次以上 的插值函数，此处的插值函数在有限单元法中一般被称为形函数。形函数的选取将会直接影响 到相关计算结果的精度。在图 1-6 中，求解域是一维的，被分成若干个线性单元，每一个单元 均采用线性形函数对场变量 u(x)在单元的内部进行插值。 x u ( x ) 节点 单元 未知的场变量函数 节点处场变量（未知） 图 1-6 单元内场变量的分段插值 如图 1-7 所示的三角形平面单元，其形函数通常采用线性的面积坐标 L。单元内及单元边

界上的位移是线性分布且由节点位移和形函数唯一确定的。由于两个相邻单元公共节点位移相 同，因此公共边上的位移也相等，即相邻单元在公共边界上位移是相协调的。 对于在边上包括中间节点的单元，通常采用高次多项式作为插值函数（形函数），如边上 有一个中间节点的单元采用二次的形函数。 x, u y, v 1 (x1, y1) (u1, v1) 2 (x2, y2) (u2, v2) 3 (x3, y3) (u3, v3) A 图 1-7 2-D 问题的三角形单元 图 1-8 是一些不同形状和阶次的平面以及实体单元，其中上面一排为一些直边的 2-D 以及 3-D 线性单元，下面一排是对应形状的边上有中间节点的二次单元。单元的阶次由插值函数来 确定。 图 1-8 各种线性单元及高阶单元 确定了插值函数后，通过在单元内应用变分原理或加权余量法即可建立单元特性方程。 （4）建立离散系统的有限元方程。 这一环节实际上就是结构的总体分析。根据相邻单元在公共节点处场变量的连续性条件 将单元特性方程集合为离散系统的总体有限元方程。集合的过程可概括为根据节点编号顺序对 号入座的过程。离散系统的总体有限元方程是奇异的线性代数方程组。 （5）边界条件的引入及求解代数方程组。 向离散系统的总体有限元方程组中引入边界条件，求解线性代数方程组，得到未知的场 变量在节点处的数值。在结构分析中，边界条件主要是位移约束。在热传导分析中，边界条件 包括恒温边界、对流边界、辐射边界等。 （6）计算导出物理量并进行可视化处理。 根据计算得到的节点处的场变量值及单元内部的插值模式得到场变量在单元内部的分 布，并通过基本解导出其他物理量。在结构分析中，可以由位移计算结果导出应变以及应力等

导出量；在热分析中，可以通过温度计算结果导出热通量。计算完成后，利用可视化的后处理 程序绘制各种变量的等值线图、路径分布图等。ANSYS Mechanical 包含了专门的计算结果后 处理程序，可以显示各种物理量的分布图形、变形过程动画，还可用于绘制变量之间的曲线关 系图，这些直观的后处理功能可以帮助分析人员对计算结果进行快速有效的检查和评价。 由上述基本处理过程可见，有限单元法具有如下特点：  思路直观，便于理解和应用。有限单元法通过单元特性集合形成结构特性的思路很直 观，分析人员可以通过直观的物理概念来理解和应用此方法。  具备复杂的几何适应性。有限单元法中通过划分单元分片进行自由度的插值，而单元 可以是各种不同形状，3D 问题可以是四面体、金字塔、六面体、三棱柱等各种形状， 因此该方法对复杂的求解域几何形状具有良好的适应性。  具备严格的理论基础。对有限单元法的理解可以建立在不同的层次上，可以通过直观 的力学概念，也可以建立在严格的理论基础之上。建立单元特性方程时采用的变分原 理或加权余量法在数学上已经被证明是微分方程的等效积分形式。  适用范围广泛。有限单元法中场变量分片插值，并未限制场变量所满足的方程形式， 因此这一算法的适用范围十分广泛，除了弹性结构力学分析，还可应用于温度场、流 场、电磁场等几乎所有物理场问题和耦合场问题的求解。  适合计算机处理。有限单元法的基本求解过程采用规范化的矩阵形式表达，最后导出 的求解方程可以统一为标准的矩阵代数计算问题，所以此方法特别适合于计算机编程 处理和执行。 正是因为上述特点，有限元方法在经历了几十年的研究和发展后，逐步成为解决各学科 数值计算问题的一种普遍性方法，被广泛地应用于结构、流体、电磁、热传导、声场等多学科 的稳态分析、瞬态分析、特征值分析以及与之相关的耦合场问题和非线性问题的分析中。有限 元方法已经成为工程计算中一种不可或缺的重要分析手段，其计算结果被作为各类工程结构设 计和工业产品性能评估的重要参考依据。

1.2 ANSYS Mechanical 的分析能力及一般分析流程

本节对 ANSYS Mechanical 的求解能力及一般分析流程进行简要的介绍。

ANSYS Mechanical 是 ANSYS 软件系列中的结构力学分析模块，该模块是基于有限单元 法的一个广义的计算力学模块，该模块可以计算的问题类型包括：  静力分析：计算结构在静力作用下的位移、应力。  特征值屈曲分析：得到结构在给定载荷作用下的失稳模式。  模态分析：得到结构在一系列特定频率下的固有振动模式。  谐响应分析：得到结构在一系列给定频率谐载荷作用下的最大稳态动力响应。  瞬态动力学分析：计算结构在任意动力载荷作用下的时间历程响应。  响应谱分析：基于模态组合方式计算结构对动力载荷谱（如地震、风载、波浪力）的 响应。  PSD 分析：基于模态组合方式计算结构对随机载荷的响应。  热传导分析：计算系统的稳态以及瞬态温度场，可以处理三种主要的热传递方式：热

传导、热对流和热辐射，支持复杂热载荷和边界条件的施加，支持非线性热特性（如 考虑与温度相关的导热系数、对流换热系数、接触传热计算等）。  声场分析：计算声波在介质中的传播问题。  耦合场分析：能处理热－结构耦合分析、声场－结构耦合、流－固耦合以及压电分析、 热－电分析、热－电－固分析等耦合问题。  转子动力学分析：转子系统的临界转速、不平衡响应及稳定性分析。  子模型分析：在关心的区域细化网格得到更加精确的应力解的结构分析技术。  子结构分析：一组单元静力凝聚为一个超单元，利用超单元减小计算自由度的结构分 析技术。子结构技术还可以应用于动力分析，即部件模态综合技术。  结构优化分析：包括参数优化和形状优化技术。  断裂力学分析：计算应力强度因子、J 积分、裂纹尖端能量释放率。 上述问题中，静力分析、瞬态分析（完全法）中可以考虑各种非线性因素。ANSYS Mechanical 具有全面的结构非线性分析功能，可以处理下列非线性力学行为：  材料非线性：可以分析非线性弹性、率无关金属塑性模型（双线性、多线性以及非线 性的等强、随动及混合硬化模型）、D-P 塑性模型、超弹性、粘弹性、粘塑性、蠕变、 混凝土材料、铸铁、垫片、层状复合材料及层间失效、形状记忆合金。支持各向异性、 与温度相关的材料特性及材料曲线拟合。  几何非线性：可以计算任意的大应变、大位移（转动）、应力刚化效应。对于大应变 问题，提供网格重划分技术以克服收敛困难问题。  接触非线性：ANSYS 具有高效的接触算法，可以分析任意的点－点、点－面、线－ 线、线－面、面－面接触问题，提供高效快捷的接触自动定义和识别技术，能处理各 种自接触、密封流体渗漏等复杂问题。此外，还提供了螺栓预紧、点焊、Joint 等部 件间连接方式。  单元非线性：施工过程模拟（材料的添加和去除），非线性连接单元（弹簧、阻尼器）、 单向受力构件（如索、仅受压的杆件）等问题。  混合非线性：以上各种非线性问题的组合，如大变形问题伴随材料塑性问题，大变形 问题伴随接触问题，或大变形、塑性和接触同时存在的问题。有几何缺陷以及力学缺 陷的受压结构的非线性屈曲分析就是一种典型的混合非线性静力分析，通常是在大变 形的同时伴随塑性。 此外，Mechanical 在温度场分析中也可以包含各种非线性因素，如与温度相关的导热系数、 对流换热系数、辐射边界条件、接触传热、摩擦生热、相变潜热等问题。 尽管 ANSYS 分析能力十分强大，可以分析的问题类型众多，但是无论进行上述哪种具体 问题的分析时，其分析流程都是类似的，即都包含前处理、加载求解、后处理三个分析环节。 下面对基于 ANSYS 有限元分析的实现过程进行简要介绍。 1．前处理：创建分析模型 前处理阶段的目标是建立分析对象的有限元模型。 在开始建模之前，要根据所分析结构的实际情况以及关心的主要问题对建模的过程进行 必要的规划，对建模过程的细节和要达到的深度做到心中有数。 建立分析模型通常是首先建立几何模型，然后对其进行 Mesh 得到分析模型。建模时根据

结构特点选用适合的单元类型，分析模型宜根据结构受力特点作必要的简化，如可以消除不必 要的细节特征、刚硬的局部通过约束方程设置为刚性体、附属子结构简化为弹簧连接等。 2．加载及计算 在前处理得到有限元模型的基础上，根据问题特点选择合适的分析类型。对于缓慢加载 的问题可进行静力计算，对于载荷作用引起的加速度不能忽略的问题则需要进行动力分析，其 中简谐载荷作用的问题可进行谐响应分析，一般动态载荷作用的问题可进行瞬态分析，随机载 荷作用可选择 PSD 谱分析。对于每一种分析类型，还需要进行相关的选项设置，如载荷步设 置、非线性设置、瞬态设置等。 另一方面，所施加的边界条件和载荷要反映结构的实际受力状态，如果简支梁采用了固 定约束，肯定无法得到正确的解答。 在选择了分析类型，施加了载荷并正确设置了分析选项后，即可进行求解。 3．后处理 计算完成后，通过后处理程序对结果进行列表、图形显示、直观动画显示等，对结果进 行分析，验证结果的正确性。

1.3 ANSYS 结构计算的理论基础

1.3.1 ANSYS 静力分析的理论基础 本节从虚功原理出发，以一般弹性结构分析的三维 8 节点线性单元（ANSYS 的 SOLID185 单元）为例介绍弹性结构静力有限元分析的理论基础。 进行单元分析。 连续的三维实体结构被离散化为由有限个单元所组成的离散系统后，对其中任意一个单 元 e，单元内任意一点的位移、应变、应力分别用向量形式表出： { } { , , }u  u v w T { } {  _xx,_yy,_zz,_xy,_yz,_zx}T { } {  _xx,_yy,_zz,_xy,_yz,_zx}T 由弹性理论可知，单元应变与位移之间应满足几何关系： { } [ ]{ }  L u 其中，[ ]L 为三维问题的微分算子阵，其具体形式为： / 0 0 0 / 0 0 0 / [ ] / / 0 0 / / / 0 / x y z L y x z y z x                    _{   }     _{   }    _{   }   

单元应力与应变之间应满足下列物理关系： { } [ ]{ }  D  其中，三维问题的弹性矩阵 [ ]D 为： 1 0 0 0 1 1 1 0 0 0 1 1 0 0 0 (1 ) [ ] 1 2 sym 0 0 (1 )(1 2 ) 2(1 ) 1 2 0 2(1 ) 1 2 2(1 ) E D                   _{ }                   _  _   _{ }     _             _    其中， E、 分别为弹性模量和泊松比。 进行单元分析。 在离散系统的任一单元 e 的体积 Ve上应用虚功原理： * * { } { }d { } { } e T e T e V   V  u F



其中， * {ue}为虚拟的变形状态下的节点位移向量， * { } 为相应于 * {ue}的单元内虚应变， {Fe}为实际状态下的单元 e 的节点载荷向量（在单元分析中视节点力为外力），{ } 为实际状 态下的单元应力。 单元内任意点的位移（虚位移）可通过节点位移（节点虚位移）值和插值函数得到，即： 8 1 8 1 8 1 ( , , ) { } ( , , ) [ ]{ } ( , , ) i i i e i i i i i i N u u u v N v N u w N w                     _{ }   _{ }     _{ }_{ }                  







其中，N （i=1，…，8）为各节点的形函数。形函数是位置坐标的函数，需要满足两个_i 基本条件：一是各节点对应的形函数在本节点值为 1，在其他节点值为 0；二是各节点的形函 数之和为 1。前一个条件说明形函数具有很直观的意义，它表示当对应节点发生某个单位位移 分量，其余节点都不动时，整个单元内的位移分布情况。后一个条件则表示当整个单元发生刚 体位移时，单元内各点位移均相等（等于各节点的位移）。[ ]N 为形函数矩阵，对8 节点单元， 其具体形式可写为：



1 3 3 8 3 3



[ ]N  N I[ ]_ ,  ,N I[ ]_

其中，[ ]I _{3 3} 为三阶单位矩阵。 用形函数表示的单元位移代入几何关系，得到应变与节点位移向量之间的关系为： { } [ ][ ]{ } [ ]{ }  L N ue  B ue 其中，[ ]B 为应变矩阵，其形式为： 1 8 [ ] [[B  B],  ,[B ]] 其中各分块为： 3 3 / 0 0 0 / 0 0 0 / [ ] [ ][ [ ] ] / / 0 0 / / / 0 / i i i i i i i i i i i N x N y N z B L N I N y N x N z N y N z N x                      _{   }     _{   }    _{   }   _ 虚位移引起的应变也可由应变矩阵和节点虚位移表示： * * { } [ ]{  B ue} 上式代回虚功方程，可得到 * * { } [ ] [ ][ ]d { } { } { } e T e T e e T e V u

_

B D B V u  u F 两边消去节点虚位移，得到 [ ] [ ][ ]d { } { } e T e e V B D B V u  F



上式可写为简洁的形式： [Ke]{ } {ue  Fe} 式中给出了单元节点力与节点位移之间的关系，通常被称为单元刚度方程。其中，[Ke]称 为单元刚度矩阵，其表达式为： [ ] [ ] [ ][ ]d e T e V K 

_

B D B V 单元刚度方程右端的节点载荷向量{Fe}由如下几部分组成： {Fe} { Fef} { FeS} { FeC} { Fei} 其中，{Fef}为体积力{ }f 的等效节点力向量，其表达式为： { } [ ] { }d e ef T V F 

_

N f V {FeS}为表面力的等效节点载荷向量，其表达式为： { } [ ] { }d e eS T S F 

_

N T S 其中，Se_{表示单元受到表面力 T 作用的表面域。}

{FeC}为单元集中节点载荷向量，{Fei}为相邻单元对单元 e 的作用力，在后续的结构分 析时{Fei}作为单元之间的内力相互抵消。 实际计算中，ANSYS 采用等参变换技术，单元刚度矩阵及载荷向量均采用数值积分方法 计算。图 1-9 为 8 节点单元的等参变换示意图。 4(-1, 1, -1) (1, -1, 1)6 (1, -1, -1)2 1 7 5 ₈ 6 4 2 0 z y x 3 0 fsz fsy fsx 8(-1, 1, 1) 7 (1, 1, 1) (-1, -1, 1)5 (-1, -1, -1)1 3(1, 1, -1)    图 1-9 三维单元的等参变换 在自然坐标( , ,   )中，可以很方便地给出形函数的表达式： 1 (1 )(1 )(1 ) 8 i i i i N        （i=1，…，8） 所谓等参元，是指单元内各点的坐标通过与位移相同的插值函数（形函数）表示： 8 1 8 1 8 1 ( , , ) ( , , ) ( , , ) i i i i i i i i i N x x y N y z N z                     _{ }   _{ }                           







对自然坐标的导数，可通过链式微分法则表示为对总体坐标的导数的表达式： [ ] i _{i i} i i i i i i N x y z N N x x N x y z N N J y y N x y z N N z z                      _  _{    }         _{ }    _  _{    }  _    _{ }                                                  _  _{   }_{ } _ _{ } _   上式的系数矩阵[J]被称为雅克比矩阵，用自然坐标的形函数和总体坐标值表示[J]如下：

8 8 8 1 1 1 8 8 8 1 1 1 8 8 8 1 1 1 [ ] i i i i i i i i i i i i i i i i i i i i i i i i i i i N N N x y z N N N J x y z N N N x y z                                                         



















于是，上节中的应变矩阵[ ]B 中各元素对总体坐标的偏导数可以通过对自然坐标的偏导数 来表示，即： 1 [ ] i i i i i i N N x N N J y N N z              _  _{ }                    _  _          _{ } 积分体积元素按下式也进行等参变换： dV J d d d   于是，单元刚度矩阵由下式计算： 1 1 1 1 1 1 [ ] [ ] [ ][ ]d [ ] [ ][ ] d d d e T e T V K B D B V B D B J       

_



_{  }

这样，单元刚度矩阵成为关于自然坐标的标准区间内的积分。类似地，各种等效载荷向 量也可变换为这一形式的积分。 ANSYS 在实际计算中采用相关的数值积分方法（如高斯积分等）进行计算，这样单元刚 度矩阵的元素以及与单元有关的量均可通过标准数值积分来计算如下： 1 1 1 1 1 1 1 1 1 ( , , )d d d ( , , ) l m n i j k i j k i j k I    f       H H H f          

_{  }





其中，H_i、H_j、H_k为数值积分的权系数，对于 2×2×2 积分点方案，各积分点位置的自 然坐标分别取±0.577，各点权重均为 1.000。以刚度矩阵的计算为例，仅需要计算[ ]B 和J 等 在积分点处的值。 由于等参元的坐标和位移采用相同的插值函数，因此应用等参元还有一个好处，就是对 曲线边界问题的处理。在这类问题中，如果单元采用精度较高的二次形函数，则曲线边界也通 过节点坐标的二次插值来逼近，可以显著降低线性单元折线代替曲线的误差。 在单元分析的基础上，下面进行结构总体分析。 结构分析需要基于两个基本条件：节点的平衡条件、相邻单元在公共节点上的位移（自 由度）协调条件。 如图 1-10 所示为由三个刚度相等的弹簧并联而组成的弹簧系统。左端固定，要使该弹簧 系统右端发生一个公共的变形量



，所需施加的外力就必须等于系统三个弹簧弹性回复力之 和，即：

( ) 3 F KKK   K 这一并联弹簧的模型实际上已经包含了结构分析的基本思想，即各个弹性元件的刚度系 数都作为弹性系统刚度系数的一部分，由各个元件来共同抵抗外力的作用且变形保持协调（即 几何条件），各个弹簧元件的内力之和等于作用外力（即平衡条件）。 图 1-10 并联弹簧模型刚度系数计算示意图 由于相邻单元在公共节点上的位移是协调的，因此结构分析中，相邻单元在公共节点对 应位移上的刚度矩阵元素被叠加到一起，以共同抵抗公共的节点位移。这体现了上述变形协调 的原则以及节点的平衡条件。实际计算处理中，结构整体分析的过程可以形象地概括为“对号 入座”。如果与一个节点对应的单元刚度矩阵元素在总体刚度矩阵中被作为一个子块，则结构 总体刚度矩阵一共包含节点总数个子块。只要按照节点编号的顺序将各子块放到结构总体矩阵 的相应位置，即可由单元刚度矩阵组集形成总体刚度矩阵。与之相对应的，所有单元的节点载 荷向量也要按照结构中的节点编号次序送入结构载荷向量对应位置，相邻单元在公共节点的相 互作用力相互抵消，外力作用的等效载荷及节点集中载荷进行叠加，这样就形成了结构的总体 刚度方程。 [ e]{ } { ef} { eS} { eC} e e e e K U  F  F  F









其中，{ }U 为结构节点位移向量，各求和符号表示矩阵或向量元素在对应自由度上的叠加， 右端三项依次为单元等效体力、表面力载荷向量与单元节点集中载荷向量的叠加。令： [ ] [ e] e K 



K { } { ef} { eS} { eC} e e e F 



F 



F 



F 则结构总体刚度方程可简写为： [ ]{ } { }K U  F 其中，[ ]K 为总体刚度矩阵，{ }F 为总体节点载荷向量。 通过结构分析，得到的结构总体刚度方程是奇异的，外在表现为结构的各元件（单元） 连接成为一个整体后，虽然元件之间的相对位置不再变化，但仍然存在由于整体约束不足而引 起的整体刚性位移（总体平移或者转动）的任意性。 在总体刚度方程中引入位移边界条件是在求解之前必须进行的步骤。通过边界约束条件 的施加可以排除结构发生整体刚性位移的任意性，使得在一定载荷作用下的结构位移响应可以 唯一地确定。注意，必须给结构加上符合实际情况的正确的约束条件，否则会导致错误的分析 K K K F

结果。 这一步在数值实现过程中表现为对结构总体刚度矩阵的元素进行数值处理，常见处理方 式有直接代入法、改 0 置 1 法、对角元素乘大数法等。  直接代入法：将已知的节点位移（包括位移约束和强迫支座位移）相对应的自由度消 去，形成以其他剩余自由度为未知量的线性方程组。  改 0 置 1 法：仅适用于给定的零位移。当支座位移为零时，可以将总体刚度矩阵中对 应于零位移自由度相应行列的主对角元素改为 0，而主对角元素置为 1，经过这样的 处理后，在不改变原来矩阵阶数的情况下，解方程的时候可以使相应自由度为零。  对角元素乘大数法：适用于任意给定的位移数值。对于已知的节点位移，把与已知位 移自由度对应的主对角元素乘以一个较大的数 α（可取为 1010_{量级），并且在对应的} 方程右端项用改变后的主对角元素和已知位移的乘积。采用这种处理方法时方程的阶 数不变，节点位移顺序不变，编制程序十分方便。 引入了边界约束条件后的离散化结构方程就具有唯一的解，ANSYS Mechanical 提供了一 系列直接求解器和迭代求解器计算自由节点的位移，再通过总体刚度方程计算结构的支反力。 在得到节点位移向量后，单元内部任意点的位移则可通过解出的节点自由度插值得到。用各单 元的应变矩阵乘以节点位移向量即可得到单元各积分点处的应变，根据应力－应变关系（线性 分析中就是胡可定律）由积分点的应变计算积分点处的应力。对于线性分析，ANSYS 会将积 分点的应力值外插到节点，并对相邻单元公共节点的应力作平均处理后作为节点量输出。 需要指出的是，由于应变是位移的导数，由位移解导出的应变和应力解的精度将低于位 移解。对于应力梯度较大的部位（应力集中部位）的求解，需要采用适当的网格密度或采用高 阶单元来提高应力解的精度。ANSYS 提供了功能强大的结果可视化处理程序（后处理器）， 可直观地显示模型中各种物理量的分布情况。显示的量可以是单元解，也可以是节点解。如果 显示的单元应力解等值线图有明显的跳跃性间断，则表明需要加密网格以提高应力解的精度， 节点解由于是相邻单元平均化的结果，不会出现不连续现象。 以上即为 ANSYS 静力结构分析相关的理论背景。 1.3.2 ANSYS 动力分析的理论基础 对于结构动力学问题，根据达朗贝尔原理，在平衡方程中需要考虑惯性力和阻尼力。惯 性力作为体积力，即： { } {f_I  { }}u 式中的为单元的密度，各点的加速度用节点加速度值按形函数插值，即： { } [ ]{ }u  N ue 代入体积力等效节点载荷表达式，得到： { I} [ ] { { }}d [ ] [ ]{ }d e e ef T T e V V F 

_

N  u V 

_

 N N u V 采用类似的处理方式，阻尼力也作为体积力的一部分： {f_D} { c u{ }} 式中的 c 为单元阻尼系数，各点速度用节点速度值按形函数插值，即： { } [ ]{ }u  N ue

代入体积力等效节点载荷表达式，得到： { D} [ ] { { }}d [ ] [ ]{ }d e e ef T T e V V F 

_

N c u V 

_

c N N u V 上述惯性力、阻尼力作为体积力代入单元刚度方程，可得到如下单元动力学方程： [Me]{ } [ue  Ce]{ } [ue  Ke]{ } {ue  Fe} 其中，[Me]、[Ce]分别称为单元的（一致）质量矩阵和单元阻尼矩阵，其表达式为： [ ] [ ] [ ]d e e T V M 

_

 N N V [ ] [ ] [ ]d e e T V C 

_

c N N V 单元分析之后同样是结构分析，根据节点的平衡条件及相邻单元公共节点的变形协调性， 各单元矩阵元素送入总体矩阵对应自由度的位置（对号入座），单元等效节点载荷送入结构载 荷向量对应自由度位置，注意到节点载荷与时间相关，于是形成下列结构动力方程： [M]{ } [ ]{ } [ ]{ } { ( )}u  C u  K u  F t 其中： [ ] [ e] e M 



M [ ] [ e] e C 



C [ ] [ e] e K 



K { ( )} { e( )} e F t 



F t 上式为结构动力有限元分析的一般方程，其中[ ]M 、 [ ]C 和 [ ]K 分别为结构的总体质量矩 阵、总体阻尼矩阵、总体刚度矩阵，

_{ }

u 、

 

u 和

 

u 分别为节点的加速度向量、速度向量、 位移向量，{ ( )}F t 为结构节点载荷向量。式中带有下标 e 的各求和符号的意义为对所有单元进 行对应元素的叠加。 对于模态分析，仅考虑结构自身的特性，与外部作用无关，通常也不考虑阻尼，此时结 构的动力方程简化为： [M]{ } [ ]{ }u  K u 0 如果令：

 

u 

 

_i cos(_it) 代入结构自由振动有限元方程，简化得到： 2 {[ ]K _i[M]}{ }_i 0 上式为一个齐次线性方程组，其有非零解的条件为： 2 det ([ ]K _i[M])0 上式是结构频率特征值分析的基本方程，通过求解这一特征值问题可得到结构的各阶自 振频率和振型。对振型计算结果，ANSYS 程序还提供了两种归一化方法：一种方法是振型向 量最大分量归一，其他各分量成比例缩放；另一种是关于质量矩阵归一化，即满足：

{ } [_i T M]{ } 1_i  对于考虑应力刚化效应的模态分析，只需在以上频率特征值方程的刚度矩阵中增加应力 刚度项即可，依然为特征值问题。 作为结构动力分析中一类常见的特殊问题，当结构外载荷为简谐载荷时，ANSYS 提供了 谐响应分析来给出系统在简谐载荷作用下的最大稳态响应。假设外载荷的频率为 Ω，外载荷和 稳态位移响应的相位分别为 ψ 和 φ，简谐外载荷及稳态位移响应分别为：



( )





max





maxcos maxsin





1 2



i i t i t i t

F t  F e e  F  iF  e  F iF e



( )





max





maxcos maxsin





1 2



i i t i t i t u t  u e e  F iF  e  u iu e 将上两式代入结构动力有限元方程，可得： 2 1 2 1 2 (  [M] i [ ] [ ]){C  K u iu } { F iF} 求解此方程组即可求出给定加载频率 Ω 的稳态位移响应幅值和相位角。 对于一般的瞬态结构动力问题，上述结构有限元分析的一般方程实际上是一个常微分方 程组，需要引入初始条件和位移边界条件后再进行求解。瞬态问题的求解方法可分为振型叠加 法和逐步积分法两大类。振型叠加法的思路是：利用振型矩阵作为变换矩阵，将多自由度系统 原本相互耦合的振动方程组转化为等数量解耦的单自由度振动方程并分别求解，以求得的单自 由度解作为系数将结构的各阶模态进行叠加并求和，最终得出结构的瞬态响应。逐步积分法的 思路是：将原本在任意时刻都需要满足的运动方程的位移代之以只要在离散的时间点满足动力 学方程；而在一定时间间隔内，对位移、速度和加速的关系采取某种假设，这样就可由初始条 件逐步求出后续各个时间点的响应值。时域数值积分方面，ANSYS 提供了 Newmark 方法、 HHT 方法等，此处不再展开讨论。

1.4 ANSYS 热传导计算的理论基础

固体结构的热传导是有限元技术的重要应用方面，本节介绍与热传导问题相关的理论 基础。 对于一般的各向异性固体，其瞬态热传导问题的基本提法如下： 瞬态热传导方程： 0 x y z v T T T T T T c k k k Q t x x y y z z     _  _ _  _ _  _               ( , , )x y z   其中，c 为比热，k 、_x k 、_y k 分别为三个方向的导热系数，_z Q 是单位体积的热功率。_v 三类边界条件（恒温边界S 、已知热流量边界₁ S 、自然对流边界₂ S ）₃ ： 1 ( , , ) TT  x y z S 2 0 ( , , ) x x y y z z s T T T k n k n k n q x y z S x y z             0 3 ( ) 0 ( , , ) x x y y z z T T T k n k n k n h T T x y z S x y z             

其中，q 为通过边界_s S 的热通量，₂ T 为环境温度。₀ 温度场初始条件： 0 ( , , ) ( , , ) t T _ T x y z  x y z   采用有限元方法分析温度场时，对于某一个特定单元 e，其内部的温度场通过节点温度插 值得到： 1 ( , , , ) ( , , ) ( ) [ ]{ } n e i i i T x y z t N x y z T t N T  



 其中，[ ] [N  N₁... N_n]，为形函数矩阵，{Te} { ( ) ... T t₁ T t_n( )}T，为单元 e 的节点温度向量。 取近似函数 ( , , , )T x y z t 在边界S1上强制满足边界条件，在单元 e 上应用 Galekin 方法，其 域内的余量为： d e ei i x y z v V T T T T T T R N c k k k Q V t x x y y z z                  _  _{ } _{ } _{ } _               



或写为矩阵形式： { } [ ] d e T e x y z v V T T T T T T R N c k k k Q V t x x y y z z                  _  _{ } _{ } _{ } _               



对上式中的中间几项进行分部积分，得到： { } [ ] d [ ] [ ] [ ] d [ ] [ ] [ ] d [ ] d e e e e T T T T e x x y y z z V S T T T T x y z v V V T T T T R N c V N k n N k n N k n S t x y z N T N T N T k k k V N Q V x x y y z z          _   _  _    _                      











令： [ ] [ ] [ ]d [ ] [ ] [ ] [ ] [ ] [ ] [ ] d { } [ ] d e e e e T V T T T e x y z V e T Q v V C c N N V N N N N N N K k k k V x x y y z z P N Q V                        







对相关项进行改写得： { } [ ] { } [ ] { } { } [ ] [ ] [ ] d e e e e e e T T T e Q x x y y z z S T T T R C T K T P N k n N k n N k n S x y z          _   _     





_

上式的最后一项，单元 e 的表面积分项应包含各类热边界以及单元交界面，相邻单元交界 面积分项在总体求和时相互抵消，强制温度边界上不考虑此项，热流边界及对流边界条件代入 边界项中（暂不考虑辐射边界），整理即可得到单元 e 的离散导热方程，如下： [ ] {Ce Te} ([ ] K e[H] ){e Te} { P_Q}e{P_qe} { P_he} 其中：

3 2 3 0 [ ] [ ] [ ]d { } [ ] d { } [ ] d e T S e T q s S e T h S H h N N S P N q S P h N T S    







在各单元上求和： [ ] {e e} ([ ]e [ ] ){e e} ({ _Q}e { _qe} { _he}) e e e C T  K  H T  P  P  P









令： [ ] [ ] [ ] ([ ] [ ] ) { } ({ } { } { }) e e e e e e e e Q q h e C C K K H P P P P      







则有： [ ]{ } [ ]{ } [ ]C T  K T  P 上式就是离散形式的系统热传导方程，如考虑非线性因素，则改写为如下的一般形式： [ ( )]{ } [ ( )]{ } { ( , )}C T T  K T T  P T t 其中，t 是时间，{T}是节点的温度向量，[C]是系统的比热矩阵，[K]是系统的热传导矩阵， {P}是节点热载荷向量。 对于稳态问题，不考虑与时间相关的项，简化为： [ ( )]{ } { ( )}K T T  P T 以上就是固体热传导问题的基本方程及有限元方法的基本处理过程。

1.5 ANSYS 流固耦合分析的理论基础

ANSYS Mechanical 可与 ANSYS Fluent 或 ANSYS CFX 专业流体分析模块组合使用，以 完成实时双向的流固耦合动力过程分析。该分析的难点在于流场的压力作为流场中工程结构所 受的动载荷，而由于结构的振动，作为流场求解域边界的固体壁面是可动的。结构振动引起的 流场动边界对流场造成影响，流场中的压力分布的变化反作用于固体结构的表面。这类问题在 工程中经常遇到，比如充液容器的晃动、海工结构在波浪作用下的振动、水坝地震响应、桥梁 结构的风致振动、管道振动、涡轮叶片的振动等。 在流固耦合计算中，固体域的处理按照标准的位移有限元方法进行域的离散。下面简单 介绍一下流体域和流固耦合界面处理相关的基本理论。 由于液体和低速流动的空气等常见工程流体一般可视作不可压缩流体，因此本节仅介绍 不可压缩粘性流体域的基本控制方程。 0 v w u y z x         

2 2 2 2 2 2 2 2 2 2 2 2 d d d d                _   _  _   _               _   _ x y u p u u u f t x x y z v p v v v f t y x y z 2 2 2 2 2 2 d d         _   _ z w p w w w f t z x y z 对于不可压缩的粘性流体，上述以速度作为基本未知量的动量 N-S 方程，结合连续性方 程，原则上可以求解三个速度分量以及压力等共计四个未知量，这些方程构成不可压缩纯流动 问题的封闭性。 对于湍流问题，在高 Re 数情况下，其最小湍动尺度仍然远远大于分子平均自由程，因此 流体依然被视作连续介质。理论上来讲，N-S 方程也是描述湍流流场的基本方程，即湍流场中 任意位置的速度、压强、密度等瞬时值都必须满足该方程。但是，由于湍流在空间和时间上的 变化很快，基于 N-S 方程的直接数值模拟的计算量十分巨大，因此目前这种方法仅仅用于湍 流理论研究领域中。在工程湍流计算中，为减少计算量，发展了一系列实用的湍流计算模型和 算法。在数值计算中常用的湍流模型有 k-epsilon、k-omega、Reynolds Stress、DES、LES 等， 相关模型的参数请参考流体分析软件的理论手册，这里不再展开。 流固耦合界面需要满足运动学和动力学条件。 运动学条件： 流固耦合界面的流体质点与固体质点的法向速度保持连续，即满足： fn f f s f s s sn v v n v n  v n v 动力学条件： 流固耦合界面上法向应力保持连续，即满足： ij sjn ijnfj ij sjn     其中_ij和_ij分别为固体应力分量和流体应力分量。

1.6 ANSYS 结构非线性分析的基本概念与算法

在以上各节所讨论的全部结构问题中，我们采用了线性分析的基本假定，即：  材料的本构关系（应力－应变关系）是线性的，即满足 Hooke 定律。  应变与位移关系（几何方程）是线性的。 满足上述假定的问题才被看作线性问题进行分析。 在工程实践中，有很多问题符合或者近似地符合上述假定，因此结构线性分析在工程问 题的处理中得到了广泛的应用，但是也有不少问题并不符合上述假定，如果仍作线性分析，就 会造成很大的误差甚至得到不能接受的分析结果，这些情况下就涉及到非线性分析。 本节介绍非线性问题的一般概念分类以及非线性问题的一般求解方法。与 ANSYS 非线性 分析相关的求解选项及设置方法将在本书第 13 章中进行系统介绍。下面，首先通过一个简单 例子向读者说明三类非线性问题的基本概念。

（1）如图 1-11 所示的简单直杆拉伸问题中，若应变与杆端位移之间的关系是非线性的， 比如考虑位移平方项的影响： 2 1 2 d d L L   _ _   图 1-11 非线性问题概念图示 其中， d 为简单拉伸杆右端的位移，图中 P 为节点载荷，由应力与应变之间的线性关系 得到位移与载荷之间的关系： 1 2 EA d d P L L         令 ( ) 1 2 EA d k d L L    _  _  ，于是有： ( ) k d dP 其中 ( )k d 就是刚度系数，并且此刚度系数与位移是相关的。显然这种情况下节点位移与 节点载荷之间是非线性关系，图 1-11 右边的曲线就是这一关系的图示。由于应变与位移的 这种非线性关系，必然导致有限元分析的整体刚度方程为非线性方程，这种问题为几何非线 性问题。 （2）在同一简单直杆拉伸问题中，若应力与应变的关系是非线性的，比如取： (1 0.05 ) E      于是有： EA 1 0.05d d P L L         令 ( )k d EA 1 0.05d L L    _  _  ，则有： ( ) k d dP 由上式看出，应力以及应变之间的非线性关系将直接导致刚度方程的非线性，这类问题 被称为材料非线性问题。 （3）若直杆的右边放置一长为 l 的杆件，且拉伸杆与右边杆件端部之间有微小的间隙  ， 如图 1-12 所示。 变形几何关系为： 当 d≤时： EA d P L  E A L P P d ( ) Pk d d

图 1-12 状态非线性问题简例 当 d   时： 1 1 1 1_{( )} E A EA E A d d d l  L  l   可见，作用于拉伸杆右端的载荷与其右端位移 d 相关，如仍取k EA L  ，且令： 1 1 ( ) ( ) EA d d L P d E A EA d d d L l        _{    }   ≤ 这种情况下，虽然应变与位移之间、应力与应变之间的关系都是线性的，节点力与节点 位移之间的关系在分段加载中也是线性的，但是整个过程仍然表现为节点力与节点位移之间的 非线性，即 ( )k d d P，非线性关系示于图 1-13 中。 图 1-13 状态非线性载荷－位移关系 这种与结构所处的状态相关的非线性问题被称为状态非线性问题，接触问题是最为常见 的一类状态非线性问题。 通过上述简单例子，我们可以总结出各类非线性问题区别于线性问题的基本特点，即刚 度是与位移相关的变化量。 由上面的问题可见，结构非线性问题可以分为以下三大类：  几何非线性：应变与位移关系的非线性。对于几何非线性问题，多指大位移情况下， 尽管位移很大，但结构的应变依然不大，即所谓的大位移小应变问题，只表现为应变 与位移之间的非线性关系，材料的应力－应变关系依然为线性的。也有大位移大应变 的非线性问题，这一类型的问题往往伴随着材料的非线性，如超弹性材料（橡胶）的 大变形问题。 E A L 1 1 E A₁  P d

 材料非线性：材料的本构关系的非线性。对于材料非线性问题，是指材料的本构关系 （即应力－应变关系）为非线性。材料非线性问题中较常见的是结构弹塑性分析。最 常见的材料非线性问题是塑性分析，如果按线性分析得到的应力明显超过屈服强度， 则需要计算塑性响应。  状态非线性：状态变化（包括接触）引起的非线性。许多实际结构常常表现出一种与 状态相关的非线性行为，例如悬索结构中的索可能处于松散的状态，也可能处于绷紧 的状态；轴承套可能是接触的，也可能是不接触的；冻土可能是冻结的，也可能是融 化的。这些系统的刚度由于系统状态的改变在不同的值之间突然变化。状态改变可能 与载荷直接有关（如在悬索问题中），也可能由某种外部原因引起（如引起冻土状态 改变的热力学条件）。此外，结构的施工过程也是一类常见的状态非线性问题，施工 过程各阶段的结构受力情况和结构建成之后是完全不同的。 对结构进行非线性有限元分析，实际上就是求解一组非线性结构方程。非线性方程的求 解方法很多，这里简单介绍 ANSYS 程序提供的有关算法。 （1）Newton-Raphson 方法。 在 ANSYS 中采用增量形式的结构有限元列式，将整个载荷－变形过程划分为一系列增量 段，在每一级增量加载中，为了得到较为精确的载荷－变形过程，对结构的切线刚度矩阵进行 多次修正，并通过修正后的刚度进行迭代，以消除不平衡力，达到误差范围允许的相对平衡状 态后，将这一状态视作平衡状态，继续作用下一级载荷增量。这就是 Newton-Raphson 方法的 基本思路，下面简称为 N-R 方法。图 1-14 是这一方法的示意图，每一次采用当前的切线刚度 得到新的位移增量都相当于一次线性方程组的求解。图 1-14 中载荷是在一步中施加的，总的 载荷 P 也可以分成多个增量步逐级地施加，如图 1-15（a）所示，对每一级都通过多次迭代以 达到收敛容差范围内的平衡。 图 1-14 N-R 方法示意图 （2）修正的 N-R 方法。 N-R 方法可以用于几何非线性程度较高的情况，这一方法所得的计算结果可信度较高，当 然计算的时间也更长。为了避免在同一载荷增量的多次迭代中反复重新形成切线刚度矩阵，提 供了修正的 N-R 方法，即在同一级载荷增量的多次迭代中始终采用本级增量最初的刚度矩阵， 其余做法都和一般的 N-R 法相同。这样计算时间可大为缩短。因为对于非线性分析来说，系 统刚度矩阵的形成耗费的计算时间是比较多的。 （3）弧长方法。 对某些不稳定系统的非线性静态分析（如结构非线性屈曲问题），在使用 N-R 迭代的过程 P P d _d

中，切线刚度矩阵可能变为降秩短阵（奇异阵），导致很严重的收敛问题。对于这种情况，可 以使用另外一种非线性迭代方法，即弧长方法，来达到稳定的收敛。弧长方法的基本原理是： N-R 平衡迭代沿着一段弧收敛，从而即使当切线刚度矩阵的斜率为零或负值时，也可以阻止发 散。这种迭代方法以图形表示在图 1-15（b）中。 ANSYS 程序通过平衡迭代迫使在每一个载荷增量结束时，在某个容差范围内的解答能够 达到相对的平衡。如果不满足收敛准则，重新估算非平衡载荷，持续这种迭代过程直到问题收 敛。ANSYS 程序提供了一系列选项来增强问题的收敛性，如自动时间步、二分法、线性搜索、 预测器、弧长法选项、非线性稳定性技术等，这些控制选项可以在计算过程中根据问题的特点 加以应用，以提高求解的收敛性。 （a）N-R 方法 （b）弧长方法 图 1-15 传统 N-R 方法与弧长方法的比较 F a 3 F a 2 F a 1 F u _u F a 1 F 1 r 2 r 3 r 收敛解 收敛解