混合湍流模型的参数优化与分析

吴 钦，王国玉，黄 彪

（北京理工大学 机械与车辆工程学院，北京 100081）

混合湍流模型的参数优化与分析

吴 钦，王国玉，黄 彪

（北京理工大学 机械与车辆工程学院，北京 100081）

在以全流场N-S方程为研究对象的空化流动数值计算中，湍流模型在很大程度上影响了对空化流动的精确预测。为了进一步完善针对非定常空化流动数值计算构建的基于密度分域的混合湍流模型（FBDCM），利用代理模型方法对FBDCM模型参数的整体敏感度进行了分析，评价了相关经验系数的扰动对空化流场及水动力特性预测精确度的影响，建议了针对非定常空化流动计算的模型参数的具体取值。结果表明，经验系数C2通过调节FBDCM中FBM模型和DCM模型的影响比重，成为对模型预测精度影响程度最大的因素；通过代理模型优化分析确定经验系数的混合湍流模型，能更好地调整流场内的湍流粘性，精确计算了反向射流区域大尺度涡团的流动，与实验结果更为接近。

混合湍流模型；代理模型；整体敏感度；参数优化

0 引 言

空化现象是关系到水下发射、航空航天及多个工业领域面临的核心关键技术问题[1-2]。近年来，基于N-S方程框架的计算方法，已经成为空化流动数值计算的主流。在以全流场N-S方程为研究对象的空化流动数值计算中，湍流模型在很大程度上影响了对空化流动的精确预测。国内外研究表明，常用的涡粘系数模型不能很好地模拟非定常空化流动[3-5]，为了解决标准RANS模型在处理多重湍流尺度决定的动力流动结构的缺陷，Wu[6]、时素果等[7]采用基于标准k-ε模型的滤波器湍流模型（FBM）、Zhou等[8]采用密度修正模型（DCM）求解空化流动，达到了较好的效果。作者基于对实验现象的分析和总结[9-11]，建立了一种基于密度分域的混合湍流模型（FBDCM），结果表明FBDCM模型可以得到更为显著的非定常流动特征，较好地揭示了空化流场的运动特性。

然而在现有湍流模型的构建中，由于基于不同的物理假设，需要引进多个经验系数，如何确定经验系数，合理地描述空化流动的湍流多尺度和多相可压缩效应，提高非定常空化流动的数值预测精度，仍然是亟须解决的关键问题。Sobol等[12]提出了整体敏感度分析方法，通过对预测结果的合理分解，计算出各自变量的影响程度。建立合适的代理模型是进行整体敏感度分析的基础。Cho等[13]通过不同工程领域的案例研究，介绍了代理模型分析及优化的基本概念和技术方法，指出利用代理模型能有效地评估设计变量对目标的影响程度并对多重竞争性目标进行分析，基于分析结果对原问题进行了科学优化。Tseng等[14]利用代理模型分析及优化技术对基于输运方程的低温空化模型进行了研究，结果表明优化后的模型参数有效提高了低温空化模型的数值预测精度。Goel等[15]发展了多重代理模型的加权平均分析方法，并通过解析函数的代理分析结果阐述了整合代理模型的优越性。

本文利用代理模型分析与优化技术，基于绕Clark-Y翼型空化流动的实验与数值计算结果，对FBDCM模型分域函数进行分析，获得影响预测精确度的主要因素，并基于分析结果对模型参数进行优化，以期给出能够获得精确预测结果的模型参数，完善该计算模型。

1 控制方程与FBDCM模型

1.1 基本方程

采用均质平衡流模型，则Farve平均的Navier-Stokes方程为：

式中：下标i和j分别代表坐标方向，u和p分别为速度和压力，ρm为混合密度，μ和μt分别为混合介质的层流和湍流粘性系数。

1.2 FBDCM模型

标准k-ε模型[4]的控制方程为：

其中：k，ε分别为湍动能和湍流耗散率，Pt为湍动能生成项。模型常数分别为：Cε1=1.44，Cε2=1.92，σε= 1.3，σk=1.0。湍流粘性系数μt定义为：

