CN111611685A - 一种模拟抽水蓄能电站地下厂房轴流排风机工作流场的致动线方法 - Google Patents

Abstract

本发明公开了一种模拟抽水蓄能电站地下厂房轴流排风机工作流场的致动线方法，其特征是，包括如下步骤：1)建立抽水蓄能电站地下厂房排风机致动线模型；2)根据流场网格型心点与叶素点的距离关系找到离叶素点最近的网格单元，实现网格点自动识别；3)计算得到沿排风机叶片展向单位长度上的体积力，并按照笛卡尔坐标系对体积力进行分解；4)将体积力源项按照三维高斯分布的方式添加流场中去；5)计算模拟地下厂房流场。本发明可以解决因地下厂房轴流排风机实体建模而导致网格数量较多、计算时间较长的问题，适用于地下厂房排风时的流场计算，对地下厂房轴流排风机功率选择有一定的指导意义，在工程中有很好的应用前景。

Description

一种模拟抽水蓄能电站地下厂房轴流排风机工作流场的致动 线方法

技术领域

本发明属于抽水蓄能电站地下厂房排风机流场计算的技术领域，特别涉及一种模拟抽水蓄能电站地下厂房轴流排风机工作流场的致动线方法。

背景技术

近些年，大规模、大容量的抽水蓄能电站相继开工建设，抽水蓄能电站地下厂房洞室深度深、结构错综复杂，如何选取合适的大功率轴流排风机和与之对应的排风形式，使得厂房内产生的污风能够高效及时的排出，已成为人们关注的焦点。近年来，计算流体动力学方法在抽蓄电站大型排风机工作流场的研究当中得到越来越多的应用，该方法通过计算机仿真模拟的方式能够准确地描述出由抽蓄电站轴流排风机产生的排风通道复杂流场。常规的轴流排风机流场数值计算，为了保证叶片及通道流场的计算准确性，需要建立相对复杂的轴流排风机实物模型，即先将固体的实物模型用专业的三维建模软件进行建模(如SolidWorks)，建模完成后再对整个通道流场区域进行网格划分。然而，三维复杂实体模型的不规则形状和必要的边界层网格使得整个流场的网格划分成为数值计算的一大难点，而且在轴流排风机上下游较近的流场区域内网格仍要保持一定的精度，大大增加了计算域的网格数量，导致计算资源的浪费。

研究人员通过将BEM理论与常规的计算流体动力学方法相结合，开发出了致动线方法，即先用BEM理论求解轴流排风机叶片工作时产生的气动力，再将气动力作为体积力源项添加到控制方程中去，以模拟叶片对空气流的作用力。由于致动线方法中没有真实存在的叶片等的固体实物模型，因此无需考虑固壁边界的网格划分，这大大减少了流场网格划分的难度和复杂度，进而节约了计算资源。

目前致动线方法的主要研究方向是阻力型机械叶片的载荷特性及其产生的流场特性，尚未有将致动线方法应用于主动产生空气流的轴流排风机领域，尤其是抽水蓄能电站地下厂房的大功率长叶片轴流排风机，而合适的轴流排风机对于抽蓄电站排风效率、厂内热湿环境平衡具有重要意义。

发明内容

发明目的，为解决现有技术领域的缺失，本发明的目的在于提供一种模拟抽水蓄能电站地下厂房轴流排风机工作流场的致动线方法，运用致动线模型对抽蓄电站轴流排风机的排风通道进行数值计算，并与某抽蓄电站排风系统的实验模型测量结果进行对比，验证致动线模型运用在抽蓄电站轴流排风机工作流场数值仿真的可行性。

技术方案，本发明提出一种模拟抽水蓄能电站地下厂房轴流排风机工作流场的致动线方法，该方法包括如下步骤：

1)将轴流排风机叶片简化成一条线段，在该线段上施加不连续的体积力以模拟叶片对空气流的作用，形成致动线模型；

