CN101576419A - 由圆管外壁温度计算内壁温度的方法 - Google Patents

Abstract

本发明涉及测量技术领域，本发明提供了一种对于二维非稳态管内流动问题由圆管外壁温度计算内壁温度的计算方法，所述方法包括：在二维非稳态管内流动下，获得多个外壁节点的外壁温度测量值；根据所述外壁温度测量值得到多个内壁节点的内壁温度结果值。本发明可以准确的由管外壁温度计算得到管内壁温度随时间的变化，计算精度高，计算量较小，较好的解决了工业生产、科研工作中对于管内流动管内壁温度测量困难的问题，具有广泛的应用前景。

Description

由圆管外壁温度计算内壁温度的方法

技术领域

本发明涉及测量技术领域，特别涉及二维非稳态管内流动下的温度测量。

背景技术

在工业生产、科研工作中许多涉及管内流动问题的场合需要知道管壁温度随时间的变化，例如，研究热力管道中热应力对管壁材料结构的影响时。对于管外壁温度，已经有许多成熟的实验方法可以对其进行测量；而对于管内壁温度，由于管内壁测温装置的安装非常困难，而且测温装置的安装不当还会扰动管内流场，使管内壁温度的测量值偏离实际值。因此，如果能够根据容易测量的管外壁温度通过一定方法计算得到管内壁温度，即进行导热反问题的求解，则避免了测量管内壁温度给生产、科研带来的阻碍和困难。目前，国内外科研人员已做了很多导热反问题方面的研究工作，发展出了最大熵法、共轭梯度法、格林函数法、神经网络法等多种方法。对一维的稳态及非稳态问题，已有成熟的方法能够由管外壁温度求解管内壁温度；对二维、稳态管内流动问题，也有成熟的由管外壁温度求解管内壁温度反问题求解方法；但对于二维、非稳态管内流动问题，由于管内壁温度和管外壁温度随时间和空间变化，导热反问题的求解比较复杂，相关研究工作较少，因此，目前还没有用于二维、非稳态问题由管外壁温度求解管内壁温度导热反问题方法。

发明内容

本发明的目的是针对现有技术存在的上述不足，提供一种管内流动问题圆管内壁温度的计算方法，所述方法可根据由实验测量得到的圆管外壁若干点的温度随时间变化的函数，通过计算得到管内壁若干点的温度随时间变化的函数。

为实现上述目的，本发明采用如下技术方案：

在二维非稳态管内流动下，获得多个外壁节点的外壁温度测量值；

根据所述外壁温度测量值得到多个内壁节点的内壁温度结果值。

进一步地，所述根据所述外壁温度测量值得到多个内壁节点的内壁温度结果值的步骤包括：

设置内壁温度初始值；

根据所述内壁温度初始值计算得到所述多个外壁节点的外壁温度计算值；

判断所述外壁温度计算值与所述外壁温度测量值是否相符，若相符，则所述内壁温度初始值为内壁温度结果值；若不相符，通过迭代计算获得结果值。

进一步地，所述多个外壁和内壁节点由网格划分方法获得；

进一步地，所述多个外壁节点的外壁温度计算值通过求解圆管二维非稳态导热方程获得；

进一步地，所述判断所述外壁温度计算值与所述外壁温度测量值是否相符的步骤包括：

根据所述外壁温度计算值与所述外壁温度测量值得到当前的目标函数的值；

判断所述当前的目标函数的值是否小于一个预先确定的阈值，若是，则相符；若否，则不相符；

进一步地，所述目标函数为： $J=\frac{1}{2}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{[T_{\mathrm{Cal},m,n}-T_{\mathrm{Exp},m,n}]}^2,$ 其中，T_cal，m，n是外壁温度计算值，T_Exp，m，n是外壁温度测量值，M、N分别为测量值节点位置数和时间点数；

进一步地，所述通过迭代计算获得结果值的步骤之前还包括：根据迭代方程获得下一步迭代的内壁温度初始值；所述迭代方程为： $T_{k,n}^{b+1}=T_{k,n}^b-{\mathrm{\β}}^b{d}_{k,n}^b$ 其中，b为迭代次数，k为管内壁待计算节点的位置，β为迭代步长，n为时间点，d为搜索方向；

进一步地，所述迭代步长与表示内外壁节点温度变化关系的灵敏度系数有关，所述灵敏度系数由灵敏度方程获得，其中，所述灵敏度方程在圆管内壁的边界条件为： $\frac{\∂T}{{\∂T}_{k,n}}=l(s,\mathrm{\τ}),$ 并有， $\begin{array}{cc}\mathrm{if\; s}=k,\mathrm{\τ}=n.& l(s,\mathrm{\τ})=1\\ \mathrm{otherwise}& l(s,\mathrm{\τ})=0\end{array};$ 其中，

为灵敏度系数。

