一種高效率巖土體離散元數(shù)值計(jì)算顆粒搜索和運(yùn)行方法
【專利摘要】一種高效率巖土體離散元數(shù)值計(jì)算顆粒搜索和運(yùn)行方法，該方法包括基于劃定空間網(wǎng)格和顆粒窗口搜索最優(yōu)鄰近顆粒，以及通過(guò)累積顆粒位移控制和減少搜索算法運(yùn)行。對(duì)于巖土體等非流體數(shù)值模擬，使用位移矩陣控制窗口搜索算法，通常數(shù)萬(wàn)次迭代才需要運(yùn)行一次搜索算法，顯著提高了離散元數(shù)值方法的速度，可實(shí)現(xiàn)巖土體大規(guī)模顆粒系統(tǒng)的動(dòng)態(tài)模擬。
【專利說(shuō)明】一種高效率巖土體離散元數(shù)值計(jì)算顆粒搜索和運(yùn)行方法

【技術(shù)領(lǐng)域】
[0001]本發(fā)明涉及巖土體離散元數(shù)值模擬，尤其是基于窗口搜索和顆粒位移控制的高效率計(jì)算方法。

【背景技術(shù)】
[0002]離散元法是解決巖土體不連續(xù)問(wèn)題的一種主要數(shù)值模擬方法，可以動(dòng)態(tài)模擬巖體的大變形、斷裂破壞，在巖土、礦冶、環(huán)境等領(lǐng)域取得了較為廣泛的應(yīng)用，如大型高邊坡穩(wěn)定變形機(jī)理、深埋地下工程圍巖破壞、礦山崩落開采等。離散元對(duì)模型在空間上進(jìn)行分解，在時(shí)間上進(jìn)行迭代計(jì)算，其計(jì)算量非常巨大。離散元計(jì)算主要包含兩個(gè)部分:搜索每個(gè)顆粒的鄰近接觸顆粒和迭代計(jì)算顆粒運(yùn)動(dòng)。其中顆粒搜索計(jì)算需要從顆粒系統(tǒng)中找出與當(dāng)前顆粒接觸的顆粒，然后進(jìn)一步計(jì)算顆粒受力和運(yùn)動(dòng)。當(dāng)顆粒數(shù)量較多時(shí)，顆粒搜索的運(yùn)算量急劇增加，成為大規(guī)模顆粒離散元數(shù)值模擬的瓶頸。因此，需要在搜索算法本身以及運(yùn)行方法上對(duì)其進(jìn)行改進(jìn)，以滿足復(fù)雜巖土工程問(wèn)題對(duì)大規(guī)模顆粒離散元數(shù)值模擬的要求。


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

