用離散元做動(dòng)力特性的比較少，前面使用比較經(jīng)典的落球試驗(yàn)證實(shí)了離散元中模擬振動(dòng)傳遞及其衰減的可行性。本文采用樁作為波的輸入裝置，在樁底平面處布置測(cè)點(diǎn)來(lái)研究振動(dòng)的衰減。基本的模擬順序?yàn)椋撼蓸印⒆灾亍⒉鍢丁⒔o參數(shù)、振動(dòng)。

這里成樣采用的分層壓縮法，生成1*0.5的試樣。

new def par    width=1.0    height= width*0.5         y_vel=10.0        cengshu=8.0    pengzhangxishu=5.0    tbjipei=table.create("jipei")    table(tbjipei,0.004)=0.1    table(tbjipei,0.005)=0.12    table(tbjipei,0.0055)=0.14    table(tbjipei,0.006)=0.3    table(tbjipei,0.0065)=0.6    table(tbjipei,0.007)=0.8    table(tbjipei,0.0075)=0.95    table(tbjipei,0.008)=1     poro=0.12    ;dengchaxian=0.02end@pardomain extent [-height*10] [height*10]set random 10001wall generate box [-width*0.5] [width*0.5] [-height*0.5] [height*(pengzhangxishu-1)-height*0.5] expand 1.5cmat default model linear method deformability emod 100e6 kratio 1.5 def create_sample    tb_ceng=table.create("ceng_poro")    loop n(1,cengshu)        dibuzuobiao=-height*0.5+height*(n-1)/cengshu        dingbuzuobiao=-height*0.5+height*(n)/cengshu                pengzhangGaoDu=dingbuzuobiao+height*(pengzhangxishu-1)/cengshu           command            wall delete walls range id 3            wall create id 3 vertices [-width*0.5*1.5] [pengzhangGaoDu] [width*0.5*1.5] [pengzhangGaoDu]         endcommand        poroceng=poro;-(n-(cengshu+1)*0.5)*dengchaxian/(cengshu-1)        table(tb_ceng,n)=poroceng          keligeoup=table.size(tbjipei)        Vs=height*width*(1-poroceng)/cengshu         cenggName=string.build("ceng%1",n)        loop m(1,keligeoup)            keligid=keligeoup+1-m            vslocal=0            if keligid=1 then                vslocal=Vs*table.y(tbjipei,keligid)            else                            vslocal=Vs*(table.y(tbjipei,keligid)-table.y(tbjipei,keligid-1))            endif            kelinum=vslocal/(math.pi*(table.x(tbjipei,keligid)*0.5)^2)            command                 ball generate radius [table.x(tbjipei,keligid)*0.5] number [math.ceiling(kelinum)] box [-width*0.5] [width*0.5] ...                    [dibuzuobiao] [pengzhangGaoDu] group @cenggName tries 20000                ball attribute density 2.7e3 damp 0.7            endcommand        endloop                 command            wall attribute yvel [-y_vel] range id 3            solve time [(pengzhangGaoDu-dingbuzuobiao)/y_vel]                        wall attribute yvel 0 range id 3            cycle 2000        endcommand            endloopend@create_samplesolvesave sample

成樣的結(jié)果為：

第二步是加自重：

restore sampleset gravity [80*9.8/width]cycle 1 solvesave zizhongtemp1wall delete walls range id 3 cycle 1 solvesave zizhong

自重后的結(jié)果為：

    第三步為插樁，這里修改了一下之前的插樁程序，將刪除顆粒與生成wall放在一起，而且將刪除顆粒的區(qū)域設(shè)置的比樁的區(qū)域大一點(diǎn)。這樣的物理意義為是挖土插樁而不是將樁擠入土中。