附图说明

图1是本发明实施例的二维左右对称管内流动问题示意图；

图2是本发明实施例的0°、60°、90°、120°、180°5个外壁测量节点的位置图；

图3是本发明实施例的0°、60°、90°、120°、180°5个内壁测量节点的位置图；

图4是通过实验得到的本发明实施例的0°、60°、90°、120°、180°5个内壁节点处的内壁温度准实验值(温度随时间的变化关系)；

图5是通过实验得到的本发明实施例的0°、60°、90°、120°、180°5个外壁节点处的外壁温度准实验值(温度随时间的变化关系)；

图6是本发明实施例的网格划分示意图；

图7是本发明实施例的内壁0°节点处的内壁温度计算值(温度随时间的变化关系)与内壁0°节点处的温度准实验值(温度随时间的变化关系)的对比；

图8是图7中60秒至80秒时间段的放大图；

图9是本发明实施例的内壁60°节点处的内壁温度计算值(温度随时间的变化关系)与内壁60°节点处的内壁温度准实验值(温度随时间的变化关系)的对比；

图10是本发明实施例的内壁120°节点处的内壁温度计算值(温度随时间的变化关系)与内壁120°节点处的内壁温度准实验值(温度随时间的变化关系)的对比。

具体实施方式

本发明提出的圆管内壁温度的计算方法，结合附图和实施例说明如下。

如图1所示的二维左右对称圆管内流动问题，其中，20为管内壁，30为管外壁，A为管内侧流体，B为管外侧流体；

已知圆管尺寸、管外流体温度T_sir、管外侧对流换热系数h_out。

对于一组管内壁温度值T^*＝T^*(r，φ，t)，可以由二维柱坐标非稳态导热方程计算得到所述圆管外壁节点温度随时间的变化，二维柱坐标非稳态导热方程为：

$\frac{\∂T}{\∂t}=\mathrm{\α}(\frac{1}{r}\frac{\∂}{\∂r}\left(r\frac{\∂T}{\∂r}\right)+\frac{1}{r^2}\frac{{\∂}^{2}T}{\∂{\mathrm{\φ}}^2})---\left(1\right)$

式中，t为时间，r、φ为节点在柱坐标系中的坐标。

圆管外侧温度的计算为第三类边界条件下的计算，即给定对流换热系数和环境温度：

$-\mathrm{\λ}{\left(\frac{\∂T}{\∂n}\right)}_w=h_{\mathrm{out}}(T-T_{\mathrm{sir}})---\left(2\right)$

圆管内侧温度的计算为第一类边界条件，即给定节点温度值：

T＝T^*    (3)

管内壁温度由迭代的初始条件通过基于共轭梯度法的导热反问题方法通过迭代计算求解逼近真实值。

迭代的初始条件为：

T(r，φ，0)＝T₀    (4)

共轭梯度法的目标函数为：

$J=\frac{1}{2}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{[T_{\mathrm{Cal},m,n}-T_{\mathrm{Exp},m,n}]}^2---\left(5\right)$

式中，T_cal，m，n是管外壁温度计算值，T_Exp，m，n是管外壁温度实验值，下标M、N分别为测点位置数和时间点数。管内壁温度迭代式为：

$T_{k,n}^{b+1}=T_{k,n}^b-{\mathrm{\β}}^b{d}_{k,n}^b---\left(6\right)$

式中，上标b表示迭代次数，下标k表示管内壁待求点的位置，下标n表示时间点。d为搜索方向，其计算式为：

${d^b}_{k,n}={\left(\frac{\∂J}{{\∂T}_{k,n}}\right)}^b+{\mathrm{\γ}}^b{d^{b-1}}_{k,n}---\left(7\right)$

式中γ由下式计算

${\mathrm{\γ}}^b=\frac{{\mathrm{\Σ}}_{k=1}^K{\mathrm{\Σ}}_{n=1}^N\left[{\left(\frac{\∂J}{{\∂T}_{k,n}}\right)}^b\right]}{{\mathrm{\Σ}}_{k=1}^K{\mathrm{\Σ}}_{n=1}^N{\left[{\left(\frac{\∂J}{{\∂T}_{k,n}}\right)}^{b-1}\right]}^2}---\left(8\right)$

迭代步长β计算式为

${\mathrm{\β}}^b=\frac{{\mathrm{\Σ}}_{j=1}^M{\mathrm{\Σ}}_{i=1}^N(T_{i,j,\mathrm{cal}}-T_{i,j,\mathrm{mea}}){\mathrm{\Σ}}_{k=1}^K{\mathrm{\Σ}}_{n=1}^N{\left(\frac{{\∂T}_{i,j}}{\∂T_{k,n}}\right)}^b{d_{k,n}}^b}{{\mathrm{\Σ}}_{j=1}^M{\mathrm{\Σ}}_{i=1}^N{\left[{\mathrm{\Σ}}_{k=1}^K{\mathrm{\Σ}}_{n=1}^N{\left(\frac{{\∂T}_{i,j}}{{\∂T}_{k,n}}\right)}^b{d_{k,n}}^b\right]}^2}---\left(9\right)$

