張錫，楊希祥，張為華，武澤平，趙海龍，王鵬宇

國防科技大學(xué)空天科學(xué)學(xué)院，湖南長沙 410073

代理模型在空間科學(xué)、結(jié)構(gòu)分析、多學(xué)科設(shè)計優(yōu)化等諸多工程優(yōu)化問題中扮演著越來越重要的作用。計算流體力學(xué)建模中耗時持久的計算可用代理模型替代，顯著提高計算效率。此外，代理模型可集成全局尋優(yōu)、數(shù)值噪聲過濾、并行設(shè)計優(yōu)化和數(shù)值仿真于一體，提高工程設(shè)計效率【1-8】。代理模型又稱為響應(yīng)面模型、近似模型、元模型等。在諸多代理模型中，Kriging 模型因在分析和預(yù)測多維、高度非線性問題時顯示出了優(yōu)良的準確性而受到越來越多學(xué)者的青睞【9-12】?；贙riging 模型的改進模型也越來越多，CoKriging模型【13-19】、梯度增強GEKriging模型等【20-23】。在提高代理模型建模效率方面，CoKriging模型結(jié)合不同精度的樣本數(shù)據(jù)提供了一種非常有吸引力的選擇。本文對CoKriging模型進一步改進構(gòu)建了ICoKriging模型（improved CoKriging），并通過數(shù)值試驗和超聲速優(yōu)化問題檢驗該模型的有效性。

1 改進CoKriging模型

1.1 模型描述及基本假設(shè)

對一個m維優(yōu)化問題，CoKriging模型是僅設(shè)計高可信度和低可信度抽樣。這里，對此設(shè)計進行改進，

其中n1、n2、n3是高中低可信度樣本點數(shù)。（S1，y1），（S2，y2），（S3，y3）是各采樣空間的樣本點，我們將其重新表述為

1.2 ICoKriging模型預(yù)測值及均方誤差

假設(shè)最后模型的輸出可用y1，y2，y3的線性組合近似，于是ICoKriging 模型在x對y1(x)的預(yù)測值可表述為

通過矩陣求導(dǎo)方法，可得到如下方程組

作如下變換

可得

令

1.3 相關(guān)函數(shù)

在ICoKriging 模型中，需要設(shè)置相關(guān)函數(shù)R()ij的具體形式，該函數(shù)反映了不同精度樣本對模型輸出結(jié)果的影響，類似于Kriging 模型中對于任意兩點xi?，xj?，R僅僅是xi?，xj?距離的函數(shù)R=R(‖ ‖xi?-xj?)，一般可用如下形式的函數(shù)

其中εk=θk|xik-xjk|，k∈1，2，…，m. 由三次樣條函數(shù)生成的相關(guān)矩陣條件數(shù)相比高斯函數(shù)表現(xiàn)更好。

1.4 模型超參數(shù)的優(yōu)化

θ沒有解析解，可通過遺傳算法等進化算法求得θ的解。

1.5 相關(guān)矩陣求逆

2 算例分析

2.1 數(shù)值算例

為驗證ICoKriging 模型的性能，我們用一些函數(shù)進行測試，并對Kriging、CoKriging 和ICoKriging 的插值結(jié)果進行分析比較。此處引用文獻［13］提供的經(jīng)典一維測試函數(shù)y=(6x- 2)2sin(12x- 4)，x∈[0，1]，高、中、低三種可信度的樣本分別為S1=(0.0，0.6，1)，S2=(0.1，0.4，0.5)，S3=(0.3，0.8，0.9)。

