離散元方法與實(shí)踐慕課matdem-第九章2段9.2礦山斜坡

離散元方法與實(shí)踐劉春南京大學(xué)9.1礦山斜坡爆破0概述本節(jié)將介紹礦山爆破示例。該模擬基于BoxModel示例代碼，這個代碼最初用于二維和三維滑坡數(shù)值模擬。本示例需先運(yùn)行BoxModel1和BoxModel2代碼來生成斜坡，再利用Model3Exploisive代碼實(shí)現(xiàn)礦山爆破過程的數(shù)值模擬。爆破模擬的基本原理為：搜索出一定區(qū)域內(nèi)的單元，定義為爆炸點(diǎn)，增加這些單元的半徑，使單元間相互重疊，以產(chǎn)生巨大壓力和彈性應(yīng)變能。在數(shù)值模擬過程中，爆破點(diǎn)的單元在內(nèi)壓力作用下，迅速向外運(yùn)動，釋放應(yīng)變能，并模擬爆破作用。1建立斜坡模型BoxModel1與其他示例的第一步相同，在此不作介紹。本節(jié)直接使用第一步得到的二維堆積模型。在加載第一步生成的堆積模型數(shù)據(jù)后，BoxModel2的主要建模代碼分為導(dǎo)入數(shù)據(jù)并切割模型、以及設(shè)置材料和平衡模型兩個部分1.1導(dǎo)入數(shù)據(jù)并切割模型首先創(chuàng)建Tool_Cut對象C（第1行）；并采用Tool_Cut對象C和addSurf命令，讀入記事本文檔中的高程數(shù)據(jù)（第2、3行）；C=Tool_Cut(d);%cutthemodellSurf=load('slope/layersurface.txt');%loadthesurfacedataC.addSurf(lSurf);%addthesurfacestothecutC.setLayer({'sample'},[1,2,3,4]);%setlayersaccordinggeometricaldatagNames={'lefPlaten';'rigPlaten';'botPlaten';'layer1';'layer2';'layer3'};d.makeModelByGroups(gNames);為方便查看層面，這些數(shù)據(jù)記錄于“滑坡高程.xls”的BoxModel數(shù)據(jù)表中，具體數(shù)據(jù)和格式如表中所示，相應(yīng)的excel折線表如圖所示。在Excel中設(shè)定完高程數(shù)據(jù)后，再將其復(fù)制到txt文檔中，以方便MatDEM讀取。折線號surface0surface1surface2surface3坐標(biāo)xyzxyzxyzxyz高程數(shù)據(jù)1001029103310434004026403040427007022702670421000100161002010038………………………………高程數(shù)據(jù).xls文件中的折線表Excel表中用于生成折線圖的矩陣數(shù)據(jù)在圖中，我們只定義了水平方向0-80米處的高程，MatDEM會自動根據(jù)已有的高程數(shù)據(jù)向外延伸，切分出80-100米的地層（圖b中）。a高程數(shù)據(jù).xls文件中的折線表b由折線表數(shù)據(jù)生成的斜坡模型然后，通過直接設(shè)定材料性質(zhì)的方法來生成兩種材料Soil1和Soil2，得到材料數(shù)組Mats（第1-8行）。%----------setmaterialofmodelmatTxt=load('Mats\Soil1.txt');Mats{1,1}=material('Soil1',matTxt,B.ballR);Mats{1,1}.Id=1;matTxt2=load('Mats\Soil2.txt');Mats{2,1}=material('Soil2',matTxt2,B.ballR);Mats{2,1}.Id=2;d.Mats=Mats;%---------assignmaterialtolayersandbalancethemodeld.setGroupMat('layer2','Soil2');d.groupMat2Model({'sample','layer2'});d.balanceBondedModel();1.2設(shè)置材料和平衡模型如表所示，Soil2的強(qiáng)度只有Soil1的十分之一，且內(nèi)摩擦系數(shù)也較小，將作為軟弱層。由于單元默認(rèn)的材料號為1，只需將軟弱層（layer2）的材料設(shè)為Soil2，并使用groupMat2Model將材料性質(zhì)賦給單元。最后平衡模型，得到如圖所示的離散元模型，其中l(wèi)ayer2為軟弱層。材料名稱楊氏模量/MPa泊松比單軸拉伸強(qiáng)度/kPa單軸壓縮強(qiáng)度/kPa內(nèi)摩擦系數(shù)密度/kg/m3Soil1200.14202000.81900Soil2100.182200.62000Soil1和Soil2兩種材料的輸入力學(xué)性質(zhì)