标准k-ε模型的主要问题在于，它是基于定常平均流动，粘性系数由湍动能和耗散率求解，这暗示着模型反馈的流场主要湍流尺度很大。Johansen[16]提出了基于标准k-ε模型的滤波器湍流模型（FBM），将湍流粘性系数修正为：

其中：fFBM为滤波函数，由滤波尺度（λ）和湍流特征尺度的比值大小决定。

考虑汽液混相的可压缩特性对湍流粘性系数的影响，Coutier-Delgosha[17]提出了基于密度修正的湍流模型（DCM），将湍流系数修正为：

其中：fDCM为密度函数，根据相关文献[17]，n取为10。

针对FBM模型和DCM模型在空化流动模拟应用中的特点，作者提出的FBDCM模型[9]对空化流场基于混合密度的分布进行分域，在不同区域采用不同的湍流粘性修正方式，形成一种基于混合密度分域的湍流模式，充分发挥FBM模型和DCM模型的优势，以捕捉湍流和空化之间的交互作用和动态行为。为了保证湍流粘性系数的光滑过渡，两种湍流粘性系数通过下面分域函数进行桥接，将湍流系数修正为：

其中：C1、C2、C3和C4为模型参数，通过调节模型参数的取值，可影响混合模型（FBDCM）中FBM模型和DCM模型在混相介质中的过渡区间和作用区域。

1.3 几何模型与网格

计算采用与实验[18]一致的Clark-Y翼型空化流动条件，图1给出计算区域及其边界条件。为了较好地匹配翼型头部形状，翼型前端区域采用C型结构化网格划分，并在翼型近壁区域进行了网格加密，如图2所示。采用速度入口、平均静压出口边界条件，流动边界为无滑移固壁，数值模拟工况与实验[18]一致，即σ=0.80，Re=7×105。具体计算细节参见文献[9]。

图1 计算区域与边界条件设置示意图Fig.1 Computational domain and boundary conditions

图2 水翼周围网格示意图Fig.2 Computational grids around the hydrofoil

1.4 分析方法与评价指标

由（10）式可知，混合模型对于空化流场的调节作用对分域函数χ（ρm/ρl）中经验系数（C1、C2等）的取值存在较大的依赖性，目前还难以通过对实验数据的模化来确定不同空化流场特征区域的过渡区间。本文利用基于方差的无量纲整体敏感度分析方法对模型进行评价，说明分域函数χ（ρm/ρl）中相关经验系数的扰动对空化流场及水动力特性预测精确度的影响。

基于代理模型的整体敏感度分析方法，其本质是将目标函数分解成自变量以及反映自变量间相互作用的函数，关键在于建立一个足够精度的数学模型。采用一阶交叉验证方法对不同代理模型的近似精确度进行评价，选取PRESS（Predition error sum of squares）作为评价指标，其表达式如下：

其中：Ns为建立代理模型选取的点的个数，yi为xi点的数值计算结果，为利用除xi点外其他样本点建立代理模型在xi点的拟合值。

敏感度指标定义如下：

其中：Vi表示单个自变量的偏方差，Viz表示多个自变量组合的偏方差，总方差V=ΣVi+Vij+…+V1…N。Si为局部敏感度，反映独立自变量对结果的影响；为全局敏感度，综合考虑了自变量的局部敏感度和其他自变量对其的相互作用。

2 FBDCM模型优化结果与讨论

2.1 FBDCM模型的参数分析

综合考虑时空复杂度，将分域函数χ（ρm/ρl）中经验系数C1、C2、C3和C4作为模型评价的自变量。为实现分域函数χ（ρm/ρl）在混相介质区间内的合理过渡，根据初步的试算分析确定自变量取值范围如表1所示，代理模型的分析目标定义为云状空化工况下，水翼的时均升力系数（CL）和阻力系数（f）实验值与数值预测值之差，分别记为 CL_diff和 f_diff。

