物理求解的伴随

求解建模物理现象伴随的根本任务是计算成本函数相对网格的灵敏度 $\frac{d L^{T}}{d X}$ 。

成本函数相对网格的灵敏度是 Eqn. (5067) 中最右侧的两项之积：

图 1. EQUATION_DISPLAY

\frac{d L^{T}}{d X} = \frac{\partial L^{T}}{\partial X} + \frac{\partial Q^{T}}{\partial X} \frac{\partial L^{T}}{\partial Q}

(5075)

其中，此灵敏度的计算分为两步：

迭代求解伴随流体求解，得出报告目标相对控制方程残差的灵敏度， $\frac{d L^{T}}{d R}$ 。
计算相对网格的残差 $\frac{\partial R^{T}}{\partial X}$ 和伴随求解之积，并将其添加到目标的显式网格依赖性 $\frac{\partial L^{T}}{\partial X}$ 。收敛伴随求解后，第二步作为一次性操作执行。

这些步骤需要对原始求解过程进行微分，包括流体方程、能量方程和湍流方程的求解（适用时）。下面介绍了耦合流体模型求导。耦合固体能量模型求导遵循类似的模式。

原始求解

Simcenter STAR-CCM+ 内的迭代求解器的目标是推进一组变量，使它们可以满足一些离散控制方程组。例如，耦合流体模型的控制方程是纳维-斯托克斯方程 Eqn. (5076)。

预处理纳维-斯托克斯方程表示如下：

图 2. EQUATION_DISPLAY

Γ \frac{\partial}{\partial t} \int_{V} Q d V + \oint [F - G] \cdot d A = \int_{V} H d V

(5076)

其中， $Γ$ 为预处理矩阵， $F$ 和 $G$ 分别为无粘性通量和粘性通量， $H$ 为体积力。

对于一般性，假设有 $N$ 个不同的迭代求解器，每个迭代求解器涉及推进一个变量子集，使它们满足某些控制方程。例如对于湍流问题，这些求解器将是耦合隐式求解器、湍流求解器和湍流粘度求解器。对于 $i^{t h}$ 求解器，变量集称为 $Q_{i}$ ，并且可以采用残差形式将离散控制方程表示为 $R_{i}$ 。由于求解器是耦合的，因此该残差可以是模拟中所有变量的函数。使用此方法，原始求解将归结为求解以下一系列耦合非线性方程：

图 3. EQUATION_DISPLAY

\begin{matrix} R_{1} (Q_{1}, Q_{2} \dots, Q_{i}, \dots, Q_{N}, X) = 0 \\ \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} R_{2} (Q_{1}, Q_{2} \dots, Q_{i}, \dots, Q_{N}, X) = 0 \\ ⋮ \end{matrix} \\ R_{i} (Q_{1}, Q_{2} \dots, Q_{i}, \dots, Q_{N}, X) = 0 \end{matrix} \\ ⋮ \end{matrix} \\ R_{N} (Q_{1}, Q_{2} \dots, Q_{i}, \dots, Q_{N}, X) = 0 \end{matrix} \end{matrix}

(5077)

对于稳态问题，此过程采用非线性块高斯-赛德尔方法，每个求解器都基于最新的变量值构建其残差，并基于局部预解算子 $[P_{i}]$ 推进其变量集。对于求解器 $i$ ，非线性迭代 $k$ 的更新给定为：

图 4. EQUATION_DISPLAY

\begin{matrix} R_{i}^{k} = R_{i} (Q_{1}^{k + 1}, Q_{2}^{k + 1} \dots, Q_{i}^{k}, \dots, Q_{N}^{k}) \\ Q_{i}^{k + 1} = Q_{i}^{k} - {[P_{i}]}^{- 1} R_{i}^{k} \end{matrix}

(5078)

对于耦合流体模型，预设条件是相对 (P, U, V, W, T) 的流体变量的纳维-斯托克斯方程的近似线性化，其定义如下：

图 5. EQUATION_DISPLAY

[P] = \frac{V}{δ t} Γ + \sum_{j}^{N_{f}} (\frac{\partial F_{j}}{\partial Q} - \frac{\partial G_{j}}{\partial Q}) A_{j}

(5079)

伴随求解

要派生稳态伴随求解，必须对原始求解中定义的耦合非线性方程进行微分和转置。要确定相对于 $X$ 的导数，链规则会应用于每个原始残差方程：

图 6. EQUATION_DISPLAY

