龍文佳,李躍新

(湖北大學(xué) 計(jì)算機(jī)與信息工程系，武漢 430000)

?

三維多物質(zhì)歐拉界面處理的并行算法研究

龍文佳,李躍新

(湖北大學(xué) 計(jì)算機(jī)與信息工程系，武漢 430000)

針對現(xiàn)有的算法無法對爆炸沖擊進(jìn)行有效的研究這一問題，在改進(jìn)了Youngs界面技術(shù)的基礎(chǔ)上設(shè)計(jì)了三維多物質(zhì)Euler界面處理的并行算法;然后采用MPI標(biāo)準(zhǔn)進(jìn)行了算法的程序設(shè)計(jì)并對程序進(jìn)行了并行性能測試和不同分區(qū)的算例測試，測試結(jié)果表明并行算法對加速比有較大的提高，而不同分區(qū)對并行算法的計(jì)算結(jié)果并沒有影響;最后應(yīng)用所編寫的程序?qū)罩斜ê途勰苌淞鬟M(jìn)行了仿真，結(jié)果表明程序的模擬結(jié)果與前人的研究數(shù)據(jù)是符合的，說明了程序的有效性。

爆炸與沖擊;數(shù)值模擬;歐拉方法;并行計(jì)算;Youngs界面技術(shù)

0 引言

爆炸沖擊涉及大量的非線性瞬態(tài)動(dòng)力學(xué)問題，屬于復(fù)雜大變形的流體彈塑性力學(xué)范疇，而Euler算法[1-2]對于解決該類問題具有很大的優(yōu)勢，但其很難描述各介質(zhì)的界面，所以界面處理的精度問題一直是Euler算法的核心問題。Youngs界面技術(shù)能夠?qū)煞N物質(zhì)界面進(jìn)行精確的處理，但在爆炸過程中，會(huì)涉及到多種物質(zhì)的混合，只有各物質(zhì)界面能實(shí)現(xiàn)網(wǎng)格類型的精確區(qū)分，才能得到正確的計(jì)算結(jié)果，故需對現(xiàn)有的Youngs界面技術(shù)進(jìn)行改進(jìn)以解決其無法對兩種以上物質(zhì)混合進(jìn)行界面處理的問題。

針對上述問題，本文通過使用網(wǎng)格中量最大的物質(zhì)來代替量最少的物質(zhì)對Youngs界面技術(shù)進(jìn)行改進(jìn)，使其能夠處理多物質(zhì)混合界面，然后將其與Euler算法結(jié)合用來處理爆炸沖擊問題。最后對改進(jìn)的算法進(jìn)行了并行設(shè)計(jì)及程序編寫，并對空中爆炸和聚能射流進(jìn)行了數(shù)值模擬，證明了算法及程序的正確性。

1 改進(jìn)的多物質(zhì)歐拉界面處理算法概述

在笛卡爾坐標(biāo)系下，Euler方法對控制方程組采用算子分裂算進(jìn)行離散。一方面，按物理效應(yīng)分裂為Lagrange步和Euler步，Lagrange步計(jì)算偏應(yīng)力、壓力對內(nèi)能和速度的影響，Euler步計(jì)算各網(wǎng)格之間質(zhì)量、體積、動(dòng)量和能量輸運(yùn)。另一方面，通過交替輸運(yùn)方法將多維輸運(yùn)問題轉(zhuǎn)化為多個(gè)一維問題，同時(shí)能保證按網(wǎng)格間物質(zhì)流通的方向完成輸運(yùn)。

在Euler算法中，當(dāng)混合網(wǎng)格中有3種物質(zhì)的時(shí)候，無法直接使用Youngs界面技術(shù)進(jìn)行界面處理。為了使計(jì)算繼續(xù)進(jìn)行，將Youngs界面技術(shù)進(jìn)行如下改進(jìn)：首先按物質(zhì)的體積份額進(jìn)行從大打小排序，用最大的物質(zhì)替代最少的物質(zhì)，如此就將3種物質(zhì)就轉(zhuǎn)化為兩種物質(zhì)了，從而就可以繼續(xù)應(yīng)用Youngs界面技術(shù)進(jìn)行界面處理。

