【PFC6.0】子*弹穿墙模拟

lobby 2022年9月13日浏览：1899 评论：3 收藏：2

0 引言

最近刚好刷到《战狼1》的短视频，里面有一个经典的子*弹穿墙镜头，想着或许可以使用离散元的方法进行模拟。对于弹道分析、穿墙速度损失等方面应该会有一定的指导意义。案例中使用到了很多技术方法，比如子*弹的形状是用圆锥二次化后生成的，墙体的边界处理等。

【PFC6.0】子*弹穿墙模拟的图1

1 成样

破碎方面一般只在被击中的砖周边，所以我们只建了三行三列的砖墙，其余的边界用有限元的zone模拟。砖头采用240*115*90mm尺寸，粒径在2-3mm之间，生成了20w左右的颗粒。由于后面是想用力清零的方式来形成胶结物，所以这里对模型的平衡度要求不高，ratio到1e-2停止。

model newdef zhuan_par    zhuan_length=240e-3    zhuan_width=115e-3    zhuan_height=90e-3        rdmin=2e-3    rdmax=3e-3end@zhuan_parmodel domain extent [-zhuan_length*2] [zhuan_length*2] 

wall generate  box [-zhuan_length*1.5] [zhuan_length*1.5] ...            [-zhuan_width*0.5]  [zhuan_width*0.5] [-zhuan_height*1.5] [zhuan_height*1.5]            ball distribute porosity 0.3 radius @rdmin @rdmax  box [-zhuan_length*1.5] [zhuan_length*1.5] ...            [-zhuan_width*0.5]  [zhuan_width*0.5] [-zhuan_height*1.5] [zhuan_height*1.5]ball attribute density 3.8e3 damp 0.7contact cmat default type ball-ball model linear method deform emod 1e9 kratio 1.5contact cmat default type ball-facet model linear method deform emod 10e9 kratio 1.5

model mechanic timestep scalemodel cycle 1model solve ratio-average 1e-2

model save "sample"

2 砖墙边界

模拟混凝土墙的话，使用一整个式样即可。对于砖墙的话一般是错开排列的，使用水泥等进行粘结。这里首先建立每块砖边界上的形状，然后导入墙体进行边界平整，然后赋予参数成样。

model restore "sample"

geometry set "fenjie"

geometry polygon create  by-positions  [-zhuan_length*1.5] [-zhuan_width*0.5] [zhuan_height*0.5] ...                        [-zhuan_length*1.5] [zhuan_width*0.5] [zhuan_height*0.5] ...                        [zhuan_length*1.5] [zhuan_width*0.5] [zhuan_height*0.5] ...                        [zhuan_length*1.5] [-zhuan_width*0.5] [zhuan_height*0.5] geometry polygon create  by-positions [-zhuan_length*1.5] [-zhuan_width*0.5] [-zhuan_height*0.5] ...                        [-zhuan_length*1.5] [zhuan_width*0.5] [-zhuan_height*0.5] ...                        [zhuan_length*1.5] [zhuan_width*0.5] [-zhuan_height*0.5] ...                        [zhuan_length*1.5] [-zhuan_width*0.5] [-zhuan_height*0.5] 



def gen_on_layer(x_pos_qi,z_pos_qi,num)    loop n(1,num)        x_pos=x_pos_qi+(n-1)*zhuan_length        command           geometry polygon create  by-positions @x_pos [-zhuan_width*0.5] [z_pos_qi] ...                        @x_pos [zhuan_width*0.5] [z_pos_qi] ...                        @x_pos [zhuan_width*0.5] [z_pos_qi+zhuan_height] ...                        @x_pos [-zhuan_width*0.5] [z_pos_qi+zhuan_height]        endcommand    endloopend@gen_on_layer([-zhuan_length*1.0],[zhuan_height*0.5],3)@gen_on_layer([-zhuan_length*0.5],[-zhuan_height*0.5],2)@gen_on_layer([-zhuan_length*1.0],[-zhuan_height*1.5],3)



def create_wall_from_geo    geo_pt=geom.set.find("fenjie")    loop foreach pg geom.poly.list(geo_pt)        n1= geom.poly.node(pg,1)        n2= geom.poly.node(pg,2)        n3= geom.poly.node(pg,3)        n4= geom.poly.node(pg,4)        command            wall generate  group "fenjie" polygon [geom.node.pos(n1)] ...                                [geom.node.pos(n2)] ...                                [geom.node.pos(n3)] ...                                [geom.node.pos(n4)]        endcommand    endloopend@create_wall_from_geomodel cycle 1model solve ratio-average 1e-2