restore zizhong[chicunbi=80/width][pile_pos_x=-width*0.5*0.5][pile_length=28/chicunbi][pile_bizhi=0.9][pile_rad=0.5/chicunbi]def get_max_pos    tu_height=1e-100    loop foreach bp ball.list        ball_shang_pos_y=ball.pos.y(bp)+ball.radius(bp)        if ball_shang_pos_y>tu_height then            tu_height=ball_shang_pos_y        endif    endloopend@get_max_poswall generate box [pile_pos_x-pile_rad] [pile_pos_x+pile_rad] [tu_height-pile_length*0.9] [tu_height+pile_length*0.1] onewallball delete range x [pile_pos_x-pile_rad*1.3] [pile_pos_x+pile_rad*1.3] y [tu_height-pile_length*0.9] [tu_height+pile_length*0.1] cycle 1solvesave addpile

插入樁的狀態(tài)為：

最后一步為樁的振動(dòng)激勵(lì)：

restore canshuball attribute displacement multiply 0set mech age 0[vel_max=-2e-3][wp_pile=wall.find(5)]def addwallvel    whilestepping    pinlv=1+199*mech.age    y_vel=vel_max*math.sin(2*math.pi*pinlv*mech.age+math.pi)    wall.vel.y(wp_pile)=y_velend@addwallvel[bp1=ball.near(pile_pos_x+5/chicunbi,tu_height-pile_length*0.9)][bp2=ball.near(pile_pos_x+13.5/chicunbi,tu_height-pile_length*0.9)][bp3=ball.near(pile_pos_x+25/chicunbi,tu_height-pile_length*0.9)][bp4=ball.near(pile_pos_x+39/chicunbi,tu_height-pile_length*0.9)][ball.group(bp1,2)="jiance"][ball.group(bp2,2)="jiance"][ball.group(bp3,2)="jiance"][ball.group(bp4,2)="jiance"]def jiancewhilestepping    time=mech.age    weiyi1_y=ball.disp.y(bp1)    weiyi2_y=ball.disp.y(bp2)    weiyi3_y=ball.disp.y(bp3)    weiyi4_y=ball.disp.y(bp4)     sudu1_y=ball.vel.y(bp1)    sudu2_y=ball.vel.y(bp2)    sudu3_y=ball.vel.y(bp3)    sudu4_y=ball.vel.y(bp4)endhistory id 1 @timehistory id 2 @y_velhistory id 3 @weiyi1_yhistory id 4 @weiyi2_yhistory id 5 @weiyi3_yhistory id 6 @weiyi4_yhistory id 7 @sudu1_yhistory id 8 @sudu2_yhistory id 9 @sudu3_yhistory id 10 @sudu4_y[baocunpinlv=0.1] ;;這個(gè)是多長(zhǎng)時(shí)間保存一個(gè)save文件[time_record=mech.age-100][count=1]def savefile        if mech.age-time_record > baocunpinlv then        filename=string.build("jieguo%1",count)        command            save @filename                    endcommand        time_record=mech.age        count +=1    endif    endset fish callback -1.0 @savefilesolve time 2save result

    

    這里給樁施加一個(gè)余弦速度，可以設(shè)置峰值，這里以?huà)哳l的方式輸入，即頻率在一秒內(nèi)從1HZ增加到200HZ。

下面為0.3s內(nèi)速度的變化趨勢(shì)，可以看出還是滿(mǎn)足我們的需求的。

測(cè)點(diǎn)布置為：

    先看看1s內(nèi)的時(shí)程曲線(xiàn)，可以看出距離越遠(yuǎn)，其振動(dòng)是越弱的。而且傳遞的滯后性也很明顯。

    但是到0.3s左右的時(shí)候，隨著樁頻率的增加，土體的振動(dòng)也增加，而且最后的速度幅值超過(guò)了輸入幅值

    后面算到2s左右，振動(dòng)逐漸變的類(lèi)似于交通荷載振動(dòng)的感覺(jué)，不知道其內(nèi)在機(jī)理是什么，如果有朋友懂的話(huà)可以在公眾號(hào)留言交流一下。

看一下速度場(chǎng)：

