2.3　数学模型的求解方法

2.3.1　有限差分法

2.3.1.1　有限差分法简介

有限差分法（FDM）是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法通过Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。

有限差分法在材料成形领域的应用较为普遍，是材料成形计算机模拟技术领域中最主要的数值分析方法之一。目前材料加工中的传热分析（如铸造过程中的传热凝固、塑性成形过程中的传热、焊接过程中的传热等）和流动分析（如铸件成型过程，焊接熔池的产生、移动等），都可以用有限差分法进行模拟。与有限元法相比，有限差分法在流场分析方面优势明显。

对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑，可分为中心格式和逆风格式。考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，构造差分的方法有多种形式，包括Taylor级数展开法、多项式拟合法、控制容积积分法和平衡法，目前主要采用的是Taylor级数展开方法。其基本的差分表达式主要有四种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。

2.3.1.2　有限差分法数学基础

设有自变量为x的解析函数y=f（x），则y对x的导数为：

　　        （2-34）

式中，dy、dx分别为函数及自变量的微分；为函数对自变量的导数，又称微商；Δy、Δx分别为函数及自变量的差分；为函数对自变量的差商。

在导数的定义中，Δx是以任意方式趋近于0的，即Δx是可正可负的，而在差分方法中，Δx总为正数，则差分可以有以下三种形式。

向前差分：

Δy=f（x+Δx）-f（x）　　        （2-35）

向后差分：

Δy=f（x）-f（x-Δx）　　        （2-36）

中心差分：

　　        （2-37）

上述对应于一阶导数的差分称为一阶差分，相应地，对应于二阶导数的差分称为二阶差分，二阶差分是在一阶差分的基础上再作一阶差分，标记为Δ²y，二阶差分的三种形式分别表示为如下形式。

向前差分：

Δ²y=Δ（Δy）=Δ［f（x+Δx）-f（x）］=Δf（x+Δx）-Δf（x）

=［f（x+2Δx）-f（x+Δx）］-［f（x+Δx）-f（x）］

=f（x+2Δx）-2f（x+Δx）+f（x）　　        （2-38）

向后差分：

Δ²y=Δ（Δy）=Δ［f（x）-f（x-Δx）］=Δf（x）-Δf（x-Δx）

=［f（x）-f（x-Δx）］-［f（x-Δx）-f（x-2Δx）］

=f（x-2Δx）-2f（x-Δx）+f（x）　　        （2-39）

中心差分：

　　        （2-40）

任何阶差分都可由比其低一阶的差分再作一阶差分得到。例如n阶向前差分为：

Δⁿy=Δ（Δ^n-1y）=Δ［Δ（Δ^n-2y）］=…=Δ{Δ…［Δ（Δy）］}
=Δ{Δ…［（Δf（x+Δx）-f（x））］}　　        （2-41）

相应地，n阶的向后差分以及中心差分分别表示如下。

向后差分：

Δⁿy=Δ{Δ…［Δ（f（x）-f（x-Δx））］}　　        （2-42）

中心差分：

　　        （2-43）

函数的差分与自变量的差分之比，称为函数对自变量的差商。一阶差商的三种形式分别如下。

向前差商：

　　        （2-44）

向后差商：

　　        （2-45）

中心差商：

　　        （2-46）

或

　　        （2-47）

相应地，二阶差商的三种形式如下。

向前差商：

　　        （2-48）

向后差商：

　　        （2-49）

中心差商：

　　        （2-50）

上述均为一元函数的差分及差商，多元函数的差分与差商也可以用类似的方法得到，如对于多元函数的一阶向前差分为：

　　        （2-51）

　　        （2-52）

有限差分法的数学基础是用差分代替微分，用差商代替微商，而用差商代替微商的几何意义是用函数在某区域内的平均变化率来代替函数的真实变化率。对于一阶微商，存在以下三种典型的差分形式。

向前差商：

　　        （2-53）

向后差商：

　　        （2-54）

中心差商：

　　        （2-55）

根据泰勒级数式，可以计算出上述三种差分形式的误差，分别为：

　　        （2-56）

