1、PFC2D (Particle Follow Code 2 Dimen)即二維顆粒流程序，是通過離散單元方法來模擬圓形顆粒介質(zhì)的運動及其相互作用。最初，這種方法是顆粒介質(zhì)特性的一種工具，它采用數(shù)值方法將物體分為有代表性的數(shù)百個顆粒單元，期望利用這種局部的模擬結(jié)果來邊值間題連續(xù)計算的本構(gòu)模型。以下兩種促使 PFC2D 方法產(chǎn)生與發(fā)展:(1)通過現(xiàn)場實驗來得到顆粒介質(zhì)本構(gòu)模型相當(dāng):(2)隨著微機功能的逐步增強，用顆粒模型模擬整個問題成為可能，一些本構(gòu)特性可以在模型中自動形成。因此，PFC2D 便成為用來模擬固體力學(xué)和顆粒流問題的一種有效。2、顆粒流方法的基本假設(shè) ：顆粒流方法在模擬過程中作了如下假

2、設(shè): 1)顆粒單元為剛性體;2)接觸發(fā)生在很小的范圍內(nèi)，即點接觸;3)接觸特性為柔性接觸，接觸處允許有一定的“”量;4) “”量的大小與接觸力有關(guān)，與顆粒大小相比，“”量很小;接觸處有特殊的連接強度;顆粒單元為圓盤形(或球形)。3、顆粒流方法的特點：PFC2D 可以直接模擬圓形顆粒的運動和相互作用問題。顆料可以代表材料中的個別顆粒，例如砂粒，也可以代表粘結(jié)在一起的固體材料，例如混凝土或巖石。當(dāng)粘結(jié)以漸進(jìn)的方式破壞時，它能夠破裂。粘結(jié)在一起的集合體可以是各向同性，也可以被分成一些離散的區(qū)域或塊體。這類物理系統(tǒng)可以用處理角狀塊體的離散單元程序 UDEC 和 3DEC 來模擬。 PFC2D 有三個優(yōu)

3、點：第一、它有潛在的高效率。因為圓形物體間的接觸探測比角狀物體間的更簡單。第二、對可以模擬的位移大小實質(zhì)上沒有限制。第三、由于它們是由粘結(jié)的粒子組成，塊體可以破裂，不象 UDEC 和 3DEC 模擬的塊體不能破裂。用 PFC2D 模擬塊體化系統(tǒng)的缺點是，塊體的邊界不是平的，用戶必須接受不平的邊界以換取 PFC2D 提供的優(yōu)點。PFC2D 能模擬任意大小圓形粒子集合體的動態(tài)力學(xué)行為。粒子分布規(guī)律分布。根據(jù)粒子的指定分布規(guī)律自動概率地生成。粒子半徑按均勻分布或按初始孔隙度一般比較高，但通過控制粒子半徑的擴(kuò)大可以獲得密度壓實。在任何階段任何都可以改變半徑。所以不需反復(fù)試驗就可以獲得指定孔隙度的壓實狀

4、態(tài)。屬性與各個粒子或接觸有關(guān)，而不是與“類型號”有關(guān)。因此，可以指定屬性和半徑的連續(xù)變化梯度?！肮?jié)理”用來修改沿指定軌跡線的接觸特性。假定這些線疊加在顆粒集合體上。用這種方法，模型可以被成組的弱面，如巖石節(jié)理切割。粒子顏色也是一種屬性，用戶可以指定各種標(biāo)記方案。PFC2D 模型中為了保證數(shù)據(jù)長期不漂移，精度數(shù)據(jù)坐標(biāo)和半徑。接觸的相對位移直接根據(jù)坐標(biāo)而不是位移增量計算。接觸性質(zhì)由下列單元組成：1）線性彈簧或簡化的 Hertz-Mindlin 準(zhǔn)則；2）滑塊；3）粘結(jié)類型：粘結(jié)接觸可承受拉力，粘結(jié)存在有限的抗拉和抗剪強度?？稍O(shè)定兩種類型的粘結(jié)，接觸粘結(jié)和平行粘結(jié)。這兩種類型粘結(jié)對應(yīng)兩種可能的物理接

5、觸：接觸粘結(jié)再現(xiàn)了作用在接觸點一個很小區(qū)域上的附著作用；平行粘結(jié)再現(xiàn)了粒子接觸后澆注其它材料的作用（剛度。泥灌漿）。平行粘結(jié)中附加材料的有效剛度具有接觸點的塊體邏輯支持附屬粒子組或塊體的創(chuàng)建，促進(jìn)了程序的推廣普及。塊體內(nèi)粒子可以任意程度的，作為剛性體具有可變形邊界的每一個塊體，可作為一般形狀的超級粒子。通過指定墻的速度、混合的粒子速度、施加外力和重力來給系統(tǒng)加載?！皵U(kuò)展的 FISH 庫”提供了在集合體內(nèi)設(shè)置指定應(yīng)力場或施加應(yīng)力邊界條件的函數(shù)。時步計算是自動的，包括因為 Hertz接觸模型剛度變化的影響。模擬過程中，根據(jù)每個粒子周圍接觸數(shù)目和瞬間剛度值，時步也在變化?；诠烙嫷牧Ｗ訑?shù)，單元來適應(yīng)

