PFC双轴案例

已于 2023-08-06 21:21:16 修改
阅读量4.3k 收藏 40
点赞数 10
分类专栏： PFC 文章标签： 开发语言
于 2023-04-28 20:32:23 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_52546971/article/details/130433099
版权

目录

双轴压缩

伺服原理

参数化成样

应力计算并预压

加围压

重定义函数与加载

双轴压缩

伺服原理

        通过位移给墙体施加一个力。（墙体没有质量，只能根据重叠量来受力）

        使用墙给一个颗粒施加力

new

domain extent -10 10

ball create position 0 0 radius 1.0

ball attribute density 2e3 damp 0.7
ball fix vel spin
cmat default model linear property kn 1e7

wall create vertices -1 1 1 1
def cal_chongdie
    chongdieliang=1e4/1e7
end
@cal_chongdie
；计算重叠量
wall attribute yvel -1.0 range id 1 
；给一个y方向的速度

set timestep fix 1e-5
solve time [chongdieliang]

        使用墙给多个颗粒施加力

new

domain extent -10 10

ball create position -2 0 radius 0.5
ball create position -1 0 radius 0.5
ball create position -0 0 radius 0.5
ball create position 1 0 radius 0.5
ball create position 2 0 radius 0.5

ball attribute density 2e3 damp 0.7
ball fix vel spin
cmat default model linear property kn 1e7

wall create vertices -3 0.5 3 0.5
def cal_chongdie
    chongdieliang=1e4/(1e7*5)
end
@cal_chongdie
；重叠量为力除以刚度，墙体的力为1e4N
；一般不给试样施加力，施加应力
wall attribute yvel -1.0 range id 1 

set timestep fix 1e-5
solve time [chongdieliang]

        对上面的加载进行优化，使得cycle 1 步就能得到目标值。

new

domain extent -10 10

ball create position -2 0 radius 0.5
ball create position -1 0 radius 0.5
ball create position -0 0 radius 0.5
ball create position 1 0 radius 0.5
ball create position 2 0 radius 0.5

ball attribute density 2e3 damp 0.7
ball fix vel spin
cmat default model linear property kn 1e7

wall create vertices -3 0.5 3 0.5

def cal_chongdie
    chongdieliang=(1e4-0)/(1e7*5)
    sudu=chongdieliang/global.timestep    
    command        
        wall attribute yvel [-sudu] range id 1 
    endcommand
end
；对墙底施加速度
set fish callback 1.0 @cal_chongdie
；在时步计算之后进行回调
cycle 1

        球体串联情况

new

domain extent -10 10

ball create position 0 0 radius 1.0
ball create position 0 2 radius 1.0
ball attribute density 2e3 damp 0.7
ball fix vel spin range id 1
cmat default model linear property kn 1e7

wall create vertices -1 3 1 3

;k=k1*k2/(k1+k2)  
;k=k1/2
；串联时两个球体刚度计算
def cal_chongdie
    chongdieliang=1e4/(1e7*0.5)
end
@cal_chongdie

wall attribute yvel -1.0 range id 1 

set timestep fix 1e-5
solve time [chongdieliang]
wall attribute yvel -0.0 range id 1 
solve

        三角形分布球体伺服

        whilestepping，stop_me

new

domain extent -10 10

ball create position -1 0 radius 1.0
ball create position 1 0 radius 1.0

ball create position 0 [math.sqrt(3)] radius 1.0
；生成三个颗粒，根号3

ball attribute density 2e3 damp 0.7
ball fix vel spin range id 1
ball fix vel spin range id 2
cmat default model linear property kn 1e7

wall create vertices -1 [math.sqrt(3)+1] 1 [math.sqrt(3)+1]

;1、力小的时候，下压
;2、力大了，上浮
[wp=wall.find(1)]
；定义墙体指针
def get_force
    whilestepping
    wall_force=wall.force.contact.y(wp)
end
；whilestepping的效果和回调函数一样，即函数结果将会随着计算的更新而更新

def apply_vel
    whilestepping
    if wall_force<1e4 then
        wall.vel.y(wp)=-1.0
    else
        wall.vel.y(wp)=1.0
    endif    
