一种求解区间热对流扩散问题的高精度数值方法

摘要

本发明公开了一种求解区间热对流扩散问题的高精度数值方法，步骤如下：建立热对流扩散问题的传热控制方程；引入区间变量对传热模型中输入参数的不确定性进行表征，建立热对流扩散问题的区间控制方程；利用勒让德多项式对区间控制方程中的温度响应进行近似表示；根据张量积运算，确定高维空间中的配点集合；利用有限元程序计算所有配点处的温度响应，建立关于温度响应近似表达式中展开系数的线性方程组，并采用最小二乘法进行求解；基于多项式函数的光滑性，确定温度响应近似表达式的极值点，得到区间温度响应的上下界。本发明可系统化解决含有区间不确定参数的热对流扩散问题，有效提高了区间数值方法的计算精度。

说明书

技术领域

本发明属于机械工程领域，具体涉及一种求解区间热对流扩散问题的高精度数值 方法。

背景技术

在自然界和各种生产技术领域中，通过流体流动来进行热量的传递是一种极其普 遍的物理现象。尤其是在航空航天、能源化工等工业装备以及电子器件等精密产品的生产 制造过程中，如何有效的实现热量传递，已成为产品设计的一个重要方面。现有关于热分析 的许多研究都是针对确定性模型进行的，没有考虑模型输入参数的不确定性。实际工程中， 由于制造工艺的限制、测量误差以及认知的局限，结构的材料属性、外部载荷和边界条件等 物理参数不可避免的受到多种不确定因素的影响，使得系统的温度响应也表现出一定的不 确定性。

基于概率理论和数理统计的随机建模与计算方法在不确定系统的分析、试验和设 计方面已经取得了不少研究成果，成功解决了许多工程实际问题。但用随机理论求解问题 时，往往假定不确定因素为随机变量或随机过程，这就事先需要大量的试验信息确定其概 率分布规律。然而，多数实际情况下，获取充足的试验数据往往代价昂贵。如此一来，信息的 缺乏使得概率模型不能真实反映客观实际，也就使得随机不确定分析方法失去意义。实际 工程中的多数情况，设计者往往只关注某些响应的变化幅度，而获得不确定模型输入参数 的取值范围相对概率密度函数来说要容易的多，且所需的不确定信息也大大减少，这就有 效拓宽了不确定分析的应用范围。近些年来，利用区间理论和分析方法来处理工程中的不 确定因素受到学者们越来越多的重视。目前，将区间理论与有限元计算方法相结合衍生出 来的区间有限元法在不确定结构的静、动力特性分析方面已经取得了不少研究成果，但对 具有区间参数的传热问题分析的文献却十分少见。另外，传统区间分析方法因区间运算所 引起的区间扩张问题还比较严重，计算精度还亟待提高。因此，如何建立高精度的区间分析 方法对不确定传热问题进行数值求解，是目前学术领域的一个研究热点，对于弥补现有传 热数值计算方法的不足，具有重要的工程应用价值。

发明内容

本发明所要解决的技术问题为：克服现有技术在热对流扩散问题求解中存在的不 足，充分考虑传热问题中的区间不确定因素，基于多项式逼近思想和配点分析理论，提出了 一种预测温度响应区间变化范围的高精度数值计算方法，可系统化解决含有区间不确定参 数的温度场预测问题，在保证计算精度的同时，进一步降低了传统抽样方法的计算耗费。

本发明为解决上述技术问题采用的技术方案为：一种求解区间热对流扩散问题的 高精度数值方法，包括以下步骤：

步骤一：根据传热模型建立热对流扩散问题的传热控制方程；

步骤二：引入区间变量对传热模型中输入参数的不确定性进行表征，根据步骤一 中的传热控制方程建立热对流扩散问题的区间控制方程；

步骤三：利用勒让德多项式对步骤二区间控制方程中的温度响应进行逼近，得到 区间温度响应的近似表达式；

