PFC2D学习笔记之边界条件与初始条件.doc

PFC2D学习笔记之边界条件与初始条件边界条件墙边界一般，颗粒组是在一组墙内创建并压缩这些墙也可以当作边界约束，按一定速度运动监测其反力；或保持反力一定，伺服控制其速度但是不能在墙上直接施加力当墙为多段线时，要考虑到墙与球段转角处的接触当转角处为凹时，两段墙可能会同时与一个球接触，但PFC2D只允许每个墙与球的接触最多一个，因此要在此处将墙断开，设置成两个墙；当转角为凸时，两段墙不肯能同时与一个球接触，因此不需要将墙断开墙也可以是圆形或弧或点颗粒边界可以创建一个颗粒串，并使用这些颗粒串作为边界条件固定速度的颗粒边界用FISH函数得到边界颗粒，并固定其平动自由度，然后删除墙，并施加加速度从而可以得到内部颗粒的速度边界颗粒的速度保持初始值不变；该模拟对应与应变控制式试验有两种方法可以得到边界颗粒一种是，遍历每个球的接触列表并探测与墙接触的球，作为边界；另一种是，使用RANGE，将落在一个空间范围内的所有颗粒作为边界颗粒第一种方法的核心程序为：bp=ball_headloopwhilebp#null;scanallballssectioncp=b_clist(bp)loopwhilecp#null;scanball'scontactsifc_nforce(cp)#0.0thenb2=c_ball2(cp)ifpointer_type(b2)=101then;b2isawallb_xfix(bp)=1;fixoriginalballinx,yb_yfix(bp)=1b_color(bp)=1;identifybycolorexitsection;alldoneforthisballend_ifend_ififc_ball1(cp)=bpcp=c_b1clist(cp)elsecp=c_b2clist(cp)end_ifend_loopend_sectionbp=b_next(bp)end_loop第二种方法参考RANGE的用法，较为简单，此处略去。

在颗粒边界上施加力先固定所有边界颗粒，并移除墙，然后CYCLE一次，这样每个边界颗粒上的不平衡力就与之前的墙反力相等或相反接着，用FISH函数在边界颗粒上施加与不平衡力相反的力，这样，边界颗粒受力平衡然后撤除所有颗粒的约束，对边界颗粒施加加速度但是为了防止由于体力的未平衡导致一些边界颗粒移动位置，就需要使用各种稳定方法比如，固定边界颗粒的旋转自由度，防止颗粒从平衡位置旋转开也可以使用混合边界条件与固定速度的颗粒边界不同的是，这种边界条件要使用FISH函数施加力，平衡掉颗粒边界上的不平衡力，其核心程序为：bp=ball_headloopwhilebp#nullifb_xfix(bp)=1b_xfap(bp)=-b_xfob(bp)b_yfap(bp)=-b_yfob(bp)end_ifbp=b_next(bp)end_loop混合边界条件垂直于X轴方向的两个边界固定住并施加速度边界；其他边界施加力当几何边界改变时,这种双轴试验就不能继续了，因为边界力将不可用应该使用FISH定义一个可以随几何边界更新的边界力，参考《双轴试验》部分核心程序与前面相似，不再重复应用命令freexyspinrangex=0.2,9.8,只释放上下边界和内部颗粒的自由度。

初始条件为了得到所需要的初始应力状态，需要调整颗粒组状态，因此必须在一定程度上改变孔隙比，其结果是，孔隙比的很小的变化都会引起应力的很大变化这里的应力是指测量区域内的平均应力，有两种计算方法：一种是，将墙的反力之和除以墙的面积(长度与厚度之积)；另一种是使用一个或多个测量圆获得各向同性应力状态1、生成所需要的紧密颗粒组，用半径扩展法或颗粒排斥法其核心程序参考《颗粒生成细节》2、定义一个测量圆，用于测量区域内平均应力测量圆的大小将会影响计算结果measid=1x5y2.5rad2.03、用内置的MP=FIND_MEAS(1)函数获取测量圆1的内存地址4、设置所要求的应力大小及其精度，并适当CYCLE达到初始应力状态SETreq_isostr=-4e5req_isostr_tol=0.01cyc205、根据理论部分所提供的计算方法，计算为了达到所要求应力大小及精度所需要放缩半径的比例其核心程序如下：_lambda=dim_sum=0.0bp=ball_headloopwhilebp#nullcp=b_clist(bp)loopwhilecp#nullifc_nforce(cp)#0.0then;notavirtualcontact_rcp=sqrt((c_x(cp)-b_x(bp))?2+(c_y(cp)-b_y(bp))?2)if_rcp<=b_rad(bp)then;just-touchingoroverlappingifc_ball1(cp)=bpthenbp_other=c_ball2(cp)elsebp_other=c_ball1(cp)end_ififpointer_type(bp_other)=101then;otherentityisawall_phi=b_rad(bp)else_phi=b_rad(bp)+b_rad(bp_other)end_if_sum=_sum+_rcp*c_kn(cp)*_phiend_ifend_ififc_ball1(cp)=bpthen;getnextcontactof[bp]cp=c_b1clist(cp)elsecp=c_b2clist(cp)end_ifend_loopbp=b_next(bp)end_loop_alpha=-1.0*_lambda*tot_vol*_diso/_sumbp=ball_headloopwhilebp#nullb_rad(bp)=(1.0+_alpha)*b_rad(bp)bp=b_next(bp)end_loop6、适当CYCLE,达到平衡，计算测量圆MP内的平均应力。

若不满足精度，则回到步骤5计算方法见手册理论部分其核心程序如下：loopwhile1#0_diso=req_isostr-isostrifabs(_diso/req_isostr)<=req_isostr_tolthenexitend_ifoo=out('Currentisotropicstress='+string(isostr))expand_radiicommandsolveend_commandend_loop获得非各向同性应力状态使用FISH函数实现，理论复杂，暂时略去。