读者对象仍然为初步接触求解器的研发人员，本文主要侧重于多物理场耦合，求解器算法和相关工具介绍1. GMRES：Generalized Minimal Residual2. BiCGSTAB：BiConjugate Gradient StabilizedJacobian-Free Newton-Krylov MethodJFNK核心思想是将牛顿迭代法与Krylov子空间方法相结合。传统的Newton-Krylov（NK）方法是通过利用牛顿迭代法求解非线性系统，并使用Krylov子空间法来逼近残差向量和牛顿方向的点乘，求出每次牛顿迭代的改进向量。然而，在传统的NK方法中，需要计算Jacobi矩阵，导致计算困难和时间复杂度较高。JFNK算法则是在NK方法的基础上进行了优化。它使用不同的策略来逐步逼近Jacobi矩阵，从而避免了昂贵的数值计算成本。特别地，JFNK算法直接构建雅各比矩阵的Krylov 近似矩阵，这样就能从迭代过程中得到更为精确的切线模型。因此，JFNK算法较之前的Newton-Krylov 更加高效，能够用较少的迭代步骤获得更高的数值精度与计算性能5. AMG：Algebraic Multigrid
AMG算法是一种基于代数方法的多重网格算法，它通过逐步将线性方程组转换为更简单的问题，从而快速求解大型稀疏的线性方程组。AMG算法是通过构造一个层次结构来处理线性方程组的解。这个层次结构包含多个较粗和较细的网格，并将原始的线性方程组分解成许多小的子问题，在每个子问题中使用预处理器或直接方法进行求解。通常情况下，该算法会优先考虑在当前较粗网格上直接求解小规模问题，然后再插值回到较细的网格中进行迭代求解。参见一篇文章入门多重网格方法
以上常用的几种方法除了JFNK，其它在以下书里有介绍，之前有推荐过：
METIS是一种常用的图分割程序，主要用于基于分治的图划分和元素重排等问题求解。在矩阵计算中，METIS被广泛应用于稀疏矩阵代数方程组的求解，主要步骤如下：1. 以稀疏矩阵的行或列为节点，构建相应的图结构； 2. 使用METIS进行图分割，将整个图分割成若干个不相交的子图； 3. 利用某种并行求解算法对每个子图进行求解，得到每个子图的解向量； 4. 通过合并各个子图的解向量，获取原始方程组的解向量。 METIS采用的基本策略是贪心策略，即从某个起点开始依次选择邻居节点来构建连通块，并对连通块进行规模划分，METIS通常采用四种划分策略：k-way、recursive k-way、refinement、和 time-limit 限制策略；METIS适用于稀疏、规则性数据较强的场合，而对于高度粗糙数据其效果可能不佳。此外，还要根据实际问题进行参数选择和优化，比如选定合适的划分方法、块数，需要一定的经验。 线性方程组求解的敏感度分析主要研究输入数据、参数误差对方程组解的影响。一般而言，可通过以下步骤实现： 1. 对输入数据/参数进行微小改变，如增加或减小一个较小的扰动 2. 计算新方程组的解，并与原来的解进行比较 3. 评估解的相对变化与参数扰动之间的关系
强耦合是指直接求解多个物理场的偏微分方程组，软件研发中的多物理场耦合仍然以弱耦合为主，也就是各种物理场单独计算，然后把单物理场的计算结果作为输入，传给另外一个物理场。传导又分为单向传导和双向传导。以电路为例，通电后导线温度会升高，变热的导线电阻会增加，进一步导致电流变化，从而影响温度，最终达到一个平衡，实际物理场是同步发生的。不同物理场在计算的时候会有一定的时间间隔，这是一个典型的双向传导。考虑到多物理场的复杂性，一般建议尽量使用单向传导。TSI：Thermal-Solid Interaction
热固耦合是多物理场耦合中最常见也最简单的耦合。通常是首先计算热场，将计算出的温度作为荷载加载到结构上，再计算结构的应力应变。且结构和热的网格兼容性好，使用相同的二次网格单元可以解决常见问题。温度网格单元的节点为标量，自由度为1，形函数物理含义明确，单刚和总刚计算都比较简单，不像电热耦合双向传递，从热场到结构是单向传递，也符合由易到难的顺序。参见一篇文章入门热学有限元
FSI：Fluid-Structure Interaction1.固体和流体部分或者完全重合，叠加在一起，流固两种介质难以分开。描述介质物理现象的本构方程受耦合效应影响，需要针对具有物理现象建立专门的偏微分方程(Partial Differential Equation)，比如多孔介质中的渗流，水凝胶等软物质的流固耦合。2.固体和流体在界面接触。流体载荷影响固体运动，固体运动又反过来影响流场的运动，这种相互影响是通过耦合界面的能量和信息交换实现。在流固耦合界面上，如果流体和固体的运动都未知，则要对整个耦合系统控制方程求解。气动弹性力学分析(气弹分析)是比较典型也是早期研究最多的流固耦合问题。之前求解气弹问题使用的Nastran，算是标准程序。在一篇文章入门多物理场有限元（全篇）一文中介绍了塔科马海峡吊桥垮塌事件，这是流固耦合计算发展史上一个重要节点。多物理场研发的内容，包括单个物理场模块封装和不同物理场之间的数据传递。有时候不同物理场的模型之间不是同步生成，所以在重合和边界的地方需要单独处理。MpCCI（Mesh-based parallel Code Coupling Interface）是由德国圣奥古斯丁SCAI研究中心（全称Fraunhofer Institute for Algorithms and Scientific Computing）开发出来的多物理场耦合商用工具。MpCCI接口软件可以实现不同模拟软件耦合区域的网格量的数据交换。由于耦合区域网格通常属于不同模拟程序，一般而言这些网格是不匹配的，MpCCI在实现网格值的数据交换前，先执行节点值之间的插值，求解器研发可以作为参考 
下图是之前介绍过的通用求解器框架，在这个结构中，将常用的业务流程封装为Solver Flow，在Solver Flow的基础上封装成一次单物理场迭代计算模块(One Interation Module)， 方便上层多物理耦合模块调用，这种设计目的也是为方便其它功能可扩展。封装是求解器开发一个重要指导原则
我们的物理世界是三维的，正常情况下，建模求解都采用三维模型很容易理解，但很多时候采用简化模型不仅求解效率更高，计算结果反而更准确。1. 3维简化到1维。整体结构分析中(比如建筑)往往很少使用三维实体单元(四面体/六面体/多面体)，只有局部分析时才会使用。结构分析中最常见的梁单元，如果导入或者建模时生成的是三维实体(比如长方体)，需要从几何上将其变成一条直线，然后计算截面惯性矩，形心矩等几何特征，将其作为属性赋值到梁单元上。所以有些专业建筑软件几何建模时，直接创建直线几何作为仿真单元，避免三维实体建模，反而简化了流程。结构分析中的杆单元，质点单元类似2. 3维简化到2维。结构分析中，针对薄板，一般采用壳单元计算，就是把一个三维实体在厚度方向简化成一个面，然后把厚度值作为属性赋到壳单元上。前处理软件中往往有抽中面功能，主要目的就是为了壳单元计算。电路分析中，针对元器件在横截面无变化的结构，可以将三维转成二维结构计算，流体也有类似的简化操作。二维分析相比三维分析网格数更少，更容易计算，也更容易验证求解器原型
3. EDA中的2.5维分析。2.5D仿真上世纪80年代由James C.Rautio博士提出，适合EDA电路中的层状结构分析，即使用三维全波公式，使用边界元/矩量法，考虑Z方向的结构厚度，不考虑Z方向的电流磁场变化。在某些场合下，2.5D计算结果优于3D，缺点是不能处理非层状，比如Bondwire结构，且对于边缘效应，介质精确建模等效率不高4. 对称结构。当整体结构出现两面，四面对称，可以采用对称单元，如两面对称可以只计算1/2模型，大幅减少计算量，需要额外处理的是对称面的边界和工况5. 模型清理。在前处理阶段会进行几何修复和清理，但也有在计算阶段发现问题，需要对计算模型修复清理。这个要求加强网格处理和求解器之间的关联处理6. 计算光学主要分为几何光学和波动光学：几何光学应用于物体尺度大于波长的光学现象，它是通过追踪光线来描述光在光学系统中传播的方法，适用于解决如透镜成像、折射、反射等问题；波动光学应用于物体尺度小于波长的光学现象，它是通过求解麦克斯韦方程组来描述光的传播，适用于解决如衍射、干涉、光栅作用等问题7. 结构中的模态分析主要用于确定模型的振动特性，即固有频率，振型，阻尼等参数，是结构受动态荷载分析中的重要参数，也是瞬态动力学，谐响应分析，谱分析等仿真内容的基础8. 在进行复杂的CFD分析前，通常会进行一些相对简单的分析，比如不考虑湍流，采用低雷诺系数，简化边界设置，其目的是对计算模型输入输出参数有一个大概了解
9. 高频电磁仿真中要获得频段范围内的S参数曲线，需要进行扫频，即对每个频点的S参数进行计算，然后拟合，如果单频点的计算误差较大，会导致最后的S参数曲线误差更大
可以看出，求解器模型的计算内容和简化方法，取决于模型特点和具体仿真内容。在求解器开发中，有一个原则：就是需要保证简单，基础的计算功能要精准可靠。这些计算功能要能进行封装，作为模块稳定的提供给上层或第三方使用。相比于CFD，EDA等领域的求解器，CAE尤其是结构领域，单元选择比较突出，通用商用软件的单元类型普遍都有几百种，而且还会扩展。对于通用结构求解器开发，要花一定的时间精力进行单元类型设计，这个关系到后续的功能扩展维护，是需要从软件架构层面考虑的内容。
求解器里涉及到多种迭代，包括解线性方程组迭代法，瞬态和时间相关的迭代，以及物理场自身算法迭代，还有网格自适应迭代。线性方程组迭代收敛涉及到预条件方法，初始值，迭代参数，网格质量等因素，也是求解器开发工作中比较麻烦的问题。瞬态仿真中有很多和迭代时间步长相关的设置。步长太小，浪费计算资源，过大，结果不正确。一般是求解器结合物理模型，网格尺寸，前一步计算结果，收敛情况给算出一个步长，比如CFD中的CFL数，热分析中的Biot和Fourier数，实际中很多时候是半自动设置，即用户可以给出一些经验参数。如果还在纠结开发语言，那就选C++，不是说C++本身有多好，而是C/C++在工业软件行业生态比较好。对于求解器开发，不仅要会常规的C++应用，还需要深入研究C++高性能计算，涉及到CUDA，OpenMPI，OpenMI，SIMD等领域，以及内存管理，CPU寄存器使用，Cache，寻址，磁盘网络性能等和硬件相关知识，这也是一个需要长期积累的过程。
|
