自主工业CAE软件，关键在于“自主”二字。

真正掌握了核心算法技术，才算是“子弹上膛”，

从而在国际赛场上占有一席竞争之地。

若是仅仅打着“自主技术”的幌子，

却没有应对技术洪流和市场危机的自主能力，

那就是对企业，甚至社会的不负责。

本文作者郭志鹏先生，会从求解器的开发过程开始，详细讲解适创科技是如何基于LBM算法开发了超快速、大尺度的流体力学求解器，从而完成了压铸领域模拟仿真云计算平台的研发。

我们的自主化技术，就是最核心的国际竞争力。

以下正文：

对于工业CAE软件来讲，其最核心的技术不是人-机交互界面、不是图形展示、不是数据库，而是其最底层的物理、数学算法架构，也就是我们一般称作的求解器。

所谓求解器，指的是针对特定场景（比如液体流动、温度传播以及结构变形），用程序编码的方式实现的对物理规律、数学原理的客观还原。

求解器在很大程度上类似于我们生活中使用的计算器，当我们在其屏幕上输入若干数字和运算符后，计算器会自动计算并输出结果。所有和计算相关的内部运算操作都集成在计算器的内部，使用者无法直观地获得对运算过程的观察和监控——就像个“黑盒子”一样。可以这么认为，生活中计算器涉及到的运算实际上是广义求解器的一个基础和组成部分。与计算器相比，CAE软件内部的求解器往往运算量更大，输入参数更多，对问题的还原更为全面。

求解器的开发过程

以流体力学求解器为例，其计算结果往往是一个三维区域内部的流体特性，包括速度、压力、形态等，随着时间的演化过程，每一个小的局部位置在特定时间点上都涉及到大量的代数运算，而这些代数运算其实就是我们生活计算器的主要工作。与实际人手操作计算器进行计算不同，CAE软件求解器可以根据物理规律演化的特点自动匹配当前涉及到的代数运算，这个就是我们俗称的“求解过程”。

从计算器-求解器的关联上看，开发求解器最主要的工作是把求解器的物理规律、数学算法通过编码的方式实现。这个过程看似清晰、简单，但却充满了挑战。

依然以流体力学求解器为例，我们来看一下这个开发的过程。

首先，我们需要深入理解描述流动过程的物理原理。对于液体来说，主导或者说决定其流动的底层物理规律是动量、质量守恒（不考虑热传导引起的能量耗散），而通过人类几个世纪的研究和总结，我们可以将整个流动过程用两个核心方程来描述，一个是Navier-Stokes方程（NS方程），描述动量守恒；一个是连续性方程，描述质量守恒。人们对于物理世界的理解往往是基于三维空间的，因此，将NS方程对笛卡尔坐标系进行拆分，将速度矢量沿着每一个坐标轴（X、Y、Z）进行分解，我们就会得到三个速度变量（U、V、W），这样的话一个NS方程就会变成针对离散速度变量的三个偏微分方程，再加上连续性方程，我们就可以形成一套完备的方程组来描述真实世界中不可压缩流体（水、金属液）的流动过程。

有了描述流体流动过程的方程组，下一步就是用数学的方法来求解这个方程组。数学方法一般分为两大类，一类是通过理论推导直接获得精确的方程组的解，另一类则是通过数值算法，“近似”地获得方程组的解。在现有的数学、物理框架下，或者说基于人类现有的科学发展状态，绝大多数的问题是无法获得精确解的，而解决实际问题的最主要途径是数值计算。

数值计算一般分为三个环节：模型离散、条件设定和求解过程本身。

“模型离散”指的是将真实的物理模型打碎成一个个的网格点，模型和网格点的关系特别像人体和细胞的关系；“条件设定”一般指的是正确地设定计算的输入参数，包括材料热物性参数、初始条件和边界条件等，正确的边界条件设定对于产生准确的计算结果是极其关键的，要不然就是“垃圾进-垃圾出（Trash in，trash out）”了；“求解过程本身”则是编码的数值算法在计算机硬件上的逐步、系统地运行。

数值计算的三个环节中，只有条件设定是完全依托于使用者的，模型离散和求解过程都可以且应该实现自动化。

我们把以上两个过程（描述物理规律的偏微分方程和数值算法）联系起来，通过计算机编码（高效算法编程一般使用C、C++以及fortran语言）实现成计算机可以识别的语言，也即计算机程序，最终通过特定计算机硬件平台和操作系统就可以实现对物理问题的求解了。

