Fluent多孔介质与多孔压力阶跃面设置方法一文说清

Fluent多孔介质与多孔压力阶跃面设置方法一文说清

本文详细介绍了在Fluent软件中设置多孔介质（porous zone）和多孔压力阶跃面（porous jump）的方法，重点讲解了如何通过Excel进行数据拟合以获取惯性阻力系数和黏性阻力系数，并在Fluent中进行具体设置。内容专业且深入，适合对Fluent有一定基础的读者阅读。

01 引言

对于换热器、多孔格栅等在采用CFD方法计算时，划分网格存在困难，但又不可忽略其引起的压力损失，所以在Fluent中一般采用多孔介质porous zone计算散热器，采用压力阶跃面porous jump计算多孔薄板。（若不考虑换热，仅计算流场，也可以利用压力阶跃面计算散热器，将其简化为一个面。）

无论是多孔介质还是压力阶跃面，其理论基础相同，只是在Fluent中的设置略有差别。本文中我们假设所有的散热器都可以视为均匀的多孔介质。下面我们分多孔介质和压力阶跃面两部分分别说明如何在Fluent中设置。

02 散热器等采用多孔介质计算

1、风速压降数据准备

首先我们来看利用多孔介质计算散热器的Fluent设置方法，以下为某个散热器在不同风速下的压降值（该值一般通过台架试验测得）：

散热器厚度为0.1m，空气密度1.17kg/m^3，黏度为0.0000185kg/(m·s)。

2、理论基础

对于均匀多孔介质，可以将动量源项简化为下式：

式中ΔP为两侧压降值，单位为Pa；μ为流体黏度，单位kg/(m·s)；ρ为流体密度，单位为kg/m^3；v为流体的流速，单位为m/s；Δm为散热器或多孔板的厚度，单位为米m；C2为惯性阻力系数，1/α为黏性阻力系数；

我们将上表数据拟合为一元二次且截距为0的方程Y=ax^2+bx，a和b分别对应上式中的1/2ρΔmC2和μΔm/α两系数，然后就可以将黏性阻力系数1/α和惯性阻力系数C2算出来了。

3、Excel数据拟合

3.1 将以上数据复制到Excel，选中数据，点击插入折线图，并将横坐标改为对应的风速，如下图所示。

3.2 右键折线，点击添加趋势线，并将右侧设置为2阶多项式，将截距设为0，并勾选显示公式。

3.3 显示趋势线与公式如下，我们将风速2m/s代入公式计算得风阻应该为162Pa，显然与趋势线不符，公式是错的，这是因为显示的公式并非以2、5、8、10为x值进行拟合得，而是1、2、3、4m/s。那么该如何纠正呢？

3.4 我们右键横坐标的风速，点击设置坐标轴格式，选中日期坐标轴，此时折线图和公式会相应的发生变化，再次将风速代入验算，发现此时公式已经没问题了。

3.5 根据拟合公式系数计算C2及1/α将拟合出来的公式y=7x^2-5x与下式对应，可以确定1/2ρΔmC2=7和μ/α=-5，而根据散热器厚度Δm为0.1m，空气密度ρ为1.17kg/m^3，黏度μ为0.0000185kg/(m·s)可以确定惯性阻力系数C2=119.7，1/α=-2702702.7。

4、Fluent设置

4.1 利用scdm生成下图流体域，流体域中间位置为一方形换热器，左侧为流体域入口，右侧为出口，进行网格划分并按照常规方法在Fluent中进行流体域设置，我们着重看散热器区域的设置方法，如果需要了解其他部分的设置方法，可在文末获取源文件或评论区催更。

4.2 流体域沿Y轴正方向穿过换热器流动，双击Cell Zone Conditions中的heatexchanger流体域进行散热器的设置。

4.3 将散热器流体域勾选Laminar Zone，内部流通视为层流。并勾选Porous Zone，将散热器视为均匀多孔介质区域。将Porous Zone栏中的方向1设置为流体流通方向，即Y方向设为1，将方向2设为任意另一个方向。并且在方向1上将黏性阻力系数Viscous Resistance设置为上文算得的-2702702.7，将另外两个方向设置为极大值，可以乘以1000或10000（设置太大可能会报错）。然后将惯性阻力系数Inertial Resistance主流方向设为119.7，另外两方向同样可乘以1000，即流动及流通阻力都在主流方向上，而另外两个方向阻力极大几乎不可流通。

4.4 湍流方程、进出口条件等设置完毕，初始化并计算会发现无法计算，直接报错，这就是此前公众号（HVAC小李工）文章链接：解决报错Divergence detected in AMGsolver pressure correction/coupled所述的问题，即惯性阻力系数和黏性阻力系数不能为负值，因此我们需要改变一下Excel拟合公式的方式。