[0003]為了克服大規(guī)模離散元數(shù)值計(jì)算搜索算法計(jì)算量大，計(jì)算效率低的問(wèn)題。本方法提出一種高效率巖土體離散元數(shù)值計(jì)算顆粒搜索和運(yùn)行方法。通過(guò)建立空間網(wǎng)格和劃定顆粒窗口，迅速查找最優(yōu)鄰近顆粒。該方法采用位移累積矩陣減少搜索算法的運(yùn)行頻率，能有效地提高離散元法的計(jì)算效率，實(shí)現(xiàn)大規(guī)模顆粒離散元數(shù)值計(jì)算模擬和復(fù)雜巖土工程問(wèn)題快速模擬。
[0004]本發(fā)明為解決離散元快速計(jì)算問(wèn)題提出的技術(shù)方案是，一種高效率巖土體離散元數(shù)值計(jì)算顆粒搜索和運(yùn)行方法，其步驟包括:
[0005](I)建立初步鄰近顆粒數(shù)組；依據(jù)顆粒坐標(biāo)，建立具有特定單元邊長(zhǎng)的空間立方單元網(wǎng)格，以覆蓋所有顆粒，并建立相應(yīng)的網(wǎng)格編號(hào)矩陣；計(jì)算顆粒在網(wǎng)格位置，將每個(gè)顆粒的鄰近網(wǎng)格中的顆粒存儲(chǔ)為初步鄰近顆粒數(shù)組，每個(gè)顆粒相對(duì)于其鄰近顆粒為中心顆粒;
[0006](2)采用窗口法計(jì)算最優(yōu)鄰近顆粒；以每個(gè)顆粒為中心(即中心顆粒)建立特定大小的窗口，搜索初步鄰近顆粒位于窗口內(nèi)的顆粒，并存儲(chǔ)于最優(yōu)鄰近顆粒數(shù)組；而窗口外顆粒至少沿坐標(biāo)軸方向移動(dòng)dS距離才能與中心顆粒接觸；
[0007](3)利用位移累積矩陣控制搜索算法；在離散元迭代運(yùn)算中，累積計(jì)算每個(gè)顆粒位移矩陣D [X，Y, Z]，當(dāng)D某一位移分量超過(guò)dS/2時(shí)，意味著可能有新顆粒與中心顆粒接觸，運(yùn)行步驟1-2重新計(jì)算最優(yōu)鄰近顆粒數(shù)組，并重置顆粒位移矩陣D為O。
[0008]步驟(2)中應(yīng)用窗口法計(jì)算最優(yōu)鄰近顆粒，得到最優(yōu)鄰近顆粒數(shù)組，并通過(guò)位移矩陣來(lái)控制窗口搜索算法的運(yùn)行；數(shù)組也是一維矩陣；
[0009]步驟10構(gòu)建初始的顆粒離散元堆積模型，顆?？倲?shù)為m，顆粒由I開始連續(xù)編號(hào)到m ；
[0010]步驟11通過(guò)空間網(wǎng)格劃分計(jì)算得到初步鄰近顆粒；網(wǎng)格單元的邊長(zhǎng)為wg，所有顆粒均包含在網(wǎng)格中；
[0011]步驟12通過(guò)建立半徑為Rw顆粒窗口確定最優(yōu)鄰近顆粒；由此獲得顆粒窗口內(nèi)最優(yōu)鄰近顆粒數(shù)組，將顆粒位移矩陣D [X，Y, Z]全部數(shù)值設(shè)為零，重新開始累積顆粒位移；(步驟31里窗口半徑為:RW = l+dS+R.，這個(gè)半徑保證了窗口外顆粒運(yùn)動(dòng)dS距離才可能與窗口內(nèi)顆粒接觸，那么通過(guò)位移矩陣?yán)鄯e顆粒位移，當(dāng)顆粒位移超過(guò)dS/2時(shí)，才需要運(yùn)行搜索算法。)
[0012]步驟13利用最優(yōu)鄰近顆粒數(shù)組完成離散元數(shù)值計(jì)算；在每次迭代中，僅需判斷中心顆粒與最優(yōu)鄰近顆粒矩陣中的顆粒接觸關(guān)系，并計(jì)算相互作用力；一次迭代計(jì)算完成后，得到中心顆粒的位移量d[X，Y，Z]；
[0013]步驟14累積中心顆粒位移矩陣D[X，Y, Z] = D[X, Y, Z]+d[X, Y, Z];通過(guò)累積每次迭代中的中心顆粒位移，位移矩陣D記錄了從前一次最優(yōu)鄰近顆粒搜索后的顆粒位移；
[0014]步驟15判斷max |D [X，Y，Z] I是否大于dS/2 ;即判斷顆粒位移在某一坐標(biāo)軸方向上是否超過(guò)dS/2。如果其大于dS/2，返回步驟11重新計(jì)算最優(yōu)鄰近顆粒；反之，進(jìn)入步驟16 ；
[0015]步驟16根據(jù)迭代次數(shù)等限定，判斷是否繼續(xù)迭代。如果是，進(jìn)入步驟13 ;反之進(jìn)入步驟17。步驟15在顆粒位移達(dá)到一定量時(shí)才進(jìn)入步驟11重算最優(yōu)鄰近顆粒矩陣，再進(jìn)行迭代計(jì)算；而步驟16則可直接進(jìn)入步驟13進(jìn)行迭代計(jì)算，而無(wú)需運(yùn)行搜索算法。由于巖土體數(shù)值模擬中顆粒相對(duì)靜態(tài)，而dS—般取0.1?0.5R(顆粒半徑)。通常數(shù)萬(wàn)次以上迭代才需運(yùn)行一次搜索算法；
[0016]步驟17返回?cái)?shù)值模擬結(jié)果。
[0017]離散元亦是指每個(gè)顆粒位移矩陣數(shù)據(jù)是離散的。
[0018]迭代方法是利用計(jì)算機(jī)不斷用變量的舊值遞推新值的過(guò)程。
[0019]通過(guò)空間網(wǎng)格計(jì)算初步鄰近顆粒矩陣的具體實(shí)現(xiàn)過(guò)程；
[0020]步驟20輸入顆粒坐標(biāo)P [X，Y, Z]和顆粒半徑數(shù)組R、采用m*l矩陣；
[0021]步驟21根據(jù)顆粒坐標(biāo)數(shù)據(jù)，建立空間網(wǎng)格并包含中心和鄰近顆粒所有顆粒；網(wǎng)格立方單元的邊長(zhǎng)取:Wg = 2*Rmax+dS
[0022]其中Rmax為最大顆粒半徑，dS為窗口半徑增量，通常取0.1-0.5R(顆粒半徑)；
[0023]步驟22根據(jù)顆粒的坐標(biāo)，計(jì)算顆粒在三維網(wǎng)格中的位置索引矩陣[IX，IY, IZ]、采用大小m*3矩陣，即顆粒在xyz坐標(biāo)方向上的網(wǎng)格單元序號(hào)；
[0024]步驟23將網(wǎng)格內(nèi)的顆粒編號(hào)記錄在對(duì)應(yīng)的四維矩陣中，即第一至第三維分別對(duì)應(yīng)顆粒在xyz方向上的網(wǎng)格編號(hào)，第四維上記錄存儲(chǔ)顆粒編號(hào)；
[0025]步驟24找出每個(gè)顆粒所在的中心網(wǎng)格單元索引[ix，iy，iz]，將鄰近網(wǎng)格中的顆粒定義為初步鄰近顆粒。從網(wǎng)格四維矩陣中將鄰近網(wǎng)格內(nèi)顆粒編號(hào)合并，建立初步鄰近顆粒數(shù)組Pntl ；
[0026]窗口法計(jì)算鄰近顆粒，以顆粒為中心劃定特定大小窗口，并將窗口內(nèi)顆粒定義為最優(yōu)鄰近顆粒。附圖3a給出了二維情況下的窗口示意圖。暗色顆粒PO為中心顆粒，以PO為中心劃定半徑為Rw的窗口，窗口內(nèi)灰色顆粒為最優(yōu)鄰近顆粒。
[0027]步驟30輸入各顆粒初步鄰近顆粒數(shù)組PnQ，以及顆粒坐標(biāo)P [X，Y, Z]和半徑R。
[0028]步驟31以每個(gè)顆粒為中心建立窗口，窗口半徑(Rw)為:
[0029]Rw = RfdS+R.其中Rtl為中心顆粒半徑；dS為窗口半徑增量；Rmax為最大顆粒半徑。此公式保證窗口外所有顆粒至少須在坐標(biāo)軸方向上運(yùn)動(dòng)dS距離才可能與中心顆粒接觸。
[0030]步驟32計(jì)算xyz坐標(biāo)方向上初步鄰近顆粒與中心顆粒的距離分量數(shù)組dPX, dPY,dPZ ；
[0031]步驟33將dPX，dPY，dPZ各距離分量均小于dS的顆粒定義為最優(yōu)鄰近顆粒，并存儲(chǔ)于最優(yōu)鄰近顆粒數(shù)組Pu;
[0032]步驟34重置顆粒位移矩陣D[X，Y，Z]為O ;在進(jìn)一步離散元迭代運(yùn)算中即可使用最優(yōu)鄰近顆粒數(shù)組進(jìn)行受力和運(yùn)動(dòng)計(jì)算。
[0033]每一個(gè)顆粒計(jì)算其鄰近顆粒時(shí)，均稱為中心顆粒。特定單元邊長(zhǎng)的設(shè)定與顆粒的粒徑有關(guān)(1-2倍顆粒直徑的長(zhǎng)度)。
[0034]本發(fā)明的有益效果是:結(jié)合使用窗口搜索算法和位移累積矩陣，一方面窗口搜索算法可以迅速得到最小化的鄰近顆粒數(shù)組，加快了迭代計(jì)算效率；另一方面通過(guò)位移累積矩陣控制搜索算法，極大地減少了搜索算法的運(yùn)行頻率。本方法避免了傳統(tǒng)離散元計(jì)算方法搜索算法效率低，且需要在迭代運(yùn)算中不斷運(yùn)行等缺點(diǎn)，而是通過(guò)顆粒位移矩陣記錄顆粒在各坐標(biāo)方向上位移，當(dāng)其位移超過(guò)一定值時(shí)運(yùn)行窗口搜索算法，以減少搜索算法運(yùn)行頻率。對(duì)巖土體等非流體模擬，使用位移累積矩陣控制窗口搜索算法，通常數(shù)萬(wàn)次迭代才需要運(yùn)行一次搜索算法，顯著提高了離散元數(shù)值計(jì)算的速度，可實(shí)現(xiàn)巖土體大規(guī)模顆粒系統(tǒng)的動(dòng)態(tài)模擬。

