徐振東，段宇軒，徐華松，楊 帆，李 鐵

（上海機電工程研究所，上海 201109）

0 引 言

流體數(shù)值模擬中涉及到復雜區(qū)域，需要生成貼合物面邊界的曲線網(wǎng)格，即貼體方法，通過變換坐標系來構(gòu)造貼體網(wǎng)格，由于其近壁處理的準確性已被廣泛使用。然而，為復雜結(jié)構(gòu)生成貼體網(wǎng)格可能非常耗時。模擬動邊界流動時所需的瞬態(tài)網(wǎng)格重新劃分進一步增加了計算成本和算法復雜性。浸沒邊界法能夠簡化網(wǎng)格生成。此外，當使用IB 方法時，高效的笛卡爾求解器可以直接應用于復雜問題。

Peskin開創(chuàng)的IB方法最初用于研究心血管循環(huán)中的流固耦合相互作用。到目前為止，為了提高該方法的穩(wěn)定性和適用性，眾多學者已經(jīng)提出了許多修正。在Iaccarino和Verzicco以及Mittal和Iaccarino的文章中可以找到對不同技術(shù)的詳細討論。盡管IB方法在模擬層流方面取得了非凡的成功，但其在湍流中的應用并不廣泛。大渦模擬（LES）是一種經(jīng)濟有效的湍流模擬方法，其對控制方程進行空間過濾以解決大尺度的動力學問題，并且僅對“通用”小尺度進行建模。與LES 相比，雷諾平均NS 方程（Reynoldsaveraged Navier-Stokes equations，RANS）模擬無法正確捕捉由湍流各向異性引起的二次流，直接數(shù)值模擬（direct numerical simulation，DNS）對于工程中的復雜湍流計算密集且不切實際。隨著雷諾數(shù)（）的增加，需要對近壁流進行高分辨率模擬。IB 方法的一個主要限制是笛卡爾網(wǎng)格無法將節(jié)點聚集到實體邊界，即使進行了網(wǎng)格局部細化，用于模擬湍流的笛卡爾網(wǎng)格的節(jié)點總數(shù)仍然非常大，對于三維流動這種情況尤其突出。在LES 環(huán)境中，必須明確求解近壁區(qū)域中最具能量的流動結(jié)構(gòu)，并且網(wǎng)格分辨率與DNS的分辨率相當。這使得在IB 方法中進行中高雷諾數(shù)的LES 是不可行的?？朔@一困難的一種方法是使用壁面模型求解近壁面流動，該模型通常以壁面剪應力的形式為外部LES 提供近似邊界條件。因此，在壁面?；腖ES 中，網(wǎng)格間隔尺度與黏性壁單元無關(guān)，而是與邊界層厚度有關(guān)。

Shu等提出了一種IB方法，稱為自由域離散化（DFD）方法，用于求解不規(guī)則域上的偏微分方程（partial differential equation，PDE）。在DFD 方法中，在IB 緊鄰的內(nèi)部節(jié)點處的偏微分方程的離散形式可能涉及一些外部節(jié)點，這些外部節(jié)點起到實現(xiàn)邊界條件的作用。在DFD方法的早期應用中，外部相關(guān)節(jié)點處的函數(shù)值沿僅涉及兩個邊界點的整個網(wǎng)格線進行評估，這種方式不適用于復雜的領域。為了使該方法更通用，開發(fā)了一種當?shù)谼FD 方法，外部相關(guān)節(jié)點處的數(shù)值通過沿垂直于壁的方向結(jié)合邊界條件進行適當?shù)木植客馔?，在每個時間步進行更新。

在當?shù)谼FD 方法中，應在壁法線上和解域內(nèi)構(gòu)造一些點，這些參考點通常不是網(wǎng)格節(jié)點，也被稱為虛擬點。虛擬點離壁面的距離約為一個局部網(wǎng)格間隔。到目前為止，當?shù)谼FD方法已成功應用于模擬各種無黏性或?qū)恿鳌hou通過引入壁面?；夹g(shù)將該方法擴展到高雷諾數(shù)湍流的RANS 模擬，以減輕對湍流邊界層的網(wǎng)格分辨率的要求。

本文在LES 框架中實現(xiàn)了當?shù)谼FD 方法，動態(tài)亞網(wǎng)格尺度（SGS）模型用于湍流閉合，并且將由Meneveau 等提出的拉格朗日平均法應用于SGS 模型系數(shù)的動態(tài)計算。為了減少靠近壁面區(qū)域的網(wǎng)格間隔，引入了基于湍流邊界層方程（turbulent boundary-layer equations，TBLEs）的壁面模型。

本文詳細討論了外部相關(guān)節(jié)點處的流量變量的處理，并進行數(shù)值實驗以驗證當前LES-DFD 方法的能力。