其中，为表示内外壁节点的温度变化关系的灵敏度系数，为得到该灵敏度系数，需要求解灵敏度方程。灵敏度方程由非稳态导热方程对T_k，n求微分得到：

$\frac{\∂}{\∂t}\left(\frac{\∂T}{{\∂T}_{k,n}}\right)=\mathrm{\α}[\frac{{\∂}^2}{{\∂x}^2}\left(\frac{\∂T}{{\∂T}_{k,n}}\right)+\frac{{\∂}^2}{{\∂y}^2}\left(\frac{\∂T}{{\∂T}_{k,n}}\right)]---\left(10\right)$

灵敏度方程的边界条件和初始条件为：

圆管外侧

$-\mathrm{\λ}{\frac{\∂T}{\∂n}\left(\frac{\∂T}{{\∂T}_{k,n}}\right)}_w=h_{\mathrm{out}}\left(\frac{\∂T}{{\∂T}_{k,n}}\right)---\left(11\right)$

圆管内侧

$\frac{\∂T}{{\∂T}_{k,n}}=l(s,t)---\left(12\right)$

函数l(s，t)的值采用如下方式计算：

当s＝k，t＝n时，l(s，t)＝1

其它情况，l(s，t)＝0    (13)

初始条件

$\frac{\∂T(x,y,0)}{T_{k,n}}=0---\left(14\right)$

计算γ时需要知道

其计算式为：

$\frac{\∂J}{{\∂T}_{k,n}}=\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(T_{i,j,\mathrm{cal}}-T_{i,j,\mathrm{mea}})\frac{{\∂T}_{i,j}}{{\∂T}_{k,n}}---\left(15\right)$

共轭梯度法的收敛目标是使

$J\left(T_{k,n}^b\right)<\mathrm{\&epsiv;}---\left(16\right)$

其中，ε为趋向于0的正数.

具体地，共轭梯度法的求解步骤如下：

步骤S01：由初始温度T⁰求解灵敏度方程得到灵敏度系数

步骤S02：求解导热方程(1)得到管外壁温度计算值；

步骤S03：根据管外壁温度计算值和测量值，通过公式(16)检查收敛目标是否达到；如果达到收敛目标，则停止迭代，否则，转到步骤S04；

步骤S04：分别由公式(7)(8)(9)(15)计算得到迭代参数d，γ，β，

步骤S05：由公式(6)计算得到下一迭代步的管内壁温度T^b+1，返回到步骤S02。

下面将列举本发明的具体实施例。

图1所示为二维左右对称管内流动问题示意图，为了能够在通过本发明的技术方案由管外壁温度计算内壁温度的同时检验本发明技术方案中共轭梯度法的实施效果，首先可以通过实验测量得到设定节点处管内壁节点温度准实验值和管外壁节点温度准实验值，本实施例选取0°、60°、90°、120°、180°5个外壁测温点，由于流动为左右对称，只需考虑半个圆周；管内壁测温点位置选择如图2所示，管外壁测温点位置如图3所示；根据这两组准实验值得到的管内壁和外壁温度随时间的变化分别如图4、图5所示。根据已有实验测量经验，获得一组管内壁温度初始值，利用本发明的方法，计算得到管内壁温度计算值，再把管内壁温度准实验值与根据本发明技术方案的管内壁温度计算值相比较，获得其相符程度，从而实现对本发明技术方案的检验。其中，求解导热正问题和求解导热反问题的计算过程都可以利用FORTRAN语言编制程序或其他可行的方式进行。

本实施例中作为测量对象的二维圆管的管外径为60.5mm，管内径为43.1mm。由于所述二维左右对称管内流动问题为左右对称，只需考虑半个圆周即可。

本实施例采用有限容积法对半个圆管进行计算。采用有限容积法进行计算时需要对计算区域进行网格划分，本实施例中网格划分采用四边形网格，沿圆周向划分19个网格节点，沿径向6个网格节点，由此半圆面上总共存在19×6＝114个节点，如图6所示。

首先通过迭代求解计算得到所述114个节点的温度随时间的变化：其中，内壁19个节点的温度通过共轭梯度法计算得到，其它114-19＝95个节点的温度通过求解导热方程(公式(1)(2)(3))得到；然后，根据外壁节点计算值与实验测量值，通过公式(5)计算目标函数值，判断其是否达到收敛要求。

其中，由于本实施例仅在外壁设置了5个实验测量点，而网格划分得到的外壁节点为19个，其中所述5个实验测量点与5个网格划分得到的节点重合，剩余的19-5＝14个节点的值可以由实验测点进行线性插值得到，插值得到的温度值可以被视为准实验值。