【專利附圖】

【附圖說(shuō)明】
[0035]圖1顆粒搜索和運(yùn)行方法總體結(jié)構(gòu)和工作原理流程圖；
[0036]圖2劃分空間網(wǎng)格計(jì)算初步鄰近顆粒流程圖；
[0037]圖3窗口法計(jì)算最優(yōu)鄰近顆粒示意圖。

【具體實(shí)施方式】
[0038]圖1為本方法實(shí)施的總體結(jié)構(gòu)和工作原理，其特點(diǎn)是應(yīng)用窗口法得到最優(yōu)鄰近顆粒矩陣，并通過(guò)位移矩陣來(lái)控制窗口搜索算法的運(yùn)行。
[0039]步驟10構(gòu)建初始的顆粒離散元堆積模型，顆粒總數(shù)為m，顆粒由I開始連續(xù)編號(hào)到m0
[0040]步驟11通過(guò)空間網(wǎng)格劃分計(jì)算得到初步鄰近顆粒。相應(yīng)附圖1a為二維情況下的網(wǎng)格示意圖，網(wǎng)格單元的邊長(zhǎng)為wg，所有顆粒均包含在網(wǎng)格中。以暗色顆粒PO為中心顆粒，黑色虛線范圍為鄰近網(wǎng)格，灰色顆粒為PO的初步鄰近顆粒，共11個(gè)。
[0041]步驟12通過(guò)建立顆粒窗口確定最優(yōu)鄰近顆粒。相應(yīng)附圖1b為二維情況下的窗口示意圖，窗口內(nèi)的顆粒可能不與中心顆粒接觸(如P5)，而窗口外的顆粒需要一個(gè)方向上(如X方向)運(yùn)動(dòng)至少dS距離，才可能與中心顆粒接觸(如P6和P7)。圖中灰色的最優(yōu)鄰近顆粒共有5個(gè)，通過(guò)窗口搜索有效減少了可能接觸顆粒數(shù)目。獲得新的最優(yōu)鄰近顆粒數(shù)組后，將顆粒位移矩陣D [X，Y, Z]全部數(shù)值設(shè)為零，重新開始累積顆粒位移。
[0042]步驟13利用最優(yōu)鄰近顆粒數(shù)組完成離散元數(shù)值計(jì)算。在每次迭代中，僅需判斷中心顆粒與最優(yōu)鄰近顆粒矩陣中的顆粒接觸關(guān)系，并計(jì)算相互作用力，因此計(jì)算量大大減少。一次迭代計(jì)算完成后，得到顆粒的位移量d[x，Y, Z]。
[0043]步驟14累積顆粒位移矩陣D [X，Y, Z] =D [X，Y, Z] +d [X，Y, Z]。通過(guò)累積每次迭代中的顆粒位移，位移矩陣D記錄了從前一次最優(yōu)鄰近顆粒搜索后的顆粒位移。
[0044]步驟15判斷111&1|0[乂，￥，2]|是否大于(15/2。即判斷顆粒位移在某一坐標(biāo)軸方向上是否超過(guò)dS/2。如果其大于dS/2，返回步驟11重新計(jì)算最優(yōu)鄰近顆粒；反之，進(jìn)入步驟
16。由于中心顆粒和鄰近顆粒均可運(yùn)動(dòng)，所以最大位移取dS的一半。
[0045]步驟16判斷是否繼續(xù)迭代。如果是，進(jìn)入步驟13 ;反之進(jìn)入步驟17。步驟15和16的區(qū)別在于:步驟15在顆粒位移達(dá)到一定量時(shí)才進(jìn)入步驟11重算最優(yōu)鄰近顆粒矩陣，再進(jìn)行迭代計(jì)算；而步驟16則可直接進(jìn)入步驟13進(jìn)行迭代計(jì)算，而無(wú)需運(yùn)行搜索算法。由于巖土體數(shù)值模擬中顆粒相對(duì)靜態(tài)，而dS —般取0.1?0.5R(顆粒半徑)。通常數(shù)萬(wàn)次以上迭代才需運(yùn)行一次搜索算法，因此最大限度地減少了搜索算法的運(yùn)行時(shí)間。
[0046]步驟17返回?cái)?shù)值模擬結(jié)果。
[0047]圖2表示通過(guò)空間網(wǎng)格計(jì)算初步鄰近顆粒矩陣的具體實(shí)現(xiàn)過(guò)程(步驟11的細(xì)化)，其特點(diǎn)是網(wǎng)格邊長(zhǎng)取特定值，以保證進(jìn)一步的窗口法可以正確有效運(yùn)行，以及通過(guò)多維矩陣來(lái)迅速完成計(jì)算。
[0048]步驟20輸入顆粒坐標(biāo)P [X，Y, Z]和顆粒半徑數(shù)組R (m*l矩陣)。
[0049]步驟21根據(jù)顆粒坐標(biāo)數(shù)據(jù)，建立空間網(wǎng)格并包含所有顆粒。網(wǎng)格立方單元的邊長(zhǎng)取:
[0050]Wg = 2*R眶+dS
[0051]其中Rmax為最大顆粒半徑，dS為窗口半徑增量，通常取0.1-0.5R(顆粒半徑)。此式可保證所有的窗口顆粒(最優(yōu)鄰近顆粒)都包含在初步鄰近顆粒矩陣中。附圖2a給出了二維情況下的網(wǎng)格劃分示意圖，暗色顆粒為中心顆粒。
[0052]步驟22根據(jù)顆粒的坐標(biāo)，計(jì)算顆粒在三維網(wǎng)格中的位置索引矩陣[IX，IY, IZ](大小m*3)，即顆粒在xyz坐標(biāo)方向上的網(wǎng)格單元序號(hào)。例如，在x方向的索引(IX)可以通過(guò)下式計(jì)算得到:
[0053]IX = floor ((PX-min (PX)) /Wg) +1 ； (Matlab 矩陣指令)
[0054]其中PX為顆粒X坐標(biāo)數(shù)組。
[0055]步驟23將網(wǎng)格內(nèi)的顆粒編號(hào)記錄在對(duì)應(yīng)的四維矩陣中，即第一至第三維分別對(duì)應(yīng)顆粒在xyz方向上的網(wǎng)格編號(hào)，第四維上記錄存儲(chǔ)顆粒編號(hào)。如附圖2a 二維網(wǎng)格中，網(wǎng)格[ix, iy] = [2，4]中包含Pl和P8兩個(gè)顆粒。
[0056]步驟24找出每個(gè)顆粒所在的中心網(wǎng)格單元索引[ix，iy，iz]，將鄰近網(wǎng)格中的顆粒定義為初步鄰近顆粒。從網(wǎng)格四維矩陣中將鄰近網(wǎng)格內(nèi)顆粒編號(hào)合并，建立初步鄰近顆粒數(shù)組Pu。。圖2a中虛線范圍內(nèi)為鄰近網(wǎng)格，其包含9個(gè)網(wǎng)格單元，網(wǎng)格內(nèi)灰色顆粒為鄰近顆粒。在三維模型情況下，鄰近網(wǎng)格包含27個(gè)立方單元。
[0057]圖3為窗口法計(jì)算鄰近顆粒示意圖(步驟12的細(xì)化)，其特點(diǎn)是以顆粒為中心劃定特定大小窗口，并將窗口內(nèi)顆粒定義為最優(yōu)鄰近顆粒。附圖3a給出了二維情況下的窗口不意圖。暗色顆粒PO為中心顆粒，以PO為中心劃定半徑為Rw的窗口，窗口內(nèi)灰色顆粒為最優(yōu)鄰近顆粒。
[0058]步驟30輸入各顆粒初步鄰近顆粒數(shù)組PnQ，以及顆粒坐標(biāo)P [X，Y, Z]和半徑R。
[0059]步驟31以每個(gè)顆粒為中心建立窗口，窗口半徑(Rw)為:
[0060]Rw = R0+dS+Rmax
[0061]其中Rtl為中心顆粒半徑；dS為窗口半徑增量；Rmax為最大顆粒半徑。此公式保證窗口外所有顆粒至少須在坐標(biāo)軸方向上運(yùn)動(dòng)dS距離才可能與中心顆粒接觸。
[0062]步驟32計(jì)算xyz坐標(biāo)方向上初步鄰近顆粒與中心顆粒的距離分量數(shù)組dPX, dPY,dPZ。例如X方向上的距離分量dPX為:
[0063]dPX = PX (Pn0) -Px0 ； (Matlab 矩陣指令)
[0064]其中PX(m*l矩陣)為顆粒x坐標(biāo)；公式右邊第一項(xiàng)得到初步鄰近顆粒x坐標(biāo)數(shù)值^xci為中心顆粒坐標(biāo)。
[0065]步驟33將dPX，dPY，dPZ各距離分量均小于dS的顆粒定義為最優(yōu)鄰近顆粒，并存儲(chǔ)于最優(yōu)鄰近顆粒數(shù)組Pn。
[0066]步驟34重置顆粒位移矩陣D[X，Y, Z]為O。在進(jìn)一步離散元迭代運(yùn)算中即可使用最優(yōu)鄰近顆粒數(shù)組進(jìn)行受力和運(yùn)動(dòng)計(jì)算，并且可以通過(guò)顆粒位移矩陣來(lái)控制和減少顆粒搜索算法的運(yùn)行，提高計(jì)算效率。
【權(quán)利要求】
1.一種高效率巖土體離散元數(shù)值計(jì)算顆粒搜索和運(yùn)行方法，其步驟包括: (1)建立初步鄰近顆粒數(shù)組；依據(jù)顆粒坐標(biāo)，建立具有特定單元邊長(zhǎng)的空間立方單元網(wǎng)格，以覆蓋所有顆粒，并建立相應(yīng)的網(wǎng)格編號(hào)矩陣；計(jì)算顆粒在網(wǎng)格位置，將每個(gè)顆粒的鄰近網(wǎng)格中的顆粒存儲(chǔ)為初步鄰近顆粒數(shù)組，每個(gè)顆粒相對(duì)于其鄰近顆粒為中心顆粒； (2)采用窗口法計(jì)算最優(yōu)鄰近顆粒；以每個(gè)顆粒為中心(即中心顆粒)建立特定大小的窗口，搜索初步鄰近顆粒位于窗口內(nèi)的顆粒，并存儲(chǔ)于最優(yōu)鄰近顆粒數(shù)組；而窗口外顆粒至少沿坐標(biāo)軸方向移動(dòng)dS距離才能與中心顆粒接觸； (3)利用位移累積矩陣控制搜索算法；在離散元迭代運(yùn)算中，累積計(jì)算每個(gè)顆粒位移矩陣D[X，Y, Z]，當(dāng)D某一位移分量超過(guò)dS/2時(shí)，意味著可能有新顆粒與中心顆粒接觸，運(yùn)行步驟1-2重新計(jì)算最優(yōu)鄰近顆粒數(shù)組，并重置顆粒位移矩陣D為O。
2.根據(jù)權(quán)利要求1所述的高效率巖土體離散元數(shù)值計(jì)算顆粒搜索和運(yùn)行方法，其特征是步驟(2)中應(yīng)用窗口法計(jì)算最優(yōu)鄰近顆粒，得到最優(yōu)鄰近顆粒數(shù)組，并通過(guò)位移矩陣來(lái)控制窗口搜索算法的運(yùn)行； 步驟10構(gòu)建初始的顆粒離散元堆積模型，顆粒總數(shù)為m，顆粒由I開始連續(xù)編號(hào)到m ；步驟11通過(guò)空間網(wǎng)格劃分計(jì)算得到初步鄰近顆粒；網(wǎng)格單元的邊長(zhǎng)為Wg，所有顆粒均包含在網(wǎng)格中； 步驟12通過(guò)建立半徑為Rw顆粒窗口確定最優(yōu)鄰近顆粒；由此獲得顆粒窗口內(nèi)最優(yōu)鄰近顆粒數(shù)組，將顆粒位移矩陣D [X，Y，Z]全部數(shù)值設(shè)為零，重新開始累積顆粒位移； 步驟13利用最優(yōu)鄰近顆粒數(shù)組完成離散元數(shù)值計(jì)算；在每次迭代中，僅需判斷中心顆粒與最優(yōu)鄰近顆粒矩陣中的顆粒接觸關(guān)系，并計(jì)算相互作用力；一次迭代計(jì)算完成后，得到中心顆粒的位移量d[X，Y，Z]； 步驟14累積中心顆粒位移矩陣D[X，Y, Z] = D[X, Y, Z]+d[X, Y, Z];通過(guò)累積每次迭代中的中心顆粒位移，位移矩陣D記錄了從前一次最優(yōu)鄰近顆粒搜索后的顆粒位移； 步驟15判斷max |D [X，Y，Z] |是否大于dS/2 ;即判斷顆粒位移在某一坐標(biāo)軸方向上是否超過(guò)dS/2。如果其大于dS/2，返回步驟11重新計(jì)算最優(yōu)鄰近顆粒；反之，進(jìn)入步驟16 ；步驟16根據(jù)迭代次數(shù)等限定，判斷是否繼續(xù)迭代。如果是，進(jìn)入步驟13 ;反之進(jìn)入步驟17。步驟15在顆粒位移達(dá)到一定量時(shí)才進(jìn)入步驟11重算最優(yōu)鄰近顆粒矩陣，再進(jìn)行迭代計(jì)算；而步驟16則可直接進(jìn)入步驟13進(jìn)行迭代計(jì)算，而無(wú)需運(yùn)行搜索算法。由于巖土體數(shù)值模擬中顆粒相對(duì)靜態(tài)，而dS —般取0.1?0.5R、R為顆粒半徑；通常數(shù)萬(wàn)次以上迭代才需運(yùn)行一次搜索算法； 步驟17返回?cái)?shù)值模擬結(jié)果。
3.根據(jù)權(quán)利要求1所述的高效率巖土體離散元數(shù)值計(jì)算顆粒搜索和運(yùn)行方法，其特征是通過(guò)空間網(wǎng)格計(jì)算初步鄰近顆粒矩陣的具體實(shí)現(xiàn)過(guò)程； 步驟20輸入顆粒坐標(biāo)P [X，Y, Z]和顆粒半徑數(shù)組R、采用m*l矩陣； 步驟21根據(jù)顆粒坐標(biāo)數(shù)據(jù)，建立空間網(wǎng)格并包含中心和鄰近顆粒所有顆粒；網(wǎng)格立方單元的邊長(zhǎng)取:
Wg = 2*Rfflax+dS 其中R_為最大顆粒半徑，dS為窗口半徑增量，通常取0.1-0.5R(顆粒半徑)； 步驟22根據(jù)顆粒的坐標(biāo)，計(jì)算顆粒在三維網(wǎng)格中的位置索引矩陣[IX，IY, IZ]、采用大小m*3矩陣，即顆粒在xyz坐標(biāo)方向上的網(wǎng)格單元序號(hào)； 步驟23將網(wǎng)格內(nèi)的顆粒編號(hào)記錄在對(duì)應(yīng)的四維矩陣中，即第一至第三維分別對(duì)應(yīng)顆粒在xyz方向上的網(wǎng)格編號(hào)，第四維上記錄存儲(chǔ)顆粒編號(hào)； 步驟24找出每個(gè)顆粒所在的中心網(wǎng)格單元索引[ix，iy，iz]，將鄰近網(wǎng)格中的顆粒定義為初步鄰近顆粒。從網(wǎng)格四維矩陣中將鄰近網(wǎng)格內(nèi)顆粒編號(hào)合并，建立初步鄰近顆粒數(shù)組 Pn0 ； 窗口法計(jì)算鄰近顆粒，以顆粒為中心劃定特定大小窗口，并將窗口內(nèi)顆粒定義為最優(yōu)鄰近顆粒。附圖3a給出了二維情況下的窗口示意圖。暗色顆粒PO為中心顆粒，以PO為中心劃定半徑為Rw的窗口，窗口內(nèi)灰色顆粒為最優(yōu)鄰近顆粒。 步驟30輸入各顆粒初步鄰近顆粒數(shù)組Prv以及顆粒坐標(biāo)P [X，Y, Z]和半徑R。 步驟31以每個(gè)顆粒為中心建立窗口，窗口半徑(Rw)為:
Rw = R0+dS+Rmax 其中Rtl為中心顆粒半徑；dS為窗口半徑增量；Rmax為最大顆粒半徑。此公式保證窗口外所有顆粒至少須在坐標(biāo)軸方向上運(yùn)動(dòng)dS距離才可能與中心顆粒接觸。 步驟32計(jì)算xyz坐標(biāo)方向上初步鄰近顆粒與中心顆粒的距離分量數(shù)組dPX,dPY,dPZ ；步驟33將dPX，dPY，dPZ各距離分量均小于dS的顆粒定義為最優(yōu)鄰近顆粒，并存儲(chǔ)于最優(yōu)鄰近顆粒數(shù)組Pn ； 步驟34重置顆粒位移矩陣D [X，Y，Z]為O ;在進(jìn)一步離散元迭代運(yùn)算中即可使用最優(yōu)鄰近顆粒數(shù)組進(jìn)行受力和運(yùn)動(dòng)計(jì)算。
【文檔編號(hào)】G06F17/30GK104239414SQ201410405178
【公開日】2014年12月24日 申請(qǐng)日期:2014年8月15日 優(yōu)先權(quán)日:2014年8月15日
【發(fā)明者】劉春 , 施斌, 王寶軍, 張丹, 索文斌, 梁修雨, 張馳 申請(qǐng)人:南京大學(xué)