首页/文章/ 详情

PFC模拟三维单剪实验

lobby

2年前浏览4259

二维单剪是比较简单的，在双轴的基础上改是比较好实现的。但是三维单剪有比较多的细节问题需要解决，而且对于结果的分析都是比较困难的事情。本文主要针对于三维单剪的建模过程和应力分析进行讲解。

一、单剪实验

大家接触比较多的可能是直剪实验，上下两个剪切盒横向移动，在剪切面上产生剪切力使得式样发生破坏。而单剪实验相当于很多个剪切盒堆在一起进行剪切，相对于直剪实验，更加符合土体的变形特性。

二、单剪实验建模

1、成样

这一步和常规三轴或者巴西劈裂一样，我们需要一个圆柱形的式样，注意这里的是一个扁圆柱样。

new def chicun_par         banjing=0.3    sample_hight=banjing*4/7.0     keli_rdmin=0.006    keli_rdmax=0.009end@chicun_par





domain extent [-banjing*1.5] [banjing*1.5] [-banjing*1.5] [banjing*1.5] ...    [-sample_hight*0.5*1.5] [sample_hight*0.5*1.5]

[n=1.4]wall generate cylinder base 0 0 [-sample_hight*0.5*n] axis 0 0 1 ...                    height [sample_hight*n] radius [banjing] cap false false

wall generate plane position 0 0 [sample_hight*0.5] dip 0 ddir 0wall generate plane position 0 0 [-sample_hight*0.5] dip 0 ddir 0

ball distribute radius [keli_rdmin] [keli_rdmax] porosity 0.28 ...    range cylinder end1 0 0 [sample_hight*0.5-keli_rdmin] ...    end2 0 0 [-sample_hight*0.5 keli_rdmin] radius [banjing-keli_rdmin]

cmat default model linear method deform emod 100e6 kratio 1.5 property fric 0.1

ball attribute density 2.7e3 damp 0.7cycle 2000 calm 50

solve

ball delete range cylinder end1 0 0 [sample_hight*0.5] ...    end2 0 0 [-sample_hight*0.5] radius [banjing] not

save sample

2、加叠环

第二步算得上是这个技术的核心部分了，在PFC中加叠环。PFC中是没有叠环这个元素的，所以我们需要在外部绘制一个导入进来。这里是用3DMax绘制的空心圆环，绘制的时候并没有注意位置和尺寸，只是限定了外环直径是内环的2倍，这样方便调整。（圆环形状文件在文末）

首先指定了一下希望式样范围内圆环的个数为10，为了防止式样漏出，在上下多生成了8个圆环。

圆环的图形变化为：

移动到原点---moveToOrigin
在原点进行缩放---scaleInOrigin
移动到式样的顶端，墙体导入----moveBy
每次都往下移动一个圆环的高度，再进行墙体导入

restore sample[num_yuanhuan=10]

wall deletedomain extent [-banjing*5.0] [banjing*5.0] [-banjing*5.0] [banjing*5.0] ...    [-sample_hight*0.5*3] [sample_hight*0.5*3]geometry import yuanhuan.stl call geo_tools@moveToOrigin("yuanhuan")[x_fac=banjing*4.0/(x_pos_max-x_pos_min)][y_fac=banjing*4.0/(y_pos_max-y_pos_min)][z_fac=sample_hight/float(num_yuanhuan)/(z_pos_max-z_pos_min)]@scaleInOrigin("yuanhuan",@x_fac,@y_fac,@z_fac)@moveBy("yuanhuan",0,0,[sample_hight*0.5 sample_hight*3.5/float(num_yuanhuan)])def shengchengqiang(number)    loop n(1,number)        wallname=string.build('yuanhuan_%1',n)                command            wall import geometry yuanhuan name @wallname         endcommand        moveBy("yuanhuan",0,0,[-sample_hight/float(num_yuanhuan)])    endloopend@shengchengqiang([num_yuanhuan 8])save yuanhuan

geo_tools代码

def get_min_and_max(geo)     local gs = geom.set.find(geo)    x_pos_min=1e100    y_pos_min=1e100    z_pos_min=1e100    x_pos_max=-1e100    y_pos_max=-1e100    z_pos_max=-1e100    loop foreach local gn geom.node.list(gs)        if geom.node.pos.x(gn)>x_pos_max then            x_pos_max=geom.node.pos.x(gn)        endif         if geom.node.pos.y(gn)>y_pos_max then            y_pos_max=geom.node.pos.y(gn)        endif         if geom.node.pos.z(gn)>z_pos_max then            z_pos_max=geom.node.pos.z(gn)        endif         if geom.node.pos.x(gn)<x_pos_min then            x_pos_min=geom.node.pos.x(gn)        endif         if geom.node.pos.y(gn)<y_pos_min then            y_pos_min=geom.node.pos.y(gn)        endif         if geom.node.pos.z(gn)<z_pos_min then            z_pos_min=geom.node.pos.z(gn)        endif    endloopend

