书名：有限元基础与COMSOL案例分析
ISBN：978-7-115-62710-0
本书由人民邮电出版社发行数字版。版权所有，侵权必究。
您购买的人民邮电出版社电子书仅供您个人使用，未经授权，不得以任何方式复制和传播本书内容。
我们愿意相信读者具有这样的良知和觉悟，与我们共同保护知识产权。
如果购买者有侵权行为，我们可能对该用户实施包括但不限于关闭该帐号等维权措施，并可能追究法律责任。
著 江 帆 温锦锋 谢智铭 叶宇星
责任编辑 吴晋瑜
人民邮电出版社出版发行 北京市丰台区成寿寺路11号
邮编 100164 电子邮件 315@ptpress.com.cn
网址 http://www.ptpress.com.cn
读者服务热线：(010)81055410
反盗版热线：(010)81055315
本书主要介绍有限元法基础知识及COMSOL在弹性力学、流体力学、电磁学、电化学、多物理场耦合等方面的应用。全书先介绍有限元法的基础知识，然后介绍COMSOL的界面组成与基本操作和网格划分的方法与实例，最后给出了结构力学分析实例、流体力学分析实例、电磁学分析实例、电化学分析实例和多物理场耦合分析实例，即以实例方式介绍COMSOL各方面应用分析的详细操作过程及一些需要注意的问题，多数案例有明确的工程应用背景，部分案例有实验对比结果，具有较强的实用性。
本书可作为机械、材料、水利、土木、暖通、动力、能源、化工、航空、冶金、环境、交通、电力电子、建筑等领域的科研与工程技术人员使用COMSOL软件进行CAE/CFD分析的参考书，也可作为这些专业的本科生和研究生学习有限元及COMSOL软件的教学用书。
有限单元法（简称有限元法，Finite Element Method，FEM）是计算机辅助工程（CAE）的重要手段，它把计算区域离散划分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数（或形函数、插值函数），用单元基函数的线性组合来逼近单元中的真解。若整个计算区域上总体的基函数可以看作是由每个单元基函数组成的，则整个计算区域内的解可以看作是由所有单元上的近似解构成的。有限元法是一种化整为零、集零为整、化未知为已知的方法。有限元分析软件可以帮助企业减少产品设计、优化、控制环节中原型测试的数量和次数。有限元分析软件很多，其中COMSOL Multiphysics具有高效的计算性能和独特的多物理场全耦合分析能力，可以保证数值仿真的高度精确，帮助科研人员与工程师模拟真实世界中的多物理场现象，探索新规律，开发出更好的技术和产品。该软件具有力学、电磁学、流体、传热、化工、MEMS、声学、光学等不同应用领域的专业模块，在机械工程、能源化工、流动与传热、土木工程、航天航空、海洋工程、环境工程、光学、生物医药、微流控、MEMS、电力电子、交通运输、数字孪生等方面都有着广泛的应用。
为了让读者掌握有限元法基础知识与COMSOL各方面的应用操作方法，我们根据COMSOL软件的新版本编写了本书，介绍了弹性力学与有限元法、流体力学、传热、电磁学、化学工程、多孔介质等基础知识及COMSOL在各个工程领域的实例分析。实例涵盖网格划分、结构力学分析、流体力学分析、电磁学分析、电化学分析、耦合分析等方面的工程，尽可能多地涉及了COMSOL的各个应用领域。实例操作步骤详细，并配有演示视频及练习题，便于读者进一步研究。
全书共分8章，第1章介绍有限元法的基础知识，第2章介绍COMSOL的界面组成与基本操作，第3章介绍网格划分的方法与实例，第4章给出了结构力学分析实例，第5章给出了流体力学分析实例，第6章给出了电磁学分析实例，第7章给出了电化学分析实例，第8章给出了多物理场耦合分析实例。本书通过基础知识的介绍，让读者了解软件中一些物理量及其设置的具体意义；利用实例促进软件学习，结合工程实例介绍COMSOL在各个领域应用分析的详细操作过程及一些注意事项，让读者在真实的工程实例分析中明确如何建模，如何进行参数设置，以及如何分析计算结果。
本书由江帆、温锦锋、谢智铭、叶宇星著，江帆负责构思全书框架与实例筛选，并编写第1章，温锦锋编写第5章和第8章，谢智铭编写第2～4章，叶宇星编写第6章和第7章，全书由江帆统稿。许钟松、张冥聪等参与了部分文字编辑工作。
本书部分材料选自百度、COMSOL博客、知乎、哔哩哔哩等网络资源，感谢这些网友所做的工作。本书编写得到广州大学研究生优秀教材资助项目（2022YJS-9）、广州市科技计划项目（202102010386）的支持，得到广州大学机械与电气工程学院老师的支持，得到人民邮电出版社的大力支持，还得到作者的家人所给予的理解与大力支持，在此一并致以深深的谢意！
本书既是工程技术人员利用COMSOL软件进行工程应用分析的指导书，又可作为高等院校相关专业本科和硕士研究生的CAE分析教学、复杂力学问题（特别是多物理场耦合问题）研究的教材或参考书。本书相关实例的模型文件和MPH文件已上传至人民邮电出版社异步社区，扫描书后二维码即可获得相关资料。限于笔者的水平，书中难免有不当之处，还请广大读者给予指正，请致信jiangfan2008@126.com或2908551877@qq.com，不胜感激。
江帆
2023年于广州
有限元法（Finite Element Method，FEM）最初是一种结构分析的方法，现已发展成为一种近似的（除杆件体系结构静力分析外）数值分析方法。其基本思想是将连续的求解区域离散为由有限个单元组成的、按一定方式相互连接在一起的单元集合体，在单元分析中，用近似函数表征待求场函数，建立相关物理量之间的相互联系，然后依据单元之间在结点处的联系，将各单元组装成整体，从而获得整体待求场函数。这种先化整为零，然后集零为整和化未知为已知的研究方法，对求解复杂力学等科学与工程问题具有重要作用。有限元法已经广泛应用于机械工程、热能动力、土木水利、铁道、汽车、船舶工业、航空航天、石油化工、微流控、环境工程、生物医疗、电力等领域。
1.1 弹性力学与有限元在本节中，我们主要介绍弹性力学中的基本概念，包括体力、面力、应力、应变、位移、主应力、相当应力和主应变。
体力是一种外载荷，是随体积分布的力，如重力和惯性力。当材料具有磁性或者分布有非自由电荷，这时磁力和静电力也是体力。体力的单位为N/m3。体力在三个坐标轴上的投影为X、Y、Z，它们的总体可用体力列阵表示：
。
面力也是一种外载荷，是作用在物体表面的力，如接触力和流体压力，其单位为N/m2。面力在三个坐标轴上的投影为
、
、
，它们的总体可用面力列阵表示：
。集中力也是一种面力，它作用在物体表面，忽略其作用面积，认定只作用在一点，单位为N。
应力是单位截面面积的内力（或内力的分布集度），是表示物体内某位置、沿某截面分布内力的大小和方向的物理量。物体由于外力或湿度、温度改变，其内部将产生内力。
在物体内的同一点，不同方向截面上的应力是不同的。过某点，各截面上应力的大小和方向的总和称为该点的应力状态。为了研究某点的应力状态，围绕该点取一个微小单元体，通常用与坐标面平行的平面截出微小的平行六面体，如图1-1所示，单元体三个方向上的尺寸都是非常小的，分别是dx、dy、dz。单元体每个面上的应力分解为一个正应力、两个切应力，分别与坐标轴平行。正应力的作用面用下标表示，如法线平行于x轴的面上正应力记作
；切应力有两个下标，第一个下标为作用面，第二个下标为切应力方向，如
表示其作用面垂直于x轴，方向平行于y轴。正应力与切应力正负号规定如下：正面上与坐标轴正方向一致为正，相反为负；负面上与坐标轴正方向一致为负，相反为正。正应力也简单记作拉为正、压为负。图1-1的单元体上，各面上的应力分量均为正。根据切应力互等定理，6个切应力有3组互等关系，即
，
，
。
图1-1 单元体及其各面上的应力分量
物体中某点的应力状态完全可以由、
、
、、
、
6个应力分量确定，用应力阵列表示为
。
应变反映弹性体在外力作用下，其内部各部分发生变形的程度。变形可以归结为长度的改变与角度的改变，故有线应变与角应变。各线段单位长度的伸缩称为正应变或线应变；各线段之间直角的改变，用弧度表示，称为切应变或剪应变。正应变用
表示，切应变用
表示，用下标表示应变方向，如
表示x方向上的线应变，
表示x、y两方向的线段之间的直角改变量。正应变以伸长为正，缩短为负；切应变以直角变小为正，变大为负。正应变与切应变都是无量纲的量。由切应力互等和胡克定律，得到切应变也是两两互等的，即
，
，
。
物体中某点的变形状态可用
、
、
、
、
、
这6个应变分量确定，用应变列阵表示为
。
位移是指物体受力过程中，物体上各点位置的变化量（位置的移动）。物体内一点（微元体）的位移由刚性位移和弹性位移两部分组成，刚性位移是由其他点的形变引起的位移，弹性位移是本身弹性变形产生的位移，与应变有着确定的几何关系。位移是一个矢量，用
表示，它在空间直角坐标系中，三个坐标方向的位移分量用u、v、w表示。位移分量以沿坐标轴正方向为正，沿坐标轴负方向为负。位移及其分量的单位为m，用一个位移列阵表示为
。
主应力指的是物体内某一点的微面积元上剪应力为零时的法向应力。这时，法向量n=(n1，n2，n3)的方向称为这一点的应力主方向。
弹性体在外力作用下是否会破坏要通过应力来判断，有限元计算的直接应力结果是各点的6个应力分量，其中3个主应力是判断该点材料是否破坏的主要参数。对于不同的失效形式，适用不同的强度理论。根据这些强度理论求得某点的相当应力，用以判断该点强度是否足够。下面介绍几种常用的强度理论及其相当应力。
（1）第一强度理论（最大拉应力理论）。当材料发生断裂且受力弹性体内的某点有拉应力存在，即σ１大于零时，可以按照该点最大拉应力σ１是否小于许用应力来判断该点强度是否足够。第一强度理论的相当应力为
（1-1）
（2）第二强度理论（最大拉应变理论）。当材料发生断裂且受力弹性体内的某点没有拉应力存在，即σ１小于零时，可以按照该点最大拉应变是否小于许用值来判断该点强度是否足够。经过变换得到用主应力表示的第二强度理论的相当应力为
（1-2）
（3）第三强度理论（最大切应力理论）。当材料发生屈服，受力弹性体内的某点最大切应力大于某一定值时，根据第三强度理论，经过变换得到用主应力表示的相当应力为
（1-3）
（4）第四强度理论（最大形状改变比能理论）。当材料发生屈服，受力弹性体内的某点形状改变比能大于某一定值时，根据第四强度理论，经过变换得到用主应力表示的相当应力为
（1-4）
（5）莫尔强度理论。对于一些材料，如铸铁、混凝土等，它们的抗拉能力和抗压能力不同，当它们受到剪切作用、发生剪切破坏时，不仅与切应力大小有关，还与剪切面上的正应力有关，遵从莫尔强度理论，其相当应力为
（1-5）
根据这些强度理论，只要判断某点的相当应力是否小于相应材料的许用应力即可判断该点强度是否足够。
由单元体6个应变分量εx、εy、εz、γxy、γyz、γzx，可以求出过该点任意方向的线应变和任意两线段之间角度的改变，表达式如下。
（1-6）
（1-7）
其中，l、m、n为过物体内某点P的线段PN的方向余弦，l1、m1、n1为过P点与PN成θ角的线段PN1的方向余弦，θ'为物体受力变形后线段PN与PN1的夹角，如图1-2所示。
图1-2 过物体内某点P的线段PN和PN1
进一步分析还可知，物体内任意一点，一定存在3个相互垂直的应变主向，这3个方向的应变称为主应变。3个主应变中最大的一个就是该点的最大线应变，3个主应变中最小的一个就是该点的最小线应变。3个应变主方向与3个应力主方向是重合的，在线弹性范围内，主应力、主应变服从胡克定律（见后续的物理方程）。
弹性力学研究弹性体受外力作用或由于温度变化等原因而发生的应力、应变和位移。弹性体占有三维空间，描述弹性体受力和变形的应力、应变、位移等物理量都是三维坐标的函数。
弹性力学基本方程的导出可从三方面分析：静力学方面，建立应力、体力和面力之间的关系；几何学方面，建立应变、位移和边界位移之间的关系；物理学方面，建立应变与应力之间的关系。通过分析分别得到平衡微分方程、几何方程和物理方程，统称为弹性力学基本方程。
围绕物体内任意一点，取如图1-1所示的一个微小平行六面体，它的3组面平行于3个坐标面，各边长度都是微量dx、dy、dz。外力作用下物体处于静力平衡状态，物体内任意一点也处于静力平衡状态，单元体各面上所受应力及单元体受到的体力满足平衡方程。3个力的平衡方程为
（1-8）
如图1-3所示，可以列出平面问题中的几何方程：
（1-9）
同理可得空间问题的几何方程：
（1-10）
图1-3 平面应变与位移
物理方程体现了应力与应变的关系，也称为本构方程，即
（1-11）
其中，
为体积应变。
上述各方程均可用矩阵方程表示，如式（1-11）可用矩阵方程表示为
（1-12）
简写成σ=Dε。其中，D称为弹性矩阵，它完全由弹性常数E和泊松比μ决定。
弹性力学基本方程共15个，由于平衡方程和几何方程都是微分方程，求解定解还需要边界条件。根据边界条件的不同，弹性力学问题分为位移边界问题、应力边界问题和混合边界问题。
在位移边界问题中，物体在全部边界上的位移是已知的，即
其中，
、
、
在边界上是坐标的已知函数，这就是位移边界条件。
在应力边界问题中，物体在全部边界上的面力分量是已知的。根据面力分量和应力分量之间的关系，可以把面力已知的条件转换成应力方面的已知条件，这就是所谓的应力边界条件，即
（1-13）
其中，面力分量
、
、
在边界上是坐标的已知函数，l、m、n为边界面外法线方向的方向余弦。
在混合边界问题中，物体的一部分边界具有已知位移，即具有位移边界条件，另一部分边界则具有已知面力，即具有应力边界条件。例如，图1-4（a）所示的固定铰支和可动铰支处为位移边界条件，DC边界上分布面力大小为q，其他边界上应力为零，为面力边界条件，整个问题为混合边界问题。图1-4（b）中同一边界存在两种边界条件：x方向位移u|x=a=0；y方向切应力τxy|x=a=0。
图1-4 混合边界问题
任何一个弹性体都是一个空间物体，任何一个实际的弹性力学问题都是空间问题。如果研究的弹性体具有某种特殊形状，并且所受的外力满足一定的条件，就可以把空间问题简化成平面问题。这样处理，分析和计算的工作量将大大减少，而所得的结果仍能满足工程精度要求。
几何特征：一个方向尺寸比另外两个方向尺寸小得多，即t << a，t << b。
受力特征：假设有很薄的等厚平板，在板边上受有平行于板面且不沿厚度变化的面力，同时体力也平行于板面且不沿厚度变化，此类问题就是平面应力问题。
设薄板厚度为t，以薄板的中面（平分板厚的平面）为xy面，z轴垂直于中面，如图1-5所示，因为板面上不受力，所以有
图1-5 平面应力问题
因为板很薄，外力不沿厚度变化，故可以认为整个薄板的所有点都有
由切应力互等关系得
，
。这样只剩下平行于xy面的3个应力分量、、
非零，所以称为平面应力问题。
因为板很薄，3个应力分量、3个应变分量和两个位移分量都可以认为不沿厚度变化，即它们只是x和y的函数，与z无关。
几何特征：一个方向尺寸比另外两个方向尺寸大得多，且沿该方向截面尺寸和形状不变。
受力特征：设有很长的柱形体，在柱侧面上受有平行于横截面且不沿长度变化的面力，同时体力也平行于横截面且不沿长度变化，此类问题就是平面应变问题。
假想该柱体无限长，以任意一横截面为xy面，z轴垂直于xy面，如图1-6所示，则所有应力分量、应变分量和位移分量都不沿z方向变化，它们只是x和y的函数。此外，在这一情况下，由于柱体无限长，任意一横截面都可看作对称面，所有点都只会沿x和y方向移动，而不会有z方向位移，即w=0，εz=γzx=γzy=0，不为零的应变分量只有εx、εy、γxy，所以称为平面应变问题。
图1-6 平面应变问题
平面应力问题和平面应变问题都只有8个独立的未知量σx、σy、τxy、εx、εy、γxy、u、v，它们只是x和y的函数，因此统称平面问题。
平面问题的平衡微分方程为
（1-14）
平面问题中的几何方程为
（1-15）
平面应力问题中的物理方程为
（1-16）
写成矩阵形式为
（1-17）
记作σ=Dε。其中，D为弹性矩阵。
平面应变问题中的物理方程为
（1-18）
写成矩阵形式为
（1-19）
同样记作σ=Dε。其中，D为弹性矩阵。
平面问题的位移边界条件为
（1-20）
平面问题的应力边界条件为
（1-21）
对于弹性力学基本方程，只要给出边界条件，理论上完全可以解出空间问题15个未知量、平面问题8个未知量。这种问题在数学上称为微分方程的边值问题。通常有3种基本解法，即按应力求解、按位移求解和混合求解。按应力求解以应力分量为基本未知函数，先求应力分量，再求其他未知量，是超静定问题，需要补充变形协调条件。由于位移边界条件不能改用应力分量表达，按应力求解时，弹性力学问题只能包含应力边界条件。按位移求解则以位移分量为基本未知函数，此时应通过物理方程和几何方程将平衡微分方程改用位移分量表达。应力边界条件也可以用位移分量表达，按位移求解时，弹性力学问题可以包括位移边界条件和应力边界条件。混合法就是以一部分应力分量为基本未知量，再以一部分位移分量为基本未知量，既建立变形协调方程，又建立内力平衡方程，最后加以求解。不管用哪种方法，工程实际中提出的弹性力学问题，能求得解析解的极其有限，多数还要用数值方法求解。
弹性力学的变分解法属于能量法，是与微分方程边值问题完全等价的方法，将弹性力学问题归结为能量的极值问题。能量表达成位移分量的函数，位移本身又是坐标的函数，因此能量是函数的函数，称为泛函。变分法就是研究泛函的极值问题。
弹性体在外力作用下发生变形，外力对弹性体做功，若不考虑变形中的热量损失、弹性体的动能以及外界阻尼，则外力功将全部转换为储存于弹性体内的位能——应变能。把虚功原理应用于连续弹性体，则可表述为：弹性体在外力作用下处于平衡状态，外力在弹性体所能发生的任何一组虚位移上所做虚功的代数和等于弹性体所储存的虚应变能。
弹性体某位置处在外力作用下实际发生的位移分量u、v、w，既满足位移分量表达的平衡微分方程，又满足边界条件以及用位移分量表达的应力边界条件。假想这些位移分量发生了边界条件所允许的微小改变，即所谓虚位移或位移变分δu、δv、δw成为
，
，
，则外力在虚位移上所做的虚功为
当弹性体发生虚位移后，虚应变能
为应力在虚位移引起的虚应变上所做的虚功。
假定在发生虚位移的过程中，没有其他形式的能量损失，依据能量守恒定理，应变能的增加等于外力在虚位移上所做的功，即虚应变能等于外力虚功。这就是连续弹性体的虚功原理（或称虚位移原理）。
（1-22）
虚功原理（虚功方程）可具体表示为
（1-23）
也可以写成矩阵形式，即
（1-24）
其中，[δ*]为虚位移列阵，[Fb]为体力列阵，[FN]为面力列阵，[ε*]为虚应变列阵，[σ]为应力列阵。
外力从位移状态
退回到无位移的初始状态时所做的功，称为外力势能，记为E。弹性体在这个退回过程中，内部产生变形势能（应变能），记为U。
（1-25）
（1-26）
变形势能和外力势能的总和称为总势能
，
。从前面的虚功原理看到，位移状态d为真实位移状态的充分必要条件是：对应位移d的总势能一阶变分为零，即对应位移d的总势能取驻值。满足位移边界条件的所有位移中，实际发生的位移使弹性体的势能最小，这就是最小势能原理。
（1-27）
对比虚功原理与最小势能原理，可知二者是完全等价的，一个用功的形式表达，另一个以能的形式表达。通过运算，可以由它们导出平衡微分方程和应力边界条件。
弹性力学中的三大类变量为位移、应变和应力，三大类方程是平衡方程、几何方程和物理方程。一般求解弹性力学问题的方法包括以位移为基本未知量的位移法，以应力为基本未知量，通过假设应力函数进行求解的逆解法和半逆解法。一般来说，能够直接进行解析求解的弹性力学问题相当少，而绝大多数弹性力学问题的求解需要借助于数值解法。有限元法为偏微分方程（组）提供了有效的数值近似求解手段，是数值求解弹性力学问题的重要途径。目前，大多数商业有限元软件在对弹性力学问题的分析中采用了按位移求解的方法。
有限元分析的基本步骤主要包括前处理、求解和后处理三个阶段。前处理阶段是将问题的求解域离散成有限个结点和单元，以结点的某些物理量作为基本未知量（在弹性力学问题中，一般是结点的位移），对单元进行分析，构造描述单元物理属性的形函数，描述每个单元的解答并建立起单元刚度矩阵，组装单元，形成总体刚度矩阵，并施加载荷、边界条件和初值条件；求解阶段一般是求解大型稀疏线性方程组，得到各个结点的位移值；后处理阶段是在得到结点的位移值后，进一步计算应力、应变、主应力、相当应力等，例如考虑屈服的强度问题中所需的von Mises应力等。下面通过一个简单的实例向读者逐步说明有限元求解的各个具体步骤。
假设有一个如图1-7所示的变截面直杆，其一端承受10kN集中力，上端的横截面积为100mm2，下端的横截面积为50 mm2，杆的长度为1m，弹性模量为E=200GPa，利用有限元法分析沿杆长度方向上不同点变形的大小及其轴向应力（忽略杆件的自重）。
图1-7 变截面直杆
（1）单元离散。离散化为有限个结点和单元，简单起见，我们将变截面杆离散化为4个单元、5个结点，如图1-8（a）所示（说明：理论上离散的结点和单元数越多，结果可能越精确，但带来的计算量也随之增大）。
（a） （b） （c）
图1-8 将杆离散化成单元
（2）材料力学描述。假定一个近似描述单元特性的解。对于横截面积为A，长度为l，弹性模量为E的等截面直杆，在受到轴向力F的作用下且材料处于线弹性阶段时，由材料力学可知：
（1-28）
（1-29）
其中，σ是杆件横截面的正应力，ε是杆件横截面的正应变，
是杆件的伸长量，FN是横截面上的轴向内力。
类比线性弹簧的控制方程F=kx，将方程式（1-28）和式（1-29）合并为
（1-30）
取等效刚度
。由于直杆在y方向上横截面积是有变化的，作为近似，可以将该杆件看作一系列受到轴向载荷作用的具有不同横截面积的阶梯杆，如图1-8（b）所示，亦可以看作由4根具有不同等效刚度的弹簧串联而成的模型，如图1-8（c）所示。
设第i个结点的位移是u，结点1的约束力为FR，对各个结点进行受力分析，如图1-9所示。建立平衡方程如下：
图1-9 结点受力分析
（1-31）
其中，k1、k2、k3、k4分别表示4个单元的等效刚度。
将式（1-31）重组，并改写成矩阵形式，分离出载荷和约束力，则有
（1-32）
为了不同时考虑未知的约束力和位移，我们将已知结点1的零位移，即用u1=0取代方程式（1-32）的第一行，将式（1-32）改写为
（1-33）
式（1-33）也可表示为
KU=F （1-34）
其中，K是系统的整体刚度矩阵，U是结点位移矩阵，F是载荷矩阵。式（1-33）中的5个方程包含5个未知量，可以求出全部的结点位移，然后可以利用式（1-32）求出约束力。
（3）单元刚度矩阵。前述方法直接求出了系统的整体刚度矩阵，若分析其中的一个单元，则会发现这种单元只有轴向力的作用，对任一结点只有一个轴向位移。考虑弹簧单元内力和位移的关系，根据图1-10所示第i个结点和第i+1个结点的单元内力fi和fi+1，有
（1-35）
图1-10 弹簧单元
列向量
称为单元的结点力向量，矩阵
称为单元刚度矩阵，列向量
称为单元的结点位移列阵。这三者亦有如下关系：
（1-36）
（4）单元组装。将式（1-35）所描述的单元刚度方程应用到所有单元，并将它们组合成整体刚度矩阵K。以单元②为例，该单元连接结点2和结点3，因此
。用k(2G)表示单元②进行扩展后的矩阵，用以进行单元刚度矩阵的组装，也能够清晰地看出单元②在整体刚度矩阵中的位置：
通过在单元刚度矩阵的右边列上结点位移，可以帮助观察该结点对临近单元的影响。同样可以写出单元③扩展后的矩阵：
将各个单元在整体刚度矩阵中的位置进行组合，即对扩展后的矩阵相加，得到最后的整体刚度矩阵K：
（1-37）
其中，n表示单元的总数，这里n=4。
因此，整体刚度矩阵可显式地写为
（1-38）
此式与式（1-32）的刚度矩阵一样。有了整体刚度矩阵以后，我们就可以施加边界条件，并进一步根据式（1-34）进行求解。
假定各单元的长度相等，并以单元的平均横截面积作为单元的计算截面面积，计算等效刚度
。各单元属性见表1-1。
表1-1 单元属性
| 单元号 |
结点号 |
平均横截面积Ai/mm2 |
长度li/mm |
等效刚度keqi/N·mm−1 |
|---|---|---|---|---|
| 1 |
1、2 |
93.75 |
250 |
75 000 |
| 2 |
2、3 |
81.25 |
250 |
65 000 |
| 3 |
3、4 |
68.75 |
250 |
55 000 |
| 4 |
4、5 |
56.25 |
250 |
45 000 |
根据式（1-38）算出整体刚度矩阵，有
引入边界条件u1=0，F=10×103N，根据式（1-33），有
划去第一行和第一列，只需要求解4×4的矩阵方程：
解得：u2=0.133 333 mm，u3=0.287 179 mm，u4=0.468 998 mm，u5=0.691 220 mm。
每个单元的平均应力可按式（1-39）计算：
（1-39）
代入位移计算结果后，求得
无论在何处截断杆件，截面的内力都是FN=10 000 N，因此对于每个单元，可以利用已知材料力学公式
进行验证。结果表明，通过位移信息计算得到的单元应力和采用材料力学公式计算得到的单元应力完全相同，因此就本问题而言，位移计算是正确的。同样也可将结果代入式（1-32），求出约束力FR=10 000 N=10kN，显然该结果与整体平衡条件相符。
说明：这里计算的应力是单元的平均应力，并不是指结点的应力，在有限元法中，经常把与结点相关联的单元的应力进行平均作为结点的应力值。
杆、梁单元一般可以根据其构件之间的连接进行“自然离散化”，而对于连续变形体，需要在对象的几何域上按照分析的需要，采用“人工布置”结点和划分单元的方式进行有限元建模。这种离散方法称为“逼近性离散”，如图1-11所示。
图1-11 逼近性离散
在单元的划分过程中，要求单元之间仅在结点处连接，外载荷将被等效作用到结点上，同时假设单元内部的几何和物理特性是均匀的，这样就把原先连续体的无限自由度问题转变成了近似的有限多个自由度问题。对于弹性力学平面问题而言，用一个全局的位移函数表示整个结构的变形，对于复杂结构几乎是不可能的。但是对于每个单元，采用其结点位移进行插值来构造比较简单的函数，近似表达单元的真实位移，这是简单易行的。当结构被划分成足够多的单元后，把各单元的位移连接起来，就可以近似地表示整个区域的真实位移。
在平面问题的有限元处理中，最常见的单元是3结点三角形平面单元，如图1-12所示。这种单元共有3个结点、6个自由度，对应的有6个结点力。单元的结点位移列阵可以表示为uT=[uix uiy ujx ujy ukx uky]T，单元的结点力列阵可以表示为f T=[fix fiy fjx fjy fkx fky]T。若单元承受分布载荷，可将其等效到结点上。
图1-12 3结点三角形平面单元
从单元的结点位移可以看出，x、y方向上的位移场u（x，y）、v（x，y）可分别由3个结点的x和y方向位移确定，因此分别假设单元中各个方向的位移模式为
（1-40）
其中，a1～a6是待定系数。单元的结点条件为
（1-41）
当已知各结点的位移时，上面的6个方程可以求解6个未知量，计算出a1～a6，并将其回代入式（1-40），经整理，可写成形函数和结点位移向量的乘积形式：
（1-42）
其中，形函数为
（1-43）
系数ai、bi、ci分别按式（1-44）计算：
（1-44）
假设A是单元的面积，有
（1-45）
为了保证面积为正，我们规定结点i、j、k的次序按逆时针排列，如图1-12所示。将式（1-42）写成矩阵形式，有
（1-46）
其中，N是形函数矩阵，I是二阶单位矩阵，ue是单元结点位移向量。对于形函数，同样具有单位分解的特性，即单元中任意点的形函数之和为1；形函数在其对应的结点上值为1，在其他点为0。例如，在结点i上，Ni=1；在结点j、k上，Ni=0。对于Nj和Nk也有同样的性质，这是由插值函数的基本性质所决定的。
在有限元法中，位移模式决定计算误差。载荷的移置以及应力矩阵、刚度矩阵的建立都依赖于位移模式。为了正确反映弹性体中的真实位移形态，位移模式需要满足完备性条件和连续性条件。完备性条件要求位移模式能够反映单元的刚体位移。每个单元的位移不仅包含本单元变形引起的位移，还包括研究对象移动、转动引起的和变形无关的刚体位移，因此位移模式中必须要反映单元的刚体位移。另外，完备性条件还要求位移模式能反映单元的常量应变，每个单元的应变包括与该单元中各点坐标位置无关的应变，即所有点都相同的常量应变。随着单元尺寸的变小，各个点的应变趋于相等，也就意味着单元的变形趋于均匀，常量应变就成为应变的主要部分，这就能保证单元划分逐步增加的情况下，结果趋于真实解。连续性条件则是要求相邻单元在它们的公共结点、公共边上具有相同的位移。
有了单元位移模式，根据几何方程可求出应变：
（1-47）
把形函数表示的位移模式代入式（1-47），则有
（1-48）
其中，矩阵L是算子矩阵，而矩阵B是应变位移矩阵。式（1-43）表明形函数是关于x、y的一次函数，在求一次偏导以后，其结果就为常数，容易算出：
（1-49）
由式（1-44）可知，对于位置确定的单元，其面积A以及各个系数均只和结点的坐标位置有关，所以应变矩阵中的各个元素都是常量。可见，应变ɛ的各分量也是常量。所以，3结点三角形平面单元内部各点的应变均相等，它是常应变单元。
根据物理方程，考虑平面应力问题，有
（1-50）
其中，μ是材料的泊松比。
对于平面应变问题，只需要将弹性模量E替换为
，将泊松比μ替换为
即可。根据式（1-48），应力矩阵可进一步写为
（1-51）
其中，S矩阵即为应力位移矩阵。对于某种弹性模量和泊松比确定的材料，D是常量矩阵，那么S矩阵也是常量矩阵，所以在每一个单元中应力分量也是常量。相邻单元一般不具有相同的应力，在它们的公共边上应力具有突变。但是，随着单元的逐步减小，这种突变将急剧变小，并不妨碍有限元法的解答收敛于正确解答。
通过虚功原理，可以推出3结点三角形平面单元的单元刚度矩阵为
（1-52）
其中，t是单元的厚度，其他符号含义同前。式中的
是2×2的矩阵，可通过式（1-53）计算：
（1-53）
对于平面应变问题，只需要将弹性模量E替换为
，将泊松比μ替换为
即可。
同样也有单元刚度方程
，将其展开，有
（1-54）
单元刚度矩阵具有如下性质。
（1）刚性，单元刚度矩阵描述了单元的刚度特性，即描述了单元受到外力作用时的应变和应力的关系。单元刚度矩阵越大，能够承受越大的力和扭矩。
（2）对称性，即
。
（3）奇异性，即
。单元刚度矩阵的奇异性表明，给定了单元结点载荷列阵，不能得出单元位移列阵，因为即使它满足平衡条件，单元还可以有任意的刚体位移。
（4）主元素恒为正，即单元刚度矩阵的对角线上的元素恒为正。
其他类型的单元刚度矩阵也具有上述性质，对于线性位移模式的三角形单元，可以证明两个相似三角形单元具有相同的单元刚度矩阵；单元水平或竖向移动不会改变刚度矩阵的数值。这种性质使得将求解域划分成相似三角形的单元，只需要计算一次单元刚度矩阵即可。
对于作用在单元内部的体力和边界上的分布面力，需要将它们移置到结点上成为结点载荷。这种移置必须按照静力等效的原则来进行。所谓静力等效，是指原载荷与结点载荷在任何虚位移上的虚功相等。在一定的位移模式下，这种移置的结果是唯一的。对于线性位移模式的3结点三角形平面单元，静力等效意味着原载荷（体力、面力）与结点载荷向任意一点简化得到的主矢和主矩都相等。
设单元在坐标为（x，y）的任意一点处受到集中力P的作用，如图1-13所示，可以将P表示为
，等效结点载荷用
表示，根据静力等效，可导出
（1-55）
图1-13 三角形单元载荷移置
其中，N是形函数矩阵。将其展开，有
（1-56）
若单元上受到分布体力
的作用，则单元结点载荷列阵：
（1-57）
其中，
是单元所覆盖的求解域。
若上述单元的ij边上有分布面力
作用，则
（1-58）
有限元网格中任取一个结点i，如图1-14所示，该结点受到环绕该结点的单元对它的作用力（内力）fi，这些作用力与各单元的结点力等值反向。另外，该结点上还有从环绕该结点的那些单元上移置过来的等效结点荷载（外力）Ri。根据平衡关系，有
图1-14 典型结点
（1-59）
对于所有结点均可以建立这样的方程，若有n个结点，平面问题就有2n个这样的方程。代入结点力公式，可得关于结点位移u的2n个线性代数方程组，即整体的有限元支配方程：
（1-60）
其中，K是整体刚度矩阵，其拼装组合方式与杆、梁单元类似，u为整体结点位移列阵，F是整体载荷列阵。整体刚度矩阵有对称性、奇异性、稀疏性的特点。如果结点、单元合理编号，整体刚度矩阵就具有非零元素带状分布的特点，即非零元素集中在以主对角线为中心的一条带状区域。每行的第一个非零元素到主元素之间的元素个数称为半带宽。半带宽越小，计算机求解的效率越高。
悬臂深梁如图1-15（a）所示，已知右端面作用有均布拉力，其合力为P。若按照图1-15（b）的方式离散化为2个单元、4个结点，设泊松比μ=1/3，厚度为t，求结点位移。
图1-15 3结点三角形平面单元实例
考虑单元①，与结点1、2、3 关联，取i、j、k分别为1、2、3。求出各系数
，
。
单元面积A=1m2。
本实例属于平面应力问题，根据式（1-52）、式（1-53），算出单元刚度矩阵中的各个元素，可得
同理，单元②与结点 1、3、4 关联，取i、j、k分别为1、3、4。算出单元刚度矩阵，即
拼装为整体刚度矩阵，即
则载荷阵列为
其中，F1x、F1y、F4x、F4y实际上是作用在结点1和结点4上的约束力分量。
根据整体有限元支配方程Ku=F，有
考虑边界条件
，划去支配方程的第1、2、7、8 行和第1、2、7、8列，处理后得到
解此矩阵方程，可得
3结点三角形平面单元是有限元法中最早提出的单元，其适应边界能力强，但由于其是常应力单元，因此单元内的应变和应力都是常量，精度相对较低。
图1-16所示的矩形单元也是最基本的单元，它有4个结点、8个自由度。考虑单元的结点位移，x方向的位移场
可以由该方向上的4个结点位移
来确定，而y方向的位移场
可以由该方向上的4个结点位移
来确定。其位移模式为
图1-16 4结点矩形单元
（1-61）
在单元边界上，当x或y是一个常量时，对于单元的每一条边，位移是线性变化的，所以也称为双线性单元。在单元内部，由于存在xy的乘积项，位移是非线性变化的。根据平面问题的几何方程，应变是位移的一阶导数，如
；应力和应变是线性关系，所以4结点矩形单元的应变模式和应力模式是一次线性变化的，其单元内部的应力不再是常量。虽然相邻的矩形单元在公共边界处的应力也有差异，但这种差异较小。在整理应力结果时，采用绕结点平均法，即将环绕某一结点的各单元在该结点处的应力求平均，用来代表该结点处的应力，这种方法的表征性较好。但是，矩形单元也存在比较明显的缺陷，它不能适应曲线边界或斜线边界，也不便在不同部位采用不同大小的单元。为了弥补这种缺陷，我们可以混合使用矩形单元和三角形单元。
在三角形单元的3条边上各增设一个结点，这样每个单元就有6个结点、12个自由度，可以采用二次完全多项式的位移模式。由于应变是位移的一阶导数，而应力与应变是线性关系，因此单元中的应力按照线性变化，能够更好地反映弹性体中应力的变化。在结点数目大致相同的情况下，其精度远高于3结点三角形单元。但是6结点三角形单元对于非均匀性及曲线边界的适应性却不如简单三角形单元，而且由于一个结点的平衡方程涉及较多的结点位移，整体刚度矩阵的带宽较大。
限于篇幅，这里不对4结点矩形单元和6结点三角形单元的单元构造做具体的推导和描述，有兴趣的读者可参阅相关教材和参考书。
在柱坐标下，轴对称问题中的一些非对称力学变量都为0，其三大类力学变量如下所示。
位移：
应变：
应力：
在轴对称问题中，以上这些变量都只和坐标r、z有关，与θ无关。对于轴对称问题的有限元离散，在每一个截面，它的单元形状与一般的平面问题相同，但需要注意的是，轴对称问题的单元都是环形单元，如图1-17所示。
图1-17 轴对称单元
三角形单元能够应用于曲折的几何边界，但其精度较低；矩形单元虽然精度较高，但其适应性较差，不便使用在曲线边界和非正交的直线边界。因此，需要采用一些具有斜边的四边形单元。任意四边形单元可以通过已有的4结点矩形单元进行坐标映射的方式来获得，我们将这种通过坐标映射的方式构造的单元称为参数单元。
在全局坐标系Oxy下的4结点四边形单元的坐标称为物理坐标；其映射母单元是4结点矩形单元，称为基准单元；采用坐标系
，称为基准坐标，如图1-18所示。
图1-18 单元映射
设两个坐标系的坐标映射为
（1-62）
基准坐标中的一点对应于物理坐标中的一个相应点，4个角点的结点映射条件为
（1-63）
每一个方向有4个结点条件，可利用包含4个待定系数的多项式来建立映射关系，即
（1-64）
将式（1-63）代入式（1-64），可以求出待定系数
。将各系数回代入式（1-64），并参照位移模式以形函数和结点位移乘积的模式进行改写，有
（1-65）
其中，
称为几何形状插值函数，具体为
（1-66）
将位移模式表达为基于结点位移的显式表达，即
（1-67）
其中，位移插值函数N(x，y)也可根据式（1-64）改写为
。如果几何形状插值函数与位移插值函数N的阶次相同，则这种单元称为等参元；若的阶次小于N，则称为亚参元。等参元和亚参元能保证单元的收敛。若的阶次大于N，则称为超参元，超参元不能保证单元的收敛。
限于篇幅，这里仅介绍参数单元的基本思路，构造参数单元的单元刚度矩阵、载荷列阵需要用到坐标系之间的偏导数变换，还要利用坐标系之间的雅可比矩阵，同时在计算的过程中需要用到勒让德-高斯求积。具体推导请读者参阅相关教材和参考书。
动力学分析是用来确定惯性和阻尼起作用时结构或构件动力学特性的响应技术，其特性一般包括如下几个方面：①振动特性，包括结构的振动方式与振动频率；②载荷随周期性变化的响应效应，主要是指施加周期性变化的载荷时结构的位移和应力的响应情况；③周期振动或者随机载荷的效应，主要是指结构受周期性载荷或者随机载荷时的变化规律。
结构动力学分析有如下类型。
在设计工程结构或机器部件的振动特性时，设计人员需要进行模态分析，即确定承受动态载荷结构设计中的重要参数（固有频率和振型），同时也可以此作为瞬态动力学分析、谐响应分析等其他动力学分析的基础。
模态分析最终目标是识别系统的模态参数，为结构系统的振动特性分析、振动故障诊断和预报以及结构动力学特性的优化设计提供依据。
模态分析通常求解如下方程：
（1-68）
其中，K为刚度矩阵，M为质量矩阵，u为固有模态位移矢量，
是角频率。
模态问题求解的方法有雅克比法、Givens法、Householder法、对分法、逆迭代法、QR法、子控件法、兰索斯法等。
谐响应分析是用于确定线性结构在承受随时间按正弦（简谐）规律变化的载荷时稳态响应的一种技术，旨在计算结构在几种频率下的响应，并得到一些响应值对频率的曲线。该技术只计算结构的稳态受迫振动，不考虑结构激励开始时的瞬态振动。通过谐响应分析，设计人员能预测结构的持续动力特性，从而验证其设计是否能够克服疲劳、共振及其他受迫振动引起的有害因素。
瞬态动力学分析（也称时间历程分析）是用于确定承受任意的随时间变化载荷的动力学响应的一种求解问题的方法，可用于分析确定结构在稳态载荷、瞬态载荷和简谐载荷的任意组合作用下随时间变化的位移、应变、应力及力。
瞬态动力学分析通常求解如下方程：
（1-69）
其中，K为刚度矩阵，M为质量矩阵，u为位移矢量，C为阻尼系数矩阵（通常表示为比例阻尼，即质量矩阵和刚度矩阵的比例，
）。
瞬态动力学问题的求解方法有中心差分法、纽马克法等。
有限元模型是有限元程序可以处理的对象，是对实际结构的合理模拟。有限元模型一方面要保证力学的完整性，另一方面要保证计算的有效性。也就是说，建立的有限元模型既要承载完整的力学信息，尽可能真实地反映实际情况，又要保证计算机可以快速计算。
在开始建立有限元模型之前，设计人员需要对要解决的问题有透彻的认识，理解问题的力学本质，弄清楚结构几何特征、所受载荷性质、结构材料特性，初步估计响应情况。此外，设计人员需要根据力学概念，分析判断研究对象属于哪一类性质的问题，是线性问题还是非线性问题，是静力问题还是动力问题。
线性问题是指受力与变形成正比关系，例如，做出连续性、均匀性、线弹性、各向同性、小变形等基本假设，最后得出的结构刚度为常量，这就属于线性问题。非线性问题是指受力与变形不成正比关系，引起非线性的因素主要有3种。一是材料本身的刚度随变形而变化，如一般低碳钢受力比较小时，受力与变形呈线弹性关系，受力较大发生屈服时刚度很小，进入强化阶段具有一定刚度但小于线弹性阶段，这类问题称为材料非线性问题。二是结构在载荷作用下发生过大变形，结构形状影响刚度，这类问题称为几何非线性问题。当物体变形的大小与物体某个几何尺寸可以相比拟时，应按大变形来处理；当应变量大于0.3时，应按大应变问题处理。大变形、大应变问题都属于几何非线性问题。三是状态非线性，如结构中两零件为接触关系，随着受力变形，接触位置、接触面积都可能发生变化，接触刚度当然也会发生变化，即状态改变影响刚度，这类问题称为状态非线性问题。
静力问题是指所有变量和关系式都与时间无关，只要任一变量或关系式与时间有关，即属于动力问题。作为动力问题，若要计算结构的固有特性，则属于模态分析；若要计算在随时间变化的载荷的作用下结构各结点随时间变化的位移、速度、加速度及应力等，则属于瞬态响应分析；若要计算在随时间按照正弦或余弦规律变化的载荷的作用下结构的稳态响应，则属于谐响应分析，分析目的在于得到结构在不同频率简谐载荷的作用下响应与频率的关系；若要计算结构在某载荷谱作用下的位移、应力等，则属于谱响应分析，载荷谱可以是位移、速度、加速度或力随频率变化的关系图，谱分析主要用于确定结构对随机载荷或随时间变化载荷（如地震、风载荷等）的动力响应情况，可代替瞬态响应分析。
建模前，设计人员还要根据物体的形状和受力情况，判断是一维问题、二维问题还是三维问题，有无对称性、反对称性或周期对称等可利用。这是一个把工程问题进行力学建模的过程。有了准确的力学模型，有限元建模就有了基础。有限元建模过程包括选择单元类型、确定单元的尺寸大小、保证网格划分质量、定义材料和单元特性、处理载荷和边界条件、确定计算方法和控制参数、求解、输出结果等。如果计算结果揭示了事物的内在规律，说明有限元模型和实际物理模型的力学性能是一致的，那么所建的模型就是个好模型。如果计算结果不符合实际情况或者计算进行不下去，那么可能出现单元类型不对、网格数量太少、材料模型错误、约束和载荷的施加方式不对、接触定义有问题、网格质量差、计算方法不对等问题，建模过程中的每个因素都可能造成计算结果错误或计算困难。针对工程问题，利用有限元分析方法，要想得到一个实用的计算结果，准确建模是关键。
有限元建模的准则是根据工程分析的精度要求，建立合适的、能模拟实际结构的有限元模型。要使分析结果有足够的精度，所建立的有限元模型必须在能量上与原连续系统等价。有限元模型具体应满足下述准则。
（1）满足平衡条件。结构的整体和任意一单元在结点上都必须保持静力平衡。
（2）满足变形协调条件。交汇于一个结点上的各单元在受力变形后也必须保持交汇于同一结点。
（3）满足边界条件和材料的本构关系。边界条件包括整个结构的边界条件和单元间的边界条件。
（4）符合刚度等价原则。有限元模型的抗拉压、抗弯曲、抗扭转、抗剪切刚度应尽可能与原来结构等价。
（5）认真选取单元，包括单元类型、形状、阶次，使其能够很好地模拟几何形状、反映受力和变形情况。单元类型有杆单元、梁单元、平面单元、板单元、空间单元等，空间块体又分四面体块单元或六面体块单元，六面体块单元又分八结点六面体或二十结点六面体等。选取单元时应综合考虑结构的类型、形状特征、应力和变形特点、精度要求和硬件条件等因素。
（6）应根据结构特点、应力分布情况、单元的性质、精度要求及其计算量的大小等仔细划分计算网格。
（7）在几何上要尽可能地逼近真实的结构体，其中特别要注意曲线与曲面的逼近问题。
（8）仔细处理载荷模型，正确生成结点力，同时载荷的简化不应该跨越主要的受力构件。
（9）质量的堆积应该满足质心及惯性矩等效要求。
（10）超单元的划分尽可能单级化并使剩余结构最小。
对于基于位移模式的有限元法，在结构的边界上必须严格满足已知的位移约束条件。例如，弹性体某位置处有固定支撑，这些边界上的位移、转角等于零，如图1-19（a）所示，
，图1-19（b）中，
；或者弹性体某位置处位移或转角有已知值，如图1-19（c）中，
，计算模型必须让它能实现这一点。
图1-19 各种边界条件示例
有时边界支撑不是沿坐标方向，称为斜支撑，如图1-19（d）所示，可以设定与整体坐标不一致的结点坐标方向来实现。还有的约束是单向的，例如，绳索只能承受拉力，光滑支撑面只能提供压力，这就需要按非线性对待。
当边界与另一弹性体相连，构成弹性边界时，我们可分两种情况处理。如果弹性体对边界点的支撑刚度已知，那么可将它的作用简化成弹簧，在此结点上加一弹簧单元，如图1-19（e）所示。如果弹性体对边界点的支撑刚度不清楚，那么可将此弹性体的一部分划出来和结构连在一起进行分析，所划区域的大小视其有影响的区域大小而定，如图1-19（f）所示。
如果整个结构存在刚体位移，就无法进行静力分析、动力分析。为此，我们必须根据实际结构的边界位移约束情况，对模型的某些结点施加约束，消除结构的刚体位移影响。对于平面问题，应消去结构的两个平移自由度、一个转动自由度；对于三维问题，须消去3个平移自由度、3个转动自由度。此外，要保证这些消除模型刚体位移的约束施加得当，如果不恰当，就会产生不真实的支反力，改变了原结构的受力状态和边界条件，从而得到错误的结果。例如，在图1-19（g）中根据对称性，C点两方向位移均为零，因此对C点施加约束是适当的。若把点A、B、D的两方向位移指定为零，则与实际情况不符。
复杂结构通常由杆、梁、板、壳及二维体、三维体等多种形式的构件组成。由于构件中各组件之间的自由度个数不匹配，因此在梁和二维体、板、壳和三维体的连接处必须妥善加以处理，否则模型会失真，得不到正确的计算结果。
在复杂结构中，我们还能遇到各种各样其他的连接关系，只要将这些连接关系彻底弄清楚，就能写出相应的位移约束关系式，这些关系式为构件间复杂的连接条件，需要在计算中使程序严格满足这些条件。
应当指出，在不少的实用结构有限元分析程序中，已为用户提供输入连接条件的接口，用户只需严格遵守用户使用规定，程序将自动处理自由度之间的用户所规定的位移约束条件。
有限元的计算时间和结点数的多少有很大关系。在保证计算精度的条件下，用户应采用各种手段减少结点数，以节约计算时间和成本。
如果计算对象的结构具有对称性，就可利用这个特点减少参加计算的结点数。所谓结构的对称性，是指结构的几何形状和支撑条件对某轴（面）对称，同时截面和材料性质也对称于此轴（面）。也就是说，结构绕对称轴对折后，左、右两部分完全重合。
如果对称结构上有对称载荷作用，则变形和应力也是对称的。只需取一半的结构建模即可，对称轴上的结点给出对称边界条件，算完后还可以根据对称性扩展出另一半结果。这样解题规模可减小一半。
如果作用在对称结构上的载荷是反对称的，即将结构绕对称轴对折后，两载荷的作用点重合，载荷大小相同，但载荷方向相反。根据结构力学可知，在反对称载荷作用下，结构的位移及应力都将反对称于对称轴。
在常用的商业有限元软件中，只要用户给出对称或反对称条件，程序会自动加上相应的位移约束。
有些结构可以划分为若干形状完全相同的子结构，当任意一子结构绕对称轴旋转一定角度时，该子结构的形状将与其他子结构完全重合，这种结构称为循环对称结构或周期对称结构。工程中常见的风扇叶片、花键、螺旋桨、齿轮、法兰等都是周期对称结构。如果结构所受载荷和位移的约束也是周期对称的，且各子结构材料和物理特性也完全相同，则应力和变形关于同一轴周期对称。若所受载荷不是周期对称的，如齿轮、法兰，则不属于周期对称问题。对于周期对称问题，计算时可以只取一个子结构进行分析。注意：在取一个子结构时，应使应力集中区域在子结构内部而不在边界。
另外，还有一种周期对称结构可以看作由一个子结构沿某一方向多次重复得到，称为重复对称结构。如果结构所受载荷和约束同样满足重复对称条件，与循环对称类似，那么只需要模拟和分析一个子结构即可。
对于一个复杂的工程构件，设计人员可以根据其在几何学、力学或传热学上的特点，进行降维处理。对于一个三维物体，如果可以忽略某些几何上的细节或次要因素，就能按照二维问题来处理。例如螺纹连接结构中，由于螺纹升角很小，可认为螺纹牙的受力在周向是相同的，从而近似看作轴对称结构。一个二维问题，若能近似地看作一维问题，就尽量当作一维问题计算。维数降低，计算量将大大降低。例如，齿轮、连杆、径向轴承等许多零件的结构计算都可以近似作为平面问题。在复杂的结构计算中，应尽量减少按三维问题处理的部分。
此外，某些零件上会有许多小圆孔、小圆角、小凸台、浅沟槽等几何细节，细节的存在将影响网格的大小、数量及分布。因为在自动划分网格时，一段直线或曲线至少划分一个单元边，一个平面或曲面至少划分一个单元面，一个圆最少也要3个单元边来离散，所以细节将限制网格的大小。另外，单元由密到疏应该平缓过渡，这也会影响整个模型的网格数量和分布。但细节的取舍要遵循两条原则：一是细节处应力的大小，只要这些不是位于应力峰值区域中分析的要害部位，根据圣维南原理就可以将其忽略；二是与分析的内容也有关系，一般情况下，由于细节会影响应力的大小及分布，静应力、动应力计算中要注意细节的影响，而结构的固有频率和模态振型主要取决于质量分布和刚度，因此计算固有特性时就可以少考虑细节。
现代机械设计中进行力学计算的目的，往往在于求出结构最大承载能力或结构最薄弱区域，这些处理方法虽然会带来一定的误差，但一般都能满足工程上的设计要求，而计算成本却能大大降低。如果对个别部位分析后不能满意，则可将这块单独取出再作细致分析。
当计算的结构比较复杂时，整体刚度矩阵的阶数往往会很大，从而超出计算机容量，这时可以考虑一小块一小块地来计算，最后再将各子块边界结点归结在一起，这就是子结构分析法。
子结构方法还可以用在需要局部精确分析的场合，如应力集中位置、局部发生塑性变形需要进行非线性分析的地方、设计可能改变的局部等，设计人员可以只重复计算部分结构，节约计算时间和计算成本。
子结构分析法的基本思路：①几何分割；②子结构离散；③定义边界自由度；④凝聚内部自由度；⑤子结构集成；⑥求解整体模型；⑦回代。现有大型有限元程序一般包括子结构法的内容，用户根据需要调用即可。
在工程上，我们常常将一些呈微弱非线性的问题当作线性问题来处理，所得结果既能满足工程要求，又可降低成本。例如，许多混凝土结构（水坝、高层建筑、冷却塔、桥梁、大型机电设备地基等）实际上都是非线性结构，其非线性现象较弱，初步分析时，常看作线性结构。只有当分析其破坏形态时，才按非线性考虑。
有时我们要对一个结构进行多种载荷工况的分析，如果每一种工况都作为一个新的问题重新分析一次，则计算量会很大，也没有必要。对于上述情况，用户可以将每一种载荷矢量
合并成载荷矩阵R，一起进行求解。方程系数只需进行一次三角分解，计算量就会大大降低。对于线性问题，用户还可以先解出某些标准载荷模式
、
、
下的解
、
、
，若其他载荷模式可以写成这些载荷的线性组合
，则它对应的解也是这些解的线性组合，即
，其中a、b、c为线性组合系数。
有限元计算需要输入和存储大量信息，计算量也非常大，这就要求编程人员探索减少存储量、减少运算次数的方法。
有限元算法计算量大致与总体刚度带宽的平方成正比。为了减少存储量，我们可以根据总体刚度矩阵具有对称性、稀疏性、带状分布等特点，采用一种名为“半带存储”的技术，即只存储总体刚度矩阵中沿主对角线非零元素的一半的数据。对同一个模型，如果按不同的次序对各结点进行编号，那么得到的总体刚度矩阵形式是不同的，半带宽也不一样，相应存储量也就不同。
1.2 流体力学基础本节主要介绍流体力学中的一些基本概念。
（1）流体质点。几何尺寸同流动空间相比是极小量，又含有大量分子的微元体。
（2）连续介质。质点连续地充满所占空间的流体或固体。
（3）连续介质模型。连续介质的所有物理量都是空间坐标和时间的连续函数的一种假设模型，可表示为u =u(t,x,y,z)。
（1）惯性。流体不受外力作用时，保持其原有运动状态的属性。惯性与质量有关，质量越大，惯性就越大。单位体积流体的质量称为密度，以
表示，单位是kg/m3。对于均质流体，设其体积为V，质量为m，则密度为
。对于非均质流体，密度随点而异。
（2）压缩性。作用在流体上的压力变化可引起流体的体积变化或密度变化，这一现象称为流体的可压缩性。压缩性可用体积压缩率k来量度，即
（1-70）
其中，p为外部压强。
在研究流体流动过程中，若考虑流体的压缩性，则称之为可压缩流动，该流体称为可压缩流体，如高速流动的气体。若不考虑流体的压缩性，则称之为不可压缩流动，该流体称为不可压缩流体，如水、油等。
（3）黏性。在运动的状态下，流体所产生的抵抗剪切变形的性质。黏性大小用黏度来量度。流体的黏度是由流动流体的内聚力和分子的动量交换所引起的。黏度有动力黏度μ和运动黏度
之分。动力黏度由牛顿内摩擦定律导出，即
（1-71）
其中，τ表示切应力，单位为Pa；μ表示动力黏度，单位为Pa·s；du/dy表示流体的速度梯度。
运动黏度与动力黏度的关系为
（1-72）
其中，为运动黏度，单位为m2/s。
在研究流体流动过程中，若考虑流体的黏性，则称之为黏性流动，相应的流体称为黏性流体；若不考虑流体的黏性，则称之为理想流体的流动，相应的流体称为理想流体。
根据是否满足牛顿内摩擦定律，流体可以分为牛顿流体和非牛顿流体。牛顿流体严格满足牛顿内摩擦定律且μ保持为常数。非牛顿流体的切应力与速度梯度不成正比，一般又分为塑性流体、假塑性流体和胀塑性流体。
（1）质量力。质量力是与流体微团质量大小有关并且集中在微团质量中心的力。在重力场中有重力mg；直线运动时，有惯性力ma。质量力是一个矢量，一般用单位质量所具有的质量力来表示，其形式如下：
f=fxi + fy j + fzk （1-73）
其中，i、j、k分别为x、y、z轴方向的单位矢量。
（2）表面力。表面力是大小与表面面积有关而且分布作用在流体表面上的力。表面力按其作用方向可以分为两种：一是沿表面内法线方向的压力，称为正压力；另一种是沿表面切向的摩擦力，称为切向力。
对于理想流体的流动，流体质点只受到正压力，没有切向力。对于黏性流体流动，流体质点所受到的作用力既有正压力，也有切向力。
作用在静止流体上的表面力只有沿表面内法线方向的正压力。单位面积上所受到的表面力称为这一点处的静压强。静压强具有两个特征：①静压强的方向垂直指向作用面；②流场内一点处静压强的大小与方向无关。
（3）液体的表面张力。在液体表面，界面上液体间的相互作用力为张力，液体表面有自动收缩的趋势，收缩的液面存在与该处液面相切的拉力。正是这种力的存在，使弯曲液面内外出现压强差以及常见的毛细现象等。
实验表明，表面张力大小T与液面的截线长度L成正比，即
（1-74）
其中，
称为表面张力系数，它表示液面上单位长度截线上的表面张力，其大小由液体性质与接触相温度、压力等决定，其单位为N/m。
（4）相对压强、绝对压强及真空度。标准大气压的压强是101 325Pa（760mmHg），它是压强的一个单位，记作atm。若压强大于大气压，则以此压强为计算基准得到的压强称为相对压强，也称表压强，通常用pr表示。若压强小于大气压，则压强低于大气压的值就称为真空度，通常用
表示。如以压强0Pa为计算的基准，则这个压强就称为绝对压强，通常用
表示。这三者的关系如下：
（1-75）
（1-76）
在流体力学中，压强都用符号p表示，但一般有一个约定：对于液体，压强用相对压强；对于气体，特别是马赫数大于0.1的流动，应视为可压缩流，压强用绝对压强。
压强的单位较多，一般用Pa，也可用单位bar，还可以用毫米汞柱、毫米水柱，这些单位换算如下：1Pa＝1N/m2；1bar＝105Pa；1atm＝760mmHg＝10.33mH2O＝101 325Pa。
（5）静压强、动压强和总压强。对于静止状态下的流体，只有静压强。对于流动状态的流体，有静压强、动压强、测压管压强和总压强之分，我们可以从伯努利方程中分析它们的意义。
伯努利方程的物理意义是一条流线上流体质点的机械能守恒。对于理想流体的不可压缩流动，其表达式为
（1-77）
其中，
称为压强水头，也是压能项，为静压强；
称为速度水头，也是动能项；
称为位置水头，也是重力势能项：这三项之和就是流体质点的总机械能。H称为总的水头高。
将式（1-77）两边同时乘以
，则有
（1-78）
其中，
称为静压强，简称静压；
称为动压强，简称动压；
称为总压强，简称总压。对于不考虑重力的流动，总压就是静压和动压之和。
根据流体流动过程以及流动过程中流体的物理参数是否与时间相关，流动可以分为定常流动与非定常流动。
（1）定常流动。流体流动过程中各物理量均与时间无关。
（2）非定常流动。流体流动过程中某个或某些物理量与时间有关。
常用迹线和流线来描述流体的流动。
（1）迹线。随着时间的变化，空间某一点处的流体质点在流动过程中留下的痕迹称为迹线。在t＝0时刻，位于空间坐标(a,b,c)处的流体质点，其迹线方程为
（1-79）
其中，u、v、w分别为流体质点速度的3个分量；x、y、z为在t时刻此流体质点的空间位置。
（2）流线。在同一个时刻，由不同的无数多个流体质点组成的一条曲线，曲线上每一点处的切线与该质点处流体质点的运动方向平行。流场在某一时刻t的流线方程为
（1-80）
对于定常流动，流线的形状不随时间变化，而且流体质点的迹线与流线重合。在实际流场中，除驻点或奇点外，流线不能相交，不能突然转折。
（1）流量。单位时间内流过某一控制面的流体体积称为该控制面的流量Q，其单位为m3/s。若单位时间内流过的流体是以质量计算，则称为质量流量
。若不加说明，则“流量”一词泛指体积流量。在曲面控制面上有
（1-81）
（2）净通量。在流场中取整个封闭曲面作为控制面A，封闭曲面内的空间称为控制体。流体经一部分控制面流入控制体，同时也有流体经另一部分控制面从控制体中流出。此时流出的流体减去流入的流体，所得流体称为流过全部封闭控制面A的净通量（或净流量），通过式（1-82）计算：
（1-82）
对于不可压缩流体来说，流过任意封闭控制面的净通量等于0。
由速度分解定理，流体质点的运动可以分解为随同其他质点的平动、自身的旋转运动和自身的变形运动（拉伸变形和剪切变形）。
在流动过程中，若流体质点自身做无旋运动，则称流动是无旋的，也就是有势的，否则就称流动是有旋流动。流体质点的旋度是一个矢量，通常用ω表示，其大小为
（1-83）
若ω＝0，则称流动为无旋流动，即有势流动，否则就是有旋流动。
ω与流体的流线或迹线形状无关。黏性流动一般为有旋流动。对于无旋流动，伯努利方程适用于流场中任意两点之间。对于无旋流动（也称为有势流动），存在一个势函数
，满足：
（1-84）
即 
（1-85）
流体的流动分为层流流动和湍流流动。层流流动中流体层与层之间相互没有任何干扰，层与层之间既没有质量的传递，也没有动量的传递；而湍流流动中层与层之间相互有干扰，而且干扰的力度还会随着流动而加大，层与层之间既有质量的传递，又有动量的传递。
判断流动是层流还是湍流，需要看其雷诺数是否超过临界雷诺数。雷诺数的定义为
（1-86）
其中，V为截面的平均速度；L为特征长度；
为流体的运动黏度。
对于圆形管内流动，特征长度L取圆管的直径d，即
（1-87）
一般认为临界雷诺数为2320。当Re < 2320时，管中是层流；当Re > 2320时，管中是湍流。
对于异型管道内的流动，特征长度取水力直径dH，则雷诺数的表达式为
（1-88）
异型管道水力直径的定义为
（1-89）
其中，A为过流断面的面积；S为过流断面上流体与固体接触的周长。
1.3 CFD基本模型流体流动所遵循的物理定律是建立流体运动基本方程组的依据。这些定律主要包括质量守恒、动量守恒、动量矩守恒、能量守恒、热力学第二定律，加上状态方程、本构方程。在实际计算时，还要考虑不同的流态，如层流与湍流。
在流体力学中，系统是指某一确定流体质点集合的总体。系统以外的环境称为外界。分割系统与外界的界面称为系统的边界。系统通常是研究的对象，外界则用来区别于系统。系统将随系统内质点一起运动，系统内的质点始终包含在系统内，系统边界的形状和所围空间的大小，则可随运动而变化。系统与外界无质量交换，但可以有力的相互作用及能量（热和功）交换。
控制体是指在流体所在的空间中，以假想或真实流体边界包围，固定不动、形状任意的空间体积。包围这个空间体积的边界面称为控制面。控制体的形状与大小不变，并相对于某坐标系固定不动。控制体内的流体质点组成并非不变的。控制体既可通过控制面与外界有质量和能量交换，也可与控制体外的环境有力的相互作用。
本书将用到如下一些数学运算符。
梯度 
（1-90）
散度 
（1-91）
旋度 
（1-92）
其中，
，称为那勃勒算子。
（1-93）
（1-94）
（1-95）
其中，
称为拉普拉斯算子。
冒号运算符 
。 （1-96）
在流场中，流体通过控制面
流入控制体，同时也会通过另一部分控制面A2流出控制体，在这期间控制体内部的流体质量会发生变化。按照质量守恒定律，流入的质量与流出的质量之差，应该等于控制体内部流体质量的增量，由此可导出流体流动连续性方程。
（1-97）
其中，
表示密度，u表示速度矢量。
动量守恒是流体运动时应遵循的另一个普遍定律，描述为：在一给定的流体系统，其动量的时间变化率等于作用于其上的外力总和，其数学表达式即为动量守恒方程，也称为运动方程，或N-S方程，其表达式为
（1-98）
其中，是压力，I是单位矩阵，K为黏性应力张量，F是体积力矢量。
动量守恒方程在实际应用中有许多表达形式，需要根据实际计算情况来选择使用。
将热力学第一定律应用于流体运动，把流体相对运动方程中的各项用有关的流体物理量表示出来，即得到能量守恒方程。
（1-99）
其中，Cp是恒压比热容；T是绝对温度；q是热通量矢量；Q为热源；S为应变率张量，
；K为黏性应力张量，
。同样，该式在实际应用中有许多表达形式，需要根据实际计算情况来选择使用。
湍流流动是自然界广泛存在的现象，其核心特征是其在物理上近乎于无穷多的尺度和数学上强烈的非线性，这使得人们无论是通过理论分析、实验研究还是计算机模拟来彻底认识湍流都非常困难，因此研究湍流机理，建立相应的模式，并进行适当的模拟仍是解决湍流问题的重要途径。COMSOL提供的湍流模型包括Spalart-Allmaras模型、L-VEL模型、代数yPlus模型、标准
模型、可实现的（Realizable）模型、低雷诺数模型、
模型、SST模型、v2-f模型等。
选取湍流模型时，需要考虑的因素包括流体是否可压、针对特定问题的习惯解法、精度的要求、计算机的计算能力和时间的限制。COMSOL还有壁面函数、自动壁面处理、湍流模型间自动切换等方式和方法帮助用户解决湍流求解问题。
Spalart-Allmaras 模型增加了一个额外的无衰减运动学涡流黏度变量。它是一个低雷诺数模型，可求解实体壁之内的整个流场。这个模型最初针对空气动力学应用而开发，优势在于相对稳健，分辨率要求不高，内存需求小，具有良好的收敛性，不使用壁面函数使可精确计算力（升力与曳力）、流量（传热与传质）。该模型不能精确计算包含剪切流、分离流或衰减湍流的流场。
L-VEL和代数yPlus湍流模型仅基于局部流速和与最近壁面的距离来计算湍流黏度，它们不求解附加变量。这两种模型的鲁棒性好，且计算强度低。虽然它们是精度较低的模型，但对内部流动却有很好的近似，尤其是在电子冷却应用中。
标准k-ε模型求解了两个变量：湍流动能k和湍流动能耗散率ε。本模型使用了壁面函数，但壁附近的解不够精确。标准k-ε模型稳定，具有很好的收敛速率和相对较低的内存要求，在工业领域应用广泛。标准k-ε模型可以在壁附近使用较粗的网格，对于复杂几何形状外部流动问题的求解效果很好，如标准k-ε模型可用于求解钝体周围的气流。但它不能精确地计算流动或射流中的逆压梯度和强曲率的流场。标准k-ε模型的湍流动能k和耗散率ε的方程为
（1-100）
（1-101）
其中，
为湍流黏度，
；常数
，
，
，
，
；
为产生项，表达式如下：
（1-102）
可实现的k-ε模型与标准k-ε模型相比，有两个主要不同点：①可实现的k-ε模型为湍流黏性增加了一个公式。②可实现的k-ε模型为耗散率增加了新的传输方程。除强旋流过程无法精确预测外，其他流动都可以使用此模型来模拟，包括有旋均匀剪切流、自由流（射流和混合层）、腔道流动和边界层流动。
低雷诺数k-ε模型类似于标准k-ε模型，但没有使用壁面函数。它求解了每个位置的流动，是对标准k-ε模型的合理补充，拥有和后者一样的优势，但通常要求网格更加密集；它的低雷诺数属性不仅表现在壁面上，而是在各处都发挥作用，使湍流衰减。该模型有两种常用的方法：一种方法是首先使用标准k-ε模型计算出一个良好的初始条件，然后用它求解低雷诺数k-ε模型；另一种方法是使用自动壁面处理功能，先利用粗化的边界层网格来获取壁面函数，然后对所需壁面处的边界层进行细化，进而获得低雷诺数k-ε模型。
低雷诺数k-ε模型可以计算升力和曳力，而且热通量的建模精度远远大于标准k-ε模型。在许多情况下，它表现出了卓越的预测分离和黏附的能力。
k-ω模型通过两个输运方程求解
与
。对于有界壁面和低雷诺数的可压缩性和剪切流动，该模型能取得较好的模拟效果，尤其适合处理圆柱绕流、放射状喷射、混合流动等问题，它包含转捩、自由剪切和压缩性选项。
SST模型结合了自由流中的k-ε模型和近壁的k-ω模型。它是一个低雷诺数模型，在工业应用中是一个“万能”模型。在对分辨率的要求方面，该模型与k-ω模型和低雷诺数k-ε模型相似，但消除了k-ω模型和低雷诺数k-ε模型表现出的一些弱点。
在接近壁面边界的地方，平行方向上的速度脉动通常会远远大于垂直于壁面的方向，速度脉动被认为是各向异性的。在远离壁面的地方，所有方向的脉动大小均相同，速度脉动变为各向同性。
除了使用两个分别描述湍流动能k和耗散率ε的方程，v2-f湍流模型使用了两个新方程来描述湍流边界层中湍流强度的各向异性：第一个方程描述了垂直于流线的湍流速度脉动的传递；第二个方程解释了非局部效应，例如由壁面引起的、垂直和平行方向之间的湍流动能的再分配阻尼。
大涡模拟（Large Eddy Simulation，LES）用于解析较大的三维非定常湍流涡，而小涡流的影响则通过近似方法表示。这项技术与边界层网格划分一起使用时，可以精确描述瞬态流场以及边界上的精确通量和力。COMSOL中提供的 LES 模型包括“基于残差的变分多尺度”（RBVM）、“基于残差的黏性变分多尺度”（RBVMWV）和 Smagorinsky 模块。
在流体动力学计算中，初始条件和边界条件的正确设置是关键的一步。COMSOL软件提供了流动的初始条件和边界条件。
初始条件是计算初始给定的参数，即
时给出各未知量的函数分布。初始条件需要根据实际情况来设置。当流体运动定常时，无初始条件问题。
边界条件是流体力学方程组在求解域的边界上流体物理量应满足的条件。例如，流体被固壁所限，流体就不应有穿过固壁的速度分量；在水面边界上，大气压强认为是常数（一般在距离不大的范围内可如此）；在流体与外界无热传导的边界上，流体与边界之间无温差等。虽然各种具体问题不同，但边界条件一般要保持恰当：①保持在物理上是正确的；②要在数学上不多不少，刚好能用来确定微分方程中的积分常数，而不是矛盾的或有随意性。
COMSOL软件的常用流动分析的初始条件与边界条件设置详见5.1节。
1.4 多相流多相流通常包括气-液、液-液、液-固、气-固、气-液-液、气-液-固、气-液-液-固混合物的流动，在机械、能源、动力、核能、石油、化工、冶金、制冷、运输、环境保护及航天技术等许多领域都有其踪迹。多相流的建模计算通常比较难，尤其是追踪流体与流体之间的交界面更难，这里简单介绍COMSOL中多相流的相关内容。
在COMSOL中，多相流模型分为两大类：分离多相流模型与分散多相流模型。
在较小尺度上，我们可以对相边界的形状进行详细建模，并把这种模型称为分离多相流模型。分离多相流模型有清晰的相界面，对此我们一般用相场
描述相界面，
是一相，而
是另一相。分离多相流模型主要用于气泡、液滴和颗粒流的模拟分析，这些相界的尺度量级与流场尺度相当，而且数量较少；也可用于微流体中的多相流、宏观流场中的单相流的自由液面。通常，我们使用表面追踪法来描述此类模型。
在较大尺度上，如果仍详细描述相边界，则模型方程无法求解。这时我们可以使用体积分数场描述不同的相，并把这种模型称为分散多相流模型。在分散多相流模型方程中，相间效应（例如表面张力、浮力和跨越相边界的传递）被视为源和汇。分散多相流模型的相界面不太清晰，为此我们用体积分数场描述相间的关系，其中
。分散多相流模型主要用于气泡数量多而气泡体积小的气泡流，也可用于乳液和气溶胶的模拟，还可用于流场中有大量固体颗粒的情况以及宏观的多相流分析。
如图1-20所示，分离多相流模型详细描述了相边界，分散多相流模型则只考虑分散在连续相中的一个相的体积分数。在分离多相流模型中，不同相之间相互排斥，并存在一个清晰的相边界，在此边界上相场函数发生突变。除了追踪相边界的位置，相场函数没有任何物理意义。在分散多相流模型中，函数描述了气相（分散相）和液相（连续相）的局部平均体积分数。通过平均体积分数可以在该区域的任一点顺利地找到介于0 和1之间的值，这预示着在其他均质域中是存在少量还是大量气泡。也就是说，在分散多相流模型中，可以在同一时间和空间点上定义气相和液相；而在分离多相流模型中，在给定的时间和空间点上，只能定义气相或液相。
（a）分离多相流模型
（b）分散多相流模型
图1-20 多相流模型
对于分离多相流模拟，COMSOL软件提供了3种不同的界面追踪方法：相场法、水平集法与移动网格法。
相场法和水平集法都是基于场的方法，其中相之间的界面代表相场或水平集函数的等值面。与上述两种方法完全不同，移动网格法将相界面模拟为分隔两个域的几何表面，每个域对应不同的相。
基于场的问题通常是在固定的网格上解决，而移动网格问题要在移动的网格上解决。在相场法和水平集法中，有限元网格不必与两个相的边界一致，如图1-21所示的搅拌自由液面模拟。
（a）相场法
（b）水平集法
图1-21 相场法与水平集法
对于移动网格，网格与相边界的形状保持一致，并且网格边缘与相边界重合，如图1-22所示的搅拌自由液面模拟，是在单相流中用移动网格模拟自由液面。但是，移动网格模型也有缺点，即目前无法处理拓扑变化（例如界面分离等），而相场法与水平集法不存在这个缺点，可以处理相边界形状的任何变化。
图1-22 移动网格法
对于给定的网格，移动网格法具有更高的精度。基于这一优势，我们可以直接在相边界上施加力和通量。基于相场的方法需要围绕相边界表面建立密集网格，以解析该表面的等值面。由于很难定义一个精确贴合等值面的自适应网格，因此通常必须在等值面周围建立大量密集网格。在具有相同精度的情况下，与移动网格相比，这样做会降低基于场的方法的效率。
对于不希望发生拓扑变化的微流体系统，通常首选移动网格法；如果需要拓扑变化，则必须使用相场法。如果表面张力的影响较大，则首选相场法；如果可以忽略表面张力，则首选水平集法。
在湍流模型中，由于仅解析平均速度和压力，流体的细节会丢失。从这一点来看，表面张力效应在流体的宏观描述中也变得不那么重要。由于湍流表面的流动比较剧烈，几乎不可能避免拓扑变化，因此对于湍流模型和分离多相流模型的组合，最好使用水平集法。
在COMSOL软件中，所有湍流模型都可以与相场法和水平集法相结合来模拟两相流，例如将水平集法与k-ε湍流模型相结合来模拟反应堆中水和空气的两相流，如图1-23所示。
图1-23 反应堆中水和空气的两相流
如果多相流的相边界过于复杂而无法解析，则必须使用分散多相流模型。COMSOL软件的CFD模块提供了4种不同的分散多相流模型（原理上）：①气泡流模型，适合高密度相中包含较小体积分数的低密度相；②混合物模型，适合连续相中包含较小体积分数的分散相（或几个分散相），其密度与一个或多个分散相相近；③欧拉-欧拉模型，适用于任何类型的多相流，可以处理气体中有密集颗粒的多相流，例如流化床；④欧拉-拉格朗日模型，适合包含相对较少（成千上万，而不是数十亿）的气泡、液滴或悬浮颗粒流体，也适合气泡、颗粒、液滴或使用方程模拟的颗粒，该方程假定流体中每个颗粒的力平衡。
（1）气泡流模型显然适用于液体中的气泡。由于忽略了分散相的动量贡献，因此该模型仅在分散相的密度比连续相小几个数量级时才有效。
（2）混合物模型与气泡流模型相似，但考虑了分散相的动量贡献，通常用于模拟分散在液相中的气泡或固体颗粒。混合物模型还可以处理任意数量的分散相。混合物模型和气泡流模型均假设分散相与连续相处于平衡状态，即分散相不能相对于连续相加速。因此，混合物模型无法处理分散在气体中的大固体颗粒。
当多相流混合物被迫通过孔口时，用混合物模型模拟了5种不同大小的气泡，流动中的剪切力导致较大的气泡破裂成较小的气泡，如图1-24所示。
图1-24 混合物模型
（3）欧拉-欧拉模型是最精确的分散多相流模型，也是用途最多的分散多相流模型。该模型可以处理任何类型的分散多相流，允许分散相加速，对不同相的体积分数也没有限制，但是它为每个相定义了一组Navier-Stokes方程。
在实践中，欧拉-欧拉模型仅适用于两相流，并且其计算成本（CPU时间和内存）较高。正因如此，该模型使用起来也相对困难，并且需要良好的初始条件才能在数值解中收敛。使用欧拉-欧拉多相流模型模拟流化床中固体颗粒的体积分数分布如图1-25所示。
图1-25 欧拉-欧拉模型的模拟
（4）如果连续流体中悬浮有一些（成千上万，但不是数十亿）非常小的气泡、液滴或颗粒，则可以使用欧拉-拉格朗日模型模拟多相流系统。该方法的优点是计算成本相对较低。从数值的角度来看，这个方法通常也不错。因此，如果连续流体中分散相的颗粒数量相对较少，那么优选欧拉-拉格朗日模型，其模拟效果如图1-26所示。
图1-26 欧拉-拉格朗日模型的模拟
还有一些方法可以使用欧拉-拉格朗日模型来模拟大量粒子，它们使用的相互作用项和体积分数可以模拟具有数十亿个粒子的系统。这些方法可以在COMSOL软件中实现，但在预定义的物理接口中无法实现，需要自定义或多模块组合来实现，如在COMSOL软件中用附加的CFD模块和粒子追踪模块可实现欧拉-拉格朗日多相流模型。
混合物模型能够处理任何相的组合，并且计算成本较低。在大多数情况下，我们可以使用此模型模拟。对于流化床（具有高密度和高体积分数的大颗粒分散相）之类的系统，只能使用欧拉-欧拉模型模拟。
各种分散多相流模型本质上是近似的，并且也与近似的湍流模型非常吻合。可以在分散相和连续相之间以及在分散相中的气泡、液滴和颗粒之间引入相互作用。这些相互作用的起源可以是用湍流模型模拟的湍流。气泡流、混合物流和欧拉-拉格朗日多相流模型可以与 COMSOL软件中的所有湍流模型结合使用。
水平集法是Osher和Sethian于1988年提出的一种用于界面追踪的数值方法。在水平集法中，界面被看作零水平集的光滑函数。由于水平集函数的对流本身很光滑，因此可以代替界面处对流引起的物性陡变梯度。虽然水平集法不像其他某些方法拥有守恒属性，但它的优势在于能够轻易计算界面曲率。水平集法采用连续逼近方法，将表面张力和交界面局部曲率表示为体积力，这简化了在计算中捕捉由表面张力变化引起的拓扑结构变化过程。
在水平集法中，用水平集平滑函数来描述两相交界面。在连续相中水平集函数始终为正，在分散相中始终为负，而相间交界表面是由水平集函数为零的点构成，即
（1-103）
从上面的描述可知，交界面上的单位法线由分散相指向连续相，交界面的曲率可以用水平集函数表示为
（1-104）
（1-105）
交界面的运动可以通过水平集函数的对流来捕获：
（1-106）
流速和压力的控制方程可以用不可压缩N-S方程表示：
（1-107）
（1-108）
其中，F是体积力，包括重力和由于对交界面应力进行水平集处理引入的表面张力项。F的两个分量可以表示为
（1-109）
（1-110）
其中，为张力系数，
为界面曲率，
函数用来处理交界面处的表面张力项，可以有很多个液-液交界面来标定分散相。描述物性急剧变化的Heaviside函数可以用水平集函数表示为
（1-111）
计算区域流体的密度与黏度可以表示为
（1-112）
（1-113）
其中，
为第一相的密度，
为第二相的密度，
为第一相的黏度，
为第二相的黏度。
在COMSOL软件微流动的不混溶的两相流模拟中，水平集描述的两相流界面输运满足如下方程
（1-114）
其中，第三项为保持数值稳定性所需的项；
为界面厚度控制参数，大多数情况下使用默认值
，
为流域典型网格大小；
为水平集函数重新初始化或稳定性的参数，通常情况下是采用流体流动的最大速度。太小可能会使相界面厚度不再保持恒定，也可能会由于数值不稳定性引起的振荡，而太大会导致相界面的移动不正确。
相场法是基于Cahn-Hilliard和Ginzburg-Landau方程的改进数值方法，其相场变量由Cahn-Hilliard扩散方程决定，可以用以下两个二阶偏微分方程表示。
（1-115）
（1-116）
其中，为相场变量，取值为（‒1, 1）；λ为混合能密度；ε为界面厚度控制参数，用于评价相界面的厚度，λ与ε两个参数与表面张力系数相关，有
；γ 为迁移率，γ 足够大时，可使界面厚度保持恒定，γ 足够小时，将使得对流项不会被过分抑制，γ=χε2；χ为迁移调节参数；ψ为相场助变量；
为外部自由能（多数情况下为0）。
相场法求解时，主场方程中流体属性的控制方程为
（1-117）
（1-118）
（1-119）
相场法中的参数设置：迁移调节参数
的默认值为1，对于大多数模型来说是一个不错的初始值；界面厚度控制参数大多数情况下使用默认的，为流域典型网格大小。
混合物模型是一种宏观两相流模型，在许多方面类似于气泡流模型。该模型跟踪平均相浓度或体积分数，并求解混合物速度的单个动量方程，适用于由浸没在液体中的固体颗粒或液滴组成的混合物。
在混合物模型中，颗粒与流体组合被视为具有宏观特性（如密度与黏度）的单个连续流动体，一般是由分散相与连续相组成的两相流。例如，连续相为液体，分散相为固体颗粒、液滴或气泡。然而对于液体中的气泡，气泡流模型更合适些。混合物模型依赖于以下假设：①每相的密度是近似常数；②两相共享相同的压力场；③颗粒的松弛时间少于宏观流动的时间尺度。
混合物模型可以求解连续性方程、动量方程、分散相体积分数输运方程、质量传输方程、湍流方程、滑移速度方程等，这里简单介绍一下前面3个方程，其他的请查阅参考文献[16]。
混合物模型的连续性方程为
（1-120）
其中，是混合密度，计算式为
（1-121）
ս是混合速度，计算式为
（1-122）
其中，
分别是连续相、分散相的体积分数；
分别是连续相、分散相的密度；
分别是连续相、分散相的速度。
混合物模型的动量方程可以通过对所有相各自的动量方程求和来获得，可表示为
（1-123）
其中，u为速度矢量；为压力；为减少的密度差，
；
为两相的滑移速度矢量；
为滑移通量，
；
为黏性与湍流应力之和；
为湍流耗散系数；
为分散相到连续相的质量传输率；F为任意外部体积力。
分散相体积分数输运方程为
（1-124）
欧拉-欧拉模型界面基于在体积上平均每个当前相的Navier-Stokes方程，该体积与计算区域相比较小，但与分散相（颗粒、液滴或气泡）相比较大。
欧拉-欧拉模型包含连续性方程、动量方程、黏度方程、相间动量输送方程、固体压力方程等，这里简单介绍其中的一些方程，详细内容请查阅参考文献[16]。
欧拉-欧拉连续性方程适用于连续相与分散相：
（1-125）
（1-126）
其中，是相体积分数，
，表示密度，u为速度，各项中的下标c和d分别表示连续相与分散相；为分散相到连续相的质量传输率。
连续相和分散相的动量方程使用非保守形式，即
（1-127）
（1-128）
假设上述方程中的流体相为牛顿流体，黏性应力张量定义为
（1-129）
（1-130）
其中，是混合流体压力，
是每相的黏性应力张量，g为重力加速度，
为相间动量传递项（一相被其他相施加的体积力），F是任何其他体积力，
是相间速度，
为动力黏度，各项中的下标c和d分别表示连续相与分散相。
欧拉-欧拉模型中使用混合物黏度的表达式，两个互穿相的动力黏度默认值为
（1-131）
简单的混合物黏度（能覆盖整个颗粒浓度范围）方程可以用Krieger方程表示为
（1-132）
其中，
为最大填充限制，对于固体颗粒，其默认值为0.62。
双流体欧拉-欧拉模型是两相流体流动的一般宏观模型，它将两相视为互穿介质，跟踪相的平均浓度。其中的速度场与每个相场相关联，动量方程和连续性方程描述每个相的动力学过程。气泡流模型是双流体模型的简化，它依赖于以下假设：①与液体密度相比，气体密度可以忽略不计；②气泡相对于液体的运动由黏性阻力和压力之间的平衡决定；③两个阶段共享相同的压力场。
基于这些假设，列出两相流的动量方程和连续性方程，结合气相输运方程，就可以跟踪气泡的体积分数。
（1-133）
（1-134）
上面的方程中，是相体积分数，表示密度，u为速度，p是压力，
是每相的黏性应力张量，
为重力加速度，F是任何额外体积力，
是液体的动力黏度，
为湍流黏度，各项中的下标l和g分别表示液相和气相。
（1-135）
其中，mgl为质气体到液体的质量传输。
在实际应用时，上述模型可能根据实际情况进行变化，限于篇幅，这里不一一说明，详细内容请查阅参考文献[16]。
对于层流两相流，当关注界面的精确位置时，移动网格模型可用于模拟两种不同的不混溶流体的流动。界面位置由移动网格跟踪，边界条件考虑了表面张力和润湿以及界面上的质量传输。两相流移动网格界面是单相流界面和移动网格界面之间的预定义物理界面耦合。在对应于各个相的区域内，流体流动使用Navier-Stokes方程求解。
我们通过流体流动域内的网格变形，来说明两种流体之间的界面。软件会扰动网格结点，使其与移动界面以及模型中的其他移动或静止边界一致。边界位移在整个区域中传播，以获得平滑的网格变形。这是通过求解网格位移方程（拉普拉斯方程、温斯洛方程或超弹性平滑方程）实现的。通常情况下，在“移动边界平滑”选项中根据式（1-136）平滑法向网格速度。
（1-136）
其中，X为x坐标变化量，n为单位矢量，
为理想的法向网格速度，
是平滑速度，
，
是移动边界平滑调整参数（无量纲），
是网格尺寸（单位：m），
是平均表面曲率（单位：1/m）。
以二维为例，变形网格中的一个位置坐标（x，y）可以与其在原始未变形网格中的坐标（X，Y）相关，写成函数形式为
（1-137）
原始未变形的网格称为材料框架（或参考框架），而变形网格称为空间框架。COMSOL还定义了几何体和网格框架，这些框架与该物理界面的材料框架一致。流体流动方程（以及其他耦合方程，如电场或化学物质传输方程）在网格被扰动的空间框架中求解。因此，在这些界面中考虑了相边界的移动。
用网格位移和流体流动的特定边界条件跟踪两相之间的界面。有两个选项可用：自由表面和流体-流体界面。当外部流体的黏度与内部流体的黏度相比可以忽略时，自由表面边界条件是合适的。在这种情况下，外部流体的压力是建模流体所需的唯一参数，并且在外部流体中不求解流动。对于流体-流体界面，对两个相的流动进行求解。
两种不混溶流体（流体1和流体2）界面处的边界条件为
（1-138）
（1-139）
（1-140）
其中，
和
分别是流体1和流体2的速度，
是两流体间界面网格速度，
为界面法向（见图1-27），
和
分别为区域1与区域2的总应力张量，
为由界面张力引起的单位面积力，
，
是表面梯度算子，
是穿过界面的质量通量。
图1-27 流体1与流体2界面法向的定义
式（1-139）的切向分量在边界处的流体之间施加无滑移条件。在没有穿过边界的传质的情况下，式（1-138）和式（1-140）确保垂直于边界的流体速度等于界面的速度。当发生传质时，这些方程是质量守恒的结果，在边界静止的框架中很容易导出。
总应力张量的分量
表示垂直于
方向的每单位面积力的第
分量。因此，
（使用求和约定）被解释为作用在边界上的每单位面积的力，通常这不是边界的法线。因此，式（1-139）表示了两种流体之间界面上的力平衡。
通常情况下，流体1的黏度显著大于流体2的黏度（例如对于气液界面）。在这种情况下，流体2的总应力中的黏度项可以忽略，式（1-139）变为
（1-141）
外部流体（流体2）仅通过压力项进入方程系统，并且系统可以用仅由流体1组成的域表示，该域在流体2的域中具有外部压力
的表达式（或恒定值）。
1.5 其他物理场分析模型传热一般包括传导、对流、辐射3种方式，相应地有3种传热问题。固体的传热方程为
（1-142）
其中，
是恒定应力下的比热容，表示密度，
为平移运动速度，
是绝对温度， q是传导热通量，qr为辐射热通量，
是热膨胀系数，S是第二Piola-Kirchhoff应力张量，
为额外的热源。
流体的传热方程为
（1-143）
其中，
是黏性应力张量，
是流体的速度，是恒定应力下的比热容，表示密度，T是绝对温度，q是传导热通量，qr为辐射热通量，
是热膨胀系数，
，p是压力，Q为黏性耗散以外的热源。
宏观层面上的电磁分析问题是在一定边界条件下求解麦克斯韦（Maxwell）方程组。麦克斯韦方程组描述了基本电磁量之间的关系，其中的主要物理量为电场强度E、电位移或电通量密度D、磁场强度H、磁通密度B、电流密度J、电荷密度。麦克斯韦方程可以用微分形式或积分形式表示，采用微分形式便于有限元法处理。对于一般时变场，麦克斯韦方程可以写成：
（1-144）
（1-145）
（1-146）
（1-147）
上述方程的前两个分别称为麦克斯韦-安培定律和法拉第定律。第三个方程和第四个方程分别是高斯定律的两种形式：电形式和磁形式。
另一个基本方程是连续性方程，如下所示。
（1-148）
上述5个方程中，只有3个是独立的。前两个方程与高斯定律的电形式或连续性方程结合可以形成一个封闭的系统。
为了获得封闭系统，需要包括描述介质宏观性质的本构关系，如下所示。
（1-149）
其中，
为真空介电常数，
为真空磁导率，为电导率，P为极化强度，M为磁化强度。在国际单位制中，真空磁导率与无量纲精细结构常数成正比，其值为4π×10−7 H/m。真空中电磁波的速度为
，约为3×108m/s，真空介电常数
，约为
。
在某些情况下，用标量电势V和矢量磁势A来描述问题有利于数值求解，表达式如下。其中矢量磁势的定义方程由磁形式的高斯定律直接给出，电势由法拉第定律产生。
（1-150）
将式（1-149）应用到安培环路定律和电形式的高斯定律中，经推导，分别得到以下磁场偏微分方程和电场偏微分方程。
（1-151）
（1-152）
要全面描述电磁问题，必须在材料界面和物理界面处指定边界条件。在两种介质之间的界面处，边界条件可以表示为
（1-153）
其中，
和
分别表示表面电荷密度和表面电流密度，
是介质2的外法线。这些条件中只有两个是独立的，这是一个超定方程组，因此需要简化。先选择方程一或方程四，然后选择方程二或方程三，这些选择一起形成一组独立的条件。根据这些关系，我们可导出电流密度的界面条件：
（1-154）
时谐场量与其相量之间的关系为
（1-155）
其中，
是一个相量，它包含场的振幅和相位信息，但与t无关。
电磁力公式可以表示为
（1-156）
通过有限元法，我们可以求解全波形式的麦克斯韦方程。假设角频率已知，为
，电磁场随时间呈正弦变化，且材料的所有属性相对于场强呈线性变化，则三维麦克斯韦控制方程可简化为
（1-157）
其中，k0表示波数，
表示相对磁导率，
表示相对介电常数，表示电导率。已知真空中光速为，则可在整个模拟域内对电场E=E(x, y, z)求解上述方程，其中E为矢量，可用其分量表示为
。其他诸如磁场强度、功率、电流等物理量都可从电场推导出。
在使用“电磁波，波束包络”接口时，我们可以从“电磁波，波束包络”设置窗口中查看该接口的控制方程：
（1-158）
其中，
是包络函数，为求解的因变量。
在场的相量表示中，对应于振幅，
代表相，即
（1-159）
式（1-158）为波束包络接口的控制方程，可以通过将式（1-159）代入亥姆霍兹方程（电磁波方程，
，k为波数，A为振幅）中导出。假设已知，是唯一未知量，这样就可以求解，因此，在使用此方法时，需要提前知道波矢或相函数。
光波也是电磁波谱的一种，在使用“几何光学”接口时，可以从“几何光学”设置窗口中查看该接口的内置方程：
（1-160）
其中，q为光线位置，k为波矢。
电磁场分析理论涉及电磁学的较多知识，上述仅简单介绍电磁学的基础知识，详细内容可查阅参考文献[8]。
多孔材料由固体结构（多孔基质）和填充有液体或气体的孔隙（空洞）组成。多孔材料有各种尺寸和广泛的应用——从纳米材料到多孔反应器，从电子元件的冷却到大规模的岩土工程应用。它们的共同点是，材料的总尺寸远大于平均孔径，因此必须使用宏观方法建立模型。
通常用孔隙率和渗透率两个参数表征多孔材料。孔隙率
描述了孔隙或空洞体积与总体积之比，
。渗透率
表征了流体通过多孔材料的能力。
描述多孔材料中液体流动的基本定律是达西定律。它描述了速度场u（m/s）和压力梯度p（Pa）之间的线性关系，此式仅用于速度很低（Re < 10）的情况。
（1-161）
在流速相对较快（Re>10）或克努森数相对较高（Kn>0.1）的情况下，达西定律不再有效。因此，引入了不同的渗透率模型来捕捉这些影响。
压力梯度与速度的非线性关系的一般形式可以写成：
（1-162）
其中，β是取决于多孔介质特性的常数。
通过填充床的流态可由床的雷诺数确定。通常情况下，对于雷诺数Re < 10的情况，可以用科泽尼-卡曼方程（达西流）来描述流动，
，
为颗粒平均直径。对于10 < Re < 1000（有时称为过渡区）的情况，流动由埃尔根方程更好地描述，
。对于Re > 1000的情况，埃尔根方程可由湍流的伯克-普卢默方程近似，
，L为填充床的长度。
（1-163）
其中，
是多孔介质每单位体积的质量源（不是每单位孔隙体积）。
（1-164）
其中，K是黏性应力张量。
标准声学问题涉及求解固定背景压力
之上的小声压变化。从数学角度来讲，这代表了围绕固定静态值的线性化（小参数扩展）。
通过对动量方程（欧拉方程）和连续性方程的变换，可以得到无损介质中声波的波动方程。
（1-165）
其中，
是体积模量，表示密度，是压力，
为偶极子声源，为单极子声源。
这里用简短的实例说明反应工程和化学程序中如何处理物质平衡方程中的平衡反应。对于组分A与B的如下反应：
（1-166）
组分A到B的反应率为
（1-167）
其中，
和
分别是A和B的浓度，kf和kr分别是正向和反向速率常数。
假定组分A与B的反应是平衡的，则反应率为
，也有
，
。反应工程程序能够定义平衡系统的质量平衡，而无须反应速率表达式。求解的方程组如下。
（1-168）
（1-169）
其中，
表示正向和反向反应速率之间的关系。
一般而言，对于贡献k个质量平衡且j个反应处于平衡的反应系统，要求解的简化方程组由k–j个质量平衡式和j个反应平衡式组成。COMSOL生成上述方程组的消除过程是自动化的，允许对化学平衡反应以及不可逆或可逆反应进行简单建模。
COMSOL软件涉及电化学的3个物理场分别是“一次电流分布”“二次电流分布”和“三次电流分布，Nernst，Planck”。
在分析“一次电流分布”时，忽略了电极动力学和浓度依赖性效应造成的损耗，假定电解液中的电荷转移遵守欧姆定律，仅考虑几何因素的影响，其控制方程为
（1-170）
（1-171）
（1-172）
其中，i为电流密度矢量，Q为一般电流源项，
为电势，为电导率，各项中的下标l代表电解质，s代表电极。
为反应的平衡电位。
在“二次电流分布”分析中，忽略了浓度极化时的电流分布，考虑了电极动力学的影响，也假定电解液中的电荷转移遵守欧姆定律。它的域方程为式（1-170）和式（1-171），而电极和电解质界面上的电位方程为
（1-173）
其中，
为活化过电位，即实际电位差和平衡电位的差值。
在“三次电流分布，Nernst，Planck”分析中，考虑电解质组成和离子强度的变化对电化学过程的影响，以及溶液电阻和电极动力学的影响，利用Nernst-Planck方程来描述电解质中化学物质的传递，不再假定电解液中的电荷转移遵守欧姆定律。该方法考虑的因素较多，会导致模型过于复杂、求解时间变长。“三次电流分布，Nernst，Planck”方法包括“三次分布，电中性”“三次分布，水基电中性”“三次分布，支持电解质”三个子方法。下面的方程为“三次分布，电中性”子方法的控制方程，其余两个子方法的控制方程可查阅参考文献[10]或从软件系统的具体设置窗口中查看。
（1-174）
（1-175）
（1-176）
（1-177）
练习题1．简述有限元法的基本思想及其分析的流程。
2．说明有限元建模的准则。
3．简述减小有限元计算规模的措施。
4．图1-28所示的矩形薄板受表面拉力作用，板厚t = 20mm，E = 210GPa，μ = 0.3，试用有限元法确定节点位移和单元应力（划分单元数量自定）。
图1-28 矩形薄板受表面拉力作用
5．边长为2a的正方形薄板如图1-29所示，它的厚度为t，两侧边固定，上边受均布载荷q作用，试利用有限元法对该薄板进行应力分析。
图1-29 正方形薄板的受力与约束
6．图1-30所示的区域内有牛顿流体的流动，流体黏度为10Pa·s，密度为1100kg/m3。流动区域为100mm×40mm的矩形区域，入口给定压力1000Pa，出口敞开。在考虑惯性项影响的情况下，计算区域内的速度-压力分布，并分析惯性项对出口流量的影响。
图1-30 带有小入口和出口的矩形计算区域
7．在图1-31所示的二维物体中，左边的温度保持在40℃，顶部和底部绝热，右边有对流，对流系数
，自由流温度T∞=10℃，导热系数
，板厚1m，其余尺寸均标注在图中，试求该二维物体的温度分布。
图1-31 有温度变化和对流的二维物体