1 數(shù)值近似的控制方程和基本方案

本文考慮等密度和等黏性系數(shù)的紊流。過濾形式的連續(xù)性方程和動量方程表達為

對于本文中使用的四面體網(wǎng)格，過濾尺度與網(wǎng)格間距相關(guān)，通過單元體積的立方根計算。為了確定系數(shù)C，引入了Meveneau等提出的拉格朗日平均法。

1.1 控制方程的空間近似

使用式（3）對τ的定義，控制方程（1）和（2）可以改寫為以下無量綱形式：

式中：、f、g分別是流動變量、對流通量和黏性通量矢量；I是修正單位矩陣。

式（6）中：為無量綱層流黏性系數(shù)；=/為雷諾數(shù)，為特征速度，為特征長度。

I是修正單位矩陣，從連續(xù)性方程中消除壓強的時間導數(shù)

假設?R是一個包含一個實體的連通開集，由四面體網(wǎng)格離散化。將Mavriplis和Belov提出的Galerkin有限元方法應用于可壓縮Navier-Stokes方程和集總質(zhì)量矩陣的概念，可以得到節(jié)點處控制方程（5）的半離散形式為

式中：表示非黏性通量張量；表示黏性通量張量；n是共享頂點的所有四面體的數(shù)量；Ω是這些四面體的體積總和。如圖1所示，表示與相對的三角形面的有向（向外法線）面積。F、F和F是三個頂點處的非黏性通量張量，G是每個四面體上的恒定黏性通量。

圖1 四面體影響域邊界面的有向面積Fig.1 Directed area of influence-domain-boundaryface of a tetrahedron

1.2 時間離散化

雙時間步長用于控制方程的時間離散化。

空間離散化將式（5）變換為式（9）所示的耦合常微分方程組。

式中：是計算節(jié)點的總數(shù)；是時間步長序號；R()是離散對流通量()和黏性通量()的總和。

將關(guān)于虛擬時間的推導添加到式（9）中，得到

在每個物理時間步，式（10）的解是通過在虛擬時間內(nèi)進入穩(wěn)態(tài)而得到。

根據(jù)Belov 等的工作，在殘差R()中添加了人工可壓縮性項，以減少聲波和對流波的差異。引入當?shù)仡A處理矩陣，可產(chǎn)生

2 IB方法

根據(jù)第1.1 節(jié)中描述的空間離散化，任何網(wǎng)格節(jié)點的解都取決于通過四面體的邊與該節(jié)點相連的相鄰節(jié)點的解。對于與固體邊界相交的單元邊緣，內(nèi)部節(jié)點（位于流體域中）的流動變量通過求解控制方程來計算，邊的另一端是該內(nèi)部節(jié)點的外部相關(guān)點。如圖2所示，節(jié)點是內(nèi)部點、和的外部相關(guān)點。為了使控制方程的離散形式封閉，外部相關(guān)點（點）上的流動變量值通過邊界附近解的近似形式進行重構(gòu)。解的近似形式的構(gòu)建通過沿著物面法向的外插實現(xiàn)，插值過程中嵌入邊界條件，這是當?shù)谼FD 方法的基本思想。

為了求解外部相關(guān)點上的流動變量，需要構(gòu)造一個虛擬點。虛擬點在三個維度上的定義描述如圖2所示。對于外部相關(guān)點，首先計算法線與壁面的交點，以及法線與頂點都是內(nèi)部節(jié)點的離墻最近的三角形面的交點。和之間的距離用δ表示，然后可以得到一個常數(shù)=maxδ，其中N是外部相關(guān)點的總數(shù)。法線上的點與主體截距點的距離為，被定義為外部相關(guān)節(jié)點的虛擬點。處的流量變量可以通過其主四面體上的線性插值獲得（圖2中的）。

圖2 浸沒邊界模型Fig.2 Immersed boundary model

2.1 層流外部相關(guān)點處流動變量的計算

在層流的當?shù)谼FD 方法中，外部相關(guān)點（圖2中的）處的笛卡爾速度分量通過使用虛擬點處的確定值并施加無滑移邊界條件確定。為了考慮體加速度的影響，在計算節(jié)點處的壓強時，考慮了以下垂直于壁面方向的簡化無黏動量方程：

式中：是流體的速度矢量；是垂直于壁面的單位外向矢量。使用無穿透條件，可以得到處的壓強為

2.2 湍流外部相關(guān)節(jié)點處過濾流動變量的評估

在這項工作中，引入了一種基于簡化TBLEs的壁模型，其中僅保留擴散項，即Tessicini 等提出的平衡應力模型，用于為LES 提供近似壁邊界條件。在確定外部相關(guān)節(jié)點處的速度分量時，不采用無滑移條件，而是采用無穿透條件和壁模型產(chǎn)生的壁剪應力。本文采用平衡壁面模型。