圖1展示了不同插值模型對樣本的擬合效果以及插值結(jié)果和真實函數(shù)的比較。圖1中實心紅色方塊為高精度樣本S1，藍色實心三角形為中精度樣本S2，綠色實心三角形為低精度樣本S3。Kriging模型僅使用三個高精度的樣本S1，插值見圖1（a）；CoKriging 模型使用了S1、S2，插值結(jié)果見圖1（b）；ICoKriging 模型則綜合了S1、S2、S3三種樣本，插值結(jié)果見圖1（c）部分。圖1（d）是三種插值模型與真實函數(shù)的比較。從圖中可以看出在x∈[0，0.6]和x∈[0.6，1]上，Kriging 模型與真實函數(shù)y=(6x- 2)2sin(12x- 4)在絕對誤差和變化趨勢上均存在較大誤差；CoKriging 模型使用S1、S2后插值結(jié)果與真實函數(shù)仍有不小誤差，加入樣本S2后得到改善的是在x∈[0，0.6]和x∈[0.6，1]上其變化趨勢與真實函數(shù)相比基本相同，但是絕對誤差仍然較大；ICoKriging 模型綜合樣本S1、S2、S3后，插值結(jié)果與真實函數(shù)的絕對誤差進一步降低，整個設(shè)計空間上的插值結(jié)果與真實函數(shù)變化趨勢基本相同。

圖1 Kriging，CoKriging和ICoKriging插值方法的比較Fig.1 Comparison of Kriging,CoKriging and ICoKriging interpolation methods

2.2 超聲速減阻問題

減阻是超音速飛行器的重要問題之一。減阻桿（Aerospike）和減阻盤（Aerodisk）搭配使用作為減阻的重要方法之一，其減阻效果好，目前已在一些裝備中得到應(yīng)用。對不同減阻桿長度和平面圓柱減阻盤半徑與半球形頭部（Hemispherical Body）的組合在超聲速流場中的流動情況進行數(shù)值模擬，分析減阻桿減阻機理以及減阻盤半徑和減阻桿長度參數(shù)對減阻效果的影響。研究發(fā)現(xiàn)，減阻盤半徑和減阻桿長度在一定范圍內(nèi)的增大有利于減阻。

2.2.1 計算模型本文選取的幾何模型采用圓球頭部加減阻桿的外形，頭部半徑R= 100 mm。減阻桿的結(jié)構(gòu)形式采用桿加頂端減阻盤，主要設(shè)計參數(shù)為減阻桿桿長d和減阻盤的半徑r，減阻盤厚度固定為t= 10 mm。由于整個幾何外形為軸對稱的，因此在計算中采取了軸對稱計算假設(shè)。研究表明，對于流動較為穩(wěn)定的情況下，雖然在再循環(huán)區(qū)出現(xiàn)了剪切層不對稱的情況，但整個流場體現(xiàn)出了明顯的軸對稱現(xiàn)象，因此使用軸對稱假設(shè)在流動穩(wěn)定的情況下是可以接受的。半球形頭部的縱向截面外形如圖2。圓球頭部及減阻桿的網(wǎng)格劃分如圖3。鑒于各類湍流模型對減阻桿飛行器數(shù)值計算結(jié)果影響的分析，采用了SST k - w 湍流模型［25］，邊界條件采用壓力遠場條件、理想氣體假設(shè)，氣體參數(shù)設(shè)置為標準大氣海平面參數(shù)，其他物理條件馬赫數(shù)Ma= 2，攻角α= 0，壓強P= 101325 Pa，溫度T= 288.15 K。

圖2 計算域和幾何模型：半球形頭部Fig.2 Computational domain and geometric model:hemispherical head

圖3 流場網(wǎng)格劃分Fig.3 Meshing of flow field

目前，減阻盤半徑與彈頭部半徑比一般為10%～60%，減阻桿長度與彈頭部半徑比一般為100%～400%，則設(shè)計減阻盤半徑10 ≤r≤60，減阻桿長度100 ≤d≤400。解如下優(yōu)化問題

其中Cd是減阻桿、減阻盤和彈身頭部總的阻力系數(shù)。

流體問題的數(shù)值仿真結(jié)果受多種因素的影響，其中網(wǎng)格數(shù)量是重要的影響因素之一。對文中的減阻問題，在80 000、120 000、240 000、360 000、700 000 網(wǎng)格數(shù)下計算阻力系數(shù)，迭代收斂情況如圖4所示。