6、粒子缺失和指定的新對象。單元目線性增加，而不是二次方增加。策略采用最佳的單元數(shù)目，自動調(diào)整單元的外部尺寸方案支持接觸探測算法以保證求解時間隨粒子數(shù)類似于 FLAC，PFC 提供了局部無粘性阻尼。這種阻尼形式有以下優(yōu)點：對于勻速運動，體力接近于零，只有阻尼系數(shù)是無因次的；運動時才有阻尼；3）因阻尼系數(shù)不隨頻率變化，集合體中具有不同自然周期的區(qū)域被同等阻尼，采用同樣的阻尼系數(shù)。PFC2D 可以在半靜態(tài)模式下運行以保證迅速收斂到靜態(tài)解，或者在完全動態(tài)模式下運行。PFC2D 包含功能強大的內(nèi)嵌式程序語言 FISH，允許用戶定義新的變量和函數(shù)使數(shù)值模型適合用戶的特殊需求。例如，用戶可以定義特殊材料的模型

7、和性質(zhì)、加載方式、實驗條件的伺服控制、模擬的順序以及繪圖和打印用戶定義的變量等。2.FISH 語言簡介FISH 是一種內(nèi)置于 Itasca內(nèi)的編程語言，使用 FISH用戶可以定義新的變量和函數(shù)，從而使得這些函數(shù)被用來擴(kuò)展 Itasca的用法或者增加用戶自定義的特性。例如，用戶可以繪制(PLOT)或打印(PR)新的變量，也能夠改進(jìn)特殊的顆粒體模型(網(wǎng)格)，可以對數(shù)值試驗進(jìn)行伺服控制，可以設(shè)定一些性質(zhì)的特殊分布以及自動化參數(shù)。Itasca 已經(jīng)寫了一些簡單但非常有用的 FISH 函數(shù)作為庫文件包含在各個具體的中，這些 FISH 函數(shù)一方面方便了一些沒有編程經(jīng)驗的用戶寫一些簡單的 FISH 函數(shù)，另

8、一方面也能使用戶在這些提供的簡單函數(shù)的基礎(chǔ)上作進(jìn)一步的改進(jìn)。不過，F(xiàn)ISH 語言象其它編程語言一樣，也可以編制出非常復(fù)雜的程序利用 FISH 語言進(jìn)行編程，應(yīng)該首先編一些簡單的函數(shù)，然后仔細(xì)檢查函數(shù)的功能，測試是否有錯誤。如果沒有發(fā)現(xiàn)錯誤，再逐漸增加其功能，增加一項功能檢查一下，直至發(fā)展到最后比較復(fù)雜的程序。這是因為雖然 FISH 是一種編譯型語言，但它沒有自己獨立的編譯器，不象 VC+或 VB 能夠?qū)崟r全面地檢查錯誤，F(xiàn)ISH 檢查錯誤的能力很差，因此在使用他們到真實的應(yīng)用之前，一定要用一些簡單的數(shù)據(jù)（假如可能的話）來檢查所有定義的函數(shù)。FISH 函數(shù)內(nèi)置于標(biāo)準(zhǔn)的 Itasca的數(shù)據(jù)文件中，

9、函數(shù)的格式必須以 DEFINE 開始，以END 結(jié)束。函 數(shù)可以嵌套調(diào)用，但定義函數(shù)的次序沒有關(guān)系，只要在使用之前全部定義就行。由于 FISH 函數(shù)的編譯格式數(shù)以及相關(guān)變量的當(dāng)前值。在 Itasca的內(nèi)存中，因此可以用 SAVE 命令保存函FISH 也可以用來改進(jìn)用戶寫的本構(gòu)模型，如例 1：DEF abcabc=22*3+5ENDPr對上例子稍作改進(jìn)(例 2)：abcnewdef abchh=22abc=hh*3+5end的人可以看出，執(zhí)行上面的例子（PR1).稍有編程abc），其結(jié)果與例 1 相同：abc=71.在這個函數(shù)中，首先把 22 賦值給變量 hh，然后把這個變量帶入 abc 的表達(dá)

10、式中，因此二者的結(jié)果相同。2).FISH 的執(zhí)行過程如下：當(dāng)在程序命令中使用一個 FISH 符號名時（例如執(zhí)行 PR符號名），如果符號名也是一個函數(shù)名，那么執(zhí)行這個函數(shù)（例如 abc）；如果符號名不是函數(shù)名，那么使用符號目前的值(例如 hh)。3).在輸入完例 2 的各行后，如果執(zhí)行命令：PRhh，此時 hh=0,因為在這個時候沒有執(zhí)行 FISH 函數(shù)，因此 hh 的初始值為 0；接著執(zhí)行 PRabc，結(jié)果顯示 abc=71;再次執(zhí)行 PR被賦值。hh,此時結(jié)果為 hh=22，這是因為首先運行了 abc 函數(shù)，在這個過程中 hh 已4). 下面的試驗將進(jìn)一步解釋函數(shù)與變量之間的差別。注意：It