根據(jù)壁面截點處無穿透的條件，線性外推法給出處的法向速度

如果網(wǎng)格足夠細，則可以使用無滑移邊界條件和線性外推計算處的切向速度，如圖3（a）所示。為了降低對近壁面網(wǎng)格分辨率的要求，采用壁面建模技術(shù)，并通過規(guī)定的壁面剪應力計算處的切向速度。如圖3（b）所示，壁面模型將無滑移邊界條件轉(zhuǎn)化為了滑移邊界條件（點的切向速度非零）。

圖3 浸沒邊界的壁面模型Fig.3 Wall modeling for immersed boundary

壁面剪應力τ可以近似為

式（16）中的τ是通過壁面?；夹g(shù)計算的，采用的壁面模型是基于湍流邊界層方程假設的。

在已知法向和切向分量的情況下，可以直接獲得外部相關(guān)節(jié)點處速度的笛卡爾分量。利用式（13）計算壓強。最后，在IB附近的一個節(jié)點處的離散形式是封閉的。

3 數(shù)值試驗

后臺階流動的幾何外形簡單，但是包含分離、再附著、分離泡等復雜流動現(xiàn)象，適用于LES 方法的可靠性驗證。

圖4為后臺階的幾何機構(gòu)和計算域的示意圖，計算域的流向、橫向、展向尺寸分別為：L/=40、L/=6、L/=3，擴張比為L/(L-)=1.2，臺階前方的長度為L/=40。雷諾數(shù)=/為5 000。入口條件為速度條件，速度型取Spalart獲得的平板湍流邊界層速度型；出口條件為對流邊界條件，即=L面上的Cartesian速度分量u滿足

圖4 后臺階流動的計算域示意圖Fig.4 Flow over back step

式中：為出口的平均流向速度。

計算域的上邊界采用對稱邊界條件；展向采用周期性邊界條件；下壁面上的無滑移邊界條件利用非平衡解析壁面模型轉(zhuǎn)化為無穿透邊界條件和壁面切應力條件，邊界網(wǎng)格與后臺階物面不重合。

圖5分別為采用平衡型壁面模型和無壁面模型計算得到的平均流場的流線。Le等的DNS 結(jié)果顯示在臺階拐角處存在一個二次分離泡，無壁面模型計算得到的二次分離泡的大小與DNS 結(jié)果相差較大。

圖5 平均流場的流線Fig.5 The streamlines of the mean flow field

圖6為流向位置/=0.5 和/=1 處的時間平均流向速度分布，并與DNS 結(jié)果進行對比。圖中紅色實線表示采用壁面模型的計算結(jié)果，綠色虛線表示不采用壁面模型的計算結(jié)果，黑色虛線表示DNS 的計算結(jié)果。由于不采用壁面模型得到的二次分離泡明顯偏小，因此不采用壁面模型得到的速度型在壁面附近與DNS 結(jié)果偏差較大，而采用了壁面模型得到的速度型與DNS 結(jié)果更接近。

圖6 時間平均流向速度Fig.6 Mean flow velocities

圖7為采用平衡壁面模型計算得到的雷諾應力（、、），并與無壁面模型和Jovic、Driver 的實驗結(jié)果進行對比。圖中紅色實線表示采用平衡壁面模型得到的計算結(jié)果，藍色虛線表示采用無壁面模型得到的計算結(jié)果，方形符號表示實驗結(jié)果?？梢钥闯?，平衡壁面模型的計算結(jié)果與參考的實驗結(jié)果吻合良好，無壁面模型的計算結(jié)果在近壁面與試驗結(jié)果誤差相對更大，該算例的計算結(jié)果一定程度上說明了LES-DFD 方法與平衡壁面模型結(jié)合的可靠性。

圖7 不同流向位置處的雷諾應力Fig.7 Reynolds stress at different flow directions

4 結(jié)束語

本文將當?shù)谼FD 方法推廣應用于湍流LES 模擬。采用基于簡化的湍流邊界層方程的壁面模型對浸沒邊界進行處理，將無滑移邊界條件轉(zhuǎn)化為無穿透邊界條件和壁面剪切應力，以計算外部相關(guān)節(jié)點上的速度分量。

通過數(shù)值試驗將后臺階流動的平均流向速度和雷諾應力與已發(fā)表的實驗數(shù)據(jù)和數(shù)值結(jié)果進行了比較。結(jié)果表明，將壁面模型引入DFD方法可以有效地降低近壁面網(wǎng)格分辨率的要求，大大降低了計算成本。采用壁面模型得到的計算結(jié)果比采用無壁面模型得到的計算結(jié)果更接近參考結(jié)果，驗證了引入壁面模型的可靠性和高效性。