湍流边界层广泛存在于航空、航天、建筑、生物等领域，其主要的计算方法包括直接数值模拟(DNS)、大涡模拟(LES)和雷诺平均Navier-Stokes(RANS)模拟。直接数值模拟求解非定常流动中的湍流扰动，大涡模拟湍流小尺度结构而分辨湍流大尺度结构，具有较好的预测精度，但DNS和LES所需的计算网格规模与雷诺数的指数幂成正比，在高雷诺数下的计算成本巨大，使得利用DNS或LES开展高雷诺数流动仿真存在困难^[1]。RANS通过引入湍流模型方程模化所有的时空多尺度湍流脉动，通过显式表达式^[2]或偏微分方程^[3]求得雷诺应力，其计算量与雷诺数呈弱相关关系，计算效率高，在工程中得到广泛应用。

然而，RANS的预测精度与DNS或LES存在差异，且在实际应用中湍流模型的选择依赖于设计者的经验, 通过调整或优化湍流模型中的个别系数难以提高其预测精度，且不具有普适性。RANS计算最大的不确定性来源于湍流模型方程的固有形式^[4-5]，近年来，利用数据挖掘和机器学习方法修正RANS湍流模型显示出巨大的潜力和前景^[6-7]。

Duraisamy等^[8-9]利用流场反演获得湍流模型修正项，进而使用神经网络建立流场特征与修正项之间的映射关系，获得嵌入式的湍流模型，经重构后湍流模型计算结果更接近实验数据。Ling和Templeton^[10]利用DNS和LES的计算结果，发展了机器学习修正的RANS湍流模型，比较了不同机器学习方法，如支持向量机、决策树以及随机森林对流场不确定性的预测效果，结果表明这些方法在与训练样本相异的流场中均能准确地捕捉流动特征。Xiao等^[11]使用贝叶斯方法量化RANS计算中模型形式的不确定性，并加入经验知识约束，结果表明该方法利用稀疏的先验数据能较好地捕获基准模型的误差，获得更准确的后验平均速度分布。Wang等^[12]的研究结果表明，雷诺应力的预测能力是影响RANS计算精度的重要因素，由此建立了基于雷诺应力张量的随机森林机器学习重构方法，并针对湍流方管及典型分离流动进行了模型重构。Wu等^[13]采用相似的方法构造了随机森林模型，给出了RANS计算置信度的先验估计。Zhang等^[14]在翼型的近壁区、尾迹区等不同区域分别构造了径向基神经网络，模化了高雷诺数湍流边界层中的涡黏项，并将其实时地引入流场求解过程，提高了计算准确度，且在翼型或流动条件变化时具有一定通用性。在湍流计算领域，基于大数据的机器学习方法具有广阔的应用前景，文献[15]对数据驱动式湍流建模理论和方法进行了综述分析。

数据驱动式建模方法依赖于样本数据的个数和质量，小样本或样本冗余误差会导致模型的失真或失效。为克服数据驱动式建模对数据样本个数和质量的依赖关系，需要引入基于物理知识约束的先验条件，在建模过程中约束模型的边界和搜索空间，有利于避免非物理现象的产生，提高模型的预测精度。近年来，基于物理知识的数据驱动式建模成为研究热点，Raissi等^[16]以自由振动的圆柱为研究对象，将有限的速度散点数据作为训练样本，利用神经网络重构的过程中，引入流体控制方程作为先验的物理知识约束，对流动参数进行了预测，获得准确的升/阻力预测结果。Raissi等^[17]构造了基于物理信息的神经网络，用于求解非线性偏微分方程，利用少量的训练数据，成功捕捉了系统的非线性特征。针对湍流模拟的RANS方程，尚未发现基于物理知识约束的数据驱动式建模方法。

准确预测湍流摩擦阻力具有重要的工程价值，本文作者在利用DNS数据研究湍流摩擦阻力的分解中，发现湍流摩擦阻力的分解项在湍流边界层的内层具有统一的比拟关系^[18]，与雷诺数大小无关，与可压缩性无关^[19]，并在槽道和平板湍流边界层中得到了验证^[18-20]，结果表明该比拟关系是一种物理的、普适的规律，是一种先验的物理知识。本文基于该物理知识开展了针对S-A湍流模型的数据驱动式建模及修正，在建模过程中引入物理知识约束，对比了有/无物理知识约束的修正对槽道湍流摩擦阻力预测精度的影响。

1 基于物理知识约束的数据驱动式建模方法 1.1 经典S-A湍流模型