结合面心立方法（FCCD）和拉丁超立方法（LHS），在自变量空间内共生成70个样本点，然后对各样本点进行数值实验，根据数值实验结果，在标准化的变量空间分别采用多项式响应面方法（PRS）、Kriging方法（KRG）、辐射神经元网络方法（RBNN）以及加权平均值近似方法（WAS）[19]建立目标函数与自变量之间的近似关系，即代理模型。表2给出了由不同代理模型拟合得到的PRESS指标分布情况。综合四种模型拟合结果可以看出，WAS模型结果能获得最小的误差分布，拟合效果明显优于其它模型，PRS模型次之。因此，选择WAS模型的拟合结果对目标变量做进一步分析。

表1 自变量取值范围Tab.1 The range of the variables

表2 不同代理模型拟合误差分布Tab.2 The error estimates for different surrogate models

利用WAS模型对绕Clark-Y翼型云状空化流动的数值计算结果进行拟合，图3反映了自变量对目标变量的敏感度分布情况。从图中不难发现，相同自变量对同一目标变量的局部敏感度与全局敏感度分布趋势基本一致，但局部敏感度指标所占比重的具体数值明显小于全局敏感度，说明自变量间的交互作用不容忽视。

图4给出了自变量对目标变量的整体敏感度指标分布。在对目标变量CL_diff和f_diff的敏感度分析中，各自变量敏感程度分布趋势基本一致，由大到小依次为C2、C3、C4和C1，不同目标变量的整体敏感度分析结果均说明模型预测精度主要受经验系数C2的影响，受经验系数C1的影响相对较小。

图3 自变量对目标变量的敏感度分布图Fig.3 The distribution of the Main Sensitivity and the Global Sensitivity

图4 整体敏感度分析饼状图Fig.4 Pie-chart of the Global Sensitivity analysis

图5 目标变量随自变量C1、C2的变化趋势图Fig.5 The variation of the objectives with design variables

考虑到自变量间相互作用对敏感性分析结果的影响，通过目标变量的变化趋势进一步分析其敏感程度。图5给出了目标变量随自变量C1、C2的变化情况。从图5可以看出，目标变量CL_diff和f_diff沿Y轴方向随C2的增加分别有减小和增大的趋势，而沿X轴方向随着C1的增大，变化趋势不明显，且目标变量f_diff的变化范围较CL_diff大，说明目标变量CL_diff和f_diff均对C2的变化较敏感，且C2对f_diff的影响程度较对CL_diff的影响程度大，这与图4所述敏感度分布情况一致。

2.2 FBDCM模型的参数优化

为了更直观地体现与自变量空间相对应的目标函数分布情况，在自变量空间中随机生成均匀分布的104个点，WAS模型的拟合结果在目标空间的分布如图6所示。图中红色区域为最优解，即CL_diff和f_diff同时取最小值。由于设计变量与目标函数是一一对应的，因此由最优解可缩小自变量的取值范围，从而实现对设计空间的优化。表3给出了自变量取值范围的优化结果。其中，空间优化程度定义为优化后取值范围的减小值与优化前取值范围之比。不难发现，经过代理模型优化后，设计空间大大缩小，从而减小了进一步分析的时空复杂度。

表3 自变量取值范围优化分析Tab.3 The optimization results of the design variables

由前述分析可知，在当前计算工况下，模型参数C2对模型预测精度的影响最大，其它参数在敏感度分析中仅起到参考作用，而非作为优化参数。根据自变量的优化分析结果，固定C1=4、C3=0.6和C4=0.2，图7给出采用不同C2值（C2=0.1、0.2和0.3）条件下，混合模型分域函数χ（ρm/ρl）的分布图。结合（9）、（10）式可知，通过改变分域函数χ（ρm/ρl）中经验系数C2的取值可以调节混合模型（FBDCM）中FBM模型和DCM模型的影响比重，即相同含汽量区域，当χ（ρm/ρl）较小时，对湍流粘性采用密度函数的修正方式（DCM），当χ（ρm/ρl）较大时，则采用空间滤波的修正方式（FBM）。

