基于高性能計算的離散介質(zhì)沖擊過程

吳志平，劉波平，王 康，李石濱，胡畢煒，胡必偉，游 杰

(江西省科技基礎條件平臺中心，江西 南昌 330003)

0 引 言

散體介質(zhì)廣泛分布于沙漠、海洋、礦山等自然界[1-3]，也被應用于化學和生物領域[4-6]或工業(yè)裝備的減震隔振[7-9]等領域。其分布甚至不局限于地球。在行星探測領域，月球表層土壤、火星土壤等均為散體介質(zhì)。隨著我國對月球及火星探索工作的開展，在探測器降落月球或者火星表面的過程中，無論是發(fā)動機吹起的表面土壤的過程，還是探測器支撐腳與火星或月球土壤之間的沖擊行為，都存在與散體介質(zhì)的相互作用，也涉及各種離散介質(zhì)間的沖擊過程，因此對離散介質(zhì)沖擊問題的研究具有重要意義[10-11]。對于大規(guī)模顆粒介質(zhì)的模擬，例如：對流化床[12-14]、風沙流[15]的模擬最常用的方法是基于高性能計算平臺的離散元模擬，在上述模擬過程中，顆粒的恢復系數(shù)是一個重要的參數(shù)。

為了對離散介質(zhì)沖擊過程進行精確描述，本文對離散介質(zhì)沖擊過程的速度變化過程進行模擬，研究離散介質(zhì)沖擊過程恢復系數(shù)的變化規(guī)律。關于這個課題已有相當多的研究，對離散介質(zhì)沖擊過程中動力學行為給出了較為細致的分析，如模型選取[16-17]和鏈長度對結果的影響[18-19]等。Liu等[18]在定義了接觸模型的基礎上，將對時間的積分轉換為對沖量的積分，解決了因時間步長的選取造成的計算量過大以及精度問題。此外，Gharib等[19]發(fā)現(xiàn)，一組不同粒徑的鏈在不同排布下，其沖擊過程會產(chǎn)生顯著差別。但在上述研究中，很少關注不同直徑的離散介質(zhì)沖擊過程的宏觀統(tǒng)計規(guī)律，如恢復系數(shù)的統(tǒng)計規(guī)律對離散介質(zhì)沖擊過程的影響。

基于高性能計算平臺的仿真模擬是科學研究的重要方向，特別是對離散介質(zhì)這類數(shù)量十分龐大的物質(zhì)，其模擬研究必須借助超級計算機等高性能計算平臺進行。李勇俊[20]使用高性能計算平臺對離散介質(zhì)堆積生成進行研究，Gan等人[21]對工業(yè)生產(chǎn)中的大顆粒搬運、Herman等人[22]對離散介質(zhì)攪拌等領域使用高性能計算平臺進行了模擬。Yue等人[23]使用高性能計算平臺對顆粒填充過程進行了模擬。劉麗芝等人[24]使用高性能計算平臺對離散介質(zhì)的沖擊過程進行模擬。此外，流體與離散介質(zhì)的耦合過程也是研究的重點，高性能計算平臺在該領域也得到了重要的應用[25]。

但是，使用高性能計算平臺對顆粒鏈沖擊過程進行模擬較為少見，特別是對多分散顆粒鏈的統(tǒng)計研究更為罕見。因此，本文通過使用高性能計算平臺進行并行計算，研究單分散鏈、多分散鏈與入射介質(zhì)末速度、恢復系數(shù)與壓縮階段沖量的變化規(guī)律，以及與傳統(tǒng)恢復系數(shù)的對比。

1 接觸模型和計算模擬方法

1.1 沖擊過程的動力學模型

圖1所示為一維離散顆粒鏈，由n個相同材料的均質(zhì)顆粒緊密排列構成，由左至右依次編號為1～n，所有顆粒的質(zhì)心在一條直線上。初始狀態(tài)下，顆粒之間恰好相互接觸但無相互作用力。編號為1的顆粒為碰撞顆粒，即入射顆粒，編號為2～n的顆粒組成顆粒鏈的其余部分。這里顆粒鏈的半徑考慮為單分散和高斯分布多分散2種情況，當考慮顆粒鏈為單分散分布時，入射顆粒的半徑與顆粒鏈其余顆粒的半徑相同。入射顆粒的初始速度為v1=v0，其余顆粒初始狀態(tài)下都靜止。