需要说明的是，网格划分得到的网格节点也可以不与实验测点重合，即网格节点的温度值可以全部由插值得到。但由于插值会带来误差，而网格点与实验点重合的划分方式可减小误差，因此，优选地，采用网格点与实验点重合的方式。

其中，管外壁温度的计算在第三类边界条件下进行，即已知管外侧环境流体温度T_sir＝20℃，对流换热系数h_out＝1w/(m·k)，圆管的热扩散系数α＝4.46×10^-6m²/s，导热系数λ＝19.35w/(m·k)。由设定的管内壁温度初始值通过由程序实现的导热反问题的迭代求解得到管内壁温度计算值；本实施例中，设置收敛目标值ε＝0.05，为达到该收敛目标值，迭代共进行了1160次，在奔腾IV 1.8G CPU/1G内存的个人电脑上耗时约12小时。图7为0°节点位置处利用本发明的方法得到的管内壁温度计算值与准实验值的比较，从图中可以看到，共轭梯度法计算值与准实验值相符较好；图8为图7中60秒至80秒数据的放大图，如图中所示，共轭梯度法不仅能准确的捕捉到管内壁温度值的较小波动，对于60秒至80秒之间的较大波动也能精确的捕捉到。通过比较其它几个节点计算值与准实验值，可以发现，共轭梯度法还能准确捕捉到其它位置处温度的波动，且计算值与准实验值相符较好。图9和图10分别为管内壁60°和120°节点位置处共轭梯度法计算值与准实验值的比较，如图中所示，根据本发明的计算方法的计算值与准实验值相符较好。

从以上计算和比较结果可以看到，本发明可以准确的由管外壁温度计算得到管内壁温度随时间的变化，计算精度高，计算量较小，较好的解决了工业生产、科研工作中对于管内流动管内壁温度测量困难的问题，具有广泛的应用前景。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (8)

1、一种由圆管外壁温度计算内壁温度的方法，其特征在于，所述方法包括：

在二维非稳态管内流动下，获得多个外壁节点的外壁温度测量值；

根据所述外壁温度测量值得到多个内壁节点的内壁温度结果值。

2、如权利要求1所述的方法，其特征在于，所述根据所述外壁温度测量值得到多个内壁节点的内壁温度结果值的步骤包括：

设置内壁温度初始值；

根据所述内壁温度初始值计算得到所述多个外壁节点的外壁温度计算值；

判断所述外壁温度计算值与所述外壁温度测量值是否相符，若相符，则所述内壁温度初始值为内壁温度结果值；若不相符，通过迭代计算获得结果值。

3、如权利要求1所述的方法，其特征在于，所述多个外壁和内壁节点由网格划分方法获得。

4、如权利要求2所述的方法，其特征在于，所述多个外壁节点的外壁温度计算值通过求解圆管二维非稳态导热方程获得。

5、如权利要求2所述的方法，其中，所述判断所述外壁温度计算值与所述外壁温度测量值是否相符的步骤进一步包括：

根据所述外壁温度计算值与所述外壁温度测量值得到当前的目标函数的值；

判断所述当前的目标函数的值是否小于一个预先确定的阈值，若是，则相符；若否，则不相符。

6、如权利要求5所述的方法，其特征在于，所述目标函数为： $J=\frac{1}{2}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{[T_{\mathrm{Cal},m,n}-T_{\mathrm{Exp},m,n}]}^2,$ 其中，T_cal，m，n是外壁温度计算值，T_Exp，m，n是外壁温度测量值，M、N分别为测量值节点位置数和时间点数。

7、如权利要求2所述的方法，其特征在于，所述通过迭代计算获得结果值的步骤之前还包括：根据迭代方程获得下一步迭代的内壁温度初始值；所述迭代方程为： $T_{k,n}^{b+1}=T_{k,n}^b-{\mathrm{\&beta;}}^b{d}_{k,n}^b$ 其中，b为迭代次数，k为管内壁待计算节点的位置，β为迭代步长，n为时间点，d为搜索方向。

8、如权利要求7所述的方法，其特征在于，所述迭代步长与表示内外壁节点温度变化关系的灵敏度系数有关，所述灵敏度系数由灵敏度方程获得，其中，所述灵敏度方程在圆管内壁的边界条件为： $\frac{\&PartialD;T}{\&PartialD;T_{k,n}}=l(s,\mathrm{\&tau;}),$ 并有， $\begin{array}{cc}\mathrm{if\; s}=k,\mathrm{\&tau;}=n.& l(s,\mathrm{\&tau;})=1\\ \mathrm{otherwise}& l(s,\mathrm{\&tau;})=0\end{array};$ 其中，

为灵敏度系数。

Images