光热发电吸热塔对飞机颠簸影响的数值分析

Numerical Analysis of the Effect of Photothermal Power Generation Absorber Tower on Aircraft Bumpiness

摘要: 光热发电吸热塔在运行条件下对周围空气形成加热，易导致吸热塔周围出现热对流现象。热对流对临近区域飞行运行安全可能造成一定影响。本文基于计算流体方法建立了光热发电吸热塔三维计算模型，分析了吸热塔顶部不同高度的温度场、湍流强度以及涡流耗散率的大小，并基于涡流耗散率的大小评估了吸热塔对飞机颠簸的影响。结果显示，在当前条件下吸热塔所形成的局部热湍流影响的范围较小，飞机在距离吸热塔高50 m以上及吸热塔后方600 m范围以外运行基本没有影响。

Abstract: In solar thermal power generation systems, the heat absorption tower heats the surrounding air during operation, which can induce thermal convection in its vicinity. This thermal convection may pose potential risks to the flight safety of aircraft operating in adjacent airspace. To investigate this phenomenon, a three-dimensional computational model of the heat absorption tower was developed using computational fluid dynamics (CFD). The model was employed to analyze the temperature field, turbulence intensity, and eddy dissipation rate at various heights near the top of the tower. Based on the eddy dissipation rate, the impact of the heat absorption tower on aircraft turbulence was evaluated. The results indicate that, under the given conditions, the spatial extent of local thermal turbulence generated by the heat absorption tower is limited. Specifically, the influence on aircraft operations becomes negligible at altitudes exceeding 50 meters above the tower and at distances greater than 600 meters downstream of the tower.

1. 引言

光热发电吸热塔对飞行的影响主要是两个方面，一个方面是吸热塔的高度较高，此高度对航路运行可能会作为障碍物被考虑来分析其对飞行的影响；另一个方面，吸热塔中的吸热器表面聚集的大量热量会通过对流换热和辐射换热的方式释放到周围环境中，使得吸热塔附近区域的空气密度和温度形成明显的差异，从而导致热对流的形成，即对流型湍流。本文关注由热对流而形成的对流型湍流对飞行的影响。热气流通常不会对飞机造成严重问题、但在飞行过程可能很艰难。轻型的飞机可能会经历颠簸引起乘客的不适，而大型飞机由于质量较大，惯性较大，反而可以抵消部分额外的颠簸[1]。湍流是大气边界层中最重要的物理过程之一，湍流的出现通常会引起飞机的颠簸，颠簸现场可能会导致乘客的不适感出现，提高对湍流的理解可以更好地为飞行运行提供指导以及为飞行安全提供保障[2]。

普遍认为有了较强的湍流就必然有较强的飞机颠簸。这种观点导致经常出现预报结果与事实不符的现象[3]。事实上有了较强的湍流，有时有较剧烈的飞机颠簸，而有时却是轻微的颠簸，甚至只是略有振颤。另一方面，在湍流不强时，有时却有一定强度的颠簸出现。湍流强度分为稳定飞行、弱湍流、中度湍流以及强湍流四种等级[4]。为了呈现飞机颠簸的程度，可以通过计算湍流强度和涡流消散率EDR (Eddy Dissipation Rate, m^2/3s⁻¹)等效表达。

目前，大量研究使用EDR来评价由环境大气湍流引起的飞机涡流。世界气象组织从2002年开始在飞机气象数据中加入湍流监测信息，其中EDR和等效垂直阵风速度都作为测量指数报告[5]。EDR表示扰动能量从大尺度漩涡向小尺度漩涡的传递率，EDR越大，湍流就越强。当前关于大气湍流的研究主要基于雷达或者快速存储器记录的数据来分析。Huang [6]基于快速存储器记录数据分析了航路中可能出现的大气湍流对飞行的影响，并提出了新的计算方法。Sarpkaya [7]提出了新的模型去预测涡流耗散率的大小，通过与测量的数据比较显示了模型的可靠性。Sharman [8]执行了一系列的试验研究了非对流型湍流的能量耗散流动状态。另外，Chan [9]以及Takacs [10]等人基于雷达测量的数据分析了机场周围的湍流流动情况，分析了不同类型湍流对航空器运行的影响。

本文结合光热吸热器的运行条件，通过数值计算的方法来评估光热吸热塔周围形成热湍流的状态[11]，并评估其对飞行安全的影响。根据吸热器的额定功率以及吸热效率可以估算出吸热器的热量损失，综合考虑多方面的因素本项目热损失功率约为30 MW。本文通过计算吸热塔周围的湍流强度以及涡流耗散率来评估热湍流对飞机颠簸的影响。

