目录

COMSOL与偏微分方程

COMSOL仿真步骤

COMSOL中的有限元方法

1. 区域的离散或划分

2. 插值函数的选择

3. 方程组的建立

4. 方程组的求解

COMSOL中常见的数学概念

1. 梯度

2. 散度

3. 旋度

4. 通量

5. 绝对误差

6. 相对误差

7. 数值稳定性

8. 病态问题

COMSOL 的数值模拟是基于求解联立的偏微分方程(PDE)，而大多数PDEs是从质量守恒、动量守恒和能量守恒这些守恒律推导而来，这些PDEs采用积分形式的方程来描述物理变量在任意求解域内的行为。

对于界面的处理，可以通过高斯散度定理将表面积分转化为体积分，然而这就对方程自身的连续性、可导性提出了近乎苛刻的要求，一旦PDEs中涉及的变参不可导，则所有物理场的数值求解过程会十分不稳定。

为了应对这种挑战，COMSOL采用了将PDEs统一转化为弱形式（本质是将研究的物理变量变换为 积分形式，从而将不可导的风险转移到试函数项中)再进行求解的方法。

这种方法非常适合求解非线性较高的多物理场耦合问题。

另一方面，COMSOL 提供了一系列丰富的预定义物理场接口，用于模拟各种物理现象，其中囊括了多个交叉学科的物理耦合效应，譬如电热耦合、流固耦合和压电耦合等。

物理场接口是专门针对特定学科或工业领域问题建模的用户界面，用户可以在其中通过任意自定义的方程、函数、变量来建立模型。

COMSOL 作为一个集成的仿真平台，其工作流程包涵了仿真中涉及的所有步骤：

同时，软件内置有流体力学、电磁学、结构力学、化工等4大领域共计29个物理模块，用户可以通过上百个物理场接口任意组合创建研究所需的多场耦合模型。

有限元方法最早形成于工程力学，后来逐渐扩展到其他领域，并不断证明其优越性和可靠性。

该方法的理论基础是最小势能原理和最小余能原理。其基本思想是用许多规则形状的连续子域来近似代表整个求解域，这些子域称为“网格单元”。

在这些单元顶点处，物理量都精确满足原控制方程，而在单元内部任一点处的物理量，则是依据单元点处的值使用插值法求得。

这样，原控制方程(或偏微分方程)的连续性、可导性要求都被弱化，同时其精确性要求，由原来的 需要在整个求解域内处处满足控制方程(或偏微分方程)，被弱化为只要求在各个单元顶点处满足控制方程(或偏微分方程)， 而在单元内部，则通过插值方法求得其近似解。

这就是有限元“弱形式”的基本思想。

实际上，在COMSOL 中进行有限元计算时，系统总是先自动将控制方程转化为其弱形式，然后构造有限元方程组进行计算。

通过网格剖分，使得无限个自由度的原边值问题被转化成了有限个自由度的问题（即单元顶点处的方程)。

然后，用里兹变分法或伽辽金法得到一组代数方程（即方程组）。

最后，通过求解方程组得到这些单元顶点处方程的近似解，然后通过插值法求得单元内部任一点处的解。

所以，边值问题的有限元分析应包括下列基本步骤：

这一步即网格剖分，可以在COMSOL 中手动或者自动实现。

COMSOL 拥有强大的网格剖分功能，可以对复杂的求解域进行结构化和非结构化的网格剖分。将复杂的整体求解域离散为形状规则的网格单元之后，就可以在每一个单元上分别对原来的问题进行求解。

网格剖分是为构造插值函数和方程组做好准备。

在进行网格剖分之后，对于每个单元， 其内部任意一点处的物理量，是根据其顶点处求得的物理量的值，通过插值方法得到的。

在有限元方法中，根据所选用的插值函数的最高次数，确定插值函数的阶次，例如线性插值函数对应的是线性单元，二次插值函数对应的是二次单元，依此类推。

在对求解域进行网格剖分之后，各单元顶点处的物理量必须精确满足控制方程，以此得到一个方程组。

由于该方程组的每一个方程都是偏微分方程或者常微分方程，求解仍非常困难。因此需要通过里兹变分法或伽辽金法，将其近似变换为一个对应的代数方程组。

这一步从控制方程组到代数方程组的变换，是在系统内部自动完成的。

通过使用数值方法对该代数方程组进行求解，即得到控制方程在单元顶点处的近似解。

当近似解的误差小于指定的最大容许误差（即相对容差)时，即认为该近似解是准确的，可以作为该方程的解。

当得到代数方程组之后，点击求解器中的“计算”按钮，求解器就会以指定的相对容差为目标，使用数值方法求解方程组。由于方程数量较大、未知数较多，在进行数值求解时，首先将其转化为矩阵方程，然后通过直接消元法或者迭代法等多种数值方法对其进行求解，以得到代数方程组的近似解。