图6 目标变量分布图Fig.6 The distribution of the objectives

2.3 FBDCM模型的计算结果分析

在代理模型的优化基础上，对绕水翼的非定常云状空化发展过程进行进一步模拟分析。图8给出了应用不同C2取值的FBDCM模型计算所得空穴形态与实验观测结果的对比。由实验结果观察可知，附着空穴的发展经过了透明状空穴线性增长至最大长度和反向射流向水翼前缘发展两个阶段，当反向射流达到水翼前缘时，空穴突然断裂，附着在翼型上的空穴长度骤减。由数值结果观察可知，在空穴发展的前一阶段，不同C2取值的数值计算结果都合理地预测了附着空穴的非定常变化，其中，C2=0.3时混合模型预测得到的空化区域相对较大，最大空穴长度也较长；在空穴发展的后一阶段，C2=0.1和0.2时混合模型能捕捉到空泡的大尺度脱落现象，与实验较为一致，而C2=0.3时混合模型过早地预测了空穴的断裂和脱落过程。这是由于在云状空化中，反向射流的产生和发展会直接导致空泡团的断裂和脱落现象，在这部分区域内，利用滤波尺度的湍流粘性修正方式应起到主导作用，而反向射流往往发生在含汽量较小的近壁水汽混相区域。结合图7所示，当C2= 0.3时，相同含汽量区域分域函数χ（ρm/ρl）值较小，即混合模型（FBDCM）中FBM模型比重较小，因此低估了多尺度涡团脱落的湍流信息。

表4给出了不同C2取值条件下数值预测得到的水翼升力系数、升力特征频率与实验测量值的对比。对比发现，当C2=0.3时，模型预测的升力特征频率偏大，当C2=0.1和0.2时，模型预测的结果基本一致。与图8所示结论一致，C2=0.1和0.2时分域函数χ（ρm/ρl）值较C2=0.3时大，即在近壁混合密度较大的反向射流区域，主要采用滤波函数而非密度函数对湍流粘性进行修正，因此对大尺度涡团的流动捕捉更精确，由此预测得到的空泡团非定常发展情况及升力的变化情况与实验结果更为接近。

图7 不同C2取值的混合模型分域函数分布示意图Fig.7 Distribution of hybrid function χ with different values of C2for the hybrid FBDCM model

图8 不同C2取值的模型预测得到的空穴形态随时间的变化过程Fig.8 Time evolution of cavity shape obtained from experiment and calculation

表4 不同湍流模型预测得到的水翼升力系数的时均值和特征频率Tab.4 Comparisons of predicted and measured mean value of lift coefficient(Cl)and frequency f

为了进一步说明不同C2取值对空化流场结构预测精度的影响，提取图9所示的一组特征截面，图10给出了绕翼型的时均流场分布情况，并与实验结果进行了对比。其中，横坐标为相应特征截面上沿主流方向的时均速度分布，纵坐标为沿Y轴方向的位置，c为水翼弦长，Uref为参考速度。由图观察可知，在x/c=0.2与x/c=0.4位置处，不同C2取值的模型预测得到的时均速度变化趋势基本一致，空穴外部的平均速度大小接近主流平均速度，空穴内部速度值迅速减小，呈现出很强的剪切层结构；在x/ c=0.6与x/c=0.8位置处，不同C2取值的模型预测得到的时均速度分布均能有效地反映近壁区域的反向射流现象，速度剪切层区域较前一位置更为显著，其中C2=0.1时模型预测得到的时均速度值小于C2= 0.2和0.3时的预测结果；在x/c=1.0与x/c=1.2的翼型尾部流场中，不同C2取值的模型预测得到的时均速度呈大幅变化趋势，翼型尾端汽液两相区的大规模非定常运动加强了空穴区域的动量交换，产生了空泡团的旋涡分离现象，从而对该区域内的时均速度分布产生了较大的影响，其中C2=0.1和0.2时模型预测得到的反向射流区域时均速度明显大于C2=0.3时的预测结果，与实验结果更为接近，这是由于前者通过增大相同含汽量区域的分域函数χ（ρm/ρl）值，增大了混合模型中FBM模型的比重，采用滤波函数对湍流粘性进行修正，对水翼尾缘的湍流运动采用直接求解的方法，精确计算了大尺度涡团的流动，使得时均速度分布具有更好的空间分辨率。综合以上分析，采用C2=0.2作为混合模型的分域函数χ（ρm/ρl）经验系数。