11、asca的 SET 命令可以用來設(shè)置任何用戶定義的 FISH 符號的值，與在 FISH 中使用的符號無關(guān)。下面的例 3 建立在例 2 的基礎(chǔ)之上，不使用 NEW 命令來清除內(nèi)存中的值，因為想繼續(xù)使用那些值:set abc=0 hh=0pr prprhh abchh5).在這個例子中，首先把 abc 和 hh 都賦值為 0，由于 hh 是一個變量，第一個 Pr命令顯示當(dāng)前 hh 的值，hh=0;第二個 Pr命令由于 abc 是一個函數(shù)名，因此執(zhí)行 abc 函數(shù)，先前定義的 abc=0 不起作用，重新計算了 hh 和 abc 的值，因此第三個 Pr命令顯示的值是它在 abc 函數(shù)內(nèi)指定的值，即 h

12、h=22,例 4 是這個試驗完整令。象其它高級編程語言一樣，F(xiàn)ISH 有執(zhí)行循環(huán)命令的功能，標(biāo)準(zhǔn)的格式如下：LOOPvar(expr1,expr2)END_LOOP其中 LOOP 和 END_LOOP 是 FISH 語句，符號 var 代表循環(huán)變量，expr1 和 pxpr2 代表表達(dá)式或者單個變量，下面的例子用循環(huán)命令計算從 1 到 10 的和以及乘積,見下例：new def abcsum= 0prod = 1loop n (1,10)sum= sum + n prod = prod * nend_loopend abcprsum, prod在這個例子中，首先給兩個變量賦于初始值，sum 用

13、來保存和的結(jié)果，prod 用來保存積的結(jié)果，然后執(zhí)行循環(huán)，最后分別打印出這兩個變量的最后結(jié)果。循環(huán)變量 n(1,10)表示從 1 開始，連續(xù)計算到 10 結(jié)束。關(guān)于 LOOP 的注意事項：1).FISH 接受 END_LOOP 和 ENDLOOP 的寫法，但不接受 END LOOP 這樣中間有空格的寫法，其它類似令有著同樣的規(guī)則，如 END_IF，MAND 等命令；2).在上面的例子中，如果執(zhí)行 Prn 或 Prfish 命令，你會看到 n=11 而不是 10，注意：這不是 FISH 的錯誤，這是一個基本的計算機指令規(guī)則，當(dāng)循環(huán)結(jié)束后，計數(shù)器的值保存的是 n+1 而不是 n,所有的高級編程語言

14、有著相同的規(guī)則。 1.DEFINEfunctionENDCASEOFexpr Case nendcaseIFexpr1 test expr2 THEN ELSEENDIF4. LOOPvar (expr1, expr2)ENDLOOP5. LOOP WHILE expr1 test expr2 ENDLOOP6.DMAND7. HISTORY varPRvarSET var value PLOT add .sh fname另外，在 FISH 中還有許多其它的預(yù)定義對象,其中一類是尺度變量(scalar variables),它們是單個的數(shù)字,下面是總的尺度變量:clock時鐘時間,unbal最

15、大不平衡力pi圓周率是秒的 100 倍.step目前的時步數(shù)目urand0.0-1.0 之間均勻分布的隨量這僅是其中的一小部分,完全的列表以后再述另一類非常有用的內(nèi)置對象是固有函數(shù)(rinsic functions),這些函數(shù)能在 FISH 內(nèi)進(jìn)行一些比較高級的數(shù)算,完整的列表見 FISH 手冊,下面給出其中的一部分:abs(a)a 的絕對值cos(a)a 的(a 為弧度)log(a)a 的底數(shù)為 10 的對數(shù)max(a,b)返回 a,b 中的最大值sqrt(a)a 的平方根3.PFC2D 計算模型的生成方法有兩個命令可用于生成顆粒流模型:BALL 和 GENER-ATE，其中，BALL 命令

16、是生成單個的顆粒，該命令生成的顆?？膳c已存在的顆粒，而 GENERATE 可生成一系列指定數(shù)目的顆粒流，該命令生成的顆粒是不允許的。PFC2D 里主要有兩種類型的顆粒流：規(guī)則排列的和無規(guī)則排列的。一系列規(guī)則排列的顆粒流可以用來模擬模擬結(jié)構(gòu)部分，如梁，而不規(guī)則排列的顆粒流可用來模擬實體或內(nèi)部結(jié)構(gòu)無規(guī)則的顆粒材料，如巖石內(nèi)部所包含的膠結(jié)顆粒。盡管顆粒的排列是隨機的，但在顆粒模型生成后，整個模型的結(jié)構(gòu)特性還是可能會受影響的，比如弱的結(jié)構(gòu)面或各向異性。對于無規(guī)律排列的顆粒流模型，一般不可能去描述它的初始接觸力的量級大小，這必須在后期要經(jīng)過一個壓縮的過程才可能給予較好的評價。3.1 規(guī)律排列顆粒流：生成

17、規(guī)則排列的顆粒流，主要采用 FISH 語言配合 BALL 命令，循環(huán)生成一系列的顆粒，如下例:Newdef hexxc = x0 yc = y0rc = radius idc = id_startr2 = 2.0 * radiusyinc = radius * sqrt(3.0) loop row (1,n_row)loop col (1,n_col)dball id=idc x=xc y=yc rad=rc mandidc = idc + 1 xc = xc + r2 end_loopyc = yc + yincxc = x0 + radius * (row - (row/2) * 2) e

