桩激励振动在砂土中的传播衰减

    用离散元做动力特性的比较少，前面使用比较经典的落球试验证实了离散元中模拟振动传递及其衰减的可行性。本文采用桩作为波的输入装置，在桩底平面处布置测点来研究振动的衰减。基本的模拟顺序为：成样、自重、插桩、给参数、振动。

这里成样采用的分层压缩法，生成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

成样的结果为：

第二步是加自重：

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

自重后的结果为：

    第三步为插桩，这里修改了一下之前的插桩程序，将删除颗粒与生成wall放在一起，而且将删除颗粒的区域设置的比桩的区域大一点。这样的物理意义为是挖土插桩而不是将桩挤入土中。

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

插入桩的状态为：

最后一步为桩的振动激励：

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] ;;这个是多长时间保存一个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

    

    这里给桩施加一个余弦速度，可以设置峰值，这里以扫频的方式输入，即频率在一秒内从1HZ增加到200HZ。

下面为0.3s内速度的变化趋势，可以看出还是满足我们的需求的。

测点布置为：

    先看看1s内的时程曲线，可以看出距离越远，其振动是越弱的。而且传递的滞后性也很明显。

    但是到0.3s左右的时候，随着桩频率的增加，土体的振动也增加，而且最后的速度幅值超过了输入幅值

    后面算到2s左右，振动逐渐变的类似于交通荷载振动的感觉，不知道其内在机理是什么，如果有朋友懂的话可以在公众号留言交流一下。

看一下速度场：

0.2s左右还是比较正常的波传递的图。

0.5s的时候，激励的频率逐渐提高，很明显会在模型中出现波谷了。

算到2s的时候基本上就比较乱了，我的专业知识也不足以帮助我分析了。

将时程图做傅里叶变化得到频域图。

下图为桩输入的频域曲线

下图为测点4的频域曲线

    可以发现、。。。。。。嗯。。。什么也看不出来，分析不出啥玩意，求动力学专家指导交流。

阻尼在离散元系统中有两种，一种是介质阻尼，一种是接触阻尼。介质阻尼的意思是你的模型是放在真空中，还是放在空气中，或者是放在水中，不同介质对模型的影响就是介质阻尼。介质阻尼的参数靠damp赋予，ball attribute damp就可以设置模型的介质阻尼。接触阻尼和介质无关，这个和材料有关系，材料受力会出现徐变或者蠕变，其实这个大部分是由于接触阻尼相关的，当然蠕变或者徐变在岩石中还有劣化的影响。接触阻尼的设置通过dp_nratio和dp_sratio实现。

    上一篇文章在介质阻尼设置为0，其实相当于将模型放在真空中，这里重新将介质阻尼设置为0.7，这个特性可能是放在水里面。这个阻尼系数可能需要一些实验区标定，这里先随意设定一下。所有东西都是和之前一样，在参数文件中改阻尼参数：

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

这里看一下最终的结果：

    下面为速度曲线，可以看出相当于上一个工况，这个工况的振动衰减现象还是比较明显的，而且从曲线的振幅可以看出土层的自振频率。当输入的频率比较小的时候，其随着频率增大，振动越大。当频率比较大的时候，频率越大，振动越小。

    对时程曲线进行傅里叶变化得到频域曲线，可以比较明显的看出高频波的衰减以及卓越频率。可以看出高频的衰减是比较明显的，但是低频的衰减不是特别好。

    将横坐标的区域改一下，可以看到峰值大概在100左右，除以尺寸比80为1.25Hz左右，也比较符合实际。

    所以在动力学分析中damp和dp的取值可能还是需要标定的，常见的标定方式可能是标定阻尼比，或者做共振柱试验，后面有时间或许做一下共振柱试验的模拟。