一種經(jīng)gpu加速的浸沒邊界-格子玻爾茲曼流固耦合模擬方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明屬于計(jì)算流體動(dòng)力學(xué)和計(jì)算固體力學(xué)的交叉領(lǐng)域，特別設(shè)及高性能計(jì)算在 上述領(lǐng)域的應(yīng)用。
【背景技術(shù)】
[0002] 流體和可變形的結(jié)構(gòu)之間的相互作用在自然界普遍存在，特別在生物組織和器官 中更是隨處可見，典型的例子包括昆蟲翅膀、魚罐、人類屯、臟瓣膜和聲帶。雖然上述結(jié)構(gòu)屬 于解創(chuàng)學(xué)或生理學(xué)領(lǐng)域的研究對(duì)象，但是他們?cè)?維空間的運(yùn)動(dòng)學(xué)特征和隨時(shí)變化的空間 形狀對(duì)生物體完成指定的功能至關(guān)重要。因?yàn)榻Y(jié)構(gòu)所經(jīng)歷的劇烈變形和周圍流場(chǎng)的內(nèi)在復(fù) 雜性，所W迄今為止正確模擬=維條件下的流固禪合問題仍然是學(xué)術(shù)界的一大挑戰(zhàn)之一。
[0003] 傳統(tǒng)的用來模擬流固禪合問題的方法通?；谶m體網(wǎng)格，典型的例子包括任意拉 格朗日一歐拉法和有限元法。上述方法的劣勢(shì)在于底層流體網(wǎng)格需要根據(jù)結(jié)構(gòu)位置來進(jìn)行 調(diào)整，而且為了防止網(wǎng)格的劇烈變形需要進(jìn)行特殊處理?？偟膩碚f網(wǎng)格生成是代價(jià)高昂的 過程，特別是當(dāng)結(jié)構(gòu)形狀復(fù)雜時(shí)為了更新網(wǎng)格拓?fù)浣Y(jié)構(gòu)需要大量耗時(shí)的計(jì)算，使得將傳統(tǒng) 方法推向?qū)嶋H應(yīng)用存在困難。
[0004] 因此，非常有必要針對(duì)效率問題開發(fā)一種新的流固禪合模擬方法，能夠在不W降 低計(jì)算精度為代價(jià)的前提下，縮短計(jì)算時(shí)間，為實(shí)際工程提供可靠的指導(dǎo)建議。
[0005] 浸沒邊界法作為一種?？谀M流固禪合問題的數(shù)值算法自提出W來已經(jīng)得到了 廣泛和成功的應(yīng)用。在此法中，流體和結(jié)構(gòu)體分別采用歐拉法和拉格朗日法描述，且分別用 來離散流體和結(jié)構(gòu)體的正交笛卡爾網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格沒有必要重合在一起，而該正是浸 沒邊界法相比于傳統(tǒng)方法的一大優(yōu)勢(shì)，解決了每一時(shí)間步都需要根據(jù)結(jié)構(gòu)體位置和形態(tài)重 新構(gòu)造網(wǎng)格的難點(diǎn)。流體和結(jié)構(gòu)體的相互作用是通過狄拉克函數(shù)實(shí)現(xiàn)的，利用此函數(shù)能夠 將定義在結(jié)構(gòu)體網(wǎng)格節(jié)點(diǎn)上的作用力擴(kuò)散至周邊流體網(wǎng)格點(diǎn)，并且也能夠?qū)⒘黧w速度插值 到結(jié)構(gòu)體網(wǎng)格點(diǎn)上，從而更新結(jié)構(gòu)體的位置。除了上述步驟外，無需額外的附加處理即可非 常簡(jiǎn)潔地完成流固禪合的模擬。
[0006] 在浸沒邊界法中，最基本的一點(diǎn)在于如何求解粘性不可壓縮NS方程。該方程的求 解方法非常豐富，比如說快速傅里葉變換、投影方法或者是有限體積法。此外新近發(fā)展起來 的非常受歡迎的格子玻爾茲曼方法也是很好的選擇。與傳統(tǒng)方法求解的NS方法不一樣，采 用介觀視角的格子玻爾茲曼方法求解的是分布函數(shù)所滿足的玻爾茲曼方程。在過去的十年 間，格子玻爾茲曼方法已經(jīng)被確認(rèn)為一種非常強(qiáng)大的計(jì)算流體動(dòng)力學(xué)工具，并且在其分支 領(lǐng)域都已得到了成功應(yīng)用。和傳統(tǒng)方法比起來，格子玻爾茲曼方法的易于并行W及額外物 理效應(yīng)的方便引入等優(yōu)點(diǎn)得到越來越多學(xué)者的關(guān)注。
[0007] 自從2003年W來，半導(dǎo)體工業(yè)的微處理器設(shè)計(jì)開始朝兩個(gè)主要方向發(fā)展。首先多 核路徑依靠強(qiáng)化串行代碼的執(zhí)行效率和多核技術(shù)提高性能。多核開始于兩核模型，并且隨 著每一次更新?lián)Q代核的數(shù)量加倍，比如Intel的Core i7處理器有四核。與此相反，眾核路 徑則將關(guān)注的重點(diǎn)放在了提高并行代碼執(zhí)行的吞吐量上面，和中央處理器比較，眾核處理 器將包含數(shù)量更多但是小得多的核屯、，比如NVIDIA公司的Ge化rceGTX280圖形處理器包 含240個(gè)核。作為比較，2009年圖形處理器和中央處理器的峰值運(yùn)算性能的比值是10:1， 2014年NVIDIA公司的Geforce780Ti所提供的峰值吞吐量為5500GFL0PS，而Intel公司 的Ivy化idge架構(gòu)的數(shù)據(jù)為700G化OPS。如此大的運(yùn)算性能差距已經(jīng)讓越來越多的科研人 員將它們程序中的計(jì)算密集部分轉(zhuǎn)移到圖形處理單元來執(zhí)行。
[000引因?yàn)楸热缦窕痚n化和Direct3D技術(shù)需要直接針對(duì)巧片編程，而直到2006年程 序員必須通過應(yīng)用程序編程接口來訪問其核屯、，所W圖形處理器的使用非常困難。該些應(yīng) 用程序編程接口限制了能夠?yàn)榍善帉懙膽?yīng)用程序的類型，所W只有小部分程序員掌握了 利用圖形處理器運(yùn)行特定程序的必要技術(shù)，進(jìn)而限制了該項(xiàng)技術(shù)的廣泛傳播。一切隨著 NVIDIA公司在2007年推出CUDA技術(shù)而改變。因?yàn)橐?guī)避了傳統(tǒng)的應(yīng)用程序編程接口，通過 CUDA，現(xiàn)在程序員們能夠用他們所熟悉的類C/C++語言來直接針對(duì)巧片編程，從而使得并 行編程技術(shù)即便對(duì)于初學(xué)者來說也非常容易上手，大大地推進(jìn)了該項(xiàng)技術(shù)在科學(xué)領(lǐng)域的應(yīng) 用。
[0009] 因?yàn)楦褡硬柶澛椒ㄊ嵌x在正交網(wǎng)格上的顯格式，并且只需要交換相鄰節(jié)點(diǎn) 的信息，所W它非常適合于在圖形處理器上并行化。發(fā)展至今，發(fā)表的關(guān)于格子玻爾茲曼方 法并行化的文章已非常豐富，并且相比于中央處理器，運(yùn)行速度一般能提高兩個(gè)數(shù)量級(jí)。雖 然格子玻爾茲曼方法的并行化已經(jīng)日臻完善，但是將其與浸沒邊界法結(jié)合起來模擬流固禪 合問題的研究仍然沒有起步，更沒有全面透徹的分析浸沒邊界法潛在數(shù)據(jù)并行性的文章發(fā) 表。