图9 特征截面位置示意图Fig.9 Profile location

图10 不同C2取值的模型预测得到的流场时均速度分布与实验结果的对比Fig.10 Time-average u-velocity of experimental and simulation results

3 结 论

本文针对绕Clark-Y翼型云状空化流动进行的数值计算结果和实验结果，应用代理模型分析及优化方法对基于密度分域的混合湍流模型的参数取值进行了分析，评价了相关经验系数对空化流场及水动力特性预测精确度的影响。主要结论如下：

（1）混合模型对于空化流场的调节作用对分域函数χ（ρm/ρl）中经验系数的取值存在较大的依赖性，其预测精度主要受到经验系数C2的影响，且经验系数C2通过调节混合模型（FBDCM）中的FBM模型和DCM模型的影响比重，即对相同含汽量区域不同湍流粘性的修正方式影响模型的预测精度；

（2）C2=0.2时混合湍流模型（FBDCM）中的分域函数χ（ρm/ρl）值较C2=0.3时大，后者低估了云状空化发展过程中多尺度涡团脱落的湍流信息，而前者在近壁混合密度较大的反向射流区域，主要采用滤波函数对湍流粘性进行修正，更精确地计算了大尺度涡团的流动，预测得到的空泡团非定常发展情况及升力的变化情况与实验结果更为接近，时均速度分布具有更好的空间分辨率；

（3）通过代理模型优化分析，确定混合模型分域函数χ（ρm/ρl）中经验系数的取值分别为C1=4、C2= 0.2、C3=0.6和C4=0.2。

[1]Yu X X,Wang Y W,Huang C G,Wei Y P,Fang X,Du T Z,Wu X C.Experiment and simulation on air layer drag reduction of high-speed underwater axisymmetric projectile[J].European Journal of Mechanics B-Fluids,2015,52:45-54.

[2]Ji B,Luo X,Wang X,Peng X,Wu Y,Xu H.Unsteady numerical simulation of cavitating turbulent flow around a highly skewed model marine propeller[J].Journal of Fluids Engineering,2011,133:1-8.

[3]Wang G,Ostoja-Starzewski M.Large eddy simulation of a sheet/cloud cavitation on a NAVA0015 hydrofoil[J].Applied Mathematical Modeling,2007,31(3):417-447.

[4]Launder B E,Spalding D B.The numerical computation of turbulent flows[J].Computational Methods in Applied Mechanics and Engineering,1974,3:269-289.

[5]Ji B,Luo X W,Arndt R E A,Peng X X,Wu Y L.Large eddy simulation and theoretical investigation of the transient cavitating vortical flow structure around a NACA66 hydrofoil[J].International Journal of Multiphase Flow,2015,68:121-134.

[6]Wu J,Wang G,Shyy W.Time-dependent turbulent cavitating flow computations with interfacial transport and filter based models[J].International Journal for Numerical Methods for Fluids,2005,49:739-761.

[7]时素果,王国玉,余志毅,张敏弟.FBM湍流模型在非定常通气超空化流动计算中的评价与应用[J].船舶力学,2012, 16(10):1099-1106. Shi S G,Wang G Y,Yu Z Y,Zhang M D.Evaluation of the filter-based turbulence model(FBM)for computation of unsteady ventilated-supercavitating flows[J].Journal of Ship Mechanics,2012,16(10):1099-1106.