end
；wall_force向上，为正值
；力小的话则墙体下压，否则墙体上浮

[stop_me=0]
def stop_me
    if math.abs(wall_force-1e4)<1 then
        stop_me=1
    endif
end
；stop_me函数判断截止条件，被solve调用，solve fishhalt @stop_me调用了停止函数
；当值为0时继续运行，当为1时停止运行


set timestep fix 1e-5
；时步不能太大，否则容易使得力来回变动
;solve fishhalt @stop_me  
cycle 1
；使用cycle 1 来打破平衡
solve
；两个要求，一施加力达到要求，二达到平衡

;cycle 2.2w

        伺服优化（速度+伺服系数）

new

domain extent -10 10

ball create position -1 0 radius 1.0
ball create position 1 0 radius 1.0

ball create position 0 [math.sqrt(3)] radius 1.0

ball attribute density 2e3 damp 0.7
ball fix vel spin range id 1
ball fix vel spin range id 2
cmat default model linear property kn 1e7

wall create vertices -1 [math.sqrt(3)+1] 1 [math.sqrt(3)+1]

;1、力小的时候，下压
;2、力大了，上浮
[wp=wall.find(1)]
def get_force
    whilestepping
    wall_force=wall.force.contact.y(wp)
end

[g=5e-4]
def apply_vel
    whilestepping
    vel=g*math.abs(wall_force-1e4) ；速度=伺服系数乘以力的差值
    if wall_force<1e4 then
        wall.vel.y(wp)=-vel
    else
        wall.vel.y(wp)=vel
    endif    
end

[stop_me=0]
def stop_me
    if math.abs(wall_force-1e4)<1 then
        stop_me=1
    endif
end


set timestep fix 1e-5
;solve fishhalt @stop_me 
cycle 1
solve

;g=1e-3 cycle 5k8
;g=1e-2 波动 不平衡 伺服系数g不能太大
;g=5e-3 cycle 1w4
;g=2e-3 cycle 2w2
;g=5e-4 cycle 1w4
；伺服系数需要试算

        试算伺服系数

new

domain extent -10 10

ball create position -1 0 radius 1.0
ball create position 1 0 radius 1.0

ball create position 0 [math.sqrt(3)] radius 1.0

ball attribute density 2e3 damp 0.7
ball fix vel spin range id 1
ball fix vel spin range id 2
cmat default model linear property kn 1e7

wall create vertices -1 [math.sqrt(3)+1] 1 [math.sqrt(3)+1]

;1、力小的时候，下压
;2、力大了，上浮
[wp=wall.find(1)]
def get_force
    whilestepping
    wall_force=wall.force.contact.y(wp)
end

；重叠量=力除以刚度
；速度=重叠量除以时步
[servo_factor=0.8]
def get_g
    ct_wp=contact.find("ball-facet",1)   
    gangdu=contact.prop(ct_wp,"kn")
    g=1.0*servo_factor/(gangdu*global.timestep)  
end
set fish callback 1.0 @get_g

def apply_vel
    whilestepping
    vel=g*math.abs(wall_force-1e4)
    if wall_force<1e4 then
        wall.vel.y(wp)=-vel
    else
        wall.vel.y(wp)=vel
    endif    
end

[stop_me=0]
def stop_me
    if math.abs(wall_force-1e4)<1 then
        stop_me=1
    endif
end


set timestep fix 1e-5
;solve fishhalt @stop_me 
cycle 10
solve

参数化成样

        加载板和试样之间的摩擦力为0，而颗粒与颗粒之间为0.5。（两者摩擦力取最小值）

new
def par
    width=0.4
    height=width*2
    
    rdmin=0.006
    rdmax=0.009
    
    poro=0.12
end
@par
；给出需要的参数

domain extent [-width*2.0] [width*2.0] ...
        [-height*2.0] [height*2.0]
set random 10001
；给定一个random的数值，使得每次运行得到的结果相同

wall generate box [-width*0.5] [width*0.5] ...
            [-height*0.5] [height*0.5] expand 1.5
；把边缘坐标进行放大（把加载板稍微加大一点）
ball distribute porosity @poro radius [rdmin] [rdmax] box ...
    [-width*0.5] [width*0.5] [-height*0.5] [height*0.5] 