　　        （2-57）

　　        （2-58）

从以上三式可以看出，用不同的方法定义的差商代替微商所引起的误差是不同的。用向前差商或向后差商代替微商，其截断误差是O（Δx），是Δx一次方的数量级；用中心差商代替微商，其截断误差是O（Δx）²，是Δx二次方的数量级，即用中心差商代替微商比用向前差商或向后差商代替微商的误差小一个数量级。

同样，对于二阶差商，其差分形式一般采用中心式：

　　        （2-59）

其截断误差为：

　　        （2-60）

从上面的分析可以看出，用差商代替微商必然会带来截断误差。相应地，用差分方程代替微分方程也会带来误差，因此，在应用有限差分法进行计算的时候，必须注意差分方程的形式、建立方法及由此产生的误差。

2.3.1.3　有限差分法解题基本步骤

有限差分法的主要解题步骤如下。

（1）建立微分方程　根据问题的性质选择计算区域，建立微分方程式，写出初始条件和边界条件。

（2）构建差分格式　首先对求解区域进行离散化，确定计算节点，选择网格布局、差分形式和步长；然后以有限差分代替无限微分，以差商代替微商，以差分方程代替微分方程及边界条件。

（3）求解差分方程　差分方程通常是一组数较多的线性代数方程，其求解方法主要包括两种：精确法和近似法。其中精确法又称直接法，主要包括矩阵法、Gauss消元法及主元素消元法等；近似法又称间接法，以迭代法为主，主要包括直接迭代法、间接迭代法以及超松弛迭代法。

（4）精度分析和检验　对所得到的数值解进行精度与收敛性分析和检验。

2.3.1.4　有限差分法解题示例

【例2-5】　设有一炉墙，厚度为δ，炉墙的内壁温度T₀=900℃，外壁温度T_m=100℃，求炉墙沿厚度方向上的温度分布。

【解】　这是一个一维稳态热传导问题，其边界条件为T₀=900℃、T_m=100℃，可以用有限差分方法求得沿炉墙厚度方向上的若干个节点的温度值。

（1）建立微分方程。根据热力学知识，对于常物性、一维、稳态热传导的微分方程为：

　　        （2-61）

（2）构建差分格式。首先确定计算区域并将其离散化。对于稳态热传导问题，只需将空间离散化。如图2-3所示，把需求解的空间区域0～δ以某一定间距划分为m等份，这些等分线成为网格线。以每一网格线为中心，取宽度为1组成一系列的子区间，成为单元体（图中阴影部分）。单元体的中心点成为节点，节点依次标记为0，1，…，m。在计算过程中，将节点的温度作为单元体的平均温度，如将节点i温度作为单元体i的平均温度，记为T_i；边界节点温度则为半个单元体的平均温度，记为T₀和T_m。在此计算区域内构件差分格式，根据式（2-58），可得：

　　        （2-62）

当m=4时建立差分方程如下：

T₀=900

T₂-2T₁+T₀=0

T₃-2T₂+T₁=0

T₄-2T₃+T₂=0

T₄=100

（3）求解差分方程。利用Gauss消元法可解除上述的线性方程组，得到炉墙特定点的温度分布，见表2-4。

（4）求解结果分析与检验。根据热力学知识可知，炉墙的温度分布应与其厚度呈线性变化关系，求解结果符合之一规律。同时，通过解析法求解微分方程（2-47），得到的解析解为：，将、、分别代入后可得到相应的温度值为700℃、500℃和300℃，这与表2-4中的计算结果是一致的。

【例2-6】　利用差分法解Laplace方程第一边值问题（要求画出差分网络及写出差分方程组）。

　　        （2-63）

【解】　采用正方形网格剖分，内节点按如图2-4所示编号。设内节点总数为N，对于每一个（x_i，y_j）∈D₀，利用数值微分公式得

　　        （2-64）

　　        （2-65）

式中，h₁、h₂分别为沿x、y轴方向的步长。

将式（2-64）和式（2-65）代入式（2-63）得

　　        （2-66）

其中

　　        （2-67）