2)建立抽蓄电站地下厂房排风通道物理模型，并进行流场网格划分，包括轴流排风机的邻域网格划分(以下简称加密区)和其它流域网格划分，并且邻域网格划分比其它流域网格划分精细；

3)将轴流排风机叶片沿展长方向等分成若干个叶素，定义各叶素的中心点为叶素点，利用轴流排风机稳定工作时的旋转角速度ω和叶素点与排风机中心距离r确定在任一时刻下各叶素点的三维坐标，再根据流场加密区网格遍历的方法自动识别离叶素点最近的网格单元，将网格单元里存储的空气流速度信息赋给该叶素点；

4)通过叶素理论对致动线模型各叶素所在翼型计算其入流角

和迎角α大小，再根据二维翼型气动数据查表获得每个翼型的升力系数和阻力系数；计算得到叶片展长方向上各叶素点处空气流的体积力，包括轴向力和切向力，然后根据力的分解法则将轴向力和切向力按照笛卡尔坐标系进行三维分解；

5)将各叶素点处空气流计算得到的在笛卡尔坐标系下的三维体积力按照三维高斯分布的方式以源项的形式光顺地添加到流场中去；

6)计算模拟流场：叶片以某一角速度进行稳定旋转后，在每个时间步长上对流场加密区网格重复以上步骤，自动识别离每个叶素点最近的网格单元，并完成体积力的计算和光顺。

进一步地，所述步骤2)中建立抽蓄电站地下厂房流场计算物理模型，包括地下厂房、排风下平洞、排风竖井和排风机房流场模型。

在排风机房中，以轴流排风机叶轮平面为中心，将其所处的邻域单独进行网格加密该邻域为长方体，长、高一般为3倍的叶片半径，与风轮轴线方向平行的宽度一般大于等于6倍的光顺参数，加密区中的网格单元尺寸与叶素展长相等，而流场中其它区域采用渐变渐疏的蜂巢网格进行划分。

进一步地，所述步骤3)中对流场加密区网格遍历，找到离某叶素点最近的网格单元的判断方法是：当某个叶素点在某个网格单元内部时则认为该网格单元离叶素点最近，即网格单元型心到叶素点的x、y、z三个方向的距离均小于等于1/2倍的加密区网格单元对应边长。

进一步地，所述步骤3)中具体网格自动识别过程如下：

对于一个距离排风机中心为r，所在叶片的编号为i的叶素点j而言，其在当前瞬态时间步下的三维笛卡尔坐标是：

x＝x_wt

式中，x_wt、y_wt、z_wt为轴流排风机旋转平面中心的三维坐标，PHI0为叶片启动相位角，B为叶片数，θ为当前时间步下叶片已经旋转的角度，θ＝ωt，ω为轴流排风机旋转角速度，R为叶片半径，N为叶片分段数。

对加密区内所有的网格单元进行遍历，如果某网格单元型心坐标满足

式中，x、y、z为叶素点坐标；x_cell、y_cell、z_cell为网格单元型心坐标；DX、DY、DZ为加密区内网格单元对应边长，则认为该网格单元距离该叶素点最近，此时将该网格单元存储的空气流速度u_cell按速度梯度的方式赋给叶素点：

式中，u_ij为计算得到的第i个叶片上第j个叶素点处的空气流速度；

为该叶素点到网格单元型心处的流场速度梯度；dS为该叶素点到最近网格型心的距离。

进一步地，所述步骤4)中具体计算过程如下：

相对于某叶素点的空气流合速度

式中，U为轴向速度，

式中

为叶素点空气流三维速度矢量；

为轴流排风机旋转的法向量；V_rot为径向速度，

式中，

为叶素点展长向量，

为叶片旋转角速度向量；旋转平面的入流角

根据叶素点所在的扭角β得到叶素翼型迎角

叶片展长方向上单位长度上空气流受到的体积力(由轴向力和切向力组成)由式

和式

确定，式中，T为轴向力；F_N为切向力；ρ为抽蓄电站轴流排风机所在地的空气密度；V_rel为叶素翼型上空气流的合速度；c为叶素翼型弦长；

