做CFD的都知道算不可壓縮流的時候最讓人頭疼的就是壓力和速度之間的糾纏關系。動量方程里有壓力梯度項但壓力本身卻沒有獨立的控制方程只能通過連續(xù)性方程來約束速度場。這就導致我們在求解時不能直接一次性得到答案必須走預測-修正的迭代路線。而在這一堆修正算法里PISO和SIMPLEC是我個人用得最多的兩個也恰恰是很多新手容易混淆的兩個。這篇就從原理到實操把這兩兄弟徹底掰開揉碎講清楚。很多人一開始接觸PISO和SIMPLEC都是因為OpenFOAM或者Fluent里的求解器設置選項。但說實話如果不知道它們背后的修正邏輯你連殘差曲線都看不懂更別說針對自己的算例去調算法參數(shù)了。這篇文章我會先用大白話解釋為什么不可壓縮流必須要用這類算法再分別拆解SIMPLEC和PISO的修正步驟最后結合我在實際項目中遇到的收斂問題給你一套可以直接用的選型和排錯思路。1. 壓力-速度耦合為什么SIMPLE家族如此重要1.1 壓力沒有自己的方程先回憶一下不可壓縮Navier-Stokes方程。連續(xù)性方程是速度散度等于零動量方程是那個包含壓力梯度、對流項和擴散項的偏微分方程。如果直接離散這個方程組你會發(fā)現(xiàn)矩陣里壓力項出現(xiàn)在動量方程中但是壓力本身沒有任何時間導數(shù)項也沒有獨立的演化方程。這就意味著壓力更像是一個拉格朗日乘子它的作用是在每一步強制速度場滿足連續(xù)性約束。如果我們把壓力項去掉速度場就會變成純對流擴散問題解出來根本滿足不了質量守恒。所以壓力場的本質是通過連續(xù)性約束反算出來的。這就帶來一個麻煩離散后的代數(shù)方程是耦合的而且壓力矩陣常常是奇異或半正定的直接做聯(lián)立求解在大多數(shù)非結構網(wǎng)格上都非常困難內存和時間開銷也大得離譜。因此工程上幾乎沒有人去做完全的隱式聯(lián)立求解大家都走分離式求解這條路。分離式求解的核心思想就是先猜一個壓力用動量方程算出速度再修正壓力讓速度滿足連續(xù)方程。這個猜-算-修的循環(huán)就是SIMPLE家族的共同基因。PISO和SIMPLEC的區(qū)別不過是在這個循環(huán)里多修幾步、少修幾步、用什么近似來修正而已。1.2 分離式求解的通用框架所有SIMPLE類算法都可以總結成一個大框架用當前的壓力場求解動量方程得到一個不滿足連續(xù)性的預測速度。通過壓力修正方程算出一個修正量。用修正量更新壓力和速度。判斷收斂不收斂就重復。這個框架里最關鍵的一步是建立速度修正量和壓力修正量之間的線性關系。SIMPLE算法在這里做了一個近似把鄰居單元的速度修正直接丟掉只保留本單元的修正。這個近似對對角線占優(yōu)有利但也帶來了收斂慢的問題。SIMPLEC和PISO的差異本質上就是對這個近似的處理方式不同一個是從系數(shù)上補償一個是從迭代結構上補償。1.3 為什么不能多迭代幾次就完事有人會問既然都是迭代那我用最簡單的SIMPLE算法多迭代幾次效果不也一樣嗎理論上對于穩(wěn)態(tài)問題多迭代確實能收斂到同一結果。但實際問題講究的是收斂速度和穩(wěn)定性。SIMPLE算法的這個近似在網(wǎng)格長寬比很大的邊界層網(wǎng)格上會讓壓力修正量算得過猛或過偏導致殘差震蕩甚至發(fā)散。這時候你加再多迭代步數(shù)都沒用需要從算法層面去改進這也就是SIMPLEC和PISO存在的意義。2. SIMPLEC算法拆解2.1 從SIMPLE到SIMPLEC的一步之遙SIMPLEC的全稱是SIMPLE-Consistent意思是一致的SIMPLE。它和SIMPLE的唯一區(qū)別就在處理鄰居速度修正的時候不再完全丟掉了。SIMPLE算法假設鄰居單元的速度修正為零這樣壓力修正方程的形式很簡單但代價是每輪修正量偏小需要多迭代。SIMPLEC的改進思路是既然鄰居速度修正不能精確算出來那我就不假設它為零而是把它近似地表達成本單元速度修正的某種比例。具體來說在推導壓力修正方程時把鄰居修正項移項合并到本單元的系數(shù)里這樣一來求解出的壓力修正量代表了整個單元塊的平均修正趨勢更接近真實解。雖然壓力修正方程的非對角項變復雜了但在很多算例里一次修正就能達到SIMPLE兩三次修正的效果。2.2 SIMPLEC的修正方程怎么說這里我用一個簡化的二維不可壓縮流動來說明。離散后的動量方程可以寫成a_P * u_P Σ a_N * u_N - (?p/?x)_P * V_P b_P其中a_P是本單元系數(shù)a_N是鄰居系數(shù)u_N是鄰居速度b_P是源項。引入壓力修正 p 和速度修正 u 之后速度修正滿足類似的關系a_P * u_P Σ a_N * u_N - (?p/?x)_P * V_PSIMPLE的做法是忽略Σ a_N * u_N這一項于是u_P - (V_P / a_P) * (?p/?x)_PSIMPLEC不忽略鄰居項而是用一個近似把Σ a_N * u_N表示成 u_P 乘以某種系數(shù)加權。這樣處理之后速度修正u_P的表達式分母會從 a_P 變成 a_P - Σ a_N這個修正量在數(shù)值上比SIMPLE更大修正效果也更直接。如果你去看OpenFOAM的源碼會發(fā)現(xiàn)SIMPLEC的壓力修正方程系數(shù)里面確實多了一個求和項就是干這個用的。2.3 SIMPLEC適合哪些情況在我的經(jīng)驗里SIMPLEC在穩(wěn)態(tài)問題里的表現(xiàn)往往比SIMPLE好尤其當網(wǎng)格質量一般、存在大長寬比單元的時候。因為它的修正路徑更直接不容易因為鄰居速度被忽略而出現(xiàn)修正滯后。對于不可壓縮穩(wěn)態(tài)流、傳熱共軛問題、多孔介質流動這類計算SIMPLEC是很好的默認選擇。不過也要注意SIMPLEC對松弛因子的取值比SIMPLE更敏感。因為它修正得比較猛如果壓力松弛因子給得太大很容易在迭代初期產(chǎn)生壓力場的過沖。我常用的做法是壓力松弛因子取0.7到1.0之間速度松弛因子取0.7到0.8之間具體還要看算例的非線性強度。記得先小后大等殘差有下降趨勢了再逐步提升。3. PISO算法拆解3.1 兩步修正的來歷PISO的全稱是Pressure-Implicit with Splitting of Operators最早是為了模擬非穩(wěn)態(tài)可壓縮流提出的后來被廣泛用在不可壓縮瞬態(tài)計算中。PISO和SIMPLECAC最直觀的區(qū)別是它在每輪迭代里做了兩次壓力修正而不是一次。PISO的計算流程是這樣預測步用當前壓力求解動量方程得到預測速度。第一次修正構造并求解壓力修正方程更新壓力和速度。第二次修正重新計算速度修正與壓力修正之間的耦合關系再次求解壓力修正方程。這里的第二次修正并不是簡單重復而是把第一次修正后仍然存在的速度散度殘余再處理一遍。因為第一次修正方程本身也是基于近似得到的所以總會有殘差。PISO的第二步修正能把這些殘差再壓下去一截得到更滿足連續(xù)性的速度場。3.2 瞬態(tài)計算中的時間步優(yōu)勢PISO在瞬態(tài)計算里的優(yōu)勢主要體現(xiàn)在它可以配合更大的庫朗數(shù)Courant number。由于每一步內做了兩次修正速度場在每個時間步內更接近滿足連續(xù)性約束因此即使對流項比較強也不容易因為壓力-速度失配而產(chǎn)生震蕩。我用它算過圓柱繞流和機翼附近的非定常流動在庫朗數(shù)接近0.8-1.0的時候PISO依然能維持穩(wěn)定的殘差下降而SIMPLEC如果庫朗數(shù)超過0.5就有點危險了。當然PISO也不是萬能的。它的兩次修正會讓每一步的計算量比SIMPLEC大一些但這個代價換來了更大的穩(wěn)定時間步在瞬態(tài)問題里往往非常劃算。如果你的算例時間依賴性很強比如渦街脫落、顫振這類現(xiàn)象PISO基本是首選。3.3 PISO的邊界條件敏感性PISO算法對邊界條件的處理比SIMPLEC更敏感尤其是出口邊界。因為第二次修正會重新計算壓力場和速度場的耦合關系如果出口邊界的壓力設定不合理很容易出現(xiàn)回流時壓力場的振蕩。我在實際項目中就遇到過出口設置為固定壓力但計算域出口附近有強回流的情況PISO算出來的出口附近壓力出現(xiàn)了鋸齒狀分布后來把出口段加長作為緩沖域才把這個問題壓下去。還有一點PISO在面對非正交網(wǎng)格時非正交修正項的處理要格外小心。OpenFOAM里用PISO必須注意nonOrthogonalCorrectors的設置如果網(wǎng)格非正交度高建議把它設為3到5次否則會看到壓力殘差先降后升怎么都壓不下去。這一點在很多新手教程里都沒提我踩過坑之后就特別注意了。4. PISO與SIMPLEC的橫向對比4.1 修正邏輯的對比SIMPLEC和PISO雖然都是為了克服SIMPLE收斂慢的問題但它們走的是兩條完全不同的路SIMPLEC是從修正量計算入手通過保留鄰居速度修正項來讓每輪修正更準確從而減少迭代輪數(shù)PISO是從修正流程入手通過一輪預測加兩輪修正讓同一預測步內的壓力-速度耦合更緊密從而減少時間步內誤差。打個通俗的比方SIMPLEC像是一個判斷力很強的醫(yī)生一次診斷就基本找準病灶PISO像是一個嚴謹?shù)尼t(yī)生看一次不放心復查一次再確認。對于慢性病穩(wěn)態(tài)問題SIMPLEC效率更高對于急性發(fā)作瞬態(tài)問題PISO的復查機制更保險。4.2 參數(shù)與適用場景速查這里我整理了一個在實際選型時經(jīng)常參考的對比表對比維度SIMPLECPISO修正次數(shù)單次修正兩次修正主要適用場景穩(wěn)態(tài)問題瞬態(tài)問題對大長寬比網(wǎng)格的容忍度高中高對庫朗數(shù)的承受能力較低建議小于0.5較高可到0.8-1.0每步計算量小中壓力松弛因子建議0.7-1.01.0一般不用松弛速度松弛因子建議0.7-0.81.0對出口邊界條件敏感性中高需要注意的是現(xiàn)在的商業(yè)CFD軟件里Fluent對PISO也做了很多細節(jié)優(yōu)化比如自動開啟偏斜修正和鄰近修正所以它的適用性比教科書里說的更廣。但在OpenFOAM這類開源工具里你還是需要手動去配置這些修正項的這時候理解原理就特別重要。4.3 選型決策先問自己三個問題如果你現(xiàn)在要開始一個新的算例不知道怎么選建議先問自己三個問題第一是穩(wěn)態(tài)還是瞬態(tài)穩(wěn)態(tài)首選SIMPLEC瞬態(tài)首選PISO。第二網(wǎng)格質量如何如果網(wǎng)格長寬比大、非正交度高SIMPLEC更能抗干擾PISO就需要多開非正交修正。第三計算資源緊張嗎資源緊張、時間步數(shù)多PISO的兩次修正雖然每步貴一點但穩(wěn)定時間步可以放大總體算下來往往更省。我用這套思路處理過很多通風管道、換熱器、閥門內部流場的項目大部分情況下都能在第一次試算時就選對算法省去了反復試錯的時間。5. 實戰(zhàn)在OpenFOAM和Fluent中的設置與調參5.1 OpenFOAM中的算法選擇在OpenFOAM里選擇算法本質上就是選擇用哪個求解器icoFoam和pimpleFoam里你可以通過fvSolution文件中的solver和PISO關鍵字來控制算法。舉個例子對于瞬態(tài)不可壓縮流動如果要用PISO你會在fvSolution里看到類似這樣的配置PISO { nCorrectors 2; nNonOrthogonalCorrectors 2; pRefCell 0; pRefValue 0; }nCorrectors對應的是每一步內壓力修正的次數(shù)PISO一般設2到3次設太大了收益也不明顯。nNonOrthogonalCorrectors對應非正交修正的次數(shù)網(wǎng)格越歪這個值越要大但每次非正交修正都要額外求解一次壓力泊松方程計算量會增加所以網(wǎng)格質量好的話保持1到2就夠。如果你想用SIMPLEC做穩(wěn)態(tài)計算可以用simpleFoam自帶的SIMPLE算法也可以手動調整相關參數(shù)達到SIMPLEC的效果。simpleFoam里的SIMPLE子字典中有一個consistent選項OpenFOAM早期版本可以通過設置這個參數(shù)來啟用SIMPLEC模式新版則通過不同的求解類名來實現(xiàn)。如果你用的版本支持直接SIMPLE { nNonOrthogonalCorrectors 0; consistent yes; }5.2 Fluent中的面板配置在Fluent里操作就更圖形化了。打開Solver設置面板在Solution Methods頁你可以看到Scheme下拉菜單里有SIMPLE、SIMPLEC和PISO三個選項。對于穩(wěn)態(tài)計算我一般在Fluent里選SIMPLEC并把Skewness Correction和Skewness-Neighbor Coupling兩個選項打開能顯著改善非結構網(wǎng)格上的收斂性。注意面板上有一個Skewness Correction下拉菜單默認可能是0建議根據(jù)網(wǎng)格偏斜程度選1或2。對于瞬態(tài)計算選PISO后會出現(xiàn)兩個子選項Skewness Correction和Neighbor Correction。如果網(wǎng)格質量不錯Neighbor Correction保持1就夠如果局部網(wǎng)格有畸變Skewness Correction可以開到2或3。我遇到過一個復雜裝配體的流動局部四面體網(wǎng)格質量很好但偏斜不小開Skewness Correction 3之后壓力殘差直接從振蕩變?yōu)槠巾樝陆?。在Fluent里還要注意選定算法后亞松弛因子并不總是在面板上直接暴露。SIMPLEC算法下如果你在Solution Controls里看到Pressure和Momentum的松弛因子建議初始值設為0.7和0.7不要貿(mào)然調到1.0。雖然從原理上SIMPLEC可以不用松弛但實際算例里非線性源項的影響還是需要一點過渡的。5.3 離散格式的配合算法選好了離散格式也不能隨便配。我的經(jīng)驗是算法與離散格式的匹配優(yōu)先級是先保證收斂再追求精度。對SIMPLEC動量方程建議先用一階迎風格式迭代幾百步讓殘差降下來再切到二階迎風格式繼續(xù)算。否則一開始就用高階格式流場還在劇烈調整中很容易在激波或回流區(qū)出現(xiàn)振蕩。對PISO如果庫朗數(shù)已經(jīng)比較高空間格式也盡量用有界格式比如Fluent里的Third-Order MUSCL或OpenFOAM里的linearUpwind能減少數(shù)值振蕩。壓力插值格式也很關鍵。OpenFOAM里默認的動量插值就是為壓力-速度耦合設計的不要隨便改成線性插值。Fluent里壓力插值一般選Second Order但如果你的算例旋轉或強體積力主導可以考慮PRESTO格式。5.4 松弛因子和時間步長的聯(lián)調最后說說松弛因子和時間步長的關系。穩(wěn)態(tài)算法SIMPLEC沒有時間步的概念所以只能靠松弛因子控制收斂速率。瞬態(tài)算法PISO有時間步松弛因子一般就保持1.0不能靠降低松弛因子來救場只能縮小時間步。這里有一個我常用的經(jīng)驗公式先估算你物理尺度下的對流時間尺度初始時間步設為它的百分之一比如特征長度0.1米、流速1米每秒對流時間尺度是0.1秒那初始時間步就取0.001秒。算起來之后看庫朗數(shù)如果最大值超過1.0就把時間步減半。這種保守起步、動態(tài)調整的方式能在穩(wěn)定性和計算效率之間找到很好的平衡。6. 收斂問題排查實錄6.1 壓力場出現(xiàn)棋盤振蕩如果殘差曲線一路下降但壓力場看起來像跳棋盤一樣一格高一格低這大概率是壓力-速度失耦問題。這個問題在SIMPLEC算法中更常見尤其是使用同位網(wǎng)格collocated grid的時候。雖然現(xiàn)在的代碼都用Rhie-Chow插值來抑制壓力振蕩但如果插值系數(shù)設置不當還是會出現(xiàn)spurious oscillation。排查思路是這樣先檢查壓力松弛因子是不是偏高SIMPLEC下控制在0.7左右再檢查動量方程離散格式如果是一階格式可以切到二階試試最后檢查網(wǎng)格質量那些長寬比特別大的單元往往是振蕩的發(fā)源地。我在一次管殼式換熱器模擬中就是因為局部網(wǎng)格長寬比超過2000壓力場總是振蕩后來重新切網(wǎng)格把長寬比壓到500以下問題才徹底解決。6.2 殘差停滯在某個水平降不下去這種情況在SIMPLEC算法里很常見通常表現(xiàn)為殘差在前幾百步大幅下降然后卡在1e-3就不再動了。我的第一反應是檢查通量平衡。在Fluent里看報告里的Net Imbalance如果進出口質量流量差超過千分之一說明問題出在邊界條件或網(wǎng)格上而不是算法本身。第二反應是檢查有沒有回流穿透出口邊界。如果出口設置成壓力出口但回流區(qū)域很大SIMPLEC求解器會在邊界上來回調整壓力導致殘差波動。解決辦法是加長出口段或者改用outflow邊界條件讓流動充分發(fā)展后再出域。還有一招對穩(wěn)態(tài)問題可以在迭代過程中切換算法。先用SIMPLEC或PISO先跑幾百步把流場大致建立起來再切回SIMPLE繼續(xù)迭代到收斂。這招對很多強非線性算例有奇效因為PISO和SIMPLEC的前期修正能力更強能更快把流場拖出初始猜測的誤區(qū)。6.3 PISO瞬態(tài)計算中庫朗數(shù)過大引起發(fā)散PISO可以承受更高的庫朗數(shù)但也不是無上限。如果你發(fā)現(xiàn)某個時間步內殘差突然翻倍或發(fā)散十有八九是局部網(wǎng)格太細導致庫朗數(shù)局部爆表。PISO雖然做了兩次修正但如果庫朗數(shù)超過2很多情況下也壓不住。我的排查經(jīng)驗是先用后處理看一下每個網(wǎng)格單元的庫朗數(shù)云圖找到超標的局部區(qū)域。如果只是少數(shù)小網(wǎng)格引起最簡單的辦法是局部加密這些區(qū)域的網(wǎng)格把庫朗數(shù)降下來比全局縮小時間步劃算得多。如果整場庫朗數(shù)都高那就老老實實時間步減半。還有一次我遇到過一種奇怪情況PISO算法在時間推進幾千步后壓力殘差突然上升。最后發(fā)現(xiàn)是物性參數(shù)在某個溫度區(qū)間內突變導致密度不可壓縮假設失效。這時候選什么算法都沒用得先確認你的流動前提是否成立。不可壓縮求解器算到接近可壓縮的區(qū)域必然發(fā)散這不是算法的問題是模型選擇的問題。6.4 非正交修正次數(shù)該怎么調在OpenFOAM里用PISO時nNonOrthogonalCorrectors設置不當是很常見的問題。這個參數(shù)設少了非正交網(wǎng)格上的壓力梯度計算不準殘差會反復跳設多了每一步要多解幾次壓力方程計算時間明顯增加。我的調試方法是先看網(wǎng)格的非正交性報告OpenFOAM運行checkMesh會給出最大非正交角度。如果最大非正交角小于60度設1次就夠60到75度設2到3次超過75度建議先修網(wǎng)格不要光靠調算法參數(shù)硬扛。因為算法修正只能緩解問題不能根治網(wǎng)格質量問題。如果你發(fā)現(xiàn)壓力方程迭代次數(shù)總是設置次數(shù)上限且殘差還降不到設定值就要考慮是不是非正交修正次數(shù)不夠。試著逐步加一次看殘差曲線是否有明顯改善有改善就繼續(xù)加沒改善就維持現(xiàn)狀不要盲目加太多。這類問題排查其實挺考驗耐心的但每次通過調參和網(wǎng)格優(yōu)化把殘差曲線從震蕩拉到平滑下降我都會對算法本身的特性多一些具體認識。PISO和SIMPLEC這兩個算法表面上只是幾個字母的差別實際用起來卻有完全不一樣的脾氣。根據(jù)我個人經(jīng)驗新手最忌諱的是上來就追求高階格式和大庫朗數(shù)那等于把算法的所有缺點都暴露出來了。反而是先用一階格式、小步長、保守松弛因子把流程跑通再逐步釋放參數(shù)才能穩(wěn)扎穩(wěn)打地把問題收斂到底。希望這篇能幫你少走一些我當年走過的彎路。