为封闭雷诺平均的Navier-Stokes方程，Spalart和Allmaras^[21]构造了一方程湍流模型方程，通过求解标量运输方程来计算湍流涡黏系数并获得雷诺应力，具体形式为

$ \frac{\partial \hat{\nu}}{\partial t}+u_{j}\left(\frac{\partial \hat{\nu}}{\partial x_{j}}\right)=P-D+T $	（1）

式中：$ \hat{\nu}$为湍流涡黏系数求解变量；P为生成项；D为“破坏”项；T为耗散项，具体形式为

$ P=C_{\mathrm{b} 1}\left(1-f_{\mathrm{t} 2}\right) \hat{S}{\hat{\nu}} $	（2）

$ D=\left(C_{\mathrm{w} 1} f_{\mathrm{w}}-\frac{C_{\mathrm{b} 1}}{\kappa^{2}} f_{\mathrm{t} 2}\right)\left(\frac{\hat{\nu}}{d}\right)^{2} $	（3）

$ T=\frac{1}{\sigma}\left[\frac{\partial}{\partial x_{j}}\left((\nu+\hat{\nu}) \frac{\partial \hat{\nu}}{\partial x_{j}}\right)+C_{\mathrm{b} 2} \frac{\partial \hat{\nu}}{\partial x_{i}} \cdot \frac{\partial \hat{\nu}}{\partial x_{i}}\right] $	（4）

其中：参数C_b1、f_t2、C_w1、f_w、κ、d、σ、ν、C_b2和$\hat S $的定义与文献[21]一致，不再赘述。求解RANS方程的数值方法可参见文献[22]，对流项使用三阶MUSCL格式，界面通量计算方法为Roe格式，时间推进采用LU-SGS格式。

1.2 源于湍流摩擦阻力分解的物理知识

固壁与流体相对运动而产生摩擦，摩擦阻力可以表征为固壁对流体做功的一种形式，对于槽道湍流而言，摩擦阻力做功，一部分能量通过黏性直接耗散转换成热能，一部分能量则通过湍动能生成项维持湍流态。基于这一思想，Renard和Deck^[23]推导了湍流摩擦阻力系数C_f的分解公式：

$ C_{f}=\underbrace{\frac{2}{u_{\mathrm{b}}^{3}} \int_{0}^{h} \nu\left(\frac{\partial U}{\partial y}\right)^{2} \mathrm{d} y}_{C_{f{1}, \mathrm{RD}}}+\underbrace{\frac{2}{u_{\mathrm{b}}^{3}} \int_{0}^{h}\left(-\left\langle u^{\prime} v^{\prime}\right\rangle\right) \frac{\partial U}{\partial y} \mathrm{d} y}_{C_{f{2}, \mathrm{RD}}} $	（5）

式中：U为平均流向速度；u′和v′分别为流向和法向脉动速度；〈u′v′〉为雷诺应力；体速度$u_{\mathrm{b}}=\frac{1}{h} \int_{0}^{h} U \mathrm{d} y, h$为槽道半高，y为壁面距离；ν为运动黏度系数。式(5)适用于不可压缩槽道湍流，可压缩流动中的湍流摩擦阻力分解方法可参考文献[19-20]。

文献[18]分析了不同雷诺数槽道湍流的摩擦阻力分解，发现直接黏性耗散项$\mathit{ \Omega}=y^{+}\left(\frac{\partial U^{+}}{\partial y^{+}}\right)^{2}$ (其中上角标+表示以黏性尺度无量纲后的变量)主要集中在近壁区，其法向分布规律与雷诺数无关，仅与法向高度相关，在y⁺≈6.5处达到峰值。直接黏性耗散项的法向分布规律是一种先验的物理知识，该物理知识中没有耦合参数依赖关系，满足自相似性，可作为约束条件用于数据驱动式建模。

图 1给出了不同雷诺数Re_τ条件下直接黏性耗散项的法向分布曲线，曲线下的面积表征了$\sqrt{2} C_{f 1, \mathrm{RD}} / C_{f}^{1, 5} $预乘积分的数值结果，揭示直接黏性耗散分量在总摩擦阻力中贡献的占比情况。为保持通用性，对直接黏性耗散项的法向分布规律进行数值拟合，拟合公式为

$ \begin{array}{c} \mathit{\Omega}\left(y^{+}\right)=a_{1} \exp \left(-\left(\frac{\lg \left(y^{+}\right)-b_{1}}{c_{1}}\right)^{2}\right)+ \\ a_{2} \exp \left(-\left(\frac{\lg \left(y^{+}\right)-b_{2}}{c_{2}}\right)^{2}\right) \end{array} $	（6）

