再簡單說說CFD

在之前那篇《從賽車風洞說起》後本來就計劃著寫這一篇,但是當時越寫越覺得對內容沒有把握,且又趕上開始讀博,於是就想著先暫時擱下,等自己充實一些後再動筆,沒想到一擱就是四年。一轉眼,自己也即將入職F1,於是想趁著入職前可以稍微少些顧慮地寫寫相關科普文章。在這裡先聲明本文所有的內容都源自公開信息,文中出現的所有觀點均為我個人觀點,不代表任何車隊立場。

另外在開頭需要特殊說明的是,這篇文章面向的對象是對氣動特別是賽車氣動有興趣的愛好者,因此重點會放在如何通俗易懂地介紹CFD的工作原理。換而言之對於從事相關行業的人來說,內容肯定會過於淺顯,希望大佬們輕噴。

0. CFD概述

CFD,即Computational Fluid Dynamics,中文名稱為計算流體力學,是一種使用計算機對於流體力學的控制方程(governing equation)進行數值求解的方法。數值求解的本質是放棄尋找方程的解析解,轉而通過將整個流場分為若幹個小單元並對每個小單元進行數值計算。

理想情況下,我們都希望能夠求得流場中各項參數與位置和時間的函數表達式。然而對控制方程中用來描述運動的納維-斯托克斯方程(Navier-Stokes equations)而言,其解析解是否存在作為千禧年大獎問題之一,至今仍未被解決。因此,目前對於流場的求解,受限於數學工具,隻能對某些簡化過的模型(例如二維層流管流, laminar pipe flow)計算解析解。因此數值解作為一種折中,是目前最主要的計算方法。

除CFD以外,我們還會使用一些其它的方法對於翼片或幾何體的氣動特性進行計算,例如Lumped Vortex Method、Panel Method等。與CFD對控制方程進行求解不同,這些方法多基於勢流(potential flow),即前提條件是無旋(irrotational)、無黏(inviscid)、穩態(steady)。基於勢流的計算方法主要在上世紀六、七十年代提出,雖然為瞭簡化計算對模型進行一定程度的妥協,但其具有迭代速度快的特點,因此在當下的氣動設計中,依然在預研階段扮演著一定角色。

對於CFD軟件而言,在實際運用過程中一般要分為以下幾個階段:

  1. CAD準備;
  2. 湍流模型選擇;
  3. 網格及邊界條件設計;
  4. 計算及後處理。

下面就根據上述的幾個階段做一下詳細的介紹。

1. CAD準備

在使用CFD軟件之前,首先要做的是使用CAD繪制幾何體,並對其進行表面處理。表面處理主要考慮的是兩方面問題:即繪圖質量的問題和幾何體本身的問題。

繪圖質量的問題指的是由於制圖者的繪圖習慣,導致的面與面之間形成幹涉,或產生孤立的邊或平面。具體

(a)面與面之間發生幹涉;(b)沒有厚度的自由邊;(c)三面共用一邊;(d)面面之間僅有點相連 [1]

另一種是由於幾何體本身形狀的關系,會導致貼體生成的網格形狀不好。一個最簡單的例子就是賽車的輪胎。由於在仿真軟件中,模型默認為剛體,因此當我們把輪胎放在地面上的時候,實際上輪胎與地面之間會產生一段非常細長的空間(如下圖中左側紅圈所示),在這段空間內繪制網格會導致網格畸形。一般而言,我們認為網格的長寬比不應過大或過小(實踐中當然要結合具體流場變化速率具體分析),因此這樣空間內生成的網格是不合格的。

【示意圖】輪胎與地面之間的狹長區域(左側)及其解決方法(右側)

為瞭解決這一問題,最簡單的方法是在輪胎下墊一層“臺階”,通過這樣一個結構,可以使得輪胎下方的尖銳空間被填充,從而避免瞭畸形網格的生成。當然在實際F1的模擬中,輪胎並非這樣簡單進行模擬的,很多時候還要考慮輪胎接地之後的形變等問題。

輪胎下的“臺階”

除輪胎以外,翼片的後緣(trailing edge)有時也會面臨類似問題。對於尖銳的後緣,通常使用的方法是將尖銳部分截去(如下圖所示)。看上去截去之後會帶來一個很誇張的鈍緣(blunt edge),但實際上這部分的厚度一般在幾個毫米,與翼片的弦長(100mm以上)相較而言並不十分誇張。

截去尖銳後緣後的翼片末端

2. 湍流模型