2 算法的并行化及實(shí)現(xiàn)

2.1 Euler方法的并行計(jì)算理論分析

在Euler方法數(shù)值模擬的過程中，網(wǎng)格固定，物質(zhì)在網(wǎng)格間流進(jìn)流出。各網(wǎng)格的物質(zhì)輸運(yùn)會(huì)導(dǎo)致其周圍網(wǎng)格部分物理量的變化，而在輸運(yùn)算法中隱含的子區(qū)域間的關(guān)聯(lián)性和數(shù)據(jù)相關(guān)性不容易被發(fā)現(xiàn)，而這兩個(gè)問題僅存在于Euler的輸運(yùn)步中，是Euler方法并行計(jì)算[3-6]實(shí)現(xiàn)的關(guān)鍵。

2.1.1 數(shù)據(jù)子區(qū)域的劃分和特點(diǎn)

并行計(jì)算需要將計(jì)算域劃分為若干塊來獨(dú)立進(jìn)行計(jì)算，那么各子區(qū)域間必然會(huì)存在重疊和相關(guān)的部分，具體如圖1所示。

圖1 重疊子區(qū)域

網(wǎng)格物理量的改變量等于該網(wǎng)格3個(gè)方向上的流進(jìn)量減去流出量的差值，故Euler方法的輸運(yùn)歩里各網(wǎng)格的更新應(yīng)包含該網(wǎng)格與鄰近網(wǎng)格的輸運(yùn)。

以圖2中的單方向輸運(yùn)為例，流場速度為正，則單方向k網(wǎng)格一次物理量的改變包含k-1和k網(wǎng)格的輸運(yùn)，其值為k-1網(wǎng)格流進(jìn)量減去k網(wǎng)格流出量的差值。

圖2 單向輸運(yùn)

因此對于圖1中的區(qū)域4而言，需在其下邊和左邊分別增加一層網(wǎng)格，具體如圖3所示。

圖3 輸運(yùn)相關(guān)性

此外，在實(shí)際運(yùn)算中，k網(wǎng)格物理量的數(shù)值更新也需要知道k，k+1網(wǎng)格物理量的值，因此4號(hào)區(qū)域的右邊和上邊也要增加一層網(wǎng)格，但所增加的網(wǎng)絡(luò)為虛擬網(wǎng)格，無需參與數(shù)據(jù)更新。在實(shí)際的并行算法設(shè)計(jì)中，要盡量避免在虛擬網(wǎng)格上進(jìn)行通信，以免出錯(cuò)。

2.1.2 數(shù)據(jù)更新與通信的分析

由圖2的分析可知，單方向輸運(yùn)的時(shí)候，網(wǎng)格物理量的輸運(yùn)變化會(huì)導(dǎo)致自身網(wǎng)格和其下階段網(wǎng)格的物理量產(chǎn)生變化，但在實(shí)際計(jì)算中，由于CFL條件的限制，各網(wǎng)格的輸運(yùn)變化只會(huì)對其周圍26個(gè)網(wǎng)格的物理量產(chǎn)生更新和影響，即網(wǎng)格輸運(yùn)的數(shù)據(jù)相關(guān)性。

可以使用原本未被更新的物理量來消除網(wǎng)格輸運(yùn)的數(shù)據(jù)相關(guān)性，但這樣又會(huì)出現(xiàn)過量輸運(yùn)的問題。故綜合考慮，將3個(gè)方向的輸運(yùn)轉(zhuǎn)變?yōu)槿匮h(huán)下的單一方向的輸運(yùn)，這是本文進(jìn)行并行算法設(shè)計(jì)的基礎(chǔ)。

2.1.3 變量數(shù)據(jù)的定義