參照離散單元法(DEM)，顆粒只與其兩側與之接觸的顆粒產(chǎn)生相互作用。設任一個顆粒i的質(zhì)量、粒徑、速度和位置分別為mi、Ri、vi以及xi，相鄰顆粒對顆粒i的接觸力為Fi-1,i和Fi,i+1，其中入射顆粒只受到顆粒鏈中顆粒2向左的作用力，而顆粒鏈最右端的顆粒n只受到它左側顆粒n-1向右的作用力。因此顆粒鏈中每一個顆粒的離散牛頓動力學方程組為：

(1)

(2)

(3)

對于顆粒與顆粒之間的接觸力，存在多種通過壓縮量、相對速度等參數(shù)定義的接觸力模型[26]，這里采用Thornton模型[27]計算：

(4)

δp,q=Rp+Rq-|xp-xq|

(5)

vp,q=vq-vp

(6)

以顆粒鏈左端為坐標原點建立坐標系，因此各顆粒的初始條件為：

(7)

(8)

當顆粒鏈中的每一個顆粒均與其余顆粒不發(fā)生接觸，并且所有顆粒之間相互遠離：

δi,i+1<0,vi

(9)

即所有顆粒都在相互遠離，便可認為顆粒鏈的沖擊過程結束。通過上述方程和條件，就可以模擬顆粒與一維顆粒鏈的碰撞過程。

對于恢復系數(shù)，有2種常用的定義，一種是基于碰撞前后的相對速度之比，其定義如下：

(10)

(11)

其中，pc、pf分別為碰撞過程中，2個顆粒壓縮階段和恢復階段的沖量。對于2個顆粒對心正碰，eN和eP完全等價。

1.2 計算模擬方法

以聚四氟乙烯(PTFE)為例，采用數(shù)值模擬研究不同入射顆粒的初速度、入射顆粒與顆粒鏈顆粒半徑的粒徑比、顆粒鏈長度等情況下，入射顆粒的回彈速度、顆粒鏈顆粒的末速度、顆粒之間的碰撞恢復系數(shù)、能量耗散率等結果的變化規(guī)律。該材料的參數(shù)如表1所示。

表1 PTFE的部分力學參數(shù)

這里假設聚四氟乙烯PTFE是一種完全的粘彈性材料，遵循完全的彈塑性壓縮和恢復過程。為確保足夠的精度，本文所有模擬中的時間步長均為1.0×10-9s。

在確定了時間步長、材料參數(shù)和顆粒的粒徑分布形式，接觸模型以及其他必須變量之后，便可以對整個碰撞過程進行模擬。圖2給出了本文模擬過程的流程圖。

碰撞顆粒的半徑是影響入射顆粒碰撞后速度的關鍵因素。實際顆粒系統(tǒng)中，即便顆粒鏈中顆粒個數(shù)和大小不變，但是因為不同大小顆粒的排列不同，導致與入射顆粒碰撞的顆粒的半徑不同，因此入射顆粒碰撞后的速度有所不同。為了描述這種情況導致的入射顆粒碰撞速度的變化，下面模擬不同種顆粒排列下入射顆粒碰撞速度分布和對不同粒徑分布下的沖擊過程。由于計算量較大，基于現(xiàn)有超算平臺的多個高性能計算節(jié)點以及MPI并行算法庫，對顆粒鏈沖擊過程的算法進行并行化處理，將計算任務分配到各個從節(jié)點上，并通過MPI算法庫將所有從節(jié)點進程的計算結果匯總到主節(jié)點上，最后對數(shù)據(jù)結果進行統(tǒng)計。本文所使用模擬程序的偽代碼段如算法1所示。

算法1 模擬過程的偽代碼。

Algorithm: Discrete Material Chain Impact

Require: equation (7), equation (8) and necessary parameters

// Initialization of MPI

MPI_Init()

MPI_Comm_size()

MPI_Comm_rank()

do i=myid,n,numprocs

// n is the times of simulation

initialization of the chain

sp←true// the sign of the end of the impact

while sp is true

//the main procedure of the impact in different process

