[发明专利]GIS设备内部多金属颗粒流动特性的三维多场耦合仿真方法

权利要求书

1.GIS设备内部多金属颗粒流动特性的三维多场耦合仿真方法，其特征在于：该方法包括以下步骤：

步骤1：以球形的拉格朗日粒子为基础的数学模型单元，设定好构成每一个金属颗粒的拉格朗日粒子中的每一个粒子的密度ρ_ii、半径r_ii、杨氏模量Y_ii、泊松比B_ii，构造出需要的颗粒模型，其中下标ii表示构成第i个颗粒的第i个拉格朗日粒子的序号；对于球形颗粒，ρ_i＝ρ_ii，r_i＝r_ii，Y_i＝Y_ii，B_i＝B_ii，ρ_i是第i个颗粒的密度，r_i是第i个颗粒半径，Y_i是第i个颗粒的杨氏模量，B_i是第i个颗粒的泊松比，对于线型颗粒可以利用刚性键连接拉格朗日粒子的方法生成；

步骤2：设定计算域，计算域的大小取待研究的GIS设备的几何模型的尺寸；

步骤3：设定边界壁的条件，边界的位置被计算域包含在内，GIS设备表面材料根据实际情况设置为固定的壁面，材料根据实际需要自行设置其密度ρ_b、杨氏模量Y_b和泊松比B_b，当金属颗粒在与壁面碰撞时，颗粒会发生反弹，同时会发生电荷转移；

步骤4：设定颗粒的碰撞模型hertz tangential，对于多颗粒系统，补充设定颗粒间的相互作用力库仑力f_q,ij和弹性碰撞力f_c,ij，下标ij表示第i和第j个颗粒之间的相互作用，因此步骤3中的颗粒碰撞力和电荷转移的过程采用如下公式进行计算：

对于颗粒发生碰撞电荷转移，其可能性有两个：

Case1：颗粒间的相互碰撞：

其中，Δq_i表示第i个颗粒表面发生的电荷转移的量，Δq_j表示第j个颗粒表面发生的电荷转移的量，C_i是第i个颗粒的电容，C_j是第j个颗粒的电容，第i个颗粒的电容C_i＝4πεr_i，第j个颗粒的电容C_i＝4πεr_j，ε是背景流体的绝对介电常数，r_i是第i个颗粒的半径，r_j是第j个颗粒的半径；

Case2：颗粒与壁/电极的相互碰撞，颗粒与壁/电极碰撞时，其原有的电荷会中和在壁/电极上，并重新荷载新的电荷，采用下式求得q_i：

其中，n是从壁或电极表面指向外侧的法向量，S_i是第i个颗粒的表面积，E是颗粒与电极碰撞时刻的瞬时电场强度；

颗粒间的弹性碰撞力：

f_c,ij＝(k_nδ_ijn_ij-γ_nv_n,ij)+(k_tt_ij-γ_tv_t,ij)

其中，δ_ij＝|x_i-x_j|-r_i-r_j，为颗粒i和j之间重叠的距离，x_j为颗粒j的位移；n_ij＝(x_i-x_j)/|x_i-x_j|为颗粒i和j碰撞的碰撞法向量；v_n,ij＝[(v_i-v_j)·n_ij]n_ij为颗粒i和j碰撞的相对速度的法向分量；v_t,ij＝v_i-v_j-v_n,ij-(ω_i+ω_j)(x_i-x_j)/2为颗粒i和j碰撞的相对速度的切向分量，ω_i为颗粒i和j发生碰撞之后的颗粒i的旋转角速度，ω_j为颗粒i和j发生碰撞之后的颗粒j的旋转角速度；为切向位移，t_i是冲击时间，t是碰撞接触的总时间，Δt是DEM的时间步长；对于颗粒i碰撞壁面的情况，则δ_ij＝x_i-r_i；

最后，采用下式计算出hertz tangential模型中的非线性弹性常数k_n、k_t、γ_n、γ_t：

其中，Y_eff为颗粒i和j发生碰撞时的等效的杨氏模量，G_eff为颗粒i和j发生碰撞时的等效剪切模量，r_eff为颗粒i和j发生碰撞时的颗粒半径，m_eff为颗粒i和j发生碰撞时的颗粒质量，β_ij为颗粒i和j发生碰撞时的阻尼比；

采用下式计算出Y_eff、G_eff、r_eff、m_eff和β_ij：

其中，Y_i为第i个颗粒的杨氏模量，B_i为第i个颗粒的泊松比，Y_j为第j个颗粒的杨氏模量，B_j为第j个颗粒的泊松比，m_i为第i个颗粒的质量，m_j为第j个颗粒的质量，e_ij为碰撞恢复系数；

步骤5：对每一个颗粒进行受力分析，并根据牛顿第二定律计算每一个颗粒的速度和位置，采用以下公式计算出每一个金属颗粒的转动惯量：

其中，是第i个颗粒的速度对时间的导数，g是重力加速度，V_i是第i个颗粒的体积，ρ_f是背景流体的密度，P是气体的压强，τ是流体粘性应力张量，f_dep,i是第i个颗粒的电场梯度力，f_drag,i是第i个颗粒的拖曳力，f_vm,i是第i个颗粒的虚拟质量力，f_mag,i是第i个颗粒的Magnus力，I_i是第i个颗粒的转动惯量矩阵，是第i个颗粒角速度的一阶导数；

上式中的参数采用以下公式求得：

其中，C_d为拖曳力系数，U为流体的速度，为流体的速度对时间的一阶导数，C_l为经验系数，C_d的取值采用下式求得：

lgC_d＝1.51exp[-0.2865(lgR_e,i-1.56)]-0.2exp[-0.1(lgR_e,i-0.4127)²](lgR_e,i-0.4127)-1.04

其中，R_e,i是第i个颗粒的雷诺数，R_e,i＝2ρ_fr_i|U-v_i|/μ，μ是流体的粘性系数；

步骤6：采用下式计算背景流场的速度和压强：

其中，α是体积分数，V_cell是FVM单元的体积，t是时间；

步骤7：根据步骤6求得的流场的速度U计算步骤5中每一个颗粒的f_drag,i、f_mag,i和f_vm,i，对于初始时刻，即t＝0时刻来说，按U＝0来计算上述力；

步骤8：采用下式计算由颗粒上的电荷所产生的泊松电场E_p：

E_p＝∑E_pm

其中，E_pm是第m个FVM单元体心上的由金属颗粒带电而产生的泊松电场，E_p是整个GIS被划分成的所有的FVM单元上的E_pm的总和，为第m个FVM单元体心上的由金属颗粒带电而产生的电位，且Q_pm是第m个FVM单元内包含的所有的金属颗粒所带的电荷量的总和：

其中，V_cell-m是第m个FVM单元的体积，q_pm是第m个FVM单元内包含的颗粒所带电荷的总量，α_m则是第m个FVM单元的体积分数；最后将E_p与背景电场E_b相叠加，得到E：

U_c＝U_msin(2πf_act+θ)

其中，U_c是交流电压，U_m是U_c的幅值，f_ac是频率，当计算工频交流电工况时，f_ac取50Hz，当计算直流电工况时，f_ac取0，θ是初相位，取0，r是颗粒到电极的径向距离，R₂是GIS模型的外径，R₁是GIS模型的内径，采用公式E＝E_b+E_p计算出E；

步骤9：按照上述步骤编写基于Linux系统下的LIGGGHTS代码和OpenFOAM代码及其算例脚本文件，编写编译完成之后按照CFDEM耦合的流程进行并行计算；

步骤10：利用Paraview软件将计算结果后处理生成颗粒流动的实时动态分布的图像。