步骤四：根据张量积运算，由一维空间的配点集合确定整个高维空间的配点集合； 所述高维是指多于一维，维数等于区间变量的个数；

步骤五：利用有限元程序计算步骤四配点集合中所有配点处的温度响应，建立关 于步骤三区间温度响应近似表达式中展开系数的线性方程组，并采用最小二乘法对此线性 方程组进行求解，得到展开系数的一组值；

步骤六：将步骤五中得到的展开系数的一组值代回到步骤三温度响应的近似表达 式中，基于多项式函数的光滑性，确定此近似表达式的极值点，进而得到区间温度响应的上 下界。

其中，所述步骤三中利用勒让德多项式对温度响应进行近似表示，多项式的截断 阶数并不是固定不变的，根据逼近精度要求选取，截断阶数越高，逼近精度就越高。

其中，所述步骤四中配点方案的选取并不是固定不变的，根据计算耗费和计算精 度的要求来选取配点数量，配点数量越多，计算精度就越高，而计算耗费就越大。

其中，所述步骤四中的高维是指多于一维，维数等于区间变量的个数。

上述各步骤具体包括以下过程：

步骤一：根据传热模型建立热对流扩散问题的传热控制方程：

$\rho cu\frac{\partial T\left(x\right)}{\partial x}=k\frac{{\partial }^{2}T\left(x\right)}{\partial x^2}+Q\left(x\right)$

其中x表示物理坐标，T(x)表示温度响应，ρ,c,k分别表示材料的密度、比热容和热 传导系数，u为传热流体的流动速度，Q(x)表示系统的热源强度。

步骤二：引入n个区间变量对传热模型中输入参数的不确定性进行表 征，并记为向量的形式其中上标I是区间符号，和表示 区间变量的下界和上界，和称作区间变量的中点和半径， 为标准区间变量根据步骤一中的控制方程建立热对流扩散问题的区间控制方 程：

$\rho \left({\alpha }^I\right)c\left({\alpha }^I\right)u\left({\alpha }^I\right)\frac{\partial T(x,{\alpha }^I)}{\partial x}=k\left({\alpha }^I\right)\frac{{\partial }^{2}T(x,{\alpha }^I)}{\partial x^2}+Q(x,{\alpha }^I)$

步骤三：利用有限阶的勒让德多项式对步骤二区间控制方程中的温度响应T(x, α^I)进行逼近，得到区间温度响应的近似表达式：

$T(x,{\alpha }^I)\approx T_N(x,{\alpha }^I)=\underset{\left|i\right|\le N}{\Sigma }{w}_i\left(x\right){\Phi }_i\left({\alpha }^I\right)$

其中Φ_i(α^I)为事先选定的勒让德正交多项式基底，w_i(x)为对应的展开系数，i＝ (i₁,i₂,...,i_n)表示多维指标，且满足|i|＝i₁+i₂+...+i_n，N为此多项式的截断阶数。根据多 项式理论，上述近似表达式中展开项的个数可用变量数n和截断阶数N表示为

步骤四：根据张量积运算，由一维空间的配点集合确定整个高维空间的配点集合。 首先，对于一维区间变量来说，设定其配点数量为m_i，则各个配点的具体位置为：

其次，用点集表示一维区间变量内所有配点组成的集 合，那么对于n个变量组成的高维空间而言，直接利用张量积运算，可得其配点集合Θ

$\Theta ={\Theta }_1^{m_1}\times {\Theta }_2^{m_2}\times ...\times {\Theta }_n^{m_n}$

而配点总数M为：

$M=\underset{i=1}{\overset{n}{\Pi }}{m}_i=m_1\times m_2\times ...\times m_n$

因此将Θ改写为的形式，用来表示高维空间所有的配点

步骤五：利用有限元程序计算步骤四配点集合中所有配点处的温度响应，建立关 于步骤三区间温度响应近似表达式中展开系数的线性方程组，并采用最小二乘法对此线性 方程组进行求解，得到展开系数的一组值。首先，步骤二中所建立的区间控制方程在配点 处可改写为：