use the equation (5) and equation (6) to get the compressions and relative velocities

get the contact force by equation (4)

use equation (1)-equation (3) to update the locations and the velocities

use the equation (9) to determine whether the sign sp is true

end while

aggregate the data to the main process (process0) using MPI_Gather()

MPI_Barrier()

//interprocess synchronization

end do

If myid=0 // all of statistics were on main process

count the data and output results in process0

end if

MPI_Finalize( )

基于高性能計算平臺，本文采用獨立開發(fā)的離散介質(zhì)沖擊過程模擬程序，對整個沖擊過程進行模擬。

2 結果與討論

下面，研究給定不同入射初速度、顆粒鏈粒徑分布、顆粒鏈長度等情況下，由PTFE組成的顆粒與一維顆粒鏈碰撞后顆粒鏈中顆粒的速度，研究入射顆粒碰撞恢復系數(shù)、顆粒鏈能量留存率。顆粒的質(zhì)量根據(jù)材料密度和顆粒的半徑計算。

首先為了驗證使用的程序的準確性，本文對比了不同入射速度下，2個顆粒對心碰撞恢復系數(shù)的模擬結果與理論結果。如圖3(a)的對比結果可知，本文所使用的模擬程序有著較高的精確性，其誤差在1%以下，能準確模擬整個彈塑性Thornton接觸過程。之后，本文模擬了顆粒與顆粒鏈碰撞后，各顆粒碰撞速度隨時間變化，如圖3(b)所示。顆粒鏈的長度為10，單分散顆粒鏈顆粒的半徑為150 μm；多分散顆粒鏈顆粒的半徑服從正態(tài)分布，均值為150 μm，標準差為30 μm。如圖3所示，碰撞結束后，顆粒鏈中每一個顆粒的速度最終都會達到穩(wěn)定，并且顆粒的速度在達到穩(wěn)定后都呈現(xiàn)依序號遞增的狀態(tài)，該結果與Liu等人[18]的模擬結果相似，驗證了本文算法的準確性。

(12)

圖3顯示，顆粒與顆粒鏈碰撞后速度最終達到穩(wěn)定，但由于碰撞過程中能量損失，顆粒鏈長度對入射顆粒末速度有一定影響。同時，將入射顆粒的泊松恢復系數(shù)隨顆粒鏈長度變化繪制于圖4(a)，顯示入射顆粒的碰撞恢復系數(shù)隨著顆粒鏈長度增加而增大。當顆粒鏈長度大于7時，入射顆粒的泊松恢復系數(shù)不再改變，保持穩(wěn)定，并且入射顆粒在壓縮階段的沖量pc不隨顆粒鏈長度的變化而變化，見圖4(b)，其中Δv=pc/m1為壓縮階段的速度變化量。

對于不同粒徑分布下的情況，分別將統(tǒng)計結果繪制于圖5。從圖5(a)可以看出，入射顆粒與多分散顆粒鏈碰撞后的速度應該與顆粒鏈不同大小的顆粒排列相關，因為不同大小顆粒隨機排列，導致入射顆粒碰撞后速度的方向可能與原來速度相同，也可能相反。其次，入射顆粒碰撞后的速度分布較寬，類似于正態(tài)分布，隨著模擬次數(shù)的提高，這個分布已經(jīng)穩(wěn)定，因此后續(xù)關于多分散顆粒鏈碰撞恢復系數(shù)等研究，都要通過多次碰撞模擬，計算入射顆粒的平均碰撞速度和標準差。圖5(b)～圖5(c)分別是多分散顆粒鏈中入射顆粒碰撞速度分布受初始速度的影響以及受入射顆粒與被碰撞顆粒半徑比的影響。由圖5(b)可知，初始速度越大，入射顆粒碰撞后速度分布越寬且平均速度越大。給定入射顆粒初始速度以及粒徑比R2/R1，碰撞結束時入射顆粒的末速度分布近似于正態(tài)分布。從圖5(c)可以很明顯地看出，當粒徑比大于1.4時，碰撞后的入射顆粒末速度分布集中于很小的一個范圍內(nèi)。由此可推測，當顆粒鏈顆?？偙慌鲎差w粒尺寸遠大于入射顆粒尺寸時，即R2/R1?1.0，入射顆粒的末速度將近似于2個顆粒碰撞的情形。