在這裡我們主要關註的是基於有限體積法和納維斯托克方程的湍流模型,有限差分、有限元、格子玻爾茲曼、光滑粒子等方法不在此討論中。目前在對湍流的處理上,主流方法(不局限工業界)主要分為三類——直接數值模擬、大渦模擬和雷諾平均方程模擬。

2.1 直接數值模擬(Direct Numerical Simulation,DNS)

理論上,我們可以直接通過邊界條件和時間不長,對NS方程進行求解,這就是所謂的“直接數值模擬”。這種方法需要能夠模擬所有的湍流尺度——直到Kolmogorov scale,因此對於計算能力有極高的要求。一般認為DNS的計算力需求隨雷諾數呈三次方增長(~mathrm{Re}^3

2.2 大渦模擬(Large Eddy Simulation,LES)

相較於DNS,LES是使用濾波函數將大尺度的渦和小尺度的渦分離開,對大尺度的渦直接模擬,並對小尺度的渦使用模型來封閉的一種方法。LES建立在Kolmogorov假設的基礎上——即湍流在耗散過程中,首先通過黏性在積分紊流尺度(integral length scale,最大的湍流尺度)上將平流中的動能提取出來,並向更小Taylor microscale傳遞,在這一尺度上湍流能量並不會耗散,而是直接傳遞給最小的Kolmogorov length scale,並在這一尺度上化為內能耗散掉。一般認為在Taylor microscale這一尺度一下,湍流一般會表現出相似性,即在統計學上遵循相似的模式,不隨流場具體工況變化而變化,因此對於這一部分的湍流的計算可以用模型替代。LES對於計算力的需求小於DNS,但是要大於RANS,因此當前階段在工業界中運用也相對較少,但其與RANS的混合方法(hybrid RANS/DES method)近幾年在工業界中運用較為廣泛。

2.3. 雷諾平均方程模擬(Reynolds-averaged Navier-Stokes equations,RANS)

RANS是三種主流方法中對於計算力需求最小的一個,它通過對NS方程在時間上進行平均,避免瞭對小尺度脈沖進行計算,從而大大降低瞭對於計算力的需求。不過,由於在平均過程中引入瞭雷諾應力從而產生瞭關於方程封閉性的問題(即方程數小於未知數),因此使用湍流模型對雷諾應力進行模擬。在F1中,RANS是最經常被用到方法,不過相較於標準的RANS模型,各個車隊一般不僅會對RANS重新標定,更會使用各種延伸的改進型(如lag模型)

3. 網格

在完成CAD準備的基礎上,就要考慮網格的問題瞭。網格分類的方式有很多,總的來說分為兩類——結構網格與非結構網格。結構網格是指網格中每個節點都連接著四個相鄰的網格,而非結構網格即為不滿足上述條件的網格(如下圖)。

二維(a)結構網格,與(b)非結構網格 [2]

而對於任何一種網格,在具體設計中都要考慮兩部分:面網格(surface mesh)和體網格(volume mesh)。面網格是指如何用網格劃分幾何體的表面,一定程度上表示瞭計算機眼中的幾何體長什麼樣子。因此對面網格而言,需要仔細選取合適的解析度,既不能太稀疏,因為會喪失幾何體的特性(比如用8個點畫球就變成瞭立方體);也不能太密,因為會讓後續的體網格過大(比如正36邊形在很多情況下就和圓形差不多瞭)。

體網格(volume mesh)的設計一般分為兩部分:棱柱層(prism layer),以及核心網格(core mesh)。註意這裡“核心網格”是西門子軟件中的術語,用來指代體網格中棱柱層以外的部分,不同情況可能有不同稱呼。

棱柱層的用意使用規整的網格來包裹需要準確計算的邊界層(通常在賽車車身和靠近車身的地面),以期更加準確地對分離等現象進行預測。關於邊界層是什麼,完全沒有概念的同學可以參考我之前那篇文章《從賽車風洞說起》的第二部分,裡面有一些比較直觀的介紹。由於流場在這一區域,隨著與幾何體邊界的距離增大而劇烈變化,因此用網格捕捉這一部分的變化尤其重要。一般設計當中需要考慮幾個參數:第一層網格高度(即y_1^+),棱柱層總高度,以及這一區間內的網格增長率。 y_1^+基本取決於湍流模型的選擇:如果使用壁函數(wall function)則需要將第一層網格置於對數區(log-law region, y_1^+>30 )或粘性底層(viscous sublayer, y_1^+<5 );如果不使用壁函數,則一般需要另 y_1^+approx1 。棱柱層厚度則一般取決於邊界層的厚度,通常需要包裹整個邊界層。而增長率則是用來確保相鄰網格的大小變化不會過大,以保證有足夠的精度來計算速度梯度極高的邊界層。

prism layer 示意圖 (註意一般最外層與體網格的體積比不會如圖所示這麼大)邊界層中速度隨壁面距離變化情況 [3]

在核心網格(core mesh)部分,主要需要考慮的是在何處進行網格加密(mesh refinement)。受到規則限制(關於CFD和風洞測試規則部分,如果感興趣的人多我可以再另外寫一篇),普遍方法是對半輛車(鏡面對稱)做穩態RANS仿真,一般網格大小大致在1億(100 million)左右 [4]。加密區的設計則更多是經驗導向的,並無一定之規。在設計時需要對氣流走向有一個大致的判斷——氣流會在哪裡分離,哪裡的流動會比較劇烈。一般而言需要對速度梯度大的地方(例如剪切層,shear layer)進行加密,以便更好地捕捉流場中小尺度的律動——這是由於在速度梯度大的地方一般有更高的湍流強度(湍流的生成,production,與速度梯度正相關),也因此會有更加復雜的流場結構。這裡順便提一句,對於仿真感興趣的同學可以參考這篇文章,作者是一位即將到梅奔做Graduate Aerodynamist的人寫的。雖然技術上還有很多值得討論的地方,但是對圍場外的大傢是個不錯的參考。

4. 計算及後處理

在F1當中,本文所述的網格和湍流模型部分都是由CFD部門通過代碼預先設置好的,空氣動力學工程師們(比如 @Next 大神)一般隻需要提交需要進行仿真的幾何體,程序就會自動判定需要如何劃分網格,使用多少核進行並行計算,並自動運行後處理內容(計算包括但不限於各種力的系數、表面壓力分佈、各種切面的雲圖、以及各種等值面)。這裡稍微列舉幾個常見的參數:

  • 力的系數:升力/下壓力(down force)、阻力(drag)和側向力(side force)。一般使用公式 C_F=frac{F}{frac{1}{2}rho U_infty^2A} 計算,其中 F 是力、 rho 是密度、 A 是參考面積(一般定義為幾何體的正投影,在實際使用中可能會固定大小以比較力的大小)
  • 速度場:一般會關註速度大小(velocity magnitude, U_mathrm{mag}=sqrt{U_x^2+U_y^2+U_z^2} ),和流向速度 U_x
  • 壓力場:壓力一般分為動壓(dynamic pressure)、靜壓(static pressure)和總壓(total pressure),初學者們很容易分不清三者的區別。總壓是靜壓和動壓之和;靜壓是流體的內稟屬性,源自分子的無規則熱運動,在實際測量中一般需要將管口垂直於流體流向;而動壓是流體由於運動產生的壓力,在實際測量中很難直接進行測量,通常是使用皮托靜壓管測量總壓和靜壓後相減得到。
  • 湍流動能:又稱TKE(turbulence kinetic energy),一般分為模擬(modelled)和解析(resolved)兩部分。在RANS模擬,特別是steady RANS中,通常隻有模擬的部分。湍流動能和流場結構息息相關,因此在CFD中非常重要如何準確計算湍流動能十分重要。
  • 渦:渦(vorticity)一般是指流體沿某一軸旋轉的運動。在數學上,渦量一般定義為 mathbf{omega}=nablatimesmathbf{U} . 但在實際運用中,由於所有帶旋轉的運動都會反應在渦量上,因此很多時候無法將渦與剪切層(shear layer)區分開,所以在實際應用上很多時候會使用 lambda_2Qmathrm{-criterion}

5. 總結

不知不覺寫瞭快5000字,但是內容還是比較淺顯而且感覺和F1的關聯並不大(笑),更多是籠統的CFD介紹,以期能讓愛好者們有個大致的瞭解。事實上這裡的每一步都有細化深挖的空間:比如網格究竟需要多密;湍流模型具體是怎麼結構、如何工作的;自動化CFD程序需要如何編排;並行計算要如何設置等等。如果後續有機會可以單獨就每一部分詳談,同時也歡迎大傢評論交流。

[1] CD-adapco STAR-CCM+ 用戶手冊 (侵刪)(PDF) Localized Coarsening of Conforming All-Hexahedral Meshes

[2]

[3] What is y+ (yplus)?

[4]

赞(0)