SPHysics浮体模型及对应边界处理

SPH方法的特点是采用粒子求和来近似场变量。有效的粒子求和近似依赖空间划分、邻近粒子搜索等算法，如前文介绍。本文主要学习SPHysics中的浮体处理方法。

 

首先需要了解SPH与结构的耦合（刚体、弹性体）都是通过界面粒子实现的。界面粒子当然也属于一类边界粒子，文献中称为“free moving boundary particles”(FM)。相对于固定在壁面上为防止流体粒子四处乱跑而设计的Monaghan repulsive边界粒子(固定边界粒子)，FMBP具有双重作用，既约束流体粒子“四处乱跑”又能将流体对结构的力传递到结构求解器（刚体求解器或有限元求解器）。故程序中作为两类不同的粒子，通过内存位置区别，各自需要的计算自然就不同，详细如下：

 

Interaction pairs

 

1 流体粒子-固定边界粒子  将边界粒子对流体粒子的斥力加入流体粒子加速度，更新流体粒子；固定边界粒子不做更新；

 

2 流体粒子-自由边界粒子  同上完成流体粒子更新后，将边界对流体的反力加入(自由)边界粒子的加速度，(之后由外部结构求解器调用)

 

3 自由边界粒子-自由边界粒子 两粒子不属于同一浮体的话，类似2分别更新这两个自由边界粒子

 

4 自由边界粒子-固定边界粒子 类似1， 只更新自由边界粒子

 

!"碰撞对“为流体粒子-固定边界粒子

if(i.gt.nb.and.j.le.nb) thenax(i)=ax(i)+fxbMON if(j.gt.nbfm) then ! "碰撞对"是 流体粒子-自由边界粒子ax(j)=ax(j)-fxbMONendifelse if(i.lt.nbp1.and.j.lt.nbp1) then !bp-bpif(i.gt.nbfm) then !i is free moving BPif(j.gt.nbfm) then ! "碰撞对”是两自由边界粒子

id_num_FB_i=i_Fb_pointer_info(i) id_num_FB_j=j_fb_pointer_info(j)if(id_num_fb_i.ne.id_num_fb_j) thenax(i)=ax(i)+fxbMON*pm(j)ax(j)=ax(j)-fxbMON*pm(j) endifelse  ! "碰撞对"是 自由边界粒子-固定边界粒子ax(i)=ax(i)+fxbMON*pm(j)endif  ! if(j.gt.nbfm)

 

Monaghan Repulsive

上述可见应该设计两大类的粒子碰撞：流体粒子-边界粒子  及 边界粒子-边界粒子 的碰撞。对于前者，其斥力模型SPH相关文献介绍很多；对于后者，估计是刚体-刚体接触力学，考虑切向摩擦、正压力等等。

 

刚体运动

do i_num_fb=1, num_fb        omegaYdot(i_num_fb)=0.0bigUdot_init = x_friction(i_num_fb)  bigWdot_init = Z_friction(i_num_fb)! big_前缀的量，表示刚体质心的变量（欧拉运动定律）do i=i_start, nb_fb(i_num_fb) bigUdot(i_num_fb)=bigUdot(i_num_fb)+ax(i) ! 质心加速度=sum(界粒子的加速度)bigWdot(i_num_fb)=bigWdot(i_num_fb)+az(i) omegaYdot(i_num_fb)=OmegaYdot(i_num_fb)+ &(zp_minus_R0(i)*ax(i)-xp_minus_R0(i)*az(i)) !角加速度enddo!体力(重力加速度)bigUdot(i_num_fb)=bigUdot(i_num_fb)+bigMass(i_num_fb)*grxbigWdot(i_num_fb)=bigWdot(i_num_fb)+bigMass(i_num_fb)*grz! forces per unit massbigUdot(i_num_fb)=bigUdot(i_num_fb)/bigMass(i_num_fb)bigWdot(i_num_fb)=bigWdot(i_num_fb)/bigMass(i_num_fb)! OmegaYdot!质心位置更新box_xc(i_num_Fb)=box_xc_old(i_num_Fb)+dt2*bigU(i_num_Fb)box_zc(i_num_Fb)=box_zc_old(i_num_Fb)+dt2*bigW(i_num_Fb)! 质心速度更新bigU(i_num_fb)=bigU_old(i_num_fb)+dt2*bigUdot(i_num_Fb)bigW(i_num_fb)=bigW_old(i_num_fb)+dt2*bigWdot(i_num_Fb)do i=i_start, nb_fb(i_num_fb)xp(i)=xp(i)+dt2*up(i) !界面粒子位置更新zp(i)=zp(i)+dt2*wp(i)xp_minus_RO(i)=xp(i)-box_xc(i_num_fb) !转轴距离zp_minus_RO(i)=zp(i)-box_zc(i_num_fb)up(i)=bigU(i_num_fb)+bigOmegaY(i_num_FB)*zp_minus_R0(i) !粒子速度更新wp(i)=bigw(i_num_fb)-bigOmegaY(i_num_FB)*xp_minus_R0(i)enddoenddo 

最近翻了SPHysics 2013.4.24