18、nd_loopendset echo off？set x0=0.2 y0=0.4 radius=0.1set id_start=100 n_col=7 n_row=8 hexset echo onplot set cap size 20 plot add axes black plot add ball yellow plot show3.2 不規(guī)則排列；無規(guī)則排列，即：對一個給定空隙率的區(qū)域，采用顆粒來充填其中需要進(jìn)行填充的空隙，并確保整個模型保持平衡。對于所能被填充的模型的初始空隙率，是有一個限制值，不能任意小。對于某些空隙率的模型，顆粒的填充可以無接觸地排列，對于其它情況的空隙率，顆粒又

19、可以排列目前沒有一個普遍的方法來將緊密充填的圓形顆粒體限制在一個任意封閉的面積內(nèi)。這里來介紹兩種方法如何充填特定半徑的顆粒體達(dá)到所需要的空隙率。第法，首先建立封閉區(qū)域的邊界(簡稱墻體)，然后在封閉區(qū)域內(nèi)任意生成一系列無接觸的顆粒，最后移動區(qū)域的限制墻體，至所需要的空隙率。這種方法有三個缺點：1.區(qū)域的幾何形狀改變；2.收斂速度慢；3.最終的分布趨勢是不均勻的第二種方法：運用 GENERATE 命令生成顆粒體，同時配合分配，即指定顆粒體半徑的上下限，然后相應(yīng)分配一個標(biāo)準(zhǔn)差，同時配合 FISH 函數(shù)來選擇顆粒半徑，最終生成所需要的模型。下面來介紹采用這種方法是如何進(jìn)行操作的，重點介紹 2 種方法。

20、3.2.1 半徑擴(kuò)展法首先需要介紹一個概念：半徑放大系數(shù) m，它會引起模型空隙率的改變。空隙率的定義：顆粒體的總面積； 封閉區(qū)域的總面積因此：定義初始空隙率為，新生成的空隙率為,那么，則有實例：區(qū)域?qū)挘?0區(qū)域高：5目標(biāo)空隙率：0.12顆粒體數(shù)目：300最大最小半徑比：1.5初始假設(shè)一個 m 之值，可以求出初始的為從而得到平均半徑大小為：由初始給定的顆粒體最大直徑和最小直徑的比例 r到顆粒體半徑的上下限：通過上面的計算，同時編寫 FISH 函數(shù)來生成滿足需要的空隙率的模型，并能自動計算空隙率的大小，以便用戶進(jìn)行檢驗。newdef expand;- 輸入數(shù)據(jù) - n_stiff = 1e8 s_

