邊界層理論及壁面方法
1. 問題
邊界層處理是所有CFD工程師都面臨的一個難題，
第一層取多厚？
y+值是多少？
做多少層？
層與層間的增長比例？
邊界層總厚度多少？
與試驗值不符合，要不要再加密一下……種種問題不一而足。
邊界層理論在所有流體力學著作里都有涉及，最專業的著作當屬德國航空專家H. Schlichting的《Boundary-Layer Theory》，這本專著已經出到第八版了，足見其影響力，目前亞馬遜有售。本文主要內容都是筆者的工程經驗，理論部分主要取自上述專著，難免有疏漏或不當之處，望方家斧正。
2. 理論
邊界層的概念是1904年德國著名的力學家普朗特在海德爾堡第三屆國際數學家學會上宣讀的“關于摩擦極小的流體運動”的論文中首先提出的。他根據理論研究和實際觀察，證實了對于水和空氣等粘性系數很小的流體，在大雷諾數下繞物體流動時，粘性對流動的影響僅限于緊貼物體壁面的薄層中，而在這一薄層外粘性的影響很小，完全可以忽略不計。普朗特把這薄層稱為邊界層，或稱附面層。從邊界層內的流動過渡到外部流動是漸變的，所以邊界層的厚度δ通常定義為從物面到約等于99%的外部流動速度處的垂直距離，它隨著離物體前緣的距離增加而增大。根據雷諾數的大小，邊界層內的流動有層流與湍流兩種形態。一般上游為層流邊界層，下游從某處以后轉變為湍流，且邊界層急劇增厚。層流和湍流之間有一過渡區，稱為轉捩（Transition）現象，NASA AmesCenter曾在導彈風洞試驗中復現這一現象，如圖1所示。
圖1 風洞試驗中的層流-湍流轉捩現象
圖1很直觀的給大家展現了層流到湍流的轉捩過程，我們對湍流區域的邊界層進一步放大，可以歸納出另一張經典的湍流邊界層分區，如圖2所示。
圖2 湍流邊界層分區
圖2使用歸一化坐標而非物理坐標y，這是因為，使用物理坐標時，隨著雷諾數Re的變化，其速度值是變化的，而研究中更需要尋找出不變性，通過歸一化處理后，消除了Re顯示變化影響（吸收到摩擦速度和粘性長度的定義里），因此得到了普適的結果。
下面對圖2的物理量進行解釋，橫坐標y+的定義為
縱坐標
其中y是離壁面的距離，ut可以理解為壁面附近流體的剪切速度，v是運動粘度， u代表主流速度。細心的讀者應該可以注意到，這y+其實是個雷諾數的公式，所以y+是y處漩渦的典型雷諾數，也反映了粘性影響隨y的變化。根據y+值的大小，可以將邊界層內分割為三個區域：
1) 粘性底層（Viscous sublayer）：??+ < 5，粘性底層內，u+ = y+，由于橫坐標y+采用對數坐標，所以在圖2上顯示為曲線；
2) 過渡區（Buffer layer）：5 ≤ ??+ ≤ 30；
3) 對數區（Log-law region）：30 < ??+ < ??(??e)，對數區內u+ = 2.5*ln（y+）+ 5.45，在對數坐標下為一直線。該區域的上限即邊界層總厚度，隨Re變化，一般認為是99%的外部流動速度處，定義為δtot，對于平板流動，Schlichting給出了經驗估算方法，
3. STAR-CCM+壁面模型
根據壁面y+值及邊界層網格分布，發展出了不同的壁面模型，如圖3所示。
圖3 不同y+壁面算法示意
STAR-CCM+壁面模型的特點是自動及全面，如圖4所示，具體應用技巧如下：
1) Low y+ walltreatment：y+值小于1, 一般推薦25+層，其中3+層位于粘性底層內，5+層位于過渡區，5+層位于對數區內，直接求解近壁面速度分布，主要用于外氣動計算中的流動分離，轉捩現象；
2) High y+wall treatment：y+值位于30至150之間，5層左右邊界層，主要用于無流動分離的情況；
3) 默認選項為All y+ wall treatment，顧名思義，不需要用戶選擇，可以自動根據來流速度，邊界網格計算y+值，并自動切換壁面算法，使用時注意y+<5或30<y+<150，盡量避開過渡區。
圖4 STAR-CCM+壁面模型
錯誤的試用壁面模型可能使計算結果嚴重失真，圖5為一管道流動案例，分別采用y+=1（紅色實線），y+=20（藍色實線）獲得的計算結果與DNS（Direct Numerical Simulation，可認為真值，圖中圓點）結果對比，可以看到，在橫坐標y+值小于200的區間內，y+=1的U+與軟流粘性系數計算結果與DNS結果吻合度極高，而藍色的y+=20位于過渡區內，現有的壁面模型處理該區間都有缺陷，因此在計算中應盡量規避。
圖5 不同y+模型計算結果對比
4. STAR-CCM+邊界層網格處理
STAR-CCM+邊界層網格生成方法非常簡單，默認設置如圖6所示，邊界層定義方法采用Stretch factor，需要定義的參數有：Numberof Prism Layer （邊界層總層數），Prism layer stretching（相鄰層間的比例），Prism LayerTotal Thickness（邊界層總厚度），邊界第一層網格厚度通過上述三個參數計算得出。
圖6 Stretch Factor方法
如果想直接定義底層厚度，只需要在Prism Layer Mesher的屬性欄里更改生成方法，相應的參數變更為Prism Layer Near WallThickness（底層厚度），層間增長率通過計算得出。
圖7 wallThickness方法
以NASA CRM飛機模型巡航狀態工況為例（Ma=0.85，Re=5000000，特征弦長Lc=7m），若采用圖7的wallThickness方法，需要確定三個主要參數：
1. Prism Layer Near WallThickness（底層網格厚度）
可通過下列公式計算y+=1條件下的底層網格厚度，為0.037mm
2. Prism Layer TotalThickness（邊界層總厚度）
利用
計算得出
3. Number of PrismLayer （邊界層總層數）
外氣動計算一般推薦層間增長率在1.2~1.3之間，此處取增長率1.2，可計算得到大約需要35層邊界層。
文章最后展示一個很常見的誤區，許多工程師熟悉邊界層理論，也按照要求做好了網格，計算結果卻不理想，很多時候都是由圖8所示的問題引起，忽視了邊界層最外層與主流網格的過渡，導致從一個很小的邊界層網格直接過渡到主流區的大網格，體積變化率（Volume Change，衡量網格質量標準之一）可達上千。正確的做法是右圖所示，通過計算，設置邊界層總厚度為面最小網格尺寸1~2倍，即可獲得適當的過渡比例。
圖8 邊界層與主流網格過渡
文末福利，關注公眾號：新能源汽車熱管理仿真技術，回復1，領取新能源熱管理資料。回復“熱管理”領取大額優惠券。
同時本人也在技術鄰平臺更新新能源動力電池熱管理仿真和設計課程如下：如需購買可添加VX：fxy33186375領取優惠券。
1、 基于starccm+在動力電池熱管理仿真技術應用、
工程師必備
- 項目客服
- 培訓客服
- 平臺客服
TOP
