CN113901665A - 一种基于条件穿越率的时变可靠度准确分析方法 - Google Patents

一种基于条件穿越率的时变可靠度准确分析方法 Download PDF

Info

Publication number
CN113901665A
CN113901665A CN202111221590.7A CN202111221590A CN113901665A CN 113901665 A CN113901665 A CN 113901665A CN 202111221590 A CN202111221590 A CN 202111221590A CN 113901665 A CN113901665 A CN 113901665A
Authority
CN
China
Prior art keywords
time
tau
crossing rate
probability
calculating
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.)
Pending
Application number
CN202111221590.7A
Other languages
English (en)
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.)
Beijing University of Technology
Original Assignee
Beijing University of Technology
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 Beijing University of Technology filed Critical Beijing University of Technology
Priority to CN202111221590.7A priority Critical patent/CN113901665A/zh
Publication of CN113901665A publication Critical patent/CN113901665A/zh
Pending legal-status Critical Current

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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • 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/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Biology (AREA)
  • Computer Hardware Design (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Evolutionary Computation (AREA)
  • Algebra (AREA)
  • Geometry (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Complex Calculations (AREA)

Abstract

本发明公开了一种基于条件穿越率的结构时变可靠度准确分析方法，采用基于条件穿越率v+(t|x)的累积失效概率Pf,c(0,T)进行结构的时变可靠度准确分析；现有的方法仅可以计算累积失效概率的上限值而非准确值，因此现有方法将会低估结构或结构体系的可靠水平，造成不必要的养护维修花费。本发明通过定义一种条件穿越率，提出基于条件穿越率的累积失效概率数值算法，实现对结构服役期内累积失效概率准确值的估算，从而可以避免由于现有基于穿越率的时变可靠度分析方法仅能评估累积失效概率上限值造成的结构运营成本上的浪费。

Description

一种基于条件穿越率的时变可靠度准确分析方法
技术领域
本发明涉及工程结构可靠度分析领域，具体涉及一种基于条件穿越率的时变可靠度准确分析方法。
背景技术
结构可靠性分析旨在计算一个结构或结构体系在其整个服役周期内与其预期功能相关的累积失效概率，即在服役期内无法完成预期功能的概率。实际工程中，由于结构的材料性能(如钢筋屈服强度、混凝土抗压强度、受力钢筋截面尺寸、不同结构层间界面内聚力等)会随时间退化且外部载荷作用(如风荷载、地震荷载、温度作用等)具有明显的时变特性，使得结构的服役可靠性本质上是时变的。显然，与经典的时不变可靠性分析相比，时变可靠性分析由于引入了时间维度而变得更加复杂。
结构时变可靠度分析的首要任务是计算结构在服役期[0，T]内的累积失效概率Pf,c(0，T)。理论上讲，结构在[0，T]时间段内失效定义为其极限状态函数G(X,Y(t),t)在[0,T]内出现至少一次负值，即失效事件Ef＝{t∈[0，T]:G(X,Y(t),t)≤0}。其中，X＝{X1,...,Xi,...,Xn}表示n维随机变量向量，Y(t)＝{Y1(t),...,Yi(t),...,Ym(t)}表示m维存在时间变异性的随机过程向量，t表示时刻。蒙特卡洛方法(Monte Carlo Simulation，MCS)是目前公认最具有普适性的时变可靠度分析方法。该方法的流程包括随机生成输入样本(即X和Y(t)的样本)，计算不同样本下极限状态函数G(X,Y(t),t)的取值，通过统计G(X,Y(t),t)小于零的次数来估计结构在服役期内的累积失效概率。为得到较为准确的估算值，MCS方法需要生成大量随机样本(往往大于105)，因此计算效率较低，主要在科学研究中作为准确性的度量，难以应用于实际工程。
近年来，基于穿越事件的时变可靠度分析方法得到了众多学者的关注，此类方法将一次失效事件表示为等价一次穿越事件，在时刻t发生穿越事件定义为极限状态函数G(X，Y(t)，t)在该时刻由正值转变为负值。相应的，累积失效概率则表示为穿越率v+(t)的函数。
Figure BDA0003312838760000011
其中，Pf,i(0)为初始时刻的时不变失效概率；v+(t)定义为在t时刻发生一次穿越事件的概率。目前，针对基于穿越事件时变可靠度分析方法的研究均集中于穿越率v+(t)的求解。通过将前后时刻的极限状态函数看作并联体系的子系统，经典PHI2时变可靠度分析方法将穿越率v+(t)表示为两个连续时刻时不变可靠指标的二维正态分布函数除以时间间隔Δt的形式。时不变可靠指标通常采用经典一阶时不变可靠度分析方法(First orderreliability method，FORM)计算。针对PHI2方法某些条件下对时间间隔Δt不稳定的问题，PHI2方法的发明者进一步提出了更稳定的PHI2+方法。虽然现有研究对穿越率v+(t)的求解。通过将前后时刻的极限状态函数看作并联体系的子系统，经典PHI2时变可靠的计算方法提出了改进，但现有方法均只能计算[0，T]时间段内累积失效概率Pf,c(0，T)的上限值Pf,c_upper(0，T)，并以此为依据对结构或结构体系的可靠水平进行评估。
由于现有基于穿越率的时变可靠度分析方法实际上仅可以计算累积失效概率的上限值Pf,c_upper(0，T)而非准确值，因此现有方法将会低估结构或结构体系的可靠水平(即高估结构或结构体系的风险)，造成不必要的养护维修花费。
发明内容
本发明针对现有技术存在的问题提供一种基于条件穿越率的时变可靠度准确分析方法。
本发明采用的技术方案是：
一种基于条件穿越率的结构时变可靠度准确分析方法，采用基于条件穿越率v+(t|x)的累积失效概率Pf,c(0,T)进行结构的时变可靠度准确分析；
其中累积失效概率计算方法如下：
Figure BDA0003312838760000021
其中：Ps,i(0)＝1-Pf,i(0)为初始时刻的时不变生存概率，X为随机变量，x为X取的固定值，nx为X的维数，L0为条件累积失效概率
Figure BDA0003312838760000022
在均值点处的估计值，Li为
Figure BDA0003312838760000023
在单个随机变量向量xi0(l)处的估计值，Lij为
Figure BDA0003312838760000024
在双随机变量向量xij0(7l+q-7)处的估计值，f指代失效，c表征累积概率，s指代生存，i表征瞬时概率，Pf,i(0)为瞬时失效概率；其中，
Figure BDA0003312838760000025
Figure BDA0003312838760000026
Figure BDA0003312838760000027
Figure BDA0003312838760000028
式中，μ为X的均值向量，xi0(l)为将μ中第i个元素替换为IN[ur(l)]的向量，xij0(7l+q-7)为分别将μ中的第i和第j个元素替换为IN[ur(l)]和IN[ur(q)]的向量，ωij0(7l+q-7)＝ωr(l)ωr(q)，IN[y]＝F-1[Φ(y)]为逆正态变换，F[·]为随机变量的累积分布函数，Φ[·]为标准正态累积分布函数，y为随机变量，l表征向量中元素的位置，q表征向量中元素的位置；v+(τT(l)|x)为条件穿越概率，ωT(l)为四阶向量，τT(l)为四阶向量。
进一步的，所述初始时刻的时不变生存概率采用FORM时不变可靠度分析方法计算得到，计算过程如下：
S11：极限状态函数为G(X,Y(0),0)，随机变量为{X,Y(0)}；采用YT表示Y(0)对应的随机变量，给定初始验算点{X0,YT0}＝{μx,μY}；
S12：计算正态空间内的验算点u0；
S13：计算初始验算点处的可靠指标
Figure BDA0003312838760000031
S14：计算极限状态函数在验算点u0处的梯度
Figure BDA0003312838760000032
S15：将
Figure BDA0003312838760000033
转换为正态空间内的梯度
Figure BDA0003312838760000034
J为对角矩阵；
S16：生成正态空间内的新验算点u1；
Figure BDA0003312838760000035
S17：计算新验算点处的可靠指标
Figure BDA0003312838760000036
S18：判别β1和β0间相对误差是否满足允许限值要求，如果误差满足要求则结束循环；否则对验算点u1进行逆正态变换得到正常空间内新的验算点，重复步骤S12～S17，直至β1和β0间相对误差满足允许限值要求。
进一步的，采用双重嵌套循环计算Lij，
初始设定l＝0，Lij＝0；
对于给定的l，开始内部循环：
S21：初始设定k＝0，Ev[T|xij(l)]＝0，其中Ev[T|xij(l)]为X＝xij(l)条件下穿越率在时域内的积分；
S22：令X＝xij(l)，t＝τT(k)，建立条件时不变极限状态函数Gc(xij(l)，此时，条件极限状态函数中的随机变量为Y(τT(k))；
S23：计算条件穿越率v+[τT(k)|xij(l)]：
将τT(k)和τT(k)+Δt时刻的随机过程表示为两组相关随机变量Y1和Y2；采用Cholesky分解，将相关的Y1和Y2表示为两组相互独立的随机变量，建立相应的极限状态函数；
计算τT(k)和τT(k)+Δt时刻的瞬时可靠指标β[τT(k)|xij(l)]和β[τT(k)+Δt|xij(l)]，获得相应的验算点u*[τT(k)|xij(l)]和u*[τT(k)+Δt|xij(l)]；计算单位法向向量α[τT(k)|xij(l)]＝u*[τT(k)|xij(l)]/β[τT(k)|xij(l)]和α[τT(k)+Δt|xij(l)]＝u*[τT(k)+Δt|xij(l)]/β[τT(k)+Δt|xij(l)]；
Figure BDA0003312838760000041
式中：a＝u*[τT(k)|xij(l)]·u*[τT(k)+Δt|xij(l)]-β[τT(k)+Δt|xij(l)]β[τT(k)|xij(l)]；ε(·)为单位阶跃函数，θ1和θ2为计算参数，
Figure BDA0003312838760000042
和Λ(θi)为函数；
Figure BDA0003312838760000043
Figure BDA0003312838760000044
式中：β0为计算参数，y为用以描述
Figure BDA0003312838760000045
形式的参数；
Figure BDA0003312838760000046
式中：ρG[τT(k),τT(k)+Δt|xij(l)]＝α[τT(k)|xij(l)]·α[τT(k)+Δt|xij(l)]为τT(k)时刻和τT(k)+Δt时刻条件极限状态函数在正态空间内验算点处线性展开函数间的相关系数。
S24：计算Ev[T|xij(l)]＝Ev[T|xij(l)]+v+[τT(k)|xij(l)]ωT(k)；
S25：判断k是否等于设定值，若等于则结束内循环，否则更新k＝k+1，重复步骤S21～S24；
S26：计算
Figure BDA0003312838760000047
S27：判断l是否等于设定值，若等于则结束外循环，否则更新l＝l+1，重复步骤S22～S23。
进一步的，采用双重嵌套循环计算Li，
Figure BDA0003312838760000048
初始设定l＝0，Li＝0；
对于给定的l，开始内部循环：
S31：初始设定k＝0，Ev[T|xi(l)]＝0，其中Ev[T|xi(l)]为X＝xi(l)条件下穿越率在时域内的积分；
S32：令X＝xi(l)，t＝τT(k)，建立条件时不变极限状态函数Gc[xi(l)]，此时，条件极限状态函数中的随机变量为Y(τT(k))；
S33：计算条件穿越率v+[τT(k)|xi(l)]：
将τT(k)和τT(k)+Δt时刻的随机过程表示为两组相关随机变量Y1和Y2；采用Cholesky分解，将相关的Y1和Y2表示为两组相互独立的随机变量，建立相应的极限状态函数；
计算τT(k)和τT(k)+Δt时刻的瞬时可靠指标β[τT(k)|xi(l)]和β[τT(k)+Δt|xi(l)]，获得相应的验算点u*[τT(k)|xi(l)]和u*[τT(k)+Δt|xi(l)]；计算单位法向向量α[τT(k)|xi(l)]＝u*[τT(k)|xi(l)]/β[τT(k)|xi(l)]和α[τT(k)+Δt|xi(l)]＝u*[τT(k)+Δt|xi(l)]/β[τT(k)+Δt|xi(l)]；
条件穿越率计算方法如下：
Figure BDA0003312838760000051
式中：a＝u*[τT(k)|xi(l)]·u*[τT(k)+Δt|xi(l)]-β[τT(k)+Δt|xi(l)]β[τT(k)|xi(l)]；ε(·)为单位阶跃函数，θ1和θ2为计算参数，
Figure BDA0003312838760000052
和Λ(θi)为函数；
Figure BDA0003312838760000053
Figure BDA0003312838760000054
Figure BDA0003312838760000055
式中：β0为计算参数，y为用以描述
Figure BDA0003312838760000056
形式的参数；
Figure BDA0003312838760000057
式中：ρG[τT(k),τT(k)+Δt|xi(l)]＝α[τT(k)|xi(l)]·α[τT(k)+Δt|xi(l)]为τT(k)时刻和τT(k)+Δt时刻条件极限状态函数在正态空间内验算点处线性展开函数间的相关系数。
S34：计算Ev[T|xi(l)]＝Ev[T|xi(l)]+v+[τT(k)|xi(l)]ωT(k)；
S35：判断k是否等于设定值，若等于则结束内循环，否则更新k＝k+1，重复步骤S21～S24；
S36：计算
Figure BDA0003312838760000058
S37：判断l是否等于设定值，若等于则结束外循环，否则更新l＝l+1，重复步骤S32～S33。
本发明的有益效果是：
(1)本发明提出的基于条件穿越率的累积失效概率数值算法，可以实现对结构服役期内累积失效概率准确值的估算；
(2)本发明数值算法可以避免由于现有基于穿越率的时变可靠度分析方法仅能评估累积失效概率上限值造成的结构运营成本上的浪费；
(3)本发明数值算法计算得到的累积失效概率的准确度更高，可以优化结构的养护维修策略，降低运营成本。
附图说明
图1为本发明数值算法流程示意图。
图2为实施例中锈蚀简支钢梁模型图。
图3为实施例中不同方法得到的累积失效概率对比示意图。
图4为现有的基于穿越率的时变可靠度分析方法流程示意图。
具体实施方式
下面结合附图和具体实施例对本发明做进一步说明。
如图1所示，一种基于条件穿越率的结构时变可靠度准确分析方法，采用基于条件穿越率v+(t|x)的累积失效概率Pf,c(0,T)进行结构的时变可靠度准确分析；
其中累积失效概率计算方法如下：
Figure BDA0003312838760000061
其中：Ps,i(0)＝1-Pf,i(0)为初始时刻的时不变生存概率，X为随机变量，x为X取的固定值，nx为X的维数，L0为条件累积失效概率
Figure BDA0003312838760000062
在均值点处的估计值，Li为
Figure BDA0003312838760000063
在单个随机变量向量xi0(l)处的估计值，Lij为
Figure BDA0003312838760000064
在双随机变量向量xij0(7l+q-7)处的估计值，f指代失效，c表征累积概率，s指代生存(不失效)，i表征瞬时(时不变)概率，Pf,i(0)为瞬时失效概率；
其中，
Figure BDA0003312838760000065
Figure BDA0003312838760000066
Figure BDA0003312838760000067
Figure BDA0003312838760000071
式中，μ为X的均值向量，xi0(l)为将μ中第i个元素替换为IN[ur(l)]的向量，xij0(7l+q-7)为分别将μ中的第i和第j个元素替换为IN[ur(l)]和IN[ur(q)]的向量，ωij0(7l+q-7)＝ωr(l)ωr(q)，IN[y]＝F-1[Φ(y)]为逆正态变换，F[·]为随机变量的累积分布函数，Φ[·]为标准正态累积分布函数，y为随机变量，l表征向量中元素的位置，q表征向量中元素的位置；v+(τT(l)|x)为条件穿越概率，ωT(l)为四阶向量，τT(l)为四阶向量。IN[y]采用四阶矩逆正态变换模型定义，ur＝{±3.750439,±2.366759,±1.154405,0}；ωi0＝ωr＝{5.4826×10-4,5.4826×10-4,3.0757×10-2,3.0757×10-2,0.240123,0.240123,0.457142}；τT＝T/2τ0+T/2,τ0＝{±0.3400，±0.8611},ωT＝{0.6521,0.6521,0.3479,0.3479}。
为了反映随机变量X在整个服役期内不随时间变化的特征，定义一种在以X取固定值x为前提的条件穿越事件。在t时刻发生条件穿越事件意味着条件极限状态函数Gc(x,Y(t),t)在该时刻由正值变为负值。相应的，在X＝x前提下，t时刻发生一次条件穿越事件的概率定义为条件穿越率v+(t|x)。基于条件穿越概率v+(t|x)的累积失效概率Pf,c(0,T)定义为：
Figure BDA0003312838760000072
Figure BDA0003312838760000073
其中，Ps,i(0)＝1-Pf,i(0)为初始时刻的时不变生存概率；Ωx为X的定义域；fx(x)为X的联合概率密度函数；Pc(T|x)定义为(0,T]时间段内的条件累积失效概率。条件穿越概率v+(t|x)的表达式为：
Figure BDA0003312838760000074
Figure BDA0003312838760000075
Figure BDA0003312838760000076
其中，a＝u*(t|x)·u*(t+Δt|x)-β(t+Δt|x)β(t|x)；u*(t|x)和u*(t+Δt|x)分别为t时刻和t+Δt时刻条件极限状态函数在正态空间内的验算点；β(t|x)和β(t+Δt|x)分别为t时刻和t+Δt时刻条件极限状态函数的时不变可靠指标，采用FORM计算；ε(·)为单位阶跃函数，即当括号内数值为正数时ε(·)＝1，反之ε(·)＝0；ρG(t,t+Δt|x)＝u*(t|x)·u*(t+Δt|x)/β(t+Δt|x)/β(t|x)为t时刻和t+Δt时刻条件极限状态函数在正态空间内验算点处线性展开函数间的相关系数。
本发明计算方法的流程如下：
1、采用FORM时不变可靠度分析方法计算初始时刻的时不变生存概率Ps,i(0)，此时考虑的极限状态函数为G(X,Y(0),0)，随机变量为{X,Y(0)}。
FORM时不变可靠度分析方法计算过程如下：
S11：极限状态函数为G(X,Y(0),0)，随机变量为{X,Y(0)}；采用YT表示Y(0)对应的随机变量，给定初始验算点{X0,YT0}＝{μx,μY}；
S12：计算正态空间内的验算点u0；对于分布已知的变量，采用Rosenblatt变换计算验算点在正态空间内的映射u0；对于分布未知的变量，采用四阶矩正态变换计算验算点在正态空间内的映射u0。
S13：计算初始验算点处的可靠指标
Figure BDA0003312838760000081
S14：计算极限状态函数在验算点u0处的梯度
Figure BDA0003312838760000082
S15：将
Figure BDA0003312838760000083
转换为正态空间内的梯度
Figure BDA0003312838760000084
J为对角矩阵；其对角线上的第i个元素定义为φ(ui)/f(xi)；xi为{X0,YT0}的第i个元素，f(xi)表示第xi的概率密度函数；ui为u0的第i个元素。
S16：生成正态空间内的新验算点u1；
Figure BDA0003312838760000085
S17：计算新验算点处的可靠指标
Figure BDA0003312838760000086
S18：判别β1和β0间相对误差是否满足允许限值要求，如果误差满足要求则结束循环，FORM计算得到的可靠指标为β1；否则对验算点u1进行逆正态变换得到正常空间内新的验算点{X0,YT0}＝IN(u1)，重复步骤S12～S17，直至β1和β0间相对误差满足允许限值要求。其中IN(·)为逆正态运算，对于分布已知的变量，采用Rosenblatt变换定义；对于分布未知的变量，采用四阶矩逆正态变换定义。
2、生成用于数值积分的各类矩阵，49(nx-1)nx/2×nx维矩阵xij是所有xij0的集合(i＝1,...,nx；j＝i+1,...,nx)；49(nx-1)nx/2维向量ωij为所有ωij0(i＝1,...,nx；j＝i+1,...,nx)的集合；7nx×nx维矩阵xi为所有xi0(i＝1,...,nx)的集合；7nx维向量ωi所有ωi0(i＝1,...,nx)的集合；τT＝T·τ0/2+T/2。
3、采用双重嵌套循环计算Lij，
初始设定l＝0，Lij＝0；
对于给定的l，开始内部循环：
S21：初始设定k＝0，Ev[T|xij(l)]＝0，其中Ev[T|xij(l)]为X＝xij(l)条件下穿越率在时域内的积分；
S22：令X＝xij(l)，t＝τT(k)，建立条件时不变极限状态函数Gc(xij(l)，此时，条件极限状态函数中的随机变量为Y(τT(k))；
S23：计算条件穿越率v+[τT(k)|xij(l)]：
将τT(k)和τT(k)+Δt时刻的随机过程表示为两组相关随机变量Y1和Y2；采用Cholesky分解，将相关的Y1和Y2表示为两组相互独立的随机变量，建立相应的极限状态函数；
采用FORM方法计算τT(k)和τT(k)+Δt时刻的瞬时可靠指标β[τT(k)]和β[τT(k)+Δt]，获得相应的验算点u*[τT(k)]和u*[τT(k)+Δt]；计算单位法向向量α[τT(k)]＝u*[τT(k)]/β[τT(k)]和α[τT(k)+Δt]＝u*[τT(k)+Δt]/β[τT(k)+Δt]；
Figure BDA0003312838760000091
式中：a＝u*(t|x)·u*(t+Δt|x)-β(t+Δt|x)β(t|x)；u*(t|x)和u*(t+Δt|x)分别为t时刻和t+Δt时刻条件极限状态函数在正态空间内的验算点；β(t|x)和β(t+Δt|x)分别为t时刻和t+Δt时刻条件极限状态函数的时不变可靠指标；ε(·)为单位阶跃函数，θ1和θ2为计算参数，
Figure BDA0003312838760000092
和Λ(θi)为函数；
Figure BDA0003312838760000093
Figure BDA0003312838760000094
式中：β0为计算参数，y为用以描述
Figure BDA0003312838760000095
形式的参数；
Figure BDA0003312838760000096
式中：ρG(t,t+Δt|x)＝α(t|x)·α(t+Δt|x)为t时刻和t+Δt时刻条件极限状态函数在正态空间内验算点处线性展开函数间的相关系数。
S24：计算Ev[T|xij(l)]＝Ev[T|xij(l)]+v+[τT(k)|xij(l)]ωT(k)；
S25：判断k是否等于设定值此处设定值为4，若等于则结束内循环，否则更新k＝k+1，重复步骤S21～S24；
S26：计算
Figure BDA0003312838760000101
S27：判断l是否等于设定值49nx(nx-1)/2，若等于则结束外循环，否则更新l＝l+1，重复步骤S22～S23。
4、采用双重嵌套循环计算Li，
Figure BDA0003312838760000102
初始设定l＝0，Li＝0；
对于给定的l，开始内部循环：
S31：初始设定k＝0，Ev[T|xi(l)]＝0，其中Ev[T|xi(l)]为X＝xi(l)条件下穿越率在时域内的积分；
S32：令X＝xi(l)，t＝τT(k)，建立条件时不变极限状态函数Gc[xi(l)]，此时，条件极限状态函数中的随机变量为Y[τT(k)]；
S33：计算条件穿越率v+[τT(k)|xi(l)]：
将τT(k)和τT(k)+Δt时刻的随机过程表示为两组相关随机变量Y1和Y2；采用Cholesky分解，将相关的Y1和Y2表示为两组相互独立的随机变量，建立相应的极限状态函数；
采用FORM方法计算τT(k)和τT(k)+Δt时刻的瞬时可靠指标β[τT(k)]和β[τT(k)+Δt]，获得相应的验算点u*[τT(k)]和u*[τT(k)+Δt]；计算单位法向向量α[τT(k)]＝u*[τT(k)]/β[τT(k)]和α[τT(k)+Δt]＝u*[τT(k)+Δt]/β[τT(k)+Δt]；
Figure BDA0003312838760000103
式中：a＝u*(t|x)·u*(t+Δt|x)-β(t+Δt|x)β(t|x)；u*(t|x)和u*(t+Δt|x)分别为t时刻和t+Δt时刻条件极限状态函数在正态空间内的验算点；β(t|x)和β(t+Δt|x)分别为t时刻和t+Δt时刻条件极限状态函数的时不变可靠指标；ε(·)为单位阶跃函数，θ1和θ2为计算参数，
Figure BDA0003312838760000104
和Λ(θi)为函数；
Figure BDA0003312838760000105
Figure BDA0003312838760000106
式中：β0为计算参数，y为用以描述
Figure BDA0003312838760000107
形式的参数；
Figure BDA0003312838760000108
式中：ρG(t,t+Δt|x)＝α(t|x)·α(t+Δt|x)为τT(k)时刻和τT(k)+Δt时刻条件极限状态函数在正态空间内验算点处线性展开函数间的相关系数。
S34：计算Ev[T|xi(l)]＝Ev[T|xi(l)]+v+[τT(k)|xi(l)]ωT(k)；
S35：判断k是否等于设定值，若等于则结束内循环，否则更新k＝k+1，重复步骤S21～S24；
S36：计算
Figure BDA0003312838760000111
S37：判断l是否等于设定值，若等于则结束外循环，否则更新l＝l+1，重复步骤S32～S33。
5、计算
Figure BDA0003312838760000112
6、计算
Figure BDA0003312838760000113
目前基于穿越率的时变可靠度分析方法主要包括PHI2及PHI2+时变可靠度分析方法。两者的计算流程完全一致，仅在计算穿越率v+(t)时采用不同的公式。在PHI2方法中，穿越率表示为：
Figure BDA0003312838760000114
式中：Φ2(·,·；·)为二维正态累积分布函数；β(t)和β(t+Δt)表示极限状态函数G(X,Y(t),t)在t和t+Δt时刻对应的瞬时可靠指标，采用FORM计算；ρG(t,t+Δt)＝-α(t)·α(t+Δt)为G(X,Y(t),t)和-G(X,Y(t+Δt),t+Δt)在正态空间内验算点处线性展开函数的相关系数；α(t)＝u*(t)/β(t)和α(t+Δt)＝u*(t+Δt)/β(t+Δt)分别为t和t+Δt时刻标准正态空间内极限状态面在验算点处的单位法向向量；u*(t)和u*(t+Δt)分别为t时刻和t+Δt时刻的验算点。为满足精度要求，PHI2方法要求Δt的选取应能够使得Y(t)和Y(t+Δt)的自相关系数在[0.99,0.995]之间。
PHI2+方法中，穿越率表示为：
Figure BDA0003312838760000115
式中，||·||表示矩阵的秩；
Figure BDA0003312838760000116
表示标准正态分布的概率密度函数；
Figure BDA0003312838760000117
表示一种函数运算，其中Φ[y]表示标准正态分布累积分布函数在y处的取值。在PHI2+方法中，Δt一般取为Y(t)自相关长度的1％左右。
v+(t)的确定依赖于采用FORM分析的时不变可靠指标。由于FORM方法为基于求导迭代的优化运算，基于FORM的v+(t)没有解析解。因此，v+(t)在时域内的积分通常采用数值方法近似计算。通过将时间段[0,T]均匀离散为Nt个时刻，PHI2/PHI2+方法将累积失效概率的上限值近似表示为：
PHI2方法：
Figure BDA0003312838760000121
PHI2+方法：
Figure BDA0003312838760000122
其中，tj表示第i个离散时刻，dt为离散时间段采用的时间间隔。
基于FORM的时变可靠度方法分析流程如图3所示，具体过程如下：
1、根据结构的极限状态建立相应的极限状态函数G(X,Y(t),t)。
2、基于FORM方法计算初始时刻的瞬时失效概率Pf,i(0)。FORM时不变可靠度分析方法的流程为：
(a)给定初始验算点{X0,Y0}＝{X,Y}。
(b)计算正态空间内的验算点。对于分布已知的变量，采用Rosenblatt变换计算验算点在正态空间内的映射u0；对于分布未知的变量，采用四阶矩正态变换(Fourth momentnormal transformation,FMNT)计算验算点在正态空间内的映射u0。
(c)计算初始验算点处的可靠指标
Figure BDA0003312838760000123
(d)计算极限状态函数在验算点{X0,Y0}处的梯度
Figure BDA0003312838760000124
(e)基于Jacobian矩阵J将
Figure BDA0003312838760000125
转换为正态空间内的梯度
Figure BDA0003312838760000126
其中，J为对角矩阵，其对角线上的第i个元素定义为φ(ui)/f(xi)；xi为{X0，Y0}的第i个元素，f(xi)表示第xi的概率密度函数；ui为u0的第i个元素。
(f)生成正态空间内的新验算点u1。
Figure BDA0003312838760000127
(g)计算新验算点处的可靠指标
Figure BDA0003312838760000128
(h)判别β1和β0间相对误差是否满足允许限值要求，如果误差满足要求则结束循环，FORM计算得到的可靠指标为β1；否则，对u1进行逆正态变换得到正常空间内新的验算点{X0,Y0}＝IN(u1)，重复步骤(b)-(g)直至β1和β0间相对误差满足允许限值要求。其中IN(·)为逆正态运算，对于分布已知的变量，采用Rosenblatt变换定义；对于分布未知的变量，采用四阶矩逆正态变换定义。
3、设置用于离散[0，T]的时间间隔dt以及用于分解Y(t)的时间增量Δt。一般情况下，dt越小越准确(以使Nt≥100为宜)。当采用PHI2方法时，Δt取为使Y(t)及Y(t+Δt)的相关系数在[0.99,0.995]范围内；采用PHI2+方法时，Δt取为Y(t)自相关长度的1％。
4、设置时域内数值积分的初始值t＝dt及Pf,c(0,T)＝Pf,i(0)。
5、将随机过程中的时间变量固定，将随机过程表示为两组相关随机变量Y(t)和Y(t+Δt)。随后，采用Cholesky分解，将相关的Y(t)和Y(t+Δt)表示为两组相互独立的随机变量，建立相应的极限状态函数。
6、采用FORM方法计算t和t+Δt时刻的瞬时可靠指标β(t)和β(t+Δt)，获得相应的验算点u*(t)和u*(t+Δt)。
7、计算单位法向向量α(t)＝u*(t)/β(t)和α(t+Δt)＝u*(t+Δt)/β(t+Δt)。
8、采用PHI2方法或PHI2+方法计算t时刻的穿越率v+(t)，更新Pf,c_upper(0,T)＝Pf,c_upper(0,T)+v+(t)dt。
9、更新t＝t+dt，重复步骤3-8直至t＝T。
现有上述方法仅能估计累积失效概率的上限值，无法得到累积失效概率的准确值。因此，这种方法将会过高的估计结构的风险，低估结构的服役可靠性，进而增加养护维修次数，提高结构的运营成本。
下面以如图2所示的锈蚀简支钢梁模型分析时变可靠度，考虑的极限状态由跨中截面的极限抗弯能力决定，相应的极限状态函数由下式给出：
Figure BDA0003312838760000131
式中，fy为钢材屈服强度；b0为梁截面长度；h0为梁截面宽度；L＝5m为钢梁跨度；ρst＝78.5kN/m3为钢材密度；Fi(t)(i＝1,...,10)为集中荷载；k＝0.03mm/年为钢材锈蚀速率。集中荷载F(t)为高斯过程，其余变量均为相互独立的随机变量。随机过程及随机变量统计信息列于表1。
表1.示例中的随机变量和随机过程统计特征
Figure BDA0003312838760000141
表2.示例中的随机变量和随机过程统计特征
Figure BDA0003312838760000142
本实施例中分别采用PHI2方法、PHI2+方法、MCS方法以及本发明方法分析示例的时变可靠度，以MCS方法得到的结果作为验证不同方法准确性的标准。其中PHI2和PHI2+以0.1年为时间间隔均匀划分时间段，而MCS以0.05年为时间间隔均匀划分时间段且每个时刻采用107个样本进行分析。预测时间段[0,T]的长度变化范围为1年到30年，即T∈[1年，30年]。不同方法得到的累积失效概率如图4所示，其代表值如表2所示。
可以看出现有的PHI2方法和PHI2+方法得到的累积失效概率均明显大于MCS得到的准确值，而本发明方法的计算结果与MCS得到的准确值完全一致。这说明采用现有方法分析累积失效概率将会高估结构面临的风险，可能导致结构运营过程中开展不必要的养护维修，提高结构运营成本。而本发明方法可以准确估计结构的累积失效概率，因此可以更好的优化结构的养护维修策略，进而降低运营成本。

Claims (4)

1.一种基于条件穿越率的结构时变可靠度准确分析方法，其特征在于，采用基于条件穿越率v+(t|x)的累积失效概率Pf，c(0，T)进行结构的时变可靠度准确分析；
其中累积失效概率计算方法如下：
Figure FDA0003312838750000011
其中：Ps，i(0)＝1-Pf，i(0)为初始时刻的时不变生存概率，X为随机变量，x为X取的固定值，nx为X的维数，L0为条件累积失效概率
Figure FDA0003312838750000012
在均值点处的估计值，Li为
Figure FDA0003312838750000013
在单个随机变量向量xi0(l)处的估计值，Lij为
Figure FDA0003312838750000014
在双随机变量向量xij0(7l+q-7)处的估计值，f指代失效，c表征累积概率，s指代生存，i表征瞬时概率，Pf，i(0)为瞬时失效概率；其中，
Figure FDA0003312838750000015
Figure FDA0003312838750000016
Figure FDA0003312838750000017
Figure FDA0003312838750000018
式中，μ为X的均值向量，xi0(l)为将μ中第i个元素替换为IN[ur(l)]的向量，xij0(7l+q-7)为分别将μ中的第i和第j个元素替换为IN[ur(l)]和IN[ur(q)]的向量，ωij0(7l+q-7)＝ωr(l)ωr(q)，IN[y]＝F-1[Φ(y)]为逆正态变换，F[·]为随机变量的累积分布函数，Φ[·]为标准正态累积分布函数，y为随机变量，l表征向量中元素的位置，q表征向量中元素的位置；v+(τT(l)|x)为条件穿越概率，ωT(l)为四阶向量，τT(l)为四阶向量。
2.根据权利要求1所述的一种基于条件穿越率的结构时变可靠度准确分析方法，其特征在于，所述初始时刻的时不变生存概率采用FORM时不变可靠度分析方法计算得到，计算过程如下：
S11：极限状态函数为G(X，Y(0)，0)，随机变量为{X，Y(0)}；采用YT表示Y(0)对应的随机变量，给定初始验算点{X0，YT0}＝{μX，μY}；
S12：计算正态空间内的验算点u0；
S13：计算初始验算点处的可靠指标
Figure FDA0003312838750000019
S14：计算极限状态函数在验算点u0处的梯度
Figure FDA00033128387500000110
S15：将
Figure FDA0003312838750000021
转换为正态空间内的梯度
Figure FDA0003312838750000022
J为对角矩阵；
S16：生成正态空间内的新验算点u1；
Figure FDA0003312838750000023
S17：计算新验算点处的可靠指标
Figure FDA0003312838750000024
S18：判别β1和β0间相对误差是否满足允许限值要求，如果误差满足要求则结束循环；否则对验算点u1进行逆正态变换得到正常空间内新的验算点，重复步骤S12～S17，直至β1和β0间相对误差满足允许限值要求。
3.根据权利要求1所述的一种基于条件穿越率的结构时变可靠度准确分析方法，其特征在于，采用双重嵌套循环计算Lij，
初始设定l＝0，Lij＝0；
对于给定的l，开始内部循环：
S21：初始设定k＝0，Ev[T|xij(l)]＝0，其中Ev[T|xij(l)]为X＝xij(l)条件下穿越率在时域内的积分；
S22：令X＝xij(l)，t＝τT(k)，建立条件时不变极限状态函数Gc(xij(l)，此时，条件极限状态函数中的随机变量为Y(τT(k))；
S23：计算条件穿越率v+[τT(k)|xij(l)]：
将τT(k)和τT(k)+Δt时刻的随机过程表示为两组相关随机变量Y1和Y2；采用Cholesky分解，将相关的Y1和Y2表示为两组相互独立的随机变量，建立相应的极限状态函数；
计算τT(k)和τT(k)+Δt时刻的瞬时可靠指标β[τT(k)|xij(l)]和β[τT(k)+Δt|xij(l)]，获得相应的验算点u＊[τT(k)|xij(l)]和u＊[τT(k)+Δt|xij(l)]；计算单位法向向量α[τT(k)|xij(l)]＝u＊[τT(k)|xij(l)]/β[τT(k)|xij(l)]和α[τT(k)+Δt|xij(l)]＝u＊[τT(k)+Δt|xij(l)]/β[τT(k)+Δt|xij(l)]；
Figure FDA0003312838750000025
式中：a＝u＊[τT(k)|xij(l)]·u＊[τT(k)+Δt|xij(l)]-β[τT(k)+Δt|xij(l)]β[τT(k)|xij(l)]；ε(·)为单位阶跃函数，θ1和θ2为计算参数，
Figure FDA0003312838750000026
和Λ(θi)为函数；
Figure FDA0003312838750000027
Figure FDA0003312838750000028
式中：β0为计算参数，y为用以描述
Figure FDA0003312838750000031
形式的参数；
Figure FDA0003312838750000032
式中：ρG[τT(k)，τT(k)+Δt|xij(l)]＝α[τT(k)|xij(l)]·α[τT(k)+Δt|xij(l)]为τT(k)时刻和τT(k)+Δt时刻条件极限状态函数在正态空间内验算点处线性展开函数间的相关系数。
S24：计算Ev[T|xij(l)]＝Ev[T|xij(l)]+v+[τT(k)|xij(l)]ωT(k)；
S25：判断k是否等于设定值，若等于则结束内循环，否则更新k＝k+1，重复步骤S21～S24；
S26：计算
Figure FDA0003312838750000033
S27：判断l是否等于设定值，若等于则结束外循环，否则更新l＝l+1，重复步骤S22～S23。
4.根据权利要求1所述的一种基于条件穿越率的结构时变可靠度准确分析方法，其特征在于，采用双重嵌套循环计算Li，
Figure FDA0003312838750000034
初始设定l＝0，Li＝0；
对于给定的l，开始内部循环：
S31：初始设定k＝0，Ev[T|xi(l)]＝0，其中Ev[T|xi(l)]为X＝xi(l)条件下穿越率在时域内的积分；
S32：令X＝xi(l)，t＝τT(k)，建立条件时不变极限状态函数Gc[xi(l)]，此时，条件极限状态函数中的随机变量为Y(τT(k))；
S33：计算条件穿越率v+[τT(k)|xi(l)]：
将τT(k)和τT(k)+Δt时刻的随机过程表示为两组相关随机变量Y1和Y2；采用Cholesky分解，将相关的Y1和Y2表示为两组相互独立的随机变量，建立相应的极限状态函数；
计算τT(k)和τT(k)+Δt时刻的瞬时可靠指标β[τT(k)|xi(l)]和β[τT(k)+Δt|xi(l)]，获得相应的验算点u＊[τT(k)|xi(l)]和u＊[τT(k)+Δt|xi(l)]；计算单位法向向量α[τT(k)|xi(l)]＝u＊[τT(k)|xi(l)]/β[τT(k)|xi(l)]和α[τT(k)+Δt|xi(l)]＝u＊[τT(k)+Δt|xi(l)]/β[τT(k)+Δt|xi(l)]；
条件穿越率计算方法如下：
Figure FDA0003312838750000041
式中：a＝u*[τT(k)|xi(l)]·u＊[τT(k)+Δt|xi(l)]-β[τT(k)+Δt|xi(l)]β[τT(k)|xi(l)]；ε(·)为单位阶跃函数，θ1和θ2为计算参数，
Figure FDA0003312838750000042
和Λ(θi)为函数；
Figure FDA0003312838750000043
Figure FDA0003312838750000044
Figure FDA0003312838750000045
式中：β0为计算参数，y为用以描述
Figure FDA0003312838750000046
形式的参数；
Figure FDA0003312838750000047
式中：ρG[τT(k)，τT(k)+Δt|xi(l)]＝α[τT(k)|xi(l)]·α[τT(k)+Δt|xi(l)]为τT(k)时刻和τT(k)+Δt时刻条件极限状态函数在正态空间内验算点处线性展开函数间的相关系数。
S34：计算Ev[T|xi(l)]＝Ev[T|xi(l)]+v+[τT(k)|xi(l)]ωT(k)；
S35：判断k是否等于设定值，若等于则结束内循环，否则更新k＝k+1，重复步骤S21～S24；
S36：计算
Figure FDA0003312838750000048
S37：判断l是否等于设定值，若等于则结束外循环，否则更新l＝l+1，重复步骤S32～S33。
CN202111221590.7A 2021-10-20 2021-10-20 一种基于条件穿越率的时变可靠度准确分析方法 Pending CN113901665A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111221590.7A CN113901665A (zh) 2021-10-20 2021-10-20 一种基于条件穿越率的时变可靠度准确分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111221590.7A CN113901665A (zh) 2021-10-20 2021-10-20 一种基于条件穿越率的时变可靠度准确分析方法

Publications (1)

Publication Number Publication Date
CN113901665A true CN113901665A (zh) 2022-01-07

Family

ID=79192930

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111221590.7A Pending CN113901665A (zh) 2021-10-20 2021-10-20 一种基于条件穿越率的时变可靠度准确分析方法

Country Status (1)

Country Link
CN (1) CN113901665A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114117873A (zh) * 2022-01-25 2022-03-01 浙江大学 基于重要性采样代理模型的复杂装备时变可靠性分析方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114117873A (zh) * 2022-01-25 2022-03-01 浙江大学 基于重要性采样代理模型的复杂装备时变可靠性分析方法
CN114117873B (zh) * 2022-01-25 2022-05-03 浙江大学 基于重要性采样代理模型的复杂装备时变可靠性分析方法

Similar Documents

Publication Publication Date Title
Karve et al. Digital twin approach for damage-tolerant mission planning under uncertainty
CN108052770B (zh) 一种考虑时变效应的大跨桥梁主梁性能预警方法
Ding et al. Structural damage detection using artificial bee colony algorithm with hybrid search strategy
Hombal et al. Surrogate modeling of 3D crack growth
CN106372278A (zh) 一种联合考虑输入参数不确定性和代理模型不确定性的灵敏度分析方法
Chowdhury et al. Probabilistic fracture mechanics by using Monte Carlo simulation and the scaled boundary finite element method
CN114547974B (zh) 基于输入变量选择与lstm神经网络的动态软测量建模方法
CN103324798B (zh) 基于区间响应面模型的随机模型修正方法
CN113065702B (zh) 基于st-seep分段法和时空arma模型的滑坡位移多线性预测方法
CN110096805B (zh) 基于改进自助法的桥梁结构参数不确定性量化及传递方法
JP5707230B2 (ja) プロセスの状態予測方法
He et al. Frequency modification of continuous beam bridge based on co-integration analysis considering the effect of temperature and humidity
Song et al. BUAK-AIS: Efficient Bayesian updating with active learning Kriging-based adaptive importance sampling
CA2689246C (en) Monitoring methods and apparatus
CN112435054B (zh) 基于广义最大相关熵准则的核极限学习机售电量预测方法
CN113010954A (zh) 桥梁结构损伤识别方法、装置及终端设备
CN113901665A (zh) 一种基于条件穿越率的时变可靠度准确分析方法
Ouyang et al. Non-probabilistic uncertain inverse problem method considering correlations for structural parameter identification
Balu et al. Inverse structural reliability analysis under mixed uncertainties using high dimensional model representation and fast Fourier transform
Tang et al. Novel solution framework for inverse problem considering interval uncertainty
CN112241832B (zh) 一种产品质量分级评价标准设计方法及系统
CN116187153B (zh) 基于层次贝叶斯的水工结构数字孪生模型更新方法
Yang et al. A method for time‐varying reliability calculation of RC bridges considering random deterioration process
Wang et al. A combined method of autoregressive model and matrix factorization for recovery and forecasting of cyclic structural health monitoring data
Glashier et al. An iterative regression-based thermal response prediction methodology for instrumented civil infrastructure

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