一种乏燃料贮存格架附加质量和附加阻尼计算方法

1.本发明涉及一种方法及其应用，更具体地，涉及一种乏燃料贮存格架附加质量和附加阻尼计算方法，尤其是涉及一种基于谱方法和多孔介质方法的乏燃料贮存格架附加质量和附加阻尼计算方法，属于流固耦合领域。


背景技术：

2.第三代压水堆广泛采用了多孔式格架设计(如图1所示，装有乏燃料组件的格架密集的排列在乏燃料水池内)，传统格架的流固耦合作用研究是将格架等效为同心方管，并采用经验公式或解析解计算附加质量。这种方法计算便捷，但是，如果采用这种方法计算管端连接式格架，实际上忽略了小室间窄缝流动的影响，使得计算得到的附加质量偏大，附加阻尼偏小。目前，现有对此类管端连接式格架的流固耦合效应研究已经取得一些成果。刘雨等(《振动条件下平行孔板间流体的作用力实验》[j].原子能科学技术,2017,51(8):1419
–
1424.)开展机理性实验，将相邻格架的壁面视为孔板，研究了平行孔板和平行板的流体作用力，发现孔板的流固耦合作用力显著小于平板。li和yin(li w,lu d,liu y.study on the fluid added mass of cap1400 spent fuel storage rack under vibration condition[j].nuclear engineering and design,elsevier,2018,337(july):439
–
449.in t,lu d,huang y等.study on added mass model for fluid-structure interaction of spent fuel storage rack[j].nuclear engineering and design,elsevier,2019,353(july):110224)等针对cap1400格架开展模型实验和cfd数值模拟，发现此类管端连接式格架的附加质量显著小于传统方法计算的附加质量。然而，实验测量的成本过高，而基于精细建模的瞬态cfd流固耦合参数计算方法耗时过大，而且阻尼的计算结果与实验尚存在一定误差。近年来，研究人员采用dbd模型，研究了浸没多孔球体(felderhof b u.velocity relaxation of a porous sphere immersed in a viscous incompressible fluid[j].journal of chemical physics,2014,140(13).ollila s t t,ala-nissila t,denniston c.hydrodynamic forces on steady and oscillating porous particles[j].journal of fluid mechanics,2012,709(1948):123
–
148.vainshtein p,shapiro m.forces on a porous particle in an oscillating flow[j].journal of colloid and interface science,elsevier inc.,2009,330(1):149
–
155.)的流体作用力，并取得了一些进展，kang等(kang s j,dehdashti e,vandadi v等.optimal viscous damping of vibrating porous cylinders[j].journal of fluid mechanics,2019,874:339
–
358.)利用bdb模型，研究了多孔圆柱的流固耦合现象，发现相比于实心圆柱，多孔圆柱的附加阻尼在一些情况下要高几倍。而上述模型是圆柱形的，具有中心对称和轴对称特征，本质上是一维问题。而在方柱情况下，各个方向流场并无完全相同，圆柱模型无法解决此方面的问题，因此，需要提出方柱模型研究在不同方向上的流场差异。为研究多孔式格架的流固耦合作用，也需要方柱来精确模拟。为提升多孔式格架的流固耦合参数计算效率，本程序针对流体浸没的多孔方柱，结合谱方法，计算得到相对精确的附加
阻尼和附加质量，为工程设计提供参考。


技术实现要素：

