基于GPU的原油液體仿真效率優化研究

劉賢梅閆沖

摘要：針對有桿抽油系統的采油過程，原油液體在油管內受力而不斷運動。為了將此過程形象、逼真的展現出來以便于相關工作人員觀察和分析，需建立一個高精度、高效率、易使用的可視化環境。本文在VC++6.0環境下，運用粒子系統技術模擬原油液體的運動過程，并采用通過GPU加速和非GPU加速對比的方式優化了系統效率。通過數據對比表明：隨著粒子數的增加，通過采用GPU加速運算方式的系統效率大大提高。

關鍵詞：可視化；粒子系統；GPU

中圖分類號：TP393文獻標識碼：A

1引言

隨著油田鉆井技術的發展，有桿采油系統中的可視化問題已經成為了油田行業的一個重大問題，特別是在采油作業的過程中，需要時刻定性、定量的了解井筒中設備的運動狀況以及原油液體狀態變化，方便相關技術人員觀察與分析，以確保采油過程順利進行。

GPU的并行處理技術已經在各個領域廣泛應用。1999年，文獻1[1]在求解輻射度方程時利用顏色來標示物體序號，利用深度緩沖來輔助計算物體可見性。2002年，文獻22]等利用圖形硬件進行實時的運動規劃，即用顏色來標記多邊形。2004年，文獻3[3]等利用圖形硬件的多邊形光柵化過程和深度緩沖檢測來計算二維和三維的通用Voronoi圖。2006年，文獻4[4]等采用均勻網格來計算場景中的紋理，并建立三角形鏈表紋理，整個計算過程都放在GPU上進行，提高了效率。自此，GPU的并行計算技術在碰撞檢測、數學運算、像素模擬等方面得到了廣泛的應用。

本文充分考慮圖形處理技術的并行性，以VC++6.0為開發平臺，采用GPU加速的方法，運用粒子系統技術模擬了原油液體在油管內的運動狀態，建立一個高精度、高效率、易使用的可視化環境。本系統對液體仿真過程進行優化，是GPU技術在液體仿真和油田領域的一次新的突破。最終通過數據對比表明：隨著粒子數的增加，通過采用GPU加速運算的方法使系統效率大大提高。

2粒子系統

2.1粒子系統簡介

迄今為止，粒子系統作為一種有效的圖形生成算法通常被用來模擬不規則的模糊物體，其基本思想就是將無數微小粒子作為最基本的元素來表示不規則物體[5-6]。根據粒子系統的應用，大致分為四種類型：

1）隨機粒子系統。在隨機粒子系統中，粒子的屬性隨機變化，例如煙花的綻放過程、流水過程等。

2）流體粒子系統。流體粒子系統中的粒子運動線路和運動軌跡受流體力學的影響。例如，工廠里的煙霧擴散情況、瀑布、以及高溫高壓下的巖漿流動情況。

3）方向粒子系統。方向粒子系統中的粒子在活動過程中相互作用、相互影響。在方向粒子系統中，所有粒子的屬性具有關聯性，主要用于綿軟不規則的物體，例如紡織物等。

4）結構化粒子系統。結構化粒子系統是由具有一定組織結構的較小粒子組成的系統。其基本單元在一定的時間內變化量較小，例如微風下的花草樹木。

2.2粒子的基本性質

1）粒子的物質性。在粒子系統中，在任何時間和空間上，粒子與粒子、粒子與其他物體之間都不會發生重疊，具有絕對的物質性。

計算技術與自動化2015年12月

第34卷第4期劉賢梅等：基于GPU的原油液體仿真效率優化研究

2）粒子的基本屬性。粒子系統中的每一個粒子都有其基本屬性，包括運動屬性和靜態屬性兩種。靜態屬性是指屬性物體的屬性不會隨時間的變化而變化，例如大小、顏色、生命周期、質量等。運動屬性是指物體的屬性會隨時間的變化而產生變化，這個變化可以是線性的，例如速度、位移等。