21、stiff = 1e8 width = 10.0height = 5.0; 法向連接剛度; 剪切連接剛度; 區(qū)域?qū)? 區(qū)域高tot_vol = width*heightporos = 0.12num = 300rat = 1.5導(dǎo)出所需數(shù)據(jù)mult = 1.6; 最終目標(biāo)空隙率; 顆粒體數(shù)目; 最大最小半徑比-; 初始半徑放大系數(shù);-n0 = 1.0 - (1.0 - poros) / mult2r0 = sqrt(height*width*(1.0 - n0)/(pi*num) rlo = 2.0 * r0 / (1.0 + rat)rhi = rat * rlo_x1 = width*(1

22、.0 + extend)_y1 = 0.0dwall id=1 ks=s_stiff kn=n_stiff nodes (_x0,_y0) (_x1,_y1) mand_x0 = width_y0 = -extend*height_x1 = width_y1 = height*(1.0 + extend) dwall id=2 ks=s_stiff kn=n_stiff nodes (_x0,_y0) (_x1,_y1) mand_x0 = width*(1.0 + extend)_y0 = height_x1 = -extend*width_y1 = height dwall id=3 ks

23、=s_stiff kn=n_stiff nodes (_x0,_y0) (_x1,_y1)mand_x0 = 0.0_y0 = height*(1.0 + extend)_x1 = 0.0_y1 = -extend*height dwall id=4 ks=s_stiff kn=n_stiff nodes (_x0,_y0) (_x1,_y1) mand;- generate the balls and give them their properties dgen id=1,num rad=rlo,rhi x=0,width y=0,height prop dens=1000 ks=s_st

24、iff kn=n_stiffmandget_porosmult = sqrt(1.0 - poros) / (1.0 - pmeas) dini rad mul mult cycle 1000prop fric 0.2cycle 2000mandenddef get_porossum = 0.0bp = ball_headloop whip # nullsum = sum + pi * b_rad(bp)2bp = b_next(bp) end_looppmeas = 1.0 - sum / (width * height)end expand get_poros plot wall ball

25、plot showprpmeassave expand.SAV3.2.2 擠壓排斥法：第 2 種方法：指定顆粒體的半徑，不限制顆粒的數(shù)目，使足夠多的顆粒產(chǎn)生來達(dá)到所需要的空隙率。但這種方法所 帶來的缺點是可能在局部區(qū)域造成大面積的顆粒，這將會產(chǎn)生很大的擠壓力，從而給予顆粒較大的初始速度，這就可能使得顆粒體脫離墻體的限制。為避免此情況的發(fā)生，可通過初始的有限步循環(huán)計算將顆粒的動能減至零，然后再計算至平衡態(tài)。newset random ; 設(shè)定顆粒體數(shù)目def exploden_stiff = 1e8 s_stiff = 1e8 width = 10.0height = 5.0poros = 0.

26、12n_max = 1000rlo = 0.174rhi = 0.261 d模式wall id 1 ks=s_stiff kn=n_stiff nodes (0,0) (width,0)wall id 2 ks=s_stiff kn=n_stiff nodes (width,0) (width,height) wall id 3 ks=s_stiff kn=n_stiff nodes (width, height) (0,height) wall id 4 ks=s_stiff kn=n_stiff nodes (0,height) (0,0)mand pvol_sum = 0.0tot_vo

27、l = width * height count = 0sectionloop n (1,n_max)r_ball = rlo + urand * (rhi - rlo) pvol_new = pvol_sum + pi * r_ball2if (1.0 - pvol_new / tot_vol) poros thenexit section ; (new porosity will be = id1 thenif b_id(bp) = id2 thensum = sum + pi * b_rad(bp)2 end_ifend_ifbp = b_next(bp) end_looppmeas =

28、 1.0 - sum / tot_volenddef final_porostot_vol = width * height id1 = 1id2 = 1200get_poros final_poros = pmeasendset x1=0.0 x2=5.0 y1=0.0 y2=5.0 id1=1 id2=50make_blockset mult_a=multset x1=5.0 x2=10.0 y1=0.0 y2=5.0 id1=1001 id2=1200make_blockset mult_b=multini rad mul=mult_a c_index 0 range id 1,50in

29、i rad mul=mult_b c_index 1 range id 1001,1200 plo create the_assemblyplot add ball lgreen lorange plot add wall blackplot show cycle 1000prop fric 0.2cycle 500prfinal_porossave expand2.SAV3.4 運用模型生成“過濾器”有些情況需要建立復(fù)雜區(qū)域形狀的顆粒流模型，如右圖，此時可運用模型生成過濾器來獲得所需要的模型，即 filter 命令，其后由用戶定義 FISH 函數(shù)來控制,其中,顆粒的半徑通過 fc_arg(0

30、)進(jìn)行檢驗,x 和y的坐標(biāo)位置分別通過 fc_arg(1) 和 fc_arg(2)進(jìn)行檢驗。如果顆粒滿足要求，則 FISH 函數(shù)值設(shè)為0，否則為 1。詳見下例。newde_rect;用戶定義生成過濾器生成方形環(huán)狀顆粒流模型; 中心 (ff_x, ff_y), 內(nèi)徑 ff_r1 and 外徑_brad = fc_arg(0)_bx = fc_arg(1)_by = fc_arg(2)_skip = 0_rx = abs( _bx - ff_x ) - _brad_ry = abs( _by - ff_y ) - _brad if _rx ff_r1 thenif _ry ff_r1 then_s

31、kip = 1 end_ifend_ifff_rect = _skipff_r2.enddef gen_balls_xlo = ff_x - ff_r2_xhi = ff_x + ff_r2_ylo = ff_y - ff_r2_yhi = ff_y + ff_r2dgenerate x=(_xlo, _xhi) y=(_ylo, _yhi) & rad=(0.09, 0.11) &filter=ff_rect & id=(1,250)mandendset ff_x=1.0 ff_y=1.0 ff_r1=2.0 ff_r2=3.0gen_ballsproperty dens=1000 kn=1

32、e8 ks=1e8wall id=1 nodes (-1.0,-1.0) (-1.0, 3.0) (3.0,3.0) (3.0,-1.0) closewall id=2 nodes (-2.0,-2.0) ( 4.0,-2.0)wall id=3 nodes (4.0,-2.0) ( 4.0, 4.0)wall id=4 nodes ( 4.0, 4.0) (-2.0, 4.0)wall id=5 nodes (-2.0, 4.0) (-2.0,-2.0) wall id=1 kn=1e8 ks=1e8wall id=2 kn=1e8 ks=1e8 wall id=3 kn=1e8 ks=1e

33、8wall id=4 kn=1e8 ks=1e8wall id=5 kn=1e8 ks=1e8 plot create the_viewplot add ball yellow plot add axes blackplot add wall blue id=on plot showpauseproperty rad mul 1.5 plot add cf green cycle 5004.邊界條件：PFC2D 中有三種邊界條件，分別是：墻體邊界、顆粒體邊界和混合邊界。其中顆粒體邊界又分為速度邊界和受力邊界。1).墻體邊界建模過程中，墻體可作為顆粒體的生成范圍約束，但同時也可以將墻體作為邊界來

34、施加約束。對于墻體，只能施度約束，而不能直接對其施加外力，因為運動定律對墻體是不適用的。其速度由以下三個參數(shù)控制:線速度、角速度和旋轉(zhuǎn)中心。墻的運動是通過不斷更新定義墻的基點的位置來描述。采用 WALL 命令設(shè)置，如：wall id=1 x=1.0 y=1.0 spin=10.0(2)顆粒體邊界PFC2D 中的模型可以將一連串的顆粒體作為邊界條件?；痉椒ǎ涸谀Ｐ途o密壓縮至平衡后，通過 FISH 函數(shù)將與墻體相接觸的顆粒體逐個提取，將這一系列的顆粒體采用共同的邊界條件限制，最后刪除初始的限制墻體，即實現(xiàn)了以顆粒體代替墻體來作為邊界條件。顆粒體邊界條件分速度邊界和外力邊界，下例為速度邊界程序?qū)崿F(xiàn)

