基于OpenFOAM的猪体空气阻力系数模拟计算

吴雪飞 李 浩* 施正香

(1.中国农业大学 水利与土木工程学院，北京 100083；2.农业农村部设施农业工程重点实验室，北京 100083；3.北京市畜禽健康养殖环境工程技术研究中心，北京 100083)

集约化生产中，畜禽动物对舍内气流的阻挡影响了动物活动区域甚至整个舍内的温度、湿度和污染物浓度等环境参数[1-5]，同时也阻碍通风系统的调控以及舍内气流的组织。因此研究畜禽动物气流阻力变化规律对畜禽舍内环境参数的调控和通风系统的运行有重要意义[6-7]。

气流流速和流向、猪的体重和姿势等对猪舍内气流传热的影响已经得到研究[8-11]，但目前关于气流流速和流向等因素对猪体气流阻力影响的研究仍然较少。同时由于畜禽动物几何形状不规则及动物活动性的影响，使得影响猪群的气流阻力的因素较多，畜禽动物周围的气流变得非常复杂[8]。而单只猪的气流阻力影响因素相对较少，以单只猪气流阻力的研究作为猪群阻力研究的基础。

随着计算机的快速发展，基于计算流体力学(Computational fluid dynamics，CFD)的数值模拟技术在畜禽舍通风系统中得到越来越多的应用[12]。这些CFD软件分为商业软件与开源软件2种：商业CFD软件，如ANSYS和Star CCM+，具有便利的图形化界面，计算精度和稳定性较高，得到广泛使用，但是商业软件闭源封装计算代码，限制了使用者二次开发，同时高额的证书费用增加了使用成本；开源CFD软件，如OpenFOAM，因其免费的证书、灵活的程序构架、良好的拓展性等优势也得到了一定程度的推广应用[13-14]。近年来，OpenFOAM在农业通风相关的CFD模拟应用中也有涉及[15-16]。然而，OpenFOAM计算稳定性低，软件操作复杂(直接使用C++语言指令进行操作)，在农业行业的应用依然非常有限[17]， 尤其是针对农业环境领域开发准确性更高的数值模拟算法的潜力没有得到发挥。因此本研究使用OpenFOAM软件数值模拟单只猪对气流的阻挡作用，探讨开源软件OpenFOAM在农业通风研究中的可行性。

综上，本研究旨在基于开源CFD软件OpenFOAM数值模拟技术，研究单只猪对气流的阻挡作用，探究气流流速和流向对猪体空气阻力系数的影响规律，以期为猪群阻力研究提供基础理论依据，为开源CFD软件在农业建筑通风领域的深化应用提供支撑。

1 研究方法

1.1 几何模型与计算域

1.1.1单只猪模型和圆柱体的构建

本研究参照体重约为90 kg的育肥猪(图1(a))进行建模，为了方便建模和减轻网格计算负担，略去猪的四肢、耳朵、尾巴，得到简化猪体模型(图1(b))，该简化猪体模型尺寸约为：长1.17 m、高0.48 m、宽0.27 m，体表面积1.15 m2，其中简化猪体模型代替全尺寸猪体模型进行模拟的准确性已经得到证实[10，18-19]。由于活体猪数据采集困难，虚拟风洞中模拟圆柱体作为模型验证的方法简单快速[20-21]，因此本研究利用圆柱体阻力系数半经验公式计算所得的圆柱体阻力系数与模拟所得的阻力系数对比从而验证模型的准确性。为了确保圆柱体与简化猪体模型空气阻力系数的可比性，圆柱体与简化猪体模型体积相同；同时为保障圆柱体阻力系数半经验公式的适用性，圆柱的长度与直径之比为4.0[22]，因此构建直径为0.37 m，高度1.48 m的圆柱体(图1(c))。

图1 实际猪模型，简化猪体模型和圆柱体关系示意图Fig.1 Diagram of the relationship between the actual pig model, simplified pig model and cylinder

1.1.2计算域及模型位置确定