圖5(d)比較了2個顆粒碰撞與顆粒鏈碰撞的eP隨半徑比變化情況。從圖5(d)可以看出，R2/R1≤1.2時，2個顆粒碰撞情形與一維顆粒鏈碰撞情形下的eP差異較大，而當R2/R1>1.4時，2種情形下eP趨于一致。這表明，當粒徑比R2/R1>1.4時，入射顆粒與顆粒鏈碰撞可近似為2個顆粒碰撞。

接下來，給定相同的初始速度，對比單分散和多分散鏈的沖擊過程。入射顆粒eP隨初始速度的變化規(guī)律，如圖6所示。同時與無限大半空間的碰撞結果也繪制于圖6。由圖6(a)可以看出，相較于傳統(tǒng)恢復系數(shù)即顆粒與無限大半空間的沖擊后的速度比，單分散鏈與多分散鏈的泊松恢復系數(shù)eP差別不大，特別是在入射速度增加時，差別可以忽略，且大于傳統(tǒng)沖擊過程的恢復系數(shù)。如圖6(b)所示，入射顆粒的末速度遠小于與無限大半空間沖擊的情況。這說明鏈對入射顆粒的能量耗散要顯著強于與無限大半空間碰撞。其原因是入射顆粒與鏈的沖擊過程中，所受到的沖量大于與無限大球的碰撞過程，因而能量耗散更高。這從理論上驗證了離散介質(zhì)在抗振減振方面的顯著效果。

基于高性能計算平臺，本文研究算法1所給出的并行算法的效率。圖7對比了不同顆粒鏈長度下，模擬1000次在不同核數(shù)下的加速比和并行效率。圖7(a)展示了上述條件下的加速比，隨著核數(shù)的增加，本文所使用的并行算法的加速比基本上呈線性變化。在核數(shù)增加到20核情況下，如圖7(b)所示，本文并行算法的并行效率保持在88%左右，顯示了較好的并行效率。

3 結束語

基于高性能計算平臺，本文模擬了牛頓擺，即單個顆粒沖擊無邊界顆粒鏈的過程，對不同入射速度、顆粒鏈尺度分布條件下，顆粒鏈中顆粒的速度變化規(guī)律特別是入射顆粒的碰撞后速度進行了詳細分析。發(fā)現(xiàn)入射顆粒與顆粒鏈碰撞后，顆粒鏈中各個顆粒的速度最終將保持不變；泊松恢復系數(shù)在顆粒鏈長度較小時，隨顆粒鏈長度增長而增大，最后趨于穩(wěn)定。而入射顆粒壓縮階段的沖量不隨顆粒鏈長度的增長而變化，只與入射速度線性相關。對于多分散顆粒鏈，當被碰撞顆粒與入射顆粒的半徑比大于某一個值時，入射顆粒沖擊顆粒鏈的泊松恢復系數(shù)與2個顆粒的碰撞基本相等，可以忽略顆粒鏈后續(xù)顆粒對入射顆粒碰撞速度的影響，這個值與沖擊速度和材料參數(shù)相關，比如沖擊速度3 m/s的PTFE顆粒，這個值是1.4。對于顆粒與多分散顆粒鏈的沖擊過程，入射顆粒的末速度分布與初速度和粒徑比相關，且分布呈現(xiàn)正態(tài)分布形式。本文使用高性能計算平臺對顆粒鏈的沖擊過程進行模擬和統(tǒng)計分析，并給出了不同核數(shù)下并行效率和加速比的變化關系，結果顯示，本文并行算法具有較好的并行效率。

本文驗證了對離散介質(zhì)的沖擊過程，其厚度只需要7～8層，之后的顆粒不會對沖擊過程產(chǎn)生顯著影響，可以大幅減少諸如流化床在內(nèi)的各種離散介質(zhì)的計算量。并且，根據(jù)本文的計算結果，對于多分散一維離散介質(zhì)的沖擊過程，可以基于高斯分布的概率分布，而非通過傳統(tǒng)的數(shù)值模擬描述其沖擊過程，極大地簡化了模擬離散介質(zhì)沖擊過程所需要的計算量。