无论是 DGT 靶还是 HLM 靶,其基本原理都是通过流动将沉积在束流反应区的热量进行移除,并携带到换热装置中进行热量的交换。因此对于靶系统的设计,对于流动特性的掌握是换热问题的关键所在。
流体的最基本组成是微观尺度的分子,但是当着眼于宏观尺度的问题时,一般是通过以包含大量分子的流体元为最基本单位,可以将分子的行为在一定程度上忽略,同时认为系统的适用连续性假设,认为流体的宏观特征量在空间中是连续均匀的确定值。对于液态靶中水或者 HLM 工质来说,分子自由程非常小,能够充分满足流体作为连续介质的前提条件。此外对于这两种流体,在研究的状态范围内,也非常充分的满足不可压缩流体的近似条件即工质的密度不随压强而改变,这些流体的条件性质不仅简化了问题的计算和设计,同时也为对于理解靶系统中的问题提供了重要基础。而颗粒体系的力学特性则要复杂的多,尽管曾有研究者试图以类流体模型对颗粒体系运动进行处理,但目前还没有好的普适方法,而且即使采用类流体模型,颗粒体系的流变特性与通常可以被认为是牛顿流体的水、液态金属也是截然不同的。对于通常状态下的水和液态金属,其由于流动速度的梯度所产生的切向应力与速度梯度成固定正比,这一比例即使流体的粘性。
对于靶问题中的多数物理问题,无论是在模拟与计算过程当中,还是在通过实验对体系物理量的测量方面,通常使用的都是欧拉方法,即关注固定空间坐标位置的物理量变化。确定流体运动的方程一般有连续方程和动量方程组成:
连续性方程:
式中(u, v, w)为速度的 x,y,z 方向分量,Fb为体积力,下标表示 x,y,z 方向分量,Pij为体积元所受应力张量中的对应分量。式中第一项为随体导数:
对于不可压缩流体:
因此,有:
描述流体运动的动量方程为:
应力张量的定义法向是朝向控制体外的,因而考虑主应力作用时,需要进行反号。考虑到本构关系方程后,靶材料这样的粘性不可压缩流体最常见采用的特殊化方程即是 Navier-Stocks 方程:
在直角坐标系下的分量形式:
由于最后一项粘性项的存在,导致了这一方程是非线性方程,也就决定了流体体系属于典型的非线性复杂系统,即使对于粘性为一定值的牛顿流体,也难以通过计算手段给出准确的系统状态,通常只能使用简化方法进行处理,或者采用数值求解方法。
一般来说,在确定的粘性和初始流动速度的情况下,方程存在定常和周期性的解,这个解一般也是比较光滑的,这正是流体中的层流现象。但当边界条件发生扰动时,正是由于粘性非线性项的存在,使得方程出现两种截然不同的性质:一是通过粘性耗散,初时的扰动逐渐衰减,流动能够恢复成稳定的状态;二是由于非线性作用,完全失去稳定,这时就成为了通常所说的湍流状态。由于粘性和流动速度间关系的重要作用,定义了无量纲系数雷诺数。
通常当雷诺数超过一定的阈值,在外界的扰动下就会发生湍流现象,当湍流发生时不仅导致压力速度的波动,甚至还会诱发特定的震动现象,最典型的即卡尔曼涡街。
除了上述非线性系统有关的不稳定性现象,水和液态金属作为不可压缩流体还会产生长距离的冲击波效应,最常见的是水击现象,冲击波的瞬时压强可以超过正常值数十倍,造成振动、噪声,严重的还会对于回路中存在缺陷和脆弱的部分产生巨大的破坏性效果,这一作用在管道的水力学研究中一维化的可以用以下方程表达:
连续性方程:
运动方程:
式中 H 为测压水头,v 为液流速度,s 为延管道长度,c 为水击波速, 为管道倾斜角度,g 为重力加速度,D 为管道直径。方程中可以看出,在较粗管道中,当管道系统沿程阻力不大时,水击波的运动可以表现出反复振荡的特点,同时在长距离上的衰减速度很慢。而计算研究和实验结果都表明了以上的现象,这样系统由于泵、阀门的突然剧烈的工况变化,会导致系统严重的振动冲击,对连接处和工作设备造成系统性的巨大破坏。
相比于流体,颗粒体系的基本构成是单个的颗粒,其尺度通常为宏观水平,而将颗粒体系作为整体研究时,既可以从颗粒之间两两作用力为出发点进行离散方法的研究,也可采取类似流体的等效连续模型从而避开体系内巨量的相互作用。使用离散的方法思路比较直观,其基本方程就是两体之间平动转动的运动方程。同时实际中的两体间非弹性相互作用模型也有利于理解和分析作用力的耗散。但其复杂之处一正是在于由于颗粒的材料、外形特性导致的复杂作用力形式,二在于体系内大量的颗粒关系对计算量的需求。离散计算模型中最为常见的是 HertzMindlin 接触模型,在这一模型中,法向作用由 hertz 模型提供:
切向作用由 Mindlin 模型以增量形式给出:
上两式中, 为剪切模量, 为泊松比, 为等效半径定义为R1 R2 /(R1 + R2), δ为形变重叠量。
当体系的规模超过离散处理能力的时候,采用连续性模型在一定条件下也是可行的。问题正在于体系所处的条件状况适用于何种连续力学模型,尤其是材料流变性质和耗散性质。相比于一般的连续性介质,其本身能够承受一定的应力这是与液体不同的;而受到应力作用时,其形变又是非线性和不可逆的,这一点又与固体不同。波在穿过固体堆积颗粒的过程中,不仅表现出波的向前传播,更重要的还同时存在波的反射、衍射等现象,对于高密度的堆积体,振波难以向前穿透,这样一来,局部产生的振动影响,主要都局限在振动产生的部位,不会产生系统性的长距离和持续效应。
颗粒中力学波的传播是极其复杂的形式,目前对这一问题仍缺乏全面的了解。一种常见的思路是基于上面介绍的 Hertz-Mindlin 模型,以平均场近似将其应用到颗粒系统整体,常见的等效介质理论 EMT(Effective Medium Theory)认为颗粒体系的性质如下:
上式中 K 为弹性模量, 为剪切模量,φ为堆积率,Z 为平均接触数,P 为堆积体挤压力,其它符号与 Hertz-Mindlin 模型中字母含义相同。
由于颗粒难以传播拉伸作用,只能讨论压缩作用和剪切作用,其波速分别以Cp ,Cs 表示:
因此,波速应当与压缩作用力的 1/6 次方成比例关系。但实验表明这一结论仅针对特定体系成立。但无论如何,通过正相关关系表明,对于并不处于完全紧密接触作用的体系,作用力减小,对于体系的波传递是更为困难。即是针对沙丘鸣沙表面振动的研究中,除了波速远小于气固声波速度外,其向系统内衰减也是非常迅速的。
应当说在不同的条件下,颗粒体系呈现出的特性是介于固体和液体之间的,但又与二者不同。设计过程中应当充分注意到这一效应,尤其是当颗粒体系的运动介于一定范围内时,呈现出固体液体转变的特征常被研究者称为玻璃态或jamming 状态。在临界范围内,颗粒流动也经常出现间歇和不稳定的情况,是设计过程中应当注意的因素。
尽管在数值模拟方法和一些设计的具体问题上,经常会采用离散模型的观点,但通过连续介质的观点来看更有利于将 DEM 靶和 HLM 靶进行对比。针对颗粒流动所建立的连续介质力学有时也被称为颗粒固体流体动力学,对于运动速度不快、堆积系数较大的颗粒体系当中低频长波动力行为物理图景的理解有着指导意义。当采用类似流体力学的观点看待颗粒体系时,其基本连续性方程和运动方程形式上是相似的,但由于连续性和内部作用力不同,因此在关键项上有区别。
连续性方程:
式中左侧第一项密度变化显然在颗粒体系中不为 0。因此,颗粒体系的连续方程需要使用:
运动方程同样满足,但是由于本构关系复杂,却无法做出化简,需要根据实际情况确定应力张量,除了基本的摩擦关系用于处理主应力和其他应力的关系,其它目前常见的流变模型还有 Mohr 环方法,或者以类似宾汉流体的处理,还有一些新方法试图采用高分子力学当中的模型,这些模型都在不同情况下取得了一定的成功,但在通用模型上还没有很好的解决办法。
但对于颗粒流动体系而言,由于其基本单元是有具有确定运动状态的宏观物体集合,因此也会以统计力学的方法借用热力学观点,对这一体系进行处理,例如体系中能量的耗散过程。
另外作为一种固相与气相的混合系统,对于固体气体的相互作用关系也应当在予以考虑,尤其是当固体颗粒尺寸密度较小而气体流速较大的时候,这一作用尤为显著,可以产生流化效应,甚至是堵塞情况产生,但针对靶设计环境而言,这一作用可通过系统材料、运行状态的设计予以避免。
就具体以数值模拟方法而言,HLM 靶的设计方法是 CFD 方法。在对流体行为的模拟中,总是需要通过网格进行离散便于进行数值模拟的。由于方程属于非线性方程,求解困难,对于这一问题有不同的解决思路:1.直接求解非线性 N-S方程,称为直接数值模拟方法,这一特点用于细密的时空离散的时候能够给出信息丰富的多尺度湍流信息,缺点则是对于这样的离散需求的计算量非常巨大,一旦网格粗糙,收敛性就会严重变差。2.采用时均-脉动的观点处理湍流,即雷诺平均方法,并且认为存在一个湍流造成的粘性系数即波西尼克假设,这样一来,运动方程得到了极大的简化,于物理过程而言忽视了就是湍流的细节结构,主要关注平均量。这样一来体系能够处理更为粗糙的时空离散和更大的系统。Fluent 中主要采用的即这一方法,其内设了多个不同的湍流模式可供使用。也是最为广泛采用的方法。3.介于上述思路之间的大涡模拟方法,这一方法对于网格尺度以下的湍流涡行为使用湍流模式的方法,而网格尺度以上的则类似采用直接数值模拟方法,得到的结果自然也介于两者之间。
对于 HLM 靶而言,以雷诺平均的方法进行计算时所采用的模型对于计算结果无疑对于流动传热有着重要的影响,MYRRHA 项目曾经在 CFD 模型方面做出了大量的工作,使用了数种 CFD 的软件和多种模型对于非 Detachedflow 靶的流动状态进行了模拟,并进行了 benchmark 实验。最终模拟中采用了多相模型(以 VOF 为主要手段并有其他如动网格或欧拉-欧拉法方法)、湍流模型(以k-ε 模型为主,LES 等其他手段进行辅助验证对比)、汽蚀模型(cavitation 模型为主)等。系统性研究方面,则采用了 Relap5 对额定状况开机与稳定运行,主流速度的小幅波动,主泵在束流关闭或未关闭情况下停转,束流开关,热交换器束流关闭或未关闭情况下失效,feeder drag 波动等情况进行了模拟与分析。
上海交通大学团队也曾经了对于非 Detached flow 的流动情况采用了不同模型进行模拟并与实验结果进行比较并给出他们的推荐方案。这些努力无疑是有益的,尤其是在大涡模拟的效果方面,但是客观的说,具有自由界面的特点使得非 Detached flow 本身的回流区的变化受到系统运行状态影响非常严重,几乎找不到一个稳定状态用来准确的进行比较。由于实验本身的波动不稳定一方面由于靶区流体流动本身造成,另一方面也与系统内其它组件导致的波动有关。这种不稳定的情况是 HLM 靶中由于物理机制自然存在而无法避免的现象,这一点无论在本工作中建立的水回路、上海交通大学水回路、还是 MYRRHA 项目的实验结果都是同样的,甚至通过三维模拟也能发现类似的不稳定现象。
计算模型方面三维模型能够更好的符合物理事实,表现流动的不对称性,但相比于二维算法,三维模型存在计算时间长,网格数量多或网格尺寸大,不易收敛等特点。如果在研究过程中针对靶段流动的参数影响,需要扫描参量,而对界面的瞬态波动关注较低,优先选择二维模型仍然是很有意义的,同时,经过与实验的对比,二维模型在流动表现形态方面是能够满足要求的:
目前模拟状态对于界面较稳定的部分的模拟结果基本还是准确的。在Detached Flow 状态下由于回流区区域显著的减小了,对于剩余较稳定界面部分的模拟价值无疑就更高。可以作为进一步 HLM 靶的研究方向之一。对于小回流区的情况,相比于弄清楚对于流体过程,也可以从工程角度通过环形束流更为方便的避开回流区的热量沉积。
因此,本工作在进行 MYRRHA 构型无窗靶的进一步参数研究之前,采用Detached Flow 状态的靶流动对可能采用到的模型、离散、计算方法共同进行了比较。模拟中在瞬态/稳态模拟、piso/ simple 算法、k-ε/k-Ω 湍流模型、网格有效性方面进行了对比研究。对比中其它条件为 VOF 多相流模型、cavitation 相变模型、液态金属入口流速 1.0m/s、出口压强 18000pa、相变临界压强 1pa、参考压强 0pa、重力 9.8m/s2、计算步长 0.00005S、计算总时长 5s。基础模型如下图所示,其尺寸按照 MYRRHA 项目 V0.10 靶型进行设置,共有计算网格 63980 个。模型采用 2D 轴对称构建,上部内管设为压强入口,模拟加速器束流管道真空环境;内管与外壁形成部分为速度入口,模拟以确定流量流入靶区的 HLM 工质;下部为压力出口,模拟流出靶区的工质的限制条件。模拟过程中虽然对入口部分的进行了增加阻力设置,但是主要评估研究区域为锥形汇流段。
在各个算例当中,HLM 工质在靶区的流动均形成了预期的自由液面,其相图和压力分布如下所示。
上图的结果均整体上与文献结果呈现基本一致的特征,尽管在局部的有边界条件造成的区别。由于在针对靶区回流现象和自由喷射流进一步的研究中,这一结构沿对称轴的流动是非常重要的,因此,将算例中沿对称轴的流速特性进行对比,如下图所示
结果表明,几种设置在这一问题的模拟上所得到的主要结果基本是一致的。另一方面,网格相关性所产生的结果变化也不显著,总体来说,速度方面,采用各个方法所得到的结果在作为研究重点液相区域区别不大,考虑到压强分布,尤其是发生汇流作用的区域的不同,认为二维轴对称下网格尺寸为 1mm 的计算模型和 K-ε 湍流模型能够起到良好的计算效率与计算质量的折衷,同时,采用稳态计算研究液相区域的流动,也是可以接受的方案。在这些算例在计算中所花费的四核并行计算时长如下表所列:
进一步地,使用稳态计算对计算模型中的离散方式进行了对比,对比算例的设置如下表所示,计算迭代为 10 万步,并且由于高阶算法的收敛问题,将动量松弛压因子调为 0.2。
除了 MP 组合导致结果不收敛,VB 组合导致显著不符合物理事实的结果外,采用不同的离散方法时,束流作用界面的形状和位置基本不变,主要发生变化的是位于下方的第二液面位置和两相界面的清晰程度,采用高阶相体积分数算法的算例中,两相界面明显更为清晰;在体积分数离散方法上 Compressive 产生的的两相界面最为清晰。得到的轴向压强分布如下:
图 5.9 对不同网格、压力、动量、体积分数、湍动的设置计算得到的轴向压强
图中,采用高阶动量算法,得出的压强值明显降低,而采用高阶湍流算方法在产生压力峰值的汇流点处的压强值明显较高。0.2m 以下的部分出现了不同的压强分布,采用不同网格离散的结果差距不大。总体来说,基于不同 1 阶精度算法得到的压强值在压力峰值处为(19895.75±4679.45)pa,变化幅度 23.5%,位置为(0.2530±0.002)m。
由于采用高阶算法的准确度高于一阶算法,因此,在稳定性可靠的情况下,尽量考虑在二阶算法中选择,提高计算精度。从中组合选取高阶计算方法进行比较,在比较过程中,由于高阶动量算法产生的差异最大,因此,以高阶动量算法为主,测试两种高阶算法共同使用时产生的效果,采用不同离散方法所产生的靶区气液两相分布、轴向流动速度分布如下所示:
第二界面的位置不同,反映了在这一区域的计算过程中不同的算法计算得到的压力不相同,为进一步量化这一因素,提取沿对称轴的压力分布进行比较,如下图所示:
与图中 SD 曲线比较可见,当多种离散方程复合使用的时候,计算值均较一阶时降低,各算例间的差距也都相应减小,可见高阶动量方程仍起到了主要作用;总体来说高阶湍流方程表现出与低阶时计算结果的不同,而高阶湍流算法之间的差异较小,如峰值压强变化减小为(14735.55±2038.75)pa 变化幅度 13.8%,压力峰值的位置变化减小为(0.25654±0.0015)m,另外还需要注意,在 0.1m-0.15m 处的压力分布也有着较大的不同,最大变化(2538.945±1197.135)pa,变化幅度为47.1%,这一位置附近是上方向下流动的液态金属注入第二液面的作用区域;高阶相体积分数离散方程也表现出与湍流离散方程相似的特征,为进一步验证,选取多个高阶离散,采用二阶动量下不同湍流、相体积分数离散方程进行对比。计算设置及计算结果相图如下:
上图中,随着更多的采用高阶算法,计算结果基本区域一致,其压强峰值波动 为 (15501.75±234.25)pa , 波 动 范 围 1.5% , 压 力 峰 值 的 位 置 变 化 为 (0.25754±0.0005)m,并且 0.1-0.15m 处的压力变化趋于一致。
综上,在进一步的计算中,选择二阶动量方程、二阶湍流方程、可压缩体积分数方程,这样的设置组合能够保证:针对靶区不同构型的两相流动状态保持较好的收敛性;提供较好的精度与计算速度的折衷;在该物理模型下,采用不同的算法导致的结果偏差不大,在预计误差 30%的情况下,可以涵盖 Fluent 不同算法导致的误差范围。
对该设置进行计算迭代步数结果偏差的测试和稳态与暂态算法的比较:
该模型下,到达 10000 步时,其各项残差曲线如下图:
各项残差随计算时间保持收敛,该模型下,随迭代步数的增加,得到的轴线压力分布如下图所示。
自 20000 步起,从其结果来看已经可以认为处于基本稳定。
而在颗粒靶的设计过程中采用的 DEM 方法,相比于 CFD 计算,由于基本作用原理简单,它的计算方法中使用的近似则要少得多。做一个不一定很恰当的类比,DEM 方法实际上就是与 CFD 方法中的直接数值模拟类似,而主要问题一方面是颗粒体之间的作用力形式问题,另一方面是计算规模问题。
对于前一问题,在 DGT 靶设计使用前述的 DEM 软件中接触力学采用广泛接受的 Hertz-Mindli 模型,并加入阻尼作用项,该模型用于计算颗粒间的受力,仅需输入颗粒的材料参数,如直径、密度、杨氏模量、恢复系数、摩擦系数就可计算任意两个接触颗粒间的受力情况。应当说,这种方法比之流体计算的方法要“第一性”的多。
针对后一问题,巨大的计算量可以采用大规模并行计算来解决,由于颗粒在不断的运动,其周围接触的颗粒、颗粒的运动性质都是在时刻变化的。在实际计算中通过通过采用元胞列表法和邻居列表法结合对于颗粒内部空间的巧妙分块,使得每个颗粒在寻找相互作用的其他颗粒时仅在其临近的域内寻找,对于其中计算量最大最为巨大的是计算每个颗粒的加速度和角加速度计算,利用 GPU 的众多轻量处理器实现,可有效降低计算量并可极大的提高并行速度。