伴随态格子玻尔兹曼方法：当成本函数在边界上包含宏观量时的必要修正

《Computers & Fluids》：Adjoint-state Lattice Boltzmann Method: A necessary correction when cost function incorporates macroscopic quantities on boundaries

　　摘要：流体流动的拓扑优化需要计算成本函数梯度。基于描述物理现象的偏微分方程，高效的计算方法为伴随状态法。本文指出，在LBM框架中，当成本函数包含边界处的宏观数据时，成本梯度往往难以准确计算。这是因为玻尔兹曼方程以介观尺度上的分布函数为基础，而成本函数则涉及宏观尺度上的变量。为此，本文提出了两种提高梯度计算精度的方法：调整成本函数的积分域以及采用聚合策略，其中聚合策略被证明能显著提升计算精度。

引言：在过去二十年中，基于格子玻尔兹曼方法的伴随状态法逐渐受到更多关注。这两种方法的结合通常被称为伴随状态格子玻尔兹曼方法。以下是关于在LBM基础上发展伴随状态方法的简略历史概述。最早涉及ALBM的文章可能是Tekitek等人[1]为解决多松弛时间背景下的参数估计问题所写的文章。在那篇具有开创性的论文中，他们考虑了完整的离散格子玻尔兹曼方案，即包含平流和碰撞两个步骤，并针对这两个步骤推导了伴随状态。随后，Pingen等人[2]将这种（完全离散的）方法应用于流体流动问题的拓扑与形状优化，重点研究了反弹边界条件的处理方法。他们采用了Spaid和Phelan[3]的孔隙率模型，通过全局应用的浸入边界条件来表示格子玻尔兹曼边界条件。之后，Krause等人[4]在三维空间中实现了ALBM方法，用于最优控制问题。他们强调了所谓的“先优化再离散化策略”，该策略使得伴随状态方程的推导可以在连续状态下进行，从而得到结构与LB算法相似的算法，这一特点对于实现高效并行计算尤为重要。Vergnault和Sagaut[5]将二维ALBM方法应用于噪声控制问题，采用了完全离散的格子玻尔兹曼方案，尤其是通过复微分与格子玻尔兹曼方法中的伴随问题相结合来实现牛顿下降法中的线搜索步骤。同时，Hekmat和Mirzaei[6]比较了基于LBM的连续型和离散型伴随状态方法在空气动力学中的应用，从数学推导和算法结果两方面进行了分析，他们指出伴随状态方程具备原始LB方程的所有优点（如简洁性、可并行性等），因此预计其计算成本与LB方程相近。Liu等人[7]将完全离散的ALBM方法应用于二维流体流动问题的拓扑优化，其中引入了体力来抵消流体流动效应，研究了扩散器、单弯管和双弯管等多种测试案例。Yaji等人[8]将水平集方法与LBM结合，用于隐式表示结构配置，他们在半离散方案上进行了伴随状态推导，即为了便于推导，玻尔兹曼方程在时间和空间上保持连续形式，研究了二维和三维空间中的单弯管和多弯管问题，还将该方法应用于换热器和氧化还原液流电池的拓扑优化[9]、[10]。?aniewski-Wo??k和Rokicki[11]在多GPU架构上实现了ALBM方法，将其应用于三维换热器和混合器的拓扑优化。Dugast等人[12]、[13]开发了一种用于包含热现象和反应-扩散现象的二维流体流动问题的伴随状态玻尔兹曼求解器，他们引入了指标置换的概念，结合分部积分法和格林定理，简化了伴随状态的推导过程，同时还证明了成本函数存在某些局部最小值，以及控制变量初始化对最优拓扑结果的影响。除上述研究外，还有许多其他工作涉及化学趋化系统的最优控制[14]、空气动力学应用中的形状优化[15]、域识别[16]以及拓扑优化领域的其他研究，详见[17]、[18]、[19]等文献。

本文旨在回归基础问题，特别是探讨在结合LBM推导伴随状态时如何处理边界条件。实际上，本文将证明，在LBM框架下，当成本函数在域边界上计算且包含宏观变量时，如果不进行特殊修正，其梯度就无法正确计算，因为计算得到的成本函数值可能与实际值存在巨大差异。这是因为玻尔兹曼方程以介观尺度上的分布函数为基础，而成本函数则涉及宏观尺度上的变量。实际上，玻尔兹曼问题仅需定义向内的边界条件，但宏观量会同时影响向内和向外的分布函数流动，这就是导致成本函数梯度计算出现误差的原因。

