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

（上海機(jī)電工程研究所，上海 201109）

0 引 言

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

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

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

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

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

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

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

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

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

1.1 控制方程的空間近似

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

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

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

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

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

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

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

1.2 時(shí)間離散化

雙時(shí)間步長(zhǎng)用于控制方程的時(shí)間離散化。

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

式中：是計(jì)算節(jié)點(diǎn)的總數(shù)；是時(shí)間步長(zhǎng)序號(hào)；R()是離散對(duì)流通量()和黏性通量()的總和。

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

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

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

2 IB方法

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

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

圖2 浸沒(méi)邊界模型Fig.2 Immersed boundary model

2.1 層流外部相關(guān)點(diǎn)處流動(dòng)變量的計(jì)算

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

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

2.2 湍流外部相關(guān)節(jié)點(diǎn)處過(guò)濾流動(dòng)變量的評(píng)估

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

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

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

圖3 浸沒(méi)邊界的壁面模型Fig.3 Wall modeling for immersed boundary

壁面剪應(yīng)力τ可以近似為

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

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

3 數(shù)值試驗(yàn)

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

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

圖4 后臺(tái)階流動(dòng)的計(jì)算域示意圖Fig.4 Flow over back step

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

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

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

圖5 平均流場(chǎng)的流線Fig.5 The streamlines of the mean flow field

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

圖6 時(shí)間平均流向速度Fig.6 Mean flow velocities

圖7為采用平衡壁面模型計(jì)算得到的雷諾應(yīng)力（、、），并與無(wú)壁面模型和Jovic、Driver 的實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。圖中紅色實(shí)線表示采用平衡壁面模型得到的計(jì)算結(jié)果，藍(lán)色虛線表示采用無(wú)壁面模型得到的計(jì)算結(jié)果，方形符號(hào)表示實(shí)驗(yàn)結(jié)果。可以看出，平衡壁面模型的計(jì)算結(jié)果與參考的實(shí)驗(yàn)結(jié)果吻合良好，無(wú)壁面模型的計(jì)算結(jié)果在近壁面與試驗(yàn)結(jié)果誤差相對(duì)更大，該算例的計(jì)算結(jié)果一定程度上說(shuō)明了LES-DFD 方法與平衡壁面模型結(jié)合的可靠性。

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

4 結(jié)束語(yǔ)

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

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