[8]Zhou L J,Wang Z W.Numerical simulation of cavitation around a hydrofoil and evaluation of a RNG k-ε model[J].Journal of Fluids Engineering,2008,130(1):1-7.

[9]Zhao Y,Wang G,Huang B.A curvature correction turbulent model for computations of cloud cavitating flows[J].Engineering Computations,2016,33(1):202-216.

[10]Hu C,Wang G,Chen G,Huang B.A modified PANS model for computations of unsteady turbulence cavitating flows[J]. Science China Physics,Mechanics&Astronomy,2014,57(10):1967-1976.

[11]Huang Biao,Wang Guoyu.A density modify based cavitation model for time dependent turbulent cavitating flow computations[J].Chinese Science Bulletin,2011,56(19):1985-1992.

[12]Sobol I M.Sensitivity analysis for nonlinear mathematical models[J].Mathematical Modeling&Computational Experiment,1993,1(4):407.

[13]Cho Y C,Du W,Gupta A,Tseng C,Sastry A M,Shyy W.Surrogate-based modeling and dimension-reduction techniques for thermo-fluid&energy systems[J].Proceeding of the ASME/JSME 2011 8th Thermal Engineering Joint Conference, 2011,3.

[14]Tseng C,Shyy W.Modeling for isothermal and cryogenic cavitations[J].Int.J Heat Mass Transfer,2010,53:513-525.

[15]Goel T,Haftka RT,Shyy W,et al.Ensemble of surrogates[J].Journal of Structural and Multidisciplinary Optimization, 2007,33(3):199-216.

[16]Johansen S T,Wu J,Shyy W.Filter-based unsteady RANS computations[J].International Journal of Heat and Fluid Flow, 2004,25(1):1-21.

[17]Coutier-Delgosha O,Fortes-Patella R,Reboud J L.Evaluation of the turbulence model influence on the numerical simulations of unsteady cavitation[J].Journal of Fluids Engineering,2003,125(1):38-45.

[18]张 博,王国玉,张淑丽.修正的RNG k-ε模型在云状空化流动计算中的应用评价[J].北京理工大学学报,2008,28 (12):1065-1069. Zhang B,Wang G Y,Zhang S L.Evaluation of a modified RNG k-ε model for computations of cloud cavitating flows[J]. Transctions of Beijing Institute of Technology,2008,28(12):1065-1069.

[19]Shyy W,Cho Y C,Du W,et al.Surrogate-based modeling and dimension reduction techniques for multi-scale mechanics problems[J].Acta Mechanica Sinica,2011,27(6):845-865.

Parameter optimization and analysis of a Filter-Based Density Correction model

WU Qin,WANG Guo-yu,HUANG Biao
(School of Mechanical and Vehicular,Beijing Institute of Technology,Beijing 100081,China)

Turbulence model has great effect on the prediction precision of cavitating flows.A systematic investigation based on the surrogate modeling techniques is conducted to assess and improve the performance of a filter-based density correction model(FBDCM),which is suggested for the computations of unsteady cavitating flows.Based on the surrogate model,global sensitivity analysis is conducted to assess the role of model parameters.The results show that,the performance of the filter-based density correction model is affected more by model parameter C2,which is used to regulating the proportion of FBM model and DCM model based on local mixture density,so that it will help to limit the over-prediction of the turbulent eddy viscosity in the cavitating regions on the foil wall and in the wake.Based on the recommended model parameter values,better prediction of the cavitating flows can be attained.

Filter-Based Density Correction Model;surrogate model; global sensitivity;parameter optimization

O35

：Adoi:10.3969/j.issn.1007-7294.2016.07.001

1007-7294（2016）07-0789-10

2016-01-08

国家自然科学基金项目资助（50979004）；国家重点自然科学基金项目资助（51239005）

吴 钦（1989-），女，博士研究生，E-mail：wuqin919@163.com；王国玉（1961-），男，教授，博士生导师。