圖4 不同網(wǎng)格數(shù)下阻力系數(shù)Fig.4 Drag coefficient of different grid numbers

由圖4可見，網(wǎng)格數(shù)量較小時，計算過程可快速收斂，但其仿真結(jié)果不太理想。隨著網(wǎng)格數(shù)的逐漸增多，仿真結(jié)果逐漸接近，50 000 步以后120 000 個以上網(wǎng)格數(shù)的仿真模型計算結(jié)果基本重合。為建立ICoKriging 模型，需要設(shè)計三種不同可信度的樣本。劃分80 000、120 000、240 000 個網(wǎng)格，使用CFD 模型計算不同減阻桿長d和減阻盤半徑r對應(yīng)的阻力系數(shù)Cd，三種網(wǎng)格分別對應(yīng)高精度（S1）、中精度（S2）、低精度（S3），三種樣本S1，S2，S3樣本量分別為25、36、49。S1的Cd-d（固定r）、Cd-r（固定d）樣本圖見圖5。

圖5 阻力系數(shù)隨減阻桿長度、減阻盤半徑的變化趨勢Fig.5 The relationship between drag coefficient and drag-reducing rod length，drag-reducing disk radius

2.2.2 計算結(jié)果將減阻桿長d、減阻盤半徑r、阻力系數(shù)Cd歸一化后分別建立模代理模型ICoKriging、CoKriging、Kriging，使用遺傳算法對代理模型尋優(yōu)、加點、建模、再次尋優(yōu)，直至滿足收斂條件。數(shù)值計算結(jié)果見表1。

表1 基于代理模型的優(yōu)化結(jié)果Table 1 Optimization results based on Surrogate model

基于以上三種優(yōu)化結(jié)果，重新設(shè)置桿長和半徑，求解CFD模型得到頭部阻力系數(shù)，發(fā)現(xiàn)ICoKriging模型與CFD 模型求解結(jié)果最為吻合，速度分布與壓力分布見圖6。試驗證明基于ICoKriging 模型的優(yōu)化結(jié)果更準確。頭部阻力系數(shù)與減阻桿長、圓盤半徑有如下關(guān)系：半徑較小時，隨著半徑增加，阻力系數(shù)減少；當(dāng)半徑進一步增大，減阻率達到飽和，再將半徑增加會致使阻力系數(shù)增加；桿長的作用效果與半徑類似，在一定范圍內(nèi)桿長增加，阻力系數(shù)減少；當(dāng)桿長達到一定值，再將減阻桿桿長增加，阻力系數(shù)會增加。

圖6 速度和壓力分布Fig.6 Speed and pressure distribution

3 結(jié) 論

本文提出了一種新的建立變可信度代理模型的方法。在CoKriging 模型的基礎(chǔ)上將CoKriging 模型中的高精度樣本、低精度樣本擴展至高、中、低三種精度不同的樣本，建立ICoKriging模型。該模型在CoKrig?ing 模型基礎(chǔ)上重新定義了插值權(quán)重系數(shù)，并構(gòu)建了新的相關(guān)矩陣。通過重定義權(quán)重系數(shù)可以減少ICoKriging模型符號表達復(fù)雜性，模型預(yù)測值在表達形式上與CoKriging保持一致，在理解上和具體實現(xiàn)上并無難處。同Kriging 模型、CoKriging 模型，建立ICoKriging 前，需要求解模型超參數(shù)，這本身即是一個優(yōu)化問題，文中對超參數(shù)的求解提出了簡單可行的算法。最后通過解析函數(shù)和超聲速飛行器頭部減阻問題對Kriging、CoKriging、ICoKriging 三種模型進行了比較，試驗表明ICoKriging 模型預(yù)測更為準確，優(yōu)化能力更強，在氣動數(shù)據(jù)校驗和外形優(yōu)化方面具有更多的應(yīng)用前景。