并行計(jì)算設(shè)計(jì)中，涉及大量變量的定義和MPI庫函數(shù)的調(diào)用。在變量的定義過程中，需注意各進(jìn)程內(nèi)私有變量的數(shù)據(jù)更新無法對其它進(jìn)程中該變量的數(shù)據(jù)進(jìn)行同步更新，這樣就會(huì)出現(xiàn)錯(cuò)誤，所以當(dāng)私有變量的數(shù)據(jù)更新會(huì)對周圍網(wǎng)格產(chǎn)生影響的時(shí)候，必須要將更新后的變量數(shù)據(jù)向相關(guān)進(jìn)程進(jìn)行傳遞。全局變量的數(shù)據(jù)一經(jīng)改變，所有進(jìn)程中該變量的數(shù)據(jù)都可以得到相應(yīng)的更新。

2.2 并行前處理

2.2.1 創(chuàng)建可以任意劃分的計(jì)算域

設(shè)計(jì)的并行程序能夠完成計(jì)算域3個(gè)方向任意維度的劃分，劃分原則如下：

①計(jì)算域最外層的循環(huán)劃分為較多的子域，內(nèi)層循環(huán)盡量少劃分子域；

②盡量將計(jì)算域劃分為正方體，以使通信面積最??；

③在劃分過程中，若3個(gè)方向的網(wǎng)格劃分?jǐn)?shù)目不同，則使外層網(wǎng)格最多。

設(shè)npx，npy，npz為X，Y，Z3個(gè)方向的分區(qū)數(shù)量(都為正整數(shù))，那么計(jì)算域子區(qū)域的數(shù)目np=npx×npy×npz，子區(qū)域的編號(hào)myid可通過MPI庫函數(shù)中的MPI_Comm_rank來獲取，計(jì)算機(jī)的進(jìn)程數(shù)目可通過MPI_Comm_size來獲取,具體實(shí)現(xiàn)代碼分別為：

MPI _Comm_ rank(MPI _Comm_world,myid, ierr)

MPI_Comm_size(MPI_omm_world,numprocs, ierr)

Mepx，Mepy，Mepz表示3個(gè)方向的進(jìn)程坐標(biāo)，具體如式(1)所示：

(1)

對于距計(jì)算邊界較近的子區(qū)域，若其周圍沒有其它子區(qū)域，可使用MPI 的虛擬化子區(qū)域MPI_Proc_null。

用來確定與目標(biāo)子區(qū)域相鄰子區(qū)域的編號(hào)的公式如式(2)所示：

(2)

假設(shè)取X，Y，Z3個(gè)方向的分區(qū)數(shù)量分別為 2、3、2，那么共有12個(gè)子區(qū)域，其中一種劃分方法如圖4所示。

圖4 三維拓?fù)浣Y(jié)構(gòu)示意圖

將X方向分為2塊，Y方向分為3塊，Z方向分為2塊，各進(jìn)程的編號(hào)myid的值為0,1,2，…，11。

2.2.2 任意劃分區(qū)域的實(shí)現(xiàn)

為了實(shí)現(xiàn)計(jì)算域能夠在3個(gè)方向上的任意劃分，最關(guān)鍵的就是要解決網(wǎng)格在3個(gè)方向的整除問題。所以需要在無法被整除的方向上附加一些網(wǎng)格，以X向?yàn)槔郊泳W(wǎng)格的數(shù)量如式3所示：

(3)

附加網(wǎng)格不參與數(shù)據(jù)計(jì)算，只進(jìn)行存儲(chǔ)和通信。但在進(jìn)行計(jì)算域劃分的時(shí)候，還是應(yīng)盡量保證整除，因?yàn)檫@樣設(shè)計(jì)的程序會(huì)較為簡單。

2.2.3 子區(qū)域的通信

本文設(shè)計(jì)的并行化算法需要實(shí)現(xiàn)計(jì)算域3個(gè)方向上的任意劃分，但其僅能保證一個(gè)方向上的數(shù)據(jù)連續(xù)通信，其它兩個(gè)方向的數(shù)據(jù)通過MPI函數(shù)庫的MPI_TYPE_HVECTOR的自定義數(shù)據(jù)類型來處理。其調(diào)用形式如下：

MPI_TYPE_HVECTOR(COUNT,BLOCKLENGTH,STRIDE,OLDTYPE,NEWTYPE,IERROR)