2. 数值方法及计算模型

2.1. 数学模型

任何流体的运动都必须满足质量守恒、动量守恒、能量守恒，对于牛顿流体可以使用纳维–斯托克斯方程(N-S方程)来描述：

连续性方程

$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho v\right)=0$ (1)

动量方程

$\frac{\partial }{\partial t}{\displaystyle \underset{{\Omega }_{CV}}{\int }\rho vd\Omega }+{\displaystyle \underset{S_{CV}}{\int }\rho vv\cdot ndS}={\displaystyle \underset{S_{CV}}{\int }T\cdot ndS}+{\displaystyle \underset{{\Omega }_{CV}}{\int }\rho bd\Omega }$ (2)

能量方程

$\frac{\partial }{\partial t}{\displaystyle \underset{{\Omega }_{CV}}{\int }\rho hd\Omega }+{\displaystyle \underset{S_{CV}}{\int }\rho hv\cdot ndS}={\displaystyle \underset{S_{CV}}{\int }k\nabla T\cdot ndS}+{\displaystyle \underset{{\Omega }_{CV}}{\int }\left(v\cdot \nabla p+S:\nabla v\right)d\Omega }+\frac{\partial }{\partial t}{\displaystyle \underset{{\Omega }_{CV}}{\int }pd\Omega }$ (3)

湍流方程

$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho k{u}_i\right)=\frac{\partial }{\partial x_j}\left({\Gamma }_k\frac{\partial k}{\partial x_j}\right)+\overline{G_k}-Y_k+S_k$ (4)

$\frac{\partial }{\partial t}\left(\rho \omega \right)+\frac{\partial }{\partial x_j}\left(\rho \omega u_j\right)=\frac{\partial }{\partial x_j}\left({\Gamma }_{\omega }\frac{\partial \omega }{\partial x_j}\right)+G_{\omega }-Y_{\omega }+D_{\omega }+S_{\omega }$ (5)

其中，$\overline{G_{\text{k}}}$项表示由于平均速度梯度产生的湍流动能，由$G_{\text{k}}$ 项计算得到。$G_{\omega }$ 是由ω方程产生的。${\Gamma }_k$ 和${\Gamma }_{\omega }$ 分别表示k方程和ω方程的有效扩散系数。$Y_{\text{k}}$ 和$Y_{\omega }$ 分别表示由于湍流产生的k方程和ω方程的耗散率。$D_{\omega }$ 表示交叉扩散项。$S_{\text{k}}$ 和$S_{\omega }$ 表示自定义源项。

辐射模型

模拟选取DO模型来计算辐射带来的影响。对于吸收、发射以及散射性介质的$\overrightarrow{r}$ 位置$\overrightarrow{s}$ 方向的辐射传输方程为：

$\frac{dI\left(\overrightarrow{r},\overrightarrow{s}\right)}{ds}+\left(\alpha +{\sigma }_s\right)I\left(\overrightarrow{r},\overrightarrow{s}\right)=\alpha n^2\frac{\sigma T^4}{\pi }+\frac{{\sigma }_s}{4\pi }{\displaystyle {\int }_0^{4\pi }I\left(\overrightarrow{r},{\overrightarrow{s}}^{/}\right)\Phi \left(\overrightarrow{s},{\overrightarrow{s}}^{/}\right)}d{\Omega }^{/}$ (6)

2.2. 几何模型及计算网格

如图1(a)所示给出了吸热塔的热物理特性计算简化模型，主要由吸热器以及吸热器上下两端的连接部分组成，吸热器的高度为14 m，外径为15 m，吸热器距离地面高度为210 m。考虑到除吸热器以外的其他区域热量损失相比吸热器而言较小，因此在分析吸热塔的热量传递时，仅仅考虑吸热器的热量交换。

如图1(b)给出了流场计算所需要的计算域，由于整个外界环境有微风吹过，并且吸热塔的部分热量会散发到空气中，从而导致吸热塔的周围出现热对流现象。综合考虑湍流尺度效应以及定日镜的安装范围，选取如图1(b)所示的计算域，长、宽、高为1600 m × 1600 m × 1422 m。

在前处理中，需要根据几何模型的结构选择合适的网格划分软件，本研究网格划分利用网格生成工具Fluent Meshing。Fluent Meshing能够处理复杂结构及区域的网格，可以建立高质量非结构化网格。吸热塔流体域网格为Ployhedra体网格，如图2所示，其体网格总数量为800 W左右，根据文献显示，同样网格密度对计算结果无较大影响。

(a) (b)

Figure 1. Geometry model of absorber tower

图1. 吸热塔计算几何模型