值得注意的是，描述物理问题的方程组经过离散之后的维度会急剧增加，这是因为这些偏微分方程组不仅需要建立在整体物理模型上，而且必须要在每一个网格单元上都得满足。因此，如果说一个物理模型经过网格离散获得了1000万个网格单元，那么相应的描述物理现象的方程组就得有1000万套，也就是说在每一个演化步中（一般是指时间步），我们要同时求解1000万套方程组。仍然以流动过程为例，对于每一个单元，在每一个时间步中，我们需要求解4个核心方程，包括描述U、V和W的三个速度变量的控制方程（NS方程）和描述质量守恒的连续性方程。从这个角度上讲，精确求解NS方程和连续性方程就是流动求解器需要解决的最关键问题。

现有的大多数工程应用涉及到的气体、水以及金属液等的流动都属于粘性不可压缩流动，针对这种流动，现有的流动求解器通常使用两种方法来实现对控制方程组的求解：基于压力迭代的算法和基于人造压缩特性的算法。这两种算法同时也是现有绝大部分商业软件采用的方法。

基于压力迭代的算法又可以细分为MAC算法和分裂算子算法，详细的算法原理可以参考相关的计算流体力学的书籍和文献，这里就不再赘述。MAC和分裂算子算法都可以通过数学变换最终转换成一个聚焦于求解泊松方程的算法，而求解泊松方程本身则需要采用所谓的隐式算法不断迭代才能完成。

事实上，求解泊松方程是近几十年来所有流体力学求解器最难啃的“骨头”。从数学上讲，泊松方程属于椭圆方程，每一个网格单元中心的压力项会直接和其周围的“邻居”单元中心的压力项形成关系，也即方程组，因此，对每一个单元来说，这个方程都会至少含有7个变量（当前单元自己的压力和周围单元的六个压力）。采用线性方程来进行描述就是：AX = b，其中，A是一个N行和N列的系数矩阵，X是含有N个压力项的变量，而b则是一个含有N个常量的向量。所以求解泊松方程就转换成了求解这个方程组。事实上，A是一个稀疏矩阵，因为其组元大部分都为零，在每一行中，只有7个（或多个，取决于采用的离散格式）位置的值不为零，可以想像，如果N=10,000,000，也即含有10^7个网格，那么这个矩阵就是一个极其稀疏的矩阵了。求解大型稀疏矩阵线性方程组是数值分析的主要任务，也是当代应用数学正在攻克的难点之一。

对于求解大型稀疏矩阵线性方程组，有兴趣的读者可以阅读一些和数值分析或者高等数值分析相关的一些书籍，经典的算法在其中都有提及，包括共轭梯度、SVD分解、幂指数迭代等等，我这里不做一一赘述。我想提两个算法，一个是公认的现有最强的矢量算法，叫GMRES，另一个则是多重网格（Multigrid）算法。

GMRES算法的核心在于建立正交坐标系，然后采用坐标系的坐标来表达最终的解。用一个小维度的例子来说明一下GMRES算法，比如我们之前AX=b的方程组，如果N=10的话，那么A就是一个10x10的系数矩阵，如果这个方程组有解的话，其解X必将是一个10x1的向量。从数学的角度理解AX=b的意思就是采用一个10维的系统将X映射到b，而这个过程其实是等效于建立正交坐标系。

GMRES在计算的过程中，绝大部分时间其实是在建立正交坐标系，因此，建立正交坐标系的时间也就最能反应GMRES算法的计算效率。在实际的工程问题上，我们发现，当计算维度不太大的时候，比方说百万、千万网格，GMRES算法的计算效率是非常高效的，但是当网格数量继续增加的时候，其效率就会急剧下降，这其实是受限于计算机浮点数本身的精度的，一个正交坐标系，如果每一个向量的值都有些许误差的话，到最后整个所谓的正交坐标系就不那么正交了，也就是相互之间不再垂直了。

而多重网格算法则没有这个问题，所谓的多重网格指的是这个算法在求解实际问题的时候采用的是不同尺寸的网格系，而非单层的网格结构。为什么要这么做呢，那是因为不同尺寸的网格对于消除“计算缺陷”的贡献是不同的，大尺寸网格消除较为粗糙的缺陷，而小尺寸网格消除的缺陷则更为“精细”，这有点像电视屏幕分辨率的感觉，当电视像素分辨率较低时，我们只能看到粗糙的图像，反之则可以看到更加清晰的画面。多重网格算法非常巧妙地利用了这点，在求解线性方程组的时候以更加高效、精准的方式来滤除计算误差，从而极大地加快了算法本身的求解速度。

虽然GMRES和多重网格算法是极其高效的，但是在求解实际工程问题中仍然面临极大的挑战。首先，不管是GMRES、多重网格还是其他各种算法，其本身对方程组的奇异性非常敏感，当系数矩阵A存在一定奇异性的时候，这些方法就会很快失效。什么是方程组的奇异性呢，其实就是矩阵A的特征值差异太大，我们大可以不用去理会什么是特征值，可以这么认为：矩阵的特征值反应了矩阵本身的“强度”，一个系数矩阵，当其某一个特征值在绝对值上远远小于其他特征值的时候，那么这个矩阵就是奇异的。