【發(fā)明內(nèi)容】

[0010] 針對(duì)現(xiàn)有技術(shù)存在的不足，本發(fā)明提出了一種基于GPU的高效流固禪合模擬方 法。
[0011] 為了解決上述技術(shù)問題，本發(fā)明采用如下的技術(shù)方案：
[0012] 本發(fā)明的具體操作步驟如下：
[0013] 1)在專業(yè)前處理軟件中建立待測(cè)結(jié)構(gòu)體及其內(nèi)部流體的浸沒邊界體型，并劃分邊 界網(wǎng)格W供程序讀?。?br>[0014] 2)程序讀取邊界體型數(shù)據(jù)，設(shè)置計(jì)算區(qū)域的邊界條件和程序模型參數(shù)；
[0015] 3)根據(jù)邊界當(dāng)前時(shí)刻構(gòu)型和參考構(gòu)型利用有限元法計(jì)算邊界施加在流體上的作 用力；
[0016] 4)將邊界的作用力經(jīng)由狄拉克函數(shù)擴(kuò)散至周邊流體；
[0017] 5)利用格子玻爾茲曼方法求解攜帶外力項(xiàng)的納維-斯托克斯方程；
[0018] 6)同樣經(jīng)由狄拉克函數(shù)將流體速度插值到邊界上得到邊界的移動(dòng)速度，進(jìn)而更新 邊界位置；
[0019] 7)重復(fù)步驟3-6,直至達(dá)到計(jì)算終點(diǎn)。
[0020] 所述專業(yè)前處理軟件為Gambit或ICEM或ProE。
[0021] 本發(fā)明中，詳細(xì)過程如下：
[0022] 一、采用浸沒邊界法處理流體和結(jié)構(gòu)體之間的相互作用。
[0023] 具體來說包含W下幾方面的特征：
[0024] 1)二維條件下，采用一維線網(wǎng)格離散結(jié)構(gòu)體。S維條件下，采用非結(jié)構(gòu)化S角形網(wǎng) 格離散結(jié)構(gòu)體。
[0025] 2)二維條件下，同時(shí)考慮結(jié)構(gòu)體抵抗拉伸和彎曲的能力，定義在節(jié)點(diǎn)上的抗拉和 抗彎力密度分別按下式計(jì)算：
【主權(quán)項(xiàng)】
1. 一種經(jīng)GPU加速的浸沒邊界-格子玻爾茲曼流固耦合模擬方法，其特征在于：步驟 如下： (1) 在專業(yè)前處理軟件中建立待測(cè)結(jié)構(gòu)體及其內(nèi)部流體的浸沒邊界體型，并劃分邊界 網(wǎng)格以供程序讀取； (2) 程序讀取邊界體型數(shù)據(jù)，設(shè)置計(jì)算區(qū)域的邊界條件和程序模型參數(shù)； (3) 根據(jù)邊界當(dāng)前時(shí)刻構(gòu)型和參考構(gòu)型利用有限元法計(jì)算邊界施加在流體上的作用 力； (4) 將邊界的作用力經(jīng)由狄拉克函數(shù)擴(kuò)散至周邊流體； (5) 利用格子玻爾茲曼方法求解攜帶外力項(xiàng)的納維-斯托克斯方程； (6) 同樣經(jīng)由狄拉克函數(shù)將流體速度插值到邊界上得到邊界的移動(dòng)速度，進(jìn)而更新邊 界位置； (7) 重復(fù)步驟3-6,直至達(dá)到計(jì)算終點(diǎn)。
2. 如權(quán)利要求1所述的經(jīng)GPU加速的浸沒邊界-格子玻爾茲曼流固耦合模擬方法，其 特征在于：所述專業(yè)前處理軟件為Gambit或ICEM或ProE。
3. 如權(quán)利要求1所述的經(jīng)GPU加速的浸沒邊界-格子玻爾茲曼流固耦合模擬方法，其 特征在于：所述步驟3中： (1) 當(dāng)浸沒邊界體型為二維時(shí)，采用一維線網(wǎng)格離散結(jié)構(gòu)體，同時(shí)考慮結(jié)構(gòu)體抵抗拉伸 和彎曲的能力，定義在節(jié)點(diǎn)上的抗拉和抗彎力密度分別按下式計(jì)算：
上式中的（Fs)1和（Fb)別表示定義在編號(hào)為1節(jié)點(diǎn)上的抗拉和抗彎作用力線密度； KJP Kb分別為抗拉系數(shù)和彎曲剛度；Λ a 1為編號(hào)1的節(jié)點(diǎn)的原始弧長(zhǎng)；X 和Xm分別 為編號(hào)m+1、m-Ι和m的節(jié)點(diǎn)坐標(biāo)；δ ml為克羅內(nèi)克符號(hào)，當(dāng)m = 1時(shí)其取值為1，其余情況 下取值為〇 ;nf為用來離散結(jié)構(gòu)體的總節(jié)點(diǎn)數(shù)； (2) 當(dāng)浸沒邊界體型為三維時(shí)，采用非結(jié)構(gòu)化三角形網(wǎng)格離散結(jié)構(gòu)體， 結(jié)構(gòu)體在編號(hào)為m節(jié)點(diǎn)上的作用力密度計(jì)算公式如下：
其中：W為結(jié)構(gòu)的總能量，3為偏微分符號(hào)，Xm為節(jié)點(diǎn)m的坐標(biāo)；結(jié)構(gòu)總能量的計(jì)算式 為： W = / wsdA 式中的dA為面元面積，而ws為結(jié)構(gòu)的拉能密度，其計(jì)算式為：
，其中E為彈性模量，IjP 12分別為第一和第二應(yīng)變不變量，可以 證明應(yīng)變不變量能夠按下式計(jì)算：
其中λ Jp λ 2是法平面內(nèi)沿一對(duì)垂直方向的拉伸比，而a、b、c的計(jì)算式為：
上式中幾何參數(shù)的定義如下：在結(jié)構(gòu)體參考構(gòu)型和變形后構(gòu)型上選取兩塊相應(yīng)的三角 形面元，分別在上述三角形中任意選取兩條邊，在參考構(gòu)型上其邊長(zhǎng)分別為^和Γ C1，而在 變形后構(gòu)型上其邊長(zhǎng)分別為1和Γ，灼和爐分別為上述兩條邊在參考構(gòu)型和變形后構(gòu)型 上的夾角。
4. 如權(quán)利要求1所述的經(jīng)GPU加速的浸沒邊界-格子玻爾茲曼流固耦合模擬方法，其 特征在于：所述步驟4和步驟6中的作用力擴(kuò)散方程為：
上式中的f(x，t)為t時(shí)刻點(diǎn)X處流體受到的作用力密度，F(xiàn)mS編號(hào)為m節(jié)點(diǎn)上的作用 力密度，分別針對(duì)二維和三維條件，△ Sk分別是定義在編號(hào)為k的結(jié)構(gòu)體節(jié)點(diǎn)上的弧長(zhǎng)和面 積，而S (x-Xm)即為狄拉克函數(shù)，一般采取如下具體形式：
為了能夠得到結(jié)構(gòu)體形狀的變化過程，我們需要根據(jù)結(jié)構(gòu)體網(wǎng)格點(diǎn)的速度更新結(jié)構(gòu)體 的構(gòu)型，在已知流場(chǎng)速度分布條件下，結(jié)構(gòu)體速度同樣將利用狄拉克函數(shù)通過從流場(chǎng)插值 得到：
上式中的Uni分別為結(jié)構(gòu)體網(wǎng)格點(diǎn)的坐標(biāo)和速度，u(x，t)是t時(shí)刻點(diǎn)X處流體的 速度，Ah是網(wǎng)格步長(zhǎng)。
5. 如權(quán)利要求1所述的經(jīng)GPU加速的浸沒邊界-格子玻爾茲曼流固耦合模擬方法，其 特征在于：所述步驟5中，格子玻爾茲曼方程為：
式中，ga (x，t)表示t時(shí)刻X處沿離散速度方向ξ α的分布函數(shù)，而g α (χ+ξ α St，t+ St)則為t+St時(shí)刻χ+ξ aSt處沿離散速度方向ξ α的分布函數(shù)，g】q(x，0是t時(shí)刻x處 沿離散速度方向ξ a的平衡態(tài)分布函數(shù)，τ是無量綱的松弛時(shí)間，與分子間的粘度系數(shù)有 關(guān)，St為時(shí)間步長(zhǎng)。
【專利摘要】本發(fā)明提供一種經(jīng)GPU加速的浸沒邊界-格子玻爾茲曼流固耦合模擬方法，步驟如下：在專業(yè)前處理軟件中建立待測(cè)結(jié)構(gòu)體及其內(nèi)部流體的浸沒邊界體型，并劃分邊界網(wǎng)格以供程序讀取；程序讀取邊界體型數(shù)據(jù)，設(shè)置所需測(cè)量點(diǎn)的邊界條件和程序模型參數(shù)；根據(jù)邊界當(dāng)前時(shí)刻構(gòu)型和參考構(gòu)型利用有限元法計(jì)算邊界施加在流體上的作用力；將邊界的作用力經(jīng)由狄拉克函數(shù)擴(kuò)散至周邊流體；利用格子玻爾茲曼方法求解攜帶外力項(xiàng)的納維-斯托克斯方程；同樣經(jīng)由狄拉克函數(shù)將流體速度插值到邊界上得到邊界的移動(dòng)速度，進(jìn)而更新邊界位置；重復(fù)步驟3-6，直至達(dá)到計(jì)算終點(diǎn)。本發(fā)明能夠處理任意形狀的結(jié)構(gòu)體，計(jì)算過程中不需要重構(gòu)底層流體網(wǎng)格，計(jì)算效率高。
【IPC分類】G06F17-50
【公開號(hào)】CN104866695
【申請(qǐng)?zhí)枴緾N201510355061
【發(fā)明人】吳家陽, 程永光, 張春澤, 刁偉
【申請(qǐng)人】武漢大學(xué)
【公開日】2015年8月26日
【申請(qǐng)日】2015年6月24日