基于随机生长四参数法的多孔介质模型.pdf

能源研究与管理2023,15(2)能源研究与管理2023,15(2)收稿日期：2023-02-20基金项目：核技术应用教育部工程研究中心开放基金（HJSYB2022-1）第一作者：王召（1997），女，硕士研究生，主要从事物理与数值模拟研究。E-mail：*通信作者：杜俊杰（1987），男，讲师，博士，主要从事超导磁体及常规磁体技术研究。E-mail：引用格式：王召,杜俊杰,刘俊雄.基于随机生长四参数法的多孔介质模型J.能源研究与管理,2023,15(2):88-94.摘要：为了研究所列因素的改变对生成孔隙空间以及程序运行时间的影响，使用向量乘法建立四参数随机生长方法（QSGS）模型，并讨论模型参数对生成孔隙结构及计算机生成孔隙模型耗时的影响，使用Matlab软件针对不同孔隙度、固相核分布和生长概率参数建立二维和三维孔隙模型。结果表明：孔隙率与固相核占比呈正相关，孔隙率选取过大会生成无效孔隙，过小会导致固相核孤立。在二维和三维模型中，模孔隙率参数越大生成孔隙模型时间越长，而生长概率与孔隙分布参数改变对模型生成时间影响较小。在三维模型情况下由于消耗更多计算资源，因此研究计算区域网格数对孔隙模型生成时间影响，且随网格数增加模型生成时间显著增加。关键词：随机生长四参数法；多孔介质；土壤；孔隙模型中图分类号：X502文献标志码：A文章编号：20967705（2023）02008807WANG Zhao1,2,DU Junjie1,*,LIU Junxiong2（1.Engineering Research Center of Ministry of Education for Nuclear Technology Application,East China Institute of Technology,Nanchang 330013,China;2.School of Nuclear Science andEngineering,East China Institute of Technology,Nanchang 330013,China）To investigate the impact of changes in the listed factors on the generation of pore space and program runtime.Multiplication of vectors is used to the establish of four-parameter stochastic growth method(QSGS)model.The influence ofmodel parameters on the generation of pore structure and the time consumption of computer generated pore model isdiscussed.Using Matlab software,two-dimensional and three-dimensional pore models are established for different porosity,solid core distribution and growth probability parameters.Numerical simulation and analysis are carried out for differentporosity,distribution probability of solid phase nuclei and growth probability,and the changes of program running timecaused by parameter changes are counted.In two dimensions,the porosity is positively correlated with the proportion of thesolid phase.High porosity generates invalid pores and low porosity leads to the isolation of solid phase nuclei.The smaller isthe distribution probability,the fewer is the number of initial solid phase nuclei.The growth probability affects the poredistribution in precision.In three dimensions,the effect of different parameters on the pore distribution is similar to that of 2D,but the running time is significantly increased.In addition,it is found that in three dimensions the running time increasesDOI：10.16056/j.2096-7705.2023.02.013基于随机生长四参数法的多孔介质模型王召1,2，杜俊杰1,*，刘俊雄2（1.东华理工大学核技术应用教育部工程研究中心，南昌330013；2.东华理工大学核科学与工程学院，南昌330013）研究与探讨88窑窑能源研究与管理2023,15(2)能源研究与管理2023,15(2)significantly when the number of grids increases.quartet structure generation set(QSGS);porous medium;soil;heavy metal speciation引言土壤孔隙指的是土壤颗粒之间或聚团物之间存在的间隙。土壤的孔隙作为容纳水、空气和溶质的空间，对于提高作物产量以及治理环境污染等问题起到非常重要的作用。例如在土壤孔隙中含水、空气和溶质容量和种类不同会直接影响到土壤里的动物以及微生物的数量以及质量，对农作物产量产生不同程度影响。同时在我国重金属污染土壤的环境问题研究中1-4，张荣辉等5通过高氧官能团的多孔碳材料为铁氧化物结合土壤中的重金属离子提供了有效的表面配体，可以协同铁氧化物修复重金属污染土。因此孔隙结构对于相关问题的研究具有重要的意义。但在理论计算孔隙流，水迁移和污染物扩散等问题中，由于孔隙结构复杂往往不能使用传统方法构造孔隙通道，并且随着计算机技术不断更新发展，因此使用计算机生成孔隙结构成为具有成本更低、效率更高且更直观方法之一。传统的测量土壤孔隙率的方法有直接测出密度和比重来计算及压汞法6，但这2种方法对于研究孔隙的具体形态效果较差。随着图像处理技术的快速发展，土壤切片技术7、X射线断层扫描8-10以及核磁共振6等技术能更为直观地观测孔隙，但也有技术复杂或有成本高等缺陷。简化孔隙结构建立模型，通过计算机进行模拟分析，既能解决成本问题，又能解决其他方法样本量小代表性不强等问题是很好的选择。选择更优的建模方法，建立更准确实用的模型成为了新课题。首先，由WANG等11提出了可人为控制的多孔介质微观结构重构方法，即四参数随机生长方法（quartet structure generation set，QSGS）。与其他随机生长算法相比，QSGS融入更多物理意义明确的参数，使得多相之间的接触连接问题得到了较好解决，同时颗粒分布更加随机，统计特性能够得到更效提取。然而该模型中颗粒尺寸单一，与自然多孔介质中颗粒尺寸多变特性相比仍存在一定差别，不能够很好地表征多孔介质的尺度不变属性。张季如等12对QSGS方法的微观孔隙结构的分析特征进行了研究，为参数的合理化取值提供参考依据。罗齐彬等13及张朔等14利用QSGS方法建立二维微观孔隙模型，分别使用有限元和格子玻尔兹曼方法（lattice Boltzmann method，LBM）对孔隙内流体进行模拟，得出了空隙度、缩放覆盖率等变化对孔隙分布的影响结果。李滔等15使用QSGS方法建立的三维孔隙模型利用LBM模型对三维孔隙结构模拟并对孔隙流体进行研究。基于前人的研究以及深入的讨论，本文利用Matlab来实现QSGS方法，选取孔隙度、固相核分布概率s、生长概率g与网格数作为变量，建立二维和三维孔隙模型，进一步利用LBM模型对孔隙流体的数值模拟。1模型建立1.1二维随机生长模拟为得到最终的多孔介质模型，QSGS算法将土壤基质分为固相以及非固相（孔隙）。在固定区域内首先随机构建固相生长内核，再在选择方向上让此内核按照给定生长概率连续随机生长。土的团聚体由细小的矿物颗粒经凝聚和胶结作用而形成3,16-17，是一个渐进的过程，而QSGS算法的固相生长法则类似于土的团聚体形成过程18-19。以土粒和孔隙二相构成的二维土体孔隙结构为例，若令土粒为生长相，孔隙为非生长相，构造区域的初始相全为孔隙，则构造过程如下：1）在构造区内的每个网格节点上，以分布概率d随机分布固相生长核，在最终构造出的孔隙结构中，d不能大于土粒面积的百分比，即d臆1；2）网格节点上的各生长核按照给定的方向生长概率（代表生长方向，=1,2,8），二维生长方向示意图如图1所示，从图中所示的四面八方向周围相邻网格节点生长。若8个方向的生长概率18相同，则生成各向同性的孔隙结构；3）重复步骤2）直至孔隙面积达到给定的孔隙度。应该指出的是，QSGS算法通过、d、(,)4个参数来影响和控制多相多孔介质微观结构的生成，其中(,)为概率密度，表示在方向上第相在第相上的生长概率，反映了各相之间的相互影响。对于土粒和孔隙二相构成的土体，只需采用、d和3个参数，即可获得表征土体孔隙结构的QSGS模型20-22。研究与探讨89窑窑能源研究与管理2023,15(2)能源研究与管理2023,15(2)由图1可知，在二维情况下只有8个生长方向，从正上方开始沿顺时针方向的标号分别为18，为了计算方便，将生成随机数小于固体核生成概率的点定义为(,)，沿18方向对应的位置坐标行向量为0 1 1 1 0111、列向量1 1 01110 1。由此将每次生存随机数小于生长方向概率向量坐标定位1，其他为0，这样通过简单的对应相乘便可得到生长方向的坐标。可以发现整个过程类似于随机行走，模型在这种移动模型中，移动节点随机选择一个方向和速度来从当前位置移动到新的位置。而随机行走模型也被称为布朗运动。使用以上方法另外一个难点是在固体核生长到边界时由于边界外未定义的原因使得代码不能继续进行，因此根据4边、4顶点的限制条件设定了生长边界，并且假设生长到边界时不允许“逆”生长。1.2三维随机生长模型随着维度的增加生长方向成指数式增长，正方体的共有26个生长方向如图2所示，左图为其中14个记为1，右图为其余12个记为2。其中设定了6面、8顶点和12边的限制生长条件。通过点、线、面及逆时针和由上至下的命名规则对26个生长方分别向进行标号其中18为顶点、920为边、2126为面，其对应向量分别为行向量()1 1111 1 1111 1 1 1 01 0 0 1 011 0 1 0 00、列向量()1 1 1111 1 1 111 1 0 1011 0 1 0 01 0 1 0 0、纵向量()1 1 1 11111 0 0 0 0 1 1 1 11111 0 0 0 0 11。同样假设在固相核生长到边界时不允许“逆”生长。图 1二维生长方向示意图Fig.12D growth direction diagram初始生成的随机数分布，通过与分布概率的大小比较来确定固相核的位置，然后对应每个固相核生成26个生长方向的向量，则构成了初始位置矩阵，同样与三维向量、简单的对应相乘随后便可得到随机生长方向的位置坐标。整个生长过程较为简单，依然与布朗运动情况类似，而在处理边界固相核生长时较为复杂。循环以上过程直到孔隙率达到要求为止，便可构建所需参数下的孔隙分布模型。2模拟结果与讨论2.1二维随机生长法模拟结果选用Matlab采用本文1.1部分的具体方法建立循环，成立模型。在这一模型下影响最终成像结果的变量为以下所列，通过控制变量法分别观测不同变量的变化对模拟结果的影响，以及导致的计算时间变化。孔隙度由孔隙体积占基质总体积的百分比来表示。指基质中通气孔隙与持水孔隙的总和。孔隙度可分为有效孔隙度和无效孔隙度。有效孔隙度是指那些互相连通的，在一般压力条件下，可允许流体在其中流动的孔隙体积之和与基质总体积的比值，与之相对的则是无效孔隙度。在考虑土壤实际应用价值以及污染治理时，被考虑较多的为有效孔隙度。生长概率即固相内核在选定方向生长的概率，在实际中可以选择各向同性也可选择各向异性。固相核分布概率会影响所选区域内固相核的分布。此外选取的网格数也影响模拟结果的精度以及孔隙的图 2三维生长方向示意图Fig.23D growth direction diagram（a）1=14（b）2=12研究与探讨90窑窑能源研究与管理2023,15(2)能源研究与管理2023,15(2)图 3s=0.002、g=0.001 时不同孔隙概率的孔隙分布Fig.3Pore distribution of different pore probabilities whens=0.002,g=0.001（a）=0.25（c）=0.45（b）=0.35（d）=0.55随着孔隙率的增加可能会导致孔隙之间并不连通的情况，孔隙变为无效孔隙。而当孔隙率减少，会导致固相核的孤立，不符合实际情况。因此应当选取适当的孔隙率。不同孔隙概率的孔隙分布的计算时间如图4所示，可知随选取的孔隙概率的增加，程序计算所需时间稳步增加。在相同孔隙及分布概率条件下，研究不同生长概率情况下孔隙分布情况，结果如图5所示，由图可知，各项同性下的生长概率并不会对结果产生很大的影响，但是如果选取较小的生长概率会使得计算时间显著增加。计算时间如图6所示，同样发现如果取值过大会影响孔隙率的精度。因此在各向同性且符合实际的情况下尽量选取适中的生长概率。在相同孔隙率及生长概率条件下，研究了不同固相核分布概率情况下孔隙分布情况，结果如图7图 4不同孔隙概率的孔隙分布计算时间Fig.4Calculation time of pore distribution withdifferent pore probability1020.60.40.20.8（a）g=0.001（b）g=0.000 5图 5=0.4、s=0.002 时不同生长概率的孔隙分布Fig.5Pore distribution of different growth probabilities when=0.4,s=0.002分布。对不同的孔隙度、固相核分布概率s、生长概率g进行了数值模拟。首先对孔隙度取值0.2，以0.5%为间隔至0.8，选用1 024伊1 024（伊）网格，分布概率s=0.002、生长概率g=0.001（为了计算方便采用各向同性模型），进行孔隙分布模拟，结果如图3所示。由图3可知，孔隙率与固相占比呈正相关，但研究与探讨91窑窑能源研究与管理2023,15(2)能源研究与管理2023,15(2)所示。由图可知，当分布概率越小时初始的固相核分布个数便越少，当小分布概率时就会出现如图7(a)所示的大面积孤岛情况，而选取的分布概率较大时也会使得固相核孤立分布概率增加。因此在实际过程中通常根据实际情况选取分布概率的大小。单分布概率引起的程序运行时间变化如图8所示，发现选取更大的分布概率大体上会让计算时间缩短。2.2三维随机生长法模拟结果三维情况下，选取孔隙率作为变量的结果如图图 6不同生长概率的孔隙分布计算时间Fig.6Calculation time of pore distribution ofdifferent growth probability2597510g/(103)8635459所示，三维空间（伊伊）中与二维空间进行分析的结果类似，孔隙率与固相占比仍呈正相关，孔隙率的增加可能导致孔隙不连通，孔隙率减少会导致固相核的孤立。计算时间如图10所示，由于增加1个维度，随着孔隙率的增加，计算时长更为明显的随之增加。三维情况下选取生长概率作为变量结果如图11所示，计算时间如图12所示，在三维维度下，生长概率依然对孔隙分布情况影响不大，计算时间相较二维时明显增多的同时，生长概率与计算时间的反比关系更为明显。选取分布概率作为变量结果如图13所示，计算图 7=0.4、g=0.005 时不同分布概率的孔隙分布Fig.7Pore distribution of different distribution probabilities when=0.4,g=0.005（a）s=0.000 2（c）s=0.001（b）s=0.000 4（d）s=0.002601.4g/(103)802001201600.61.01.80.2图 8不同分布概率的孔隙分布计算时间Fig.8Calculation time of pore distribution with differentdistribution probability图 9s=0.003、g=0.005 时不同孔隙概率的孔隙分布Fig.9Pore distribution of different pore probabilities whens=0.003,g=0.005（b）=0.4（a）=0.3（d）=0.6（c）=0.5研究与探讨92窑窑能源研究与管理2023,15(2)能源研究与管理2023,15(2)时间如图14所示，依然与二维维度情况类似，分布概率过大依然会造成如图13(a)所示的大面积孤岛情况，分布概率过小依然会使得固相核孤立。且计算时间依然对比低维度计算时大大增加了。同时，本文还尝试了分别计算64伊64伊64、128伊128伊128、256伊256伊256以及512伊512伊512网格下的孔隙模型，结果如图15所示，由图可以看到随网格数的增加观测到的孔隙越来越多。并且计算时长如如图16所示递增。因此在三维条件下选择合适的网格数与生长概率及分布概率对于节省生成孔隙模型时间以及提升精度具有很大帮助。3结论1）在二维情况下，孔隙率的大小会影响固相核的面积，分布概率的大小对一开始固相核分布个数图 15=0.4,s=0.003,g=0.005 时选取不同网格数生成的孔隙分布Fig.15When=0.4,s=0.003,g=0.005,the poredistribution generated by different mesh numbers is selected（b）128伊128伊128（d）512伊512伊512（a）64伊64伊64（c）256伊256伊256（b）g=0.008 0（a）g=0.004 0（d）g=0.001 4（c）g=0.010 0图 11=0.45,s=0.003 时不同生长概率的孔隙分布Fig.11 Pore distribution of different growth probabilities when=0.45,s=0.003图 12计算不同生长概率的孔隙分布所需时间Fig.12Calculate the time required for pore distribution withdifferent growth probabilities228414g/(103)6263012101418图 10计算不同孔隙概率的孔隙分布所需时间Fig.10Time required to calculate the pore distribution ofdifferent pore probabilities0.91020.70.50.3（b）s=0.001 0（a）s=0.004 0（d）g=0.002 0（c）s=0.001 4图 13=0.4,g=0.005 时不同分布概率的孔隙分布Fig.13Pore distribution of different distribution probabilitieswhen=0.4,g=0.005图 14计算不同分布概率的孔隙分布所需时间Fig.14Time required to calculate the pore distribution withdifferent distribution probabilities262.0g/(103)28301.61.20.80.4研究与探讨93窑窑能源研究与管理2023,15(2)能源研究与管理2023,15(2)有较大的影响，生长概率则是会对迭代过程中计算孔隙率步骤产生较大影响，而本身对于孔隙结构影响较小，过大则会影响生成孔隙的精度。而经过实际计算，三维情况与二维情况的这3种变量对结果的影响十分类似。2）在程序计算时间方面，无论在二维还是三维，大体上孔隙率的增加会使时间增加，生长概率与分布概率的增加会使计算时间减少。随着维度的增加所需计算时间显著增加。选取更多的网格数也会使计算时间增加，尤其是对于三维情况。3）为符合实际土壤孔隙分布谨慎的选取孔隙率、分布概率生长概率等参数是必要的。为控制时间成本，还可以选择在二维层面建模，并选取适当的网格数等。本研究仅选用了所列的几个因素的有限数据集进行了分析，然而，选取相关因素或进一步扩大数据集的规模和范围可能会对结果的稳定性和普遍性产生更多的证据。未来的研究可以考虑更多的相关因素以及扩大数据集。在方法选择方面，可以结合各种新的技术，例如，使用机器学习算法来分析数据，或者采用深度学习模型来预测结果等。这样可以进一步提高研究的准确性和可靠性。总之，本研究只是对土壤孔隙模型的初步探索，还有很多方面可以进一步研究和发展。希望未来的研究可以在这一领域做出更深入的贡献，从而推动相关领域的发展和进步。参考文献1路策,丁园.多孔碳负载铁氧化物的制备及其在重金属污染土壤修复中的应用J.能源研究与管理,2021(1):83-88.2吕贵芬,杨涛,陈院华,等.江西省土壤重金属污染治理研究进展J.能源研究与管理,2016(2):16-18.3刘亚龙,王萍,汪景宽.土壤团聚体的形成和稳定机制:研究进展与展望J/OL.(2022-08-01).http:/ B,张桃林.利用土壤切片的数字图像定量评价土壤孔隙变异度和复杂度J.土壤学报,2003(5):678-682.8杨善莲.基于CT技术的菜地土壤孔隙性质与氮磷养分含量三维分布特征研究D.合肥:安徽农业大学,2020.9张苗.基于CT扫描法探究污染场地土壤大孔隙对Cd迁移的影响D.北京:北京建筑大学,2022.10张俊.不同生育期香蕉地土壤孔隙特征及溶质运移影响研究D.桂林:桂林理工大学,2022.11 WANG Ziyan,JIN Xu,WANG Xiaoqi,et al.Pore-scale geome-try effects on gas permeability in shaleJ.Journal of NaturalGas Science and Engineering,2016,34:948-957.12张季如,钟思维.四参数随机生成法重构土体微观孔隙结构的分形特征J.水利学报,2018(7):814-822.13罗齐彬.隐伏铀矿勘探总的地气测量机理及其应用研究D.南昌:东华理工大学,2019.14张朔,刘顺喜,金毅.多相随机自相似多孔介质表征及其孔-渗关系数值分析J,中国矿业,2019(28):331-335.15李滔,李闽,荆雪琪,等.孔隙尺度各向异性与孔隙分布非均质性对多孔介质渗透率的影响机理J.石油勘探与开发,2019(3):569-579.16 SIX J,ELLIOTT E T,PAUSTIAN K.Soil macroaggregate tur-nover and microaggregate formation:A mechanism for Csequestration under no-tillage agricultureJ.Soil Biology andBiochemistry,2000,32(14):2099-2103.17 EDWARDS A P,BREMNER J M.Use of sonic vibration forseparation of soil particlesJ.Canadian Journal of Soil Sci-ence,1964,44(3):366.18 OLESCHKO K,FIGUEROA S B,MIRANDA M E,et al.Massfractal dimensions and some selected physical properties ofcontrasting soils and sediments of MexicoJ.Soil&TillageResearch,2000,55(1):43-61.19郭飞,徐绍辉,刘建立.土壤图像孔隙轮廓线分形特征及其应用J.农业工程学报,2005(7):6-10.20周潇,申林方,阮永芬,等.基于四参数随机生长法重构土体的渗流细观数值模拟J.排灌机械工程学报,2015,33(4):316-321.21申林方,王志良,李邵军.基于土体细观结构重构技术的渗流场数值模拟J.岩土力学,2015,36(11):3307-3314.22李仁民,刘松玉,方磊,等.采用随机生长四参数生成法构造黏土微观结构J.浙江大学学报(工学版),2010,44(10):1897-1901.图 16采用不同网格数计算时间变化Fig.16Uses different grid numbers to calculate time changes101500400300200100102103研究与探讨94窑窑