式中：a₁=2.644；b₁=1.895；c₁=-0.805；a₂=1.777；b₂=1.118；c₂=1.642。式(6)中仅有一个自变量y⁺，图 1表明拟合结果与DNS数据的符合良好。本文将式(6)作为一种先验物理知识约束引入到数据驱动建模中。

1.3 基于物理知识约束的湍流模型修正方法

保留S-A湍流模型方程的基本形式，对生成项P添加修正因子β，修正后的湍流模型方程为

$ \frac{\partial \hat{\nu}}{\partial t}+u_{j}\left(\frac{\partial \hat{\nu}}{\partial x_{j}}\right)=\beta P-D+T $	（7）

修正因子β对流场的影响是全局性的，其取值随着法向高度的变化而变化。本文通过添加修正因子β来校正S-A湍流模型，其本质是在S-A湍流模型中添加了源项的校正量δ，即δ=(β-1)P，其大小与生成项相关。初始时β=1，校正量δ均为0，式(7)即退化为原始的S-A方程。通过引入修正因子β，拟使得S-A湍流模型预测更准确，即校正后的S-A湍流模型计算结果更接近于真实解，真实解可取高精度的DNS、LES结果，也可利用试验测量的数据进行标定，本文中的真实解为DNS计算结果^[24]。

为获得准确可靠的修正因子，需要设定合理的目标函数，结合高效的优化方法，快速获得修正因子分布。传统方法对模型中系数进行调整优化，例如对冯卡门常数(κ)或湍流普朗特数(σ)的校正^[24-25]，但其难以克服重要流场物理特性(如雷诺数)的影响。由此，考虑到湍流摩擦阻力的预测与平均流向速度剖面紧密相关，且引入式(6)作为物理知识约束，设定目标函数为

$ \begin{array}{c} \min J_{\mathrm{ud}}=\sum\limits_{j=1}^{N_{\mathrm{c}}}\left(U_{j}-U_{j, \mathrm{DNS}}\right)^{2}+ \\ \sum\limits_{j=1}^{N_{\mathrm{c}}}\left(\mathit{\Omega}_{j}-\mathit{\Omega}\left(y^{+}\right)\right)^{2} \end{array} $	（8）

式中：(U_j－U_{j, DNS})²表征对平均流向速度剖面的逼近程度；(Ω_j－Ω(y⁺))²表征对近壁黏性耗散项分布的逼近程度；N_c为网格单元数。值得指出的是，修正因子β与直接黏性耗散和平均流向速度之间不存在直接的量化关系，仅存在隐式关联与映射。

为对比分析有/无物理知识约束对湍流摩擦阻力的预测误差，本文也采用了基于平均流向速度的目标函数：

$ \min \quad J_{\mathrm{u}}=\sum\limits_{j=1}^{N_{\mathrm{c}}}\left(U_{j}-U_{j, \mathrm{DNS}}\right)^{2} $	（9）

由于设计变量个数与网格数量相当，使用有限差分法求解梯度计算量巨大。伴随方法(Adjoint Method)是一种高效的梯度求解方法，其计算规模与设计变量的数目基本无关，可大幅减小梯度的求解时间。伴随方法以偏微分方程系统控制理论为基础，Jameson^[26]首次将伴随方法应用于气动设计中，把气动设计转换为一个满足特定约束的最优控制问题，梯度求解的计算量约为2倍的流场计算量。伴随方法包括连续伴随和离散伴随，本文采用离散伴随方法^[27]，伴随方程和梯度的表达式为

$ \left[\frac{\partial \boldsymbol{R}}{\partial \boldsymbol{Q}}\right]^{\mathrm{T}} \boldsymbol{\psi}=-\left[\frac{\partial \boldsymbol{J}}{\partial \boldsymbol{Q}}\right]^{\mathrm{T}} $	（10）

$ \frac{\mathrm{d} \boldsymbol{J}}{\mathrm{d} \beta}=\frac{\partial \boldsymbol{J}}{\partial \beta}+\boldsymbol{\psi}^{\mathrm{T}} \frac{\partial \boldsymbol{R}}{\partial \beta} $	（11）