求解矩阵最常用的方法是迭代法，在每一步迭代的同时，都会计算当前结果的相对误差。当相对误差小于指定的相对容差时，迭代计算过程就会停 止，此时得到的解就认为已经具有足够的精度。

当得到各单元顶点处的近似解之后，系统会依据指定阶次的插值函数，计算单元内部各点上物理量的值，从而得到整个求解域上的近似解。

梯度(gradient)是一个矢量，它常用来表示某个物理量的变化快慢程度。

在同一空间或者时间范围内，某个物理量的变化率越大，则说明其梯度越大。

在COMSOL 中，可以在结果中用等值线观察物理量的梯度分布，在等值线分布越密的地方，则代表其梯度越大。

在进行有限元计算时，如果某个部位的梯度远大于其他部位，则需要将该梯度较大的部位使用较小的单元尺寸进行网格剖分，以获得较精度的结果。

比如对非常复杂的几何结构进行结构力学计算时，一般需要先根据理论预测，判断结构中应力梯度较大的部位，如应力集中区域等，从而对其进行局部网格细化。

如果理论预测有困难，可以先进行初步的均匀网格剖分和试算，根据结果找出应力梯度较大的区域，然后重新剖分网格并对这些区域进行局部细化，从而在确保计算精度的基础上，节省计算资源。

散度(divergence)用于表征空间各点矢量场发散的强弱程度，物理上，散度的意义是场的有源性。

当散度大于零时，表示该点有散发通量的正源（发散源）；当散度小于 零时，表示该点有吸收通量的负源（洞或汇）；当散度等于零时，表示该点无源。

由散度的定义可知，散度表示在某点处的单位体积内散发出来的矢量的通量，所以散度描述了通量源的密度。

举例来说，假设将太空中各个点的热辐射强度向量看作一个向量场，那么某个热辐射源（比如太阳）周边的热辐射强度向量都指向外，说明太阳是不断产生新的热辐射的源头，其散度大于零。

旋度是矢量分析中的一个矢量算子，可以表示三维矢量场对某一点附近的微元造成的旋转程度。

这个矢量提供了矢量场在这一点的旋转性质。旋度矢量的方向表示矢量场在这一点附近旋转度最大的环量的旋转轴，它和矢量旋转的方向满足右手定则。

旋度矢量的大小则是绕着这个旋转轴旋转的环量与旋转路径围成的面元的面积之比。

在COMSOL 中，关于散度和旋度的直接应用较少，而是直接设置相关的边界条件。

通量是指单位时间内，流经某单位面积的某属性量，是表示某属性量输送强度的物理量，如动量通量、热通量、物质通量和流量通量等。

在传热和流体等物理场中，经常需要设置热通量或者流量通量的边界条件。

测量值与真实值之差的绝对值，绝对误差只反映计量值与实际值的差别。某个物理量的绝对误差，与该物理量具有相同的量纲。

绝对误差与真实值的比值， 相对误差以百分比表示，数值越小表示计量 的精度越高。与绝对误差不同，相对误差的量纲永远为1。

在COMSOL 中，在几何与求解阶段都要考虑误差。

在几何建模阶段，根据实际情况，可以指定结构修复的相对误差或绝对误差。

而在求解阶段，一般需要指定的是相对误差，即相对容差。

稳定性是指算法对于计算过程中的误差（舍入误差、截断误差等）不敏感。

数值稳定时，在计算过程中随着计算的进行，相对误差会逐渐减小，直到小于设定的相对容差，即得到原问题的精确解。

数值不稳定时，在计算过程中相对误差反而会逐渐增大，或者出现上下波动，从而很难达到或者达不到设定的相对容差，最终很难得到或者无法得到问题的精确解。

在COMSOL 中进行计算时，可以通过收敛图对数值稳定性进行监测。

是指当输入数据（如参数、初始值等）有微小的波动时，会引起解的大的扰动。

由于计算工具总会存在舍入误差，因而对于病态问题，用任何算法求数值解都是不稳定的。可见，病态问题是数学模型自身的问题，与算法没有关系，病态问题的病态越严重，对数值计算稳定性的影响就越大。

病态问题会将问题的误差放大，当在计算过程中出现病态问题时，一般意味着边界条件或者求解器的设置出现了错误，需要进行检查。

在COMSOL 中使用有限元方法对问题进行求解时，其本质是使用数值方法求解矩阵方程。在数值方法中，求解矩阵的方法有直接法和迭代法。

直接法适用于小规模矩阵，它通过对矩阵求逆的方法来求解矩阵；迭代法以高斯消元法为基础，它适用于求解大规模矩阵。