ball attribute density 2e3 damp 0.7
cmat default model linear method deform emod 100e6 kratio 1.5
；method可以用来定义一些宏观的参数，但emod即弹性模量，但是弹性模量也会自动被计算为刚度
cycle 2000 calm 50
；让cycle尽量大，使得颗粒散开
ball property fric 0.5
；最好给出fric
solve
save sample

应力计算并预压

        预压模拟的是土体的先前条件。伺服频率

        对于墙体和颗粒接触，模量乘以2对应kn的数值。

restore sample

[txx=1e4]
[tyy=1e4]
def wall_ini
    wpup=wall.find(3)
    wpdown=wall.find(1)
    wpleft=wall.find(4)
    wpright=wall.find(2)
end
@wall_ini
；使用指针找到墙体

def computer_chicun
    wlx=wall.pos.x(wpright)-wall.pos.x(wpleft)
    wly=wall.pos.y(wpup)-wall.pos.y(wpdown)
end
；计算尺寸
def computer_stress
    computer_chicun
    wyss=0.5*(wall.force.contact.y(wpdown)-wall.force.contact.y(wpup))/wlx
    wxss=0.5*(wall.force.contact.x(wpleft)-wall.force.contact.x(wpright))/wly
end
；计算应力，两侧应力的平均值

；伺服
[servo_factor=0.8]
def get_g

    zongKNY=100e6*2.0  ；假设初始时有一定的约束，防止颗粒过于松散，kn为0
    loop foreach ct wall.contactmap(wpup) ；遍历，ct为临时变量，它会把墙体接触依次定义
        zongKNY+=contact.prop(ct,"kn")   ；总的竖向刚度，
    endloop
    loop foreach ct wall.contactmap(wpdown)
        zongKNY+=contact.prop(ct,"kn")
    endloop
    
    zongKNX=100e6*2.0
    loop foreach ct wall.contactmap(wpleft)
        zongKNX+=contact.prop(ct,"kn")
    endloop
    loop foreach ct wall.contactmap(wpright)
        zongKNX+=contact.prop(ct,"kn")
    endloop
    
    gx=1.0*servo_factor*wly/(zongKNX*global.timestep) 
    ；wly表示对应的面积，把应力转化为力
    gy=1.0*servo_factor*wlx/(zongKNY*global.timestep)
    ；当试样特别松散的时候，刚度可能出现为0的情况
    ；时步正比于(质量/刚度)  
end

[sevro_freq=100] ；定义伺服频率，以步数作为跨度,隔一定步数计算函数
[time_record=global.step-1]
def sevro_walls
    computer_stress
    if global.step>time_record then
        get_g
        time_record=global.step+sevro_freq   ；每100步计算一次get_g
    endif    
    yvel=gy*math.abs(math.abs(wyss)-tyy) ；竖向速度
    if math.abs(wyss)<tyy then
        wall.vel.y(wpup)=-yvel
        wall.vel.y(wpdown)=yvel
    else
        wall.vel.y(wpup)=yvel
        wall.vel.y(wpdown)=-yvel
    endif  
    xvel=gx*math.abs(math.abs(wxss)-txx)    ；横向速度
    if math.abs(wxss)<txx then
        wall.vel.x(wpleft)=xvel
        wall.vel.x(wpright)=-xvel
    else
        wall.vel.x(wpleft)=-xvel
        wall.vel.x(wpright)=xvel
    endif  
end

set fish callback -1.0 @sevro_walls 
history id 1 @wxss
history id 2 @wyss
；记录x和y方向应力的变化
cycle 10
；多运行几下，打破平衡
solve
save yuya

加围压

        利用callback和fish函数。

restore yuya
；采用伺服原理
[txx=100e3]
[tyy=100e3]
；设置为100kpa
；改变txx和tyy会直接代入yuya命令部分的函数进行计算
；因为callback存在，所以运行过程中会不断更新数值

cycle 1
solve

save weiya

重定义函数与加载

        对于预压部分的伺服函数进行重定义，删掉竖向伺服，添加竖向墙体速度。同时设定hist并添加加载截止条件。