2設(shè)置爆破點(diǎn)和爆破能量模型建好后，進(jìn)入第三步代碼Model3Exploision。由于第二步默認(rèn)的材料力學(xué)性質(zhì)較弱，為保證邊坡穩(wěn)定，此處先將單元的斷裂力（aBF）和初始抗剪力（aFS0）增加為10倍（第2、3行）。…B.name=[B.name'Exploision'];d.mo.aBF=d.mo.aBF*10;d.mo.aFS0=d.mo.aFS0*10;然后通過代碼定義爆破點(diǎn)在斜坡中的位置和半徑，搜索出相應(yīng)的單元編號bombId，并將其定義為Bomb1組。centerX=15;centerZ=20;bombR=2;dX=d.mo.aX-centerX;dZ=d.mo.aZ-centerZ;bombId=find((dX.*dX+dZ.*dZ)<bombR.*bombR);%gettheIdofbombd.addGroup('Bomb1',bombId);%addanewgroupd.mo.zeroBalance();d.recordStatus();我們將爆破點(diǎn)的單元組號定義為1，利用d.show函數(shù)顯示出爆破點(diǎn)的位置（圖）；然后計算平均應(yīng)力Stress，并保存于d.data.Stress中，通過d.show(‘Stress’)顯示結(jié)果，實(shí)現(xiàn)了自定義后處理顯示。d.setData();d.data.groupId(d.GROUP.Bomb1)=1;d.show('groupId');d.data.Stress=(abs(d.data.StressXX)+abs(d.data.StressZZ));d.show('Stress');爆破點(diǎn)單元示意圖

（a）組號圖；

（b）自定義平均應(yīng)力圖

在獲得爆破單元的基礎(chǔ)上，通過增加其半徑來產(chǎn)生爆破能量。為此，先將單元半徑增加40%，并記錄新的系統(tǒng)狀態(tài)。oldKe=d.status.elasticEs(end);%recordtheoriginalenergybombExpandRate=1.4;d.mo.aR(bombId)=d.mo.aR(bombId)*bombExpandRate;%increasebombelementsized.mo.zeroBalance();d.recordStatus();newKe=d.status.elasticEs(end);dKe=newKe-oldKe;%calculatetheenergyincrementfs.disp(['Energyofthebombis'num2str(dKe)'J','~='num2str(dKe/4.2e6)'KgTNT']);可以看到，爆破點(diǎn)處單元的應(yīng)力非常大。為了計算爆破能量，利用oldKe記錄原始的彈性應(yīng)變能（第1行），用newKe記錄單元半徑增加后的彈性應(yīng)變能（第6行），計算得到由于爆破點(diǎn)單元半徑增加而增加的彈性應(yīng)變能。最后一條命令，將能量轉(zhuǎn)化為TNT當(dāng)量，以量化爆破能量。運(yùn)行代碼后，會得到以下信息Energyofthebombis17973262.3911J~=4.2793KgTNT即模擬的爆破能量相當(dāng)于4.3公斤的TNT炸藥。需要注意的是，在本示例中，為了較清晰地顯示爆破點(diǎn)的位置，采用了較大的爆破點(diǎn)半徑。在實(shí)際應(yīng)用中，可根據(jù)需要設(shè)置較小的區(qū)域半徑，并通過增加爆破點(diǎn)單元的剛度來增加爆破能量。同時，也可以通過給單元施加較大初速度的方法來增加能量，等。3迭代計算和模擬結(jié)果進(jìn)一步，可以開始迭代計算。由于本例為動力作用，需要將活動單元的阻尼系數(shù)（d.mo.mVis）設(shè)小一些，相當(dāng)于模擬空氣阻力（1、2行）。同時，運(yùn)行d.mo.setGPU(‘a(chǎn)uto’)命令，測試GPU和CPU的速度，選擇更快的一個，用于后續(xù)計算（第3行）

。visRate=0.00001;d.mo.mVis=d.mo.mVis*visRate;gpuStatus=d.mo.setGPU('auto');在設(shè)置好模型后，通過循環(huán)完成數(shù)值模擬：totalCircle=20;d.tic(totalCircle);fName=['data/step/'B.namenum2str(B.ballR)'-'num2str(B.distriRate)'loopNum'];save([fName'0.mat']);fori=1:totalCircled.mo.setGPU(gpuStatus);d.balance('Standard',0.1);d.clearData(1);save([fNamenum2str(i)'.mat']);d.calculateData();d.toc();%showthenoteoftimeendd.show('mV');這個模擬分為20次循環(huán)（totalCircle），每次循環(huán)進(jìn)行0.1次標(biāo)準(zhǔn)平衡（第7行）共模擬真實(shí)世界時間5s。圖a為0.1s時的單元速度場圖，在圖中可以清晰地看到，應(yīng)力波以爆破點(diǎn)為中心向外傳播，并在坡面形成較大的單元速度場（40m/s）。單元速度場分布