[0003]
本发明的目的是提供一种用于乏燃料贮运格架附加质量和附加阻尼的计算方法，这种方法可以较为经济的计算多孔介质下的附加质量和附加阻尼系数，为工程设计提供支持和依据，其技术方案如下：
[0004]
一种乏燃料贮存格架附加质量和附加阻尼计算方法，该方法是基于谱方法的附加质量和附加阻尼计算方法，其特征为：包括如下步骤：
[0005]
步骤1：给出模型参数、假设和近似；
[0006]
步骤2：建立流体方程；
[0007]
步骤3：给出边界条件；
[0008]
步骤4：给出格架模型渗透率；
[0009]
步骤5：利用傅里叶谱方法计算速度场和水动力载荷；
[0010]
步骤6：结果分析，给出附加质量和附加阻尼。
[0011]
本发明还公开一种将上述乏燃料贮存格架附加质量和附加阻尼计算方法应用于第三代压水堆采用了多孔式格架设计中。
[0012]
有益效果：
[0013]
根据本发明提供的一种基于谱方法的多孔介质附加质量和附加阻尼的计算方法，可以计算出流体浸没的多孔介质的附加质量和附加阻尼，并给出速度场和水动力载荷场，计算结果比已有的计算方法更加准确，和已有的试验和模拟方法相比，计算过程更加经济，可以为工程设计提供支持和依据。
附图说明
[0014]
图1为高密度乏核燃料贮存格架示意图；
[0015]
图2为谱方法计算网格示意图(浅色区域为多孔方柱；深色区域为纯流体域)；
[0016]
图3为压力场分布(纵轴为压力)；
[0017]
图4为x方向(左)和y方向(右)速度云图；
[0018]
图5为速度矢量图
[0019]
图6为涡量图及速度矢量图；
[0020]
图7为附加质量(左图)和附加阻尼(右图)随参数的变化。
具体实施方式
[0021]
下面将结合本发明中的实施例,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
[0022]
基于谱方法的附加质量和附加阻尼计算方法，其特征为：包括如下步骤：
[0023]
步骤1：给出模型参数、假设和近似；
[0024]
给出计算模型，选定高密度乏燃料格架格架作为研究模型，这种格架可以等效成
一种方形的多孔介质。假设格架在高度方向上处处相同，即不考虑格架的三维效应，假定，格架在地震作用下振动的幅值相对于格架本身的尺度来说很小。最终，得到一种理想化简化模型：中间为多孔的方柱，方柱的边长为a，方柱外为一正方形流域，流域边长为l＝4a，振动幅值为0.01a。如图2所述(浅色方形部分为多空方柱，深色部分为方形流域)。该模型内部多孔介质采用brinkman
–
debye
–
bueche(bdb)模型计算，外部采用navier-stokes(ns)方程求解流场。
[0025]
bdb模型是在嵌入多孔物质中的球形粒子模型的基础上修正而来的，是1948年提出可以用来描述多孔介质流体流动的计算模型。在bdb模型中，通过修改darcy方程，解决了边界条件冲突和渗透率偏大情况下的问题。ns方程是描述一般的不可压缩流体动量守恒的运动方程。
[0026]
假定流体不可压缩，雷诺数较小时，将内外流场的守恒方程可以写成统一的形式。
[0027]
步骤2：建立流体方程：
[0028]
内流场即多空方柱流场和外流场即方形流域流场的守恒方程如下式所示：
[0029][0030]
其中，ρ为流体密度，是流体速度矢量，是方柱的速度矢量，μ是动力粘度，为时间；为压力场；χ是二进制的几何分布函数：位于方柱内部时取1，否则取0。非稳态bdb和非稳态斯托克斯方程与连续性方程在方柱边界处耦合起来。在无穷远处近似为0。上式中，
[0031][0032]
ω是振动的角频率；κ是渗透率；r是振动雷诺数，r表征了随着某点与方柱外侧距离增大时，方柱周围流体振荡衰减的快慢；β表征了流体穿透多孔结构的深度。其比值
[0033][0034]
该比值是振动多孔结构流体动力学特性的重要参数。本文的长度、时间、速度、应力、单位长度的力和力矩分别采用a，1/ω，aω，μaω/a，μaω和μaaω无量纲化。
[0035]
步骤3：给出边界条件：
[0036]
对于中间方柱，求解u和p的解析解非常困难，因此，本文采用傅里叶谱方法求解该问题。由于谱方法要求计算域必须是满足周期性，所以，我们假定振动方柱和周围的流体也是周期性的,当相邻方柱的距离很大时，可以视为无限大流体内的多孔方柱振动问题。
[0037]
给定柱体的振动形式，可以用复数的实部表示速度和压力，即带入方程()，可以得到
[0038][0039]
其中，ξ
i2
＝ir β2,
[0040]
步骤4：给出格架模型渗透率：
[0041]
多孔结构的渗透率是模型关键参数，需要计算得到。根据多孔介质理论，简单的均质多孔介质动量源项可以写成
[0042][0043]
其中，si为动量源项，κ是渗透率，μ是液体的动力粘性系数，v是速度，c2是惯性阻力系数，ρ是液体的密度。
[0044]
沿流动方向的单位长度的压降δp/l在多孔介质中可以拟合成二次函数，即
[0045]
δp/l＝av2 bv (6)
[0046]
其中，δp/l是沿流动方向的单位长度的压降，a，b是拟合得到的系数。
[0047]
因此渗透率为
[0048][0049]
为获得本模型所需的渗透率，建立与乏燃料贮存格架相同布置的流场模型，进行cfd计算。模型出口压力设置为0，入口边界采用速度边界条件，速度采用如下的周期性函数：
[0050][0051]
其中，x是输入运动的位移函数；vs是速度输入；ω是角频率；t是时间；a是幅度。
[0052]
然后，根据计算得到的入口-出口压差曲线和入口速度曲线，选取三个周期的曲线的峰值作为被拟合数据。根据(6)式，拟合得到b＝2.13
×
104,则渗透率取κ＝4.7
×
10-8
。
[0053]
步骤5：利用傅里叶谱方法计算速度场和水动力载荷
[0054]
首先，对(1)式取离散傅里叶变换，变形后可以得到如下表达式
[0055][0056]
其中k是二维波向量，k是该向量的模，k＝|k|；算符^表示傅里叶变换后的变量和项。对上式动量方程点积ik，可得
[0057][0058]
再将上式带回到动量方程，得到
[0059][0060]
其中，f-1
表示傅里叶逆变换；
[0061]
为求解上述方程，首先需要离散上述计算域(图3)为n
×
n的正方形微元。然后采用快速傅里叶变换fft去计算(11)式在微元中心(差值点)的变换/逆变换值。然后得到以u为变量的线性方程，该方程可以通过广义最小残差法(generalized minimal residual method，即gmres)求解。由于单个计算域的速度场可以代表整个流体域空间，所以边界条件的连续性可以自动满足。最终，流体对方柱的作用力为
[0062][0063]
其中,sj,χj，uj和rj分别表示表面积、几何分布函数，第j个单元重心点的速度和位置。
[0064]
通过计算，本文设置周期域的长度为l，网格数量使得l/a≥5，同时要求最大网格尺寸小于0.1a。其中，a为方柱正方形截面的边长。不仅如此，几何分布函数，通过高斯滤波器避免快速傅里叶变换出现的吉布斯振荡。
[0065]
根据附加质量和附加阻尼的定义，方柱流体力的幅值可以写为：
[0066][0067]
步骤6：结果分析，给出附加质量和附加阻尼
[0068]
典型压力场和速度场
[0069]
本节选取一典型工况，展示计算结果。如图2所示，边长a＝1m，流体域边长l＝4a。方柱振动方向沿x方向，振动周期为2hz，幅值为0.01m，方柱最大速度无量纲后为1。渗透率取κ＝4.7
×
10-8
流体粘性系数μ＝1.0
×
10-3
，r＝1.256
×
107，δ＝1.296。程序运行时间0.46s，远小于cfd瞬态计算程序。
[0070]
计算得到的压力结果如图3所示。压力计算结果表明，在峰值时刻，多孔方柱的一侧为压力峰值，另一侧为压力谷值。图4展示了x方向和y方向的速度云图，x方向的无量纲速度最大值约为0.5，约为方柱振动速度幅值的一半，出现在多孔方柱内部；y方向的速度最大值出现在方柱的顶点处。合成的速度云图见图5，速度云图表明方柱上下两侧形成了环流，在计算域边缘的速度逐渐趋于0。图6给出了该方柱的涡量场，该涡量图表明，涡量场在方柱上下两侧呈带状分布。
[0071]
多孔方柱流固耦合参数随渗透率的变化
[0072]
振动雷诺数r和参数β是评估多孔方柱的主要无量纲参数。根据式()振动雷诺数主要由振动频率决定；β由多孔介质的渗透率决定；而参数δ是二者的比值，表征了外部流体振荡穿透深度与流体穿透多孔结构的深度的比值，δ趋近于∞时，表示外部穿透显著大于内部穿透，即流体几乎无法穿透多孔结构。
[0073]
本节在恒定的r下，选取了不同值(即不同的渗透率)，获得了附加质量
和附加阻尼参数的变化曲线(图7)和数据(表1)。由下图可知，附加阻尼随着δ的增长先增加，当δ接近1时，附加阻尼达到最大值取，随着δ进一步增大，附加阻尼逐渐降低并区域稳定。附加质量随着δ的增长先增加，当δ超过1之后，趋于稳定值0.76。
[0074]
表1计算工况及结果(r不变、δ变化)
[0075][0076]
本发明针对振动条件下多孔方柱所受流体作用力建立了dbd-ns耦合模型，结合谱方法，采用matlab求解，获得了多孔方柱的流体作用力，并得到了附加质量和附加阻尼等参数。本发明使用的谱方法可以用于求解无限大流体域内多孔方柱的流体作用力，并用于流固耦合参数的计算，计算效率远远高于基于cfd瞬态计算的方法。
[0077]
在以上的描述中阐述了很多具体细节以便于充分理解本发明。但是以上描述仅是本发明的较佳实施例而已，本发明能够以很多不同于在此描述的其它方式来实施，因此本发明不受上面公开的具体实施的限制。同时任何熟悉本领域技术人员在不脱离本发明技术方案范围情况下，都可利用上述揭示的方法和技术内容对本发明技术方案做出许多可能的变动和修饰，或修改为等同变化的等效实施例。凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均仍属于本发明技术方案保护的范围内。