model save "wall_ball"

3 加胶结

这里为了弱化边界的影响，使用了zone来模拟周边的墙体来传递应力应变。给了砖头比较强的接触，边界给了弱一倍的参数，这里是随便给的，如果是专门做的话，需要进行参数标定。

model restore "wall_ball"[emod=1e9][pb_coh=3e6]



wall delete 

zone create brick size (4,1,11) point 0 ([zhuan_length*1.5],[-zhuan_Width*0.5],[-zhuan_height*5.5])  ...                                    point 1 ([zhuan_length*5.5],[-zhuan_Width*0.5],[-zhuan_height*5.5]) ...                                  point 2 ([zhuan_length*1.5],[zhuan_Width*0.5],[-zhuan_height*5.5]) ...                                  point 3 ([zhuan_length*1.5],[-zhuan_Width*0.5],[zhuan_height*5.5])zone create brick size (4,1,11) point 0 ([-zhuan_length*5.5],[-zhuan_Width*0.5],[-zhuan_height*5.5])  ...                                    point 1 ([-zhuan_length*1.5],[-zhuan_Width*0.5],[-zhuan_height*5.5]) ...                                  point 2 ([-zhuan_length*5.5],[zhuan_Width*0.5],[-zhuan_height*5.5]) ...                                  point 3 ([-zhuan_length*5.5],[-zhuan_Width*0.5],[zhuan_height*5.5])zone create brick size (3,1,4) point 0 ([-zhuan_length*1.5],[-zhuan_Width*0.5],[zhuan_height*1.5])  ...                                    point 1 ([zhuan_length*1.5],[-zhuan_Width*0.5],[zhuan_height*1.5]) ...                                  point 2 ([-zhuan_length*1.5],[zhuan_Width*0.5],[zhuan_height*1.5]) ...                                  point 3 ([-zhuan_length*1.5],[-zhuan_Width*0.5],[zhuan_height*5.5])zone create brick size (3,1,4) point 0 ([-zhuan_length*1.5],[-zhuan_Width*0.5],[-zhuan_height*5.5])  ...                                    point 1 ([zhuan_length*1.5],[-zhuan_Width*0.5],[-zhuan_height*5.5]) ...                                  point 2 ([-zhuan_length*1.5],[zhuan_Width*0.5],[-zhuan_height*5.5]) ...                                  point 3 ([-zhuan_length*1.5],[-zhuan_Width*0.5],[-zhuan_height*1.5])zone cmodel assign elasticzone property young [emod] poisson 0.25

zone gridpoint fix velocity range pos-x [-zhuan_length*5.6]  [-zhuan_length*5.4]zone gridpoint fix velocity range pos-x [zhuan_length*5.4]  [zhuan_length*5.6]

zone gridpoint fix velocity range pos-z [-zhuan_height*5.6]  [-zhuan_height*5.4]zone gridpoint fix velocity range pos-z [zhuan_height*5.4]  [zhuan_height*5.6]

wall-zone create name 'bianjie1' skip-errors range pos-x [-zhuan_length*1.6] [-zhuan_length*1.4] ...                            pos-z [-zhuan_height*1.6] [zhuan_height*1.6] wall-zone create name 'bianjie2' skip-errors range pos-x [zhuan_length*1.4] [-zhuan_length*1.6] ...                            pos-z [-zhuan_height*1.6] [zhuan_height*1.6] wall-zone create name 'bianjie3' skip-errors range pos-x [-zhuan_length*1.6] [zhuan_length*1.6] ...                            pos-z [zhuan_height*1.4] [zhuan_height*1.6] wall-zone create name 'bianjie4' skip-errors range pos-x [-zhuan_length*1.6] [zhuan_length*1.6] ...                            pos-z [-zhuan_height*1.6] [zhuan_height*1.4] contact cmat default type ball-ball model linearpbond method deform emod @emod kratio 1.5 ...                pb_deform emod @emod kratio 1.5 property pb_coh [pb_coh] pb_ten [pb_coh*3] pb_fa 50contact cmat default type ball-facet model linearpbond method deform emod @emod kratio 1.5 ...                pb_deform emod @emod kratio 1.5 property pb_coh [pb_coh] pb_ten [pb_coh*3] pb_fa 50contact cmat add 1 model linearpbond method deform emod [emod*0.5] kratio 1.5 ...                pb_deform emod [emod*0.5] kratio 1.5 property pb_coh [pb_coh*0.5] pb_ten [pb_coh*1.5] pb_fa 30 range ...                geometry-distance "fenjie" gap [rdmin*0.2]contact cmat applymodel clean