restore weiya

ball attribute displacement multiply 0
；将位移进行清零
def sevro_walls
    computer_stress
    if global.step>time_record then
        get_g
        time_record=global.step+sevro_freq
    endif  
    ；获得g这个伺服系数
    xvel=gx*math.abs(math.abs(wxss)-txx)
    if math.abs(wxss)<txx then
        wall.vel.x(wpleft)=xvel
        wall.vel.x(wpright)=-xvel
    else
        wall.vel.x(wpleft)=-xvel
        wall.vel.x(wpright)=xvel
    endif  
    ；关掉竖向伺服，给竖向一个速度，墙体的下压与上浮
end
[strainRate=1e-2]
wall attribute yvel [strainRate*wly] range id 1
wall attribute yvel [-strainRate*wly] range id 3
；用应变率×尺寸作为速度

[Iy0=wly]
[Ix0=wlx]
；记录初始的宽度和高度
def computer_strain
    weyy=(Iy0-wly)/Iy0
    wexx=(Ix0-wlx)/Ix0
    wevol=weyy+wexx
；体变等于两个方向应变相加
end
set fish callback -1.0 @computer_strain
history delete 
；删除之前的一些历史记录然后重新定义
history id 1 @wxss
history id 2 @wyss
history id 3 @weyy
history id 4 @wexx
history id 5 @wevol

[stop_me=0]
def stop_me
    if weyy>0.2 then
        stop_me=1
    endif
