自主工業CAE軟體，關鍵在於「自主」二字。
真正掌握了核心算法技術，才算是「子彈上膛」，
從而在國際賽場上佔有一席競爭之地。
若是僅僅打著「自主技術」的幌子，
卻沒有應對技術洪流和市場危機的自主能力，
那就是對企業，甚至社會的不負責。
本文作者郭志鵬先生，會從求解器的開發過程開始，詳細講解適創科技是如何基於LBM算法開發了超快速、大尺度的流體力學求解器，從而完成了壓鑄領域模擬仿真雲計算平臺的研發。
我們的自主化技術，就是最核心的國際競爭力。
以下正文：
對於工業CAE軟體來講，其最核心的技術不是人-機交互界面、不是圖形展示、不是資料庫，而是其最底層的物理、數學算法架構，也就是我們一般稱作的求解器。
所謂求解器，指的是針對特定場景（比如液體流動、溫度傳播以及結構變形），用程序編碼的方式實現的對物理規律、數學原理的客觀還原。
圖1、工業軟體底層技術圖譜[1]
求解器在很大程度上類似於我們生活中使用的計算器，當我們在其屏幕上輸入若干數字和運算符後，計算器會自動計算並輸出結果。所有和計算相關的內部運算操作都集成在計算器的內部，使用者無法直觀地獲得對運算過程的觀察和監控——就像個「黑盒子」一樣。可以這麼認為，生活中計算器涉及到的運算實際上是廣義求解器的一個基礎和組成部分。與計算器相比，CAE軟體內部的求解器往往運算量更大，輸入參數更多，對問題的還原更為全面。
求解器的開發過程
以流體力學求解器為例，其計算結果往往是一個三維區域內部的流體特性，包括速度、壓力、形態等，隨著時間的演化過程，每一個小的局部位置在特定時間點上都涉及到大量的代數運算，而這些代數運算其實就是我們生活計算器的主要工作。與實際人手操作計算器進行計算不同，CAE軟體求解器可以根據物理規律演化的特點自動匹配當前涉及到的代數運算，這個就是我們俗稱的「求解過程」。
圖2、流體計算的核心變量：壓力，速度，溫度（熱流體）[2]
從計算器-求解器的關聯上看，開發求解器最主要的工作是把求解器的物理規律、數學算法通過編碼的方式實現。這個過程看似清晰、簡單，但卻充滿了挑戰。
依然以流體力學求解器為例，我們來看一下這個開發的過程。
首先，我們需要深入理解描述流動過程的物理原理。對於液體來說，主導或者說決定其流動的底層物理規律是動量、質量守恆（不考慮熱傳導引起的能量耗散），而通過人類幾個世紀的研究和總結，我們可以將整個流動過程用兩個核心方程來描述，一個是Navier-Stokes方程（NS方程），描述動量守恆；一個是連續性方程，描述質量守恆。人們對於物理世界的理解往往是基於三維空間的，因此，將NS方程對笛卡爾坐標系進行拆分，將速度矢量沿著每一個坐標軸（X、Y、Z）進行分解，我們就會得到三個速度變量（U、V、W），這樣的話一個NS方程就會變成針對離散速度變量的三個偏微分方程，再加上連續性方程，我們就可以形成一套完備的方程組來描述真實世界中不可壓縮流體（水、金屬液）的流動過程。
圖3、頂驅方箱流速度矢量場和標量雲圖分布[3]
有了描述流體流動過程的方程組，下一步就是用數學的方法來求解這個方程組。數學方法一般分為兩大類，一類是通過理論推導直接獲得精確的方程組的解，另一類則是通過數值算法，「近似」地獲得方程組的解。在現有的數學、物理框架下，或者說基於人類現有的科學發展狀態，絕大多數的問題是無法獲得精確解的，而解決實際問題的最主要途徑是數值計算。
數值計算一般分為三個環節：模型離散、條件設定和求解過程本身。
「模型離散」指的是將真實的物理模型打碎成一個個的網格點，模型和網格點的關係特別像人體和細胞的關係；「條件設定」一般指的是正確地設定計算的輸入參數，包括材料熱物性參數、初始條件和邊界條件等，正確的邊界條件設定對於產生準確的計算結果是極其關鍵的，要不然就是「垃圾進-垃圾出（Trash in，trash out）」了；「求解過程本身」則是編碼的數值算法在計算機硬體上的逐步、系統地運行。
數值計算的三個環節中，只有條件設定是完全依託於使用者的，模型離散和求解過程都可以且應該實現自動化。
我們把以上兩個過程（描述物理規律的偏微分方程和數值算法）聯繫起來，通過計算機編碼（高效算法編程一般使用C、C++以及fortran語言）實現成計算機可以識別的語言，也即電腦程式，最終通過特定計算機硬體平臺和作業系統就可以實現對物理問題的求解了。
值得注意的是，描述物理問題的方程組經過離散之後的維度會急劇增加，這是因為這些偏微分方程組不僅需要建立在整體物理模型上，而且必須要在每一個網格單元上都得滿足。因此，如果說一個物理模型經過網格離散獲得了1000萬個網格單元，那麼相應的描述物理現象的方程組就得有1000萬套，也就是說在每一個演化步中（一般是指時間步），我們要同時求解1000萬套方程組。仍然以流動過程為例，對於每一個單元，在每一個時間步中，我們需要求解4個核心方程，包括描述U、V和W的三個速度變量的控制方程（NS方程）和描述質量守恆的連續性方程。從這個角度上講，精確求解NS方程和連續性方程就是流動求解器需要解決的最關鍵問題。
現有的大多數工程應用涉及到的氣體、水以及金屬液等的流動都屬於粘性不可壓縮流動，針對這種流動，現有的流動求解器通常使用兩種方法來實現對控制方程組的求解：基於壓力迭代的算法和基於人造壓縮特性的算法。這兩種算法同時也是現有絕大部分商業軟體採用的方法。
圖4、基於壓力迭代的SIMPLE算法[4]
基於壓力迭代的算法又可以細分為MAC算法和分裂算子算法，詳細的算法原理可以參考相關的計算流體力學的書籍和文獻，這裡就不再贅述。MAC和分裂算子算法都可以通過數學變換最終轉換成一個聚焦於求解泊松方程的算法，而求解泊松方程本身則需要採用所謂的隱式算法不斷迭代才能完成。
事實上，求解泊松方程是近幾十年來所有流體力學求解器最難啃的「骨頭」。從數學上講，泊松方程屬於橢圓方程，每一個網格單元中心的壓力項會直接和其周圍的「鄰居」單元中心的壓力項形成關係，也即方程組，因此，對每一個單元來說，這個方程都會至少含有7個變量（當前單元自己的壓力和周圍單元的六個壓力）。採用線性方程來進行描述就是：AX = b，其中，A是一個N行和N列的係數矩陣，X是含有N個壓力項的變量，而b則是一個含有N個常量的向量。所以求解泊松方程就轉換成了求解這個方程組。事實上，A是一個稀疏矩陣，因為其組元大部分都為零，在每一行中，只有7個（或多個，取決於採用的離散格式）位置的值不為零，可以想像，如果N=10,000,000，也即含有10^7個網格，那麼這個矩陣就是一個極其稀疏的矩陣了。求解大型稀疏矩陣線性方程組是數值分析的主要任務，也是當代應用數學正在攻克的難點之一。
對於求解大型稀疏矩陣線性方程組，有興趣的讀者可以閱讀一些和數值分析或者高等數值分析相關的一些書籍，經典的算法在其中都有提及，包括共軛梯度、SVD分解、冪指數迭代等等，我這裡不做一一贅述。我想提兩個算法，一個是公認的現有最強的矢量算法，叫GMRES，另一個則是多重網格（Multigrid）算法。
GMRES算法的核心在於建立正交坐標系，然後採用坐標系的坐標來表達最終的解。用一個小維度的例子來說明一下GMRES算法，比如我們之前AX=b的方程組，如果N=10的話，那麼A就是一個10x10的係數矩陣，如果這個方程組有解的話，其解X必將是一個10x1的向量。從數學的角度理解AX=b的意思就是採用一個10維的系統將X映射到b，而這個過程其實是等效於建立正交坐標系。
GMRES在計算的過程中，絕大部分時間其實是在建立正交坐標系，因此，建立正交坐標系的時間也就最能反應GMRES算法的計算效率。在實際的工程問題上，我們發現，當計算維度不太大的時候，比方說百萬、千萬網格，GMRES算法的計算效率是非常高效的，但是當網格數量繼續增加的時候，其效率就會急劇下降，這其實是受限於計算機浮點數本身的精度的，一個正交坐標系，如果每一個向量的值都有些許誤差的話，到最後整個所謂的正交坐標系就不那么正交了，也就是相互之間不再垂直了。
而多重網格算法則沒有這個問題，所謂的多重網格指的是這個算法在求解實際問題的時候採用的是不同尺寸的網格系，而非單層的網格結構。為什麼要這麼做呢，那是因為不同尺寸的網格對於消除「計算缺陷」的貢獻是不同的，大尺寸網格消除較為粗糙的缺陷，而小尺寸網格消除的缺陷則更為「精細」，這有點像電視屏幕解析度的感覺，當電視像素解析度較低時，我們只能看到粗糙的圖像，反之則可以看到更加清晰的畫面。多重網格算法非常巧妙地利用了這點，在求解線性方程組的時候以更加高效、精準的方式來濾除計算誤差，從而極大地加快了算法本身的求解速度。
圖5、2D多重網格算法粗細層格點結構[5]
雖然GMRES和多重網格算法是極其高效的，但是在求解實際工程問題中仍然面臨極大的挑戰。首先，不管是GMRES、多重網格還是其他各種算法，其本身對方程組的奇異性非常敏感，當係數矩陣A存在一定奇異性的時候，這些方法就會很快失效。什麼是方程組的奇異性呢，其實就是矩陣A的特徵值差異太大，我們大可以不用去理會什麼是特徵值，可以這麼認為：矩陣的特徵值反應了矩陣本身的「強度」，一個係數矩陣，當其某一個特徵值在絕對值上遠遠小於其他特徵值的時候，那麼這個矩陣就是奇異的。
造成奇異的因素很多，可能是因為我們對實際問題的描述有誤，也有可能是因為邊界條件設定的不對，還有可能是這個問題本身就存在奇異性，比方說激波（shock wave）的存在就打亂了其兩側物理的連續性。那麼，當這些奇異性存在時，採用這些所謂的「隱式」迭代算法就會存在收斂性的問題，甚至在某些時候就會計算中斷，也就是我們俗稱的「發散」了。現有商業軟體開展流動過程模擬，出現計算不收斂的情況是非常常見的，當然，造成發散本身的原因也是多種多樣，有操作者自身的問題，但絕大多數其實是求解器本身沒有那麼強大，或者說魯棒性沒那麼好。
開發一個「完美」的流動求解器是極其困難的，而精確地求解Navier-Stokes方程也是當今世界7大數學難題之一（現在是6大，其中的黎曼猜想已經證實被解決了？）。這其實就是為什麼我國之前並沒有能夠開發出一款具有世界競爭力的流體力學軟體的本質原因。開發流動求解器，模擬真實自然界的各種跟流動相關的物理現象，實現真正的對工業乃至對宇宙的認識（星系爆炸、遷移和演變都涉及到流動）是現在絕大多數做應用數學、力學研究人員的終身夢想。
在真實的自然界中，跟流動相關的物理現象是非常複雜的。比如我們見到的絕大多數流動過程是涉及到多個流體的多相流，同時，當雷諾數超過2000左右的時候，流動會從層流變為紊流。多相紊流其實是自然界的常態，而真正從物理、算法的角度解決這類流動問題是極其困難的。現有商業軟體在流動求解過程中都做過大量的假設和簡化，所開發的求解器也只有在某些特定的場合才能產生相對有價值的計算和模擬結果，但是這些都不是使用者能夠確切知道的，這其實也是模擬誤差，或者說模擬精度不高的主要原因。
圖6、基於LBM算法的Xflow軟體模擬的腔體潤滑劑流動過程[6]
適創科技的自主求解器
面對這些困難，適創科技只能採用一種新的算法架構來開發流動求解器。這種新的算法就是LBM（算法本身並不新，只是相較於現有的流動求解算法來說，大規模有效使用LBM算法的時間並不長），也即格子-玻爾茲曼算法，同時採用大渦模型描述大雷諾數紊流。這兩種算法的結合可以非常有效地規避現有流動求解器遇到的核心問題（迭代求解泊松方程這一環節）。
LBM方法其實結合了有網格和無網格方法的特點。
第一，LBM採用粒子碰撞的方法從微觀來描繪基本的流動碰撞，但是這個方法不直接採用粒子碰撞來對流動狀態做描繪，而是採用密度函數的方法將粒子碰撞進行積分，在相空間裡對流態進行描述。採用密度函數的方法（而不直接採用粒子碰撞的方法）來連接微觀碰撞和宏觀流動是現代物理學的一個很重大的發現，這也是統計力學的基礎。
從這個角度講，LBM方法在基礎微觀層面採用粒子碰撞，而在宏觀流動方面採用密度函數法，其計算精度（或者說對實際物理現象的描繪）要比現有的有、無網格方法更加準確。
「有網格方法」指的是數值計算中採用網格的概念開展的算法，常見的如FDM、FEM和FVM等；而「無網格方法」則是計算中沒有採用網格這一概念，常見的如SPH、MLPG、RBF以及FPM。
但不論是哪種方法，其計算精度都依賴於「解析度」，仍然像咱們電視的像素。對於有網格的計算方法，這個解析度就是網格類型和大小，而對於無網格方法，這個解析度則是粒子的數量。顯然，我們絕對不可能用一個網格或一個粒子就把流體的流動狀態描繪出來，因此它們都有「解析度」的限制。
圖7、ARES噴射流波形幹涉模擬，4億網格，20小時
同時，為了提高計算效率，我們在數值算法上採用並行計算。LBM這種算法天生就具有極強的並行計算能力，非常適合於現有的並行計算集群，這樣的話就很好地把算法本身的優勢和現有的計算資源巧妙地結合了在一塊。值得慶幸的是，用於並行計算的集群資源越來越多，國家在這方面的投入也越來越多，因此可以極大地把計算成本壓低。
不同於大多數現有商業軟體的算法架構，適創科技開發的流動求解器，是基於LBM算法、大渦模型和並行計算為主要特點的新一代流體力學求解器，本身代表了流動計算算法的最高水平，其計算效率和解決複雜問題的能力也遠遠超越了絕大多數現有的商業軟體。
而這個求解器的誕生，也為後期適創科技拓展到其他計算應用領域奠定下了堅實的基礎，這也是我們一直引以為豪並樂以提及的：
「自主化技術，世界競爭力。」
[1] 圖片來源：
公眾號「知識自動化」：《工業軟體為什麼這麼難？》
[2] 圖片來源：
https://www.autodesk.com/solutions/simulation/cfd-fluid-flow
[3] 圖片來源：
https://www.mathworks.com/matlabcentral/fileexchange/74483-cfd101-2d-lid-driven-cavity-flow
[4] 圖片來源：
https://web.stanford.edu/class/me469b/handouts/incompressible.pdf
[5] 圖片來源：
http://web.utk.edu/~wfeng1/research.html
[6] 圖片來源：
https://scanscot.com/web/wp-content/uploads/2019/01/Media9_cropped.gif
作者：郭志鵬
畢業於清華大學（本、博），牛津大學、英國皇家學會研究會員，長期從事數位化工業方面的研究，包括高性能算法、高能X射線檢測、圖像處理以及相關工業領域的材料和核心工藝開發等，立志創造有國際競爭力的自主化CAE軟體，擺脫國際壟斷，提升和振興民族工業水平。