一種土體大變形流動(dòng)的計(jì)算模擬方法
【專利摘要】本發(fā)明公開了一種土體大變形流動(dòng)的計(jì)算模擬方法，該方法基于分子動(dòng)力學(xué)（簡(jiǎn)稱MD）方法，以牛頓運(yùn)動(dòng)方程為控制方程，引入Hertzian型摩擦公式和斯多克公式分別描述固?固顆粒間相互作用和固?液顆粒間相互作用，運(yùn)用Verlet算法求解積分過(guò)程，建立了固?液耦合的土體大變形流動(dòng)計(jì)算模型。本發(fā)明可以呈現(xiàn)土體顆粒復(fù)雜流動(dòng)的基本動(dòng)力學(xué)特征及其運(yùn)動(dòng)規(guī)律，有利于進(jìn)一步完善和提高我國(guó)防治土體大變形流動(dòng)災(zāi)害的研究水平，同時(shí)也為分子動(dòng)力學(xué)計(jì)算方法在土體大變形流動(dòng)災(zāi)害中的廣泛應(yīng)用和推廣奠定一定基礎(chǔ)。
【專利說(shuō)明】
一種土體大變形流動(dòng)的計(jì)算模擬方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明屬于巖土力學(xué)計(jì)算領(lǐng)域，涉及一種土體大變形流動(dòng)的計(jì)算模擬方法。
【背景技術(shù)】
[0002] 土體大變形流動(dòng)災(zāi)害主要表現(xiàn)為邊坡失穩(wěn)、基坑變形破壞和地面塌陷等諸多形 式，嚴(yán)重威脅人類生命與財(cái)產(chǎn)安全。土體大變形的分析對(duì)于工程設(shè)計(jì)和防災(zāi)減災(zāi)具有重要 的意義，但目前還沒(méi)有得到較好的解決，其中一個(gè)關(guān)鍵原因是數(shù)值計(jì)算并沒(méi)有完全反映土 體大變形的力學(xué)特征。土體在較大應(yīng)力狀態(tài)下的變形過(guò)程屬于非線性的大變形狀態(tài)，其變 形不再是符合小變形理論，鑒于傳統(tǒng)數(shù)值計(jì)算方法無(wú)法模擬土體大變形的完整過(guò)程。因此， 尋求一種有效的計(jì)算方法來(lái)模擬土體大變形流動(dòng)過(guò)程顯得尤為重要。

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

