基于颗粒流微波诱发损伤不规则颗粒数值模拟研究方法

摘要

本发明提供了一种基于颗粒流微波诱发损伤不规则颗粒数值模拟研究方法，其基本算法是建立规定尺寸的热力耦合数值模型，赋予平行粘接性能，利用编写的FISH语言算法生成圆形颗粒，进而形成基于骨料形状的不规则clump块体，采用“试凑法”建立宏观热力参数与细观热力参数的对应关系，为不同骨料赋予不同的热力细观参数，对不同间断比尺下的微波诱发损伤进行数值模拟研究。本专利所述的研究方法对于采用采用离散单元法对微波诱发损伤的不规则颗粒数值模拟提供了一种新的研究思路。

说明书

技术领域

本发明涉及一种基于离散单元法的数值模拟方法，尤其是一种用于微波诱发损伤不规则颗粒的数值模拟研究方法。

背景技术

离散单元法是基于离散元的理论发展起来的，通过颗粒组合宏观的物质并表现出宏观的力学行为，克服了传统连续介质力学模型的宏观连续性假设，适合于模拟岩体大变形、开裂等非线性行为。其基本理论是基于牛顿第二定律与力-位移定律，通过运动定律更新颗粒-颗粒与颗粒-边界的位置；通过力-位移定律更新接触部分的接触力，达到最新的平衡。

在目前研究进程中，采用离散单元法对非连续体材料的模拟多限于力学性能的模拟研究，而针对热力耦合所产生的损伤研究较少，对于微波诱发损伤的不规则颗粒数值模拟研究更为少数。

发明内容

为解决上述问题，本发明提出一种基于离散单元法的微波诱发损伤数值模拟研究方法。

实现本发明目的的一种基于颗粒流微波诱发损伤不规则颗粒数值模拟研究方法，该方法包括以下步骤：

步骤S1、建立规定尺寸的墙体，确定数值模型，设定所模拟的颗粒数量，按照要求的颗粒级配随机分布球状颗粒，随机选中一个颗粒，按照方程式(1)：X＝X_bp+0.0003*k；Y＝Y_bp+0.0003*k生成其周围的其它颗粒；

步骤S2、在中心位置颗粒确定后，按照方程式(2)：(X/a)²ⁿ+(Y/b)²ⁿ＝1的算法生成块体；

步骤S3、当生成的颗粒在椭圆边界内部时，颗粒簇即被认为是生成的clump，按照模拟要求确定生成单个clump所需的颗粒数。在第一个椭圆生成后其根据方程式(3)：X＝X′cosθ-Y′sinθ；Y＝X′sinθ-Y′cosθ；θ＝π*k进行旋转变换；

步骤S4、将方程式(3)代入方程式(2)中即得到算法所需的椭圆旋转方程式(4)：(X′cosθ-Y′sinθ)²/a²+(X′sinθ-Y′cosθ)²/b²＝1，同时结合方程式(5)：可以得到第二个椭圆形，同时根据方程式(5)中长轴和短轴的比例来确定内部颗粒的数量。

步骤S5、根据方程式(6)：S(n)＝S(E_i)∩S(E_i-1)和方程式(7)：椭圆E_i(i≥2)和椭圆E_i-1(i≥2)之间的覆盖面积被定义为S(n)。如果在新的两个椭圆E_i(i≥2)和E_i-1(i≥2)之间有新的颗粒产生，那么颗粒数为N(p)_same；当则没有相同的颗粒，且N(p)_same＝0。

步骤S6、根据以上的算法并结合方程式(8)：构成了新的clump颗粒数CN_newi，按照这种算法依次循环直至满足方程式(9)：此时椭圆E_i+1中的颗粒完全包含在椭圆E_i中。

步骤S7、按照经验赋予数值模型初始细观参数，不同类型骨料赋予不同的细观参数，将模拟的数值结果与物理试验进行对照，使用“试凑法”反复调整参数，确定宏观热力参数与细观热力参数的关系；

步骤S8、将每个颗粒看做一个热存储器，颗粒之间通过热流管来连接，编写基于离散元的热力耦合数值模型代码，依照热传导方程式(10)：

$>\frac{\partial \left(K_{ij}\frac{\partial T}{\partial x_j}\right)}{\partial x_i}+2{\mathrm{\pi fϵ}}_0{{ϵ}_r}^{\prime \prime }{E_0}^2={\mathrm{\rho C}}_p\frac{\partial T}{\partial t}$>

可得到指定功率密度和照射时间的温度数值。数值模型中的颗粒在微波辐射作用下，其颗粒半径和颗粒之间的平行粘接力会发生变化，此变化通过相应的热应变和粘结力来调整。

步骤S9、开展数值模拟试验，运行离散元数值程序，按照编写的不规则颗粒算法进行步数迭代，调用热力耦合方程，探究微波辐射下不同间断比尺的微波损伤，并得出数值模拟结果。

所述的步骤1所述的方程式(1)中，X代表颗粒的中心点横坐标，Y代表颗粒的中心点纵坐标，X_bp代表任意颗粒的中心点横坐标，Y_bp代表任意颗粒的中心点纵坐标，k代表(0，1)范围内的随机数。

所述的步骤2所述的方程式(2)中，a代表椭圆形块体的长轴，b代表椭圆形块体的短轴。当n为1的时候，趋向于形成椭圆形块体；当n大于或等于2的时候，随着n值的增大所生成的颗粒越为规则。

所述的步骤3所述的方程式(3)中，θ代表椭圆的旋转角。

所述的步骤4所述的方程式(4)中，代表颗粒半径的最小值，代表颗粒半径的最大值，k代表范围(0，1)的随机数。