$\rho \left({\beta }_j^{node}\right)c\left({\beta }_j^{node}\right)u\left({\beta }_j^{node}\right)\frac{\partial T(x,{\beta }_j^{node})}{\partial x}=k\left({\beta }_j^{node}\right)\frac{{\partial }^{2}T(x,{\beta }_j^{node})}{\partial x^2}+Q(x,{\beta }_j^{node}),j=1,2,...,M$

利用有限元程序对上述方程进行求解，可以得到所有配点处的温度响应

其次，基于步骤三中区间温度响应近似表达式，可以建立关于所有展开系数w_i(x) 的线性方程组：

其中表示多项式基底函数Φ_i(α^I)在配点处的取值。

然后，利用最小二乘法求解上述方程组，得到展开系数w_i(x)的一组值。

步骤六：将步骤五中计算得到的展开系数w_i(x)的一组值代回到步骤三温度响应 的近似表达式中，基于多项式函数T_N(x,α^I)的光滑性，令其一阶导数为零，易确定其极值点， 连同边界点一起记为其中r为极值点和边界点的总体数量。比较这r个点处的 温度值大小，最终确定区间温度响应T(x,α^I)的下界和上界

$\underset{‾}{T}(x,{\alpha }^I)\approx {\underset{‾}{T}}_N(x,{\alpha }^I)=\underset{j=1,...,r}{min}{T}_N(x,{\alpha }_j^0)$

$\bar{T}(x,{\alpha }^I)\approx {\bar{T}}_N(x,{\alpha }^I)=\underset{j=1,...,r}{max}{T}_N(x,{\alpha }_j^0)$

本发明与现有技术相比的优点在于：

(1)与传统的热对流扩散问题分析方法相比，所提出的数值计算方法充分考虑实 际工程中材料属性、外部载荷和边界条件的区间不确定性，计算结果对温度场分析具有更 重要的指导意义。

(2)利用高阶勒让德多项式对温度响应进行近似表示，可有效提高逼近精度。同 时，利用多项式函数的光滑性，可快速确定其极值点，进而得到温度响应的上下界。

(3)基于配点理论，可以充分利用原有确定性模型的有限元计算程序而无需对其 做进一步的修改，保证了计算的可移植性。

(4)本发明操作简单，实施方便，有效提高了计算精度。

附图说明

图1为本发明的一种求解区间热对流扩散问题的高精度数值方法流程图；

图2为本发明的三维热交换器模型示意图；

图3为上部面板中心线区间温度响应示意图；

图4为流管中心线区间温度响应示意图。

具体实施方式

下面结合附图和实施例对本发明做进一步说明。

本发明可系统化解决含有区间不确定参数的热对流扩散问题，有效提高了区间数 值方法的计算精度。本发明适用于含有区间不确定参数的热对流扩散问题的温度场预测。 本发明实施方式以某三维热交换器模型为例，具体说明所述的一种求解区间热对流扩散问 题的高精度数值方法。另外，此三维热交换器模型的区间温度响应数值计算方法可以推广 到其他含有区间参数的热对流扩散问题温度场预测中。

一种求解区间热对流扩散问题的高精度数值方法的计算过程如图1所示，根据具 体的传热模型建立热对流扩散问题的传热控制方程，引入区间变量对传热模型中输入参数 的不确定性进行表征，建立热对流扩散问题的区间控制方程，利用勒让德多项式对温度响 应进行近似表示，同时根据张量积运算，确定配点集合，利用有限元程序计算所有配点处的 温度响应，对温度响应近似表达式中展开系数进行求解，并基于函数的光滑性，快速得到区 间温度响应的上下界。可分为如下几个步骤进行：