0.2s左右還是比較正常的波傳遞的圖。

0.5s的時(shí)候，激勵(lì)的頻率逐漸提高，很明顯會(huì)在模型中出現(xiàn)波谷了。

算到2s的時(shí)候基本上就比較亂了，我的專(zhuān)業(yè)知識(shí)也不足以幫助我分析了。

將時(shí)程圖做傅里葉變化得到頻域圖。

下圖為樁輸入的頻域曲線(xiàn)

下圖為測(cè)點(diǎn)4的頻域曲線(xiàn)

    可以發(fā)現(xiàn)、。。。。。。嗯。。。什么也看不出來(lái)，分析不出啥玩意，求動(dòng)力學(xué)專(zhuān)家指導(dǎo)交流。

阻尼在離散元系統(tǒng)中有兩種，一種是介質(zhì)阻尼，一種是接觸阻尼。介質(zhì)阻尼的意思是你的模型是放在真空中，還是放在空氣中，或者是放在水中，不同介質(zhì)對(duì)模型的影響就是介質(zhì)阻尼。介質(zhì)阻尼的參數(shù)靠damp賦予，ball attribute damp就可以設(shè)置模型的介質(zhì)阻尼。接觸阻尼和介質(zhì)無(wú)關(guān)，這個(gè)和材料有關(guān)系，材料受力會(huì)出現(xiàn)徐變或者蠕變，其實(shí)這個(gè)大部分是由于接觸阻尼相關(guān)的，當(dāng)然蠕變或者徐變?cè)趲r石中還有劣化的影響。接觸阻尼的設(shè)置通過(guò)dp_nratio和dp_sratio實(shí)現(xiàn)。

    上一篇文章在介質(zhì)阻尼設(shè)置為0，其實(shí)相當(dāng)于將模型放在真空中，這里重新將介質(zhì)阻尼設(shè)置為0.7，這個(gè)特性可能是放在水里面。這個(gè)阻尼系數(shù)可能需要一些實(shí)驗(yàn)區(qū)標(biāo)定，這里先隨意設(shè)定一下。所有東西都是和之前一樣，在參數(shù)文件中改阻尼參數(shù)：

restore addpile[dp_n=0.1]ball attribute damp 0.7cmat add 1  model linear method deformability emod 100e6 kratio 1.5 property fric 0.5 dp_nratio [dp_n] dp_sratio [dp_n] range contact type ball-ballcmat add 2  model linear method deformability emod 100e6 kratio 1.5 property fric 0.5 dp_nratio [dp_n] dp_sratio [dp_n] range contact type ball-facetcmat applycleancycle 1solve aratio 1e-7save canshu

這里看一下最終的結(jié)果：

    下面為速度曲線(xiàn)，可以看出相當(dāng)于上一個(gè)工況，這個(gè)工況的振動(dòng)衰減現(xiàn)象還是比較明顯的，而且從曲線(xiàn)的振幅可以看出土層的自振頻率。當(dāng)輸入的頻率比較小的時(shí)候，其隨著頻率增大，振動(dòng)越大。當(dāng)頻率比較大的時(shí)候，頻率越大，振動(dòng)越小。

    對(duì)時(shí)程曲線(xiàn)進(jìn)行傅里葉變化得到頻域曲線(xiàn)，可以比較明顯的看出高頻波的衰減以及卓越頻率。可以看出高頻的衰減是比較明顯的，但是低頻的衰減不是特別好。

    將橫坐標(biāo)的區(qū)域改一下，可以看到峰值大概在100左右，除以尺寸比80為1.25Hz左右，也比較符合實(shí)際。

    所以在動(dòng)力學(xué)分析中damp和dp的取值可能還是需要標(biāo)定的，常見(jiàn)的標(biāo)定方式可能是標(biāo)定阻尼比，或者做共振柱試驗(yàn)，后面有時(shí)間或許做一下共振柱試驗(yàn)的模擬。