其中，COUNT表示整合數(shù)據(jù)的塊數(shù)，BLOCKLENGTH表示每個(gè)子區(qū)域中所含元素的個(gè)數(shù)，STRIDE表示各子區(qū)域起始位置間的字節(jié)數(shù)，OLDTYPE表示舊的數(shù)據(jù)類型，NEWTYPE表示新的數(shù)據(jù)類型。

由于Y，Z方向的數(shù)據(jù)均不連續(xù)，以yz平面為例，用戶通過自定義一個(gè)數(shù)據(jù)類型，將Y、Z方向的不連續(xù)數(shù)據(jù)進(jìn)行打包，使他們能夠進(jìn)行正常通信。具體代碼如下：

MPI_Type_Hvector (mysizey,1, (mysizex+2)*8,

MPI_R eal8,Txx,ierr) MPI_Type_Hvector(mysizez,1,(mysizex+2)*(mysizey+2)*8,Txx,Tx,ierr)

MPI_Type_commit(Tx,ierr) MPI_Type_free(Txx,ierr)

其中，Txx表示y方向的一組數(shù)據(jù)，其為最終Tx數(shù)據(jù)的中間過渡值，在完成最終的數(shù)據(jù)類型定義后，需將中間過渡的數(shù)據(jù)類型釋放掉。

2.3 并行后處理

并行計(jì)算的后處理主要包括對關(guān)鍵點(diǎn)的數(shù)據(jù)采集、可視化軟件的接口設(shè)計(jì)方面。針對不同的需求可編制不同的后處理程序，本文編制的后處理程序可將并行計(jì)算后的結(jié)果顯示出來，滿足所需的后處理要求。

3 并行程序設(shè)計(jì)及性能分析

3.1 并行程序設(shè)計(jì)

本文設(shè)計(jì)的并行程序，由主程序、前處理初始化程序、邊界條件處理子程序和主循環(huán)計(jì)算子程序組成[7]，各部分的功能如下所示。

1)主程序：程序的核心構(gòu)成，貫穿于程序的整個(gè)執(zhí)行過程，完成各子程序的調(diào)用，采用最節(jié)簡的語言設(shè)計(jì)；

2)前處理初始化程序：主要完成人機(jī)交互，能讓用戶選擇計(jì)算的種類，并進(jìn)行有關(guān)參數(shù)的設(shè)置。同時(shí)要完成并行計(jì)算循環(huán)的起始位置設(shè)置，各進(jìn)程的位置及與其相關(guān)進(jìn)程的編號(hào)；

3)邊界條件處理子程序：處理計(jì)算模型的壓力，速度，網(wǎng)格標(biāo)志，介質(zhì)標(biāo)志等物理量的邊界條件；

4)主循環(huán)計(jì)算子程序：主循環(huán)程序包括人工粘性步、Lagrange步、Euler步、修正步和輸出子程序組成。其中，人工粘性步用來計(jì)算人工粘性、損傷度和偏應(yīng)力等；Lagrange步用來計(jì)算壓力梯度對質(zhì)團(tuán)的速度、能量、動(dòng)量和比內(nèi)能的影響；Euler步使用改進(jìn)的Youngs界面處理方法對能量、質(zhì)量、體積和動(dòng)量進(jìn)行輸運(yùn)；修正步子程序主要對計(jì)算的數(shù)值進(jìn)行修正，以保證計(jì)算能夠繼續(xù)保持下去；輸出子程序主要完成計(jì)算結(jié)果的二維、三維切片輸出。

整個(gè)并行程序的執(zhí)行流程如圖5所示。

圖5 程序執(zhí)行流程圖

3.2 性能分析

1)并行性能測試

選取三介質(zhì)聚能射流模型[8]為測試算例，處理92×92×300、151×151×200、151×151×300和200×200×300 4個(gè)不同網(wǎng)格數(shù)的模型，分別記為1，2，3，4號(hào)模型。測試共分為8個(gè)區(qū)，計(jì)算步驟為5 000步，取每100步所需的時(shí)間作為其運(yùn)行的平均速度。各模型的加速比如表1所示。

從表1中可以看出：

①同一模型中，分區(qū)越多的計(jì)算速度越快；

②不同網(wǎng)格模型中，網(wǎng)格數(shù)越多，其加速比也越高。

