網(wǎng)上經(jīng)典-38pfc匯總

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

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

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

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

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

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

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

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

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

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

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

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

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

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ù)定義對(duì)象,其中一類(lèi)是尺度變量(scalar variables),它們是單個(gè)的數(shù)字,下面是總的尺度變量:clock時(shí)鐘時(shí)間,unbal最

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

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

17、規(guī)則排列的顆粒流，主要采用 FISH 語(yǔ)言配合 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ī)則排列；無(wú)規(guī)則排列，即：對(duì)一個(gè)給定空隙率的區(qū)域，采用顆粒來(lái)充填其中需要進(jìn)行填充的空隙，并確保整個(gè)模型保持平衡。對(duì)于所能被填充的模型的初始空隙率，是有一個(gè)限制值，不能任意小。對(duì)于某些空隙率的模型，顆粒的填充可以無(wú)接觸地排列，對(duì)于其它情況的空隙率，顆粒又

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

20、3.2.1 半徑擴(kuò)展法首先需要介紹一個(gè)概念：半徑放大系數(shù) m，它會(huì)引起模型空隙率的改變?？障堵实亩x：顆粒體的總面積； 封閉區(qū)域的總面積因此：定義初始空隙率為，新生成的空隙率為,那么，則有實(shí)例：區(qū)域?qū)挘?0區(qū)域高：5目標(biāo)空隙率：0.12顆粒體數(shù)目：300最大最小半徑比：1.5初始假設(shè)一個(gè) m 之值，可以求出初始的為從而得到平均半徑大小為：由初始給定的顆粒體最大直徑和最小直徑的比例 r到顆粒體半徑的上下限：通過(guò)上面的計(jì)算，同時(shí)編寫(xiě) FISH 函數(shù)來(lái)生成滿足需要的空隙率的模型，并能自動(dòng)計(jì)算空隙率的大小，以便用戶進(jìn)行檢驗(yà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)生來(lái)達(dá)到所需要的空隙率。但這種方法所 帶來(lái)的缺點(diǎn)是可能在局部區(qū)域造成大面積的顆粒，這將會(huì)產(chǎn)生很大的擠壓力，從而給予顆粒較大的初始速度，這就可能使得顆粒體脫離墻體的限制。為避免此情況的發(fā)生，可通過(guò)初始的有限步循環(huán)計(jì)算將顆粒的動(dòng)能減至零，然后再計(jì)算至平衡態(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 運(yùn)用模型生成“過(guò)濾器”有些情況需要建立復(fù)雜區(qū)域形狀的顆粒流模型，如右圖，此時(shí)可運(yùn)用模型生成過(guò)濾器來(lái)獲得所需要的模型，即 filter 命令，其后由用戶定義 FISH 函數(shù)來(lái)控制,其中,顆粒的半徑通過(guò) fc_arg(0

30、)進(jìn)行檢驗(yàn),x 和y的坐標(biāo)位置分別通過(guò) fc_arg(1) 和 fc_arg(2)進(jìn)行檢驗(yàn)。如果顆粒滿足要求，則 FISH 函數(shù)值設(shè)為0，否則為 1。詳見(jiàn)下例。newde_rect;用戶定義生成過(guò)濾器生成方形環(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).墻體邊界建模過(guò)程中，墻體可作為顆粒體的生成范圍約束，但同時(shí)也可以將墻體作為邊界來(lái)

34、施加約束。對(duì)于墻體，只能施度約束，而不能直接對(duì)其施加外力，因?yàn)檫\(yùn)動(dòng)定律對(duì)墻體是不適用的。其速度由以下三個(gè)參數(shù)控制:線速度、角速度和旋轉(zhuǎn)中心。墻的運(yùn)動(dòng)是通過(guò)不斷更新定義墻的基點(diǎn)的位置來(lái)描述。采用 WALL 命令設(shè)置，如：wall id=1 x=1.0 y=1.0 spin=10.0(2)顆粒體邊界PFC2D 中的模型可以將一連串的顆粒體作為邊界條件?；痉椒ǎ涸谀Ｐ途o密壓縮至平衡后，通過(guò) FISH 函數(shù)將與墻體相接觸的顆粒體逐個(gè)提取，將這一系列的顆粒體采用共同的邊界條件限制，最后刪除初始的限制墻體，即實(shí)現(xiàn)了以顆粒體代替墻體來(lái)作為邊界條件。顆粒體邊界條件分速度邊界和外力邊界，下例為速度邊界程序?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顆粒體外力邊界：其前期操作與速度邊界類(lèi)似，差異處是先將邊界顆粒體的速度初始化為零，再刪除限制墻體，采用 FISH 函數(shù)反向施加顆粒體所受的不平衡力，運(yùn)行一定的計(jì)算步至平衡態(tài)。在此過(guò)程中，區(qū)域內(nèi)部的顆粒體會(huì)因輕微的擾動(dòng)而導(dǎo)致其自身產(chǎn)生運(yùn)動(dòng)而脫離邊界顆粒體，此時(shí)前面將邊界顆粒體的速度固定為零。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)混合邊界條件：混合邊界即速度邊界與外力邊界同時(shí)存在，對(duì)于雙軸壓縮試驗(yàn)，其兩側(cè)端面不允許存在明顯的幾何變形，否則其邊界力將會(huì)無(wú)效。因此，此時(shí)需要用到混合邊界，兩側(cè)端面采用外力邊界來(lái)約束，確保其不發(fā)生幾何變形，側(cè)面 Y 向采用速度邊界，使其速度線性增加，以此來(lái)模擬試塊軸向變形。5.初始條件 : 首先介紹一下應(yīng)力的概念。在連續(xù)介質(zhì)中，應(yīng)力

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