[0003] 本發(fā)明的目的在于解決現(xiàn)有數(shù)值計(jì)算方法無(wú)法有效模擬土體大變形流動(dòng)災(zāi)害的 技術(shù)問(wèn)題，進(jìn)而提出一種土體大變形流動(dòng)的計(jì)算模擬方法，本發(fā)明將Hertz ian型摩擦公式 和斯多克公式引入分子動(dòng)力學(xué)計(jì)算框架，運(yùn)用Verlet算法求解積分過(guò)程，建立了固-液耦合 的土體大變形流動(dòng)計(jì)算模型，可實(shí)現(xiàn)土體大變形流動(dòng)過(guò)程的模擬。
[0004] 為達(dá)到上述目的，本發(fā)明采用如下技術(shù)方案：
[0005] -種土體大變形流動(dòng)的計(jì)算模擬方法，其特征在于：所述土體大變形流動(dòng)的計(jì)算 模擬方法包括以下步驟：
[0006] 1)在分子動(dòng)力學(xué)數(shù)值計(jì)算模擬軟件LAMPS中輸入粒子信息文件：
[0007] 基于計(jì)算體系特征，首先設(shè)置計(jì)算體系的維數(shù)、單位、粒子類型、粒子半徑、邊界條 件以及積分算法，創(chuàng)建模擬單元;然后設(shè)置計(jì)算體系基本信息，所述計(jì)算體系基本信息包括 但不限于系綜、輸出的時(shí)間間隔、步長(zhǎng)以及運(yùn)行步數(shù);所述計(jì)算體系是要進(jìn)行數(shù)值計(jì)算模擬 的對(duì)象，所述要進(jìn)行數(shù)值計(jì)算模擬的對(duì)象是土體;所述計(jì)算體系特征是要進(jìn)行數(shù)值計(jì)算模 擬的對(duì)象的物理力學(xué)性質(zhì)及所處邊界條件；
[0008] 2)初始設(shè)置模塊啟動(dòng)，分子動(dòng)力學(xué)數(shù)值計(jì)算模擬軟件LAMMPS對(duì)計(jì)算體系的所有粒 子賦予初始位置以及初始速度，進(jìn)行初始化;輸入模型參數(shù)、粒子初始密度及粒子所受外 力;所述粒子是組成要進(jìn)行數(shù)值計(jì)算模擬的對(duì)象的分子、原子以及離子；
[0009] 3)模擬控制模塊啟動(dòng)，質(zhì)點(diǎn)鄰近搜索：
[0010]根據(jù)粒子影響半徑，運(yùn)用Verlet鄰近搜索法對(duì)計(jì)算粒子影響范圍內(nèi)所有粒子進(jìn)行 搜索，逐一判斷周圍粒子與計(jì)算粒子之間的間距與影響半徑之間的關(guān)系，確定周圍粒子與 計(jì)算粒子之間的間距小于影響半徑的周圍粒子作為計(jì)算粒子的鄰近粒子;直至確定每個(gè)粒 子的所有鄰近粒子；
[0011]所述計(jì)算粒子是正在被進(jìn)行計(jì)算的粒子；
[0012]所述周圍粒子是計(jì)算粒子周圍的粒子；
[0013]所述影響半徑是依據(jù)要進(jìn)行數(shù)值計(jì)算模擬的對(duì)象的顆粒粒徑大小設(shè)置的距離；
[0014] 4)計(jì)算每個(gè)粒子的所有鄰近粒子對(duì)這些粒子所產(chǎn)生的作用力、瞬時(shí)加速度以及瞬 時(shí)速度：
[0015] 根據(jù)步驟2)輸入的模型參數(shù)、粒子初始密度及粒子所受外力計(jì)算步驟3)所得到的 每個(gè)粒子的所有鄰近粒子對(duì)這些粒子所產(chǎn)生的作用力，然后通過(guò)矢量求和求得每個(gè)粒子上 所受合力，由牛頓運(yùn)動(dòng)方程求得瞬時(shí)加速度，采用Verlet積分算法對(duì)牛頓運(yùn)動(dòng)方程求解，得 到每個(gè)粒子的瞬時(shí)速度和位置；
[0016] 所述每個(gè)粒子的所有鄰近粒子對(duì)這些粒子所產(chǎn)生的作用力的計(jì)算方式是采用 Hertzian型摩擦公式和斯多克定律計(jì)算得到的，
[0017] 所述Hertzian型摩擦公式是：
[0018]
[0019]其中：
[0020] δ是第i個(gè)粒子和第j個(gè)粒子之間重疊的距離；
[0021] RhRj分別是第i個(gè)粒子以及第j個(gè)粒子的半徑；
[0022] kn、kt分別是第i個(gè)粒子和第j個(gè)粒子之間的法向彈性系數(shù)以及切向彈性系數(shù);所述 kt = 2/7kn；
[0023] 丫"、丫*分別是第i個(gè)粒子和第j個(gè)粒子之間的法向阻尼系數(shù)以及切向阻尼系數(shù);所 述 y t= 1/2 γ n;
[0024] nieff是第i個(gè)粒子和第j個(gè)粒子之間的有效質(zhì)量，所述meff=mimj/(mi+mj);
[0025] mi，mj分別為第i個(gè)粒子和第j個(gè)粒子的質(zhì)量；
[0026] Δ St是第i個(gè)粒子和第j個(gè)粒子之間的切向位移矢量；
[0027] 是連接第i個(gè)粒子和第j個(gè)粒子中心的單位矢量；
[0028] Vn、Vt分別是第i個(gè)粒子和第j個(gè)粒子相對(duì)速度的法向分量以及切向分量；
[0029] t是時(shí)間；
[0030] i是第i個(gè)粒子；
[0031] j是除第i個(gè)粒子外的所有粒子；
[0032]所述斯多克定律的計(jì)算公式是：
[0033] f = -3πη?ν
[0034] 其中：
[0035] η是水分子的粘度；
[0036] d是第i個(gè)粒子的粒徑；
[0037] V是第i個(gè)粒子的運(yùn)動(dòng)速度；
[0038] 所述每個(gè)粒子的瞬時(shí)速度的計(jì)算公式是：
[0039]
[0040] 其中：
[0041] i是第i個(gè)粒子，即正在被進(jìn)行計(jì)算的粒子；
[0042] t是時(shí)間；
[0043] Vi(t)是第i個(gè)粒子在t時(shí)刻的速度；
[0044] ri(t)是第i個(gè)粒子在t時(shí)刻的位置；
[0045] At是時(shí)間間隔；
[0046] ri(t_ Δ t))是第i個(gè)粒子在t_ Δ t時(shí)刻的位置；
[0047] ri(t+A t)是第i個(gè)粒子在t+Δ t時(shí)刻的位置；
[0048] 0(( At)2)是At的二階無(wú)窮?。?br>[0049] 5)重復(fù)步驟4)，對(duì)步驟4)所得到的每個(gè)粒子的所有鄰近粒子對(duì)這些粒子所產(chǎn)生的 作用力、瞬時(shí)加速度以及瞬時(shí)速度進(jìn)行更新；
[0050] 6)實(shí)時(shí)計(jì)算結(jié)果輸出模塊啟動(dòng)，統(tǒng)計(jì)所需物理量，求出模擬系統(tǒng)中粒子在相空間 的軌跡后，應(yīng)用統(tǒng)計(jì)物理原理得到系統(tǒng)的宏觀特性和結(jié)構(gòu)特點(diǎn)；
[0051] 7)判斷模擬時(shí)間是否大于設(shè)置的總時(shí)間，若模擬時(shí)間大于設(shè)置的總時(shí)間時(shí)，步驟 5)的更新過(guò)程結(jié)束，程序終止;所模擬時(shí)間小于設(shè)置的總時(shí)間時(shí)，重復(fù)步驟5)直至模擬時(shí)間 大于設(shè)置的總時(shí)間；所述模擬時(shí)間是步驟1)中步長(zhǎng)與運(yùn)行步數(shù)的乘積;所述設(shè)置的總時(shí)間 是步驟1)中步長(zhǎng)與運(yùn)行步數(shù)總和的乘積。
[0052]本發(fā)明的優(yōu)點(diǎn)如下：
[0053] (1)鑒于目前土體大變形流動(dòng)問(wèn)題數(shù)值模擬研究的現(xiàn)狀，以及傳統(tǒng)數(shù)值計(jì)算方法 的局限性，本發(fā)明提供了一種不帶近似、動(dòng)態(tài)跟蹤粒子軌跡、模擬結(jié)果精確，并且具有溝通 宏觀特性與微觀結(jié)構(gòu)作用的計(jì)算模擬方法一一分子動(dòng)力學(xué)方法，能夠有效再現(xiàn)土體大變形 流動(dòng)過(guò)程，解決了土體顆粒流動(dòng)帶來(lái)的相關(guān)軌跡動(dòng)態(tài)變化以及大變形等復(fù)雜問(wèn)題，可為防 治土體大變形流動(dòng)災(zāi)害提供相關(guān)指導(dǎo)依據(jù)。
[0054] (2)本發(fā)明首次將Hertzian型摩擦公式和斯多克公式引入MD框架中，建立了固-液 耦合的土體大變形流動(dòng)計(jì)算模型，該計(jì)算模擬方法有效彌補(bǔ)了當(dāng)前土體大變形流動(dòng)模擬精 確度不足的問(wèn)題。
[0055] (3)本發(fā)明首次應(yīng)用分子動(dòng)力學(xué)方法計(jì)算模擬土體大變形流動(dòng)問(wèn)題。分子動(dòng)力學(xué) 方法不要求模型過(guò)分簡(jiǎn)化，還能克服網(wǎng)格劃分帶來(lái)的問(wèn)題，提高了計(jì)算效率。
【附圖說(shuō)明】
[0056]圖1是本發(fā)明的流程圖；
[0057]圖2是砂土流動(dòng)最終構(gòu)型的計(jì)算模擬結(jié)果。
【具體實(shí)施方式】
[0058]下面結(jié)合附圖和具體的實(shí)施例對(duì)本發(fā)明的土體大變形流動(dòng)的計(jì)算模擬方法做進(jìn) 一步的詳細(xì)說(shuō)明：
[0059]參見(jiàn)圖1以及圖2,本發(fā)明提供了一種土體大變形流動(dòng)的計(jì)算模擬方法，該方法包 括以下步驟：
[0060] 1)在分子動(dòng)力學(xué)數(shù)值計(jì)算模擬軟件LAMPS中輸入粒子信息文件：
[0061] 基于計(jì)算體系特征，首先設(shè)置計(jì)算體系的維數(shù)、單位、粒子類型、粒子半徑、邊界條 件以及積分算法，創(chuàng)建模擬單元;然后設(shè)置計(jì)算體系基本信息，計(jì)算體系基本信息包括但不 限于系綜、輸出的時(shí)間間隔、步長(zhǎng)以及運(yùn)行步數(shù);計(jì)算體系是要進(jìn)行數(shù)值計(jì)算模擬的對(duì)象， 要進(jìn)行數(shù)值計(jì)算模擬的對(duì)象是土體;計(jì)算體系特征是要進(jìn)行數(shù)值計(jì)算模擬的對(duì)象的物理力 學(xué)性質(zhì)及所處邊界條件；
[0062] 2)初始設(shè)置模塊啟動(dòng)，分子動(dòng)力學(xué)數(shù)值計(jì)算模擬軟件LAMMPS對(duì)計(jì)算體系的所有粒 子賦予初始位置以及初始速度，進(jìn)行初始化;輸入模型參數(shù)、粒子初始密度及粒子所受外 力;粒子是組成要進(jìn)行數(shù)值計(jì)算模擬的對(duì)象的分子、原子以及離子；
[0063] 3)模擬控制模塊啟動(dòng)，質(zhì)點(diǎn)鄰近搜索：
[0064]根據(jù)粒子影響半徑，運(yùn)用Verlet鄰近搜索法對(duì)計(jì)算粒子影響范圍內(nèi)所有粒子進(jìn)行 搜索，逐一判斷周圍粒子與計(jì)算粒子之間的間距與影響半徑之間的關(guān)系，確定周圍粒子與 計(jì)算粒子之間的間距小于影響半徑的周圍粒子作為計(jì)算粒子的鄰近粒子;直至確定每個(gè)粒 子的所有鄰近粒子；
[0065] 計(jì)算粒子是正在被進(jìn)行計(jì)算的粒子；
[0066] 周圍粒子是計(jì)算粒子周圍的粒子；
[0067] 影響半徑是依據(jù)要進(jìn)行數(shù)值計(jì)算模擬的對(duì)象的顆粒粒徑大小設(shè)置的距離；
[0068] 4)計(jì)算每個(gè)粒子的所有鄰近粒子對(duì)這些粒子所產(chǎn)生的作用力、瞬時(shí)加速度以及瞬 時(shí)速度：
[0069] 根據(jù)步驟2)輸入的模型參數(shù)、粒子初始密度及粒子所受外力計(jì)算步驟3)所得到的 每個(gè)粒子的所有鄰近粒子對(duì)這些粒子所產(chǎn)生的作用力，然后通過(guò)矢量求和求得每個(gè)粒子上 所受合力，由牛頓運(yùn)動(dòng)方程求得瞬時(shí)加速度，采用Verlet積分算法對(duì)牛頓運(yùn)動(dòng)方程求解，得 到每個(gè)粒子的瞬時(shí)速度和位置；
[0070] 每個(gè)粒子的所有鄰近粒子對(duì)這些粒子所產(chǎn)生的作用力的計(jì)算方式是采用 Hertzian型摩擦公式和斯多克定律計(jì)算得到的，
[0071 ] Hertzian型摩擦公式是：
[0072]
[0073] 其中：
[0074] δ是第i個(gè)粒子和第j個(gè)粒子之間重疊的距離；
[0075] RnRj分別是第i個(gè)粒子以及第j個(gè)粒子的半徑；
[0076] kn、kt分別是第i個(gè)粒子和第j個(gè)粒子之間的法向彈性系數(shù)以及切向彈性系數(shù);kt = 2/7kn；
[0077] 丫"、丫*分別是第i個(gè)粒子和第j個(gè)粒子之間的法向阻尼系數(shù)以及切向阻尼系數(shù)； T t= 1/2 y η ；
[0078] meff是第i個(gè)粒子和第j個(gè)粒子之間的有效質(zhì)量，meff=mimj/(mi+mj);
[0079] nu，!^*別為第i個(gè)粒子和第j個(gè)粒子的質(zhì)量；
[0080] △ St是第i個(gè)粒子和第j個(gè)粒子之間的切向位移矢量；
[0081] 是連接第i個(gè)粒子和第j個(gè)粒子中心的單位矢量；
[0082] Vn、Vt分別是第i個(gè)粒子和第j個(gè)粒子相對(duì)速度的法向分量以及切向分量；
[0083] t是時(shí)間；
[0084] i是第i個(gè)粒子；
[0085] j是除第i個(gè)粒子外的所有粒子；
[0086]斯多克定律的計(jì)算公式：
[0087] f = -3 jrqd v
[0088] 其中：
[0089] η是水分子的粘度；
[0090] d是第i個(gè)粒子的粒徑；
[0091] V是第i個(gè)粒子的運(yùn)動(dòng)速度；
[0092] 每個(gè)粒子的瞬時(shí)速度的計(jì)算公式是：
[0093]
[0094] 兵干
[0095] i是第i個(gè)粒子，即正在被進(jìn)行計(jì)算的粒子；
[0096] t是時(shí)間；
[0097] Vi(t)是第i個(gè)粒子在t時(shí)刻的速度；
[0098] ri(t)是第i個(gè)粒子在t時(shí)刻的位置；
[0099] At是時(shí)間間隔；
[0100] ri(t_ Δ t))是第i個(gè)粒子在t_ Δ t時(shí)刻的位置；
[0101 ] r"t- Δ t))是第i個(gè)粒子在t+ Δ t時(shí)刻的位置；
[0102] 0(( At)2)是At的二階無(wú)窮??；
[0103] 5)重復(fù)步驟4)，對(duì)步驟4)所得到的每個(gè)粒子的所有鄰近粒子對(duì)這些粒子所產(chǎn)生的 作用力、瞬時(shí)加速度以及瞬時(shí)速度進(jìn)行更新；
[0104] 6)實(shí)時(shí)計(jì)算結(jié)果輸出模塊啟動(dòng)，統(tǒng)計(jì)所需物理量，求出模擬系統(tǒng)中粒子在相空間 的軌跡后，應(yīng)用統(tǒng)計(jì)物理原理得到系統(tǒng)的宏觀特性和結(jié)構(gòu)特點(diǎn)；
[0105] 7)判斷模擬時(shí)間是否大于設(shè)置的總時(shí)間，若模擬時(shí)間大于設(shè)置的總時(shí)間時(shí)，步驟 5)的更新過(guò)程結(jié)束，程序終止;所模擬時(shí)間小于設(shè)置的總時(shí)間時(shí)，重復(fù)步驟5)直至模擬時(shí)間 大于設(shè)置的總時(shí)間；模擬時(shí)間是步驟1)中步長(zhǎng)與運(yùn)行步數(shù)的乘積;所述設(shè)置的總時(shí)間是步 驟1)中步長(zhǎng)與運(yùn)行步數(shù)總和的乘積。
[0106] 實(shí)施例：
[0107]尺寸為40cmX IOcmX 12cm的模型箱水平放置，其中內(nèi)接箱的有效尺寸為IOcmX IOcmX 10cm，將普通中砂(粒徑0.25mm-0.5mm)緩慢均勻地裝入內(nèi)接箱直到8cm高度時(shí)停止 裝砂，快速抽掉內(nèi)接箱的擋板使砂土發(fā)生流動(dòng)。
[0108]應(yīng)用本發(fā)明的計(jì)算方法，對(duì)模型箱內(nèi)砂土流動(dòng)進(jìn)行計(jì)算模擬，具體計(jì)算步驟如下：
[0109] (1)應(yīng)用分子動(dòng)力學(xué)軟件LAMMPS計(jì)算時(shí)，設(shè)置計(jì)算區(qū)域長(zhǎng)為40cm，寬為10cm，高為 10cm，砂樣尺寸長(zhǎng)為10cm，寬為10cm，高為8cm，設(shè)置體系的單位為約化單位(無(wú)量綱），粒子 形狀視為圓球形，忽略砂粒的粒徑分布，邊界條件為固定邊界，左右兩側(cè)及下側(cè)為wall邊 界，上部與外界大氣接觸，計(jì)算直徑取粒徑期望值的3倍約為2_。
[0110] (2)砂樣初始位置與內(nèi)接箱裝好砂土?xí)r保持一致，初始速度設(shè)為零。
[0111] (3)運(yùn)用Verlet鄰近搜索法，以0.3倍間距，逐一判斷周圍粒子與計(jì)算粒子之間的 間距是否小于影響半徑，若小于則該粒子即為鄰近粒子，如此往復(fù)，逐一確定每個(gè)粒子所有 符合條件的鄰近粒子，建立鄰近列表。
[0112] (4)根據(jù)公式（1)和(2)計(jì)算體系中粒子間相互作用力，計(jì)算參數(shù)如表1所示，求出 粒子所受合力，進(jìn)而求得相應(yīng)的加速度，再根據(jù)公式(3)求得相應(yīng)的速度。
[0113] (5)重復(fù)步驟(4)，逐一更新每一步粒子的位置和速度。
[0114] (6)按照所需輸出相應(yīng)的物理量，如每一步的瞬時(shí)速度、加速度、位置、溫度等。
[0115] (7)當(dāng)模擬時(shí)間大于設(shè)置的總時(shí)間時(shí)，結(jié)束循環(huán)計(jì)算。最后輸出體系的熱力學(xué)數(shù) 據(jù):動(dòng)能、角動(dòng)量、體積、密度、時(shí)間步和計(jì)算速度等。
[0116] 表1計(jì)算參數(shù)表
[0118]最后說(shuō)明的是，以上實(shí)施例僅用以說(shuō)明本發(fā)明的技術(shù)方案而非限制，本領(lǐng)域的普 通技術(shù)人員應(yīng)當(dāng)理解，可以對(duì)本發(fā)明的技術(shù)方案進(jìn)行修改或者等同替換，而不脫離本發(fā)明 技術(shù)方案的宗旨和范圍，其均應(yīng)涵蓋在本發(fā)明的權(quán)利要求范圍當(dāng)中。
【主權(quán)項(xiàng)】
1. 一種土體大變形流動(dòng)的計(jì)算模擬方法，其特征在于:所述土體大變形流動(dòng)的計(jì)算模 擬方法包括以下步驟： 1) 在分子動(dòng)力學(xué)數(shù)值計(jì)算模擬軟件LAMMPS中輸入粒子信息文件： 基于計(jì)算體系特征，首先設(shè)置計(jì)算體系的維數(shù)、單位、粒子類型、粒子半徑、邊界條件以 及積分算法，創(chuàng)建模擬單元;然后設(shè)置計(jì)算體系基本信息，所述計(jì)算體系基本信息包括但不 限于系綜、輸出的時(shí)間間隔、步長(zhǎng)以及運(yùn)行步數(shù);所述計(jì)算體系是要進(jìn)行數(shù)值計(jì)算模擬的對(duì) 象，所述要進(jìn)行數(shù)值計(jì)算模擬的對(duì)象是土體;所述計(jì)算體系特征是要進(jìn)行數(shù)值計(jì)算模擬的 對(duì)象的物理力學(xué)性質(zhì)及所處邊界條件； 2) 初始設(shè)置模塊啟動(dòng)，分子動(dòng)力學(xué)數(shù)值計(jì)算模擬軟件LAMMPS對(duì)計(jì)算體系的所有粒子賦 予初始位置以及初始速度，進(jìn)行初始化;輸入模型參數(shù)、粒子初始密度及粒子所受外力;所 述粒子是組成要進(jìn)行數(shù)值計(jì)算模擬的對(duì)象的分子、原子以及離子； 3) 模擬控制模塊啟動(dòng)，質(zhì)點(diǎn)鄰近搜索： 根據(jù)粒子影響半徑，運(yùn)用Verlet鄰近搜索法對(duì)計(jì)算粒子影響范圍內(nèi)所有粒子進(jìn)行搜 索，逐一判斷周圍粒子與計(jì)算粒子之間的間距與影響半徑之間的關(guān)系，確定周圍粒子與計(jì) 算粒子之間的間距小于影響半徑的周圍粒子作為計(jì)算粒子的鄰近粒子;直至確定每個(gè)粒子 的所有鄰近粒子； 所述計(jì)算粒子是正在被進(jìn)行計(jì)算的粒子； 所述周圍粒子是計(jì)算粒子周圍的粒子； 所述影響半徑是依據(jù)要進(jìn)行數(shù)值計(jì)算模擬的對(duì)象的顆粒粒徑大小設(shè)置的距離;所述影 響半徑SL/2,其中L是模擬單元立方體盒子的邊長(zhǎng)； 4) 計(jì)算每個(gè)粒子的所有鄰近粒子對(duì)這些粒子所產(chǎn)生的作用力、瞬時(shí)加速度以及瞬時(shí)速 度： 根據(jù)步驟2)輸入的模型參數(shù)、粒子初始密度及粒子所受外力計(jì)算步驟3)所得到的每個(gè) 粒子的所有鄰近粒子對(duì)這些粒子所產(chǎn)生的作用力，然后通過(guò)矢量求和求得每個(gè)粒子上所受 合力，由牛頓運(yùn)動(dòng)方程求得瞬時(shí)加速度，采用Ver let積分算法對(duì)牛頓運(yùn)動(dòng)方程求解，得到每 個(gè)粒子的瞬時(shí)速度和位置； 所述每個(gè)粒子的所有鄰近粒子對(duì)這些粒子所產(chǎn)生的作用力的計(jì)算方式是采用 Hertzian型摩擦公式和斯多克定律計(jì)算得到的， 所述Hertzian型摩擦公式是：其中： S是第i個(gè)粒子和第j個(gè)粒子之間重疊的距離； 心、心分別是第i個(gè)粒子以及第j個(gè)粒子的半徑； kn、kt分別是第i個(gè)粒子和第j個(gè)粒子之間的法向彈性系數(shù)以及切向彈性系數(shù);所述kt = 2/7kn； Yn、Yt分別是第i個(gè)粒子和第j個(gè)粒子之間的法向阻尼系數(shù)以及切向阻尼系數(shù);所述 T t= 1/2 y n ； nfeff是第i個(gè)粒子和第j個(gè)粒子之間的有效質(zhì)量，所述mrff=mimj/(mi+mj); mi，mj分別為第i個(gè)粒子和第j個(gè)粒子的質(zhì)量； A st是第i個(gè)粒子和第j個(gè)粒子之間的切向位移矢量； 是連接第i個(gè)粒子和第j個(gè)粒子中心的單位矢量； vn、vt分別是第i個(gè)粒子和第j個(gè)粒子相對(duì)速度的法向分量以及切向分量； t是時(shí)間； i是第i個(gè)粒子； j是除第i個(gè)粒子外的所有粒子； 所述斯多克定律的計(jì)算公式是： f = _3Jiqdv 其中： n是水分子的粘度； d是第i個(gè)粒子的粒徑； v是第i個(gè)粒子的運(yùn)動(dòng)速度； 所述每個(gè)粒子的瞬時(shí)速度的計(jì)算公式是：其中： i是第i個(gè)粒子，即正在被進(jìn)行計(jì)算的粒子； t是時(shí)間； Vi(t)是第i個(gè)粒子在t時(shí)刻的速度； n (t)是第i個(gè)粒子在t時(shí)刻的位置； At是時(shí)間間隔； ri(t_ A t))是第i個(gè)粒子在t-A t時(shí)刻的位置； ri(t+A t)是第i個(gè)粒子在t+A t時(shí)刻的位置； 〇(( At)2)是At的二階無(wú)窮小； 5) 重復(fù)步驟4)，對(duì)步驟4)所得到的每個(gè)粒子的所有鄰近粒子對(duì)這些粒子所產(chǎn)生的作用 力、瞬時(shí)加速度以及瞬時(shí)速度進(jìn)行更新； 6) 實(shí)時(shí)計(jì)算結(jié)果輸出模塊啟動(dòng)，統(tǒng)計(jì)所需物理量，求出模擬系統(tǒng)中所有粒子在相空間 的軌跡后，應(yīng)用統(tǒng)計(jì)物理原理得到計(jì)算體系的宏觀特性和結(jié)構(gòu)特點(diǎn)； 7) 判斷模擬時(shí)間是否大于設(shè)置的總時(shí)間，若模擬時(shí)間大于設(shè)置的總時(shí)間時(shí)，步驟5)的 更新過(guò)程結(jié)束，程序終止;所模擬時(shí)間小于設(shè)置的總時(shí)間時(shí)，重復(fù)步驟5)直至模擬時(shí)間大于 設(shè)置的總時(shí)間；所述模擬時(shí)間是步驟1)中步長(zhǎng)與運(yùn)行步數(shù)的乘積;所述設(shè)置的總時(shí)間是步 驟1)中步長(zhǎng)與運(yùn)行步數(shù)總和的乘積。
【文檔編號(hào)】G06F17/50GK106055851SQ201610573468
【公開日】2016年10月26日
【申請(qǐng)日】2016年7月19日
【發(fā)明人】趙紫陽(yáng), 申俊敏, 張軍, 孫志杰, 薛曉輝, 馬林
【申請(qǐng)人】山西省交通科學(xué)研究院, 山西交科公路勘察設(shè)計(jì)院