def moveToOrigin(geo)    local gs = geom.set.find(geo)    get_min_and_max(geo)    x_zhong=(x_pos_max x_pos_min)/2.0    y_zhong=(y_pos_max y_pos_min)/2.0    z_zhong=(z_pos_max z_pos_min)/2.0    loop foreach local gn geom.node.list(gs)        geom.node.pos.x(gn)=geom.node.pos.x(gn)-x_zhong        geom.node.pos.y(gn)=geom.node.pos.y(gn)-y_zhong        geom.node.pos.z(gn)=geom.node.pos.z(gn)-z_zhong    endloopend

def scaleInOrigin(geo,x_fac,y_fac,z_fac)    local gs = geom.set.find(geo)    loop foreach local gn geom.node.list(gs)        geom.node.pos.x(gn)=geom.node.pos.x(gn)*x_fac        geom.node.pos.y(gn)=geom.node.pos.y(gn)*y_fac        geom.node.pos.z(gn)=geom.node.pos.z(gn)*z_fac    endloopend

def Tran_y_z(geo)    local gs = geom.set.find(geo)    loop foreach local gn geom.node.list(gs)       y_temp=geom.node.pos.y(gn)        geom.node.pos.y(gn)=geom.node.pos.z(gn)        geom.node.pos.z(gn)=y_temp    endloopend

def moveBy(geo,x_place,y_place,z_place)    local gs = geom.set.find(geo)    loop foreach local gn geom.node.list(gs)        geom.node.pos.x(gn)=geom.node.pos.x(gn) x_place        geom.node.pos.y(gn)=geom.node.pos.y(gn) y_place        geom.node.pos.z(gn)=geom.node.pos.z(gn) z_place    endloopend

这部分结束后，模型的状态为：

一个颜色就是一个圆环。

3、加上下板

我们还需要上下的墙体来控制围压，这里用的是两个disk。

restore yuanhuanwall generate disk position 0 0 [sample_hight*0.5 sample_hight/float(num_yuanhuan)] dip 0  ddir 0 radius [banjing*1.5]wall generate disk position 0 0 [-sample_hight*0.5-sample_hight/float(num_yuanhuan)] dip 0  ddir 0 radius [banjing*1.5]

save shangxia

到这一步为止，我们的模型算是准备好了。中间我是用clipbox切了一部分墙体出来，然后再添加一个墙体透明度调高即可。

4、加围压

按理说这一步我们的模型参数就应该定下，但是由于离散元的一些理论缺陷和颗粒数的限制，我们决定在围压的时候依然不给摩擦系数，在加载的时候给摩擦系数。

在这里我们只对上面的墙体加伺服。在这里我们伺服部分写的比较少，也比较容易理解伺服的本质了。



 restore shangxiadef wp_init    wpUp=wall.find(19)    wpDown=wall.find(20)end@wp_initwall property fric 0

define compute_wallstress  wszz = -wall.force.contact.z(wpUp) /(math.pi*(banjing)^2.0)end[tzz = -100e3]define servo_walls    gz=2e-5    compute_wallstress    zvel = gz*(wszz- tzz)    wall.vel.z(wpUp) = - zvelendset fish callback  1.0 @servo_walls

history @wszz cycle 2000solvesave weiya

由于式样比较松，所以伺服情况是一个先增大后减小到目标应力然后平衡的过程。

5、加载

最后一步就是加载了，我们把颗粒属性的摩擦系数都给上，将上部墙体的摩擦系数也给上。

之后的细节方面不去讲解，主要功能是识别上下墙体所在圆环的序号，这个理解起来应该不会特别难，进行一个位置的筛选即可。

然后我们给上下墙体之间的圆环加上一个线性变化的速度，上下墙体的水平方向无速度，上面墙体的竖向伺服依然开着。

中间生成一个比较大的测量圆，这里所有的变量都通过测量圆测得。



restore weiya

ball property fric 0.5ball attribute displacement multiply 0[Vel=banjing*0.1]def add_prop    loop n(1,18)        wp=wall.find(n)        wall.prop(wp,"fric")=0    endloop    wall.prop(wpUp,"fric")=0.5    wall.prop(wpDown,"fric")=0.5end@add_prop

