風雷軟件LES 開發設計與驗證

張子佩，趙鐘，，陳堅強，劉健，鄧小兵

1．空氣動力學國家重點實驗室，綿陽 621000

2．中國空氣動力研究與發展中心 計算空氣動力研究所，綿陽 621000

進入21 世紀以來，計算流體力學（Computational Fluid Dynamic， CFD）的快速發展和廣泛應用極大地改變了航空航天飛行器設計過程，不僅能有效減少物理試驗次數，還可模擬一些物理試驗無法開展的工況，已成為飛行器設計過程中必不可少的方法手段。

當前，計算流體力學的研究方法正從以定常、小范圍分離或附著流動為特征的RANS（Reynold-Averaged Numerical Simulation）湍流模擬，向非定常、大范圍分離和流動摻混為特征的非線性多尺度問題轉移。RANS 方法能夠很好地預測定常小范圍分離或附著流動，但對大范圍分離流動的模擬能力不足。大渦模擬（Large Eddy Simulation， LES）[1］和直接數值模擬（Direct Numerical Simulation， DNS）[2］方法被認為是解決這類問題的有效手段。DNS 方法對流動中的所有尺度進行直接求解，模擬高雷諾數流動需要的網格量達到了Re37/14L量級[3］。LES 方法只對流動中的大尺度脈動直接求解，對小尺度脈動采用模化方法，分為壁面解析的大渦模擬（Wall-Resolved LES， WRLES）和壁面模化的大渦模擬（Wall-Modelled， WMLES）。對于低雷諾數流動，兩者所需的網格量差異不大，但對于高雷諾數流動，前者需要的網格量為Re13/7L，后者需要的網格量為ReL[3］。按當前計算機的發展水平，即使對106量級雷諾數下的小展弦比機翼模擬，預計也要到2030年，WMLES 才可以在24 h 內完成[4］。

LES 的基本思想是選擇處于慣性子區內的濾波尺度，將湍流脈動量分為可以數值求解的大尺度脈動量和需要亞格子尺度（Subgrid Scale，SGS）模型模化的小尺度脈動量，由于小于濾波尺度的脈動量接近于各向同性，所以亞格子模型比RANS 模型更為容易和普適。對于大尺度脈動量所主導的流動問題，例如氣動噪聲[5］、湍流燃燒[6］、激波/邊界層干擾[7］、大范圍分離流動[8］等領域，LES 方法正在發揮著越來越重要的作用。

為了推動LES 方法盡快走向成熟應用，學者正從多個方面開展研究，主要包括：空間離散格式、時間推進格式、亞格子模型、入流生成方法、壁面模化方法、非結構網格應用、軟件工程化等。

空間離散格式方面，迎風格式的適用性是焦點。一般認為RANS 方法常用的二階迎風格式耗散性大，會使亞格子模型的截斷波數難以確定，尤其是對復雜構型而言，由于非均勻湍流的存在，截斷波數更加難以確定，因而對高階迎風格式提出更迫切的需求[9］。

時間推進格式方面，較低的時間精度會使高精度的空間離散方法無法達到效果，Choi 和Moin[10］認為時間步長需滿足u2τΔt/ν ＜1，其中uτ為壁面摩擦速度，Δt 為時間步長，ν 為黏性系數。

亞格子模型方面，Nicoud 等[11］認為渦黏性模型需要滿足4 個特性：①只利用流場的當地量，并且要保正；②滿足近壁漸進特性；③在剛性轉動、純剪切或者二維流動中，模型退化為0；④在單純軸對稱或者單純各向同性膨脹/收縮流動中，模型退化為0。

壁面模化方法方面，對于Re＞106的邊界層流動，壁面模化的大渦模擬方法要將網格量控制在O（Re），對工程應用才有實際意義。Piomelli[12］將 壁 面 模 化 方 法 劃 分 為：①壁 面 模 型，例如Schumann[13］提 出 的 平 衡 層 模 型；②在 壁 面 附近單獨求解邊界層方程的分區方法[14］；③RANS/LES 混合方法，例如DES 類方法[15］。

入流生成方法方面，在壁面附著流動中，層流邊界層和轉捩區的計算成本甚至是湍流區的10～100 倍，如果只關注湍流區域，需要準確刻畫湍流入流。目前采用的入流生成方法有“回收”方法[16-17］、人工合成湍流方法[18］等。

在非結構網格應用方面，由于生成非結構網格相對省時省力，基于非結構網格的LES 方法開始在復雜工程問題中獲得越來越多的應用，例如Forsythe 等[19］利 用 非 結 構 求 解 器 實 現 了F-15E全機外形的DES 模擬。但是傳統非結構求解器的空間離散精度只有二階，因而基于間斷伽遼金（Discontinuous Galerkin， DG）、有 限 譜 體 積（Spectral Volume， SV）等高精度求解的LES 方法是目前的研究熱點[20］。

在工程應用方面，雖然一些商業軟件中已有LES 計算模塊，但是實際應用還不成熟，還處于從“算法”向“應用”的過渡階段。由于工業設計部門缺乏對湍流非定常和多尺度效應的正確認識，在嘗試應用LES 時，誤認為僅將RANS 方法的湍流模式改為LES 方法的亞格子模型，就可以達到可觀的效果。實際上，LES 緊密依賴于所采用的數值格式、網格、初邊值條件等因素，這些因素間要互相匹配才能開展LES 模擬，不能簡單地“一鍵化”操作。

綜上所述，由于存在種種尚待解決的問題，LES 方法還難以像RANS 方法那樣成熟地開展批量工程應用，更大程度上還處于方法基礎研究、應用基礎研究階段。在商業軟件以外，大量基礎研究者更多地是采用自己編寫的課題組代碼、“in-house”代碼開展研究。

在國內，由于此前長期缺乏共享的自主CFD軟件框架，眾多研究者在開展LES 方法研究中，將大量精力耗費在開發前后處理接口、計算格式、并行計算等重復性工作上，力量難以向計算模型、流動機理等更為根本的方向聚焦。

風 雷 軟 件（PHengLEI）[21-22］，是 中 國 空 氣 動力研究與發展中心研發的、具有自主知識產權的CFD 軟件，歷經十余年發展，已成為唯一一款能同時支持結構、非結構網格的大型CFD 軟件平臺。在國家數值風洞工程（NNW）的支持下，風雷軟件致力于面向中國研究者提供高可擴展性的CFD 軟件基礎設施，共同推進中國CFD 軟件生態建設，于2020 年12 月面向全國開源，在互聯網“紅山開源平臺”開放代碼（平臺網址https：∥osredm.com）[23］。

LES 模型是風雷軟件開源框架中的重要組成部分，旨在于為國內開展LES 基礎和應用研究的學者，提供高起點工具，已在PHengLEI v2112版本中開源。主要功能包括：基于結構/非結構二階精度有限體積方法的LES 模型，基于結構高階精度有限差分方法的LES 模型，基于Fourier譜/有限差分混合方法的LES 模型，以及LES 模擬中常用的亞格子模型、初邊值條件、非定常湍流統計等方法。

本文主要介紹風雷軟件開源框架LES 模型的理論方法、軟件設計和算例驗證。首先給出LES 的理論方法，然后介紹設計的軟件框架，最后通過不可壓縮槽道湍流、亞臨界雷諾數圓柱、NACA0012 翼型臨界攻角的低頻振蕩等標準算例，給出軟件的驗證結果。

1 理論方法概述

不可壓縮流動和可壓縮流動的求解方法不同，采用投影法求解不可壓縮流動，空間離散采用Fourier 譜與高階有限差分相結合的混合方法；采用二階有限體積和高階有限差分方法求解可壓縮流動。下面對這兩類流動求解方法分別作簡要介紹。

1.1 基于Fourier 譜/有限差分的投影法

投影法（Projection Method）也稱作分步法（Fractional Step Method），由Chorin[24-25］在20 世紀60 年代末提出，最早認為在時間方向上只有一階精度，后來許多學者發展了各種改進的投影法，將時間方向提高為二階精度[26-27］。

1.1.1 投影法

無量綱后的非定常不可壓縮流動的控制方程為

對流項采用二階精度的Adams-Bashforth 格式，時間導數項和黏性項采用Crank-Nicholson格式離散[26］，控制方程離散為

2）投影步。根據Helmholtz-Hodge 分解定理，u*可分解為

也即? 要滿足泊松方程

求出? 后，壓力由式（8）求得

投影步結束后式（9）成立

3）壁面散度修正步。由式（7）求得的? 只能保證u**在壁面的法向分量為0，切向分量不能保證為0。為此，在利用式（5）求速度u**時加入無滑移邊界條件

這樣由式（10）求得的速度u**在壁面上散度不為0。為此，采用影響矩陣法對壁面散度進行修正，具體不再詳述，可參考文獻[28］。

1.1.2 Fourier 譜/有限差分混合方法

從以上的分析可看出，本文采用投影法離散控制方程后，需要求解兩種不同邊界條件類型的泊松方程，第1 種是包含Dirichlet 邊界條件的泊松方程

對于槽道湍流，流向和展向均為周期邊界，可以采用Fourier 級數展開，法向采用差分方法，用x、y、z 分別表示流向、展向、法向，則速度和壓力分別為

式中：x=(x，y，z)T；kx、ky分別為流向波數和展向波數；k=(kx，ky)T。利用Fourier 基 函數的正交性，算子?2可以離散為

這樣，對kx、ky，式（11）、式（12）就轉換為z 方向的一維二階微分方程。z 方向的離散采用Gamet 等[29］發展的非均勻網格空間導數求解方法，得到空間二階導數的一般形式為

非線性項(u·?)u 的計算采用“偽譜”方法，并利用3/2 規則來消除混淆誤差，具體可參考文獻[30］。

式（12）的求解也采用上述方法。

1.2 二階有限體積和高階有限差分方法

對三維可壓縮NS 方程采用密度加權的Favre 濾波，得到LES 控制方程為

式 中：δij為Kronecker 符 號；μl為 層 流 黏 性 系 數；Prl為層流普朗特數；Cp為定壓比熱容。式（17）中上標“SGS”表示亞格子項為黏性應力的非線性項為亞格子黏性性擴散項為熱通量的非線性項，這3 項通常都被忽略[31］。一般關注的是亞格子應力張量和亞格子熱通量項，通過采用亞格子模型使之封閉，亞格子湍流擴散項進行線性化處理[31］

二階求解器采用結構或非結構的有限體積方法[21］，無黏項采用Roe 格式計算，黏性項采用二階中心型格式。高階求解器采用基于結構網格的高精度有限差分方法[32］，無黏項采用WCNS格式[33］，黏性項采用四階中心型格式。時間推進均采用LU-SGS 隱式方法求解，非定常計算均采用雙時間步法[34］。

1.3 亞格子模型

主要采用Smagorinsky 模型、動態Smagorinsky 模型、Sigma 模型來模化亞格子項。

1）Smagorinsky 模型。Fourier 譜/有限差分方法LES 求解器采用原始不可壓縮形式的Smagorinsky 模型（Smagorinsky Model， SM）[35］

有限體積/有限差分LES 求解器采用Yoshizawa[36］提出的形式

式中：Prt為湍流普朗特數

2）動 態Smagorinsky 模 型。Fourier 譜/有限差分方法LES 求解器采用原始不可壓縮形式的動態Smagorinsky 模型（Dynamic Smagorinsky Model， DSM）[37］，模型系數在水平面平均以提高數值穩定性

有限體積/有限差分LES 求解器采用Lu等[39］提出的形式

3）Sigma 模型。Nicoud 等[11］提出的Sigma模型滿足O(y3)壁面漸進特性，在剛性轉動、純剪切、單純軸對稱、單純各向同性膨脹/收縮流動或者二維流動中可以退化為0。Fourier 譜/有限差分方法LES 求解器采用的模型方程為

式中：Cσ為模型系數，根據各向同性衰減湍流算例的標定，取為1.5；σ1～σ3為速度梯度張量gij的奇異值，并且σ1≥σ2≥σ3，σ1～σ3的具體求法可參考文獻[11］。

2 軟件架構設計

風雷軟件設計了具有良好通用性、可擴展性的體系結構和數據結構，在同一個軟件平臺上，能同時支持結構和非結構網格（算法）、二階和高階精度（算法）。風雷軟件LES 模塊提供了一個開放平臺，可使國內的LES 研究人員從前后處理接口、通量格式、時間離散格式等重復開發中解放出來，將精力聚焦于學科前沿的創新研究，為更多有志于進入這一領域的學者和團隊提供高水平的起點，為需要應用LES 方法開展基礎流動問題和工程應用問題的學者和團隊提供高效高精度的開放平臺。

風雷軟件LES 模擬同時包含松/緊耦合策略。緊耦合策略，是在二階結構求解器、二階非結構求解器、高階結構求解器上集成了LES 模型，統稱為“有限體積/有限差分方法LES 求解器”；松耦合策略，是結合Fourier 譜方法和有限差分方法，開發的求解不可壓縮槽道流動的功能模塊，稱為“Fourier 譜/有限差分方法LES 求解器”。

風雷軟件整體框架可參考文獻[21-22］，這里不做詳細介紹，僅詳細介紹與LES 方法相關的軟件框架設計。

2.1 LES 計算框架設計

對于松耦合策略，設計的Fourier 譜/有限差分方法LES 求解器計算框架如圖1 所示。計算框架由前處理、迭代流程、后處理3 個模塊構成。前處理包括網格/參數讀入、初始化物理空間/波數空間；迭代流程中每個迭代步都要遍歷所有波數；后處理包括流場變量在波數空間的統計輸出、流場變量在物理空間的統計輸出等。

圖1 Fourier 譜/有限差分LES 求解器的計算框架Fig.1 Design diagram of Fourier spectrum/finite difference solver computing framework

Fourier 譜/有限差分方法LES 求解器的功能模塊分層設計如圖2 所示，也按照內核層、算法層、功能層和應用層進行分層構建。在內核層，復用了風雷軟件的計算參數讀取，在原有面向實數的數據結構基礎上，進一步設計了復數數據結構，采用P3DFFT[40］進行并行Fourier 變換；在算法層，抽象出物理空間計算域和波數空間計算域，設計了高階緊致有限差分格式和濾波函數；功能層包括微分算子、求解控制流程、前后處理、工具類等。微分算子方面，構建了拉普拉斯算子、梯度算子、散度算子等微分算子模塊。前后處理方面，構建了流場變量在波數空間和物理空間之間的轉換功能模塊，以實現流場數據的可視化。求解控制流程方面，設計了對流步、投影步和壁面散度修正的投影法求解流程。工具類方面，構建了添加速度脈動的模塊等；在應用層，利用功能層的拉普拉斯算子、梯度算子和散度算子等，構建了泊松方程求解器。湍流求解不再作為單獨的求解器，而是封裝為亞格子模型。泊松方程求解器和亞格子模型結合，拓展為Fourier 譜/有限差分方法LES求解器。

圖2 Fourier 譜/有限差分方法LES 求解器的功能模塊分層設計示意圖Fig.2 Functional module hierarchical design of Fourier spectrum/finite difference solver

對于緊耦合策略，設計的有限體積/有限差分方法LES 求解器計算框架如圖3 所示。計算框架由前處理、迭代流程、后處理3 個模塊構成。前處理包括網格/參數讀入、網格轉換/分區、求解器加載等；迭代流程由“迭代控制器（Controller）”控制。

圖3 有限體積/有限差分LES 求解器的計算框架Fig.3 Design diagram of finite volume/finite difference solver computing framework

首先根據Zone 的信息初始化每個Zone 上的所有求解器。然后開始迭代計算，在每個迭代步內包括2 層循環遍歷：第1 層遍歷不同求解器，第2 層是在每個求解器中遍歷每個Zone；后處理包括湍流一階/二階矩的統計平均、Fourier 分析、模態分析等工具包。

對于有限體積/有限差分LES 計算，每個Zone 上加載的求解器包括NS 方程求解器NSSolver 和有限體積/有限差分LES 求解器LESSolver。對原用于進行RANS 流動模擬的NS 方程求解器NSSolver 作細微改造：將湍流黏性通量修改為反對稱項和對稱項2 部分，與SGS 應力項的對稱項和反對稱項相對應。

LESSolver 的功能模塊分層設計如圖4 所示，與風雷軟件的RANS 湍流模式方程求解器TurbSolver 類似，按照內核層、算法層、功能層、應用層進行分層構建。由于軟件現有的亞格子模型均為代數形，依賴當地的網格尺度、速度梯度張量等信息，不需要求解模型方程，因而在功能層的求解控制流程模塊刪除黏性項、對流項等求解過程，在算法層增加濾波函數、壁面函數等模塊。風雷軟件原有數據通信模塊采用異步通信模式，可以減少通信次數，但會影響數值穩定性，LESSolver 在求解SGS 應力項前增加一次速度梯度張量的計算和數據通信，以提高數值穩定性。后處理模塊增加湍流統計輸出功能。

圖4 有限體積/有限差分LES 求解器的功能模塊分層Fig.4 Functional module hierarchical design of finite volume/finite difference solver

2.2 LES 求解器設計

風雷軟件將幾何計算域Zone 與求解器Solver 作為數值計算的基礎。Zone 是框架的核心，負責存儲網格、加載求解器、存儲流場變量等，每個Zone 上可以有多個網格塊Grid，每個Grid 既可以是結構網格也可以是非結構網格。根據計算任務，每個Zone 上可以分別加載不同的Solver[22］。

將不同Solver 的共性提取出來抽象為求解器基類，從基類派生出具有差異化的求解器子類。如圖5 所示，由基類派生出的CFD 求解器類CFDSolver，其中定義了CFD 求解器的共性成員和接口，如計算網格、空間離散、時間推進、邊界條件處理等純虛函數接口。

SpecDiffHybSolver 類 是Fourier 譜/有 限 差分求解器，由于所用的空間離散和時間推進格式與傳統基于有限體積方法的CFD 求解有很大不同，因此直接從Solver 基類派生。

原始CFDSolver 派生得到NS 方程求解器類NSSolver、湍流模式方程求解器類TurbSolver。由于LES 所采用的亞格子模型多數為代數形，不需要求解模式方程，因此直接從CFDSolver 派生出有限體積/有限差分LES 求解器類LESSolver。

LESSolver 又可以派生出基于不同網格類型的求解器類：基于結構網格的LESSolverStruct和基于非結構網格的LESSolverUnstruct。NSSolverStruct 又派生出結構網格高精度求解器類NSSolverStructFD，而基于NSSolverStructFD的大渦模擬求解器類被融合在LESSolver-Struct 中。

LES 數值模擬流場數據可供氣動聲學、人工智能等領域專家使用，對此，采用CGNS（CFD General Notation System）數據底層所采用的HDF5（Hierarchical Data Format）標 準 數 據 格式，以保證數據的可擴展性。相關格式見軟件用戶手冊[41］。

2.3 并行計算模式設計

Fourier 譜/有限差分方法LES 求解器是對泊松方程式（12）、式（13）在譜空間下進行求解，需要對所有流向波數kx和展向波數ky的遍歷，任意波數對（kx，ky）的求解與其他波數無關，因而理論上所有波數對的求解是完全獨立的。然而對流項在“偽譜”方法中是通過變換到物理空間求解的，這個過程的并行計算需要并行快速Fourier變換（Fast Fourier Transforms， FFT）。因此，高效的并行快速Fourier 變換對整個求解器的并行效率和可擴展性有決定性的影響。

Fourier 譜/有限差分方法LES 求解器需要對數據在2 個維度各進行一次FFT，傳統并行方法是把數據只在其中一個維度分布到n個進程。以4 進程并行計算三維槽道為例，如圖6（a）所示，需要在x方向和z方向各進行一次FFT。首先，將數據在z方向分為4 份，分配到進程1～4，使得每個進程中的數據在x方向是連通的，這樣每個進程就可以在x方向進行FFT；然后，將上一步變換后的數據在x方向分為4 份，重新分配到進程1～4，使得每個進程中的數據在z方向是貫通的，這樣每個進程就可以在z方向進行FFT。這個方案僅能在一個方向上劃分計算域，稱為一維分解。該方案計算網格塊不能進一步細分，限制了其并行粒度，并行規模也難以增加。

圖6 并行快速Fourier 變換：一維分解和“束”分解Fig.6 1D and pencils decomposition of parallel FFT computations

近年來，加州大學開發了一款開源工具P3DFFT[40］專 門 處 理 并 行Fourier 變 換，可 將 需要進行FFT 的數據在其中一個維度分布到n1組，同時將不需要變換的維度分布到n2組。以16 進程并行計算三維槽道為例，如圖6（b）所示，需要在x方向和z方向各進行一次FFT。首先，將數據在y方向和z方向同時分為4 份，分配到進程1～16，使得每個進程中的數據在x方向是貫通的，這樣每個進程就可以在x方向進行FFT；然后，將上一步變換后的數據在x方向和y方向各分為4 份，重新分配到進程1～16，使得每個進程中的數據在z方向是貫通的，這樣每個進程就可以在z方向進行FFT。這個方案可以2 個方向上劃分計算域，稱為二維分解或者“束”分解（Pencils Decomposition）。本文采用“束”分解。

并行策略方面，有限體積/有限差分方法LES 求解器的并行計算流程與風雷軟件基本一致：當一種求解器類型的計算執行完成后，每個分區之間執行一次并行通信，同時在分區內部網格間執行一次邊界面數據交換。當求解器列表遍歷完成后，也要進行求解器之間的數據通信，以實現不同求解條件下流場或特殊邊界的數據傳遞和交換[34］。

通信模式方面，采用了基于循環遍歷所有網格塊的模式進行通信：每個進程均遍歷所有網格塊（僅僅是塊編號），當所遍歷的網格塊屬于當前進程時，則向鄰居網格塊所在進程發送數據；當所遍歷的網格塊屬于其鄰居網格塊時，則從鄰居網格塊所在進程接收數據，詳細參考文獻[21］。

3 算例驗證

采用不可壓縮槽道湍流、亞臨界雷諾數圓柱繞流、NACA0012 臨界攻角的低頻振蕩等3 個算例，對本文設計的Fourier 譜/有限差分方法LES求解器、有限體積/有限差分LES 求解器進行驗證。

3.1 不可壓縮槽道湍流

不可壓縮槽道湍流由于外形簡單且流動現象單一，是研究壁面剪切流動的標準算例，也是考量數值方法和亞格子模型的基礎算例，使用本算例可以考核本文設計的Fourier 譜/有限差分方法LES 求解器。譜方法是數值誤差最小的一種離散方法，因此是考察亞格子模型誤差的最理想方法。

圖7 是不可壓縮槽道湍流算例，計算域大小為Lx×Ly×Lz=2πh×2h×πh（h是槽道半高度），x和z方向是周期邊界條件，y方向是無滑移壁 面，雷 諾 數Reτ=uτh/ν=395（uτ為 壁 面 摩 擦速度）。

圖7 不可壓縮槽道湍流的示意圖Fig.7 Incompressible turbulent channel flow

對該算例，通常做法是在流向（x方向）、展向（z方向）采用Fourier 譜方法，在法向（y方向）采用基于Chebyshev 多項式的譜方法[42］。而本文在法向采用高階有限差分方法，這是因為，若在法向采用基于Chebyshev 多項式的譜方法，則要求網格點分布必須滿足Chebyshev 多項式關系，這極大地限制了在高雷諾數流動中的應用，而本文在法向采用高階有限差分方法，可打破這一限制，比傳統方法更具優勢。

算例A～算例D 采用不同的亞格子模型和網格 尺 度，計 算 結 果 與Moser 等[43］的DNS 結 果、Nicoud 等[11］的LES 結果對比。算例A～算例D和文獻中采用的數值方法、計算域、網格、亞格子模型如表1 所示，Nx、Ny、Nz為3 個方向的網格數目，Δx+、Δy+、Δz+為無量綱化后的網格尺度。算例A 與Moser 等[43］的DNS 算例采用完全相同的計算域和網格尺度，只是數值方法不同，用于考核Fourier 譜/有限差分方法LES 求解器的數值精度。法向網格為雙曲正切分布，流向和展向網格均勻分布，Moser 等[43］在流向和展向離散采用與算例A 相同的Fourier 譜方法，法向采用Chebyshev 多項式。算例B～算例D 采用了不同的亞格子模型，Smagorinsky 模型系數CS=0.1，引入van Driest 函數修正壁面附近的渦黏系數，Sigma模型系數Cσ=1.5。Nicoud 在所有方向均采用四階中心Taylor-Galerlin 有限元格式。算例B～算例D 流向和法向網格尺度比Nicoud 算例更小，展向網格尺度相當，足以滿足LES 模擬需求[3］，用于驗證Fourier 譜/有限差分方法LES 求解器具備考核亞格子模型的能力。

表1 本文算例和文獻采用的數值方法、計算域、網格和亞格子模型Table 1 Numerical methods， computational domain， grid and subgrid models used in this paper and references

將計算結果作時間平均和空間平均，得到平均流向速度剖面（見圖8）、平均雷諾應力分布（見圖9），橫縱坐標都采用壁面坐標。可以看到，算例A 的平均流向速度剖面、雷諾應力分布都與Moser 等[43］的基本一致，說明本文設計的Fourier譜/有 限 差 分 方 法LES 求 解 器 與Moser 等[43］的DNS 求解器數值精度相當，法向采用高階有限差分方法并沒有降低數值精度。Chebychev-tau 離散方法受限于Chebychev 配置點的布置，而高階有限差分方法對于網格點的布置沒有限制，這是本文設計的Fourier 譜/有限差分方法LES 求解器的一大優勢。

圖8 平均流向速度剖面Fig.8 Mean streamwise velocity profile

圖9 平均雷諾應力分布Fig.9 Mean Reynolds stress profile

算例B～算例D 的平均流向速度剖面與Moser 等[43］的DNS 結果基本一致。對算例C，與Nicoud 等[11］利 用DSM 模 型 得 到 的 速 度 剖 面 相比，更 接 近Moser 等[43］的DNS 結 果，這 是 由 于Nicoud 等[11］的DSM 動態模型系數是通過局部平均得到的，而算例C 的動態模型系數是在水平方向平均得到的。

算例B 的流向雷諾應力峰值略小于DNS 結果，在其他位置與DNS 結果基本一致，這可能是由于Smagorinsky 模型不滿足O(y3)的近壁特性。算例C 的流向雷諾應力峰值略小于DNS 結果，是由于LES 方法只解析部分雷諾應力，剩余部分雷諾應力由亞格子模型模化；而Nicoud 等[11］利用DSM 模型得到的流向雷諾應力峰值略大于DNS 結果，可能是由于數值誤差和模型誤差相互抵消導致的，其使用的有限差分方法數值誤差可能偏大，網格也可能略稀疏。算例B～算例D 的法向、展向雷諾應力與Moser 的DNS 結果基本一致，并且比Nicoud 的LES 結果更接近Moser 的DNS 結果，進一步說明了Nicoud 的數值誤差過大。

通過以上分析，本文設計的Fourier 譜/有限差分方法LES 求解器具有與Fourier-Chebychev方法相當的數值精度，在驗證亞格子模型時可以排除數值誤差的影響。

3.2 亞臨界雷諾數圓柱繞流

Re=3 900 條件下的圓柱繞流處于亞臨界狀態，分離前的邊界層為層流，分離后的剪切層逐漸失穩，尾跡呈現三維湍流特性，是用于考核鈍體繞流模擬能力的經典算例。采用不同精度、不同網格類型的LES 求解器和不同的亞格子模型，對該算例進行考察，主要為驗證不同類型的LES求解器和亞格子模型的可靠性。

計算域和網格如圖10 所示，x、y分別為流向、法向，r、θ分別為徑向、周向，展向計算域大小為πD（D為圓柱直徑）。圓柱表面是無滑移邊界，展向是周期邊界，其他均是遠場邊界。采用“O”型網格，網格分布為Nr×Nθ×Nz=434×264×48，周向有128 個點集中分布在尾跡區中心線兩側90°范圍內，展向網格均勻分布。

圖10 亞臨界雷諾數圓柱繞流計算域和網格Fig.10 Computation domain and its grid of flow past circular cylinder at sub-critical Reynolds number

設計了E～算例H 4 個算例（見表2），分別采用了不同網格類型、數值精度的求解器和不同的亞格子模型，以全面考核本文設計的有限體積/有限差分方法LES 求解器。

表2 算例采用的求解器和亞格子模型Table 2 Solvers and sub-grid scale models used in simulations

圖11 給出了中心線（y/D=0）平均速度（時間平均和展向平均，用來流速度Uc無量綱化，下同）沿x方向分布，統計平均周期T≈120Tvs（Tvs為渦脫落周期）與文獻[44-48］的對比。可以看到，后駐點處的平均速度為0，之后減小到某個極小值，然后逐漸恢復到接近來流速度，中間改變過一次速度方向。定義后駐點與速度改變方向處之間的距離為回流區長度（Lr），與各文獻中的回流區長度對比（見表3），可以看到回流區長度范圍為1.19D～1.67D。算例E、H 的回流區長度基本在文獻結果范圍之內，算例F、G 的回流區長度略大于文獻結果。

表3 回流區長度對比Table 3 Recirculation lengths comparision

圖11 圓柱尾跡區中心線上的平均流向速度分布Fig.11 Mean streamwise velocity in wake centerline

從算例E、F 的Q等值面圖（見圖12，采用x方向渦量ωx著色）可以看出，層流邊界層分離形成二維層流剪切層，隨著剪切層失穩，三維渦結構開始出現，流動逐漸過渡到湍流。相較于算例E，算例F 的剪切層在較長的尾跡區內保持二維層流狀態，這也是后文中平均速度、雷諾應力分布差異的主要原因。圖12（a）標注了近尾跡區的3 個流向站位分別為x/D=1.06，1.54，2.02，以下重點分析這3 個站位的平均速度、雷諾應力分布。

從近尾跡區3 個站位的平均流向速度分布（見圖13）看，本文計算結果、文獻[45-49］結果都是由U 型速度分布逐漸過渡到V 型速度分布，不同的是過渡點的位置。U 型速度分布對應于層流剪切層，V 型速度分布對應于湍流剪切層。在x/D=1.06 站位，算例E～算例H 都是U 型分布，與文獻結果一致。在x/D=1.54 站位，算例E、H 已過渡到V 型分布；而算例F、G 尚未完全過渡到V型分布，說明處于層流剪切層逐漸失穩向湍流剪切層過渡的狀態，這也可以從圖12 中看出。

圖12 圓柱尾跡區的瞬時Q 等值面圖Fig.12 Instantaneous iso-surfaces of Qcriterion in wake

圖13 圓柱近尾跡區平均流向速度分布Fig.13 Mean streamwise velocity in near wake

從近尾跡區3 個站位的平均法向速度分布（見圖14）看，算例E～算例H 的速度型都有很好的反對稱性。在x/D=1.06，2.02 站位，算例E、H 與文獻[45-49］中的實驗、計算結果符合得很好，算例F、G 的結果與文獻結果相比，峰值位置和峰值大小略有差異。在x/D=1.54 站位，文獻結果之間有很大的差異，算例E、H 基本處于不同文獻結果的中間范圍內，之所以有差異，是由于層流剪切層在此處失穩，逐漸過渡到湍流剪切層，是一個中間過渡狀態，數值噪聲和實驗的來流噪聲都會對結果有很大影響。

圖14 圓柱近尾跡區平均法向速度分布Fig.14 Mean normal velocity in near wake

圖15 是近尾跡區3 個站位的流向雷諾正應力u'u' 分布。在x/D=1.06 站位，u'u' 的峰值對應剪切層失穩的位置。在x/D=1.54 站位，主分離渦的出現使得上述峰值附近出現另2 個峰值。隨著湍流混合加強，2 個峰值逐漸融合為一個峰值。不同文獻之間的雷諾應力大小有很大的差異，本文的計算結果基本在文獻[45-49］結果的中間范圍內。在x/D=1.06，1.54 站位，算例F、G 的結果要小于算例E、H，這也是由于算例F、G 模擬的剪切層失穩較晚所致。

圖15 圓柱近尾跡區3 個站位的流向雷諾正應力分布Fig.15 Mean streamwise Reynolds stress at three locations in near wake

從以上分析可知，算例E、H 的結果基本一致，算例F、G 的結果基本一致，說明亞格子模型對結果的影響很小，數值方法對結果的影響較大。另外，對剪切層失穩的模擬對于速度、雷諾應力等的預測非常關鍵。剪切層失穩受到數值噪聲和實驗來流噪聲的影響，這使得不同的計算程序和實驗條件得到的結果相差很大，從不同文獻之間的差異可以看出這一點。對比本文的4 個算例，算例F、G 模擬的剪切層失穩位置在E、H 的上游，發展出湍流剪切層的位置更靠上游，速度型分布與文獻結果也更為接近，這可能是由于二階精度方法的數值誤差與亞格子模型誤差共同作用導致的。

3.3 NACA0012 臨界攻角的低頻振蕩

臨界攻角下的翼型繞流伴隨著流動分離、剪切層失穩、轉捩、流動再附等復雜的非定常流動現象，是體現高精度LES 方法優勢的一類流動問題。Zaman 等[50］在實驗中發現臨界攻角狀態下的流動除了渦脫落模態，還存在著一個更低頻的流動 模 態。Almutairi 和Alqadi[8］利 用LES 方 法 對NACA0012 翼型的低頻振蕩現象進行了數值模擬，認為后緣渦脫落引起的聲波激發了前緣剪切層，使剪切層失穩并轉捩，湍流剪切層再附又抑制了后緣渦脫落。這樣，后緣分離渦與前緣剪切層共同構成的反饋機制，是低頻振蕩現象的原因。Eljack 和Soria[51］的LES 計 算 結 果 發 現 前 緣 附 近存在的“三渦”結構導致了流動處于周期性的分離和附著狀態，使得流動存在著低頻振蕩現象。

選取攻角α=9.29°的NACA0012 翼型作為研究對象，馬赫數Ma=0.4，基于翼型弦長c的雷諾數Rec=50 000，計算域Lx×Ly×Lz=21c×20c×0.5c。網格如圖16 所示，壁面法向第1 層網格Δy+w＜0.35，流動分離、轉捩區的流向網格間距Δx+＜30，展向均勻布置96 個網格點，展向網格間距Δz+＜18，網格總量約為2 600 萬，滿足LES 方法對網格尺度的要求[3］。無量綱時間步長為ΔtU/c=0.002，采用400 核，模擬了20 萬個無量綱時間步長。

圖16 NACA0012 計算網格Fig.16 Computation grid of NACA0012

圖17 是模擬得到的瞬時Q等值面圖（采用x方向渦量ωx著色），圖中可以看到層流剪切層失穩、轉捩和三維渦結構產生的過程。

圖17 NACA0012 繞流瞬時Q 等值面圖Fig.17 Instantaneous iso-surfaces of Q criterion around NACA0012

圖18 是升力系數（CL）隨時間的變化曲線，圖中可以看出升力系數的振動幅值與平均值相當，流動狀態在附著流動與失速間交替。對升力系數的時間序列利用現代譜估計方法得到功率譜密度（Power Spectrum Density，PSD）[52］，如圖19所示，計算得到的升力系數功率譜的第一峰值對應的斯特勞哈爾數St=fcsinα/U=0.005 5，與Almutairi 和Alqadi[8］得 到 的 峰 值 位 置St=0.005基本一致，峰值大小也基本一致。

圖18 升力系數的時間變化曲線Fig.18 Time history and power spectrum of lift coefficient

圖19 升力系數的功率譜密度Fig.19 Power spectrum density of lift coefficient

4 結 論

基于風雷軟件開源框架，分別采用松耦合和緊耦合的模式，開發了Fourier 譜/有限差分方法LES 求解器和有限體積/有限差分方法LES 求解器。初步實現了基于結構網格和非結構網格，基于二階有限體積方法、高階有限差分方法和Fourier 譜方法的LES 功能。集成了多種亞格子模型、初邊值條件和非定常湍流統計功能。

分別采用不可壓縮槽道湍流、亞臨界雷諾數圓柱繞流、NACA0012 臨界攻角的低頻振蕩3 個算例，對所設計的開源求解器進行驗證。計算表明，Fourier 譜/有限差分方法的LES 求解器對于不可壓縮槽道湍流的模擬結果與DNS 結果一致，數值精度與Fourier-Chebychev 方法相當，在考核亞格子模型時可以排除數值誤差的影響。有限體積/有限差分方法的LES 求解器對于亞臨界雷諾數圓柱繞流算例的模擬結果與文獻中的計算和實驗結果基本一致，證明了有限體積/有限差分方法的LES 求解器和幾種亞格子模型的可靠性。高階精度有限差分方法的LES 求解器復現了NACA0012 臨界攻角的低頻振蕩現象，St數值與文獻中基本一致，體現了高階精度有限差分方法的LES 求解器對于復雜湍流模擬的能力。

以上代碼、接口和算例均已在風雷軟件PHengLEI 2112 版本中開源，可以在“紅山開源平臺”下載。LES 模塊和軟件框架的接口，可參考文獻[21-22］，或者用戶手冊[41］，受篇幅限制，這里不再述及。

下一步，將在以下方面持續開展軟件功能開發和優化：①壁面模型方面，集成WMLES、RANS/LES 混合方法等；②集成人工合成湍流功能；③發展多層次、多維度的混合計算方法，例如結構/非結構網格混合計算方法，二階精度/高階精度混合計算方法，RANS/LES 混合計算方法等。

致 謝

本項研究得到了中國空氣動力研究與發展中心計算空氣動力研究所風雷課題組和馬燕凱、閔耀兵等的大力支持，在此表示感謝。