PFC模擬三維單剪實驗

lobby

2021年10月10日 20:53

瀏覽：4327 評論：1 收藏：13

二維單剪是比較簡單的，在雙軸的基礎上改是比較好實現的。但是三維單剪有比較多的細節問題需要解決，而且對于結果的分析都是比較困難的事情。本文主要針對于三維單剪的建模過程和應力分析進行講解。

一、單剪實驗

大家接觸比較多的可能是直剪實驗，上下兩個剪切盒橫向移動，在剪切面上產生剪切力使得式樣發生破壞。而單剪實驗相當于很多個剪切盒堆在一起進行剪切，相對于直剪實驗，更加符合土體的變形特性。

PFC模擬三維單剪實驗的圖1

滑坡體變形與單剪實驗

PFC模擬三維單剪實驗的圖2

直剪實驗變形

PFC模擬三維單剪實驗的圖3

直剪實驗變形

(吳明浙江大學等)

二、單剪實驗建模

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

PFC模擬三維單剪實驗的圖4

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

這部分結束后，模型的狀態為：

PFC模擬三維單剪實驗的圖5

一個顏色就是一個圓環。

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切了一部分墻體出來，然后再添加一個墻體透明度調高即可。

PFC模擬三維單剪實驗的圖6

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

由于式樣比較松，所以伺服情況是一個先增大后減小到目標應力然后平衡的過程。

PFC模擬三維單剪實驗的圖7

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

加載后的狀態為：

PFC模擬三維單剪實驗的圖8

可以看到的是整個邊界還是如我們所想的那樣，式樣也是按照設計的方式進行變形。

PFC模擬三維單剪實驗的圖9

力鏈圖也是符合常理的

PFC模擬三維單剪實驗的圖10

測量球的位置如圖所示

PFC模擬三維單剪實驗的圖11

三、應力分析

各種研究表明，單剪實驗中式樣的應力分布是不均勻的，但是中間核心部分的式樣是符合理論解的。所以相對于邊界上面的力，直接使用測量圓測得的力可能更加有實際意義。這里的顆粒數還是不夠多，所以結果不是特別的理想。

PFC模擬三維單剪實驗的圖12

單剪實驗式樣

(吳明浙江大學等)

首先我們看一下三個方向的應力

PFC模擬三維單剪實驗的圖13

可以發現的是三個方向的應力都有不同程度的增大，z向應力的增大是最多的。但是比較有意思的是我們伺服是開著的，所以邊界上的豎向力是不變的。結合力鏈圖是可以解釋這個現象的。

式樣內部的力鏈從均勻到不均勻變化，有的地方應力減小有的地方增大，我們可以發現的是在A區力鏈是減小的，在B區是增大的。C區是B區聯系的部分，會在近乎豎直方向形成比較強的力鏈，所以豎直向的力是在變大的。

PFC模擬三維單剪實驗的圖14

我們再看一下XZ向的應變圖:

PFC模擬三維單剪實驗的圖15

這個是比較好理解的，一直在變大。

我們提取出來的三個主應力都有不同程度的變大，和xx yy zz向的正應力還是比較類似的。

PFC模擬三維單剪實驗的圖16

我們可以通過應力十字架來分析單剪實驗中的應力偏轉現象。我們可以看到在加載前十字架大小比較均勻的，因為式樣的不均勻性，略微有些偏轉。注意應力方向和大小要結合在一起分析，比如1應力是1，2應力是0.5；和1應力是1000，2應力是0.5，這兩個方向是一樣的，性質可就差多了。所以除了方向，我們一般也要觀察一下I2或者J2。加載后大主應力為加載方向，小主應力在yy向，中主應力在zz向，這也是符合實際的。

PFC模擬三維單剪實驗的圖17