所述的步骤5所述的方程式(6)中，S(E_i)代表椭圆E_i(i≥2)的面积，S(E_i-1)代表椭圆E_i-1(i≥2)的面积。

所述的步骤8所述的方程式(10)中，K_ij代表热传导率张量，ρ代表密度，C_p代表定体积比热容，T代表温度，f代表微波辐射的频率，ε₀代表真空中的介电常数，ε_r″代表介质的耗损因子，E₀代表电场强度。

本专利可以对微波诱发作用下非连续体材料在不同间断比尺照射作用下的微波诱发的损伤进行数值模拟，对裂纹的萌生、扩展、延伸、贯穿等过程进行微观研究，进而对于采用采用离散单元法对微波诱发损伤的不规则颗粒数值模拟提供了一种新的研究思路。

附图说明

图1为本发明的颗粒单元图

图2为本发明的圆的中心范围图

图3为本发明的生成第一个椭圆的示意图

图4为本发明的生成第一个clump的示意图

图5为本发明的生成第二个clump的示意图

图6为本发明的生成最后一个clump的示意图

具体实施方式

如图1至图6所示，本发明的一种基于颗粒流微波诱发损伤不规则颗粒数值模拟研究方法，该方法包括以下步骤：

步骤S1、建立规定尺寸的墙体，确定数值模型，设定所模拟的颗粒数量，按照要求的颗粒级配随机分布球状颗粒，随机选中一个颗粒，按照方程式(1)：X＝X_bp+0.0003*k；Y＝Y_bp+0.0003*k生成其周围的其它颗粒；

步骤S2、在中心位置颗粒确定后，按照方程式(2)：(X/a)²ⁿ+(Y/b)²ⁿ＝1的算法生成块体；

步骤S3、当生成的颗粒在椭圆边界内部时，颗粒簇即被认为是生成的clump，按照模拟要求确定生成单个clump所需的颗粒数。在第一个椭圆生成后其根据方程式(3)：X＝X′cosθ-Y′sinθ；Y＝X′sinθ-Y′cosθ；θ＝π*k进行旋转变换；

步骤S4、将方程式(3)代入方程式(2)中即得到算法所需的椭圆旋转方程式(4)：(X′cosθ-Y′sinθ)²/a²+(X′sinθ-Y′cosθ)²/b²＝1，同时结合方程式(5)：可以得到第二个椭圆形，同时根据方程式(5)中长轴和短轴的比例来确定内部颗粒的数量。

步骤S5、根据方程式(6)：S(n)＝S(E_i)∩S(E_i-1)和方程式(7)：椭圆E_i(i≥2)和椭圆E_i-1(i≥2)之间的覆盖面积被定义为S(n)。如果在新的两个椭圆E_i(i≥2)和E_i-1(i≥2)之间有新的颗粒产生，那么颗粒数为N(p)_same；当则没有相同的颗粒，且N(p)_same＝0。

步骤S6、根据以上的算法并结合方程式(8)：构成了新的clump颗粒数CN_newi，按照这种算法依次循环直至满足方程式(9)：此时椭圆E_i+1中的颗粒完全包含在椭圆E_i中。

步骤S7、按照经验赋予数值模型初始细观参数，不同类型骨料赋予不同的细观参数，将模拟的数值结果与物理试验进行对照，使用“试凑法”反复调整参数，确定宏观热力参数与细观热力参数的关系；

步骤S8、将每个颗粒看做一个热存储器，颗粒之间通过热流管来连接，编写基于离散元的热力耦合数值模型代码，依照热传导方程式(10)：

$>\frac{\partial \left(K_{ij}\frac{\partial T}{\partial x_j}\right)}{\partial x_i}+2{\mathrm{\pi fϵ}}_0{{ϵ}_r}^{\prime \prime }{E_0}^2={\mathrm{\rho C}}_p\frac{\partial T}{\partial t}$>

可得到指定功率密度和照射时间的温度数值。数值模型中的颗粒在微波辐射作用下，其颗粒半径和颗粒之间的平行粘接力会发生变化，此变化通过相应的热应变和粘结力来调整。

步骤S9、开展数值模拟试验，运行离散元数值程序，按照编写的不规则颗粒算法进行步数迭代，调用热力耦合方程，探究微波辐射下不同间断比尺的微波损伤，并得出数值模拟结果。

所述的步骤1所述的方程式(1)中，X代表颗粒的中心点横坐标，Y代表颗粒的中心点纵坐标，X_bp代表任意颗粒的中心点横坐标，Y_bp代表任意颗粒的中心点纵坐标，k代表(0，1)范围内的随机数。

所述的步骤2所述的方程式(2)中，a代表椭圆形块体的长轴，b代表椭圆形块体的短轴。当n为1的时候，趋向于形成椭圆形块体；当n大于或等于2的时候，随着n值的增大所生成的颗粒越为规则。

所述的步骤3所述的方程式(3)中，θ代表椭圆的旋转角。

所述的步骤4所述的方程式(4)中，代表颗粒半径的最小值，代表颗粒半径的最大值，k代表范围(0，1)的随机数。

所述的步骤5所述的方程式(6)中，S(E_i)代表椭圆E_i(i≥2)的面积，S(E_i-1)代表椭圆E_i-1(i≥2)的面积。

所述的步骤8所述的方程式(10)中，K_jj代表热传导率张量，ρ代表密度，C_p代表定体积比热容，T代表温度，f代表微波辐射的频率，ε₀代表真空中的介电常数，ε_r″代表介质的耗损因子，E₀代表电场强度。

上面所述的实施例仅仅是对本发明的优选实施方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神前提下，本领域普通工程技术人员对本发明技术方案做出的各种变形和改进，均应落入本发明的权利要求书确定的保护范围内。