为入流角；C_l＝C_l(α,Re)为升力系数，C_d＝C_d(α,Re)为阻力系数，C_l、C_d都是以迎角α和雷诺数Re为变量的函数，升力系数和阻力系数均由叶素翼型的二维气动特性曲线得到；dr为叶素展长；

按照力的分解法则将叶素空气流受到的轴向力和切向力分解为笛卡尔坐标系下的三维体积力：

f_x＝T

式中，f_x、f_y、f_z分别代表叶素空气流所受到的x、y、z三个方向上的体积力；T为轴向力；F_N为切向力；PHI0为叶片启动相位角；B为叶片数；θ为当前时间步下叶片已经旋转的角度，θ＝ωt，ω为轴流排风机旋转角速度；i为叶片编号。

进一步地，所述步骤5)中体积力的分布方式为三维高斯分布，具体计算过程如下：

根据流场加密区中网格单元的尺寸确定体积力向流场光顺集中程度的参数ε，称为光顺参数，一般为2倍的加密区网格单元尺寸；对于一个叶片编号为i，叶素编号为j的某个叶素点而言，该叶素点处的空气流所受到的任一方向体积力f_ij按照式

的形式向流场进行光顺过渡，式中，f_ε为光顺之后的体积力；d为加密区内任一网格单元型心与叶素点的距离；

由于在计算域中，每个网格单元处的空气流均不同程度的受到叶素点处的空气流体积力的影响，因此以加密区内某一网格单元型心M点为例，它所受到的体积力应是轴流排风机叶片上所有叶素点处空气流体积力光顺到此处的总和，即：

式中，f_M为网格单元型心M受到的总体积力；B为风力机叶片数；N为每个叶片的分段数。

进一步地，所述步骤6)中在计算流体力学求解器中(如Fluent)设定边界条件(包括压力出入口和无滑移壁面)，将步骤1)、3)～5)所述的致动线模型通过udf编程导入Fluent软件中作为排风机的边界条件，然后在每个时间步上计算模拟整个地下厂房的流场情况。

有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：

本发明提出的一种模拟抽水蓄能电站地下厂房轴流排风机工作流场的致动线方法，通过网格自动识别技术，可以有效地减少实体建模和数值计算时间和流域网格数量，适用于抽蓄电站轴流排风机排风工作流场计算，可对抽蓄电站轴流排风机的选型和功率选择有一定的指导意义，在抽蓄工程中有很好的应用前景。

附图说明

图1是致动线模型流程图；

图2是实验模型计算域示意图，其中(a)为计算域示意图，(b)为加密区网格划分示意图，(c)为局部“蜂巢”网格划分示意图；

图3是致动线识别示意图；

图4是轴流排风机示意图；

图5是叶素二维翼型受力分析；

图6是经计算得到轴流排风机运行时叶片x方向体积力源项示意图；

图7是流场试验测点数据与模拟结果的速度对比图。

具体实施方式

如图1-7所示，本发明所建立的是一种模拟抽水蓄能电站地下厂房轴流排风机工作流场的致动线方法，建立的步骤如下：

步骤1)是将轴流排风机叶片简化成一条线段，在该线段上施加不连续的体积力以模拟叶片对空气流的作用，形成致动线模型。

步骤2)是建立抽蓄电站地下厂房流场计算物理模型(包括地下厂房、排风下平洞、排风竖井和排风机房等流场模型)，如图2(a)所示。

然后，在排风机房中，以轴流排风机叶轮平面为中心，将其所处的邻域单独进行网格加密(该邻域为长方体，长、高一般为3倍的叶片半径，与风轮轴线方向平行的宽度一般大于等于6倍的光顺参数)，如图2(b)所示，加密区中的网格单元尺寸与叶素展长相等，而流场中其他区域采用渐变渐疏的蜂巢网格进行划分，如图2(c)所示。