为除去R_ij项后所剩下的项之和，略去余项R_ij，得

　　        （2-68）

取h₁=h₂，即取正方形网格时，差分方程为

4h_ij-［u_（i+1）j+u_（i-1）j+u_i（j+1）+u_i（j-1）］=-h²f_ij　　（x_i，y_j）∈D₀　　        （2-69）

本例中取h₁=h₂=0.125，采用正方形网格剖分，内节点按图2-4所示编号，式（2-69）得

i=1，j=1，4u₁₁-（u₂₁+u₀₁+u₁₂+u₁₀）=0

其中，u₁₁对应u₁，u₂₁对应u₂，u₀₁为0，u₁₂对应u₄，u₁₀为0，于是得4u₁-u₂-u₄=0。其余类推得差分方程：

　　        （2-70）

写成矩阵形式为：

　　        （2-71）

用Seidel迭代法求得u₁=6.25，u₂=12.5，u₃=18.75，u₄=12.50，u₅=25，u₆=37.50，u₇=18.75，u₈=37.5，u₉=56.25。

2.3.2　有限元法

有限元法（又称为有限单元法、有限元素法）是20世纪50年代初才出现的一种新的数值分析方法，最初它只应用于力学领域中，70年以来被应用于传热学计算中。与有限差分法相比较，有限元法的准确性和稳定性都比较好，且由于其单元的灵活性，使它更适应于数值求解非线性热传导问题以及具有不规则几何形状与边界，特别是要求同时得到热应力场的各种复杂导热问题。有限元法在传热学中的应用正处于开拓与发展阶段，迄今为止，其应用已波及热传导、对流传热及换热器设计与计算。

有限元法是变分法与经典有限差分法相结合的产物，它既吸收了古典变分近似解析解法——一种泛函求极值的基本原理，又采用了有限差分的离散化处理方法，突出了单元的作用及各单元的相互影响，形成了自身的独特风格。

古典变分法是要寻求定解问题的级数形式近似解析解，在这种方法中，首先构造一个与定解问题（微分方程及其边值条件）相对应的泛函，然后对此泛函求极值，从而得到满足微分方程和边值条件的近似解析解。这样一来，就把选择泛函并对泛函求极值的运算，等价于一个在数学上对微分方程及其边值条件所组成的定解问题的求解。由于这种方法首先在弹性力学中得到应用，而在弹性力学中是以最小能位原理为平衡条件加以分析的，故曾将上述泛函求极值的数学概念与最小能位原理的物理概念联系起来，从而称上述变分法为能量法（又称Ritz法）当然，这种变分法绝不只是最小能位原理的表述，得此命名只不过反映了这种方法的实用背景和历史发展过程而已。但是，遗憾的是并非所有定解问题均可找到其相对应的泛函，有些定解问题可能根本不存在其对应的泛函。

于是，人们设法直接从微分方程出发去寻找其近似级数解，从而回避了寻找泛函这一难题，这种解法就是加权余量法，其中Galerkin法是较典型的一种。由于这种方法与能量法颇相似，也要选择适当的函数代入微分方程，然后对其加权积分使其为零，故在广义上亦称为变分法。无论采用能量法，还是加权余量法，都要选择与微分方程相对应的适当的函数代入泛函或代入微分方程，再对泛函求极值或对微分方程加权积分使其为零，这种函数称为试探函数。对此函数，要求其在全区域内满足定解问题，这一要求是极苛刻的，不能适应许多工程实际中的复杂热传导问题，这就使古典变分法的应用受到了很大的限制。有限元法是对上述古典变分法的改进，也就是采取与有限差分相类似的方法，将区域离散化，以只在离散化有限小的单元内使试探函数满足定解问题要求并在单元内积分，代替在全区域内满足要求与积分的条件，消除了古典变分法的局限性。在这层意义上来说，有限元法就是有限的单元变分法。

2.3.2.1　有限元法常用术语