3）粒子的生命機制。每個粒子都有其生命周期，一般粒子要經歷產生、活動、消亡基本過程。粒子的生命閾值到達上限是會被自動刪除，在此過程中會有新的粒子產生。

4）粒子的活動機制。粒子在其生命周期內會按照預先設定好的運動規律在場景中進行運動，直到生命周期結束。

5）粒子的繪制。根據被模擬物體的基本特性設定粒子的靜態屬性，再分析被模擬物體的運動特性設定粒子的動態屬性。系統中的粒子按照產生、活動、消亡、更新的順序進行全生命周期的活動。

2.3粒子系統模型

粒子系統是按照時間順序進行全生命周期活動的系統[7]。粒子系統中的每一個粒子都要經歷產生、活動、消亡這三個過程。在進行每一次粒子繪制的時候，整個粒子系統要維持一個穩定的循環體系，對于活動中的粒子，我們要根據被模擬對象的屬性來更新粒子的屬性，其運行基本步驟主要包括五部分：

1）在指定區域隨機生成粒子，并賦予粒子隨機屬性；

2）根據被模擬物體的基本特性來歸納總結粒子的基本活動規律，設置好粒子的靜態屬性或者動態屬性，例如初始位置、顏色、大小、生命周期等。；

3）根據粒子生命周期，將粒子的當前存在時間與生命上線進行對比，判斷粒子的生命狀態，刪除生命值為零的粒子并根據需要更新粒子；

4）根據粒子的動態屬性或靜態屬性，模擬粒子的位移變換；

5）將活動粒子進行渲染，最終描繪出被模擬物體。

3可編程GPU技術

3.1GPU工作流程

GPU處理復雜和大量圖形具有很高的效率[8]，主要是由于：

1）減少了GPU與CPU之間的數據通信

當整個應用針對圖形生成的時候，不再需要在CPU與GPU之間進行多次數據交換，這樣就可以將CPU解放出來做其他的任務。這些優勢使得GPU比CPU更適用于流處理計算，因此GPU也被認為是一個SIMD的并行機或者流處理器，可以用于處理大規模數據集，使應用得到加速[9]。

2）向量運算架構

GPU本質上是一個向量計算模型架構，計算單元豐富，它的這種架構使它在處理大規模向量運算時體現出較高的性能[10]，GPU運算過程如圖1所示。

3.2GPU數據存儲與訪問

在GPU運算中，大量數據需要存儲在紋理中，通過紋理坐標來訪問紋理中的數據[11]。紋理相當于一個存儲顏色值的二維數組，紋理中每個元素都存儲著一個四維向量的顏色值。在進行GPU通用計算時，將大量的數據存儲到紋理中[12]。紋理有多種像素格式，可以根據每種數據的精度和值域來選擇合適的像素格式進行存儲。如果存儲的數據是布爾型或者是精度很低的整型，可以用低精度的8bitsARGB2222格式來存儲；如果是高精度的浮點型數據，可以選用相應的浮點型紋理格式來存儲[13]。在GPU通用計算中，浮點像素格式是比較通用的像素格式，其每個顏色通道都可以存儲一個32位的浮點數據[14-15]。

在紋理中通過紋理坐標對數據進行索引。紋理坐標是一個二維的實數向量，一般紋理坐標的值域在O到1之間，索引屬性值需要精確定位到紋理中的像素。例如在一個m×n尺寸的紋理中索引第x列第y行的像素則紋理坐標的計算方法是：

U=x/m+1/（2×m）；V=y/n+1/（2×n）。

其中U和V分別為UV兩個方向的紋理坐標值。

4系統實現與分析

4.1GPU下油管內液體粒子的數據組織

實驗描述：模擬在上沖程過程中，原油液體在泵筒中由下往上的運動狀態，通過機組實驗數據對比在粒子數變化幅度相同的情況下，CPU和GPU的效率變化。

