CN111783236A - 基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法 - Google Patents

基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法 Download PDF

Info

Publication number
CN111783236A
CN111783236A CN202010407601.XA CN202010407601A CN111783236A CN 111783236 A CN111783236 A CN 111783236A CN 202010407601 A CN202010407601 A CN 202010407601A CN 111783236 A CN111783236 A CN 111783236A
Authority
CN
China
Prior art keywords
sample
failure
function
samples
value
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.)
Granted
Application number
CN202010407601.XA
Other languages
English (en)
Other versions
CN111783236B (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.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
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 Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202010407601.XA priority Critical patent/CN111783236B/zh
Publication of CN111783236A publication Critical patent/CN111783236A/zh
Application granted granted Critical
Publication of CN111783236B publication Critical patent/CN111783236B/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/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • 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
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (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)
  • Complex Calculations (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

本发明提供一种基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法，包括：训练Kriging模型；将得到的每一层失效事件对应的条件失效概率的估算值做乘法运算，得到最终的涡轮机匣失效概率的估计值
Figure DDA0002491911290000011
得到可靠性局部灵敏度指标和可靠性全局灵敏度指标。本发明提供的基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法具有以下优点：本发明提供一种基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法，是一种有效可行的同时求解可靠性局部灵敏度指标和可靠性全局灵敏度指标的方法。可以克服现有方法所需样本量较大及耗时较长的问题，本发明提高涡轮机匣灵敏度分析的效率。

Description

基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法
技术领域
本发明属于灵敏度分析技术领域，具体涉及一种基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法。
背景技术
涡轮机匣是航空发动机的重要零件之一，是整个发动机的基座，其外形结构复杂，不同的发动机、发动机不同部位，其机匣形状各不相同，机匣零件的功能决定了机匣的形状，但其基本特征是圆筒形或圆锥形的壳体和支板组成的构件。由于机匣零件设计难度大、周期长，在整个发动机的设计中，机匣的设计占相当大的比重，提高机匣的设计效率对压缩发动机整机的设计周期具有重要意义。
为了查出对机匣失效概率影响较大的设计参数，需要对机匣进行灵敏度分析。传统的灵敏度分析方法需要多次有限元模型，在实际应用中需要过多的分析时间，且对于小失效概率问题需要更多的样本量，具有灵敏度分析效率低的问题，从而制约了其在工程实践中的应用。
发明内容
针对现有技术存在的缺陷，本发明提供一种基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法，可有效解决上述问题。
本发明采用的技术方案如下：
本发明提供一种基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法，包括以下步骤：
步骤1，建立涡轮机匣模型；对所述涡轮机匣模型进行分析，确定影响涡轮机匣灵敏度的n维随机向量x＝{x1,x2,…,xn}T；
步骤2，令k＝1；
步骤3，根据n维随机向量x的联合概率密度函数fX(x)，产生Nk个样本构建样本集合
Figure BDA0002491911270000021
其中，
Figure BDA0002491911270000022
代表样本集合Sk中的第i个样本，为n维随机向量x＝{x1,x2,…,xn}T中各维自变量取具体数值后的n维数组；
步骤4，从样本集合Sk中随机选择T1个样本，形成训练集合
Figure BDA0002491911270000023
调用涡轮机匣失效的实际功能函数y＝g(x)，分别计算训练集合T中每个样本的实际功能函数值，得到实际功能函数值集合
Figure BDA0002491911270000024
将实际功能函数值集合gt中的各个元素从大到小排序，排序结果为
Figure BDA0002491911270000025
取第
Figure BDA0002491911270000026
个实际功能函数值作为第k层失效事件Fk＝{x:g(x)≤bk}的临界值bk，即
Figure BDA0002491911270000027
其中，
Figure BDA0002491911270000028
表示向上取整符号，p0为预设已知值；
步骤5，由训练集合
Figure BDA0002491911270000029
以及相应的实际功能函数值集合
Figure BDA00024919112700000210
构建得到Kriging模型gK(x)如下：
Figure BDA00024919112700000211
其中：
f(x)＝[f1(x),f2(x),...,fp(x)]T是n维随机向量x的p个基函数形成的集合，用于提供设计空间内的全局近似模型；p是基函数的个数；
β＝[β1,β2,...,βp]为回归函数的p个系数值形成的集合；每个回归函数的系数值通过训练集合
Figure BDA00024919112700000212
以及相应的实际功能函数值集合
Figure BDA00024919112700000213
得到；
z(x)代表随机过程；
步骤6，对于样本集合Sk中任意的样本
Figure BDA00024919112700000214
代入Kriging模型gK(x)，得到对应的一个代理功能函数值
Figure BDA00024919112700000215
每个代理功能函数值
Figure BDA00024919112700000216
均服从代理功能函数均值为
Figure BDA00024919112700000217
标准差为
Figure BDA00024919112700000218
的高斯分布：
Figure BDA00024919112700000219
因此，对于任意的样本
Figure BDA00024919112700000220
可得到代理功能函数均值为
Figure BDA00024919112700000221
代理功能函数标准差
Figure BDA0002491911270000031
步骤7，根据下式，计算样本集合Sk中每个样本
Figure BDA0002491911270000032
对应的U学习函数值
Figure BDA0002491911270000033
Figure BDA0002491911270000034
步骤8，因此，由于样本集合Sk中共有Nk个样本，共计算得到Nk个U学习函数值，判断每个U学习函数值是否均大于等于2，即：判断下式是否成立：
Figure BDA0002491911270000035
如果成立，则代表当前Kriging模型已收敛，停止更新当前Kriging模型，得到训练好的Kriging模型，执行步骤9；
否则，根据下式从样本集合Sk选取新的样本，也就是选取具有最小U学习函数值的样本，表示为
Figure BDA0002491911270000036
Figure BDA0002491911270000037
计算样本
Figure BDA0002491911270000038
的实际功能函数值
Figure BDA0002491911270000039
将
Figure BDA00024919112700000310
和
Figure BDA00024919112700000311
分别加入训练集合T和实际功能函数值集合gt中，由此得到更新后的训练集合T和更新后的实际功能函数值集合gt，再返回步骤5；
步骤9，利用训练好的Kriging模型，计算样本集合Sk中每个样本
Figure BDA00024919112700000312
对应的代理功能函数均值
Figure BDA00024919112700000313
根据下述指示函数
Figure BDA00024919112700000314
从样本集合Sk中挑选可使第k层失效事件Fk＝{x:g(x)≤bk}发生的样本：
Figure BDA00024919112700000315
若
Figure BDA00024919112700000316
则样本
Figure BDA00024919112700000317
可使第k层失效事件Fk＝{x:g(x)≤bk}发生；
样本集合Sk中，所有使第k层失效事件Fk＝{x:g(x)≤bk}发生的样本组成的集合记为失效样本集合
Figure BDA0002491911270000041
其中，Mk为使第k层失效事件Fk＝{x:g(x)≤bk}发生的失效样本总数；
根据下式计算条件失效概率Pk的估算值
Figure BDA0002491911270000042
Figure BDA0002491911270000043
其中：
Figure BDA0002491911270000044
表示第k层失效事件Fk发生概率的估算值；
Figure BDA0002491911270000045
表示第k-1层失效事件Fk-1发生的条件下，第k层失效事件Fk发生概率的估算值；
P{Fk-1}表示第k-1层失效事件Fk-1发生概率；
Figure BDA0002491911270000046
表示重要抽样密度函数；
根据下式计算条件失效概率Pk的估算值
Figure BDA0002491911270000047
对
Figure BDA0002491911270000048
的偏导数
Figure BDA0002491911270000049
Figure BDA00024919112700000410
其中：
Figure BDA00024919112700000411
代表n维随机向量x＝{x1,x2,…,xn}T中的第r维自变量xr的第s个分布参数；
Pq代表第q层失效事件Fq发生的概率，q＝1,2,...,k-1；
步骤10，判断bk是否等于0，如果bk＝0，则循环过程结束，执行步骤15；否则，执行步骤11；
步骤11，令k＝k+1，在失效样本集合RSk-1中，选出联合概率密度函数值最大的样本，作为重要抽样密度函数
Figure BDA00024919112700000412
的抽样中心；
然后，根据n维随机向量x的重要抽样密度函数
Figure BDA0002491911270000051
以联合概率密度函数值最大的样本为抽样中心，产生Nk个样本构建样本集合
Figure BDA0002491911270000052
步骤12，利用训练好的Kriging模型，计算样本集合Sk中每个样本
Figure BDA0002491911270000053
对应的代理功能函数值
Figure BDA0002491911270000054
和代理功能函数均值
Figure BDA0002491911270000055
根据下述指示函数
Figure BDA0002491911270000056
从样本集合Sk中挑选可使第k层失效事件Fk＝{x:g(x)≤bk}发生的样本：
Figure BDA0002491911270000057
样本集合Sk中，所有使第k层失效事件Fk＝{x:g(x)≤bk}发生的样本组成的集合记为失效样本集合
Figure BDA0002491911270000058
其中，Mk为使第k层失效事件Fk＝{x:g(x)≤bk}发生的失效样本总数；
步骤13，利用训练好的Kriging模型，计算失效样本集合RSk中每个样本对应的代理功能函数均值；然后，将失效样本集合RSk中每个样本对应的代理功能函数均值从大到小排序，
Figure BDA0002491911270000059
取第
Figure BDA00024919112700000510
个代理功能函数均值作为第k层失效事件Fk＝{x:g(x)≤bk}的临界值bk，即
Figure BDA00024919112700000511
步骤14，判断bk是否小于0；如果小于，则令bk＝0，然后返回步骤7；如果不小于，直接返回步骤7；
步骤15，将至此得到的每一层失效事件对应的条件失效概率的估算值做乘法运算，得到最终的涡轮机匣失效概率的估计值
Figure BDA00024919112700000512
具体的，令m＝k，根据下式计算涡轮机匣失效概率的估计值
Figure BDA00024919112700000513
Figure BDA00024919112700000514
通过下式计算得到最终的n维随机向量x＝{x1,x2,…,xn}T中的第r维自变量xr的第s个分布参数的可靠性局部灵敏度指标
Figure BDA0002491911270000061
Figure BDA0002491911270000062
优选的，影响涡轮机匣灵敏度的n维随机向量中，各分量包括高温区域温度、流动冷却区域温度、冲击冷却区域温度、前段机匣长度及机匣内半径。
优选的，还包括：
步骤16，最后一层失效事件即为第m层失效事件Fm，使第m层失效事件Fm发生的样本组成的失效样本集合为RSm，设失效样本集合RSm中失效样本总数为NC；
将失效样本集合RSm中的所有失效样本放入矩阵C，则矩阵C中的样本表示为
Figure BDA0002491911270000063
其中，b＝1,2,...,NC；
步骤17，初始化
Figure BDA0002491911270000064
并令β＝1；
步骤18，构建并计算根据重要抽样密度函数产生的样本
Figure BDA0002491911270000065
向目标分布中第β+1状态zβ+1转化的接受概率
Figure BDA0002491911270000066
为：
Figure BDA0002491911270000067
其中：
Figure BDA0002491911270000068
为样本
Figure BDA0002491911270000069
的联合概率密度函数；
Figure BDA00024919112700000610
为第β状态zβ的重要抽样密度函数；
fX(zβ)为第β状态zβ的联合概率密度函数；
Figure BDA00024919112700000611
为样本
Figure BDA00024919112700000612
的重要抽样密度函数；
步骤19，根据以下判断式，判断zβ+1是否接受
Figure BDA00024919112700000613
Figure BDA00024919112700000614
其中，random[0,1]表示0与1之间的随机数；
即：如果
Figure BDA00024919112700000615
则将
Figure BDA00024919112700000616
赋值给zβ+1；否则，将zβ赋值给zβ+1；
步骤20，令β＝β+1，判断β是否大于NC，如果大于，则执行步骤21；否则，执行步骤18；
步骤21，由此得到NC个zβ，分别为
Figure BDA0002491911270000071
将
Figure BDA0002491911270000072
放入到矩阵B，矩阵B中样本服从输入变量的条件概率密度函数fX(x|F)；
步骤22，根据矩阵B中的样本估计失效域内xr条件概率密度函数
Figure BDA0002491911270000073
再结合步骤15得到的涡轮机匣失效概率的估计值
Figure BDA0002491911270000074
以及第r维自变量xr的联合概率密度函数
Figure BDA0002491911270000075
代入下式得到输入变量xr的可靠性全局灵敏度指标Sr：
Figure BDA0002491911270000076
其中：
E[·]表示期望运算。
本发明提供的基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法具有以下优点：
本发明提供一种基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法，是一种有效可行的同时求解可靠性局部灵敏度指标和可靠性全局灵敏度指标的方法。可以克服现有方法所需样本量较大及耗时较长的问题，本发明提高涡轮机匣灵敏度分析的效率。
附图说明
图1为本发明提供的基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法的流程示意图。
图2是本发明应用子集模拟引入中间事件的示意图；
图3是本发明涡轮机匣的平面简图；
图4是本发明涡轮机匣的受力示意图。
具体实施方式
为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
本发明提供一种基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法，是一种有效可行的同时求解可靠性局部灵敏度指标和可靠性全局灵敏度指标的方法。该方法旨在克服现有方法所需样本量较大及耗时较长的问题，本发明提高灵敏度分析的效率。
本发明采用的技术方案主要包含两个部分，分别为：(一)利用子集模拟法产生中间失效事件和分层抽样的样本、利用Kriging模型自适应的选点和代理功能函数；(二)基于贝叶斯公式求解失效概率的可靠性全局灵敏度指标。
参考图1，本发明提供的一种基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法，包括以下步骤：
第一部分：利用子集模拟法和自适应Kriging构造功能函数的代理模型及求解可靠性局部灵敏度指标，其中，各子集的样本由重要抽样法得到的抽样密度函数产生，Kriging模型由U学习函数进行加点并判断是否满足收敛条件。具体步骤包括步骤1-步骤14：
步骤1，建立涡轮机匣模型；对所述涡轮机匣模型进行分析，确定影响涡轮机匣灵敏度的n维随机向量x＝{x1,x2,…,xn}T；
步骤2，令k＝1；
步骤3，根据n维随机向量x的联合概率密度函数fX(x)，产生Nk个样本构建样本集合
Figure BDA0002491911270000081
其中，
Figure BDA0002491911270000082
代表样本集合Sk中的第i个样本，为n维随机向量x＝{x1,x2,…,xn}T中各维自变量取具体数值后的n维数组；
步骤4，从样本集合Sk中随机选择T1个样本，形成训练集合
Figure BDA0002491911270000083
调用涡轮机匣失效的实际功能函数y＝g(x)，分别计算训练集合T中每个样本的实际功能函数值，得到实际功能函数值集合
Figure BDA0002491911270000084
将实际功能函数值集合gt中的各个元素从大到小排序，排序结果为
Figure BDA0002491911270000085
取第
Figure BDA0002491911270000086
个实际功能函数值作为第k层失效事件Fk＝{x:g(x)≤bk}的临界值bk，即
Figure BDA0002491911270000091
其中，
Figure BDA0002491911270000092
表示向上取整符号，p0为预设已知值；
步骤5，由训练集合
Figure BDA0002491911270000093
以及相应的实际功能函数值集合
Figure BDA0002491911270000094
构建得到Kriging模型gK(x)如下：
Figure BDA0002491911270000095
其中：
f(x)＝[f1(x),f2(x),...,fp(x)]T是n维随机向量x的p个基函数形成的集合，用于提供设计空间内的全局近似模型；p是基函数的个数；
β＝[β1,β2,...,βp]为回归函数的p个系数值形成的集合；每个回归函数的系数值通过训练集合
Figure BDA0002491911270000096
以及相应的实际功能函数值集合
Figure BDA0002491911270000097
得到；
z(x)代表随机过程；
步骤6，对于样本集合Sk中任意的样本
Figure BDA0002491911270000098
代入Kriging模型gK(x)，得到对应的一个代理功能函数值
Figure BDA0002491911270000099
每个代理功能函数值
Figure BDA00024919112700000910
均服从代理功能函数均值为
Figure BDA00024919112700000911
标准差为
Figure BDA00024919112700000912
的高斯分布：
Figure BDA00024919112700000913
因此，对于任意的样本
Figure BDA00024919112700000914
可得到代理功能函数均值为
Figure BDA00024919112700000915
代理功能函数标准差
Figure BDA00024919112700000916
步骤7，根据下式，计算样本集合Sk中每个样本
Figure BDA00024919112700000917
对应的U学习函数值
Figure BDA00024919112700000918
Figure BDA00024919112700000919
步骤8，因此，由于样本集合Sk中共有Nk个样本，共计算得到Nk个U学习函数值，判断每个U学习函数值是否均大于等于2，即：判断下式是否成立：
Figure BDA00024919112700000920
如果成立，则代表当前Kriging模型已收敛，停止更新当前Kriging模型，得到训练好的Kriging模型，执行步骤9；
否则，根据下式从样本集合Sk选取新的样本，也就是选取具有最小U学习函数值的样本，表示为
Figure BDA0002491911270000101
Figure BDA0002491911270000102
计算样本
Figure BDA0002491911270000103
的实际功能函数值
Figure BDA0002491911270000104
将
Figure BDA0002491911270000105
和
Figure BDA0002491911270000106
分别加入训练集合T和实际功能函数值集合gt中，由此得到更新后的训练集合T和更新后的实际功能函数值集合gt，再返回步骤5；
步骤9，利用训练好的Kriging模型，计算样本集合Sk中每个样本
Figure BDA0002491911270000107
对应的代理功能函数均值
Figure BDA0002491911270000108
根据下述指示函数
Figure BDA0002491911270000109
从样本集合Sk中挑选可使第k层失效事件Fk＝{x:g(x)≤bk}发生的样本：
Figure BDA00024919112700001010
若
Figure BDA00024919112700001011
则样本
Figure BDA00024919112700001012
可使第k层失效事件Fk＝{x:g(x)≤bk}发生；
样本集合Sk中，所有使第k层失效事件Fk＝{x:g(x)≤bk}发生的样本组成的集合记为失效样本集合
Figure BDA00024919112700001013
其中，Mk为使第k层失效事件Fk＝{x:g(x)≤bk}发生的失效样本总数；
根据下式计算条件失效概率Pk的估算值
Figure BDA00024919112700001014
Figure BDA00024919112700001015
其中：
Figure BDA00024919112700001016
表示第k层失效事件Fk发生概率的估算值；
Figure BDA0002491911270000111
表示第k-1层失效事件Fk-1发生的条件下，第k层失效事件Fk发生概率的估算值；
P{Fk-1}表示第k-1层失效事件Fk-1发生概率；
Figure BDA0002491911270000112
表示重要抽样密度函数；
根据下式计算条件失效概率Pk的估算值
Figure BDA0002491911270000113
对
Figure BDA0002491911270000114
的偏导数
Figure BDA0002491911270000115
Figure BDA0002491911270000116
其中：
Figure BDA0002491911270000117
代表n维随机向量x＝{x1,x2,…,xn}T中的第r维自变量xr的第s个分布参数；
Pq代表第q层失效事件Fq发生的概率，q＝1,2,...,k-1；
步骤10，判断bk是否等于0，如果bk＝0，则循环过程结束，执行步骤15；否则，执行步骤11；
步骤11，令k＝k+1，在失效样本集合RSk-1中，选出联合概率密度函数值最大的样本，作为重要抽样密度函数
Figure BDA0002491911270000118
的抽样中心；
然后，根据n维随机向量x的重要抽样密度函数
Figure BDA0002491911270000119
以联合概率密度函数值最大的样本为抽样中心，产生Nk个样本构建样本集合
Figure BDA00024919112700001110
步骤12，利用训练好的Kriging模型，计算样本集合Sk中每个样本
Figure BDA00024919112700001111
对应的代理功能函数值
Figure BDA00024919112700001112
和代理功能函数均值
Figure BDA00024919112700001113
根据下述指示函数
Figure BDA00024919112700001114
从样本集合Sk中挑选可使第k层失效事件Fk＝{x:g(x)≤bk}发生的样本：
Figure BDA0002491911270000121
样本集合Sk中，所有使第k层失效事件Fk＝{x:g(x)≤bk}发生的样本组成的集合记为失效样本集合
Figure BDA0002491911270000122
其中，Mk为使第k层失效事件Fk＝{x:g(x)≤bk}发生的失效样本总数；
步骤13，利用训练好的Kriging模型，计算失效样本集合RSk中每个样本对应的代理功能函数均值；然后，将失效样本集合RSk中每个样本对应的代理功能函数均值从大到小排序，
Figure BDA0002491911270000123
取第
Figure BDA0002491911270000124
个代理功能函数均值作为第k层失效事件Fk＝{x:g(x)≤bk}的临界值bk，即
Figure BDA0002491911270000125
步骤14，判断bk是否小于0；如果小于，则令bk＝0，然后返回步骤7；如果不小于，直接返回步骤7；
步骤15，将至此得到的每一层失效事件对应的条件失效概率的估算值做乘法运算，得到最终的涡轮机匣失效概率的估计值
Figure BDA0002491911270000126
具体的，令m＝k，根据下式计算涡轮机匣失效概率的估计值
Figure BDA0002491911270000127
Figure BDA0002491911270000128
通过下式计算得到最终的n维随机向量x＝{x1,x2,…,xn}T中的第r维自变量xr的第s个分布参数的可靠性局部灵敏度指标
Figure BDA0002491911270000129
Figure BDA00024919112700001210
第二部分：基于Bayes公式
Figure BDA00024919112700001211
能建立无条件失效概率P{F}与条件失效概率P{F|Xi}之间的关系，由此可以进一步得到基于失效概率的可靠性全局灵敏度指标Si，但在步骤一中得到的失效样本均服从条件重要抽样概率密度函数
Figure BDA0002491911270000131
故为了求解Si需先利用Metropolis-Hastings准则将其转换为服从fX(x|F)的样本。具体包括以下步骤16-步骤22。
步骤16，最后一层失效事件即为第m层失效事件Fm，使第m层失效事件Fm发生的样本组成的失效样本集合为RSm，设失效样本集合RSm中失效样本总数为NC；
由于使最后一层失效事件Fm发生的失效样本较多，故利用最后一层的失效样本进行转换。
将失效样本集合RSm中的所有失效样本放入矩阵C，则矩阵C中的样本表示为
Figure BDA0002491911270000132
其中，b＝1,2,...,NC；
步骤17，初始化
Figure BDA0002491911270000133
并令β＝1；
步骤18，构建并计算根据重要抽样密度函数产生的样本
Figure BDA0002491911270000134
向目标分布中第β+1状态zβ+1转化的接受概率
Figure BDA0002491911270000135
为：
Figure BDA0002491911270000136
其中：
Figure BDA0002491911270000137
为样本
Figure BDA0002491911270000138
的联合概率密度函数；
Figure BDA0002491911270000139
为第β状态zβ的重要抽样密度函数；
fX(zβ)为第β状态zβ的联合概率密度函数；
Figure BDA00024919112700001310
为样本
Figure BDA00024919112700001311
的重要抽样密度函数；
步骤19，根据以下判断式，判断zβ+1是否接受
Figure BDA00024919112700001312
Figure BDA00024919112700001313
其中，random[0,1]表示0与1之间的随机数；
即：如果
Figure BDA00024919112700001314
则将
Figure BDA00024919112700001315
赋值给zβ+1；否则，将zβ赋值给zβ+1；
步骤20，令β＝β+1，判断β是否大于NC，如果大于，则执行步骤21；否则，执行步骤18；
步骤21，由此得到NC个zβ，分别为
Figure BDA0002491911270000141
将
Figure BDA0002491911270000142
放入到矩阵B，矩阵B中样本服从输入变量的条件概率密度函数fX(x|F)；
步骤22，根据矩阵B中的样本估计失效域内xr条件概率密度函数
Figure BDA0002491911270000143
再结合步骤15得到的涡轮机匣失效概率的估计值
Figure BDA0002491911270000144
以及第r维自变量xr的联合概率密度函数
Figure BDA0002491911270000145
代入下式得到输入变量xr的可靠性全局灵敏度指标Sr：
Figure BDA0002491911270000146
其中：
E[·]表示期望运算。
下面介绍一种具体实施例：
针对某型发动机涡轮机匣进行灵敏度分析，首先需要建立参数化的涡轮机匣部件模型。为了简化原始机匣模型，本发明忽略安装边上的孔及支座的影响，将该机匣简化为轴对称的模型，机匣外表面的加强筋不予考虑，机匣的安装边的安装孔也简化。简化后的机匣由外机匣组成，确定性设计参数见表1。假定6个随机输入变量分别为高温区域温度、流动冷却区域温度、冲击冷却区域温度、前段机匣长度及机匣内半径，随机输入变量的分布参数表见表2。之后对其利用ANSYS软件进行有限元建模，并进行参数化处理可供MATLAB调用。
表1涡轮机匣参数表
Figure BDA0002491911270000147
表2随机输入变量的分布参数
Figure BDA0002491911270000151
然后，按本发明提供的基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法，对涡轮机匣灵敏度进行分析，结果见表3：
表3可靠性灵敏度指标的计算结果
Figure BDA0002491911270000152
因此，对于6维随机向量x中的每一维分量，分别为高温区域温度T1(℃)、流动冷却区域温度T2(℃)、冲击冷却区域温度T3(℃)、压强p(Mpa)，前段机匣长度l1和机匣内半径r1(mm)，对于每一维分量，基于步骤15的可靠性局部灵敏度指标计算公式，分别取分布参数为均值和标准差，可分别计算得到对均值的可靠性局部灵敏度指标，以及对标准差的可靠性局部灵敏度指标；根据步骤22的公式，可得到每一维分量对应的可靠性全局灵敏度指标。
通过对表3进行分析，可以看出各维分量对均值的可靠性局部灵敏度指标、对标准差的可靠性局部灵敏度指标以及可靠性全局灵敏度指标的排序是一致的，即：对均值的可靠性局部灵敏度指标，高温区域温度T1(℃)、流动冷却区域温度T2(℃)、冲击冷却区域温度T3(℃)、压强p(Mpa)，前段机匣长度l1和机匣内半径r1(mm)的值逐渐增大；同样，对标准差的可靠性局部灵敏度指标，这6维分量的值逐渐增大；同样，对可靠性全局灵敏度指标，这6维分量的值逐渐增大。并且，机匣内半径r1(mm)的各个指标值均为最大，说明其对涡轮机匣部失效有最重要的影响，因此，机匣内半径r1(mm)为对涡轮机匣可靠性最为敏感的设计参数，在对涡轮机匣的可靠性优化设计过程中，应着重考虑机匣内半径的尺寸的优化。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (3)

1.一种基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法，其特征在于，包括以下步骤：
步骤1，建立涡轮机匣模型；对所述涡轮机匣模型进行分析，确定影响涡轮机匣灵敏度的n维随机向量x＝{x1,x2,…,xn}T；
步骤2，令k＝1；
步骤3，根据n维随机向量x的联合概率密度函数fX(x)，产生Nk个样本构建样本集合
Figure FDA0002491911260000011
其中，
Figure FDA0002491911260000012
代表样本集合Sk中的第i个样本，为n维随机向量x＝{x1,x2,…,xn}T中各维自变量取具体数值后的n维数组；
步骤4，从样本集合Sk中随机选择T1个样本，形成训练集合
Figure FDA0002491911260000013
调用涡轮机匣失效的实际功能函数y＝g(x)，分别计算训练集合T中每个样本的实际功能函数值，得到实际功能函数值集合
Figure FDA0002491911260000014
将实际功能函数值集合gt中的各个元素从大到小排序，排序结果为
Figure FDA0002491911260000015
取第
Figure FDA0002491911260000016
个实际功能函数值作为第k层失效事件Fk＝{x:g(x)≤bk}的临界值bk，即
Figure FDA0002491911260000017
其中，
Figure FDA0002491911260000018
表示向上取整符号，p0为预设已知值；
步骤5，由训练集合
Figure FDA0002491911260000019
以及相应的实际功能函数值集合
Figure FDA00024919112600000110
构建得到Kriging模型gK(x)如下：
Figure FDA00024919112600000111
其中：
f(x)＝[f1(x),f2(x),...,fp(x)]T是n维随机向量x的p个基函数形成的集合，用于提供设计空间内的全局近似模型；p是基函数的个数；
β＝[β1,β2,...,βp]为回归函数的p个系数值形成的集合；每个回归函数的系数值通过训练集合
Figure FDA00024919112600000112
以及相应的实际功能函数值集合
Figure FDA00024919112600000113
得到；
z(x)代表随机过程；
步骤6，对于样本集合Sk中任意的样本
Figure FDA0002491911260000021
代入Kriging模型gK(x)，得到对应的一个代理功能函数值
Figure FDA0002491911260000022
每个代理功能函数值
Figure FDA0002491911260000023
均服从代理功能函数均值为
Figure FDA0002491911260000024
标准差为
Figure FDA0002491911260000025
的高斯分布：
Figure FDA0002491911260000026
因此，对于任意的样本
Figure FDA0002491911260000027
可得到代理功能函数均值为
Figure FDA0002491911260000028
代理功能函数标准差
Figure FDA0002491911260000029
步骤7，根据下式，计算样本集合Sk中每个样本
Figure FDA00024919112600000210
对应的U学习函数值
Figure FDA00024919112600000211
Figure FDA00024919112600000212
步骤8，因此，由于样本集合Sk中共有Nk个样本，共计算得到Nk个U学习函数值，判断每个U学习函数值是否均大于等于2，即：判断下式是否成立：
Figure FDA00024919112600000213
如果成立，则代表当前Kriging模型已收敛，停止更新当前Kriging模型，得到训练好的Kriging模型，执行步骤9；
否则，根据下式从样本集合Sk选取新的样本，也就是选取具有最小U学习函数值的样本，表示为
Figure FDA00024919112600000214
Figure FDA00024919112600000215
计算样本
Figure FDA00024919112600000216
的实际功能函数值
Figure FDA00024919112600000217
将
Figure FDA00024919112600000218
和
Figure FDA00024919112600000219
分别加入训练集合T和实际功能函数值集合gt中，由此得到更新后的训练集合T和更新后的实际功能函数值集合gt，再返回步骤5；
步骤9，利用训练好的Kriging模型，计算样本集合Sk中每个样本
Figure FDA00024919112600000220
对应的代理功能函数均值
Figure FDA0002491911260000031
根据下述指示函数IFk(·)从样本集合Sk中挑选可使第k层失效事件Fk＝{x:g(x)≤bk}发生的样本：
Figure FDA0002491911260000032
若
Figure FDA0002491911260000033
则样本
Figure FDA0002491911260000034
可使第k层失效事件Fk＝{x:g(x)≤bk}发生；
样本集合Sk中，所有使第k层失效事件Fk＝{x:g(x)≤bk}发生的样本组成的集合记为失效样本集合
Figure FDA0002491911260000035
其中，Mk为使第k层失效事件Fk＝{x:g(x)≤bk}发生的失效样本总数；
根据下式计算条件失效概率Pk的估算值
Figure FDA0002491911260000036
Figure FDA0002491911260000037
其中：
Figure FDA0002491911260000038
表示第k层失效事件Fk发生概率的估算值；
Figure FDA0002491911260000039
表示第k-1层失效事件Fk-1发生的条件下，第k层失效事件Fk发生概率的估算值；
P{Fk-1}表示第k-1层失效事件Fk-1发生概率；
Figure FDA00024919112600000310
表示重要抽样密度函数；
根据下式计算条件失效概率Pk的估算值
Figure FDA00024919112600000311
对
Figure FDA00024919112600000312
的偏导数
Figure FDA00024919112600000313
Figure FDA00024919112600000314
其中：
Figure FDA0002491911260000041
代表n维随机向量x＝{x1,x2,…,xn}T中的第r维自变量xr的第s个分布参数；
Pq代表第q层失效事件Fq发生的概率，q＝1,2,...,k-1；
步骤10，判断bk是否等于0，如果bk＝0，则循环过程结束，执行步骤15；否则，执行步骤11；
步骤11，令k＝k+1，在失效样本集合RSk-1中，选出联合概率密度函数值最大的样本，作为重要抽样密度函数
Figure FDA0002491911260000042
的抽样中心；
然后，根据n维随机向量x的重要抽样密度函数
Figure FDA0002491911260000043
以联合概率密度函数值最大的样本为抽样中心，产生Nk个样本构建样本集合
Figure FDA0002491911260000044
步骤12，利用训练好的Kriging模型，计算样本集合Sk中每个样本
Figure FDA0002491911260000045
对应的代理功能函数值
Figure FDA0002491911260000046
和代理功能函数均值
Figure FDA0002491911260000047
根据下述指示函数
Figure FDA0002491911260000048
从样本集合Sk中挑选可使第k层失效事件Fk＝{x:g(x)≤bk}发生的样本：
Figure FDA0002491911260000049
样本集合Sk中，所有使第k层失效事件Fk＝{x:g(x)≤bk}发生的样本组成的集合记为失效样本集合
Figure FDA00024919112600000410
其中，Mk为使第k层失效事件Fk＝{x:g(x)≤bk}发生的失效样本总数；
步骤13，利用训练好的Kriging模型，计算失效样本集合RSk中每个样本对应的代理功能函数均值；然后，将失效样本集合RSk中每个样本对应的代理功能函数均值从大到小排序，
Figure FDA00024919112600000411
取第
Figure FDA00024919112600000412
个代理功能函数均值作为第k层失效事件Fk＝{x:g(x)≤bk}的临界值bk，即
Figure FDA00024919112600000413
步骤14，判断bk是否小于0；如果小于，则令bk＝0，然后返回步骤7；如果不小于，直接返回步骤7；
步骤15，将至此得到的每一层失效事件对应的条件失效概率的估算值做乘法运算，得到最终的涡轮机匣失效概率的估计值
Figure FDA0002491911260000051
具体的，令m＝k，根据下式计算涡轮机匣失效概率的估计值
Figure FDA0002491911260000052
Figure FDA0002491911260000053
通过下式计算得到最终的n维随机向量x＝{x1,x2,…,xn}T中的第r维自变量xr的第s个分布参数的可靠性局部灵敏度指标
Figure FDA0002491911260000054
Figure FDA0002491911260000055
2.根据权利要求1所述的基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法，其特征在于，影响涡轮机匣灵敏度的n维随机向量中，各分量包括高温区域温度、流动冷却区域温度、冲击冷却区域温度、前段机匣长度及机匣内半径。
3.根据权利要求1所述的基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法，其特征在于，还包括：
步骤16，最后一层失效事件即为第m层失效事件Fm，使第m层失效事件Fm发生的样本组成的失效样本集合为RSm，设失效样本集合RSm中失效样本总数为NC；
将失效样本集合RSm中的所有失效样本放入矩阵C，则矩阵C中的样本表示为
Figure FDA0002491911260000056
其中，b＝1,2,...,NC；
步骤17，初始化
Figure FDA0002491911260000057
并令β＝1；
步骤18，构建并计算根据重要抽样密度函数产生的样本
Figure FDA0002491911260000058
向目标分布中第β+1状态zβ+1转化的接受概率
Figure FDA0002491911260000059
为：
Figure FDA0002491911260000061
其中：
Figure FDA0002491911260000062
为样本
Figure FDA0002491911260000063
的联合概率密度函数；
Figure FDA0002491911260000064
为第β状态zβ的重要抽样密度函数；
fX(zβ)为第β状态zβ的联合概率密度函数；
Figure FDA0002491911260000065
为样本
Figure FDA0002491911260000066
的重要抽样密度函数；
步骤19，根据以下判断式，判断zβ+1是否接受
Figure FDA0002491911260000067
Figure FDA0002491911260000068
其中，random[0,1]表示0与1之间的随机数；
即：如果
Figure FDA0002491911260000069
则将
Figure FDA00024919112600000610
赋值给zβ+1；否则，将zβ赋值给zβ+1；
步骤20，令β＝β+1，判断β是否大于NC，如果大于，则执行步骤21；否则，执行步骤18；
步骤21，由此得到NC个zβ，分别为
Figure FDA00024919112600000615
将
Figure FDA00024919112600000616
放入到矩阵B，矩阵B中样本服从输入变量的条件概率密度函数fX(x|F)；
步骤22，根据矩阵B中的样本估计失效域内xr条件概率密度函数
Figure FDA00024919112600000611
再结合步骤15得到的涡轮机匣失效概率的估计值
Figure FDA00024919112600000612
以及第r维自变量xr的联合概率密度函数
Figure FDA00024919112600000613
代入下式得到输入变量xr的可靠性全局灵敏度指标Sr：
Figure FDA00024919112600000614
其中：
E[·]表示期望运算。
CN202010407601.XA 2020-05-14 2020-05-14 基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法 Active CN111783236B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010407601.XA CN111783236B (zh) 2020-05-14 2020-05-14 基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010407601.XA CN111783236B (zh) 2020-05-14 2020-05-14 基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法

Publications (2)

Publication Number Publication Date
CN111783236A true CN111783236A (zh) 2020-10-16
CN111783236B CN111783236B (zh) 2023-03-14

Family

ID=72754091

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010407601.XA Active CN111783236B (zh) 2020-05-14 2020-05-14 基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法

Country Status (1)

Country Link
CN (1) CN111783236B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112528511A (zh) * 2020-12-17 2021-03-19 中国矿业大学 一种刮板输送机链轮磨损可靠性评估方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB0916094D0 (en) * 2007-03-14 2009-10-28 Halliburton Energy Serv Inc Neural-network based surrogate model construction methods and applications thereof
US20140072185A1 (en) * 2012-09-07 2014-03-13 David D. Dunlap Biometric Identification Systems and Methods
CN107563067A (zh) * 2017-09-06 2018-01-09 电子科技大学 基于自适应代理模型的结构可靠性分析方法
CN108287808A (zh) * 2018-01-29 2018-07-17 华中科技大学 一种面向结构可靠性分析的通用动态追踪序列采样方法
CN108304679A (zh) * 2018-03-07 2018-07-20 西北工业大学 一种自适应可靠性分析方法
CN108491284A (zh) * 2018-02-13 2018-09-04 西北工业大学 多失效模式复杂机构可靠性及全局灵敏度分析方法
CN109257749A (zh) * 2018-09-05 2019-01-22 中国人民解放军陆军工程大学 面向动态拓扑的无线传感器网络自适应分层入侵检测方法
CN110532723A (zh) * 2019-09-06 2019-12-03 北京航空航天大学 一种基于egra的涡轮盘多失效模式可靠性优化方法

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB0916094D0 (en) * 2007-03-14 2009-10-28 Halliburton Energy Serv Inc Neural-network based surrogate model construction methods and applications thereof
GB2462380A (en) * 2007-03-14 2010-02-10 Halliburton Energy Serv Inc Neural-network based surrogate model construction methods and applications thereof
US20140072185A1 (en) * 2012-09-07 2014-03-13 David D. Dunlap Biometric Identification Systems and Methods
CN107563067A (zh) * 2017-09-06 2018-01-09 电子科技大学 基于自适应代理模型的结构可靠性分析方法
CN108287808A (zh) * 2018-01-29 2018-07-17 华中科技大学 一种面向结构可靠性分析的通用动态追踪序列采样方法
CN108491284A (zh) * 2018-02-13 2018-09-04 西北工业大学 多失效模式复杂机构可靠性及全局灵敏度分析方法
CN108304679A (zh) * 2018-03-07 2018-07-20 西北工业大学 一种自适应可靠性分析方法
CN109257749A (zh) * 2018-09-05 2019-01-22 中国人民解放军陆军工程大学 面向动态拓扑的无线传感器网络自适应分层入侵检测方法
CN110532723A (zh) * 2019-09-06 2019-12-03 北京航空航天大学 一种基于egra的涡轮盘多失效模式可靠性优化方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
MI XIAO 等: "An efficient Kriging-based subset simulation method for hybrid reliability analysis under random and interval variables with small failure probability", 《STRUCTURAL AND MULTIDISCIPLINARY OPTIMIZATION》 *
周苏婷 等: "可靠性全局灵敏度求解的Meta-IS-AK算法", 《航空学报HTTP://KNS.CNKI.NET/KCMS/DETAIL/11.1929.V.20191010.0931.002.HTML》 *
宋述芳 等: "基于子集模拟和重要抽样的可靠性灵敏度分析方法", 《力学学报》 *
闫阳 等: "基于Kriging模型的子集模拟方法", 《计算机仿真》 *
黄晓旭等: "基于主动学习Kriging模型和子集模拟的可靠度分析", 《应用力学学报》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112528511A (zh) * 2020-12-17 2021-03-19 中国矿业大学 一种刮板输送机链轮磨损可靠性评估方法
CN112528511B (zh) * 2020-12-17 2023-12-12 中国矿业大学 一种刮板输送机链轮磨损可靠性评估方法

Also Published As

Publication number Publication date
CN111783236B (zh) 2023-03-14

Similar Documents

Publication Publication Date Title
Toal et al. Kriging hyperparameter tuning strategies
CN111814956B (zh) 一种基于多维度二次特征提取的多任务学习的空气质量预测方法
CN108304679A (zh) 一种自适应可靠性分析方法
JP2020500420A (ja) マシンラーニング基盤の半導体製造の収率予測システム及び方法
CN112580228B (zh) 一种涡扇发动机混合结构风扇叶片结构设计优化方法
CN113591215B (zh) 基于不确定性的异常卫星组件布局检测方法
CN111797535A (zh) 一种面向多种代理模型的结构可靠性分析自适应加点方法
CN111260138A (zh) 加权和自适应并行关联规则的火电机组能耗动态优化方法
CN111783238A (zh) 涡轮轴结构可靠性分析方法、分析装置及可读存储介质
CN111783236B (zh) 基于自适应模型和子集模拟的涡轮机匣灵敏度分析方法
CN115906399A (zh) 一种改进的小样本数据下产品关键工艺质量预测方法
CN111783237A (zh) 基于Kriging模型的涡轮轴可靠性优化设计方法
J. Toal et al. Performance of an ensemble of ordinary, universal, non-stationary and limit Kriging predictors
CN114996995A (zh) 一种超材料隔振单元性能预报方法和系统
Hart et al. An advanced cost estimation methodology for engineering systems
Brooks MCMC convergence diagnosis via multivariate bounds on log-concave densities
CN113609763A (zh) 基于不确定性的卫星组件布局温度场预测方法
Han et al. A kriging-based active learning algorithm for contour estimation of integrated response with noise factors
Gal et al. Using deep neural networks and derivative free optimization to accelerate coverage closure
Hamidzadeh et al. Clustering data stream with uncertainty using belief function theory and fading function
CN116628854A (zh) 一种翼型气动特性预测方法、系统、电子设备及存储介质
CN114417653B (zh) 基于不确定性的卫星组件电路板热可靠性分析方法
McGurk et al. Data-driven Bayesian inference for stochastic model identification of nonlinear aeroelastic systems
CN114492112A (zh) 基于采样分位数回归的卫星框架结构可靠性分析方法
CN114880792A (zh) 一种基于形变预测的全方位多角度优化方法

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