\begin{matrix} \frac{d R_{1}}{d X} + \frac{\partial R_{1}}{\partial Q_{1}} \frac{\partial Q_{1}}{\partial X} + \frac{\partial R_{1}}{\partial Q_{2}} \frac{\partial Q_{2}}{\partial X} + \dots + \frac{\partial R_{1}}{\partial Q_{i}} \frac{\partial Q_{i}}{\partial X} + \dots + \frac{\partial R_{1}}{\partial Q_{N}} \frac{\partial Q_{N}}{\partial X} = 0 \\ \frac{d R_{2}}{d X} + \frac{\partial R_{2}}{\partial Q_{1}} \frac{\partial Q_{1}}{\partial X} + \frac{\partial R_{2}}{\partial Q_{2}} \frac{\partial Q_{2}}{\partial X} + \dots + \frac{\partial R_{2}}{\partial Q_{i}} \frac{\partial Q_{i}}{\partial X} + \dots + \frac{\partial R_{2}}{\partial Q_{N}} \frac{\partial Q_{N}}{\partial X} = 0 \\ \begin{matrix} \begin{matrix} \begin{matrix} ⋮ \\ \frac{d R_{i}}{d X} + \frac{\partial R_{i}}{\partial Q_{1}} \frac{\partial Q_{1}}{\partial X} + \frac{\partial R_{i}}{\partial Q_{2}} \frac{\partial Q_{2}}{\partial X} + \dots + \frac{\partial R_{i}}{\partial Q_{i}} \frac{\partial Q_{i}}{\partial X} + \dots + \frac{\partial R_{i}}{\partial Q_{N}} \frac{\partial Q_{N}}{\partial X} = 0 \end{matrix} \\ ⋮ \end{matrix} \\ \frac{d R_{N}}{d X} + \frac{\partial R_{N}}{\partial Q_{1}} \frac{\partial Q_{1}}{\partial X} + \frac{\partial R_{N}}{\partial Q_{2}} \frac{\partial Q_{2}}{\partial X} + \dots + \frac{\partial R_{N}}{\partial Q_{i}} \frac{\partial Q_{i}}{\partial X} + \dots + \frac{\partial R_{N}}{\partial Q_{N}} \frac{\partial Q_{N}}{\partial X} = 0 \end{matrix} \end{matrix}

(5080)

上述微分的结果是所有求解器之间存在的一个大的、完全耦合的线性系统：

图 7. EQUATION_DISPLAY

[\frac{\partial R}{\partial Q}] \frac{\partial Q}{\partial X} = - \frac{\partial R}{\partial X}

(5081)

与用于推进原始的预设条件矩阵不同，上述过程中的线性化是精确的，并考虑二阶项。

图 8. EQUATION_DISPLAY

\frac{\partial Q^{T}}{\partial X} = - \frac{\partial R^{T}}{\partial X} {[\frac{\partial R}{\partial Q}]}^{- T}

(5082)

由于 ${[\frac{\partial R}{\partial Q}]}^{- T}$ 的反算成本较高，因此将 Eqn. (5081) 与先前计算的 ${\frac{\partial L}{\partial Q}}^{T}$ 结果相乘（请参见求解分析）：

图 9. EQUATION_DISPLAY

{\frac{d L}{d X}}^{T} = {\frac{\partial L}{\partial X}}^{T} - {\frac{\partial R}{\partial X}}^{T} {[\frac{\partial R}{\partial Q}]}^{- T} {\frac{\partial L}{\partial Q}}^{T}

(5083)

其中， ${\frac{\partial R}{\partial X}}^{T}$ 已知，因为 $R$ 已知。最后两项定义伴随求解变量 $Λ_{Q}$ ：

图 10. EQUATION_DISPLAY

{Λ_{Q} = [\frac{\partial R}{\partial Q}]}^{- T} {\frac{\partial L}{\partial Q}}^{T}

(5084)

该变量可以解释为 ${\frac{d L}{d R}}^{T}$ 。

为了避免计算和组合大型复杂矩阵，伴随求解过程利用仅需要矩阵矢量积的算法，即重启的 GMRES 和缺陷校正。缺陷校正可以写为：

图 11. EQUATION_DISPLAY

{[P]}^{T} {δ Λ}_{Q} = \frac{\partial L}{\partial Q} - [\frac{\partial R}{\partial Q}] Λ_{Q} = - R_{Λ_{Q}} (Λ_{b f Q})

(5085)

对于 GMRES，系统直接求解，Eqn. (5085) 用作 GMRES 算法中的预设条件。要在多个求解器中实施缺陷校正扫掠，原始求解过程将以相反顺序应用，并且每个伴随求解器在整个矩阵的列上操作。每次迭代时，求解器将计算伴随残差，并使用局部预设条件的转置来推进其伴随变量。唯一偏离原始结构的情况是，伴随以类似块雅可比方式进行推进，所有求解器在计算其伴随残差时都使用上一次迭代的伴随求解。以这种方式推进求解，可以在求解多个伴随求解时节省内存。

图 12. EQUATION_DISPLAY

\begin{matrix} R_{Λ_{Q, i}}^{k} = \frac{\partial R_{1}}{\partial Q_{i}} Λ_{Q, 1}^{k} + \frac{\partial R_{2}}{\partial Q_{i}} Λ_{Q, 2}^{k} + \dots + \frac{\partial R_{i}}{\partial Q_{i}} Λ_{Q, i}^{k} + \dots + \frac{\partial R_{N}}{\partial Q_{i}} Λ_{Q, N}^{k} \\ Λ_{Q, i}^{k + 1} = Λ_{Q, i}^{k} - {[P_{i}]}^{- T} R_{Λ_{Q, i}}^{k} \end{matrix}

(5086)

上述过程表示原始过程的完全耦合伴随求解。例如，如果原始过程由使用 SA 湍流模型的流求解组成，则上述选项表示选择伴随 SA 模型时的实施。还可以为特定求解器引入冻结处理。这是在选择伴随冻结湍流模型时应用于湍流求解器的处理。当求解器冻结时，表示这些求解器残差的行和表示这些求解器变量的列将从 Eqn. (5081) 的线性系统中移除。

此移除的实际结果是，相应的伴随求解器可以排除在该过程之外，从而提高了伴随求解过程的稳健性和性能，代价是降低了精度。