M6機翼
M6是ONERA設計的一種機翼模型。該模型在跨聲速條件下進行了一系列風洞試驗。試驗馬赫數在之間,攻角區間為度,雷諾數Re(參考長度為平均氣動弦長c)約為。儘管M6機翼幾何外形簡單,但是其涉及的跨聲速流動卻十分複雜,包含局部超音速流動、激波和邊界層分離等。M6機翼具備三維可壓縮流動的典型特徵,因此被大量論文選為CFD代碼的驗證算例。本文以M6為測試算例,檢驗OpenFOAM V4在可壓縮流場模擬方面的計算效率和計算精度。
圖 1 M6機翼風洞試驗模型
M6是一種無扭曲的後掠機翼,其基本翼型為ONERA D section對稱翼型。M6機翼幾何外形和參數見圖2。試驗時,在7個展向截面上佈置了壓力感測器,測得的壓力數據可用於與計算結果進行對比。
圖 2 M6機翼幾何外形及參數
圖 3 M6機翼網格
本次計算採用NASA網站上公開發布的一種C型結構化網格。(https://www.grc.nasa.gov/www/wind/valid/m6wing/m6wing01/m6wing.x.fmt)網格包含在名為m6wing.x.fmt的PLOT3D網格文件(formatted, multi-zone, whole, three-dimensional)中。它由4個網格塊組成,表1列出了各塊網格的節點分佈,總共316932個網格點。
表 網格分區及節點分佈
Zone2的J1表面覆蓋機翼的下表面,Zone3的J1表面覆蓋機翼的上表面。所有Zone的K1表面都在對稱平面上, KMAX邊界在它們包圍翼尖時彼此耦合。Zone1的I1表面和Zone4的IMAX表面是出口邊界。所有區域的JMAX表面都在遠場邊界上。圖3顯示了機翼表面和尖端的網格。網格圍繞機翼表面聚集以解析邊界層。
本次計算需要將plot3d格式網格,轉換為OpenFOAM求解器能夠讀取的網格存儲格式。我們採用Pointwise V18.1 R1軟體進行格式轉換。具體步驟如下:
(1)用文本編輯器(推薦採用notepad++)打開m6wing.x.fmt,將其中的逗號全部替換為空格,將文件保存為m6wing.x;
(2)打開Pointwise V18.1 R1軟體,導入網格,具體操作見下圖;
(3)將求解器設置為OpenFOAM,並設置邊界條件;
(4)對網格進行縮放。
(5)輸出結果。
結合定常求解器rhoSimpleFoam和非定常求解器rhoCentralFoam計算M6機翼三維流場。首先,採用rhoSimpleFoam計算NS方程,獲得流場初始解。然後,再用rhoCentralFoam計算URANS方程(採用SST湍流模型),待流場穩定後,即獲得最終流場。
rhoSimpleFoam是OpenFOAM中一個基於壓力的可壓縮穩態求解器。所有OpenFOAM求解器設置都需要準備三種文件,分別是0文件,constant文件和system文件。0文件中包含流場中各物理量的初始值和邊界條件。constant文件中包含網格信息和物性參數。system文件中包含求解器相關參數。
(1)0文件
0文件夾中需要準備p,T,U,alphat和nut五個文件。這五個文件的內容可參考OpenFOAM自帶的相關案例。修改好的文件內容如下:
p文件如下所示。其中pOut為來流靜壓,遠場採用freestreamPressure條件,對稱面為symmetryPlane條件,壁面為zeroGradient條件。
/*--------------------------------*- C++ -*----------------------------------* | ========= | | | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \ / O peration | Version: v1806 | | \ / A nd | Web: www.OpenFOAM.com | | \/ M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; object p; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // pOut 8.353085E+04; dimensions [1 -1 -2 0 0 0 0]; internalField uniform $pOut; boundaryField { farfield { type freestreamPressure; freestreamValue uniform $pOut; } sym { type symmetryPlane; } wall { type zeroGradient; } #includeEtc "caseDicts/setConstraintTypes" }
// ************************************************************************* // T文件如下所示。其中Tinlet為來流靜溫,遠場採用inletOutlet條件,對稱面為symmetryPlane條件,壁面為zeroGradient條件。 /*--------------------------------*- C++ -*----------------------------------* | ========= | | | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \ / O peration | Version: v1806 | | \ / A nd | Web: www.OpenFOAM.com | | \/ M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; object T; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Tinlet 2.629383E+02; dimensions [0 0 0 1 0 0 0]; internalField uniform $Tinlet; boundaryField { farfield { type inletOutlet; inletValue uniform $Tinlet; value $inletValue; } sym { type symmetryPlane; } wall { type zeroGradient; } #includeEtc "caseDicts/setConstraintTypes" } // ************************************************************************* // U文件如下所示。其中Uinlet為來流速度矢量,遠場採用freestream條件,對稱面為symmetryPlane條件,壁面為noSlip條件。 /*--------------------------------*- C++ -*----------------------------------* | ========= | | | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \ / O peration | Version: v1806 | | \ / A nd | Web: www.OpenFOAM.com | | \/ M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volVectorField; object U; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Uinlet (2.724789E+02 1.456615E+01 0.000000E+00); dimensions [0 1 -1 0 0 0 0]; internalField uniform $Uinlet; boundaryField { farfield { type freestream; freestreamValue uniform $Uinlet; value uniform $Uinlet; } sym { type symmetryPlane; } wall { type noSlip; } #includeEtc "caseDicts/setConstraintTypes" }
// ************************************************************************* // alphat文件如下所示。 /*--------------------------------*- C++ -*----------------------------------* | ========= | | | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \ / O peration | Version: v1806 | | \ / A nd | Web: www.OpenFOAM.com | | \/ M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; object alphat; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [1 -1 -1 0 0 0 0]; internalField uniform 0; boundaryField { farfield { type calculated; value uniform 0; } sym { type symmetryPlane; } wall { type compressible::alphatWallFunction; value uniform 0; } #includeEtc "caseDicts/setConstraintTypes" }
// ************************************************************************* // nut文件如下所示。 /*--------------------------------*- C++ -*----------------------------------* | ========= | | | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \ / O peration | Version: v1806 | | \ / A nd | Web: www.OpenFOAM.com | | \/ M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; object nut; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 -1 0 0 0 0]; internalField uniform 0; boundaryField { farfield { type calculated; value uniform 0; } sym { type symmetryPlane; } wall { type nutkWallFunction; type nutkWallFunction; value uniform 0; } #includeEtc "caseDicts/setConstraintTypes" } // ************************************************************************* //
(2)constant文件
constant文件夾中需要準備thermophysicalProperties,turbulenceProperties和polyMesh網格文件。其中,thermophysicalProperties和turbulenceProperties的編寫可參考OpenFOAM自帶的相關案例,polyMesh網格文件直接由PointWise軟體輸出。修改好的文件內容如下:
thermophysicalProperties主要包含流體的屬性參數,一般不作修改,文件內容如下所示。
/*--------------------------------*- C++ -*----------------------------------* | ========= | | | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \ / O peration | Version: v1806 | | \ / A nd | Web: www.OpenFOAM.com | | \/ M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object thermophysicalProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // thermoType { type hePsiThermo; mixture pureMixture; transport const; thermo hConst; equationOfState perfectGas; specie specie; energy sensibleInternalEnergy; }
mixture // air at room temperature (293 K) { specie { nMoles 1; molWeight 28.9; } thermodynamics { Cp 1005; Hf 0; } transport { mu 1.82e-05; Pr 0.71; } } // ************************************************************************* // turbulenceProperties文件主要對湍流模型進行設置,採用的湍流模型為liminar,文件內容如下所示。 /*--------------------------------*- C++ -*----------------------------------* | ========= | | | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \ / O peration | Version: v1806 | | \ / A nd | Web: www.OpenFOAM.com | | \/ M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object turbulenceProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // simulationType laminar; // ************************************************************************* //
(3)system文件
system文件夾中需要準備fvSchemes,fvSolution,controlDict,fvOptions,residuals等文件。其中,前三個文件是必須提供的文件。上述文件的編寫均可參考OpenFOAM自帶的相關案例。修改好的文件內容如下:
fvSchemes文件主要對求解器採用的數值離散格式進行設定,由於M6涉及激波等超聲速流動,且rhoSimpleFoam這一求解器穩態求解器對數值格式較為敏感,對流項選用1階迎風格式。具體設置如下所示。
/*--------------------------------*- C++ -*----------------------------------* ========= | \ / F ield | OpenFOAM: The Open Source CFD Toolbox \ / O peration | Website: https://openfoam.org \ / A nd | Version: 6 \/ M anipulation | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes { default steadyState; }
gradSchemes { default Gauss linear; limited cellLimited Gauss linear 1; grad(U) $limited; grad(k) $limited; grad(omega) $limited; grad(epsilon) $limited; }
divSchemes { default none; div(phi,U) bounded Gauss upwind; turbulence bounded Gauss upwind; energy bounded Gauss upwind; div(phi,k) $turbulence; div(phi,omega) $turbulence; div(phi,epsilon) $turbulence; div(phi,e) $energy; div(phi,K) $energy; div(phi,Ekp) $energy; div(phid,p) Gauss upwind; div((phi|interpolate(rho)),p) bounded Gauss upwind; div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; }
laplacianSchemes { default Gauss linear limited corrected 0.5; }
interpolationSchemes { default linear; }
snGradSchemes { default limited corrected 0.5; }
// ************************************************************************* // fvSolution文件主要迭代過程及收斂準則等方面進行設置,設置參數如下所示。其中,pRefPoint為壓力參考點坐標位置,一般設置在遠場邊界附近。pRefValue為壓力參考值,一般設置為來流靜壓。 /*--------------------------------*- C++ -*----------------------------------* ========= | \ / F ield | OpenFOAM: The Open Source CFD Toolbox \ / O peration | Website: https://openfoam.org \ / A nd | Version: 6 \/ M anipulation | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers { p { solver PCG; preconditioner DIC; tolerance 1e-6; relTol 0.001; }
"(U|e)" { solver PBiCG; preconditioner DILU; tolerance 1e-6; relTol 0.0001; }
"(k|omega|epsilon)" { solver PBiCG; preconditioner DILU; tolerance 1e-6; relTol 0.0001; minIter 3; } }
SIMPLE { residualControl { p 2e-5; U 1e-7; "(k|omega|epsilon|e)" 1e-6; } pRefPoint (-6.528 -0.261 4.41); pRefValue 8.353085E+04; nNonOrthogonalCorrectors 0; pMinFactor 0.1; pMaxFactor 2; rhoMin rhoMin [1 -3 0 0 0] 0.001; rhoMax rhoMax [1 -3 0 0 0] 10.0; } relaxationFactors { fields { p 0.3; rho 0.01; } equations { U 0.6; e 0.6; "(k|omega|epsilon)" 0.1; } }
// ************************************************************************* // controlDict文件主要對求解過程進行設定。 /*--------------------------------*- C++ -*----------------------------------* ========= | \ / F ield | OpenFOAM: The Open Source CFD Toolbox \ / O peration | Website: https://openfoam.org \ / A nd | Version: 6 \/ M anipulation | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application rhoSimpleFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 2000; deltaT 1; writeControl timeStep; writeInterval 1000; purgeWrite 0; writeFormat binary; writePrecision 8; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable true; functions { #includeFunc MachNo #includeFunc residuals }
// ************************************************************************* // fvOptions文件中包含一些不常用的參數設定。 /*--------------------------------*- C++ -*----------------------------------* ========= | \ / F ield | OpenFOAM: The Open Source CFD Toolbox \ / O peration | Website: https://openfoam.org \ / A nd | Version: 6 \/ M anipulation | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object fvOptions; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // limitT { type limitTemperature; min 101; max 1000; selectionMode all; limitTemperatureCoeffs { selectionMode all; Tmin 101; Tmax 1000; } }
//************************************************************************** // residuals文件主要對殘差信息進行設置。 /*--------------------------------*- C++ -*----------------------------------* ========= | \ / F ield | OpenFOAM: The Open Source CFD Toolbox \ / O peration | Website: https://openfoam.org \ / A nd | Version: 6 \/ M anipulation | ------------------------------------------------------------------------------- Description For specified fields, writes out the initial residuals for the first solution of each time step; for non-scalar fields (e.g. vectors), writes the largest of the residuals for each component (e.g. x, y, z). *---------------------------------------------------------------------------*/ #includeEtc "caseDicts/postProcessing/numerical/residuals.cfg" fields (p U e k omega); // ************************************************************************* //
(4)程序執行
在cas文件所在目錄下,直接輸入rhoSimpleFoam即可。
rhoCentralFoam是OpenFOAM中一個基於密度的可壓縮求解器,適合計算跨、超聲速非定常問題。rhoCentralFoam求解器設置同樣需要準備三種文件,分別是0文件,constant文件和system文件。
直接將rhoSimpleFoam的層流計算結果拷貝過來,當作rhoCentralFoam的初始條件和邊界條件。需要拷貝的文件包括rho,p,T,U,alphat和nut六個文件。此外,還需要編寫k和omega以滿足SST湍流模型的計算。
k文件設置如下。
/*--------------------------------*- C++ -*----------------------------------* | ========= | | | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \ / O peration | Version: 4.0 | | \ / A nd | Web: www.OpenFOAM.org | | \/ M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format binary; class volScalarField; location "2000"; object k; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 -2 0 0 0 0]; internalField uniform 0.01; boundaryField { farfield { type inletOutlet; inletValue uniform 0.01; value uniform 0.01; } sym { type symmetryPlane; } wall { type kqRWallFunction; value uniform 0.01; } }
// ************************************************************************* // omega文件設置如下。 /*--------------------------------*- C++ -*----------------------------------* | ========= | | | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \ / O peration | Version: 4.0 | | \ / A nd | Web: www.OpenFOAM.org | | \/ M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format binary; class volScalarField; location "2000"; object omega; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 0 -1 0 0 0 0]; internalField uniform 10; boundaryField { farfield { type inletOutlet; inletValue uniform 10; value uniform 10; } sym { type symmetryPlane; } wall { type omegaWallFunction; Cmu 0.09; kappa 0.41; E 9.8; beta1 0.075; value uniform 10; } } // ************************************************************************* //
constant文件夾可直接複製rhoSimpleFoam案例中相關內容。僅需對turbulenceProperties文件進行適當修改,將層流模型改為SST湍流模型。
修改後的turbulenceProperties文件內容如下:
/*--------------------------------*- C++ -*----------------------------------* | ========= | | | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \ / O peration | Version: v1806 | | \ / A nd | Web: www.OpenFOAM.com | | \/ M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object turbulenceProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType RAS; RAS { RASModel kOmegaSST; turbulence on; printCoeffs on; }
// ************************************************************************* //
更換求解器後,需要對求解器的數值格式、迭代方法、求解過程等重新進行設置。具體修改如下:
由於SST湍流模型中需要得到壁面距離這一參數,fvSchemes文件增加了壁面距離wallDist的相關設置,具體內容如下。
/*--------------------------------*- C++ -*----------------------------------* | ========= | | | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \ / O peration | Version: 4.0 | | \ / A nd | Web: www.OpenFOAM.org | | \/ M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // fluxScheme Kurganov; ddtSchemes { default Euler; } gradSchemes { default Gauss linear; } divSchemes { default none; div(tauMC) Gauss linear; div(phi,k) Gauss upwind; div(phi,omega) Gauss upwind; } laplacianSchemes { default Gauss linear limited corrected 0.5; } interpolationSchemes { default linear; reconstruct(rho) Gamma 1; reconstruct(U) GammaV 1; reconstruct(T) Gamma 1; }
wallDist { method Poisson; } // ************************************************************************* //
fvSolution文件同樣增加了壁面距離yPsi的相關設置,具體內容如下。
/*--------------------------------*- C++ -*----------------------------------* | ========= | | | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \ / O peration | Version: 4.0 | | \ / A nd | Web: www.OpenFOAM.org | | \/ M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { "(rho|rhoU|rhoE)" { solver diagonal; } U { solver smoothSolver; smoother GaussSeidel; nSweeps 2; tolerance 1.0E-6; relTol 0.1; } "(e|k|omega|h)" { $U; tolerance 1.0E-4; relTol 0.1; } yPsi { solver PCG; preconditioner DIC; tolerance 1e-6; relTol 0; } } // ************************************************************************* //
controlDict文件設置如下。受格式穩定性影響,rhoCentralFoam不能太大。OpenFOAM自帶案例中一般推薦最大CFL數(maxCo)為0.2。本案例中,將maxCo設定為0.4。
/*--------------------------------*- C++ -*----------------------------------* | ========= | | | \ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \ / O peration | Version: 4.0 | | \ / A nd | Web: www.OpenFOAM.org | | \/ M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application rhoCentralFoam; startFrom latestTime; startTime 0; stopAt endTime; endTime 0.03; deltaT 0.2E-6; writeControl adjustableRunTime; writeInterval 0.0002; purgeWrite 0; writeFormat binary; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable true; adjustTimeStep yes; maxCo 0.4; maxDeltaT 1; // ************************************************************************* //
(b)
圖 4 M6機翼表面壓力分佈雲圖:(a)rhoSimpleFoam(b)rhoCentralFoam
圖 5 M6機翼表面壓力分佈計算結果(rhoSimpleFoam)與試驗結果對比
圖 6 M6機翼表面壓力分佈計算結果(rhoCentralFoam)與試驗結果對比
圖4到6分別給出了rhoSimpleFoam和rhoCentralFoam計算得到的壓力分佈曲線和雲圖。從圖中可以看出,rhoCentralFoam計算得到的壓力分佈與試驗結果相比基本吻合,而rhoSimpleFoam與試驗結果相比則存在較大差異。計算結果的差異主要在於前緣負壓區和激波附近。由於rhoSimpleFoam對數值格式極其敏感,本次計算中僅採用了一階迎風格式,從而導致壓力分佈在激波處過渡十分平滑,難以捕捉到清晰的激波結構。
(a) (b)
圖 7 M6機翼上表面激波結構:(a)rhoSimpleFoam(b)rhoCentralFoam
通過提取馬赫數等值面(Ma=1.12)的方法獲得M6機翼上表面激波結構。從圖中可以看出,rhoCentralFoam求解器獲得了清晰的λ激波結構,而rhoSimpleFoam求解器則僅捕捉到了前緣和翼尖的激波,位於機翼中部的激波則完全沒有捕捉到。
(a)
圖 8 M6機翼表面物面極限流線:(a)rhoSimpleFoam(b)rhoCentralFoam
圖8繪製了兩種求解器計算得到的機翼表面物面極限流線。結果顯示,儘管機翼上表面出現了激波/邊界層幹擾現象,但是未出現明顯的流動分離現象。
(1)rhoCentralFoam求解器是OpenFOAM中唯一的一款基於密度的求解器。跨聲速M6機翼流場計算結果顯示,該求解器獲得的壓力分佈與試驗結果基本吻合,能夠清晰地捕捉激波、局部然而,由於rhoCentralFoam採用半隱式(semi-implicit)分離(segregate)演算法,對CFL數有較大限制,實際應用中最大CFL數要求在0.2-0.4以內,導致時間推進速度十分較慢,計算效率與Runge-Kutta顯式時間推進方法接近。本算例中,最大CFL數(maxCo)為0.4,對應時間推進步長為1.4e-8秒,總推進時間約為3.4e-2秒後計算收斂,需完成約2.4e+6步非定常時間推進。
(2)rhoSimpleFoam作為一款定常求解器,能夠較為快速的獲得流場計算結果。但是由於該求解器對數值格式極其敏感,難以獲得二階及二階以上高精度收斂解。