Figure 2. Absorbing tower computational grid

图2. 吸热塔计算网格

2.3. 边界条件及离散方法

根据此光热场址典型年气象数据计算得到年平均风速为1.67 m/s，温度为8.5℃。入口边界条件为速度入口，其中湍流耗散率以及湍流粘性比分别设置为5和10。出口边界条件选定为自由流出口，定义所有出口压力为101,325 Pa。壁面为无滑移边界条件，壁面温度通过耦合方式实现。空气密度选取不可压缩理想气体模型，吸热器的热流密度为45,472.84 w/m²。

由于地表摩擦的结果，使接近地表的风速随着离地面高度的增加而增大，在离地面300~500 m以上的地方，风速将不受地表的影响，能够在气压的作用下自由流动，从而达到最大速度。平均风速沿高度的变化规律，常称为平均风速梯度，出现这种速度的高度叫梯度风高度。在400 m以上，认为风速保持不变。在400 m以下，风速按指数关系变化

$v=v_0{\left(z/z_0\right)}^{\alpha }$ (7)

其中，$v$ 为某一高度上的风速，$v_0$ 为地面上方10 m高度测得的风速，即为气象局测试的风速。$z$ 为高度，$z_0$ 为10 m，$\alpha$ 为地面粗糙系数，该值与地形地貌、环境等因素有关，对于平阔平原其值为0.167。根据提供的数据，地面10 m处平均风速为1.67 m/s。对于越平阔的地面状态，$\alpha$ 越小。外部大空间，在上风口设置为速度进口，进行UDF编程设置；其他出口面设置为，压力出口，给定为环境温度。

在数值计算的流动控制方程求解过程中，SIMPLE算法用于求解不可压缩流场，数值计算方法使用有限体积法对流体区域进行离散。控制方程的扩散项采用中心差分格式离散，而对流项则可采用二阶迎风格式。亚松驰因子分别为：密度0.8，动量0.8，能量1。

3. 结果分析与讨论

3.1. 温度场

(a) 吸热塔表面温度 (b) 吸热塔表面辐射温度

(c) X = 0 m温度剖面图

(d) X = 0 m辐射温度剖面图

(e) Z = 237 m剖面图 (f) Z = 257 m剖面图

(g) Z = 277 m剖面图 (h) Z = 297 m剖面图

Figure 3. Temperature distribution at different positions

图3. 不同位置的温度分布

在进行温度场计算时，吸热塔中的吸热器被考虑为30 MW的热源。通过离散能量方程得到吸热器表面的温度分布如图3(a)和图3(b)，吸热器表面的温度约为673℃，观察吸热塔表面辐射温度分布图可知吸热器的热量主要通过辐射的方式释放到环境中，由于空气的对流换热系数较小使得对流换热并不占主导地位。进一步地，结合图3(c)和图3(d)可以分析辐射换热的规律，通过给出的辐射温度分布规律，吸热器的温度以放射状的形态传递到周围环境中，离吸热器越远受辐射的影响较小，其影响的范围在50 m以内。

图3(e)~(h)给出了不同高度处的温度分布。从图中可以看出随着高度的增加温度逐渐降低。吸热塔顶部(237 m)相较于257 m位置，受到速度梯度的影响较大导致在吸热塔10 m范围内的温度分布规律出现比较明显的差异。一旦远离吸热塔，流动的不均匀性越小，温度的分布差异更小。在297 m位置最大温度约为9.3℃，相较于远场的温度8.5℃差异不大，表明该位置及以上高度接受的热量很小。

3.2. 湍流强度

大气湍流强度是影响飞机颠簸的重要参数。湍流强度为湍流脉动速度与平均速度的比值，是衡量湍流强弱的相对指标。湍流强度描述了风速随时间和空间变化的程度，反映脉动风速的相对强度，是描述大气湍流运动特性最重要的特征量。一般认为，I < 1%为低强度湍流，I > 10%为高强度湍流。湍流强度定义为速度脉动的均方根$\widehat{u}$ 与平均速度${\mu }_{avg}$ 的比值。其计算方式为：

$I=\frac{\widehat{u}}{u_{avg}}$ (8)

上式中，速度脉动的均方根$\widehat{u}$ 为：

$\widehat{u}=\sqrt{\left(u^2+v^2+w^2\right)/3}\approx u$ (9)

其中，μ，ν，w指速度脉动量，但因为各向同性假设，u = v = w。

(a) Z = 237 m剖面图 (b) Z = 257 m剖面图

(c) Z = 277 m剖面图 (d) Z = 297 m剖面图

Figure 4. Turbulence intensity distribution at different locations

