CN113361155B - 一种热物性参数辨识结果置信区间估计方法及系统 - Google Patents

一种热物性参数辨识结果置信区间估计方法及系统 Download PDF

Info

Publication number
CN113361155B
CN113361155B CN202110522174.4A CN202110522174A CN113361155B CN 113361155 B CN113361155 B CN 113361155B CN 202110522174 A CN202110522174 A CN 202110522174A CN 113361155 B CN113361155 B CN 113361155B
Authority
CN
China
Prior art keywords
parameter
thermophysical
calculating
confidence interval
cdf
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110522174.4A
Other languages
English (en)
Other versions
CN113361155A (zh
Inventor
窦怡彬
李煜
陆云超
李宗阳
刘陆广
任智毅
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Institute of Electromechanical Engineering
Original Assignee
Shanghai Institute of Electromechanical Engineering
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shanghai Institute of Electromechanical Engineering filed Critical Shanghai Institute of Electromechanical Engineering
Priority to CN202110522174.4A priority Critical patent/CN113361155B/zh
Publication of CN113361155A publication Critical patent/CN113361155A/zh
Application granted granted Critical
Publication of CN113361155B publication Critical patent/CN113361155B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08—Thermal analysis or thermal optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Investigating Or Analyzing Materials Using Thermal Means (AREA)

Abstract

本发明提供一种热物性参数辨识结果置信区间估计方法及系统，涉及结构热防护设计中的材料热物性参数辨识技术领域，包括：步骤S1：建立结构热传导计算模型；步骤S2：优化得到热物性参数pbest和最优卡方函数χ(pbest)；步骤S3：给定热物性参数置信度pc和增量步长dp，迭代计算直到F分布函数的连续分布函数值pcdf>pc；步骤S4：计算相应的置信区间上限
Figure DDA0003064407530000011
使得pcdf＝pc；步骤S5：计算置信区间下限
Figure DDA0003064407530000012
步骤S6：遍历所有待辨识热物性参数，重复步骤S3～步骤S5，得到热物性参数辨识结果的置信区间；步骤S7：开展热防护设计温度响应计算。本发明能够对于不服从正态分布的热物性参数估计其置信区间，置信度更高，估计出的置信区间能体现出辨识参数不确定度的非线性程度。

Description

一种热物性参数辨识结果置信区间估计方法及系统
技术领域
本发明涉及结构热防护设计中的材料热物性参数辨识技术领域，具体地，涉及一种热物性参数辨识结果置信区间估计方法及系统。
背景技术
随着导弹向着高速、远程的方向发展，气动加热及热防护设计已经成为高超声速飞行器研制过程中一个重要且必不可少的工作。在防热方案设计阶段，对弹体结构大面积部位开展传热分析，并根据得到的温度响应给出相应的防热方案。结构的温度响应结果受材料热物性参数影响，而这些热物性参数往往是无法得到或者不全面、不准确的，只能用类似材料的参数近似求解。另外在产品研制阶段，通过大量的地面热试验和飞行试验获得的温度响应数据往往只是提供给总体作为验证防热方案可行的佐证，其中蕴含的更多信息没有被充分挖掘，需要开发基于地面或者飞行试验温度测量结果的材料热物性参数辨识研究，为后续型号的防热设计和热传导计算模型修正提供依据。
基于试验测量得到的温度响应数据来辨识材料热物性参数，属于热传导逆问题研究范畴，大致可以分为两大类方法即基于梯度的方法和基于随机的方法。这两类方法在辨识数据没有噪声的情况下都可以得到较好的辨识结果，然后真实的试验数据，特别是风洞试验数据和飞行实测数据是含有测量噪声的，且试验测量数据往往只有1-2组，无法通过多次辨识取平均值来消除噪声对辨识结果的影响，因此需要对辨识结果进行置信区间分析得到一定范围内可信的热物性参数。
公开号为CN105892481A的发明专利，公开了一种三捷联惯组量化动态阈值置信区间估计方法，属于基于分位数法的阈值置信区间估计方法，首先进行分位数定义，其次将通过蒙特卡洛打靶方法生成带有三捷联惯组误差模型的弹道数据作为样本数据；之后对样本数据按照数值大小进行升序排序，将相同的数值进行合并，计算每一个数据在原样本数据中出现的次数占整体数据的比例，得到每个样本数据的出现概率；最后进行阈值置信区间的估计。
目前热物性参数辨识结果置信区间估计主要是通过计算渐进标准差来表示68％置信度下的置信区间。该方法只能对服从正态分布的热物性参数进行置信区间估计且得到的置信区间是关于辨识结果对称的，无法体现辨识参数不确定度的非线性效应。为了方便工程设计人员进一步认识测量数据噪声水平对辨识结果的影响，分析每个热物性参数不确定度的非线性程度，有必要提供一种新的热物性参数辨识结果置信区间估计方法。
发明内容
针对现有技术中的缺陷，本发明提供一种热物性参数辨识结果置信区间估计方法及系统。
根据本发明提供的一种热物性参数辨识结果置信区间估计方法及系统，所述方案如下：
第一方面，提供了一种热物性参数辨识结果置信区间估计方法，所述方法包括：
步骤S1：对几何模型进行网格划分，定义网格的单元属性和材料属性，设置边界条件，基于有限体积法建立结构的一维热传导计算模型；
步骤S2：基于目标温度响应，采用L-M算法对待辨识热物性参数进行优化，辨识得到热物性参数的最优解pbest和最优卡方函数χ(pbest)；
步骤S3：给定热物性参数置信度pc和增量步长dp，迭代计算直到F分布函数的连续分布函数值pcdf>pc；
步骤S4：采用二分法在区间[pbest,i pnew,i]内计算相应的置信区间上限
Figure GDA0003724070470000021
使得连续概率密度函数pcdf等于给定的置信度pc；
步骤S5：令dp＝-dp重复步骤S3～步骤S4，计算置信区间下限
Figure GDA0003724070470000022
步骤S6：令dp＝-dp遍历所有待辨识热物性参数，重复步骤S3～步骤S5，得到所有给定置信度下热物性参数辨识结果的置信区间；
步骤S7：取热物性参数辨识结果置信区间的中值和上下限，开展热防护设计温度响应计算，得到温度响应的上下边界供设计人员参考。
优选的，所述步骤S1包括：
所述结构的一维热传导计算模型通用表达式为：
Figure GDA0003724070470000031
式中ρ表示材料密度；cp表示比热容；λ表示热导率；
Figure GDA0003724070470000032
表示内热源，根据具体物理模型取值，结构没有内部加热源或非相变材料情况下该项等于0；T表示温度；r表示一维方向坐标；m取0、1和2分别对应笛卡尔坐标系、圆柱坐标系和球坐标系；t表示时间。
优选的，所述步骤S1还包括：
采用有限体积法对一维结构进行离散，得到如下隐式离散格式：
Figure GDA0003724070470000033
其中，
Figure GDA0003724070470000034
Figure GDA0003724070470000035
Figure GDA0003724070470000036
上式中aw、aP和aE分别为温度前的系数；W、P和E分别表示当前网格点的西侧、当前网格点和当前网格点的东侧；下标i表示第i个网格节点；上标n表示第n时刻；Sc和Sp分别表示内热源
Figure GDA0003724070470000037
线性化的常数项和线性项系数；Δt为时间步长；(·)′表示调和平均参数。变量的上标w和e分别表示变量取第i个网格单元西侧和东侧边界上的值；A表示第i个网格节点处单元横截面积；ri表示第i个网格节点处的坐标；δr表示网格节点到网格单元边界的距离。
左边界条件离散后有：
Figure GDA0003724070470000041
其中，
Figure GDA0003724070470000042
Figure GDA0003724070470000043
上式中上标1表示第一个网格点；其中
Figure GDA0003724070470000044
h、ε、σ、Tr和TAMBIENT分别表示热流密度、对流换热系数、发射率、Stefen-Bolzman常数、对流换热温度和辐射背景温度。
右边界条件离散后有：
Figure GDA0003724070470000045
其中，
Figure GDA0003724070470000046
Figure GDA0003724070470000047
上式中上标M表示第M个网格节点，即最后一个网格节点。
在实际问题中，热导率和比热容是关于温度的函数，其能够采用多项式形式表示。
优选的，所述步骤S2包括：
基于目标温度响应，采用L-M算法对带辨识参数向量进行优化，L-M算法通过迭代求解如下方程来优化参数向量p，
Figure GDA0003724070470000048
上标j表示迭代步数，上式中
Figure GDA0003724070470000049
为雅各比矩阵，Ω为对角矩阵，μ为阻尼因子，h为增量步长，W是测量误差倒数构成的协方差矩阵，
Figure GDA0003724070470000051
为实际测量到的随时间变化的温度向量。Ω矩阵的具体形式为
Ω＝diag(JTWJ)
通过对比两个迭代增量步之间目标函数的差值z(p)来判断第j个迭代增量步hj是否合适，其具体表达式为：
Figure GDA0003724070470000052
其中，tol为给定的误差，满足上式则接受hj，得到新的p，进入下一步迭代；阻尼因子根据如下方式更新：
Figure GDA0003724070470000053
对于雅各比矩阵J的计算，采用差分方法代替微分方法进行求解。计算如上步骤，得到下一个迭代步的待辨识参数pj+1，并计算相应的卡方函数χ(pj+1)，如果小于给定的收敛标准则结束迭代计算，令pbest＝pj+1即为辨识得到的热物性参数；否则更新阻尼因子继续迭代计算。
优选的，所述步骤S3包括：初始时刻令
Figure GDA0003724070470000054
其中上标k表示迭代步数；对第i个热物性参数开始迭代计算置信区间上下限；
令
Figure GDA0003724070470000055
将
Figure GDA0003724070470000056
作为初始值并采用L-M算法进行优化得到新的热物性参数
Figure GDA0003724070470000057
和新的卡方函数
Figure GDA0003724070470000058
根据卡方函数
Figure GDA0003724070470000059
和最优卡方函数χ(pbest)计算F分布函数的连续概率密度值pcdf；
如果pcdf≤pc，则k＝k+1，并重复步骤S3，如果pcdf>pc则停止计算。
优选的，所述步骤S3还包括：
初始时刻令k＝0，i＝1，
Figure GDA00037240704700000510
对第i个热物性参数开始迭代计算置信区间上下限；
令
Figure GDA00037240704700000511
并将
Figure GDA00037240704700000512
作为初始条件采用步骤S2辨识热物性参数，得到新的热物性参数
Figure GDA00037240704700000513
和新的卡方函数
Figure GDA00037240704700000514
此时L-M优化过程中
Figure GDA00037240704700000515
固定不变，根据卡方函数
Figure GDA00037240704700000516
和最优卡方函数χ(pbest)计算F分布函数的连续概率密度值pcdf，F分布的连续概率密度函数的具体表示式为：
Figure GDA0003724070470000061
其中，Г为伽马函数；v1和v2为F分布函数在分子和分母上的自由度，令v1＝1，v2＝nt-np，其中，nt为温度响应数据个数，np为热物性参数个数；x为积分上限，令
Figure GDA0003724070470000062
若pcdf≤pc，则k＝k+1，并重复步骤S3；若pcdf>pc则停止计算。
优选的，所述步骤S4包括：
对于第i个参数，使用二分法在区间[pbest,i pnew,i]内计算相应的置信区间上限
Figure GDA0003724070470000063
使得连续概率密度函数pcdf等于给定的置信度pc；
令p*＝0.5(pbest,i+pnew,i)，
Figure GDA0003724070470000064
且
Figure GDA0003724070470000065
上标0表示初始迭代步；
以
Figure GDA0003724070470000066
为初始条件调用步骤S2优化得到新的
Figure GDA0003724070470000067
并计算相应的卡方函数
Figure GDA0003724070470000068
结合最优卡方函数χ(pbest)计算F分布函数的连续概率密度值pcdf，若|pcdf-pc|<ε，则得到上限
Figure GDA0003724070470000069
否则根据pcdf-pc是否大于0，令pbest,i＝p*或者pnew,i＝p*，重复步骤S4直到p*满足要求。
第二方面，提供了一种热物性参数辨识结果置信区间估计系统，所述系统包括：
模块M1：对平板试片模型进行网格划分，定义网格的单元属性和材料属性，设置边界条件，基于有限体积法建立结构的一维热传导计算模型；
模块M2：基于目标温度响应，采用L-M算法对待辨识热物性参数进行优化，辨识得到热物性参数的最优解pbest和最优卡方函数χ(pbest)；
模块M3：给定热物性参数置信度pc和增量步长dp，迭代计算直到F分布函数的连续分布函数值pcdf>pc；
模块M4：采用二分法在区间[pbest,i pnew,i]内计算相应的置信区间上限
Figure GDA00037240704700000610
使得连续概率密度函数pcdf等于给定的置信度pc；
模块M5：令dp＝-dp重复模块M3～模块M4，计算置信区间下限
Figure GDA00037240704700000611
模块M6：令dp＝-dp遍历所有待辨识热物性参数，重复模块M3～模块M5，得到所有给定置信度下热物性参数辨识结果的置信区间。
优选的，所述模块M1包括：
所述结构的一维热传导计算模型通用表达式为：
Figure GDA0003724070470000071
式中ρ表示材料密度；cp表示比热容；λ表示热导率；
Figure GDA0003724070470000078
表示内热源，根据具体物理模型取值，结构没有内部加热源或非相变材料情况下该项等于0；T表示温度；r表示一维方向坐标；m取0、1和2分别对应笛卡尔坐标系、圆柱坐标系和球坐标系；t表示时间。
优选的，所述模块M1还包括：
采用有限体积法对一维结构进行离散，得到如下隐式离散格式：
Figure GDA0003724070470000072
其中，
Figure GDA0003724070470000073
Figure GDA0003724070470000074
Figure GDA0003724070470000075
上式中aw、aP和aE分别为温度前的系数；W、P和E分别表示当前网格点的西侧、当前网格点和当前网格点的东侧；下标i表示第i个网格节点；上标n表示第n时刻；Sc和Sp分别表示内热源
Figure GDA0003724070470000076
线性化的常数项和线性项系数；Δt为时间步长；(·)′表示调和平均参数。变量的上标w和e分别表示变量取第i个网格单元西侧和东侧边界上的值；A表示第i个网格节点处单元横截面积；ri表示第i个网格节点处的坐标；δr表示网格节点到网格单元边界的距离。
左边界条件离散后有：
Figure GDA0003724070470000077
其中，
Figure GDA0003724070470000081
Figure GDA0003724070470000082
上式中上标1表示第一个网格点；其中
Figure GDA0003724070470000083
h、ε、σ、Tr和TAMBIENT分别表示热流密度、对流换热系数、发射率、Stefen-Bolzman常数、对流换热温度和辐射背景温度。
右边界条件离散后有：
Figure GDA0003724070470000084
其中，
Figure GDA0003724070470000085
Figure GDA0003724070470000086
上式中上标M表示第M个网格节点，即最后一个网格节点。
在实际问题中，热导率和比热容是关于温度的函数，其能够采用多项式形式表示。
与现有技术相比，本发明具有如下的有益效果：
1、本发明可以对于不服从正态分布的热物性参数估计其置信区间；
2、本发明可以估计任意置信度下的热物性参数置信区间，置信度更高；
3、本发明估计出的置信区间可以体现出辨识参数不确定度的非线性程度。
附图说明
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
图1为本发明整体流程图；
图2为结构网格控制体示意图；
图3为本发明实施例TC4钛合金材料热物性参数表；
图4为本发明实施例中不同噪声水平下辨识用温度数据曲线；
图5为本发明实施例中热物性参数辨识结果和相对误差；
图6为本发明实施例中不同置信度下估计得到的热物性参数辨识结果置信区间；
图7为本发明实施例中不同置信度下辨识结果与真值的对比图；
图8为飞行器典型舱段示意图；
图9为温度响应的上下边界示意图。
具体实施方式
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。
实施例1：
本发明实施例提供了一种热物性参数辨识结果置信区间估计方法，参照图1所示，具体包括：
步骤S1：对平板试片模型进行网格划分，定义网格的单元属性和材料属性，设置边界条件，基于有限体积法建立结构的一维热传导计算模型。
步骤S2：基于目标温度响应，采用L-M算法对待辨识热物性参数进行优化，辨识得到热物性参数的最优解pbest和最优卡方函数χ(pbest)，本实施例中优化得到热物性参数pbest中，下标best表示最优解。
步骤S3：给定热物性参数置信度pc和增量步长dp，本实施例中给定热物性参数置信度pc中，下标c是置信度confidence的缩写。初始时刻令
Figure GDA0003724070470000091
对第i个热物性参数开始迭代计算置信区间上下限；令
Figure GDA0003724070470000092
将
Figure GDA0003724070470000093
作为初始值并采用L-M算法进行优化得到新的热物性参数
Figure GDA0003724070470000094
和新的卡方函数
Figure GDA0003724070470000095
根据卡方函数
Figure GDA0003724070470000096
和最优卡方函数χ(pbest)计算F分布函数的连续概率密度值pcdf；如果pcdf≤pc，则k＝k+1，并重复步骤3，如果pcdf>pc则停止计算。
步骤S4：采用二分法在区间[pbest,i pnew,i]内计算相应的置信区间上限
Figure GDA0003724070470000101
本实施例中上标up表示置信区间上限，下标i表示向量pbest中第i个元素，使得连续概率密度函数pcdf等于给定的置信度pc，下标cdf是连续分布函数cumulative distribution functions的缩写。
步骤S5：令dp＝-dp重复步骤S3～步骤S4，计算置信区间下限
Figure GDA0003724070470000102
步骤S6：令dp＝-dp遍历所有待辨识热物性参数，重复步骤S3～步骤S5，得到所有给定置信度下热物性参数辨识结果的置信区间。
步骤S7：取热物性参数辨识结果置信区间的中值和两个边界值，开展热防护设计温度响应计算，得到温度响应的上下边界供设计人员参考。
具体地，所述步骤S1包括：
所述结构的一维热传导计算模型通用表达式为：
Figure GDA0003724070470000103
式中ρ、cp、λ、
Figure GDA0003724070470000104
T和r分别为材料密度、比热容、热导率、内热源、温度和一维方向坐标；m取0、1和2分别对应笛卡尔坐标系、圆柱坐标系和球坐标系；t表示时间。对材料进行网格划分，参照图2所示。
采用有限体积法对一维结构进行离散，得到如下隐式离散格式
Figure GDA0003724070470000105
其中，
Figure GDA0003724070470000106
Figure GDA0003724070470000107
Figure GDA0003724070470000108
上式中aw、aP和aE分别为温度前的系数；W、P和E分别表示当前网格点的西侧、当前网格点和当前网格点的东侧；下标i表示第i个网格节点；上标n表示第n时刻；Sc和Sp分别表示内热源
Figure GDA0003724070470000118
线性化的常数项和线性项系数；Δt为时间步长；(·)′表示调和平均参数。变量的上标w和e分别表示变量取第i个网格单元西侧和东侧边界上的值；A表示第i个网格节点处单元横截面积；ri表示第i个网格节点处的坐标；δr表示网格节点到网格单元边界的距离。
左边界条件离散后有：
Figure GDA0003724070470000111
其中，
Figure GDA0003724070470000112
Figure GDA0003724070470000113
上式中上标1表示第一个网格点；其中
Figure GDA0003724070470000114
h、ε、σ、Tr和TAMBIENT分别表示热流密度、对流换热系数、发射率、Stefen-Bolzman常数、对流换热温度和辐射背景温度。
右边界条件离散后有：
Figure GDA0003724070470000115
其中，
Figure GDA0003724070470000116
Figure GDA0003724070470000117
上式中上标M表示第M个网格节点，即最后一个网格节点。
实际问题中，热导率和比热容是关于温度的函数，其可以采用多项式形式表示，以热导率为例，即：
λ＝a0+a1T+…+anTn
假设
Figure GDA0003724070470000121
为实际测量得到的随时间变化的温度向量，则热导率反问题的求解就是使如下卡方函数取最小值
Figure GDA0003724070470000122
其中χ(p)为实值有界函数，p为待辨识参数向量，
Figure GDA0003724070470000123
是温度测量误差，W是测量误差倒数构成的协方差矩阵。
具体地，所述步骤S2包括：
基于目标温度响应，采用L-M算法(Levenberg-Marquardt算法)对带辨识参数向量进行优化。L-M算法通过迭代求解如下方程来优化参数向量p
Figure GDA0003724070470000124
上标j表示迭代步数，上式中
Figure GDA0003724070470000125
为雅各比矩阵，Ω为对角矩阵，μ为阻尼因子，h为增量步长。Ω矩阵的具体形式为
Ω＝diag(JTWJ)
通过对比两个迭代增量步之间目标函数的差值z(p)来判断第j个迭代增量步hj是否合适，其具体表达式为：
Figure GDA0003724070470000126
其中tol为给定的误差，满足上式则接受hj，得到新的p，进入下一步迭代。阻尼因子根据如下方式更新：
Figure GDA0003724070470000127
对于雅各比矩阵J的计算，采用差分方法代替微分方法进行求解，比如
Figure GDA0003724070470000128
计算如上步骤，得到下一个迭代步的待辨识参数pj+1，并计算相应的卡方函数χ(pj +1)，如果小于给定的收敛标准则结束迭代计算，令pbest＝pj+1即为辨识得到的热物性参数；否则更新阻尼因子继续迭代计算。
具体地，所述步骤S3包括：
给定热物性参数置信度pci和增量步长dp。初始时刻令k＝0，i＝1，
Figure GDA0003724070470000131
对第i个热物性参数开始迭代计算置信区间上下限。令
Figure GDA0003724070470000132
并将
Figure GDA0003724070470000133
作为初始条件采用步骤S2辨识热物性参数，得到新的热物性参数
Figure GDA0003724070470000134
和新的卡方函数
Figure GDA0003724070470000135
注意此时L-M优化过程中
Figure GDA0003724070470000136
固定不变。根据卡方函数
Figure GDA0003724070470000137
和最优卡方函数χ(pbest)计算F分布函数的连续概率密度值pcdf，F分布的连续概率密度函数的具体表示式为
Figure GDA0003724070470000138
其中Г为伽马函数；v1和v2为F分布函数在分子和分母上的自由度，本方法中取v1＝1，v2＝nt-np，其中nt为温度响应数据个数，np为热物性参数个数；x为积分上限，本方法取
Figure GDA0003724070470000139
如果pcdf≤pc，则k＝k+1，并重复步骤S3；如果pcdf>pc则停止计算。
具体地，所述步骤S4包括：
对于第i个参数，使用二分法在区间[pbest,i pnew,i]内计算相应的置信区间上限
Figure GDA00037240704700001310
使得连续概率密度函数pcdf等于给定的置信度pc。令p*＝0.5(pbest,i+pnew,i)，
Figure GDA00037240704700001311
且p0best,i＝p*。以
Figure GDA00037240704700001312
为初始条件调用步骤S2优化得到新的
Figure GDA00037240704700001313
并计算相应的卡方函数
Figure GDA00037240704700001314
结合最优卡方函数χ(pbest)计算F分布函数的连续概率密度值pcdf，如果|pcdf-pc|<ε，则得到上限pupperi，否则根据pcdf-pc是否大于0令pbest,i＝p*或者pnew,i＝p*，重复步骤S4直到p*满足要求。
具体地，所述步骤S5包括：
令dp＝-dp重复步骤S3-步骤S4，计算置信区间下限
Figure GDA00037240704700001315
具体地，所述步骤S6包括：
令dp＝-dp遍历所有待辨识热物性参数，重复步骤S3-步骤S5，得到所有给定置信度下热物性参数辨识结果的置信区间。
具体地，所述步骤S7包括：
取热物性参数辨识结果置信区间的中值和上下限，开展热防护设计温度响应计算，得到温度响应的上下边界供设计人员参考。
实施例2：
实施例2是实施例1的变化例。
以TC4材料为例，通过查阅手册得到该材料的热物性参数如图3所示，密度ρ＝4440kg/m3。对一个4mm厚的TC4钛合金试片，对外表面施加150kW/m2热流，计算试片背温。对背温数据增加0均值白噪声，其具体表达式为：
Figure GDA0003724070470000141
其中，ξ为测量误差，η为[-1,1]内正态分布的随机数，
Figure GDA0003724070470000142
为含测量噪声的温度。计算中取ξ＝0.01/0.03/0.05，参照图4所示，给出了辨识用的温度数据。
步骤1：采用有限体积方法对上述钛合金平板厚度方向热传导问题建模，将该问题简化为一维热传导问题，在厚度方向上共划分100个单元进行空间离散。比热容按照表1给出的材料属性赋值，辨识随温度变化的热导率。
步骤2：采用L-M算法分别对不同ξ下的温度数据进行参数优化，辨识得到随温度变化的热导率，其结果参照图5所示。以ξ＝0.05为例，将该状态下的优化解作为辨识结果pbest计算相应的卡方函数χ(pbest)。
步骤3：给定置信度pc分别为68％和95％，给定增量步长dp＝0.1。以热导率第一个参数为例，在初始时刻令pknew＝[6.74+dp 7.41 8.74 9.84 10.46 11.83]，k＝0。以
Figure GDA0003724070470000143
为初始条件，仍然以
Figure GDA0003724070470000144
为辨识数据进行优化，辨识得到相应的参数
Figure GDA0003724070470000145
并重新计算卡方函数
Figure GDA0003724070470000146
注意和步骤2不同的时，此时L-M算法优化过程中
Figure GDA0003724070470000147
不变，即
Figure GDA0003724070470000148
不随优化迭代过程而变化。
根据当前的卡方函数
Figure GDA0003724070470000149
和最优卡方函数χ(pbest)，计算F分布函数连续概率密度函数的积分上限值
Figure GDA00037240704700001410
然后分子和分母自由度参数v1和v2分别取v1＝1和v2＝nt-np。带入如下函数计算pcdf。
Figure GDA0003724070470000151
如果pcdf≤0.95，则k＝k+1，令
Figure GDA0003724070470000152
重复步骤3；如果pcdf>0.95则停止计算。
步骤4：对于第1个参数，使用二分法在区间[pbest,1pnew,1]内计算相应的置信区间上限
Figure GDA0003724070470000153
使得连续概率密度函数pcdf等于给定的置信度pc。令p*＝0.5(pbest,1+pnew,1)，
Figure GDA0003724070470000154
且
Figure GDA0003724070470000155
以
Figure GDA0003724070470000156
为初始条件调用步骤2优化得到新的
Figure GDA0003724070470000157
并计算相应的卡方函数
Figure GDA0003724070470000158
Figure GDA0003724070470000159
结合最优卡方函数χ(pbest)计算F分布函数的连续概率密度值pcdf，如果|pcdf-pc|<ε，则得到上限
Figure GDA00037240704700001510
否则根据pcdf-pc是否大于0令pbest,1＝p*或者pnew,1＝p*，重复步骤4直到p*满足要求。
步骤5：令增量步长dp＝-0.1，重复步骤3～步骤4，计算置信区间下限
Figure GDA00037240704700001511
步骤6：重复步骤3～步骤5，得到给定置信度下6个热导率辨识结果的置信区间，结果参照图6和图7所示，得到的置信区间上下界不关于辨识结果对称，体现了辨识结果不确定性的非线性程度。
步骤7：取热物性参数辨识结果置信区间的中值和上下限，对飞行器典型舱段，参照图8所示，开展热防护设计温度响应计算，得到温度响应的上下边界，参照图9所示，供设计人员参考。本发明实施例提供了一种热物性参数辨识结果置信区间估计方法，可以对于不服从正态分布的热物性参数估计其置信区间；可以估计任意置信度下的热物性参数置信区间，置信度更高；本发明估计出的置信区间可以体现出辨识参数不确定度的非线性程度。
本领域技术人员知道，除了以纯计算机可读程序代码方式实现本发明提供的系统及其各个装置、模块、单元以外，完全可以通过将方法步骤进行逻辑编程来使得本发明提供的系统及其各个装置、模块、单元以逻辑门、开关、专用集成电路、可编程逻辑控制器以及嵌入式微控制器等的形式来实现相同功能。所以，本发明提供的系统及其各项装置、模块、单元可以被认为是一种硬件部件，而对其内包括的用于实现各种功能的装置、模块、单元也可以视为硬件部件内的结构；也可以将用于实现各种功能的装置、模块、单元视为既可以是实现方法的软件模块又可以是硬件部件内的结构。
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本申请的实施例和实施例中的特征可以任意相互组合。

Claims (7)

1.一种热物性参数辨识结果置信区间估计方法，其特征在于，包括：
步骤S1：对平板试片模型进行网格划分，定义网格的单元属性和材料属性，设置边界条件，基于有限体积法建立结构的一维热传导计算模型；
步骤S2：基于目标温度响应，采用L-M算法对待辨识热物性参数进行优化，辨识得到热物性参数的最优解pbest和最优卡方函数χ(pbest)；
步骤S3：给定热物性参数置信度pc和增量步长dp，迭代计算直到F分布函数的连续分布函数值pcdf>pc；
步骤S4：采用二分法在参数最优解p中第i个参数区间[pbest,i pnew,i]内计算相应的置信区间上限
Figure FDA0003724070460000011
使得连续概率密度函数pcdf等于给定的置信度pc；
步骤S5：令dp＝-dp重复步骤S3～步骤S4，计算置信区间下限
Figure FDA0003724070460000012
步骤S6：令dp＝-dp遍历所有待辨识热物性参数，重复步骤S3～步骤S5，得到所有给定置信度下热物性参数辨识结果的置信区间；
步骤S7：取热物性参数辨识结果置信区间的中值和上下限，开展热防护设计温度响应计算，得到温度响应的上下边界供设计人员参考；
所述步骤S1包括：
所述结构的一维热传导计算模型通用表达式为：
Figure FDA0003724070460000013
式中ρ表示材料密度；cp表示比热容；λ表示热导率；
Figure FDA0003724070460000014
表示内热源，根据具体物理模型取值，结构没有内部加热源或非相变材料情况下该项等于0；T表示温度；r表示一维方向坐标；m取0、1和2分别对应笛卡尔坐标系、圆柱坐标系和球坐标系；t表示时间；
所述步骤S2包括：
基于目标温度响应，采用L-M算法对带辨识参数向量进行优化，L-M算法通过迭代求解如下方程来优化参数向量p，
Figure FDA0003724070460000015
上标j表示迭代步数，上式中
Figure FDA0003724070460000016
为雅各比矩阵，Ω为对角矩阵，μ为阻尼因子，h为增量步长，W是测量误差倒数构成的协方差矩阵，
Figure FDA0003724070460000017
为实际测量到的随时间变化的温度向量；Ω矩阵的具体形式为
Ω＝diag(JTWJ)
通过对比两个迭代增量步之间目标函数的差值z(p)来判断第j个迭代增量步hj是否合适，其具体表达式为：
Figure FDA0003724070460000021
其中，tol为给定的误差，满足上式则接受hj，得到新的p，进入下一步迭代；阻尼因子根据如下方式更新：
Figure FDA0003724070460000022
对于雅各比矩阵J的计算，采用差分方法代替微分方法进行求解；计算如上步骤，得到下一个迭代步的待辨识参数pj+1，并计算相应的卡方函数χ(pj+1)，如果小于给定的收敛标准则结束迭代计算，令pbest＝pj+1即为辨识得到的热物性参数；否则更新阻尼因子继续迭代计算。
2.根据权利要求1所述的热物性参数辨识结果置信区间估计方法，其特征在于，所述步骤S1还包括：
采用有限体积法对一维结构进行离散，得到如下隐式离散格式：
Figure FDA0003724070460000023
其中，
Figure FDA0003724070460000024
Figure FDA0003724070460000025
Figure FDA0003724070460000026
上式中aw、aP和aE分别为温度前的系数；W、P和E分别表示当前网格点的西侧、当前网格点和当前网格点的东侧；下标i表示第i个网格节点；上标n表示第n时刻；Sc和Sp分别表示内热源
Figure FDA0003724070460000031
线性化的常数项和线性项系数；Δt为时间步长；(·)′表示调和平均参数；变量的上标w和e分别表示变量取第i个网格单元西侧和东侧边界上的值；A表示第i个网格节点处单元横截面积；ri表示第i个网格节点处的坐标；δr表示网格节点到网格单元边界的距离；
左边界条件离散后有：
Figure FDA0003724070460000032
其中，
Figure FDA0003724070460000033
Figure FDA0003724070460000034
上式中上标1表示第一个网格点；其中
Figure FDA0003724070460000035
h、ε、σ、Tr和TAMBIENT分别表示热流密度、对流换热系数、发射率、Stefen-Bolzman常数、对流换热温度和辐射背景温度；
右边界条件离散后有：
Figure FDA0003724070460000036
其中，
Figure FDA0003724070460000037
Figure FDA0003724070460000038
上式中上标M表示第M个网格节点，即最后一个网格节点；
在实际问题中，热导率和比热容是关于温度的函数，其能够采用多项式形式表示。
3.根据权利要求1所述的热物性参数辨识结果置信区间估计方法，其特征在于，所述步骤S3包括：初始时刻令
Figure FDA0003724070460000041
其中上标k表示迭代步数；对第i个热物性参数开始迭代计算置信区间上下限；
令
Figure FDA0003724070460000042
将
Figure FDA0003724070460000043
作为初始值并采用L-M算法进行优化得到新的热物性参数
Figure FDA0003724070460000044
和新的卡方函数
Figure FDA0003724070460000045
根据卡方函数
Figure FDA0003724070460000046
和最优卡方函数χ(pbest)计算F分布函数的连续概率密度值pcdf；
如果pcdf≤pc，则k＝k+1，并重复步骤S3，如果pcdf>pc则停止计算。
4.根据权利要求3所述的热物性参数辨识结果置信区间估计方法，其特征在于，所述步骤S3还包括：
初始时刻令k＝0，i＝1，
Figure FDA0003724070460000047
对第i个热物性参数开始迭代计算置信区间上下限；
令
Figure FDA0003724070460000048
并将
Figure FDA0003724070460000049
作为初始条件采用步骤S2辨识热物性参数，得到新的热物性参数
Figure FDA00037240704600000410
和新的卡方函数
Figure FDA00037240704600000411
此时L-M优化过程中
Figure FDA00037240704600000412
固定不变，根据卡方函数
Figure FDA00037240704600000413
和最优卡方函数χ(pbest)计算F分布函数的连续概率密度值pcdf，F分布的连续概率密度函数的具体表示式为：
Figure FDA00037240704600000414
其中，Г为伽马函数；v1和v2为F分布函数在分子和分母上的自由度，令v1＝1，v2＝nt-np，其中，nt为温度响应数据个数，np为热物性参数个数；x为积分上限，令
Figure FDA00037240704600000415
若pcdf≤pc，则k＝k+1，并重复步骤S3；若pcdf>pc则停止计算。
5.根据权利要求1所述的热物性参数辨识结果置信区间估计方法，其特征在于，所述步骤S4包括：
对于第i个参数，使用二分法在区间[pbest,i pnew,i]内计算相应的置信区间上限
Figure FDA00037240704600000416
使得连续概率密度函数pcdf等于给定的置信度pc；
令p*＝0.5(pbest,i+pnew,i)，
Figure FDA00037240704600000417
且
Figure FDA00037240704600000418
上标0表示初始迭代步；
以
Figure FDA00037240704600000419
为初始条件调用步骤S2优化得到新的
Figure FDA00037240704600000420
并计算相应的卡方函数
Figure FDA00037240704600000421
结合最优卡方函数χ(pbest)计算F分布函数的连续概率密度值pcdf，若|pcdf-pc|<ε，则得到上限
Figure FDA00037240704600000422
否则根据pcdf-pc是否大于0，令pbest,i＝p*或者pnew,i＝p*，重复步骤S4直到p*满足要求。
6.一种热物性参数辨识结果置信区间估计系统，其特征在于，包括：
模块M1：对平板试片模型进行网格划分，定义网格的单元属性和材料属性，设置边界条件，基于有限体积法建立结构的一维热传导计算模型；
模块M2：基于目标温度响应，采用L-M算法对待辨识热物性参数进行优化，辨识得到热物性参数的最优解pbest和最优卡方函数χ(pbest)；
模块M3：给定热物性参数置信度pc和增量步长dp，迭代计算直到F分布函数的连续分布函数值pcdf>pc；
模块M4：采用二分法在区间[pbest,i pnew,i]内计算相应的置信区间上限
Figure FDA0003724070460000051
使得连续概率密度函数pcdf等于给定的置信度pc；
模块M5：令dp＝-dp重复模块M3～模块M4，计算置信区间下限
Figure FDA0003724070460000052
模块M6：令dp＝-dp遍历所有待辨识热物性参数，重复模块M3～模块M5，得到所有给定置信度下热物性参数辨识结果的置信区间；
所述模块M1包括：
所述结构的一维热传导计算模型通用表达式为：
Figure FDA0003724070460000053
式中ρ表示材料密度；cp表示比热容；λ表示热导率；
Figure FDA0003724070460000054
表示内热源，根据具体物理模型取值，结构没有内部加热源或非相变材料情况下该项等于0；T表示温度；r表示一维方向坐标；m取0、1和2分别对应笛卡尔坐标系、圆柱坐标系和球坐标系；t表示时间；
所述模块M2包括：
基于目标温度响应，采用L-M算法对带辨识参数向量进行优化，L-M算法通过迭代求解如下方程来优化参数向量p，
Figure FDA0003724070460000055
上标j表示迭代步数，上式中
Figure FDA0003724070460000056
为雅各比矩阵，Ω为对角矩阵，μ为阻尼因子，h为增量步长，W是测量误差倒数构成的协方差矩阵，
Figure FDA0003724070460000057
为实际测量到的随时间变化的温度向量；Ω矩阵的具体形式为
Ω＝diag(JTWJ)
通过对比两个迭代增量步之间目标函数的差值z(p)来判断第j个迭代增量步hj是否合适，其具体表达式为：
Figure FDA0003724070460000061
其中，tol为给定的误差，满足上式则接受hj，得到新的p，进入下一步迭代；阻尼因子根据如下方式更新：
Figure FDA0003724070460000062
对于雅各比矩阵J的计算，采用差分方法代替微分方法进行求解；计算如上步骤，得到下一个迭代步的待辨识参数pj+1，并计算相应的卡方函数χ(pj+1)，如果小于给定的收敛标准则结束迭代计算，令pbest＝pj+1即为辨识得到的热物性参数；否则更新阻尼因子继续迭代计算。
7.根据权利要求6所述的热物性参数辨识结果置信区间估计系统，其特征在于，所述模块M1还包括：
采用有限体积法对一维结构进行离散，得到如下隐式离散格式：
Figure FDA0003724070460000063
其中，
Figure FDA0003724070460000064
Figure FDA0003724070460000065
Figure FDA0003724070460000066
上式中aw、aP和aE分别为温度前的系数；W、P和E分别表示当前网格点的西侧、当前网格点和当前网格点的东侧；下标i表示第i个网格节点；上标n表示第n时刻；Sc和Sp分别表示内热源
Figure FDA0003724070460000067
线性化的常数项和线性项系数；Δt为时间步长；(·)′表示调和平均参数；变量的上标w和e分别表示变量取第i个网格单元西侧和东侧边界上的值；A表示第i个网格节点处单元横截面积；ri表示第i个网格节点处的坐标；δr表示网格节点到网格单元边界的距离；
左边界条件离散后有：
Figure FDA0003724070460000071
其中，
Figure FDA0003724070460000072
Figure FDA0003724070460000073
上式中上标1表示第一个网格点；其中
Figure FDA0003724070460000074
h、ε、σ、Tr和TAMBIENT分别表示热流密度、对流换热系数、发射率、Stefen-Bolzman常数、对流换热温度和辐射背景温度；
右边界条件离散后有：
Figure FDA0003724070460000075
其中，
Figure FDA0003724070460000076
Figure FDA0003724070460000077
上式中上标M表示第M个网格节点，即最后一个网格节点；
在实际问题中，热导率和比热容是关于温度的函数，其能够采用多项式形式表示。
CN202110522174.4A 2021-05-13 2021-05-13 一种热物性参数辨识结果置信区间估计方法及系统 Active CN113361155B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110522174.4A CN113361155B (zh) 2021-05-13 2021-05-13 一种热物性参数辨识结果置信区间估计方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110522174.4A CN113361155B (zh) 2021-05-13 2021-05-13 一种热物性参数辨识结果置信区间估计方法及系统

Publications (2)

Publication Number Publication Date
CN113361155A CN113361155A (zh) 2021-09-07
CN113361155B true CN113361155B (zh) 2022-09-13

Family

ID=77526244

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110522174.4A Active CN113361155B (zh) 2021-05-13 2021-05-13 一种热物性参数辨识结果置信区间估计方法及系统

Country Status (1)

Country Link
CN (1) CN113361155B (zh)

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105868465A (zh) * 2016-03-28 2016-08-17 大连理工大学 一种随温度变化的热导率辨识的修正lm方法
CN105956344A (zh) * 2016-06-23 2016-09-21 大连理工大学 复杂结构多维瞬态非线性热传导反问题的简易快速求解方法
CN109684717B (zh) * 2018-12-23 2023-04-07 贵州大学 基于量纲分析的油炒烹饪中表面换热系数的预测方法
CN109918613B (zh) * 2019-03-06 2019-10-25 哈尔滨工业大学 一种基于阶跃响应标定的热流辨识方法
CN109900734B (zh) * 2019-04-04 2021-11-19 哈尔滨工业大学 方形锂离子电池内部热物性参数辨识方法
CN110261427B (zh) * 2019-07-04 2020-06-09 西安交通大学 基于共轭梯度法的多层复合材料导热系数测量方法
CN111125913B (zh) * 2019-12-25 2023-11-03 东北大学 一种加热炉总括热吸收率的辨识方法及装置
CN112329304B (zh) * 2020-11-04 2022-07-01 北京航空航天大学 一种连续结构动态载荷区间辨识方法

Also Published As

Publication number Publication date
CN113361155A (zh) 2021-09-07

Similar Documents

Publication Publication Date Title
Wang et al. Machine learning aided stochastic structural free vibration analysis for functionally graded bar-type structures
Baumann et al. A discontinuous hp finite element method for the Euler and Navier–Stokes equations
Choi et al. Reliability-Based Structural Optimization
Wang et al. Structural reliability optimization using an efficient safety index calculation procedure
US20200082036A1 (en) Simulation method, simulation method by mbd program, numerical analysis apparatus, numerical analysis system for mbd, numerical analysis program, and mbd program
Widhalm et al. Investigation on adjoint based gradient computations for realistic 3D aero-optimization
Yu et al. Time-variant reliability analysis via approximation of the first-crossing PDF
Li et al. Non-intrusive coupling of a 3-D Generalized Finite Element Method and Abaqus for the multiscale analysis of localized defects and structural features
Ahmed et al. Surrogate-based aerodynamic design optimization: Use of surrogates in aerodynamic design optimization
CN113987691A (zh) 激波失稳的高精度混合计算方法、装置、设备和存储介质
Li et al. Spectral stochastic isogeometric analysis for linear stability analysis of plate
Ouyang et al. A novel dynamic model updating method for composite laminate structures considering non-probabilistic uncertainties and correlations
DeCarlo et al. Bayesian calibration of aerothermal models for hypersonic air vehicles
CN113361155B (zh) 一种热物性参数辨识结果置信区间估计方法及系统
Waligura et al. Investigation of spalart-allmaras turbulence model modifications for hypersonic flows utilizing output-based grid adaptation
CN114091297A (zh) 用于热障涂层材料破坏预测的热力耦合近场动力学方法
CN112580855A (zh) 基于自适应变异pso-bp神经网络的电缆群稳态温升预测方法
CN109598059B (zh) 一种基于代理模型的热防护系统优化设计方法及设计系统
CN110096798B (zh) 一种多状态有限元模型修正的方法
Morse et al. Manufacturing cost and reliability‐based shape optimization of plate structures
Zhang et al. A novel robust aerodynamic optimization technique coupled with adjoint solvers and polynomial chaos expansion
Piqueras et al. Calculation of linear conductances for thermal lumped models by means of the CMF method
Zhuang et al. Implicit differentiation-based reliability analysis for shallow shell structures with the Boundary Element Method
Ben Ameur et al. Physics-based mesh fitting algorithms for hypersonic flows simulations
Yildiz et al. Uncertainty quantification of aeroelastic systems using active learning gaussian process

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant