基于能量守恒耗散粒子动力学方法的自然对流模拟改进研究.pdf

基于能量守恒耗散粒子动力学方法的自然对流模拟改进研究*鲁维1)陈硕1)于致远1)赵嘉毅2)张凯旋3)1)(同济大学航空航天与力学学院力学系,上海200092)2)(上海理工大学能源与动力工程学院,上海200093)3)(南开大学医学院,天津300071)(2023年 3月 31 日收到;2023年 5月 16 日收到修改稿)能量守恒耗散粒子动力学(eDPD)是一种研究热输运过程的介观尺度数值模拟方法,然而在 eDPD 系统内引入 Boussinesq 假设以研究自然对流问题时,eDPD 系统自身的热膨胀性对模拟结果的影响常常被忽略.首先研究了 eDPD 系统的热膨胀性,通过模拟获得 eDPD 系统的热膨胀系数 b;并由此模拟了不同瑞利数Ra、不同几何结构下的自然对流;利用 eDPD 系统自身的热膨胀性,在不引入 Boussinesq 假设的前提下获得了合理的温度场和速度场,与相同 Ra 数下有限体积法模拟结果相比,误差明显小于以往研究中相同条件下的对比误差.研究表明在 eDPD 系统中引入 Boussinesq 假设时,需要考虑 eDPD 系统自身热膨胀性的影响,并且进一步对 Ra 数的计算进行了修正.关键词：能量守恒耗散粒子动力学,Boussinesq 假设,自然对流,热膨胀性PACS：02.70.Ns,47.11.j,02.70.cDOI:10.7498/aps.72.202304951引言耗散粒子动力学(dissipativeparticledynam-ics,DPD)是一种介观尺度流体动力学数值模拟方法,由 Hoogerbrugge 和 Koelman1于 1992 年提出,随后 Espaol 和 Warren2构筑了该方法的理论框架.DPD 被称为“弥补宏观尺度和微观尺度之间空白的方法”3.DPD 系统中以离散粒子作为动量载体,粒子在离散的时间和空间上运动,这些粒子的运动状态等价于一个指定范围内大量分子的集中行为体现,也即 DPD 研究介观尺度上粗粒化粒子的运动状态.传统的 DPD 模拟的是等温系统,不能模拟介观尺度上粒子之间热输运问题,Avalos 和 Mackie4与 Espaol5在 DPD 的基础上提出了能量守恒耗散粒子动力学(eDPD),可描述介观尺度上的粒子间的能量传递.Ripoll 等6,7通过与一维、二维传热问题的理论对比验证了 eDPD 模型对于模拟热输运过程的可行性;Mackie 等8,9构建了 eDPD 参数与热传导系数以及热扩散率之间的关系;Lukes10通过量纲分析对 eDPD 进行了模拟参数与物理参数之间的映射研究.近年来 Homman 等11和Stoltz12分别研究了 eDPD 并行迭代算法以及算法稳定性,以谋求更大的时间步长.随着理论基础的不断完善,eDPD 被广泛应用到各个领域.Abu-Nada 完成了一系列基于 eDPD 的研究,从固体导热问题13,14到 Rayleigh-Bnard(RB)以及封闭空腔中的自然对流问题等15,16,显示 eDPD 能模拟复杂自然、受迫热对流等过程17,18.但是上述研究中当瑞利数 Ra 较小时,Abu-Nada 使用 eDPD 的模拟结果不理想,模拟方腔内自然对流时,相同条件*国家自然科学基金(批准号:11872283,12002212)和上海扬帆计划(批准号:20YF1432800)资助的课题.通信作者.E-mail:schen_2023中国物理学会ChinesePhysicalSocietyhttp:/物理学报ActaPhys.Sin.Vol.72,No.18(2023)180203180203-1下的 eDPD 模拟结果与有限体积法(finitevolumemethod,FVM)模拟结果等温线对比存在较大误差.我们分析,误差的主要原因是以往研究中没有考虑 eDPD 本身热膨胀性,这部分热膨胀性在适当的条件下也会产生相应的自然对流,如果忽略其带来的影响则会在 Ra 数量级较小时带来较大的误差.在 Abu-Nada 等的研究中,eDPD 系统被视为不可压缩系统,直接引入 Boussinesq 假设来提供模拟中所需的浮升力.然而 Mai-Duy 等19和 Pan等20的研究显示由于 DPD 采用粒子间“软”相互作用,DPD 系统本身具有一定的可压缩性,在系统温度较高或保守力足够大的情况下才可认为系统近似不可压缩.因此通常情况下 eDPD 研究自然对流问题时应当考虑系统自身热膨胀性的作用,否则会使 Ra 数的计算结果与实际的 Ra 数出现较大偏差,体现在等温线对比图上就会使相同温度的等温线无法精确拟合.本文首先通过模拟获得 eDPD 系统的热膨胀系数 b,随后在不引入 Boussinesq 假设的前提下,利用 eDPD 系统自身的热膨胀性,研究了不同 Ra数、不同几何结构下的自然对流问题,获得了合理的温度场和速度场.在相同条件下与 FVM 的模拟结果进行对比,误差要小于 Abu-Nada 等1318的模拟误差.依靠 eDPD 系统自身的热膨胀性可以形成自然对流,表明在引入 Boussinesq 假设时,需要考虑 eDPD 系统本身热膨胀性的影响.本文以 Ra为 15000 时同心圆环内自然对流的模拟为例,分析 eDPD 系统热膨胀性对不可压缩假设下采用 Bo-ussinesq 假设进行自然对流模拟所产生的影响.进一步,本文对 Ra 数的计算进行修正,结合实际模拟的 Ra 数大小,对系统自身引入 Boussinesq 假设的条件进行评估.2能量守恒耗散粒子动力学模型rc与分子动力学类似,eDPD 方法是一种在一定的截断半径 内基于粒子与相邻粒子之间相互作用的粒子方法,粒子的运动符合牛顿第二定律:dridt=vi,dvidt=Fi,(1)dvidt=(fCij+fDij+fRij).(2)viFiri这里,表示粒子速度,表示粒子受到的力,fCij,fDijfRij表示粒子的位置,和 分别为保守力、耗散力和随机力,分别表示为fCij=j=iaijc(rij)eij,(3)fDij=j=iijD(rij)(eij vij)eij,(4)fRij=j=iijR(rij)ijt12eij,(5)aijaij=A kB(Ti+Tj)/2kBTiTjijijjiijijeijvijtDCR式中,为粒子之间的保守力系数与温度有关,这里我们取 ,A=18.75,是玻尔兹曼常数,和 是粒子温度;是高斯白噪音,同时具有对称性,即 =,以此保证整体的动量守恒;和 是耗散力系数和随机力系数;是粒子 i 和粒子 j 的单位向量;是粒子 i 和粒子 j 之间的相对速度;是时间步长;耗散力 权函数与保守力权函数 和随机力权函数 之间的关系分别为D(rij)=C(rij)=(1 rijrc),rij rc,0,rij rc,(6)D(rij)=R(rij)2=(1 rijrc)s,rij rc,0,rij rc,(7)rij=|ri rj|其中 是粒子 i 和粒子 j 之间的距离.传统的 DPD 模型里能量不守恒,只能对恒温的模型进行模拟计算.因此 eDPD 在 DPD 模型的基础上引入粒子内能变量,使 DPD 方程中包含能量项,流体粒子之间的热量传递通过可以下式表示21:CVdTidt=(qviscij+qcondij+qRij),(8)CVTiqviscij,qcondij,qRij式中,为粒子的定容比热容;为粒子温度;分别为黏性热通量、碰撞热通量和随机热通量,分别为qcondij=j=iij2(rij)(1Ti1Tj),(9)qviscij=j=i12CV2(rij)ij(eij vij)22ijm ij(rij)(eij vij)ij,(10)qRij=j=iij(rij)t12eij,(11)物理学报ActaPhys.Sin.Vol.72,No.18(2023)180203180203-2ijijeijeijeji其中,和 分别表示碰撞热通量系数和随机热通量系数,为随机热流中的随机变量,具有反对称性,即 =,保证了能量守恒,权函数 随着粒子之间的距离单调递减,超过截断半径之后就为 0,如(6)式和(7)式所示.ijij耗散力系数 和随机力系数 ,它们之间满足耗散涨落定理.当此关系得以满足时,系统的温度将保持在给定值 T 附近,这一功能也被称为恒温器,由(12)(14)式来实现:ij=2ij(Ti+Tj)4kBTiTj,(12)ij=2kBij,(13)ij=C2vko(Ti+Tj)24kB,(14)kBkoC其中 是玻尔兹曼常数,而 是用来调整传导系数的参数.在 eDPD 模拟中,(7)式的 s=2,s 也可以取其他值用来调节流体的黏性.与 DPD 一样可以通过对于热能的积分得到 eDPD 中各参数与实际热物性的理论公式.其中最重要的就是热传导系数 和热扩散率 .根据Ripoll22的研究推导可得C=1d(CVko0r3C(rij),(15)CVCV=CV/kBd其中 是无量纲的定容比热,;是维度.在三维情况下可得C=16CVkor5c(s+1)(s+2)(s+3)(s+4)(s+5).(16)=CCV而热传导系数满足 的关系,从而可得 的理论值.3结果和讨论3.1 eDPD 系统物性参数设置与热膨胀系数的计算在自然对流问题的模拟中两个重要的无量纲参数分别是瑞利数 Ra 和普朗特数 Pr,可以通过下式计算得到Ra=g(TH TC)H3CC,Pr=CC.(17)CC其中 g 为重力,b 为热膨胀系数,TH为高温温度,TC为低温温度,H 为计算域的高度,为运动黏度,为热扩散率,下标 C 表示该物理量为低温温度对应的物理量.Ra 的大小能判断流体的自然对流是否能产生,以及自然对流的强烈程度;Pr 则反映了流体中能量和动量迁移过程中的相互影响.在以往研究中 eDPD 系统被视为不可压缩系统,Ra 中的热膨胀系数 b 是基于 Boussinesq 假设的输入量.然而如前所述,根据 Mai-Duy19和 Pan等20的研究 eDPD 系统不能直接视为不可压缩系统.本文研究了 eDPD 系统的热膨胀系数 b,通过以下方式计算:首先构建 3 个方向为周期性边界条件的正方体计算域,通过调整计算域体积,实现不同平均温度下的等压环境,等压条件下正方体计算域初始密度 r0=4,粒子质量 m=1,比热容 CV=27800,保守力系数 a=18.75,随机力系数 s=,介观导热系数 k0=0.000126,系统温度 T 变化范围为 11.8,所有参数均为无量纲参数.压强设置为 35,等压条件下,密度 r 随着温度 T 变化的规律如图 1 所示,拟合后得到 r=0.59279T22.681T+6.2659.图 1 也显示了压强分别为 51 和 68 时,密度与温度之间的关系,表明压强对上述 eDPD 系统的热膨胀性影响不大.1.01.21.41.61.82.03.03.23.43.63.84.0513568=0.592792-2.681+6.26259图1不同压强下 eDPD 密度随温度的变化Fig.1.DensityvariationofeDPDwithtemperatureatdif-ferentpressurelevel.根据定义,热膨胀系数 b 由下式计算得到=1VdVdT=1ddT,(18)将 r=0.59279T22.681T+6.2659 代入(18)式,得到 eDPD 系统自身热膨胀系数 b 随温度 T 变化的规律:=0.3952T+0.8937.(19)b 在温度 T=11.8 的范围内变化趋势如图 2所示.物理学报ActaPhys.Sin.Vol.72,No.18(2023)180203180203-3基于泊肃叶流动,在上述系统参数条件下可以得到系统的动力黏性系数=1.077,对应的运动黏性系数 n=/r=0.265,与相同条件下的理论解23一致.3.2 eDPD 模型验证30 基于 eDPD 系统自身的热膨胀性,模拟了方腔中 RB 问题.如图 3 所示,方腔模拟区域为 x 30,15 y 15,2 z 2,z 方向施加周期性边界条件.上板为冷板,温度 TC=1.0;下板为热板,温度 TH=1.8;左右两板为互相平行的绝热板,W/H 始终保持为 2,方腔内流体粒子受到向下的重力,g=0.02.将绝热板的 a,s,k0都设置为 0 以实现绝热条件,并施加麦克斯韦反射边界条件从而防止流体粒子穿透.模拟时间步长为 0.02,每 5000 步对区域数据进行平均统计,总步数为 100000 步.通过改变温度和重力的大小以调整 Ra,Pr 始终保持为 1.分别得到 Ra 为 1600,1800,2000,4000时的自然对流温度云图和速度矢量图,如图 4 和图 5 所示.值得注意的是,模拟中设置的无量纲重1.01.21.41.61.82.0-1.0-0.8-0.6-0.4-0.200.20.40.60.81.0图2eDPD 系统热膨胀系数随温度的变化Fig.2.VariationofthermalexpansioncoefficientofeDPDsystemwithtemperature.冷板 C热板 H绝热板dd绝热板dd图3eDPD 模拟方腔内 RB 问题的模型Fig.3.ModelofeDPDtosimulateRBprobleminsquarecavity.-30-20-100102030-15-10-5051015(a)0.99901.0421.0841.1271.1701.2121.2551.2971.340-30-20-100102030-15-10-5051015(b)0.99801.0501.1011.1521.2041.2561.3071.3581.410-30-20-100102030-15-10-5051015(c)-30-20-100102030-15-10-5051015(d)0.99801.0981.1991.2991.3991.4991.5991.7001.8000.99801.0981.1991.2991.3991.4991.5991.7001.800图4eDPD 模拟不同瑞利数下自然对流的温度云图(a)Ra=1600;(b)Ra=1800;(c)Ra=2000;(d)Ra=4000Fig.4.TemperaturecloudsofnaturalconvectionunderdifferentRayleighnumberssimulatedbyeDPD:(a)Ra=1600;(b)Ra=1800;(c)Ra=2000;(d)Ra=4000.物理学报ActaPhys.Sin.Vol.72,No.18(2023)180203180203-4力应小于 0.1,无量纲温差也须小于 1,否则会使eDPD 模型的边界产生较大的滑移,并且使系统热膨胀性降低19,20,从而导致模拟的速度场、温度场与实际结果产生较大误差.基于 eDPD 模拟中温度云图与速度矢量图随 Ra 的变化情况,可以判断产生自然对流的临界瑞利数 Ra 在 16001800 之间.当 Ra 小于 1600 时流体的传热方式以热传导为主,Ra 大于 1800 时流体的传热方式以热对流为主.此时从能量守恒耗散粒子动力学的控制方程来分析,由于 eDPD 粒子之间存在随机力,具有随机效应.当传热方式发生改变时,从介观尺度来看系统的自然对流产生并非对称发生的,而是随机地从某处开始产生微小的环流.随着 Ra 不断的增大,小的环流带动整个流体最终形成稳定较大的两个环流,这也符合 Ra 可以反映自然对流的激烈程度的特点.本文得到的临界 Ra 的范围与 Koschmieder和 Pallas24实验得到的临界瑞利数 Ra=1675 相吻合,同时也符合 Abu-Nada15得到的临界瑞利数 Ra=1707 的结果;此外 Zhang 和nskog25通过朗之万方法基于热通量变化也得到临界瑞利数 Ra 为 1708 的结论.显示,利用 eDPD 自身的热膨胀性能够有效描述自然对流的发生.3.3 eDPD 模拟方腔内 RB 问题基于 eDPD 系统自身的热膨胀性,对不同Ra条件下方腔内自然对流 RB 问题进行了验证.在 3.2 节参数设置基础上通过调节重力加速度,将自然对流模型的 Ra 降低到 3100,Pr 始终保持为 1,模拟结果归一化后得到的温度云图如图 6 所示.10(a)0-10-30-20-10010203010(b)0-10-30-20-10010203010(c)0-10-30-20-10010203010(d)0-10-30-20-100102030图5eDPD 模拟不同瑞利数下自然对流的速度矢量图(a)Ra=1600;(b)Ra=1800;(c)Ra=2000;(d)Ra=4000Fig.5.NaturalconvectionvelocityvectorsatdifferentRayleighnumbersbyeDPDsimulations:(a)Ra=1600;(b)Ra=1800;(c)Ra=2000;(d)Ra=4000.1.11.21.31.41.51.61.71.01.11.21.31.41.51.61.71.80.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.80.10.20.30.40.50.60.70.80.91.0002.0图6eDPD 模拟 Ra 为 3100 时方腔内自然对流的温度云图Fig.6.eDPDsimulationofthetemperaturecloudofnatur-alconvectioninthesquarecavitywhenRais3100.物理学报ActaPhys.Sin.Vol.72,No.18(2023)180203180203-5FVM 针对偏微分方程进行数值求解,它将求解区域划分为有限数量的体积元,并在每个体积元内计算它们与相邻体积元之间的通量.本文使用FVM 模拟相同的自然对流做对照组,左右两边同样设置为绝热边界条件,上下板温差为 0.8,下板设置为热源,计算域的尺寸、液体的密度、黏性、比热容、热扩散率、热膨胀系数、重力加速度等参数都与 eDPD 系统建模完成之后所得到的参数一致,Ra 和 Pr 都是无量纲参数,通过(17)式计算后使FVM 的 Ra 也为 3100,Pr 也为 1.在此基础上每隔 0.04 个单位提取等温线,eDPD 的等温线与相同 Ra 的 FVM 模拟结果对比如图 7 所示,可以看到 eDPD 的模拟结果与 FVM 模拟结果误差很小,整体误差小于 0.7%,边界最大误差小于 1.8%.继续调整重力大小,其他条件保持不变,将Ra 提高到 7300,Pr 保持为 1,如图 8 所示对比关系,仍获得较好的对比结果,整体误差小于 1.1%,边界处最大误差小于 2.4%.在绝热边界处,eDPD的温度值略微偏离 FVM 结果,主要是由于 eDPD中无滑移边界的不完全满足所导致的.在此基础上再进一步对比 eDPD 和 FVM 的速度场.将速度场沿着 X0.5 和 X1.5 取 VX沿着 Y 方向分布,沿着 X=0.3 和 X=1.0 方向取VY沿着 Y 方向分布,得到 VX,VY速度场度场归一化之后与 FVM 的模拟结果对比,如图 9 所示.由图 9 可知,在不同 X 坐标方位上,eDPD 与FVM 的速度场同样吻合较好,由于无滑移边界条件不完全满足,从而导致速度场在边界处存在一定的偏差.此外,由于 eDPD 中的介观热波动效应,统计的速度存在一定波动,在速度极值处也有一定差异.该现象也存在于以往研究的 eDPD 自然对流模拟之中15,16.本文中 eDPD 模拟方腔内自然对流整体误差为1.1%,相较于 Abu-Nada 在相同条件下基于 Bouss-inesq 假设15,16的模拟误差已有明显改善.综上所述,不需要引入 Boussinesq 假设,基于系统自身的热膨胀性的 eDPD 可有效模拟不同 Ra 数下的自然对流 RB 问题.3.4 eDPD 模拟同心圆环和偏心圆环内自然对流通过同心圆环自然对流问题的模拟,进一步验证了基于系统自身的热膨胀性的 eDPD 对复杂几何结构中自然对流问题模拟的可行性.图 10 为同心圆环几何模型示意图,内圆半径 Ri=10,外圆半径 Ro=30,内圆温度为高温 TH=1.8,外圆温度为低温 Tc=1.0,其余模拟参数与 3.3 节中自然对流 RB 问题中的设置相同.同心圆环的 Ra 数和Pr 数由下式得到Ra=g(TH TC)(Ro Ri)3CC;Pr=CC.(20)调整重力 g 的大小和改变 Ro大小来获得不同 Ra 下的同心圆环自然对流,Pr 始终保持为 1,可以得到 Ra 从 100010400 时 eDPD 模拟结果与 FVM 模拟结果的等温线对比图,如图 11 所示,其中图 11(a)和图 11(b)中的 Ro=30,图 11(c)和图 11(d)中 Ro=40.由等温线对比结果可知,在同心圆环自然对流问题的模拟中,eDPD 与 FVM 结果在 Ra 数较小的情况下,如图 11(a)和图 11(b)所示,对比结果00.20.40.60.81.01.21.41.61.82.000.10.20.30.40.50.60.70.80.91.0YFVM1.11.21.31.41.51.61.71.11.21.31.41.51.61.7EDPD图7Ra 为 3100 时 eDPD 和 FVM 模拟的方腔内自然对流温度等温线对比图Fig.7.Comparisonof natural convection temperature iso-thermsinthesquarecavitysimulatedbyeDPDandFVMat3100ofRa.00.20.40.60.81.01.21.41.61.82.000.10.20.30.40.50.60.70.80.91.0YFVM1.11.21.31.41.51.61.71.11.21.31.41.51.61.7EDPD图8Ra 为 7300 时 eDPD 和 FVM 模拟的方腔内自然对流温度等温线对比图Fig.8.Comparisonof natural convection temperature iso-thermsinthesquarecavitysimulatedbyeDPDandFVMat7300ofRa.物理学报ActaPhys.Sin.Vol.72,No.18(2023)180203180203-6误差较小,整体误差小于 1.4%;当 Ra 提高后,如图 11(c)和图 11(d)所示,温度在圆环左右两侧产生了较明显的偏差,不过整体误差仍小于 2.7%,低温边界处最大误差小于 3.9%,原因在于此处的流体速度较大,边界条件的影响比较显著,导致统计过程中温度受到较大影响,不过整体模拟结果的精度相比相同条件下基于 Boussinesq 假设的研究也有了明显的提高21,26.同时,通过改变内圆位置模拟了偏心圆环的自然对流,如图 12 所示.类似于同心圆环模拟,温度在速度变化较大处仍存在较大偏差,该现象证明为了实现 eDPD 对自然问题的准确模拟,不仅需要考虑热膨胀性产生的影响,同时也应引入更完善的无滑移边界条件27.上述模拟对比表明,在使用 eDPD 进行自然对流模拟研究的过程中,系统自身热膨胀性的作用不应该被忽略;考虑热膨胀性可以减小以往研究中在 Boussinesq 假设下 eDPD 与 FVM 之间的模拟误差.最后,本文将对比探讨热膨胀性在传统基于 Boussinesq 假设的 eDPD 模拟中的具体影响.3.5 eDPD 引入Boussinesq 假设的自然对流问题模拟在上述研究的基础上,进一步讨论用 eDPD模拟自然对流时引入 Boussinesq 假设的适用条件.基于不可压缩条件,通常引入 Boussinesq 假设后粒子的运动方程由下式控制:dvidt=(fCij+fDij+fRij)+B(T To),(21)BgTTo其中热膨胀系数为 ,重力向量为 ,为流体温度,为初始温度.00.250.500.751.00-0.08-0.06-0.04-0.0200.020.040.060.08(a)FVMEDPD-0.10-0.08-0.06-0.04-0.02000.250.500.751.00(c)FVMEDPD00.250.500.751.00-0.08-0.05-0.0300.030.050.08(b)FVMEDPD00.050.100.1500.250.500.751.00(d)FVMEDPD图9Ra 为 7300 时 eDPD 和 FVM 模拟的方腔内自然对流速度线对比图(a)X=0.5 时 VX沿 Y 方向变化曲线;(b)X=1.5 时VX沿 Y 方向变化曲线;(c)X=0.3 时 VY沿 Y 方向变化曲线;(d)X=1.0 时 VY沿 Y 方向变化曲线Fig.9.ComparisonofnaturalconvectionvelocityprofileinthesquarecavitysimulatedbyeDPDandFVMwhenRais7300:(a)VariationcurveofVXalongY-directionatX=0.5;(b)variationcurveofVXalongY-directionatX=1.5;(c)variationcurveofVYalongtheY-directionatX=0.3;(d)variationcurveofVYalongtheY-directionatX=1.0.oHCi图10同心圆环模型示意图Fig.10.Schematicdiagramoftheconcentricringmodel.物理学报ActaPhys.Sin.Vol.72,No.18(2023)180203180203-7由于认为流体不可压缩,流体自身热膨胀性的影响在上述算法中被忽略.为提高 eDPD 模拟自然对流的精度,本文根据所需模拟的 Ra 数量级来判断是否引入 Boussinesq 假设:当按(18)式通过调整重力、温差、计算域大小,结合 eDPD 系统自身的热膨胀系数,就能得到所需模拟的自然对流的 Ra 数时,那么只需要充分考虑系统自身的热膨胀性即可;该可压缩条件下的 eDPD 自身的热膨FVM1.11.21.31.41.51.61.7EDPD1.11.21.31.41.51.61.700.13 0.25 0.38 0.50 0.63 0.75 0.88 1.000.000.130.250.380.500.630.750.881.00(a)00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.00.10.20.30.40.50.60.70.80.9(c)01.000.13 0.25 0.38 0.50 0.63 0.75 0.88 1.000.000.130.250.380.500.630.750.881.00(b)-50-40-30-20-1001020304050-50-40-30-20-1001020304050(d)图11不同 Ra 的eDPD 和 FVM 同心圆环自然对流等温线对比图(a)Ra=1000;(b)Ra=4100;(c)Ra=7200;(d)Ra=10400Fig.11.ComparisonofnaturalconvectionisothermsofeDPDandFVMsimulationforconcentricringsatdifferentRa:(a)Ra=1000;(b)Ra=4100;(c)Ra=7200;(d)Ra=10400.FVM1.11.21.31.41.51.61.7EDPD1.11.21.31.41.51.61.7(a)00.170.330.500.670.831.0000.170.330.500.670.831.00(b)00.250.500.751.0000.250.500.751.00图12eDPD 和 FVM 偏心圆环自然对流等温线对比图(a)Ra=4100;(b)Ra=10400Fig.12.ComparisonofeDPDandFVMeccentriccircularnaturalconvectionisotherms:(a)Ra=4100;(b)Ra=10400.物理学报ActaPhys.Sin.Vol.72,No.18(2023)180203180203-8胀性产生的浮升力已经可以满足相应自然对流的模拟,3.3 和 3.4 节中方腔内 RB 问题的模拟和同心圆环内自然对流的模拟均证明模拟精度已足够准确;所模拟的自然对流 Ra 数量级在 103以下时,不需要引入 Boussinesq 假设,仅仅依靠 eDPD 自身热膨胀性可完成相应自然对流问题的模拟.当所需模拟的 Ra 数量级增大时,根据图 1 可知系统自身所能提供的热膨胀系数有限,图 2 也表明温度的取值范围有限,无量纲温差不能超过 1,设置无量纲重力应小于 0.1,否则边界会产生较大的滑移.所以由(18)式可知只能扩大计算域的大小来得到所需模拟的自然对流的 Ra,但这会大大增加计算成本,为提高计算效率、节约计算资源,可以在考虑可压缩条件对 eDPD 系统自身的热膨胀性影响的同时引入 Boussinesq 假设,以此来实现更大数量级 Ra 数的自然对流模拟.以同心圆环自然对流为例,图 13 为按照以往研究简单把 eDPD 视为不可压缩系统,直接引入Boussinesq 假设,得到理论瑞利数 Ra 为 11600 的eDPD 模拟结果与FVM 模拟结果的对比图.由图13可知,此时 eDPD 结果即使在速度较小的区域也明显偏离 FVM 的结果,整体误差为 14.6%.由此可以判断在该数量级的 Ra 下使用 eDPD 模拟自然对流时,不考虑 eDPD 自身热膨胀性的影响,误差会比较明显,这是因为 eDPD 系统自身具有一定的热膨胀性,当通过自身热膨胀性所得到的 Ra大于自然对流产生的临界点时就会有自然对流产生,此时引入 Boussinesq 假设得到的 Ra 并非是当前自然对流对应的真实 Ra,需要综合考虑系统本身和引入 Boussinesq 假设的影响,才能真实反映当前自然对流所对应的 Ra,这也是直接引入Boussinesq 假设 eDPD 模拟结果与 FVM 误差较大的原因.随后在引入 Boussinesq 假设的基础上,考虑系统热膨胀性对模拟结果的影响.由粒子运动方程(21)和浮升力公式fbouyancy=gB(T To)(22)可知,Boussinesq 假设本质上给 eDPD 粒子施加一个体力以提供实现自然对流所需的浮升力,在eDPD 中这个力的施加方式和重力一致.要保持良好的边界条件,浮升力取值范围也应该和重力一样小于 0.1,同时如要满足无量纲温差不能超过 1,那B么可以得到 Boussinesq 假设中的 取值范围为g 0.1gB(T To)0.1 B 1.(23)由此可知 Boussinesq 假设极限条件下可以提供 310 倍于系统自身的热膨胀系数.基于上述分析可知,Boussinesq 假设产生的作用与 eDPD 自身热膨胀性的作用产生了叠加影响.通过计算可知,本节中 eDPD 构建的同心圆环模型自身热膨胀性对应的瑞利数 Ra1为 3400,在此基础上按照以往研究中引入 Boussinesq 假设对应的瑞利数 Ra2为 11600,最终 Ra 总和应为 15000.在相同条件下考虑热膨胀性的 eDPD 模拟结果,与 Ra 为 15000 的 FVM 模拟结果对比如图 14 所示.结果显示误差显著减小,整体误差降低至 3.2%.上述结果证明 eDPD 此时模拟的自然对流所对应的 Ra 更接近 15000 而不是 11600,说明模拟该数量级 Ra 的自然对流时,由系统自身热膨胀性产生的浮升力在实际计算中应该被充分考虑,否则模拟所得的自然对流所对应的 Ra 与预设需要得到的Ra 会有较大偏差,从而导致模拟准确度受到较大影响.由结果还可以继续推断出当采用 eDPD 继续模拟更大数量级 Ra 的自然对流时,可在考虑自身热膨胀性的条件下,通过结合 Boussinesq 假设实现大 Ra 的模拟;并且系统热膨胀性对整体热膨胀系数影响会随着所需模拟的 Ra 增大而减小;这种情况下引入的 Boussinesq 假设可以达到提高计FVM1.11.21.31.41.51.61.71.11.21.31.41.51.61.7EDPD01.0001.000.250.500.750.250.500.75图13不考虑 eDPD 自身热膨胀性时同心圆环自然对流等温线对比图Fig.13.Comparisonofnaturalconvectionisothermsincon-centricringswithoutconsideringthethermalexpansionoftheeDPDsystem.物理学报ActaPhys.Sin.Vol.72,No.18(2023)180203180203-9算效率的目的.综上所述,应当考虑 eDPD 自身的热膨胀性对整体自然对流的影响.4结论本文基于 eDPD 系统的热膨胀性,首先获得了 eDPD 系统自身的热膨胀系数 b,在此基础上模拟了不同 Ra 数下,方腔内 RB 问题、同心圆环和偏心圆环的自然对流问题.通过对 Ra 从 16004000之间的方腔内 RB 问题温度场和速度场的变化进行分析,显示流体产生自然对流的临界瑞利数 Ra在 16001800 之间,该结果与现有研究结果一致;之后定量分析了在 eDPD 系统中 Ra 数为 3100 和7300 时方腔内 RB 问题的温度场和速度场,并与FVM 的模拟结果对比,最大误差小于 2.4%.随后模拟同心圆环自然对流问题,将 Ra 从100010400 之间同心圆环自然对流的温度场与FVM 的模拟结果进行对比,最大误差小于 3.9%,在小 Ra 数下,本文通过 eDPD 自身热膨胀性得到的自然对流误差要远小于以往研究中的误差.通过分析偏心圆环的自然对流问题模拟结果,显示由于 eDPD 中的介观热波动效应,统计的速度存在一定波动,加上边界不能完全满足无滑移条件,导致模拟结果在边界处和极值处存在一定的偏差,未来需要进一步完善.最后本文探讨了在 eDPD 系统中自身热膨胀性对于引入 Boussinesq 假设的影响,并对相应的 Ra数进行修正.按照以往的研究,以同心圆环内自然对流模拟为例,eDPD 直接引入 Boussinesq 假设得到的瑞利数 Ra 为 11600,与 Ra 为 11600 的 FVM模拟所得温度场对比,发现误差为 14.6%;而考虑自身热膨胀性后引入 Boussinesq 假设,相同条件下的同心圆环自身热膨胀性所对应的 Ra1为 3400,引入 Boussinesq 假设对应的 Ra2为 11600,两者影响相互叠加,整体的 Ra 达到 15000,eDPD 模拟结果与 Ra 为 15000 的 FVM 温度场模拟结果拟合误差大大减小,误差降低至 3.2%.进一步,本文根据所需模拟的 Ra 数量级来判断是否引入 Boussinesq 假设:当所需模拟的自然对流的 Ra 依靠 eDPD 系统自身的热膨胀系数就能达到,那么只需要充分考虑系统热膨胀性导致的浮升力带来的影响即可,此时模拟精度已足够准确,不需要引入 Boussinesq 假设;当所需模拟的 Ra数量级增大时,由于系统自身所能提供的热膨胀系数有限,只能扩大计算域的大小,但这会增加计算成本,可以在考虑系统自身热膨胀性影响的同时,引入 Boussinesq 假设实现较大 Ra 的自然对流模拟;随着所需模拟的 Ra 增大,eDPD 自身的热膨胀性对整体自然对流的影响会减弱,这种情况下引入的 Boussinesq 假设可以达到提高计算效率的目的.综上所述,应当充分考虑 eDPD 自身的热膨胀性对整体自然对流的影响,以提高自