一種基于sph算法的流體熱運動仿真系統的制作方法
【專利摘要】本發明涉及一種基于SPH算法的流體熱運動仿真系統，屬于熱力學模擬與可視化流動領域。該方法包括如下步驟：首先根據計算實例和計算要求對流體進行建模，然后基于SPH算法求解流體控制方程和修正方程，實現流體拉格朗日熱運動的仿真模擬，最后應用可視化處理技術以圖形及動畫的方式直觀的展示流體的熱運動過程。與現有網格化CFD仿真系統相比，在模擬某些高度非線性問題(如射流問題、海洋運動問題、潰壩問題、動態邊界問題)時，會得到更為合理的仿真結果；與現有無網格CFD仿真系統相比，提升了計算穩定性和計算精度，并且填補了流體熱運動仿真和湍流仿真的空白。
【專利說明】一種基于SPH算法的流體熱運動仿真系統
【技術領域】
[0001] 本發明為一種基于光滑粒子流體動力學算法（SPH)的流體熱運動仿真系統，屬于 熱力學模擬與可視化流動領域。
【背景技術】
[0002] 流體的熱運動仿真在實際工程中有著廣泛的應用背景，例如建筑物隔熱、核反應 堆安全設計、微電子設備的冷卻等。
[0003] 光滑粒子流體動力學方法是一種基于粒子的無網格數值模擬方法.該方法被應 用于解決各種流體力學和固體力學問題，包括流固耦合問題、凝固問題、射流問題、多相流 問題、熱傳導問題等.光滑粒子動力學方法作為最經典的無網格方法，其算法被不斷完善， 多種修正技術被提出并應用，并且發展了相應的前處理與后處理技術，其已經成為較為完 善的粒子型無網格方法.
[0004] 近年來，使用傳統的網格化計算流體動力學（CFD)仿真系統在解決實際流體工程 問題時遇到了許多難點，其模擬高度非線性問題（如射流問題、海洋運動問題、潰壩問題、 動態邊界問題）時往往會受到網格變形的制約，無法實現某些高度非線性問題的模擬，這 些高度非線性問題往往需要通過采用無網格CFD方法來解決；而現有的無網格CFD仿真系 統又無法解決熱力學問題。
【發明內容】
[0005] 本發明目的是為了解決現有網格化CFD仿真系統在解決某些高度非線性問題時 因為網格的變形導致結果不夠準確以及現有無網格CFD仿真系統無法對流體的熱運動進 行仿真的問題，提供一種基于SPH算法的流體熱運動仿真系統。
[0006] 本發明技術方案的思想是：根據計算實例和計算要求，首先對流體進行建模，然后 基于SPH算法求解流體控制方程和修正方程，實現流體拉格朗日熱運動的仿真模擬，最后 應用可視化處理技術以圖形及動畫的方式直觀的展示流體的熱運動過程。
[0007] 本發明是通過以下技術方案實現的：
[0008] -種基于SPH算法的流體熱運動仿真系統，包括前處理模塊、計算模塊和后處理 模塊，整個仿真系統的模擬流程如圖1所示。
[0009] 1.前處理模塊
[0010] 所述前處理模塊，用于根據實際模擬目標建立CFD仿真模型，一般CFD仿真模型需 具備以下內容：離散方程、離散方法和離散模型，本仿真系統中離散方程為流體的流動控制 方程（質量守恒方程、動量守恒方程、能量守恒方程），尚散方法為SPH算法，尚散模型為根 據用戶要求構建的粒子模型；因此，系統將前處理模塊分為以下兩個子模塊：SPH算法設置 和粒子模型建模；
[0011] 1)SPH算法設置
[0012] 前處理模塊中的SPH算法設置將決定如何使用SPH算法，其包括以下四部分內容： 計算方式設置、SPH算法格式設置、修正方法設置和SPH算法參數設置；
[0013] a)計算方式設置決定了該仿真系統的計算方式，以充分利用計算資源，其可選類 型包括：中央處理器（CPU)并行、圖形處理器（GPU)并行；其中CPU并行通過標準消息傳遞 界面（MPI)技術實現，GPU并行通過通用并行計算架構（CUDA)實現。
[0014] b)若要使用SPH算法解算流體控制方程，則需要根據模擬對象來確定如何使用 SPH算法，即需要設置SPH算法的格式；SPH算法的格式包括：時間積分格式、核函數格式、 粘性解算格式；其中時間積分格式決定了模擬時間推進的方式，可選格式包括：蛙跳法、韋 爾萊積分法、辛積分法；核函數為SPH算法的最基本函數，其決定了 SPH算法如何離散流動 控制方程，可選格式包括：三次樣條核函數、五次樣條核函數；粘性解算格式為動量守恒方 程中流體粘性的計算方式，可選格式包括：層流粘性、人工粘性、亞粒子模型湍流粘性（湍 流模型）；
[0015] C)由于SPH算法易于產生震蕩誤差，則需要在算法中引入一些修正方法，修正 方法包括：密度重構技術，人工位移技術（Xu R，Stansby P, Laurence D.Accuracy and stability in incompressible SPH(ISPH)based on the projection method and a new approach [J]· Journal of Computational Physics, 2009, 228 (18) :6703-6725);其中密度 重構技術用于修正密度計算的震蕩誤差，可選類型包括：薛帕德濾波、S -SPH ;人工位移技 術用于修正粒子位置的計算誤差。
[0016] d)SPH算法參數包括：粒子間距離、時間步長、模擬時長、核函數的影響半徑（光滑 長度）、修正方程參數、人工聲速（人工聲速為人工狀態方程的參數，人工狀態方程用于壓 力求解）。
[0017] 2)粒子模型建模
[0018] 粒子模型建模模塊用于對粒子模型進行建模，包括粒子坐標、粒子屬性（流體或 者固體）、邊界條件和初始條件；
[0019] 所述粒子坐標的生成通過輸入相應外形參數實現，根據外形參數的幾何特征及給 定的粒子間距離，便可求得該幾何形狀內所有粒子的坐標，例如：按順序輸入矩形的四個頂 點坐標，即可生成矩形粒子域；
[0020] 所述粒子屬性設置即是將邊界粒子屬性設置為固體，其他粒子屬性設置為流體的 過程；
[0021] 所述邊界條件確定了邊界固體粒子對內部流體粒子的作用方式，可設置以下兩方 面：固體壁面邊界條件、溫度邊界條件；
[0022] 計算域邊界部分粒子的屬性為固體粒子，邊界粒子位置不變，并對靠近的流體粒 子產生一定作用力，作用力施加方式可選類型包括：邊界力法、鏡像粒子法；其中邊界力法 通過施加排斥力實現固體邊界；鏡像粒子法通過施加排斥動量實現固體邊界；
[0023] 當存在具有一定溫度特性的邊界時，該部分邊界粒子需要設置溫度邊界條件，以 使其溫度特性影響其附近粒子，溫度邊界條件可選類型包括：恒溫邊界條件、絕熱邊界條 件；恒溫邊界條件情況下，該部分固體粒子溫度恒定，持續與周圍粒子進行熱交換；絕熱邊 界條件情況下，該部分固體粒子與其他粒子不發生任何熱交換。
[0024] 所述初始條件為計算域內粒子的初始物理參數，包括：速度、壓力、溫度、參考溫 度、密度、運動粘度、重力、瑞利數、普朗特數。
[0025] 當完成了前處理模塊中的SPH算法設置和粒子模型建模，便可將數據導入計算模 塊，進行流體熱運動的解算。
[0026] 2.計算模塊
[0027] 所述計算模塊用于使用SPH算法對流體的流動控制方程和修正方程進行求解，得 到各時間步的各物理參數的數值模擬結果。該模塊首先讀入前處理模塊所建立的仿真模 型，從而獲得SPH算法設置信息和粒子模型信息，進而開始SPH計算循環，在每個計算循環 中，首先SPH算法會將流體的流動控制方程（質量守恒方程、動量守恒方程、能量守恒）施 加到計算域內所有的粒子上；然后所有粒子的物理參數會按照SPH算法離散后的控制方程 產生變化；最后當所有粒子完成一次計算后，求解前處理模塊中選取的修正方程，使需要修 正的物理參數產生修正變化，由此完成所有粒子物理參數的更新，該時間步完成；當完成粒 子狀態更新后，按照時間積分格式更新時間步，并同時輸出包含所有粒子物理參數的數據 文件，然后進入下一循環，直至時間步迭代到模擬時長結束。
[0028] 3.后處理模塊
[0029] 所述后處理模塊，用于對計算模塊輸出的仿真數據進行動態可視化處理，以直觀 展現流體熱運動過程中流體的運動過程及各個參數的變化過程。該模塊的處理過程如下： 首先讀入計算模塊輸出的數據文件，進而根據每個時間步仿真數據中的位置、速度、溫度、 密度等物理參數，將粒子位置分布以坐標點的形式表示其分布，得到整個計算域內粒子分 布圖，將其他物理參數用連續變化的顏色表示其分布，得到整個計算域內物理參數分布云 圖，最后利用所有時間步的粒子分布圖及物理參數分布云圖生成整個模擬過程的粒子運動 動畫和各物理參數變化動畫。
[0030] 有益效果
[0031] 本無網格CFD仿真系統由于采用了 SPH算法，較之已有的網格化CFD仿真系統，在 模擬某些高度非線性問題（如射流問題、海洋運動問題、潰壩問題、動態邊界問題）時，會得 到更為合理的仿真結果；
[0032] 本無網格CFD仿真系統較之已有的無網格CFD仿真系統存在如下優點：
[0033] 1)添加了能量守恒方程的求解，從而使本仿真系統得以解決熱力學問題；
[0034] 2)添加了多種修正方程的求解，從而提升了系統的計算穩定性和計算精度；
[0035] 3)添加了湍流模型，從而使本仿真系統可以解決流體熱運動湍流問題。
【專利附圖】
【附圖說明】
[0036] 圖1是本發明系統的組成結構以及數據流向示意圖；
[0037] 圖2為本發明實施例的封閉方腔自然對流物理模型；
[0038] 圖3為本發明實施例的粒子模型；
[0039] 圖4為本發明實施例某一時間步所生成的流場溫度分布云圖。
【具體實施方式】
[0040] 為了使本發明的目的、技術方案和優點更加清晰明白，下面結合實施例和附圖，對 本發明實施例做進一步的詳細說明。
[0041] 下面以封閉方腔自然對流模擬為例介紹本發明基于SPH算法的流體熱運動仿真 系統的實施過程。
[0042] 按圖1所示流程實現該實施例的仿真模擬，具體實現步驟如下：
[0043] 1.前處理模塊包括：
[0044] 1) SPH算法設置
[0045] 計算方式：GPU并行。
[0046] SPH算法格式設置如下，時間積分格式：娃跳法；核函數格式：五次樣條核函數；粘 性結算格式：層流粘性。
[0047] 修正方法：人工位移技術。
[0048] SPH算法參數：粒子間距離dx = 0. 01、時間步長10'模擬時長0. 5、核函數的影響 半徑（光滑長度）h = 1. 5dx = 0. 015、人工位移技術強度參數0. 04、人工聲速400。
[0049] 2)粒子模型建模
[0050] 圖2為封閉方腔自然對流物理模型，其中L(L= 1)表示物理模型邊長，其中 $： = 〇表示該部分粒子采用了恒溫邊界條件，τΗ和t分別表示恒高溫和恒低溫，g表示重 力的方向。
[0051] 粒子坐標和粒子屬性：以圖2中正方形的四個頂點坐標為控制點，建立如圖3所示 的粒子模型，內部正方形域設置為流體粒子，正方形域外圍設置三層邊界（固體）粒子。
[0052] 邊界條件：四個壁面均采用固體壁面邊界條件，采用邊界力法；左右壁面采用恒 溫溫度邊界條件，上下壁面采用絕熱溫度邊界條件。
[0053] 初始條件：速度0、壓力0、左側壁面溫度為1，右側壁面溫度為0、參考溫度為0. 5、 密度1. 225、運動粘度v = 10'重力g = 9. 81、瑞利數Ra = 105、普朗特數Pr = 0· 71 (以 上均為無量綱形式）。
[0054] 2.計算模塊進行以下處理：
[0055] 1)讀入前處理模塊輸入的信息；
[0056] 2)根據前處理模塊所輸入的信息，設置SPH算法與粒子模型；
[0057] 3)開始循環，根據所選計算方法與參數利用SPH算法解算質量守恒方程、動量守 恒方程、能量守恒方程和人工位移技術方程；
[0058] 4)根據解算結果，更新粒子信息，包括：速度、位置、溫度、密度；
[0059] 5)輸出該時間步的數據文件；
[0060] 6)通過時間更新進入下一循環，直至時間步迭代到模擬時長循環結束。
[0061] 3.后處理模塊進行以下處理：
[0062] 1)讀入計算模塊生成的數據；
[0063] 2)生成單個時間步的粒子分布圖及物理參數分布云圖，如圖4所示為根據某一時 間步數據生成的流場溫度分布云圖，可以看出所得溫度分布較為合理；
[0064] 3)用所有時間步的粒子分布圖及物理參數分布云圖生成整個模擬過程的粒子運 動動畫和各物理參數變化動畫，以動態展示流場的變化情況。
[0065] 以上所述僅為本發明的具體實施例，并不用于限定本發明的保護范圍，凡在本發 明的精神和原則之內，所做的任何修改、等同替換、改進等，均應包含在本發明的保護范圍 之內。
【權利要求】
1. 一種基于SPH算法的流體熱運動仿真系統，其特征在于，包括前處理模塊、計算模 塊和后處理模塊；所述前處理模塊用于根據實際模擬目標建立CFD仿真模型，所述計算模 塊用于使用SPH算法對流體的流動控制方程和修正方程進行求解，得到各時間步的仿真數 據，所述后處理模塊用于對計算模塊輸出的仿真數據進行動態可視化處理，以直觀展現流 體熱運動過程中流體的運動過程及各個參數的變化過程。
2. 根據權利要求1所述的一種基于SPH算法的流體熱運動仿真系統，其特征在于，前處 理模塊包含SPH算法設置和粒子模型建模兩個子模塊； 所述SPH算法設置模塊決定計算模塊如何使用SPH算法，包括以下內容：計算方式設 置、SPH算法格式設置、修正方法設置和SPH算法參數設置； a) 計算方式設置決定了該仿真系統的計算方式，包括中央處理器（CPU)并行和圖形處 理器（GPU)并行兩種計算方式；其中CPU并行通過標準消息傳遞界面（MPI)技術實現，GPU 并行通過通用并行計算架構（CUDA)實現； b) SPH算法格式設置包括時間積分格式、核函數格式和粘性解算格式三項內容；其中 時間積分格式決定了模擬時間推進的方式，可選格式包括：蛙跳法、韋爾萊積分法、辛積分 法；核函數為SPH算法的最基本函數，其決定了 SPH算法如何離散流動控制方程，可選格式 包括：三次樣條核函數、五次樣條核函數；粘性解算格式為動量守恒方程中流體粘性的計 算方式，可選格式包括：層流粘性、人工粘性、亞粒子模型湍流粘性； c) 由于SPH算法易于產生震蕩誤差，則需要在算法中引入一些修正方法，修正方法包 括：密度重構技術，人工位移技術；其中密度重構技術用于修正密度計算的震蕩誤差，可選 類型包括：薛帕德濾波、S -SPH ;人工位移技術用于修正粒子位置的計算誤差； d) SPH算法參數包括：粒子間距離、時間步長、模擬時長、核函數的影響半徑、修正方程 參數、和人工聲速； 所述粒子模型建模模塊用于對粒子模型進行建模，其包括以下四部分內容：粒子坐標、 粒子屬性、邊界條件和初始條件； a) 粒子坐標的生成通過輸入相應外形參數實現，根據外形參數的幾何特征及給定的粒 子間距離，便可求得該幾何形狀內所有粒子的坐標； b) 粒子屬性設置即是將邊界粒子屬性設置為固體，其他粒子屬性設置為流體的過程； c) 邊界條件確定了邊界固體粒子對內部流體粒子的作用方式，包括兩方面內容：固體 壁面邊界條件、溫度邊界條件； 固體壁面邊界條件即是計算域邊界部分粒子的屬性為固體粒子，邊界粒子位置不變， 并對靠近的流體粒子產生一定作用力，作用力施加方式包括邊界力法和鏡像粒子法兩種方 式；其中邊界力法通過施加排斥力實現固體邊界；鏡像粒子法通過施加排斥動量實現固體 邊界； 溫度邊界條件即是當存在具有一定溫度特性的邊界時，該部分邊界粒子需要設置溫度 邊界條件，以使其溫度特性影響其附近流體粒子，溫度邊界條件包括恒溫邊界條件和絕熱 邊界條件兩種；恒溫邊界條件情況下，該部分固體粒子溫度恒定，持續與周圍粒子進行熱交 換；絕熱邊界條件情況下，該部分固體粒子與其他粒子不發生任何熱交換； d) 初始條件為計算域內粒子的初始物理參數，包括：速度、壓力、溫度、參考溫度、密 度、運動粘度、重力、瑞利數、普朗特數。
【文檔編號】G06F17/50GK104143027SQ201410377151
【公開日】2014年11月12日 申請日期:2014年8月1日 優先權日:2014年8月1日
【發明者】雷娟棉, 楊浩, 黃燦, 王鎖柱 申請人:北京理工大學