在上沖程過程中，抽油桿和柱塞上行，原油液體進入泵筒內部，此時液柱受到重力G（豎直向下）、推力ft（豎直向上）、摩擦阻力fm（沿著泵筒向下）。在泵筒內液柱向上運動，根據原油粒子特點，需要將以下粒子屬性儲存到紋理中：

Speed：用三維向量表示液體粒子的流動速度；

Live：粒子的生命閾值，存儲為布爾型的標量，True表示粒子處于活動中，False表示粒子消亡；

Position：用三維向量存儲粒子在世界坐標系下的位置。

由于液體粒子數量很多、大小不一，因此存放粒子時應采取紋理多樣化的方式，根據不同精度選擇合適的像素格式進行有效存儲。例如：當存儲低精度數據時，采用低精度的8 bits ARGB 2222格式，反之當存儲高精度數據時，采用高精度的32 bits ARGB 8888格式。在GPU的通用計算中，為了滿足計算的精度要求，通常浮點像素格式的每一個顏色通道都將存儲一個8位或者32位的浮點數據。在存儲數據的過程中既要考慮效率問題又要考慮存儲空間的合理利用問題，通常的像素格式都是具有RGBA四個顏色通道，粒子的速度和位置都是三維向量，這樣把粒子的速度和位置存放到紋理像素中的時候就會剩余一個通道，所以將粒子的Live屬性儲存到最后一個通道里，這樣就充分合理的利用了存儲空間。

格式選取完畢后開始創建存儲粒子數據的紋理，確保每一個粒子對應一個像素點，這樣像素點的最大值等于粒子的最大數目。每一個粒子的屬性值存儲在通過紋理中并通過各自的紋理坐標進行檢索。

4.2液體粒子屬性計算

在液體粒子產生、活動、消亡的過程中，把握并控制好粒子的屬性，準確實時計算出粒子的速度、生命值、位置是繪制液體粒子的核心。首先編寫像素著色器程序，通過GPU的通用計算將粒子屬性的計算問題具體化，計算出每一個液體粒子的屬性值并存儲到紋理中以備使用。

利用GPU通用計算來計算每一個粒子屬性的時候，粒子要經歷產生、活動。消亡三個過程。由于在GPU下不能動態為粒子分配存儲空間，所以每一個活動中的粒子的生命值由最大值降為零時，將此粒子的生命值、位置、速度進行初始化，作為新粒子重新加入到系統中。粒子系統開始運行后，粒子的數量隨著時間的推移不斷增加，粒子的屬性以二維索引的形式存儲在紋理中，在進行GPU通用計算時需要將二維索引轉化為一維索引。對于二維索引來說，需要依靠紋理坐標來將其轉換為一維索引，轉換公式為：M=（U+V×H）×W，其中M是二維索引轉換為一維索引的序列號，U和V為紋理坐標的兩個分量，W和H為兩個紋理坐標分量的大小。液體粒子屬性值計算步驟如下：

1）液體在受力中的體積和數量的計算

原油液體在油管內流動的過程中，粒子系統在每一幀都要保持產生新的粒子以維持液體的流動，這樣就要對新生成的粒子的速度、位置、生命值進行初始化。在上沖程過程中，油管內的液柱隨著柱塞的上升而升高，所有液體粒子由下部初始位置升到最高點，初始點的所有粒子位置要添加一個隨機值：P=P0+Prant，其中P是粒子的初始位置，P0是油管內液柱上升的初始位置，Prand是一個三維隨機向量，此向量決定液柱的粗細。如圖2為上沖程過程中原油液體進入泵筒內的示意圖。此時泵筒內原油液體受重力G（豎直向下）、推力ft（豎直向上）、摩擦阻力fm（沿著泵筒向下）、泵筒內壓力fp。如圖2中x表示柱塞在時刻t的位移，表示柱塞位移為x時的液面高度，Lw是尾管長度。設任意時刻活塞t運動到x的位置，進入泵筒內原油液體體積為Vt。根據受力分析，將整個泵筒內空間分為兩部分，Va為原油液體未充滿泵筒前真空部分的體積，Vb為下半部分體積，其中Va=（x-y）A，Vt=Va+Vb，Vb=（Rs（x）-Rs（0））（1-Lw）（Ls+y）ABo，其中Rs是溶解油氣比、Bo為原油體積系數與水的體積系數、A為壓力系數。假設單位粒子體積為1，則有：1×Vnum=Vt，其中Vnum為粒子隨體積變化的數量。