35、restore expand.sav def boundbp = ball_headloop whisectionp # nullcp = b_clist(bp) loop while cp # nullif c_nforce(cp) # 0.0 then b2 = c_ball2(cp)if poer_type( b2 ) = 101 then ; b2 is a wall b_xfix(bp) = 1 ; fix original ball in x,y b_yfix(bp) = 1b_color(bp) = 1exit section ; all done for this ball e

36、nd_ifend_ifif c_ball1(cp) = bpcp = c_b1clist(cp) elsecp = c_b2clist(cp) end_ifend_loop end_sectionbp = b_next(bp) end_loopend boundini xvel=0.0 grad=(-0.1, 0) ini yvel=0.0del wall 1del wall 2del wall 3del wall 4plot create assembly plot add ball red green plot add vel blackplot showplot add cf white

37、 cycle 200顆粒體外力邊界：其前期操作與速度邊界類似，差異處是先將邊界顆粒體的速度初始化為零，再刪除限制墻體，采用 FISH 函數(shù)反向施加顆粒體所受的不平衡力，運行一定的計算步至平衡態(tài)。在此過程中，區(qū)域內(nèi)部的顆粒體會因輕微的擾動而導(dǎo)致其自身產(chǎn)生運動而脫離邊界顆粒體，此時前面將邊界顆粒體的速度固定為零。ini xvel=0 yvel=0 spin=0 del wall 1del wall 2del wall 3del wall 4cycle 1def replace ; Replace unbalanced forbp = ball_headwippd forloop whip # nu

38、llif b_xfix(bp) = 1b_xfap(bp) = -b_xfob(bp) b_yfap(bp) = -b_yfob(bp)end_ifbp = b_next(bp)end_loopend replacefree x y spincycle 200(3)混合邊界條件：混合邊界即速度邊界與外力邊界同時存在，對于雙軸壓縮試驗，其兩側(cè)端面不允許存在明顯的幾何變形，否則其邊界力將會無效。因此，此時需要用到混合邊界，兩側(cè)端面采用外力邊界來約束，確保其不發(fā)生幾何變形，側(cè)面 Y 向采用速度邊界，使其速度線性增加，以此來模擬試塊軸向變形。5.初始條件 : 首先介紹一下應(yīng)力的概念。在連續(xù)介質(zhì)中，應(yīng)力

39、是續(xù)量，應(yīng)力一般定義線段上的合力。然而，為作用于面積上的力，對于二維問題，應(yīng)力則定義為作用于對于松散介質(zhì)這種定義是不合適的，因為顆粒介質(zhì)是離散的，在顆粒模型中沒有存在于每一個點上連續(xù)的應(yīng)力，而且變化幅度很大。對于應(yīng)變也存在這樣的問題。在 PFC 模型中，計算顆粒間接觸力和顆粒的位移，這些量對于在微觀上材料的特性很有意義，但不能將這些量直接聯(lián)系到連續(xù)模型中，需要經(jīng)過一個平均的過程，才能從微觀傳到連續(xù)模型中。所以,在 PFC2D 中采用平均應(yīng)力的概念來表示5.1 施加等向應(yīng)力def expand_radii_sum = 0.0bp= ball_headloop whip # null;找到與 bp

40、 相連的下一個顆粒的物理地址cp = b_clist(bp) loop while cp # nullif c_nforce(cp) # 0.0 then_rcp = sqrt(c_x(cp) - b_x(bp)2 + (c_y(cp) - b_y(bp)2) if _rcp = b_rad(bp) then; 剛好接觸或if c_ball1(cp) = bp then bp_other = c_ball2(cp)elsebp_other = c_ball1(cp) end_ifif poer_type(bp_other) = 101 then;21, 100, 101, 102, 103 o

41、r 104_phi = b_rad(bp)else; memory, ball,wall,contact,clump,circle_phi = b_rad(bp) + b_rad(bp_other) end_if_sum = _sum + _rcp * c_kn(cp) * _phi end_ifend_ifif c_ball1(cp) = bp thencp = c_b1clist(cp) elsecp = c_b2clist(cp) end_ifend_loopbp = b_next(bp) end_loop_alpha = -1.0 * _lambda * tot_vol * _diso

42、 / _sumbp = ball_headloop whip # nullb_rad(bp) = (1.0 + _alpha)*b_rad(bp) bp = b_next(bp)end_loop enddef achieve_isostrloop while 1 # 0;始終執(zhí)行，除非遇到 exit 命令_diso = req_isostr - isostrif abs(_diso/req_isostr) = req_isostr_tol then exitend_ifoo = out(Current isotropic stress = + string(isostr) expand_rad