（1）单元　有限元模型中每一个小的块体称为一个单元。根据其形状的不同，可以将单元划分为以下几种类型：线段单元、三角形单元、四边形单元、四面体单元和六面体单元等。由于单元是构成有限元模型的基础，因此单元类型对于有限元分析至关重要。一个有限元软件提供的单元种类越多，该程序功能就越强大。

（2）节点　用于确定单元形状、表述单元特征及连接相邻单元的点称为节点。节点是有限元模型中的最小构成元素。多个单元可以共用一个节点，节点起连接单元和实现数据传递的作用。

（3）荷载　工程结构所受到的外部施加的力或力矩称为荷载，包括集中力、力矩及分布力等。在不同的学科中，荷载的含义有所差别。在通常结构分析过程中，荷载为力、位移等；在温度场分析过程中，荷载是指温度等；而在电磁场分析过程中，荷载是指结构所受的电场和磁场作用。

（4）边界条件　边界条件是指结构在边界上所受到的外加约束。在有限元分析过程中，施加正确的边界条件是获得正确的分析结果和较高的分析精度的关键。

（5）初始条件　初始条件是结构响应前所施加的初始速度、初始温度及预应力等。

2.3.2.2　有限元法数学基础

（1）微分方程的等效积分形式　工程或物理学中的许多问题，通常是以未知场函数应满足的微分方程和边界条件的形式提出来的，一般可以表示为未知函数u。未知函数u应满足微分方程组：

　　        （2-72）

同时未知函数u还应满足边界条件：

　　        （2-73）

S是求解域V的边界，分为S_u、S_p两部分，如图2-5所示。

未知函数u可以是标量场或向量场。A、B表示相对于独立变量（空间坐标、时间坐标等）的微分算子。微分方程数应和未知场函数的数目相等。

由于式（2-72）和式（2-73）分别在域V内和边界S上的每一点都必须为0，所以对于任意的函数向量I和，下式恒成立：

　　        （2-74）

同样的，若要求对于所有的函数向量I和，上式恒成立，则式（2-72）和式（2-73）必须满足。因此，式（2-74）称为微分方程的等效积分形式。

（2）加权余量法　对于式（2-72）和式（2-73）所表达的物理问题，往往难以求得场函数u的精确解，此时可采用以下的近似函数来便是位置的场函数：

　　        （2-75）

其中，N为试探函数，又称形函数，a为待定参量。

通常在式（2-75）取得有限项的条件下，近似解不能精确满足式（2-72）和全部边界条件式（2-73），而会产生余量R和。

　　        （2-76）

这里用某个给定的函数W_j和来代替等效积分式（2-74）中的I和，得到：

　　        （2-77）

上式的意义是通过选择待定参量a，强迫余量在某种意义上等于0。其中W_j和成为加权函数。令式（2-76）等于0，得到一组求解方程，用以求解待定参量a，从而可以得到原问题的近似解。采用使余量的加权积分为0来求得微分方程近似解得方法成为加权余量法。任何独立的完全函数集都可以用作权函数。

（3）伽辽金法　伽辽金法是加权余量法的一种。该方法将原来的形函数作为加权函数，即：W_j=N_j，在边界S上，=-W_j=-N_j，于是式（2-76）可写为：

　　        （2-78）

伽辽金法求解方程的系数矩阵往往是对称的，因此在使用加权余量法建立有限元模型时，大多采用伽辽金法。

2.3.2.3　有限元分析基本步骤

有限元分析的基本步骤如下。

①建立求解域并将其离散化为有限单元，即将连续体问题分解成节点和单元等个体问题；

②假设代表单元物理行为的形函数，即假设代表单元解的近似连续函数；

③建立单元方程；

④构造单元整体刚度矩阵；

⑤施加边界条件、初始条件和载荷；

⑥求解线性或非线性的微分方程组，得到节点求解结果及其他重要信息。

【例2-7】　设有一炉墙，厚度为1m，炉墙的内壁温度T₀=0℃，外壁温度T_m=0℃，在炉墙内具有内热源φ，如图2-6所示。炉墙的热传导系数为1W/（m·℃），求炉墙沿厚度方向上的温度分布。

其中，