步骤一：考虑图2所示的长为300mm的热交换器传热模型，截面为40mm×20mm的矩 形，中间有30mm×15mm的方孔通过速度为u的冷却空气，入口处8的空气温度为T_s，结构的上 面板7承受密度为q_s＝20000×sin(ω)W/m²的热流载荷，在上面板中心线上选定编号为1～3 的三个点，在流管中心线上选定编号为4～6的三个点作为温度场的观测点。根据传热模型 建立热对流扩散问题的传热控制方程：

$\rho cu\left(\frac{\partial T\left(x,y,z\right)}{\partial x}+\frac{\partial T\left(x,y,z\right)}{\partial y}+\frac{\partial T\left(x,y,z\right)}{\partial z}\right)=k\left(\frac{{\partial }^{2}T\left(x,y,z\right)}{\partial x^2}+\frac{{\partial }^{2}T\left(x,y,z\right)}{\partial y^2}+\frac{{\partial }^{2}T\left(x,y,z\right)}{\partial z^2}\right)$

其中x,y,z表示三个空间方向上的物理坐标，T(x,y,z)表示温度响应，ρ,c,k分别 表示空气的密度、比热容和热传导系数，u为空气的流动速度。

步骤二：由于材料制造工艺的限制以及测量的误差，所有模型输入参数均含有一 定的区间不确定性，引入六个区间变量对不确定性进行表征ρ^I＝[1.3,1.5]kg/m³，c^I＝ [900,1100]J/(kg·℃)，k^I＝[0.023,0.029]W/(m·℃)，u^I＝[1.8,2.2]m/s，T_s^I＝[18,22] ℃，ω^I＝[1.5,2.5]。将所有区间变量统一表示为向量形式其 中上标I是区间符号，和表示区间变量的下界和上界，和称作区间变量的中点和半径，为标准区间变量根据步骤一中的控制方程建立 热对流扩散问题的区间控制方程：

$\left(\begin{array}{c}\rho \left({\alpha }^I\right)c\left({\alpha }^I\right)u\left({\alpha }^I\right)\left(\frac{\partial T(x,y,z,{\alpha }^I)}{\partial x}+\frac{\partial T(x,y,z,{\alpha }^I)}{\partial y}+\frac{\partial T(x,y,z,{\alpha }^I)}{\partial z}\right)\\ =k\left(\frac{{\partial }^{2}T(x,y,z,{\alpha }^I)}{\partial x^2}+\frac{{\partial }^{2}T(x,y,z,{\alpha }^I)}{\partial y^2}+\frac{{\partial }^{2}T(x,y,z,{\alpha }^I)}{\partial z^2}\right)\end{array}\right)$

步骤三：根据区间变量的分布特点，选用勒让德多项式对步骤二区间控制方程中 的温度响应T(x,y,z,α^I)进行逼近，截断阶数设定为N＝3，得到区间温度响应的近似表达 式：

$T(x,y,z,{\alpha }^I)\approx T_N(x,y,z,{\alpha }^I)=\underset{\left|i\right|\le N}{\Sigma }{w}_i(x,y,z){\Phi }_i\left({\alpha }^I\right)$

其中Φ_i(α^I)为事先选定的勒让德正交多项式基底，w_i(x,y,z)为对应的展开系数， i＝(i₁,i₂,...,i₆)表示多维指标，且满足|i|＝i₁+i₂+...+i₆。此时上述近似表达式中展开 项的个数为

步骤四：根据张量积运算，由一维空间的配点集合确定整个高维空间的配点集合。 首先，对于一维区间变量来说，设定其配点数量为m_i＝5，则各个配点的具体位置 为：

${\beta }_j^{\left(i\right)}={\alpha }_i^c-cos\frac{\pi (j-1)}{4}·{\mathrm{\Delta \alpha }}_i,j=1,2,...,5$

其次，用点集表示一维区间变量内所有配点组成的集 合，那么对于6个变量组成的六维空间而言，直接利用张量积运算，可得其配点集合Θ：