本研究构建一个尺寸(长×宽×高)为10.0 m×6.0 m×2.4 m的虚拟风洞，虚拟风洞设置进风口(Inlet)和出风口(Outlet)各1个，壁面4个。猪体中心距离进风口4.0 m，出风口6.0 m，满足单一物体距离风洞进风口大于物体特征长度的5倍(2.4 m)、距离出风口大于物体特征长度的10倍(4.8 m)的要求，距离壁面约为2.5 m,顶面约2.0 m,猪体的最大风洞阻塞率为1.8%，低于3.0%的设计要求[23]，距离地面0.5 m(猪只站立高度)。

本研究中圆柱体的位置因为试验目的不同而在虚拟风洞中的位置不同：为了验证模型的可靠性以及模拟结果的准确性，将圆柱体布置在虚拟风洞高度的中间位置(距离地面1.2 m)，圆柱体与其余壁面的距离与简化猪体模型相同，从而保证圆柱体空气阻力系数半经验公式的测试条件；为检验圆柱体作为极度简化猪体模型的可行性，圆柱体也布置在简化猪体模型的相同位置进行模拟，以保障其与简化猪体模型空气阻力系数的可比性。

1.2 气流流速与气流流向

研究不同流速下猪体对气流的阻挡作用，虚拟风洞进风口的气流流速v设置6个水平，分别为0.2、0.5、1.0、1.5、2.0和2.5 m/s。气流流向的改变通过猪体与气流流向的夹角θ进行表征，θ定义为猪体长轴方向与虚拟风洞进风口气流流向的夹角。本研究中，风洞气流流向不发生变化，仅改变猪体的朝向来实现θ的调节，θ设置13个水平，分别为 0°，15°，…，180°，其中θ=0°，90°和180°分别代表猪体头部垂直迎风，猪体侧面与气流流向垂直和猪体尾部垂直迎风的情况。

1.3 空气阻力系数定义及验算

物体受到的气流阻力由几何形状引起的压差阻力和流体粘性引起的摩擦阻力构成。工程上为了方便计算绕流阻力引入空气阻力系数(Drag coefficient，Cd)。本研究使用猪体表面平均空气阻力系数(Cd)量化猪体对气流阻挡作用的大小，Cd越大，说明猪对气流的阻挡作用越大。空气阻力系数计算公式为：

Cd=2F/sv2

(1)

式中：F为物体受到的空气阻力，N；s为物体迎风面积，m2；v为物体与气流的相对速度，m/s。

用于模型对比验证的圆柱体空气阻力系数半经验公式由Hongwu Tang等[24]总结已有研究的试验数据并归纳得到式(2)：

(2)

式中：Re为雷诺数，表征流体流动情况。

使用空气阻力系数相对误差(Er)评判圆柱模拟结果的准确性，计算公式为：

(3)

式中：Cd,y为模拟所得的圆柱体空气阻力系数；Cd,b为圆柱体半经验公式计算得到的圆柱体空气阻力系数。

1.4 网格处理

1.4.1网格划分

OpenFOAM网格的生成主要由BlockMesh和SnappyHexMesh算法联合处理，主要有4个步骤(图2)：首先将STL(Stereolithography)格式的三角曲面的目标物生成为OpenFOAM可以识别的曲面网格(图2(a))；其次，使用BlockMesh算法创建由计算域边界填充整个区域的初始六面体背景网格(图2(b))；然后在初始六面体背景网格的基础上由SnappyHexMesh算法进行网格细化，网格细化可以提高网格精细度，以保证网格能准确计算流场信息，同时，为节约计算资源，网格细化的范围限制在围绕猪体展开的2.5 m×2.0 m×2.0 m的长方体空间区域(图2(c))；最后检查网格质量，优化网格,得到细化后的网格(图2(d))。

1.简化猪体模型；2.六面体背景网格；3.猪模型所在位置；4.网格细化区域；5.细化网格；6.不同级别网格之间的过渡网格；7.背景网格1.Simplified pig model; 2.Hexahedron background meshes; 3.Position of pig model; 4.Mesh refinement region; 5.Refined meshes; 6.Transition meshes between different levels of meshes; 7.Background meshes图2 OpenFOAM网格生成过程Fig.2 The process of grid generation in OpenFOAM

1.4.2网格数量