5、将惯性黏性阻力系数拟合为正值的方法

5.1 首先我们需要把Excel的规划求解器调用出来，方法如下：点击Excel左上角文件→点击选项→加载项→选中规划求解加载项→点击下方的转到；勾选规划求解加载项，并点击确定，此时在Excel的数据一栏中就会出现规划求解工具；

5.2 我们把一元二次方程的自变量x与因变量y放在A和B列，然后再C列中先填入一个二次项系数预测值（比如1），D和E分别填入一次项的系数预测值和截距预测值（也可都设为1）。然后根据一元二次方程再F列中计算每个预测y值，在F16中填入=$C$16A16A16+$D$16*A16+$E$16，式中$是为了保持C16单元格不变，通过向下拖动算出四组预测值。在G列中利用F16-B16计算预测值与真实值之间的误差，向下拖拽算出4组误差值。在H16中利用=SUMSQ(G16:G19)计算误差平方和。

5.3 点击Excel中数据菜单中的规划求解，设置目标选择H16误差平方和，并选择最小值，意为使误差平方和达到最小；在通过更改可变单元格中选择C16至E16中的系数和截差，意为此三个数可以变动；在遵守约束中添加C16和D16大于等于一个大于0的数字，E16=0，意为使二次项和一次项系数均大于0，截距等于0；点击求解；

5.4 选择保留规划求解的解并确定，这样就在原位置单元格中显示求解后的满足约束的系数值，并显示了相应的预测值、误差值；

5.5 如此可以得到风速压降拟合公式为y=6.442x^2+0.001x，系数均已为正值，根据散热器厚度Δm为0.1m，空气密度ρ为1.17kg/m^3，黏度μ为0.0000185kg/(m·s)及下式可以确定惯性阻力系数C2=110.12，1/α=540.5。

5.6 按照上文方法在Fluent中将散热器对应的计算域设置多孔介质参数，然后计算即可；

03 压力阶跃面

1、对于多孔薄板，当流体流过时会存在流阻，但划分网格也会存在困难，此时就可以采用porous jump多孔压力阶跃面来考虑多孔板引起的流阻。

与多孔区域porous jump的理论基础相同，在绘制网格时，将多孔板简化为一层面网格，在Fluent中右键该面，将该面调整为porous jump面，则会弹出以下小窗口：

该窗口内存在3个数，分别为①Face Permeability面渗透率，该值为上文中的α，即黏性阻力系数的倒数；②Porous Medium Thickness多孔区域厚度，即多孔板厚度；③Pressure-Jump Coefficient压力阶跃系数，即为上文中的C2惯性阻力系数；

2、通过仿真获取风速压降数据

上文提到我们在拟合公式之前，需要首先知道风速压降数据，可是有时我们并没有试验条件，那该如何解决？

对于下面右图中规则的多孔板，我们可以将其中一块划分出来，利用仿真方法计算出几组不同风速下的压降值。但应该注意，划分出来的部分应该能够通过复制代表整块板；

3、划分模型

上图多孔板每个孔由六个孔规则地围起来，各成60度角，因此可以将其中一个孔及六边形边缘划分出来，如下图，通过中心带圆孔的六边形复制即可代表整块多孔板。

对该单个带孔六边形划分以下流体域，左侧面设为速度入口，右侧面为压力出口，透明的六边形的六个长条面设置为Fluent中的对称面，设置在流体域中的带孔多边形表面设为基本的Wall表面。

在SCDM中创建上述几何体，然后进行命名，例如用鼠标左键选中左侧六边形表面，点击左侧Groups，然后点击Create NS，在下方填入inlet，同样的方法可以命名其他表面。（选中表面后Ctrl+G快捷键效果相同）

4、Fluent Meshing划分网格将该模型放入Fluent Meshing划分体网格，此文我们主要说明方法，网格质量我们这里不作讨论；

5、Fluent求解计算我们采用SST k-omega湍流方程，和Fluent默认的壁面函数等、入口速度设置为3m/s，监测出口最大速度不再变化视为收敛，计算结果如下，速度在孔口处最大达到7.21m/s，孔洞两侧的压降为18.2Pa。同样的方法可以获取入口速度为5、9、12m/s时的压降值，如下表：

6、利用上文中的方法，我们在Excel中进行拟合，拟合出的公式为Y=1.43x^2+0.1x，此时空气密度为1.225kg/m^3，黏度为0.000017894kg/(ms)，厚度为0.002m，通过对照下式可以算出面渗透率0.00000035788，压力阶跃系数为1167.3。

将流场中的多孔板简化为一个面网格，并设置为压力阶跃面Porous jump，将以上数据填入即可在流场计算中考虑因孔板引起的风阻了，感兴趣的同学可以将以上数据代入计算，如有错误，请在评论区指正。