$\Theta ={\Theta }_1^{m_1}\times {\Theta }_2^{m_2}\times ...\times {\Theta }_6^{m_6}$

而配点总数为将Θ改写为的形式，用来表 示六维空间所有的配点

步骤五：利用有限元程序计算步骤四配点集合中所有配点处的温度响应，建立关 于步骤三区间温度响应近似表达式中展开系数的线性方程组，并采用最小二乘法对此线性 方程组进行求解，得到展开系数的一组值。首先，步骤二中所建立的区间控制方程在配点 处可改写为：

$\begin{array}{c}\rho \left({\beta }_j^{node}\right)c\left({\beta }_j^{node}\right)u\left({\beta }_j^{node}\right)\left(\frac{\partial T(x,y,z,{\beta }_j^{node})}{\partial x}+\frac{\partial T(x,y,z,{\beta }_j^{node})}{\partial y}+\frac{\partial T(x,y,z,{\beta }_j^{node})}{\partial z}\right)\\ =k\left(\frac{{\partial }^{2}T(x,y,z,{\beta }_j^{node})}{\partial x^2}+\frac{{\partial }^{2}T(x,y,z,{\beta }_j^{node})}{\partial y^2}+\frac{{\partial }^{2}T(x,y,z,{\beta }_j^{node})}{\partial z^2}\right)\end{array},j=1,2,\mathrm{...},15625$

利用软件Nastran中的有限元程序对上述传热问题进行求解，可以得到所有配点 处的温度响应

其次，基于步骤三中区间温度响应近似表达式，建立关于所有展开系数w_i(x,y,z) 的线性方程组：

其中表示多项式基底函数Φ_i(α^I)在配点处的取 值。

然后，利用最小二乘法求解上述方程组，得到展开系数w_i(x,y,z)的一组值。

步骤六：将步骤五中计算得到的展开系数w_i(x,y,z)的一组值代回到步骤三温度 响应的近似表达式中，利用多项式函数T_N(x,y,z,α^I)的光滑性，令其一阶导数为零，确定其 极值点，连同边界点一起记为其中r为极值点和边界点的总体数量。比较这r个 点处的温度值大小，最终确定区间温度响应T(x,y,z,α^I)的下界和上界

$\underset{‾}{T}(x,y,z,{\alpha }^I)\approx {\underset{‾}{T}}_N(x,y,z,{\alpha }^I)=\underset{j=1,...,r}{min}{T}_N(x,y,z,{\alpha }_j^0)$

$\bar{T}(x,y,z,{\alpha }^I)\approx {\bar{T}}_N(x,y,z,{\alpha }^I)=\underset{j=1,...,r}{max}{T}_N(x,y,z,{\alpha }_j^0)$

六个观测点处温度响应的计算结果如表1所示。与样本数为10⁶的传统蒙特卡洛抽 样方法对比可以看出，本发明方法的计算误差小于1％，计算精度完全满足工程需求。另外， 从样本数量上看，本发明方法的样本数为15625，计算耗费远远小于蒙特卡洛方法。

表1观测点处区间温度响应上下界

除了上述六个观测点外，沿x轴方向，上部面板和流管中心线区间温度响应如图3 和图4所示，横坐标表示沿x轴方向的空间位置，纵坐标表示空间位置处的温度值，实线和虚 线分别表示蒙特卡洛抽样方法和本发明方法计算得到的结果。可以看出，本发明方法计算 得到的温度响应上下界曲线与传统蒙特卡洛抽样得到的参考值吻合程度很好，计算结果真 实可信。用本发明方法可以解决含有区间不确定输入参数的热对流扩散问题，计算精度高， 此功能是一般商用软件所不能实现的。

以上所述的仅为本发明的较佳实施例而已，本发明不仅仅局限于上述实施例，凡 在本发明的精神和原则之内所作的局部改动、等同替换、改进等均应包含在本发明的保护 范围之内。