步骤3)是利用轴流排风机稳定工作时的旋转角速度ω和叶素点与排风机中心距离r确定在任一时刻下各叶素点三维坐标。对于一个与排风机中心距离为r，所在叶片的编号为i的叶素点j而言，其在当前瞬态时间步下的三维笛卡尔坐标是：

x＝x_wt

式中，x_wt、y_wt、z_wt为轴流排风机旋转平面中心的三维坐标；PHI0为叶片启动相位角；B为叶片数；θ为当前时间步下叶片已经旋转的角度，θ＝ωt，ω为轴流排风机旋转角速度；R为叶片半径；N为叶片分段数。

然后，如图3所示，对加密区内所有的网格单元进行遍历，如果某网格单元型心坐标满足

式中，x、y、z为叶素点坐标；x_cell、y_cell、z_cell为网格单元型心坐标；DX、DY、DZ为加密区内网格单元对应边长，则认为该网格单元距离该叶素点最近，此时将该网格单元存储的空气流速度u_cell按速度梯度的方式赋给叶素点：

式中，u_ij为第i个叶片上第j个叶素点处的空气流速度；

为该叶素点到网格单元型心处的流场速度梯度；dS为该叶素点到最近网格型心的距离。

步骤4)是由二维翼型数据根据叶素理论计算确定叶素翼型的气动力。如图4所示，某叶素翼型相对于空气流受到的轴向速度

式中

为叶素点空气流三维速度矢量；

为轴流排风机旋转的法向量；切向速度

式中，

为叶素点展长向量，

为叶片旋转角速度向量。

然后，如图5所示，叶素翼型相对于某叶素点的空气流合速度可由速度三角形得出

则旋转平面的入流角

进一步，根据叶素点所在的扭角β得到叶素翼型迎角

叶片展长方向上单位长度上空气流受到的体积力(由轴向力和切向力组成)由式

和式

确定，式中，T为轴向力；F_N为切向力；ρ为抽蓄电站轴流排风机所在地的空气密度；V_rel为叶素翼型上空气流的合速度；c为叶素翼型弦长；

为入流角；C_l＝C_l(α,Re)为升力系数；C_d＝C_d(α,Re)为阻力系数，两者都是以迎角α和雷诺数Re为变量的函数，升力系数和阻力系数均由叶素翼型的二维气动特性曲线得到；dr为叶素展长。

进一步，按照力的分解法则将叶素空气流受到的轴向力和切向力分解为笛卡尔坐标系下的三维体积力：

f_x＝T

式中，f_x、f_y、f_z分别代表叶素空气流所受到的x、y、z三个方向上的体积力；T为轴向力；F_N为切向力；PHI0为叶片启动相位角；B为叶片数；θ为当前时间步下叶片已经旋转的角度，θ＝ωt，ω为轴流排风机旋转角速度；i为该叶片编号。

步骤5)是叶素空气流受到的任一方向体积力按照式

的形式进行光顺过渡，式中，f_ij为叶片编号为i，叶素编号为j的某叶素点空气流所受到的任一方向体积力；f_ε为光顺之后的体积力；d为加密区内任一网格单元型心与叶素点的距离；ε为光顺参数，它与加密区网格单元尺寸相关，一般为2倍的加密区网格单元尺寸。

进一步，以加密区内某一网格单元型心M点为例，它所受到的体积力应是轴流排风机叶片上所有叶素点处空气流体积力光顺到此处的总和，即：

式中，f_M为网格单元型心M受到的总体积力；B为风力机叶片数；N为每个叶片的分段数。

步骤6)是将网格单元受到的体积力以动量源项的形式添加到流体控制方程中：

式中：

为速度；ρ为流体密度；p为压强；μ为粘度系数；

为网格单元所受到的体积力。计算模拟流场：叶片以某一角速度旋转后，

随着计算时间t和角速度ω变化，在每个时间步上对加密区网格重复以上步骤，自动识别离每个叶素点最近的网格单元，并完成体积力计算和光顺，源项添加的结果如图6所示。