43、iidsolvemandend_loopenddef isostrii= measure(mp,1); 計算應(yīng)力 1 表示應(yīng)力，2 表示應(yīng)變isostr = (m_s11(mp) + m_s22(mp) / 2.0enddef get_mp；返回 measure id=1 mp = find_meas(1)endmeas id=1 x 5 y 2.5 rad 2.0 get_mphistory id=1 isostrSET req_isostr = -4e5req_isostr_tol = 0.01 cyc 20prprmeas 1isostrachieve_isostr6.接觸本構(gòu)模型在 P

44、FC2D 中，材料的本構(gòu)特性是通過接觸本構(gòu)模型來模擬的。每一顆粒的接觸本構(gòu)模型有:1)接觸剛度模型;2)滑動模型;3)連接模型。接觸剛度模型提供了接觸力和相對位移的彈性關(guān)系，滑動模型則強調(diào)切向和法向接觸力使得接觸顆?？梢园l(fā)生相對移動，而連接模型是限制總的切向和法向力使得在連接強度范圍內(nèi)發(fā)生接觸。1).接觸剛度模型接觸剛度通過以下兩式把接觸力與相對位移聯(lián)系起來，即:式中，是法向剛度，它為割線剛度，聯(lián)系總的法向力與位移。式中，是切向剛度，它為切線剛度，通過增量形式聯(lián)系切向力與位移。在 PFC2D 中有兩種接觸剛度模型，即線性接觸剛度模型與簡化的 Hertz-Minlin 接觸剛度模型,不同的接觸剛

45、度模型有不同的接觸剛度值。.線性接觸模型:分析顆粒間接觸特性時，將這種接觸的顆粒想象為一點在顆粒中心彈性梁，受力或力矩相當(dāng)于作用于顆粒的中心。這種以下特征參數(shù)描述:1)幾何參數(shù):長度(L)、斷面積(A)、慣性矩(I);2)變形參數(shù):楊氏模量 E、泊松比 v ; 3)強度參數(shù):法向強度、切向強度 ;顆粒間接觸楊氏模量用 表示，平行連接楊氏模量表示。同時假定 PFC2D 中所有顆粒均為厚度為 t 的圓盤。若顆粒 A 和 B 接觸，則梁的半徑為:梁的長度為:式中:和分別為接觸顆粒的半徑。相互接觸顆粒之間的受力特性，相當(dāng)于彈性致，此時梁斷面各 A 與慣性矩 I 為:受純軸向或純切向荷載時的情況一，式中

46、:t 為假設(shè)的顆粒圓盤厚度。關(guān)于接觸變形，不僅指接觸連接模型的接觸連接變形，同樣適用于顆粒之間以及顆粒與墻之間的接觸變形，所以，下面介紹的計算方法即可用于接觸連接材料也可用于非接觸連接材料(砂土)。對于純軸向荷載或純切向荷載，其法向接觸剛度與切向接觸剛度分別為：；式中， 是接觸楊氏模量，不同于材料整體楊氏模量(通常大于整體楊氏模量)。對于線性接觸模型，其接觸剛度 k是假設(shè)相互接觸兩球為串聯(lián)：式中，分別代表切向和法向剛度。若兩球有相同的剛度，即：此時，通過公式轉(zhuǎn)換，在接觸連接時，接觸模量與球的剛度間的關(guān)系式:所以在描述變形微觀參數(shù)時，先給定顆粒一顆粒接觸的接觸模量 EC 以及顆粒法向剛度與切向剛

47、度比值 kn/ks。根據(jù)式 1.計算kn。再據(jù)給定比值求 ks。線性接觸模型是通過兩個接觸體(顆粒顆粒，顆粒一墻體)的切向剛度 ks 和法向剛度kn 定義的。計算線性接觸模型的接觸剛度時假定接觸體為串聯(lián)，根據(jù)式2.計算。.Hertz-Mindlin 接觸模型Hertz-Mindlin 接觸模型是 Mindlin&Deresiewich (1953)和 Cundall(1988)理論近似得出的一種非線性接觸模型。它只能嚴(yán)格用于顆粒體接觸，而且不能再現(xiàn)剪切過程中的連續(xù)非線性(特別使用初始剪切模量時)。Hertz-Mindlin 模型由兩個接觸球體的參數(shù)剪切模量 G 和泊松比 v 來定義。該模型不適

48、用于采用接觸連接的兩球體的接觸，因為這種模型沒有定義球體受張力的情況。對于顆粒顆粒相互接觸，其彈性常數(shù)為平均值，對于顆粒墻體接觸，假定墻為剛性，所以直接取球的彈性常數(shù)。接觸法向割線剛度和切向切線剛度為:式中 Un 為顆粒體接觸量，F(xiàn)in 為法向接觸力其它的系數(shù)為兩接觸體的幾何與材料特性的函數(shù)。對于顆粒一顆粒相互接觸，這些系數(shù)可表示為:對于顆粒墻相互接觸，這些系數(shù)可表示為:式中，G 為彈性剪切模量，v 為泊松比，R 為顆粒的半徑，A, B為相互接觸的兩球體。2).滑動模型滑動模型是相互接觸球體的一種固有特性，它沒有法向抗拉強度，允許顆粒在抗剪強度范圍內(nèi)發(fā)生滑動，這種模型在接觸連接模型發(fā)生作用之前

49、一直有效。同時，滑動模型也可與半行連接模型同時起作用?；瑒幽Ｐ褪峭ㄟ^兩接觸體間最小摩擦系數(shù) u 定義的，若顆粒間則令法向和切向接觸力等于零。發(fā)生滑動的判別條件為:量小于或等于零，若,則可以發(fā)生滑動，其中：3).連接模型：PFC2D 模型允許相互接觸顆粒連接在一起，有兩種連接模型:即接觸連接與平行連接模型。接觸連接認(rèn)為連接只發(fā)生在接觸點很小范圍內(nèi)，而平行連接發(fā)生在接觸顆粒間圓形或方形有限范圍內(nèi)。接觸連接只能傳遞力，而平行連接同時能傳遞力和力矩。兩種類型的接觸可以同時存在，直到超過接觸強度。連接模型只能是顆粒之間的連接，顆粒與墻之間的連接不能采用連接模型。.接觸連接模型接觸連接可以想象為一對有恒定