式中：R为原方程的残差；Q为模型的求解量；J为目标函数。通过构造偏导数矩阵$\left[\frac{\partial \boldsymbol{R}}{\partial \boldsymbol{Q}}\right] $可获得伴随方程，由于该矩阵为大型稀疏矩阵，可使用Intel函数库Math Kernel Library中的迭代求解器对伴随方程其进行求解，采用广义最小残差(GMRES)算法^[28]。求得伴随因子ψ后将其代入式(11)可得到目标函数关于设计变量β的导数。优化算法采用基于梯度的最速下降法(Steepest Descent Algorithm)，进行迭代更新并最终求得设计变量β的分布。

湍流模型修正优化流程如图 2所示，首先基于原始S-A模型完成RANS方程的求解，由流场变量及方程残差构造偏导数矩阵$\left[\frac{\partial \boldsymbol{R}}{\partial \boldsymbol{Q}}\right] $及关于目标函数的右端项$\left[\frac{\partial \boldsymbol{J}}{\partial \boldsymbol{Q}}\right] $，然后使用GMRES算法对伴随方程进行求解，获得目标函数关于设计变量β的梯度，利用最速下降法进行β的更新，重新构造新的偏导数矩阵$\left[\frac{\partial \boldsymbol{R}}{\partial \boldsymbol{Q}}\right] $及右端项$\left[\frac{\partial \boldsymbol{J}}{\partial \boldsymbol{Q}}\right] $，迭代计算直至目标函数收敛。

以二维方程^[29]来验证伴随方法所求得的梯度与有限差分法的一致性，图 3对比了有限差分法与伴随方法所得的梯度分布，其偏差小于0.57%，属于合理范围内。

2 计算结果与分析

针对槽道湍流边界层，选取3组雷诺数(Re_τ=180、550、1 000)作对比分析。采用式(9)为目标函数的算例，用J_u进行表示，采用式(8)为目标函数的算例，用J_ud表示，将计算结果与DNS数据进行比较。对S-A湍流模型中的生成项P添加修正因子β，即对原湍流模型方程添加了源项的校正量δ=(β-1)P，进行了槽道流场的修正，经过基于梯度的迭代优化，壁面法向的黏性耗散分布如图 4所示，可以看出相较于原始算例，采用J_u与J_ud为目标函数的算例均获得了一定的修正效果，其中J_ud的修正效果更为显著。图 5给出了定量的误差比较结果，定义相对误差为$ E_\mathit{ \Omega}=\frac{\mathit{ \Omega}_{\mathrm{SA}}-\mathit{ \Omega}\left(y^{+}\right)}{\mathit{ \Omega}\left(y^{+}\right)}$，在黏性耗散项非零的区域内(y⁺ < 20)，Ω_SA为S-A计算结果中的直接黏性耗散，原始S-A模型的峰值误差在3个不同雷诺数下分别达到50%、28%和33%，经J_u优化后相对误差为21%、22%和26%，优于原始模型的结果，而经J_ud优化后这3处的相对误差减小至3.01%、0.16%和1.39%以内，显著优于基于J_u的优化效果。

图 4和图 5中的比较结果表明，原始S-A湍流模型对近壁区直接黏性耗散项的预测存在明显偏差；如果通过数据驱动模型仅修正流向速度剖面，对直接黏性耗散项的预测有所改善，但误差仍然较大；在数据驱动建模的目标函数中，进一步引入关于直接黏性耗散的先验物理知识作为约束条件，预测精度显著改善，计算结果与真实结果的相对误差小于3%。

图 6展示了平均流向速度U⁺的法向分布情况，为定量描述修正效果，定义RANS和DNS平均流向速度的相对误差为$ E_{\mathrm{u}}=\frac{U_{\mathrm{SA}}-U_{\mathrm{DNS}}}{U_{\mathrm{DNS}}}$，U_SA为S-A计算结果中的平均流向速度，U_DNS为DNS计算结果中的平均流向速度，其法向分布如图 7所示，结果表明在不同雷诺数条件下，采用目标函数J_u和J_ud修正后的S-A湍流模型的预测精度均优于原始S-A湍流模型的结果，更加趋近于DNS的真实解。当Re_τ=180时，原始S-A湍流模型计算结果的最大相对误差由-8.4%减小至-1.5%以内。随着雷诺数增大，原始S-A湍流模型计算结果的误差有所减小，但基于数据的修正仍能进一步减小误差，且加入先验物理知识约束的J_ud修正具有更好的修正效果。Re_τ=550，1 000时，平均流向速度的最大相对误差位置由壁面变化到y⁺=10左右，基本与图 5中直接黏性耗散项的最大误差位置一致，说明在高雷诺条件下的平均流向速度的预测精度受到湍流变量预测精度的影响，而引入湍流变量(本文中为直接黏性耗散)的多参数修正模型，可能会取得更好的预测效果。

