【文章內(nèi)容簡(jiǎn)介】 、 球形物密度以及 沖 擊速度等因素下進(jìn)行離散元計(jì)算，分析各參數(shù)對(duì)顆粒系統(tǒng)緩沖性能的影響。 （ 2）設(shè)置試驗(yàn)方案，以高速攝像拾取圖像并進(jìn)行 Matlab 數(shù)值計(jì)算為思路， 在不同顆粒粒徑、球形物大小 、 球形物密度以及 沖擊速度 等因素下開(kāi)展了試驗(yàn)研究。 （ 3）分析離散元模擬和試驗(yàn)結(jié)果，比較不同情況下盛顆粒容器底部受力以及球形物受力，探究球形物沖擊 下顆粒物質(zhì)的緩沖性能 影響因素，沖擊力 變化趨勢(shì)和沖擊力峰值 異同。 東北大學(xué)畢業(yè)設(shè)計(jì)（論文） 第 1 章 緒論 4 東北大學(xué)畢業(yè)設(shè)計(jì)（論文） 第 2章 沖擊載荷下的接觸模型理論分析 5 第 2 章 沖擊荷載下的接觸模型理論分析 DEM 離散元 程序 流程 當(dāng)顆粒在阻尼器中運(yùn)動(dòng)時(shí)，將所有顆粒都視為彈性小球以簡(jiǎn)化計(jì)算。 整個(gè)振動(dòng)響應(yīng)過(guò)程，被分隔成若干個(gè)極小的時(shí)間步 t? ，在一個(gè)時(shí)間步內(nèi)顆粒作勻速直線運(yùn)動(dòng)，如果獲得第 i 個(gè)時(shí)間步 it? 內(nèi)顆粒的位置、速度、受力情況，便可由牛頓運(yùn)動(dòng)定理求得第 i+1 個(gè)時(shí)間步 1it?? 時(shí)顆粒的位置、速度、受力情況。通過(guò)不斷迭代，最終求得整個(gè)阻尼系統(tǒng)在一段時(shí)間內(nèi)的響應(yīng)。 程序的前處理部分主要工作是文件的讀入和初始參數(shù)的計(jì)算。 本程序的讀入文件 為。在 中，定義了顆粒尺寸、顆粒個(gè)數(shù)、材料楊氏模量、泊松比、摩擦系數(shù)、計(jì)算時(shí)間步數(shù)等必要的程序運(yùn)行參數(shù)。 讀入文件后，程序的前處理部分會(huì)自動(dòng)計(jì)算出程序循環(huán)計(jì)算部分的必要參數(shù)，如時(shí)間步時(shí)長(zhǎng) t? 、顆粒之間的法向和切向剛度 ,NSKK、整個(gè)阻尼系統(tǒng)的背景空間大小等。 除此之外，給計(jì)算需要的數(shù)據(jù)數(shù)組開(kāi)辟內(nèi)存空間、分配顆粒的初始位置、給背景空間劃分網(wǎng)格等步驟也都是在程序的前處理部分完成的。 程序的前處理部分為之后的循環(huán)計(jì)算部分做 好了鋪墊。 循環(huán)計(jì)算部分主要包括三個(gè)大塊，分別是：搜索模塊、受力模塊和運(yùn)動(dòng)模塊。 搜索模塊的主要作用是初步判斷兩顆粒之間還有顆粒與墻壁之間是否能夠接觸。分為 GRID、 SEARCH、 SEARCH_WALL 三個(gè)主要子程序。本程序的搜索模塊采用的是基于背景網(wǎng)格的搜索模式，將固定的背景空間劃分為一定數(shù)量的網(wǎng)格，網(wǎng)格大小為邊長(zhǎng)小于最小顆粒半徑 sr 的 倍，每個(gè)網(wǎng)格內(nèi)只能存在一個(gè)顆粒的球心。如果兩顆粒球心處于相近的網(wǎng)格內(nèi)，便判斷其可能接觸。 受力模塊主要包括 FORCE_SINGLE 和 FORCE_WALL 兩個(gè)子程序，主要用于計(jì)算顆粒與顆粒之間還有顆粒與墻壁之間的相互作用力，計(jì)算方式會(huì)在 節(jié)詳細(xì)說(shuō)明。 運(yùn)動(dòng)模塊主要包括 MOVE 子程序，用于計(jì)算顆粒及箱體在力的作用下所做的運(yùn)動(dòng)。 東北大學(xué)畢業(yè)設(shè)計(jì)（論文） 第 2章 沖擊載荷下的接觸模型理論分析 6 下圖為 程序流程圖示 ： 圖 程序流程圖 東北大學(xué)畢業(yè)設(shè)計(jì)（論文） 第 2章 沖擊載荷下的接觸模型理論分析 7 顆粒間接觸受力 FORCE_SINGLE 子程序 用于計(jì)算顆粒間接觸時(shí)，每個(gè) 顆粒 受到的合力、合 力矩 ?；舅枷霝椋合到y(tǒng)中存在 NLIST 個(gè)“接觸對(duì)”，計(jì)算每個(gè)“接觸對(duì)”中的 A、 B 兩球相互作用時(shí)所受到的力和力矩，然后通過(guò)對(duì)所有“接觸 對(duì)”進(jìn)行循環(huán)累加，得到每個(gè)小球所受到的合力、合力矩。 圖 顆粒間接觸受力示意圖 兩球心連線為法線方向 n ，過(guò)接觸點(diǎn)且垂直法線的為切線方向 ? ； 兩球相互作用力可分成三種：法向力，切向力和滾動(dòng)摩阻 [1]。 速度求解均采用通過(guò)全局坐標(biāo)系下三個(gè)坐標(biāo)方向的分量進(jìn)行，其小球球心速度、角速度以及它們的分量為（如圖所示） A Ax Ay Az An AB Bx By Bz Bn Bv v v v v vv v v v v v ??? ? ? ? ????? ? ? ? ???r r r r r rr r r r r r； A A x A y A zB B x B y B z? ? ? ?? ? ? ?? ? ????? ? ???r r r rr r r r 程序采用如下模型，用彈簧、粘壺、耦合器和滑動(dòng)器來(lái)模擬表示顆粒間三類(lèi)接觸作用 [10]，即法向、切向和滾動(dòng) 。 Av xyzAxvAyvAzvAnvAv?Ax?Ay?Az?A B ?n() 東北大學(xué)畢業(yè)設(shè)計(jì)（論文） 第 2章 沖擊載荷下的接觸模型理論分析 8 圖 三類(lèi)接觸示意圖 一、 法向力計(jì)算 1 o v e r la p D O T NnN n nF k c? ? ? ? ? 其中：① 第一部分 nk overlap? 為最簡(jiǎn)單的線性模型，在 overlap 很小的情況下成立，是Hertz 接觸理論的一階展開(kāi)。 overlap 為兩球重疊量，即兩球半徑之和（ DAB ）與兩球心距離（ DS ）的差值。 ② 兩球接觸的法向有效剛度系數(shù) πn RA RBkERA RB?? ?, RA 、 RB 為小球半徑， E 為材料楊氏模量。 ③法向阻尼系數(shù) 2n n nc z k A M A SS? ， AMASS 為兩小球平均質(zhì)量， ? ?? ?22ln lnn e p sz e p s??? ?，eps 為恢復(fù)系數(shù)，由 input 文件給出 。 ④ DOTN 為兩球球心速度 Av ， Bv 的法向速度差 An Bnvv? 。 二、切向力計(jì)算 1n n ns s s s sF F f F k u? ? ? ? ? ? ? 判 斷 ：若 11nns s NFF???? ，則 11nns s NFF???? 若 11nns s NFF???? ，則 1n n ns s s s s s sF F f F k u c u? ? ? ? ? ? ? ? amp。 其中： （ 1） 切向剛度 snk k RNS?? ， RNS 為法向、切向剛度比，由 input 文件給出 。 （ 2） su? 為兩球接觸點(diǎn)切向相對(duì)滑移量，通過(guò)其分 量 ,sx sy szu u u? ? ? 計(jì)算而得，以 sxu?的計(jì)算為例進(jìn)行說(shuō)明： () () 東北大學(xué)畢業(yè)設(shè)計(jì)（論文） 第 2章 沖擊載荷下的接觸模型理論分析 9 DUTsxut? ? ?? DUT 為 A、 B 球接觸點(diǎn)切向方向速度差在 x 軸上的投影 Ax B xvv??? 程序中的公式： DUT=DUDOTN*XM+VRIATXVRIBTX DU： Av ， Bv 速度差在 x 軸的投影 Ax Bxvv? DOTN： Av ， Bv 的法向速度差 An Bnvv? DOTN*XM： Av Bv 法向速度差在 x 軸上的投影 Anx Bnxvv? DUDOTN*XM： A、 B 球心平動(dòng)引起的切向速度差在 x 軸上的投影，由總量差減去法向分量差而得， ()A x B x A x B x A n x B n xv v v v v v??? ? ? ? ? VRIATX： A 球轉(zhuǎn)動(dòng)引起的切向速度在 x 軸的投影 VRIBTX： B 球轉(zhuǎn)動(dòng)引起的切向速度在 x 軸的投影 VRIATXVRIBTX：由轉(zhuǎn)動(dòng)引起的 A、 B 球切向速度差在 x 軸的投影 。 （兩球接觸點(diǎn)的相對(duì)速度根據(jù)兩球在接觸點(diǎn)的絕對(duì)速度的差值求得，每個(gè)球在接觸點(diǎn)的絕對(duì)速度由隨著球心平動(dòng)和繞著球心轉(zhuǎn)動(dòng)這兩個(gè)速度矢量疊加而得。球心（平動(dòng)）速度對(duì)法向和切向速度均有貢獻(xiàn)，而繞著球心轉(zhuǎn)動(dòng)引起的速度只對(duì)接觸點(diǎn)的切線方向有貢獻(xiàn)。將 ? 分解成三個(gè)分量，分別計(jì)算三個(gè)轉(zhuǎn)動(dòng)分量引起的速度在三個(gè)坐標(biāo)軸上的投影，然后進(jìn)行疊加。 x? 不產(chǎn)生 x 方向的速度分量，以此類(lèi)推） 。 （ 3） 靜滑動(dòng)摩擦系數(shù) s? 在本子程序內(nèi)部定義并給出數(shù)值。 （ 4） 切向阻尼系數(shù) 2s s nc z k A M A S S??， snz rzs z??， rzs 法向、切向阻尼比，由 input給出 。 （ 5） su 也通過(guò)其分量計(jì)算而得，例如： sxu DUT? 。 三、 滾動(dòng)摩阻計(jì)算： 滾動(dòng)摩擦模型： m in ( , )r r r r nM k c r F? ? ??? 判斷： 若 11nnr r NM r F???? ，則 11nnr r NM r F???? ，（此處 1nNF? 不考慮粘滯力，與上同） 。 其中： （ 1） 滾動(dòng)摩阻的剛度系數(shù) 22rsR A R BkkR A R B???? ????? () () 東北大學(xué)畢業(yè)設(shè)計(jì)（論文） 第 2章 沖擊載荷下的接觸模型理論分析 10 （ 2） 靜滾動(dòng)摩阻系數(shù) ? ? 22 2 2224rD S R B R Ar A R F R BDS???? ? ? ?， DS 為兩球球心距離，ARF 未知 。 顆粒與邊界間接觸受力 Force_WALL 子程序 計(jì)算顆粒與邊界接觸時(shí)，邊界受到的合力及顆粒受到的合力、合力距。由于圓柱的特殊性，可把圓柱桶分成三個(gè)墻壁（ NWall 3? ），即下壁（ IW1? ）、圓 柱壁（ IW2? ）以及上壁（ IW3? ） . 接觸判斷：一個(gè)球最多可能和三個(gè)墻壁同時(shí)接觸 ? ?ip = na_w ic , ic 為 接觸對(duì)粒子編號(hào) 。 ? ?iw = nb_w ic *, ic* 接觸對(duì)墻編號(hào) 。 粒子球心坐標(biāo) ? ?x(ip), y(ip), z(ip)， 接觸點(diǎn)的坐標(biāo) ? ?x _ c o n t a c t , y _ c o n t a c t , z _ c o n t a c t. 圖 顆粒與邊界接觸示意圖 Fortran 程序中 上下墻壁 ， 如果 IW1? : x_contact = x(ip) y_contact = y(ip) z_contact = dx = x_contact x(ip)=0 dy = y_contact y(ip)=0 東北大學(xué)畢業(yè)設(shè)計(jì)（論文） 第 2章 沖擊載荷下的接觸模型理論分析 11 dz = z_contact z(ip) ds = sqrt(dx*dx + dy*dy + dz*dz) 如果 IW=3 ： x_contact = x(ip) y_contact = y(ip) z_contact = zPoint(1) （ zPoint(1)為上墻壁的 z 坐標(biāo)） dx = x_contact x(ip)=0 dy = y_contact y(ip)=0 dz = z_contact z(ip) ds = sqrt(dx*dx + dy*dy + dz*dz) 如果 IW2? ： 圖 圓柱壁受力示意圖 R_Location = sqrt( x(ip)*x(ip) + y(ip)*y(ip) ) （粒子球心距原點(diǎn)距離） R_Contact = R_Cylinder （接觸點(diǎn)距原點(diǎn)距離） x_contact = x(ip)*R_Contact/R_Location y_contact = y(ip)*R_Contact/R_Location z_Contact = z(ip) dx = x_contact x(ip) dy = y_contact y(ip) 東北大學(xué)畢業(yè)設(shè)計(jì)（論文） 第 2章 沖擊載荷下的接觸模型理論分析 12 dz = z_Contact z(ip)=0 ds = sqrt(dx*dx + dy*dy + dz*dz) 當(dāng) ds R_Particle 且 ds0 時(shí)，粒子和壁之間有接觸，此時(shí)粒子與墻之間的重疊量Overlap = R_Particle – ds。 顆粒碰撞后運(yùn)動(dòng)狀態(tài) MOVE 子程序 用于計(jì)算顆粒碰撞之后的運(yùn)動(dòng)。 得到顆粒在碰撞之后的位置和運(yùn)動(dòng)參數(shù)（包括顆粒的速度 U(I)、 V(I)、 W(I)，坐標(biāo) X(I)、 Y(I)、 Z(I)，角速度 OM1(I)、 OM2(I)、OM3(I)，轉(zhuǎn)角 THX(I)、 THY(I)、 THZ(I)， 將這些參數(shù)傳遞回 中的 mon公共區(qū)中，再提供給下一個(gè)時(shí)間步，重新計(jì)算顆粒碰撞產(chǎn)生的力和力矩，從而形成顆粒的運(yùn)動(dòng)和受力之間的循環(huán)模式。 在程序中，顆粒運(yùn)動(dòng)參數(shù)的計(jì)算方法采用如下中心差分格式 ， 用 Nt 時(shí)刻的加速度 /角加速度和12Nt?時(shí)刻的速度 /角速度來(lái)更新12Nt?時(shí)刻的速度 /角速度；用12Nt?時(shí)刻的速度 /角速度和 Nt 時(shí)刻的位移 /轉(zhuǎn)角更新 1Nt? 時(shí)刻的位移 /轉(zhuǎn)角。 計(jì)算單個(gè)顆粒的質(zhì)量 ： 3432D S I Z EA M A S S D E N S? ??? ???? 其中： AMASS 為 顆粒質(zhì)量； DSIZE