表1 不同模型的加速比

2)不同分區(qū)的算例測試

選取空中爆炸模型，炸藥類型為TNT，裝藥密度為1.65 g/cm3，裝藥半徑為1 m，裝藥量為6.9 t，爆速為6.970 mm/ μs，空氣密度為1.29 × 10-3g/cm3。分區(qū)統(tǒng)一劃分為1×1×1，1×3×1，1×1×3和3×1×1四種模型，其中1×1×1模型算例等價(jià)于于串行算法，而其它模型代表3種不同的分區(qū)方法。

對算例1及算例2進(jìn)行全局監(jiān)測，對算例3采用10個(gè)關(guān)鍵點(diǎn)進(jìn)行壓力監(jiān)測。測試結(jié)果表明：算例1、算例2爆炸在各分區(qū)的檢測結(jié)果基本一致，分別如圖6、圖7所示。

圖6 算例1的二維切片結(jié)果

圖7 算例2的二維切片結(jié)果

算例3的10個(gè)關(guān)鍵點(diǎn)各分區(qū)的壓力測試結(jié)果如表2所示。

表2 算例3不同分區(qū)的沖擊波峰值超壓值

由表2中可以看出，1×3×1，1×1×3和3×1×13個(gè)分區(qū)關(guān)鍵點(diǎn)的空中爆炸沖擊波波峰數(shù)值完全一致，但卻與1×1×1的劃分區(qū)域有些許區(qū)別，主要原因如下：1)網(wǎng)格的劃分會(huì)使串行計(jì)算數(shù)值和并行計(jì)算數(shù)值有一定的誤差；2)程序設(shè)計(jì)中，很多網(wǎng)格介質(zhì)標(biāo)志的判定方法比較粗糙存在誤判問題。

從上述的算例測試可以看出，不同的分區(qū)不會(huì)對計(jì)算結(jié)果產(chǎn)生影響，證明了程序可以對計(jì)算域進(jìn)行任意劃分。

4 程序的實(shí)例仿真

4.1 空中爆炸

選取炸藥類型為TNT，裝藥密度為1.65 g/cm3，裝藥半徑為1 m，裝藥量為6.9 t，爆速為6.970 mm/μs，空氣密度為1.29×10-3g/cm3，爆轟產(chǎn)物的多方指數(shù)k0取3，爆轟產(chǎn)物充分膨脹下的多方指數(shù)k1取3.16，邊界條件為連續(xù)性邊界，計(jì)算區(qū)域大小為100 m×10 m×10 m，網(wǎng)格劃分為400×40×40。

根據(jù)前人大量的實(shí)驗(yàn)，空中爆炸沖擊波峰值的超壓計(jì)算公式如式3所示：

(3)

式中，ΔPm表示空氣沖擊波波峰超壓(0.1 Mpa)，w表示裝藥量(kg)，r表示距爆炸中心的距離(m)。根據(jù)式(3)及本文設(shè)計(jì)程序的數(shù)值模擬后的結(jié)果如表3所示。

表3 沖擊波峰值超壓數(shù)值解和經(jīng)驗(yàn)公式解的比較

從表中可以看出，本文設(shè)計(jì)的程序數(shù)值計(jì)算結(jié)果與前人的實(shí)驗(yàn)結(jié)果基本一致，說明了本文所設(shè)計(jì)算法及程序的正確性。

4.2 聚能射流

選取的計(jì)算模型藥型罩為等壁厚錐形，厚度2 mm，材料45#鋼，藥柱的錐形角120°，高33 mm，直徑40 mm，邊界條件為連續(xù)性邊界，計(jì)算域大小為46 mm×46 mm×150 mm，網(wǎng)格步長為0.5 mm。

應(yīng)用本文所設(shè)計(jì)程序的模擬結(jié)果如圖9所示。

圖9 聚能射流二維切片的程序模擬圖

實(shí)際的聚能射流如圖10所示。

圖10 聚能射流的實(shí)際效果圖

從圖9的程序模擬結(jié)果可以清楚地看出，聚能射流過程中共有3種介質(zhì)交互作用，相互之間界面清晰，說明改進(jìn)的Youngs算法的界面處理是非常精確的。