进一步地，在计算流体力学求解器中(如Fluent)设定边界条件(包括压力出入口和无滑移壁面)，将步骤1)、3)～5)所述的致动线模型通过udf编程导入Fluent软件中作为排风机的边界条件，然后在每个时间步上计算模拟整个地下厂房的流场情况。

下面通过某抽蓄电站地上式排风机房布置方式的排风系统模型试验数据对所提出的模型进行验证。模型试验所采用的轴流排风机含有7根叶片，每根叶片长度均为0.05m，额定功率为20W，额定转速为5100RPM。以该试验型轴流排风机工作在额定转速工况下，用致动线方法对排风系统流场进行了数值模拟，并将计算结果与布置在排风管道上多个测点的试验测量值进行比较。整个计算域包括水平长矩形管道(代表排风下平洞)、垂直高度的圆柱形竖井以及地上排风机房模型，其中轴流排风机致动线模型位于排风机房出口前的加密区中，具体划分如图2。

采用非结构性网格，在排风机房轴流排风机邻域进行单独结构化网格加密，其余管道流域采用“蜂巢”网格形式，以减少网格数量，在每个轴流排风机叶片上分布网格数15个左右，以保证数值计算的准确性。

进一步的，如图7，在试验模型的排风下平洞和排风竖井各设置一个空气流速度监测点，每个监测点在横截面上按照中间、1/4近壁面和近壁面分别测量三次空气流速度，用以与致动线数值模拟的结果进行比较。从图7可以看出，两个测点测量得到的风速值与致动线模型数值模拟的结果较为接近，在截面的三个位置处速度值的吻合度均较高，且展现了管道流速的对称性，模拟得到的截面速度曲线走势与实验值基本一致，证明了所提出的轴流排风机致动线模型在抽蓄电站排风系统流场数值计算中的准确性。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (5)

1.一种模拟抽水蓄能电站地下厂房轴流排风机工作流场的致动线方法，其特征是，该方法包括如下步骤：

1)将轴流排风机叶片简化成一条线段，在该线段上施加不连续的体积力以模拟叶片对空气流的作用，形成致动线模型；

2)建立抽蓄电站地下厂房排风通道物理模型，并进行流场网格划分，包括轴流排风机的邻域网格划分和其它流域网格划分，并且邻域网格划分比其它流域网格划分精细，其中，邻域网格简称加密区；

3)将轴流排风机叶片沿展长方向等分成若干个叶素，定义各叶素的中心点为叶素点，利用轴流排风机稳定工作时的旋转角速度ω和叶素点与排风机中心距离r确定在任一时刻下各叶素点的三维坐标，再根据流场加密区网格遍历的方法自动识别离叶素点最近的网格单元，将网格单元里存储的空气流速度信息赋给该叶素点；

4)通过叶素理论对致动线模型各叶素所在翼型计算其入流角

和迎角α大小，再根据二维翼型气动数据获得每个翼型的升力系数和阻力系数；计算得到叶片展长方向上各叶素点处空气流的体积力，包括轴向力和切向力，然后根据力的分解法则将轴向力和切向力按照笛卡尔坐标系进行三维分解；

5)将各叶素点处空气流计算得到的在笛卡尔坐标系下的三维体积力按照三维高斯分布的方式以源项的形式光顺地添加到流场中去；

6)计算模拟流场：叶片以某一角速度进行稳定旋转后，在每个时间步长上对流场加密区网格重复以上步骤，自动识别离每个叶素点最近的网格单元，并完成体积力的计算和光顺。

2.根据权利要求1所述的一种模拟抽水蓄能电站地下厂房轴流排风机工作流场的致动线方法，其特征是，所述步骤2)中建立抽蓄电站地下厂房流场计算物理模型包括地下厂房、排风下平洞、排风竖井和排风机房流场模型。