【解】　这是一个具有内热源的一维稳态热传导问题，边界条件已知，可采用有限元法求解其温度场分布。

（1）建立微分方程。根据热力学知识，对于常物性、一维、稳态热传导的微分方程为：

　　        （2-79）

（2）构建形函数并建立单元方程。选取傅里叶级数为近似解，即有：

　　        （2-80）

上式中，a_i为待定系数，sin（iπx）为形函数N_i。将边界条件代入可知，近似解满足边界条件，且在求解域中连续。则根据式（2-76）可得：

　　        （2-81）

对上式进行分部积分可得：

　　        （2-82）

在边界上，W_j=0，上式简化为：

　　        （2-83）

（3）构造单元整体刚度矩阵。将式（2-83）变换形式，则

Ka+F=0　　        （2-84）

上式中，F=［F₁ F₂ F₃ … F_n］^T

a=［a₁ a₂ a₃ … a_n］^T

（4）求解非线性方程组。取n=1，依据式（2-83）有：

采用伽辽金法，则有：

即：

求解得：

则该一维热传导问题的一项解为：

取n=2，依据式（2-79）有：

T=a₁sin（πx）+a₂sin（2πx）

采用伽辽金法，则有：

W₁=N₁=sin（πx），W₂=N₂=sin（2πx）

将上式代入式（2-82）可得：

即：

求解得：

则该一维热传导问题的二项解为：

本例的精确解为：

图2-7绘制了精确解和有限元模拟结果的曲线，从中可以看出，二项解和精确解比较接近，而一项解和精确解相差较大。可见，随着近似解的项数增加，解的精度将不断提高，但求解的工作量也随之增加。因此，在有限元模拟过程中，需要在求解精度和求解工做量之间做出合理选择。

2.3.2.4　有限元法的基本概念——直接刚度法

在刚提出有限元法的时候采用的是直接刚度法，它源于结构分析的刚度法。因刚度法只能处理一些比较简单的实际问题，现在已很少使用了，但它对我们理解和明确有限元法的一些物理概念是有帮助的。所以我们通过一个例子来介绍直接刚度法，同时说明有限元法求解的一般步骤。

【例2-8】　一个受到轴向荷载的变截面杆如图2-8所示。杆的一端固定，另一端承受P=1000N的载荷，杆的顶部宽w₁=2cm，底部宽w₂=1cm，厚度t=0.125cm，长度L=10cm，杆的弹性模量E=10.4×10⁶MPa。试分析该杆沿长度方向不同位置的变形情况，假设杆的质量可以忽略不计。

【解】（1）前处理过程

①将求解域离散化，先将求解的问题分解为节点和单元。为简单起见，将杆划分成五个节点和四个单元，如图2-9所示。增加节点和单元的数目，将提高近似解的精度。给定的变截面杆简化为四个独立的部分，每部分的截面面积恒定。每个单元的截面积为组成该单元的节点处的面积的平均值。

②求某一单元的近似解。为研究典型单元的变形情况，先考虑长度为l，有均一截面A的固体单元在受到外力F时的变形情况，如图2-10所示。

单元中的平均应力为：

　　        （2-85）

平均正应变为：

　　        （2-86）

在弹性范围内，应力和应变的关系由Hooke定律描述，其中，E为材料的弹性模量，结合式（2-84）和式（2-85）可得

　　        （2-87）

注意，式（2-86）与现行弹簧等式F=kx相似。因此可以用一弹簧的变形来模拟固态单元的变形，弹簧的等价刚度为：

　　        （2-88）

回到前面讨论的问题，杆的横截面积沿y轴方向变化，作为近似，将杆模型化为不同截面的等截面杆的串联，如图2-9所示，这样，杆可以由一个四个弹簧串联组成的模型来表示，每个单元的弹性行为可以用等价线性弹簧来表述

　　        （2-89）

其中，等价单元刚度为：

　　        （2-90）

A_i+1和A_i分别是节点i+1和i处的横截面积，l是中元的长度。应用上述模型，我们来考虑作用在每个节点上的力。图2-11示出了节点1到节点5的受力情况。