50、法向剛度與切向剛度的彈簧作用于顆粒接觸點處，同時，這些彈簧設(shè)定一定的抗拉與抗剪強度。只要接觸連接存在就沒有顆粒間的滑動發(fā)生。接觸連接在顆粒間量小于零時，允許出現(xiàn)張力，但法向接觸張力過接觸連接強度。在 PFC2D 中，接觸連接是有法向連接強度 Fcn 定義。當(dāng)法向抗拉接觸力大于或等于法向接觸連接強度時，連接破壞并把法向、切向接觸力賦值為零。當(dāng)切向接觸力大于或等于切向連接強度時，連接也發(fā)生破壞，但是接觸力不發(fā)生變化，并假設(shè)切向力沒有超過摩擦極限。顆粒接觸點處接觸力與相對位移的關(guān)系的本構(gòu)特性見下圖。接觸力的法向分量接觸力的切向分量.平行連接模型平行連接模型可以描述顆粒之間有限范圍內(nèi)有夾層材料的木構(gòu)特

51、性。相互連接的兩個顆粒可以看作是球體或柱體。這種連接建立顆粒間一種彈性相互作用關(guān)系，可與前面所述的滑動模型或接觸連接模型同時存在。平行連接可以想象為一組有恒定法向剛度與切向剛度的彈簧均勻分布于接觸平面內(nèi)，這些彈簧作用的本構(gòu)關(guān)系類似與點接觸彈簧模擬顆粒剛度的本構(gòu)特性。接觸的相對運動在平行連接處產(chǎn)生力和力矩，作用于相互連接的顆粒上，且與連接材料的最向、切向應(yīng)力有關(guān)。平行連接模型是由法向剛度、切向剛度、法向強度、切向強度和連接半徑五個參數(shù)定義的。與平行連接相對應(yīng)的總接觸力和力矩用和表示，根據(jù)，總的接觸力和力矩表示平行連接在顆粒 B 上的作用，如下圖所示?？梢詫⒖偟慕佑|力沿接觸面分解為切向分量和法向分

52、量：式中，表示法向分量， 表示切向分量。7.賦予材料屬性 ： PFC2D 程序里，除了顆粒體聯(lián)結(jié)相關(guān)的屬性(聯(lián)結(jié)剛度、強度以及平行聯(lián)結(jié)半徑)外，其它顆粒體的屬性賦予均采用 PROPERTY 命令來執(zhí)行。同時，在執(zhí)行命令時，對與不同區(qū)域不同位置的顆粒體，其屬性的變化可通過gradient 和group 來配合使用，達(dá)到所需設(shè)置的要求和變化。如顆粒體間接觸剪切強度 隨著 X 方向線形變化，其表達(dá)式為：可以這樣來編寫命令流文件：prop s_bond 1e6 grad -1e4, 0 range x 0 50具體求解為：x=15m 處，s_bond=1e6+(-1e4)*15同時還可將某個區(qū)域內(nèi)的顆

53、粒體均定義為一個組，用 group 來執(zhí)行，然后再給這個組賦予屬性，具體操作為：group strong_rock range x 0 10 y 5 10group weak_rock range x 0 10 y 0 5 prop fric 1.0 range group strong_rock prop fric 0.1 range group weak_rock如果顆粒體的材料屬性是一個非線形變化的量，對于顆粒間摩擦系數(shù)表達(dá)式為：就需要采用 FISH 語言來賦予屬性,如其中 d 為顆粒體距模型頂面距離，h 為模型區(qū)域的寬度和高度,此時賦予屬性時, FISH語言可這樣來表述：def var

54、y_friction bp = ball_headloop whip # nully_= b_y(bp)depth = height - abs(y_)b_fric(bp) = fric0 * (1.0 + (depth/height)2) bp = b_next(bp)end_loopendSET height=5 fric0=0.2vary_friction8.節(jié)理面的生成及屬性設(shè)置 ： PFC2D理面可以用來模擬顆粒組間的滑移和分離，還可以模擬節(jié)理、斷層和層理等，同時還能來模擬兩種不同材料間的接觸面，如基礎(chǔ)和土體。采用 JSET 命令來設(shè)置，可生成單獨的一條或一組節(jié)理面。但只有存在法向力的顆粒間的接觸，才能被設(shè)置為節(jié)理面接觸，因此，節(jié)理面命令必須在顆粒體模型生成后并達(dá)到平衡后才能執(zhí)行。dip 和 origin可配合使用來生成所需位置的節(jié)理面，dip 是指節(jié)理面的傾角，origin是指節(jié)理面所經(jīng)過的某一點的坐標(biāo)。 如若生成一條 45 度角的單節(jié)理面， JSET id=1 dip=-45.0 origin=(0.0,2.0)可以采用下面令：number 和對于多組節(jié)理面，可采用spacing 來控制。JSET id=2 dip=0 origin=(3.0,2.0) & number=2 spacing=5.0JSET i