從圖9和圖10中可以看出，在金屬射流不斷伸長的過程中，從頭部到尾部，速度逐漸減小，且圖9的模擬結(jié)果和圖10的實(shí)驗(yàn)結(jié)果基本一致，證明了算法及程序的正確性。

5 結(jié)論

本文設(shè)計(jì)了基于改進(jìn)Youngs界面處理技術(shù)的Euler并行算法，并采用MPI標(biāo)準(zhǔn)消息傳遞接口編制了流體彈塑性動(dòng)力學(xué)的并行程序，然后對程序進(jìn)行了并行性能分析和算例測試，結(jié)果表明并行算法對程序的加速比有較大的提高，而不同分區(qū)對并行算法的計(jì)算結(jié)果沒有影響。最后通過對空中爆炸和聚能射流的實(shí)例仿真，證明了并行算法及程序的正確性，對爆炸沖擊問題的研究具有重要意義。

[1] Zhang W Y, Ma T B, Ning J G. Visualized computing for explosion disasters, progress in safety science and technology (Vol. VI)[M].Science Press, 2012:1214-1218.

[2] Wang G, Wang J T, Liu K X. New numerical algorithm in SUPER CE/SE and their application in explosion mechanics[J]. Science China Physics,Mechanics& Astronomy. 2011,52(2): 237-243.

[3] 馬天寶, 費(fèi)廣磊, 張文耀. 三維多物質(zhì)彈塑性流體動(dòng)力學(xué) Euler 方法的并行算法研究及程序測試[J].高壓物理學(xué)報(bào), 2011,25(6): 508-513.

[4] 馬天寶, 郝 莉, 寧建國. Euler 多物質(zhì)流體動(dòng)力學(xué)數(shù)值方法中的界面處理算法[J]. 計(jì)算物理, 2012, 25(2): 133- 138.

[5] 夏 銳, 肖明清, 賴 根. 并行測試系統(tǒng)設(shè)計(jì)與開發(fā)[J]. 計(jì)算機(jī)測量與控制, 2015, 14(7):841-843.

[6] Ning J G, Ren H L, Li P. Mechanical behaviors and damage constitutive model of ceramics under shock compression [J]. Acta Mechanica Sinica, 2011, 24(3): 305-315.

[7] 王生武, 石秀華. 渦輪葉尖壓力邊小翼肋條對泄漏流場的數(shù)值模擬[J]. 計(jì)算機(jī)測量與控制, 2009,17(8):1527-1530.

[8] 費(fèi)廣磊, 馬天寶, 寧建國. 基于 MPI 的聚能射流問題的并行計(jì)算[A]. 中國力學(xué)學(xué)會(huì)2012 學(xué)術(shù)大會(huì)論文集[C]．2012.

3D More Material Euler Interface Processing Parallel Algorithm Research

Long Wenjia，Li Yuexin

(Computer and Information Engineering Department, Hubei University, Wuhan 430000,China)

Of explosion shock in view of the existing algorithms can not effectively study this problem, on the basis of the technology improved the Youngs interface design for more than three dimensional material Euler interface processing parallel algorithm. Then using MPI standard to the program design of the algorithm and program for the parallel performance testing and different partitions example test, the test results show that the parallel algorithm for acceleration is larger, the different partition has no effect on the calculation results of parallel algorithm. Finally application programs written by the air blast and shaped jet are simulated, the results show that the simulation results and the program is conform to the predecessors' research data, illustrates the effectiveness of the program.

explosion and impact;numerical simulation; Euler's method; parallel computing; youngs interface technology

2015-08-27；

2015-09-25。

湖北省科技支撐項(xiàng)目(2014BAA089)。

龍文佳(1978-)，女，湖北荊州人，講師，碩士，主要從事數(shù)據(jù)庫方向的研究。

李躍新(1958-),男，湖北武漢人，副教授，博士，主要從事知識(shí)工程與并行計(jì)算方向的研究。

1671-4598(2016)06-0226-04

10.16526/j.cnki.11-4762/tp.2016.06.062

TM417

A