在优化领域，关于在计算泛函梯度时应采用离散化连续方法还是离散方法的争论由来已久。Chenais在1994年针对壳体形状优化问题进行了此类比较[20]，Favennec则针对非线性抛物型问题进行了研究[21]。根据这些研究，最佳策略是在完成所有近似处理后对正向模型求导，这样成本函数与其梯度之间就能保持完全一致性，因为成本函数是在正向模型求解完成后、即所有近似处理结束后才计算出来的。Ito等人最近在LBM背景下也得出了类似结论[22]。不过，对连续模型求导的好处在于可以选用与正向求解器不同的积分方案。本文采用的就是这种离散化连续梯度方法，它有助于制定必要的修正策略，以降低成本函数梯度计算中的误差。

本文结构如下：第2节以通用方式介绍了包含多种边界条件的正向模型的构建方法，适用于时变问题和稳态问题。为证明当成本函数包含边界处的宏观量时，成本函数与其梯度之间存在差异，有必要以非常严谨精确的方式构建完整的格子玻尔兹曼求解器，采用半离散方案，即至少在空间上保持连续，这样才能定义伴随状态边界条件。尽管人们普遍认为BGK碰撞算子在许多流动情况下不够稳定，但本文出于可读性考虑仍使用了该方案，本文的结论与所采用的碰撞方案无关。第3节以通用方式介绍了需要最小化的成本函数，并概述了基于Lions开创性数学工作的伴随状态问题的构建流程[23]。第4节是对伴随状态格子玻尔兹曼问题的数学推导，这一推导基于所谓的切线问题，并遵循Lions的方法，该节暂时不讨论边界条件问题，相关内容将在下一节阐述。文中通过著名的时变Taylor–Green涡流问题作为验证案例，目标是根据最终时刻的解反推初始条件。第5节重点讨论伴随状态问题中的边界条件处理，详细推导了反弹边界条件、周期边界条件以及Zou和He边界条件，同时解释了成本函数梯度计算中出现偏差的原因，本文提出了两种减少这种偏差的方法：第一种方法是修改成本函数的积分域；第二种方法是将宏观量中非向内的分量合并到向内的分量中，后者能显著降低偏差。需要注意的是，所有用于学术验证的测试案例都被刻意设计得较为简单，以便任何人都能轻松复现。第6节总结了本文的研究成果并提出了未来展望。

片段摘录：

正向模型：半离散玻尔兹曼方程组为
Ri=?fi?t+ei??fi+1τfi?fieq=0?(i,x,t)∈Q×D×I
其中，Ri表示与分布函数fi相关的残差，fieq代表第i个离散速度方向对应的平衡分布函数，索引集定义为??Q=?0,q?1?，D和I分别表示空间域和时间域，范围分别为[0,tf]。本文中采用的平衡分布函数表达式为
fieq=ωiρ1+3ei

优化问题的建立：我们定义如下成本函数：J以状态f显式表示，另外还有一个所谓的简化成本函数j，它是控制变量ψ的显式函数。在minj(ψ)这类优化问题中，ψ即为控制变量，通常直接使用这个函数。但实际上，成本函数需要在模型方程（2）求解后，作为状态f的函数来计算。就本文而言，...

伴随状态格子玻尔兹曼模型：在本节中，出于教学目的，我们在不考虑任何边界条件的情况下推导了伴随状态格子玻尔兹曼模型（ALBM），边界条件的处理将在下一节中进行。此外，本节仅讨论时变问题（对于稳态问题，既不考虑时间微分算子，也不考虑初始条件）。伴随状态的推导基于切线格子模型；后续还将进行验证...

边界条件处理：本节探讨伴随状态边界条件的推导方法，与前一节一样，首先假设存在线性切线边界条件，本节末尾的测试案例用于验证推导结果，为简化起见，此处仅讨论稳态问题。

讨论与总结：本文强调了当成本函数包含边界处的宏观数据时，必须修正伴随状态格子玻尔兹曼问题的边界条件。之所以需要这样的修正，是因为一方面LBM处理的是受平流型偏微分方程控制的介观变量，因此边界条件仅定义在向内方向上；而另一方面，边界处却需要考虑宏观量...

作者贡献声明：
Yann Favennec：写作——审阅与编辑、写作——初稿撰写、可视化、验证、监督、软件、方法论、研究、形式分析、概念设计。
Mahmoud Assaf：写作——审阅与编辑、验证、软件、方法论。
Christophe Josset：写作——审阅与编辑、验证、软件、方法论。

利益冲突声明：作者声明没有已知的财务利益或可能影响本文研究结果的个人关系。

Yann Favennec | Mahmoud Assaf | Christophe Josset
LTEN UMR CNRS 6607 - 南特大学，La Chantrerie，南特，44300，法国