对于相同的模型，网格数量主要由网格尺寸决定。为了确定网格尺寸并检验网格独立性，本研究共设置12种网格级别(表1)，这12种网格级别具有相同的网格划分模式，仅由于背景网格尺寸不同而导致网格细化区域的网格尺寸和总网格数不同。

表1 不同网格级别的网格尺寸及网格总数Table 1 Mesh size of different mesh resolutions and the total number of meshes

1.5 边界条件

边界条件直接影响到气流流动的物理属性和模拟结果，本研究中虚拟风洞进风口、出风口和4个壁面与猪模型表面均为壁面边界(wall)，对猪模型近壁面区域和虚拟风洞地面区域使用壁面函数法进行处理(表2)。

表2 不同壁面的边界条件Table 2 Boundary conditions in different walls

1.6 计算模型与求解

由Menter[25]提出的k-ω SST两方程湍流模型适用于局部空气阻力与固体表面之间的气流模拟，并在相关研究中表现良好[26]，因此本研究使用k-ω SST模型作为模拟的湍流模型。OpenFOAM可以自由选择求解格式等参数，从而能够更加具有针对性地求解流场信息，但是容易造成求解的不稳定性，为了检验求解格式的稳定性，求解设置如下：数值模拟流体求解使用SIMPLE (Semi-implicit method for pressure-linked equation algorithms)算法，求解器使用OpenFOAM程序中的稳态求解器SimpleFoam，扩散项格式(LaplacianSchemes)设置为Guass linear corrected,表面差分格式(InterpolationSchemes)设置为Linear，梯度格式(GradSchemes)设置为Gauss linear,离散格式(DivSchemes)采用Bounded Guass upwind。

2 结果与讨论

2.1 网格划分对空气阻力系数的影响

为检验网格独立性，模拟得到v=1.0 m/s的流速下，12种网格处理模式圆柱体空气阻力系数(表3)。当背景网格尺寸为300～500 mm时，网格数量较少，空气阻力系数的相对误差和变化幅度均较大，模拟准确度和稳定性均较差；背景网格尺寸为230～280 mm时，空气阻力系数的相对误差和变化幅度均较小，模拟稳定性和准确性均较高，如本研究中第7网格级别的圆柱体模拟所得的空气阻力系数；背景网格尺寸继续减小到150～220 mm时，网格数量增加，空气阻力系数波动较大，说明在较小的网格尺寸下模拟的准确度和稳定性均较差。适当的网格尺寸保证模拟的准确度与稳定性，网格尺寸减小导致网格数量和计算成本的增加但不会提高计算精确度，Hong等[15]使用OpenFAOM模拟农业建筑自然通风时也得到了与本研究类似的结果。

表3 不同网格级别的圆柱体空气阻力系数模拟结果Table 3 Simulated results of drag coefficient ofcylinder with different mesh levels

2.2 模型验证与猪体模型对比

网格独立性检验完成后需要验证模型的准确性，分别在v=0.2、0.5、1.0、1.5、2.0和2.5 m/s流速下计算求解第7网格级别处理模式的圆柱体空气阻力系数(表4)。圆柱体空气阻力系数的相对误差波动较大,原因可能是：由于半经验公式是基于均匀来流的无限长直立圆柱体进行归纳，流体流动属于二维，但由于实际圆柱体三维流动的复杂性仍可能导致模拟结果的误差增加，尽管此次模拟的圆柱体长径比达到了使用公式所要求的4.0从而减小三维流动的影响，但是不能完全避免三维流动造成的干扰；另外一方面，OpenFOAM作为开源算法，在计算迭代过程中对结果的优化有限，特别是考虑到实际工程中湍流流动复杂，准确的测量结果的获得难度较大，所以该误差范围是合理正常的[27]。基于上述考虑，使用第7网格级别处理网格既保证网格独立性，又保证了湍流模型和OpenFOAM计算求解的可靠性，验证了模拟的准确性。

表4 不同气流流速下第7网格级别的圆柱体空气阻力系数模拟结果与相对误差Table 4 Simulated drag coefficient and the relativeerror of cylinder model with the 7th meshlevel under different air velocities