图4. 不同位置的湍流强度分布

图4给出了吸热塔上方不同位置湍流强度分布。图中可以看出，在吸热塔附近，紊流强度较大，其次是吸热塔后下风方向，气流紊流流动比较厉害，紊流强度较高。吸热塔上方10 m范围内出现高强度湍流，在10 m~60 m范围内为中度湍流强度。湍流的影响范围主要在吸热塔附近10 m，影响距离在吸热塔后方600 m范围。可能的原因在于一定速度的自然风流过吸热塔时，易在吸热塔周围形成扰流现象，并在吸热塔后方形成涡流区。涡流区内具有较大的动能，使得气流的湍动能加剧，湍流强度增大。

3.3. EDR分析

EDR值被国际民航组织认定为评价飞行中大气湍流对飞机影响强度的指标，该指标是一种独立于飞机本身的客观指标，它测量飞机周围大气状况，且不受飞机的类型或大小影响。因此，在分析吸热塔附近热对流对飞行影响时，主要参考EDR值进行评价。EDR值通过表达式EDR = ε1/3计算而来。湍流耗散率可以通过计算流场速度的梯度张量来获得。梯度张量可以用速度场的测量数据计算，通过矩阵运算来确定每个点的速度梯度张量，然后通过计算张量的迹来确定湍流耗散率：

$\epsilon =C_{\mu }^{3/4}\frac{k^{3/2}}{l}$ (10)

其中，C_u为经验常数，默认值为0.09。k为湍动能，可以通过湍流强度及平均速度进行估算：

$k=\frac{3}{2}{\left(u_{avg}I\right)}^2$ (11)

(a) Z =237 m剖面图 (b) Z = 257 m剖面图

(c) Z = 277 m剖面图 (d) Z = 297 m剖面图

Figure 5. Distribution of eddy dissipation rate at different positions

图5. 不同位置的涡流消散率(EDR)分布

根据EDR值与湍流强度之间的关系，当EDR小于0.1时认为对飞行基本没有影响，EDR在0.1~0.3之间时认为出现轻度湍流，EDR在0.3~0.5之间时认为出现中度湍流，EDR在0.5以上时认为出现重度湍流。图5为吸热塔上方不同高度处的EDR大小分布情况，可以看出随着高度增加，EDR值逐渐减小，湍流强度也逐渐降低，对飞行运行影响也更弱。在297 m高度处，最大的EDR为0.08 (小于0.1)，在该高度以上出现的湍流对飞行基本没有影响，因此本项目在分析湍流影响时最大高度就选取在297 m处。在297 m以下，根据EDR值可以看出吸热塔上方15 m范围内出现了中度湍流，甚至在吸热塔附近出现了中度湍流。在吸热塔上方15 m~50 m范围内，出现了轻度湍流。不管是中度湍流还是轻度湍流，对飞行安全都会造成一定程度的影响。

4. 结论

利用计算流体动力学方法，根据环境参数建立了储热型光热项目的吸热塔三维CFD模型，并对典型气象条件下的由于吸热塔而形成的热湍流过程进行数值计算，获得吸热塔及吸热塔附近空气的运动特征，温度场、湍流强度以及EDR分布等，为主管和决策部门分析、评价吸热塔附近的热流动对机场飞机运行的影响提供基础数据。

数值模拟结果表明，由于吸热塔中吸热器的大量热量散发到大气环境中，导致在吸热塔周围空气在热泳力的作用下形成热湍流。并且，由于受到环境风速的影响，吸热塔周围向上流动的空气流经吸热塔顶部之后，迅速向流动的下游扩散。在257 m高度处，空气基本不向上方飘动，下游最大影响范围在600 m，影响半径15 m。

目前用于定量分析热湍流对飞机运行安全影响的评价指标主要是EDR值，但是更加权威和官方的评价体系或标准尚未形成。EDR指标是一种独立于飞机本身的客观指标，它测量飞机周围大气状况，且不受飞机的类型或大小影响。EDR的计算可通过计算湍流耗散率实现，本研究通过植入相应数学公式在程序中完成分析。根据数值模拟数据显示，吸热塔所形成的局部热湍流影响的范围较小。在吸热塔50 m高度以上，最大的EDR为0.08，在该高度以上出现的湍流对飞行基本没有影响。在吸热塔上方15 m至50 m以下高度范围内出现了轻度湍流。而在吸热塔上方15 m高度以下，出现了中度湍流，吸热塔附近甚至出现了重度湍流。热湍流在吸热塔的后方影响距离约为600 m。因此，只要飞机只要在距离吸热塔高50 m、吸热塔后方600 m范围以外运行，基本没有影响。