静态平衡要求作用在每个节点上的力的总和为零，因此可得到以下方程。

节点1：R₁-k₁（u₂-u₁）=0

节点2：k₁（u₂-u₁）-k₂（u₃-u₂）=0

节点3：k₂（u₃-u₂）-k₃（u₄-u₃）=0

节点4：k₃（u₄-u₃）-k₄（u₅-u₄）=0

节点5：k₄（u₅-u₄）-P=0

　　        （2-91）

将这些方程式写成矩阵的形式，有：

　　        （2-92）

在荷载矩阵中区分施加力和反作用力也是必要的。因此矩阵式又可写为：

　　        （2-93）

一般地，可写为：

R=Ku-F　　        （2-94）

即，{反作用力矩阵}=［刚度矩阵］{位移矩阵}-{荷载矩阵}

考察以上所讨论的问题，因为杆在顶端固定，节点1的位移应为零，即R₁=0。将边界条件用于式（2-93）。得到如下矩阵

　　        （2-95）

解此方程，即可得到个节点的位移。

③确定每个单元的方程。在本问题中，每个单元有两个节点，每个节点都和一个位移相关，因此对每个节点可建立两个方程，这些方程一定包含节点位移和单元的刚度。

如图2-12所示，考虑内力f_i和f_i+1以及节点位移u_i和u_i+1根据平衡条件，不管选择哪种坐标系均有f_i+f_i+1=0。考虑到向前差分的一致性，可以选用图2-12（b）的形式。这样，节点i和i+1处的传输力可由下列方程表示：

f_i=k_eq（u_i-u_i+1）

f_i+1=k_eq（u_i+1-u_i）

写成矩阵形式为：

　　        （2-96）

④集成单元。将式（2-95）用于所有单元，并进行集成，将得到总体刚度矩阵。单元1的刚度矩阵为：

　　        （2-97）

其在总体刚度矩阵中的位置为：

　　        （2-98）

在刚度矩阵旁边给出了节点位移矩阵，有助于理解节点对相邻单元的贡献。

同样，对单元2有：

　　        （2-99）

其在总体刚度矩阵中的位置为：

　　        （2-100）

对单元3有：

　　        （2-101）

其在总体刚度矩阵中的位置为：

　　        （2-102）

对单元4有：

　　        （2-103）

其在总体刚度矩阵中的位置为：

　　        （2-104）

将组装或相加，即得到总体刚度矩阵：

　　        （2-105）

可以看到从单元刚度矩阵组合得到的总体刚度矩阵与前面通过自由体积分析得到的总体刚度矩阵完全一致。

⑤施加边界条件和载荷。将杆的顶端固定，应满足边界条件u₁=0，载荷P施加在节点5上。应用这些条件得到下式：

　　        （2-106）

对于固体力学问题，有限元分析有如下一般形式：

［刚度矩阵］{位移矩阵}={荷载矩阵}

（2）求解阶段。杆的截面沿y轴的变化可以描述为：

　　（2-107）

利用式（2-107）可以计算出各节点处截面的大小为：

A=0.25cm²，A₂=0.21875cm²，A₃=0.1875cm²，A₄=0.15625cm²，A₅=0.125cm²

下面，计算每个单元的刚度系数，因为

　　        （2-108）

据此，可计算出k₁=975×10³kN/cm，k₂=845×10³kN/cm，k₃=715×10³kN/cm，k₄=585×10³kN/cm以及单元刚度矩阵为：

　　        （2-109）

　　        （2-110）

　　        （2-111）

　　        （2-112）

将单元刚度矩阵组合，得到总体刚度矩阵：

　　        （2-113）

应用边界条件并施加荷载，得到：

　　        （2-114）

解此矩阵方程，即可得到各节点的位移值：

u₁=0，u₂=0.001026cm，u₃=0.002210cm，u₄=0.003608cm，u₅=0.005317cm

（3）后处理阶段。对分析结果进行处理，从中可以得到其他一些信息。如每个单元中的平均正应力。这些值可以从下式得到：

　　        （2-115）

因为节点的位移已知，可以从应力-应变关系直接得到下式：