end
solve fishhalt @stop_me
；人工定义截止条件，当应变大于0.2时就停止
PFC实例
kingma6666的博客
02-14 750
PFC功率因数校正方法
PFC实例源码
03-14
PFC实例代码例子，全面，对想了解PFC的朋友会有很大帮助，容易撑握。
PFC6.0——第五章 漏斗案例
最新发布
Cx330_wjy的博客
03-27 366
有些情况下会先判断当前的稳定状态，导致直接结束，一般前面加 model cycle 1，相当于晃动一下，再进行收敛计算。漏斗案例常应用于农业、机械等方面，通过此案例，可以进一步感受到项目模块化建模的重要性，并学习到多种的墙体生成方式、颗粒分组。离散元研究破坏，一般先染色，呈现更好的可视化效果。指定生成一固定box墙体。
PFC电感计算公式.PDF
09-14
开关电源的APFC电路中电感，是PFC电路中的关键点,资料分析了PFC电路原理和设计步骤.
PFC电感计算.ppt
05-10
针对PFC电路进行系统的讲解，分类单相BOOST、交错并联型，CCM、BCM、DCM，以及不同的控制方式（平均电流型，峰值电流型，滞环控制型），给出了PFC电感的计算公式和推导，较为详尽，对于学习PFC电路的工程师值得学习。
PFC2D5.0_双轴试验
qq_56811637的博客
01-18 1549
在加载过程中，墙体向土体方向移动，wlx、wly在一直发生变化，为达到指定的围压，wall servo命令指定的力必须随着wlx、wly的变化而调整，调用伺服函数和callback命令就可以实现这个目的。iso_compression函数的调用是判断围压是否达到500kPa（误差0.005），若达到，则结束。（本文施加围压500kPa）
PFC岩石双轴
weixin_52546971的博客
05-01 3514
添加胶结，胶结间距一般取rdmin*0.2。间距越大，获得的胶结越多；间距过小可能会出现悬浮颗粒。x方向和y方向分别施加0.5MPa围压。怎么看是否有胶结（监测）——重点是移除墙体伺服。可以查看胶结破坏的颗粒，会存在明显裂纹。在每个裂纹地方，都有一个dfn。（张拉裂纹或者剪切裂纹）和。：删除墙体之后仍然有接触。——调试显示应力应变曲线。：删除墙体后接触消失。
PFC 5.0岩土工程中单轴、双轴、三轴实验的颗粒流离散元算例代码包,PFC5.0颗粒流离散元算例代码集：岩土单轴、双轴及三轴测试案例全解析,PFC5.0颗粒流离散元【岩土-单轴-双轴-三轴】算例代码
02-16
PFC 5.0岩土工程中单轴、双轴、三轴实验的颗粒流离散元算例代码包,PFC5.0颗粒流离散元算例代码集：岩土单轴、双轴及三轴测试案例全解析,PFC5.0颗粒流离散元【岩土-单轴-双轴-三轴】算例代码，包含2D和3D的案例，所有...
PFC5.0颗粒流离散元【岩土-单轴-双轴-三轴】算例代码，包含2D和3D的案例，所有案例打包 本案例提供参考，可以自行修改参数
10-08
在当今工程领域，特别是在岩土力学的研究和应用中，PFC5.0颗粒流离散元软件已经成为...通过对提供的2D和3D案例代码的学习和修改，用户将能够更加深入地掌握PFC5.0的应用技巧，并将其运用到更加广泛的岩土工程实践中去。
PFC2D5.0双轴试验算例：松散与密实试样模拟及参数调整指南,PFC2D5.0颗粒流离散元模型双轴试验参数化算例解析与学习,PFC2D5.0颗粒流离散元双轴试验算例代码，算例中有松散和密实试样的
03-02
PFC2D5.0双轴试验算例：松散与密实试样模拟及参数调整指南,PFC2D5.0颗粒流离散元模型双轴试验参数化算例解析与学习,PFC2D5.0颗粒流离散元【双轴试验】算例代码，算例中有松散和密实试样的双轴试验。 可以自己更改...
"PFC5.0颗粒流离散元模拟软件：岩土单轴、双轴、三轴算例代码集，2D与3D案例全覆盖，含详细参数修改及二次加工指南，部分代码带注释便于学习参考",PFC5.0颗粒流离散元【岩土-单轴-双轴-三轴】
01-19
岩土单轴、双轴、三轴算例代码集，2D与3D案例全覆盖，含详细参数修改及二次加工指南，部分代码带注释便于学习参考",PFC5.0颗粒流离散元【岩土-单轴-双轴-三轴】算例代码，包含2D和3D的案例，所有案例打包 ...
不规则颗粒的双轴压缩试验.zip
05-06
在PFC中运用FISH语言编写程序，模拟双轴压缩模拟试验，改变的参数可以为：块石的形状、块石方向等
PFC5.0岩土单轴_双轴_三轴算例.
08-05
PFC5.0岩土单轴_双轴_三轴算例.（各种算例 2d 3d 代码）
rouxing_柔性伺服控制_pfc2d代码大全_PFC围压_
09-30
PFC2D柔性伺服代码，PFC模型的围压伺服代码 为最新5.0版本代码
PFC2D直剪试验命令流
04-04
PFC2D的直剪试验命令流 无需调用fishtank
双轴代码转换
weixin_52546971的博客
05-12 777
除此之外，因为2D的y变为了现在3D的z，一些对应的细节也需要更改。更新墙体指针，并重新定义y和z方向的尺寸/应力。关于模型的操作在前面需要加上model。→找不到对应字符的情况应该到手册中查找。cmat需要在前面加上contact。字符串需要加上双引号。
PFC6.0documentation_PFC
weixin_52546971的博客
04-11 7607
对于每一个要分配的球/块，根据用户的输入选择一个随机的大小和位置（仅对于块来说是方向），然后加入预期的球/块，直到达到指定的固体分数。一旦存在一个接触，并且窄相接触检测确定了碎片重叠（即计算了穿透矢量，并且间隙为负），计算力所需的其余信息就是接触位置。接触力学体现在采用软接触方法的粒子相互作用定律中，对于这种方法，所有的变形都发生在刚体之间的接触处。牛顿第二定律用于确定每个粒子的运动，这些运动是由作用在它身上的接触力和体力引起的，而力-位移定律则用于更新由每个接触处的相对运动引起的接触力。
PFC 3D 3.0三轴压缩建模
11-02
：颗粒流问题是众多工农业生产领域中的常见问题，然而人 们至今对颗粒流机理认识不深。目前一种基于离散元的数值模拟方 法被用以模拟和分析颗粒流问题，研究实践表明该方法是散体力学 分析的一种有效的工具，PFC 3d离散单元法在岩土工程真三轴中使用广泛。
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值