[qiang_z_width=0]@get_min_and_max("yuanhuan")[yuanhuan_z_width=(z_pos_max-z_pos_min)/2.0];====================================================================

;======================================================================def wall_group    wall_z_max_xia=wall.pos.z(wpDown) qiang_z_width    q=1    wpp_min=1e100    loop while wpp_min >  wall_z_max_xia        wpp=wall.find(q)        wpp_min=wall.pos.z(wpp)-yuanhuan_z_width         z_wall_min=q        q=q 1    endloop    wall_z_min_shang=wall.pos.z(wpUp)-qiang_z_width    q=num_yuanhuan 8    wpp_min=-1e100    loop while wpp_min <  wall_z_min_shang         wpp=wall.find(q)        wpp_min=wall.pos.z(wpp) yuanhuan_z_width        z_wall_max=q        q=q-1    endloop    end@wall_groupset mech age 0.0





def wall_vel    n=z_wall_min    loop while n >= z_wall_max        x_vel=(n-z_wall_min)*Vel/float(z_wall_max-z_wall_min)        wp=wall.find(n)        wall.vel.x(wp)=x_vel        n=n-1    endloop    loop m(1,z_wall_max)        wp=wall.find(m)        wall.vel.x(wp)=Vel    endloopend@wall_vel

measure create id 1 position 0 0 0 radius [sample_hight*0.4][mp=measure.find(1)]

def stress_strain_mea    mea_stress_xx=measure.stress.xx(mp)    mea_stress_yy=measure.stress.yy(mp)    mea_stress_zz=measure.stress.zz(mp)    mea_stress_xy=measure.stress.xy(mp)    mea_stress_yz=measure.stress.yz(mp)    mea_stress_xz=measure.stress.xz(mp)    mea_strain_xz =measure.strainrate.xz(mp)*global.timestep     stress_prin=tensor.prin(measure.stress(mp))    sigm1=comp.x(stress_prin)    sigm2=comp.y(stress_prin)    sigm3=comp.z(stress_prin)    P=(sigm1 sigm2 sigm3)/3.0    q=math.sqrt(tensor.j2(measure.stress(mp)))  endset fish callback -1.0 @stress_strain_mea

history deletehistory ncycle 200 history id 1 @mea_strain_xzhistory id 11 @mea_stress_xxhistory id 12 @mea_stress_yyhistory id 13 @mea_stress_zzhistory id 14 @mea_stress_xyhistory id 15 @mea_stress_yzhistory id 16 @mea_stress_xz

history id 21 @sigm1history id 22 @sigm2history id 23 @sigm3history id 24 @Phistory id 25 @q

def stop_me    if mea_strain_xz<-0.3 then        stop_me=1    endifend

solve fishhalt @stop_me

save simple_shear

加载后的状态为：

可以看到的是整个边界还是如我们所想的那样，式样也是按照设计的方式进行变形。

力链图也是符合常理的

测量球的位置如图所示

三、应力分析

各种研究表明，单剪实验中式样的应力分布是不均匀的，但是中间核心部分的式样是符合理论解的。所以相对于边界上面的力，直接使用测量圆测得的力可能更加有实际意义。这里的颗粒数还是不够多，所以结果不是特别的理想。

首先我们看一下三个方向的应力

可以发现的是三个方向的应力都有不同程度的增大，z向应力的增大是最多的。但是比较有意思的是我们伺服是开着的，所以边界上的竖向力是不变的。结合力链图是可以解释这个现象的。

式样内部的力链从均匀到不均匀变化，有的地方应力减小有的地方增大，我们可以发现的是在A区力链是减小的，在B区是增大的。C区是B区联系的部分，会在近乎竖直方向形成比较强的力链，所以竖直向的力是在变大的。

我们再看一下XZ向的应变图:

这个是比较好理解的，一直在变大。

我们提取出来的三个主应力都有不同程度的变大，和xx yy zz向的正应力还是比较类似的。

我们可以通过应力十字架来分析单剪实验中的应力偏转现象。我们可以看到在加载前十字架大小比较均匀的，因为式样的不均匀性，略微有些偏转。注意应力方向和大小要结合在一起分析，比如1应力是1，2应力是0.5；和1应力是1000，2应力是0.5，这两个方向是一样的，性质可就差多了。所以除了方向，我们一般也要观察一下I2或者J2。加载后大主应力为加载方向，小主应力在yy向，中主应力在zz向，这也是符合实际的。