造成奇异的因素很多，可能是因为我们对实际问题的描述有误，也有可能是因为边界条件设定的不对，还有可能是这个问题本身就存在奇异性，比方说激波（shock wave）的存在就打乱了其两侧物理的连续性。那么，当这些奇异性存在时，采用这些所谓的“隐式”迭代算法就会存在收敛性的问题，甚至在某些时候就会计算中断，也就是我们俗称的“发散”了。现有商业软件开展流动过程模拟，出现计算不收敛的情况是非常常见的，当然，造成发散本身的原因也是多种多样，有操作者自身的问题，但绝大多数其实是求解器本身没有那么强大，或者说鲁棒性没那么好。

开发一个“完美”的流动求解器是极其困难的，而精确地求解Navier-Stokes方程也是当今世界7大数学难题之一（现在是6大，其中的黎曼猜想已经证实被解决了？）。这其实就是为什么我国之前并没有能够开发出一款具有世界竞争力的流体力学软件的本质原因。开发流动求解器，模拟真实自然界的各种跟流动相关的物理现象，实现真正的对工业乃至对宇宙的认识（星系爆炸、迁移和演变都涉及到流动）是现在绝大多数做应用数学、力学研究人员的终身梦想。

在真实的自然界中，跟流动相关的物理现象是非常复杂的。比如我们见到的绝大多数流动过程是涉及到多个流体的多相流，同时，当雷诺数超过2000左右的时候，流动会从层流变为紊流。多相紊流其实是自然界的常态，而真正从物理、算法的角度解决这类流动问题是极其困难的。现有商业软件在流动求解过程中都做过大量的假设和简化，所开发的求解器也只有在某些特定的场合才能产生相对有价值的计算和模拟结果，但是这些都不是使用者能够确切知道的，这其实也是模拟误差，或者说模拟精度不高的主要原因。

适创科技的自主求解器

面对这些困难，适创科技只能采用一种新的算法架构来开发流动求解器。这种新的算法就是LBM（算法本身并不新，只是相较于现有的流动求解算法来说，大规模有效使用LBM算法的时间并不长），也即格子-玻尔兹曼算法，同时采用大涡模型描述大雷诺数紊流。这两种算法的结合可以非常有效地规避现有流动求解器遇到的核心问题（迭代求解泊松方程这一环节）。

LBM方法其实结合了有网格和无网格方法的特点。

第一，LBM采用粒子碰撞的方法从微观来描绘基本的流动碰撞，但是这个方法不直接采用粒子碰撞来对流动状态做描绘，而是采用密度函数的方法将粒子碰撞进行积分，在相空间里对流态进行描述。采用密度函数的方法（而不直接采用粒子碰撞的方法）来连接微观碰撞和宏观流动是现代物理学的一个很重大的发现，这也是统计力学的基础。

从这个角度讲，LBM方法在基础微观层面采用粒子碰撞，而在宏观流动方面采用密度函数法，其计算精度（或者说对实际物理现象的描绘）要比现有的有、无网格方法更加准确。

“有网格方法”指的是数值计算中采用网格的概念开展的算法，常见的如FDM、FEM和FVM等；而“无网格方法”则是计算中没有采用网格这一概念，常见的如SPH、MLPG、RBF以及FPM。

但不论是哪种方法，其计算精度都依赖于“分辨率”，仍然像咱们电视的像素。对于有网格的计算方法，这个分辨率就是网格类型和大小，而对于无网格方法，这个分辨率则是粒子的数量。显然，我们绝对不可能用一个网格或一个粒子就把流体的流动状态描绘出来，因此它们都有“分辨率”的限制。

同时，为了提高计算效率，我们在数值算法上采用并行计算。LBM这种算法天生就具有极强的并行计算能力，非常适合于现有的并行计算集群，这样的话就很好地把算法本身的优势和现有的计算资源巧妙地结合了在一块。值得庆幸的是，用于并行计算的集群资源越来越多，国家在这方面的投入也越来越多，因此可以极大地把计算成本压低。

不同于大多数现有商业软件的算法架构，适创科技开发的流动求解器，是基于LBM算法、大涡模型和并行计算为主要特点的新一代流体力学求解器，本身代表了流动计算算法的最高水平，其计算效率和解决复杂问题的能力也远远超越了绝大多数现有的商业软件。

而这个求解器的诞生，也为后期适创科技拓展到其他计算应用领域奠定下了坚实的基础，这也是我们一直引以为豪并乐以提及的：

“自主化技术，世界竞争力。”