2）活動粒子屬性值的計算

上沖程過程中，液體粒子受力向上運動，其位置和速度在不斷地變化著。下面分析粒子的受力狀況，根據里的平衡原理得出：fm+fp+G=ft，vm=vm-1+a×Δt，其中Vm是粒子的速度，Vm-1是上一幀粒子的速度，a是粒子的加速度，Δt是一幀持續的時間。在Y方向上，粒子的加速度由受力變化確定，在X和Z方向上粒子的速度分量不變。Pm=Pm-1+（a×Δt）×Δt，其中Pm是粒子的位置，Pm-1是粒子在上一幀的位置，Δt是一幀持續的時間。

3）消亡粒子屬性值的計算

當液體粒子在油管內達到最高點時，粒子的生命值降為零，此時的粒子屬性不在發生變化。在此粒子重新生成并再次加入粒子系統之前，既不需要計算其屬性值，也不用被渲染。

4.3液體粒子繪制

計算完液體粒子的位置和速度，下一步需要讀取粒子的屬性值并進行繪制。運用紋理獲取技術（Vertex Texture Fetch）來獲取粒子的屬性值，用tex2Dlod函數在頂點著色器讀取紋理像素中的粒子屬性值進行繪制，這樣就不需要通過CPU傳輸，大大降低了CPU與GPU之間的通信，提高了渲染速度。下面進行粒子的繪制步驟：

1）建立頂點緩沖區

選擇圓形作為繪制液體粒子的基本元素，組成每個粒子的最基本單位設定為四邊形。首先建立頂點緩沖區，頂點的個數為四倍的粒子數，然后建立頂點索引，將每四個頂點組合在一起組成四邊形，粒子以四邊形的方式被繪制出來。創建出來的頂點緩沖區需要預先設計其格式。首先給液體粒子貼上紋理（此處選擇普通水滴紋理），利用Texcoord0（紋理通道0的紋理坐標）給粒子的頂點的紋理坐標分別賦值為（0，1）（1，0）（0，1）（1，1），這樣就可以通過紋理坐標0將水滴的紋理映射在粒子上，并判斷出水滴紋理在粒子上的位置。其次要對粒子的屬性值進行索引，將粒子的紋理屬性存儲在Texcoord1（紋理通道1的紋理坐標）中，這樣就可以索引到粒子的屬性值了。

2）頂點繪制

首先讀取粒子的紋理屬性并獲取粒子在世界坐標系下的位置，由于粒子必須要正對著屏幕，所以要將粒子變換到相機坐標下，利用紋理坐標將粒子的四個頂點分別繪制在粒子的四個角上。由于GPU會將所有頂點繪制出來，即無論是活動中的頂點還是已經消亡的頂點，GPU都會無選擇性的將其全部繪制出來，因此把消亡的粒子的頂點位置全部設定到一起，這樣粒子就不會顯示出來。頂點位置的計算公式如下：

Pt=（P×Mv+（U0，V0，0）×Size×Live）×Mp，其中Pt是最終計算的頂點位置，P是屬性紋理中的粒子位置，U0和V0分別是頂點紋理坐標0的兩個分量，Size是粒子的大小，Live是粒子的生命值，Mv和Mp分別是相機矩陣和投影矩陣。當粒子消亡時，粒子的Live值為零，粒子不需要被繪制出來；反之當粒子活動時，粒子的Live值為1，粒子頂點的位置會被正確的計算出來。