3.根据权利要求1所述的一种模拟抽水蓄能电站地下厂房轴流排风机工作流场的致动线方法，其特征是，所述步骤3)中对流场加密区网格遍历，找到离某叶素点最近的网格单元的判断方法是：当某个叶素点在某个网格单元内部时则认为该网格单元离叶素点最近，即网格单元型心到叶素点的x、y、z三个方向的距离均小于等于1/2倍的加密区网格单元对应边长，具体网格自动识别过程如下：对于一个距离排风机中心为r，所在叶片的编号为i的叶素点j而言，其在当前瞬态时间步下的三维笛卡尔坐标是：

x＝x_wt

式中，x_wt、y_wt、z_wt为轴流排风机旋转平面中心的三维坐标，PHI0为叶片启动相位角，B为叶片数，θ为当前时间步下叶片已经旋转的角度，θ＝ωt，ω为轴流排风机旋转角速度，R为叶片半径，N为叶片分段数；

对加密区内所有的网格单元进行遍历，如果某网格单元型心坐标满足

式中，x、y、z为叶素点坐标；x_cell、y_cell、z_cell为网格单元型心坐标；DX、DY、DZ为加密区内网格单元对应边长，则认为该网格单元距离该叶素点最近，此时将该网格单元存储的空气流速度u_cell按速度梯度的方式赋给叶素点：

式中，u_ij为计算得到的第i个叶片上第j个叶素点处的空气流速度；

为该叶素点到网格单元型心处的流场速度梯度；dS为该叶素点到最近网格型心的距离。

4.根据权利要求1所述的一种模拟抽水蓄能电站地下厂房轴流排风机工作流场的致动线方法，其特征是，所述步骤4)中具体计算过程如下：相对于某叶素点的空气流合速度

式中，U为轴向速度，

式中，

为叶素点空气流三维速度矢量，u,v,w分别表示x，y，z方向的速度标量；

为轴流排风机旋转的法向量；V_rot为径向速度，

式中，

为叶素点展长向量，

为叶片旋转角速度向量，旋转平面的入流角

根据叶素点所在的扭角β得到叶素翼型迎角

叶片展长方向上单位长度上空气流受到的体积力，其由轴向力和切向力组成，由式

和式

确定，式中，T为轴向力；F_N为切向力；ρ为抽蓄电站轴流排风机所在地的空气密度；V_rel为叶素翼型上空气流的合速度；c为叶素翼型弦长；

为入流角；C_l＝C_l(α,Re)为升力系数，C_d＝C_d(α,Re)为阻力系数，C_l、C_d都是以迎角α和雷诺数Re为变量的函数，升力系数和阻力系数均由叶素翼型的二维气动特性曲线得到，dr为叶素展长；

按照力的分解法则将叶素空气流受到的轴向力和切向力分解为笛卡尔坐标系下的三维体积力：

f_x＝T

式中，f_x、f_y、f_z分别代表叶素空气流所受到的x、y、z三个方向上的体积力；T为轴向力；F_N为切向力；PHI0为叶片启动相位角；B为叶片数；θ为当前时间步下叶片已经旋转的角度，θ＝ωt，ω为轴流排风机旋转角速度；i为叶片编号。

5.根据权利要求1所述的一种模拟抽水蓄能电站地下厂房轴流排风机工作流场的致动线方法，其特征是，所述步骤5)中体积力的分布方式为三维高斯分布，具体计算过程如下：根据流场加密区中网格单元的尺寸确定体积力向流场光顺集中程度的参数ε，称为光顺参数，对于一个叶片编号为i，叶素编号为j的某个叶素点而言，该叶素点处的空气流所受到的任一方向体积力f_ij按照式

的形式向流场进行光顺过渡，式中，f_ε为光顺之后的体积力，d为加密区内任一网格单元型心与叶素点的距离；

加密区内某一网格单元型心M点所受到的体积力应是轴流排风机叶片上所有叶素点处空气流体积力光顺到此处的总和，即：

式中，f_M为网格单元型心M受到的总体积力；B为风力机叶片数；N为每个叶片的分段数。

Images