　　        （2-116）

应用式（2-116）可以计算出例题中各单元的平均正应力值：

σ^（1）=4268kN/cm²，σ^（2）=4925kN/cm²，σ^（3）=5816kN/cm²，σ^（4）=7109kN/cm²

此外，通过计算还可以得到一些其他信息。

2.3.2.5　有限元程序的结构和特点

有限元法的实现必须通过计算机。全部有限元法的计算原理和数值方法集中反映在有限元法的程序中，因此有限元法的程序极为重要。它应具有分析准确可靠、计算效率高、使用方便、易于扩充和修改等特点。

有限元法程序总体可分为三个组成部分：前处理部分、有限元分析本体程序、后处理部分。

有限元分析本体程序是有限元分析程序的核心，它根据离散模型的数据文件进行有限元分析。有限元分析的原理和采用的数值方法集中于此，因此它是有限元分析准确可靠的关键。选用计算方法的合理与否决定了有限元分析程序的计算效率和结果的精度及可靠性。

离散模型的数据文件主要应包括：离散模型的节点数及节点坐标、单元数及单元节点编码、载荷信息等。对于一个实际的工程问题。离散模型的数据文件十分庞大，靠人工处理和生成一般是不可能的。除工作量大不能忍受外，还会不可避免地出现数据错误，包括数据精度的不足。为了解决这一问题，有限元分析程序必须有前处理程序。前处理程序是根据使用者提供的对计算模型外形及网格要求的简单数据描述，自动或半自动地生成离散模型的数据文件，并要生成网格图供使用者检查和修改。这部分程序的功能很大程度上决定了程序使用的方便性。

同样，有限元分析程序的计算结果也是针对离散模型得到的。例如对静力平衡问题可以得到离散模型各节点的位移、各单元的应力等。输出的文本文件量很大，但却不易得到所分析对象的全貌，例如位移哪里最大、应力集中发生在什么部位以及变化趋势如何。因此，一个使用方便的有限元分析程序不仅要有可供选择输出内容的文本文件，还需有结果的图形显示，如位移图、等应力线图或截面应力分布图等。这部分程序称后处理程序，与前处理程序相似，对程序使用的方便性有举足轻重的作用。

有限元分析程序的三个组成部分对于一个较好的用于实际问题分析的有限元程序来说，前后处理的程序量常常超出有限元分析的本体程序。前后处理功能越强，程序的使用就越方便。有限元分析程序中前后处理程序一般可占全部程序条数的2/3～4/5。有的近期发展的通用程序更注重程序的“包装”和使用功能，有限元分析本体程序以外部分的比例更高。

2.3.2.6　有限元软件简介

（1）ANSYS软件　世界著名力学分析专家、匹兹堡大学教授J.Swanson创立的SASI的大型通用有限元分析软件，世界最权威的有限元产品。

（2）SAP　美国加州大学伯克利分校M. J. Wilson教授的线性静、动力学结构分析程序。

（3）IDEAS　美国SDRC公司的机械通用软件，集成化设计工程分析系统。集设计、分析、数控加工、塑料模具设计和测试数据为一体的工作站用软件。

（4）NASTRAN　美国国家航空和宇航局（NASA）的结构分析程序。

（5）ADINA　美国麻省理工学院机械工程系的自动动力增量非线性分析有限元程序。

（6）ALGOR　美国ALGOR公司在SAP5和ANINA有限元分析程序的基础上针对微机平台开发的通用有限元分析系统。

（7）ABAQUS　通用有限元分析软件。

（8）DEFORM　材料成型分析专用非线性有限元软件。

（9）AUTOFORM　薄板成形模拟软件。

（10）DYNAFORM　板料冲压成型模拟软件。

（11）MSC/MARC　非线性分析有限元软件。

（12）MSC/NASTRAN　结构分析有限元软件。

（13）VPG　专业整车分析模拟软件。

（14）SYSWELD　焊接与热处理分析软件。

（15）DIANA　荷兰TND DIANA公司开发的钢筋混凝土通用有限元分析软件。