5效果與結論

本實驗硬件環境為：Intel2.8GHz處理器，4G內存，NVIDIAGeForceGT640M顯卡，操作系統是64位Window7系統，開發工具為VC++6.0和Directx9.0。圖4.1為原油液體的流動效果圖。表1為獨立環境下CPU模擬的原油液體和GPU模擬的原油液體幀速率的比較。如圖4.2為表一數據的對比圖。

結論：針對有桿抽油系統采油過程中原油液體在油管內的受力運動，本文設計并建立了一個能滿足高精確度、高效率、易使用的可視化環境。從表中我們可以看出，當原油液體粒子數較少時，CPU實現的速度比GPU還要快，但隨著粒子數的增大，CPU實現的幀速率快速下降，而GPU實現的幀速率下降速度緩慢，當粒子數達到300萬時，GPU仍然能夠實時的繪制，而CPU實現每秒只有6.8幀，這說明GPU對大量粒子的處理能力遠遠高于CPU。

參考文獻

[1]HopfMatthias，ErtlThomas.Accelerating 3D convolutionusinggraphicshardware.[J] Proceedings of IEEE Visualization.1999，28（05）：56-70，62.

[2]FRANCISCO J Seron，Rafael Rodrigues，Eva Cerezo，AlfredoPina.AddingSupport for Highlevel Skeletal Animation[J]. IEEE Transactions onvisualization and computer graphics.2002，28（05）：120-134，122.

[3]Peter Kipfer，Mark Segal，Rüdiger Westermann.UberFlow：a GPUbased particle engine[J]. Proceedings of the ACM SIGGRAPH/EUROGRAPHICS Conference on Graphics HardwareHWWS04.2004， 28（05）：128-132，130.

[4]ARON E Lefohn，Shubhabrata Sengupta，Joe Kniss，et al.Glift：Generic，efficient，randomaccess GPU data structures[J]. ACM Transactions on Graphics.2006， 28（05）：120-134，123.

[5]田蘇昕.虛擬現實中大規模粒子系統的研究[J].沈陽航空工業學院學報，2009，30（3）：149-150，160.

[6]徐陽東.基于粒子系統的不規則景物建模研究[J].山東師范大學學報，2009， 28（5）：128-132，138.

[7]白玉.基于粒子系統的煤礦火災模擬系統[J].吉林大學學報，2010，8（2）：47-49.

[8]王濤. 基于GPU的程序分析與并行化研究[J].解放軍信息工程大學學報，2010，8（2）：47-51.

[9]胡杰. CPU-GPU異構平臺計算模型的研究與應用[J].大連理工大學學報，2011，29（9）：44-48.

[10]郭忠明. 基于CUDA的并行圖像處理性能優化[J].大連理工大學學報，2012，99（4）：49-58.

[11]方旭東. 面向大規?？茖W計算的CPUGPU異構并行技術研究[J].國防科學技術大學學報，2009，22（9）：165-169.

[12]BUCK I，FOLEY T，HORN D，SUGERMAN J，HANRAHAN P.Brook for GPUs： Stream Computing on Graphics Hardware[J].Proceedings of SIGGRAPH，2004，28（4）：55-59.

[13]Jizhu Lu，et.al.HMMerCell： High Performance Protein Profile Searching on the Cell/B.E.Processor[J].Proceeding of IEEE International Symposium on Performance Analysis of Systems and SoftwareISPASS，2008，33（7）：94-98.

[14]S.García，A.Fernández，J.Luengo et al.A study of statistical techniques and performance measures for geneticsbased machine learning： accuracy and interpretability[J].Soft Computing，2009，83（7）：104-108.

[15]Zhe Wei， Yixiong Feng， Jianrong Tan et al.. Multi-objective performance optimal design of largescale injection molding machine[J]. The International Journal of Advanced Manufacturing Technology 2009， 30（3）：166-170，160.