为检验圆柱体作为极度简化猪体模型的可行性，求解风洞中相同位置处的简化猪体模型和圆柱体在不同流速下的空气阻力系数(图3)。v为0.2～1.5 m/s时，猪的空气阻力系数Cd,p减小，v为1.5～2.5 m/s时，Cd,p增大，增值较缓(图3(a))。圆柱的空气阻力系数Cd,y随着气流流速的增大而减小。作为类圆柱体的简化猪体和圆柱体周围的气流流动属于亚临界区，空气阻力系数可以看作定值[24]，所以在v=0.5～2.5 m/s时取平均值得到Cd,p=0.943，Cd,y=1.207，二者大小差别明显。在变化规律方面，线性拟合圆柱和猪的空气阻力系数(图3(b))得到y=0.600 3x+0.223 5,R2=0.603 3，说明圆柱体模型与简化猪体模型阻力系数相关性不显著。几何形状的差异导致简化猪体模型(忽略猪体的四肢耳朵等非躯干部分的建模)相较于圆柱体更加能够表征单体猪体阻力系数变化规律。

2.3 气流流速和流向对猪体空气阻力系数的影响

不同流速和猪体与气流流向的不同夹角下猪体简化模型空气阻力系数模拟结果见图4。可见：随着猪体迎风面积的增大，空气阻力系数发生波动，但是存在增大趋势。θ=45°时，流速较高的情况下Cd,p陡增，这是因为作为类圆柱体的猪体模型在猪体与气流流向的夹角θ=45°时流场三维效应愈发明显[28-29]。相同流速下，θ=0°、15°、30°、150°、165°和180°时，猪体气流阻力较低；θ=45°、60°、120°和135°时，猪体阻力适中；θ=75°、90°和105°时，猪体阻力较大。由于猪舍生产过程中风速以不超过2.0 m/s 为宜，最佳风速在0.2 m/s左右，所以在实际生产中可以根据猪体与气流流向的夹角范围划分阻力区间范围从而评判猪体对气流的阻挡作用。

图4 猪体与气流流向夹角(θ)和气流流速(v)变化时猪体空气阻力系数(Cd,p)的模拟结果Fig.4 Simulated results of drag coefficient of pig (Cd,p) under different airflow angles (θ) and air velocities (v)

当气流流速v=0.2 m/s时，猪体相同迎风面积，头部迎风和尾部迎风的情况下Cd,p相差较小；v=1.0或2.5 m/s时，头部迎风和尾部迎风时Cd,p相差较大。以相同迎风面积，头部迎风θ=60°和尾部迎风θ=120°为例，v=0.2 m/s时，Cd,p分别为0.87 与0.86，压力云图差别较小(图5(a)、(b))；v=1.0 m/s时，Cd,p分别为0.73和0.85，压力云图差别较大(图5(c)、(d))。这是因为流速较低，表面摩擦阻力是影响阻力的关键；流速较高，形状阻力是影响阻力的关键，头部和尾部迎风，迎风形状不同则气流变化明显，空气阻力变化较大。

图5 距离风洞进风口2.4 m处不同猪体与气流流向夹角(θ)和不同气流流速(v)下的压力云图Fig.5 Pressure contour under different velocities (v) and airflow angles (θ) with 2.4 m from inlet of wind tunnel

3 结 论

本研究基于开源CFD(Computational fluid dynamics)算法OpenFOAM模拟风洞中单只猪对气流的阻挡作用，比较网格大小对计算结果稳定性和准确性的影响，分析了不同气流流速和猪体与气流流向不同的夹角下的猪体表面的空气阻力系数的变化规律，主要结论如下：

1)开源CFD软件OpenFOAM能够适用于畜禽舍空气阻力的研究，网格处理模式对计算精度和稳定性影响明显，合适的网格尺寸保证计算结果的稳定性和可靠性，减少网格数，节约计算成本。

2)流速较高时，猪体几何形状对气流流动影响较大，圆柱体作为极度简化的猪体模型不能够充分表征猪体的空气阻力系数的变化规律，二者的相关性较低。

3)气流流向对猪体空气阻力系数影响显著，可以根据猪体与气流流向的夹角范围划分阻力区间范围从而评判猪体对气流的阻挡作用的大小。