图 8比较了不同雷诺数下修正因子β的变化情况，注意原始S-A湍流模型中修正因子恒为1。当采用J_u为目标函数时，仅根据流向平均速度进行模型修正，修正因子β的法向变化趋势基本与图 7中流向平均速度误差的分布一致，说明修正因子β具有较好的活性，可根据当地的平均速度误差对S-A模型中的湍流黏性生成项进行相应调整。当采用J_ud为目标函数时，引入了直接黏性耗散项作为流向平均速度修正的补充，以Re_τ=180的槽道为例，对于复合目标函数修正因子β做出相应的调整，在y⁺=4.65处达到峰值，说明目标函数J_ud起到了复合校正效果。同时，图 8中β系数随流场位置及当地流场变量的变化而发生演变，也说明了本文建立的伴随优化方法计算可以有效地建立对方程预测精度的调整和修正。需要注意的是，修正因子β的取值范围在-10~20之间，说明原始S-A湍流模型中的生成项存在低估或者逆传递现象。修正因子的本质作用是在S-A湍流模型中引入δ=(β-1)P作为源项的修正项，该修正项的取值范围根据预测结果与真实结果的偏差而定。

最后，量化修正模型对湍流摩擦阻力的预测精度，定义湍流摩擦阻力系数的相对误差为$ E_{C_{f}}=\frac{C_{f, \mathrm{SA}}-C_{f, \text { DNS }}}{C_{f, \text { DNS }}}$，其中湍流摩擦阻力系数均由壁面平均流向速度的法向导数计算而得。修正前后的相对误差见表 1。整体而言，在所有雷诺数条件下，采用J_ud为目标函数的预测精度均优于采用J_u为目标函数的计算结果，且与原始S-A模型相比，两者的预测精度均有显著的提高。当Re_τ=180时，摩擦阻力系数的预测精度提高最为显著，误差由-8.41%减小到1.5%以内。随着雷诺数的增大，原始S-A湍流模型的预测误差相应减小，当Re_τ=1 000时，误差为0.48%，但采用J_ud修正后的S-A湍流模型相较于采用J_u的模型可将误差进一步减小至-0.07%。原始S-A湍流模型的预测偏差随着雷诺数的增大而减小，其原因可能是由于S-A湍流模型中涉及的经验常数是根据高雷诺数流动确定的。

湍流摩擦阻力是一个较小的壁面变量，目前针对大型客机的总阻力预测，利用RANS计算要确保阻力误差在1 count以内仍是巨大的挑战，其中最为困难的是湍流摩擦阻力的准确预测，其约占总阻力的50%以上。本文针对槽道湍流这一简单基础构型进行了数据驱动的湍流模型修正，在低雷诺数下其预测精度显著提高，而在较高雷诺数条件下，表 1中的数据显示无修正的湍流模型对摩擦阻力的预测误差已经小于0.48%，通过模型修正，这个误差得到了进一步的减小。该结果表明，基于数据驱动式的湍流模型修正具有较好的应用前景。

3 结论与展望

数据驱动式建模的本质是数据回归和优化问题引入物理知识的先验约束，一定程度上可降低对数据样本个数和质量的依赖，提高模型的预测精度。本文提出了一种基于物理知识的数据驱动式湍流模型修正方法，并以槽道湍流为研究对象验证了该方法的优越性，获得的结论包括：

1) 基于湍流摩擦阻力分解的先验物理知识，建立了一种基于物理知识约束的湍流模型修正方法，将物理知识显性地引入目标函数的设定中，结合伴随优化方法，高效率地获得修正因子的分布，以达到提高预测精度的目的。

2) 以槽道湍流为例，对比了有/无物理知识约束下的湍流模型预测精度，结果表明引入物理知识约束后，可提高对湍流摩擦阻力的预测精度。

本文研究存在如下不足，希望在后期研究中克服：

1) 引入的基于湍流摩擦阻力分解的先验物理知识，适用于无逆压梯度的湍流边界层，而对于有逆压梯度导致的流动分离问题，该物理知识的有效性需要进一步的验证，或引入更为普适的物理知识。

2) 仅开展了基于DNS数据的先验验证，而没有利用机器学习方法，如神经网络、决策树等，建立物理变量与目标函数之间的映射关系，以开展无数据条件下的后验验证。

3) 建立目标函数时，用到了DNS平均速度剖面数据，没有完全避免对DNS数据的依赖，在后期工作中有待进一步完善。