40、 相連的下一個(gè)顆粒的物理地址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); 計(jì)算應(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)特性是通過(guò)接觸本構(gòu)模型來(lái)模擬的。每一顆粒的接觸本構(gòu)模型有:1)接觸剛度模型;2)滑動(dòng)模型;3)連接模型。接觸剛度模型提供了接觸力和相對(duì)位移的彈性關(guān)系，滑動(dòng)模型則強(qiáng)調(diào)切向和法向接觸力使得接觸顆?？梢园l(fā)生相對(duì)移動(dòng)，而連接模型是限制總的切向和法向力使得在連接強(qiáng)度范圍內(nèi)發(fā)生接觸。1).接觸剛度模型接觸剛度通過(guò)以下兩式把接觸力與相對(duì)位移聯(lián)系起來(lái)，即:式中，是法向剛度，它為割線剛度，聯(lián)系總的法向力與位移。式中，是切向剛度，它為切線剛度，通過(guò)增量形式聯(lián)系切向力與位移。在 PFC2D 中有兩種接觸剛度模型，即線性接觸剛度模型與簡(jiǎn)化的 Hertz-Minlin 接觸剛度模型,不同的接觸剛

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

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

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

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

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

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

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

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

53、粒體均定義為一個(gè)組，用 group 來(lái)執(zhí)行，然后再給這個(gè)組賦予屬性，具體操作為：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如果顆粒體的材料屬性是一個(gè)非線形變化的量，對(duì)于顆粒間摩擦系數(shù)表達(dá)式為：就需要采用 FISH 語(yǔ)言來(lái)賦予屬性,如其中 d 為顆粒體距模型頂面距離，h 為模型區(qū)域的寬度和高度,此時(shí)賦予屬性時(shí), FISH語(yǔ)言可這樣來(lái)表述：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理面可以用來(lái)模擬顆粒組間的滑移和分離，還可以模擬節(jié)理、斷層和層理等，同時(shí)還能來(lái)模擬兩種不同材料間的接觸面，如基礎(chǔ)和土體。采用 JSET 命令來(lái)設(shè)置，可生成單獨(dú)的一條或一組節(jié)理面。但只有存在法向力的顆粒間的接觸，才能被設(shè)置為節(jié)理面接觸，因此，節(jié)理面命令必須在顆粒體模型生成后并達(dá)到平衡后才能執(zhí)行。dip 和 origin可配合使用來(lái)生成所需位置的節(jié)理面，dip 是指節(jié)理面的傾角，origin是指節(jié)理面所經(jīng)過(guò)的某一點(diǎn)的坐標(biāo)。 如若生成一條 45 度角的單節(jié)理面， JSET id=1 dip=-45.0 origin=(0.0,2.0)可以采用下面令：number 和對(duì)于多組節(jié)理面，可采用spacing 來(lái)控制。JSET id=2 dip=0 origin=(3.0,2.0) & number=2 spacing=5.0JSET i