contact method bond gap [rdmax]

contact property lin_force 0.0 0.0 0 lin_mode 1ball attribute force-contact multiply 0.0 moment-contact multiply 0.0 ball attribute damp 0.1model calm

model cycle 1model solve

model save "model_qiang"

到这里为止，砖墙模型就形成了，如下图：

【PFC6.0】子*弹穿墙模拟的图2

胶结情况如图：

【PFC6.0】子*弹穿墙模拟的图3

4 造子*弹

懒得画子*弹图形，并且参数控制也比较容易来研究子*弹形状对结果的影响。这里是首先生成一个圆锥cone作为基本图形，然后将其y坐标二次化，其实是将圆锥给不均匀拉长成一个子*弹的形状。需要注意的cone本身y向是只有一段的，所以需要将其加密才可以获得多段y坐标的cone，然后改变其y坐标。

model restore "model_qiang"

def zidan_par    zidan_zhijing=9.5e-3    zidan_changdu=67.7e-3end@zidan_pargeometry set "zidan"geometry generate cone base 0 [-zhuan_Width*0.5-zidan_zhijing*0.5] 0 ...                axis 0 1 0 ...                height [zidan_zhijing*0.5] ...                radius [zidan_zhijing*0.5] ...                cap true res 0.8

geometry refine 4def extend_geo    geo_zidan=geo.set.find("zidan")    loop foreach nd geo.node.list(geo_zidan)        y_dist=(geo.node.pos.y(nd)+zhuan_Width*0.5)*2000        y_dist_new=(-1/2000.0)* math.abs(y_dist)^(2)-zhuan_Width*0.5        geo.node.pos.y(nd)=y_dist_new    endloopend@extend_geo

def get_min_max    global x_min=1e100    global x_max=-1e100        global y_min=1e100    global y_max=-1e100    local gs = geom.set.find("zidan")    loop foreach local gn geom.node.list(gs)        local pos = geom.node.pos(gn)        if x_min > comp.x(pos)            x_min = comp.x(pos)        endif        if y_min > comp.y(pos)            y_min = comp.y(pos)        endif           if x_max < comp.x(pos)            x_max = comp.x(pos)        endif        if y_max < comp.y(pos)            y_max = comp.y(pos)        endif    endloopend@get_min_maxdef scale_geo    length_now=y_max-y_min    scale_factor=zidan_changdu/length_now    geo_zidan=geo.set.find("zidan")    loop foreach nd geo.node.list(geo_zidan)        geo.node.pos.y(nd)=-zhuan_Width*0.5+(geo.node.pos.y(nd)+zhuan_Width*0.5)*scale_factor    endloopend@scale_geo

rblock create  from-geometry 

model save "add_zidan"

子*弹用rblock来模拟，如图：

【PFC6.0】子*弹穿墙模拟的图4

5 冲击

最后一步就是给子*弹一个速度，来冲击砖墙。在子*弹冲击墙的时候，把时步设置的比较小，防止因为速度过大引起的计算问题。子*弹冲击完后，可以将时步放大，来研究墙体的破碎情况。注意动力问题阻尼要缩小。

model restore "add_zidan"model domain extent [-zhuan_length*2] [zhuan_length*2] condition destroymodel mech time-total 0contact cmat default type ball-rblock model linear method deform emod 10e9 kratio 1.5 model mechanic timestep fix 1e-8ball attribute density 2e3 damp 0.01rblock attribute density 11e3 damp 0.01rblock attribute velocity-y [800] 

[time_record=mech.time.total-100][baocunpinlv=5e-5]def savefile        if mech.time.total-time_record > baocunpinlv then        filename=string.build("jieguo_%1",count)        command            model save @filename                    endcommand        time_record=mech.time.total        count +=1    endif    endfish callback add @savefile -1.0

[rb_zidan=rblock.find(1)]def jiance    whilestepping    zidan_vel_y=rblock.vel.y(rb_zidan)endhistory deletehistory id 1 @zidan_vel_yrblock trace id 1model solve time 0.75e-3[baocunpinlv=1e-2]model mechanic timestep fix 1e-6def jianceendmodel solve time 1

model save "result"