可靠性和可靠性灵敏度分析的Monte Carlo数值模拟法
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
机械系统可靠性和可靠性灵敏度分析 的Monte Carlo 数值模拟法
授课教师：黄贤振
东北大学机械设计及理论研究所
1 机械系统
机械系统是指由多个机械基本要素（如机械零、部件， 机构，机器等）组成的，可以完成所需动作 ( 或动作过 程)，实现动作的传递和力的变换以及机械能的转化和利 用的装置。
机械系统
E I F ( x ) I F
N 1 ˆ I I (x ) P ˆ E P f F F j f N j 1
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
失效概率估计值的方差可以通过对式(1)两边求方差得如下：
x j 独立 N N 1 1 ˆ Var Var P Var I F ( x j ) I F ( x j ) 2 f N N j 1 j 1
Z g ( x ) g ( x1 , x2 , ..., xn )
则极限状态方程g(x1, x2, …, xn)=0将结构的基本随机变量空 间分为失效域和可靠区域两部分。
东北大学机械设计及理论研究所
2 机械可靠性
失效概率Pf可表示为
Pf
g ( x )0
f X ( x1 , x2 ,
3 Monte Carlo 可靠性分析
指示函数IF(x)方差的无偏估计可以进一步表达为
2 N N 1 N 2 N 1 1 2 2 Var I F ( x ) I F ( x j ) I F ( xk ) I F ( x j ) NI F N 1 j 1 N k 1 N 1 N j 1 ˆ P ˆ2) N ( P N 1 N 2 f f ˆ I F ( x j ) Pf N 1 N j 1 N 1 ？
xj与母体x独立 同分布
ˆ 1 Var I ( x ) Var P F j f N
ˆ Var P f
1 ˆ P ˆ2) (P f f N 1
ˆ ] ˆ Var[ P 1 P f f ˆ Cov( Pf ) ˆ ] ˆ E[ P ( N 1) P f f
E[.]为数学期望算子。
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
失效概率为失效域指示函数的数学期望，依据大数定律，失 效域指示函数的数学期望可以有失效域指示函数的样本均值 来近似。
N Nf 1 ˆ Pf I F ( x j ) N j 1 N
(1)
即 随 机 变 量 的 联 合 概率密度 函数 fX(x) 抽 取 N 个 样本 xj(j=1, 2, …, N), 落入失效域F内样本点的个数Nf与总样本点的个数 N之比即为失效概率的估计值。
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
ˆ P ，即为Pf的无偏估计 由上式可知， E P f f
在数值模拟的过程中，以指示函数IF(x)的样本均值 I F 近似代替 E[IF(x)], 则失效概率估计值的期望可以近似表达为
N 1 ˆ E I F ( x j ) E I F ( x j ) Pf P f N E j 1
将上式代人式(2)，可得失效概率估计值的方差估计为 1 ˆ ˆ P ˆ2) Var Pf (P f f N 1 进而得到估计值的变异系数为
ˆ ] ˆ Var[ P 1 P f f ˆ Cov( Pf ) ˆ ] ˆ E[ P ( N 1) P f f
东北大学机械设计及理论研究所
xS
f R ( xR )dxR
应力xS落在小区间dxS内，同时强度xR小于应力xS的概率为
P xR xS , xS dxS f S ( xS )dxS
xS
f R ( xR )dxR
根据可靠度的定义，对于应力xS所有的可能值强度 xR均小于应力xS的概率， 及事件(xR<xS)的边缘概率就是机械系统的可靠度(对上式的求和)
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
Monte Carlo 可靠性分析方法又称随机抽样法、概率模拟法 或统计试验法。该方法是通过随机模拟或者说统计试验来 进行结构可靠性分析的。由于它是以概率和数理统计理论 为基础的，故被无理学家以赌城Monte Carlo来命名。
东北大学机械设计及理论研究所
1 n 上式表明样本均值 n xi 是依概率收敛于母体的均值 i 1
μ 的。
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
另外，设随机事件 A 发生的概率 P(A), 在n次独立试验中， 事件A发生的频率为m，则随机事件A发生的频率W(A)=m/n, 对于任意ε>0，有
m lim P P( A) 1 n n
, xn )dx1dx2
dxn
式中，f(x1, x2, …, xn) 是基本随机变量x=[x1, x2, …, xn]T的联合
概率密度函数。若基本随机变量是相互独立的，则有
Pf
g ( x )0
f X1 ( x1 ) f X 2 ( x2 )
f X n ( xn )dx1dx2
dxn
式中, f(xi) (i=1, 2, …, n)为随机变量xi的概率密度函数。
多体系统是指由多个 刚性或柔性体通过一 定的联接方式互相关 联起来的完成预期动 作的一个整体
若干元件组成的承受 外部作用力并有特定 功能的整体
多
体 机
械
结 构
系
统
统 系
由构件和运动副组成 的机械系统
机
构
系
统
图1 机械系统的分类
东北大学机械设计及理论研究所
2 机械可靠性
机械可靠性模型 载荷 分布参数
设置抽样样本数N， 初始值m=0, j=0 j=j+1
由概率密度函数fX(x)参数随机样本点xj
将样本xj 代人功能函数，计算功能函数值g(xj)
IF(xj)=0
否 g(xj)≤0 ?
是 IF(xj)=1
N 1 ˆ Var Var P I ( x ) F j f N j 1
xj相互独立
N 1 ˆ Var P Var I F ( x j ) f N2 j 1
ˆ P ˆ2) N 1 N N (P 2 f f ˆ I F ( x j ) Pf N 1 N j 1 N 1
东北大学机械设计及理论研究所
开始
依据随机变量的分布形式和 参数，由随机样本的产生方法 产生 N 组随机样本向量的样本 xj=(xj1, xj2, …, xjn) (j=1, 2, …, N)。
将随机向量样本 xj 代人极限 状态方程，并根据状态指示函 数IF(xj)进行累加。
ˆ 求得失效概率估计值 P f ˆ 的方 估计失效概率估计值 P f 差及变异系数。
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
2. Monte Carlo 失效概率估计值的方差分析 由式(1)可以看出，失效概率估计值的随机样本xj (j=1, 2, …, N) 的函数，因此也是一个随机变量。为了对计算出的的收敛性 有一个清楚的认识，有必要对的方差进行分析。 对式(1)两边求数学期望，可得失效概率估计值 P ˆ 的期望 E P ˆ f f 如下所示：
2 机械可靠性
假定机械系统的强度随机变量 xR，应力随机变量为 xS，机械 系统极限状态函数（功能函数）为
Z g ( x R , x S ) x R xS
极限状态的定义为：整个机械系统的一部分超过某一特定状 态就不能满足设计规定的某一功能要求，此特定状态为该功 能的极限状态。
0 失效状态 Z g ( xR , xS ) 0 极限状态 0 可靠状态
东北大学机械设计及理论研究所
2 机械可靠性
f(x)
fS(xS) fR(xR) 应力xS落在附近dxS小区间内的概率为
P xS dxS f S ( xS )dxS
在应力 xS 落在小区间dxS 内的条件下，
强度xR小于应力xS的概率为
xS dxS
x
P xR xS xS dxS
3 Monte Carlo 可靠性分析
N 1 ˆ I (x ) P f F j N j 1
IF(xj)方差的无偏估计
1 N 2 2 S I ( x ) NI F j F N 1 j 1
2 2 N 1 N 2 1 N I F ( x j ) I F ( xk ) N 1 N k 1 N j 1
3 Monte Carlo 可靠性分析
3.1 Monte Carlo 模拟的理论基础 根据概率论的大数定律，若有来自同一母体且有相同分布 的n个相互独立的随机样本 x1, x2, …, xn, 他们具有相同的有 限均值 μ 和方差 σ2，则对于任意的ε >0, 有
1 n lim P xi 源自文库 1 n n i 1
东北大学机械设计及理论研究所
2 机械可靠性
Pf P ( Z 0) f R (r ) f S ( s )drds
r s xS f S ( xS ) f R ( xR )dxR dxS
FR ( xS ) f S ( xS )dxS
f (x )
机械系统 构件分布参数 fS(xS)
fR(xR)
机械系统 性能分布参数
应力分析 应力分布参数
强度分析
强度分布参数
x
机械系统 构件分布参数 应力是指对产品功能有影响的各种因素，比如机械应力、变形、磨损等。 强度是指产品承受应力的能力，比如机械强度、刚度、抗裂度等。
东北大学机械设计及理论研究所
N 1 ˆ E Pf E I F ( x j ) N j 1 样本xj与母体 N x独立同分布 1 ˆ E P E I F ( x j ) E I F ( x j ) Pf f N j 1
上式表明事件发生的频率是依概率收敛于事件发生的概率的。
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
3.2 Monte Carlo 可靠性分析的原理与计算公式 1. Monte Carlo 方法求解失效概率估计值的计算公式 失效概率可以改写成下式所示的失效域指示函数 IF(x)的数学期 望形式：
同理，也可以求得失效概率的另一种表达式
Pf P(Z 0)
f R ( xR ) f S ( xS )dxS dxR 1 FS ( xR ) f R ( xR )dxR xR
具有n设计变量的机械系统的功能函数可以表示为
xj与母体x独立同分布
1 ˆ Var Pf Var I F ( x) N
(2)
由于样本方差依概率收敛于母体的方差，所以可以用IF(.)的样
本方差
N 1 2 2 2 S I F ( x j ) NI F N 1 j 1
东北大学机械设计及理论研究所
Pf
Rn
g ( x ) 0
f X ( x1 , x2 ,
, xn )dx1dx2 , xn )dx1dx2
dxn dxn E[ I F ( x )]
I F ( x ) f X ( x1 , x2 ,
1, x F 式中， I F ( x) 为失效域指示函数；Rn为n维变量空间； 0, x F
授课教师：黄贤振
东北大学机械设计及理论研究所
1 机械系统
机械系统是指由多个机械基本要素（如机械零、部件， 机构，机器等）组成的，可以完成所需动作 ( 或动作过 程)，实现动作的传递和力的变换以及机械能的转化和利 用的装置。
机械系统
E I F ( x ) I F
N 1 ˆ I I (x ) P ˆ E P f F F j f N j 1
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
失效概率估计值的方差可以通过对式(1)两边求方差得如下：
x j 独立 N N 1 1 ˆ Var Var P Var I F ( x j ) I F ( x j ) 2 f N N j 1 j 1
Z g ( x ) g ( x1 , x2 , ..., xn )
则极限状态方程g(x1, x2, …, xn)=0将结构的基本随机变量空 间分为失效域和可靠区域两部分。
东北大学机械设计及理论研究所
2 机械可靠性
失效概率Pf可表示为
Pf
g ( x )0
f X ( x1 , x2 ,
3 Monte Carlo 可靠性分析
指示函数IF(x)方差的无偏估计可以进一步表达为
2 N N 1 N 2 N 1 1 2 2 Var I F ( x ) I F ( x j ) I F ( xk ) I F ( x j ) NI F N 1 j 1 N k 1 N 1 N j 1 ˆ P ˆ2) N ( P N 1 N 2 f f ˆ I F ( x j ) Pf N 1 N j 1 N 1 ？
xj与母体x独立 同分布
ˆ 1 Var I ( x ) Var P F j f N
ˆ Var P f
1 ˆ P ˆ2) (P f f N 1
ˆ ] ˆ Var[ P 1 P f f ˆ Cov( Pf ) ˆ ] ˆ E[ P ( N 1) P f f
E[.]为数学期望算子。
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
失效概率为失效域指示函数的数学期望，依据大数定律，失 效域指示函数的数学期望可以有失效域指示函数的样本均值 来近似。
N Nf 1 ˆ Pf I F ( x j ) N j 1 N
(1)
即 随 机 变 量 的 联 合 概率密度 函数 fX(x) 抽 取 N 个 样本 xj(j=1, 2, …, N), 落入失效域F内样本点的个数Nf与总样本点的个数 N之比即为失效概率的估计值。
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
ˆ P ，即为Pf的无偏估计 由上式可知， E P f f
在数值模拟的过程中，以指示函数IF(x)的样本均值 I F 近似代替 E[IF(x)], 则失效概率估计值的期望可以近似表达为
N 1 ˆ E I F ( x j ) E I F ( x j ) Pf P f N E j 1
将上式代人式(2)，可得失效概率估计值的方差估计为 1 ˆ ˆ P ˆ2) Var Pf (P f f N 1 进而得到估计值的变异系数为
ˆ ] ˆ Var[ P 1 P f f ˆ Cov( Pf ) ˆ ] ˆ E[ P ( N 1) P f f
东北大学机械设计及理论研究所
xS
f R ( xR )dxR
应力xS落在小区间dxS内，同时强度xR小于应力xS的概率为
P xR xS , xS dxS f S ( xS )dxS
xS
f R ( xR )dxR
根据可靠度的定义，对于应力xS所有的可能值强度 xR均小于应力xS的概率， 及事件(xR<xS)的边缘概率就是机械系统的可靠度(对上式的求和)
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
Monte Carlo 可靠性分析方法又称随机抽样法、概率模拟法 或统计试验法。该方法是通过随机模拟或者说统计试验来 进行结构可靠性分析的。由于它是以概率和数理统计理论 为基础的，故被无理学家以赌城Monte Carlo来命名。
东北大学机械设计及理论研究所
1 n 上式表明样本均值 n xi 是依概率收敛于母体的均值 i 1
μ 的。
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
另外，设随机事件 A 发生的概率 P(A), 在n次独立试验中， 事件A发生的频率为m，则随机事件A发生的频率W(A)=m/n, 对于任意ε>0，有
m lim P P( A) 1 n n
, xn )dx1dx2
dxn
式中，f(x1, x2, …, xn) 是基本随机变量x=[x1, x2, …, xn]T的联合
概率密度函数。若基本随机变量是相互独立的，则有
Pf
g ( x )0
f X1 ( x1 ) f X 2 ( x2 )
f X n ( xn )dx1dx2
dxn
式中, f(xi) (i=1, 2, …, n)为随机变量xi的概率密度函数。
多体系统是指由多个 刚性或柔性体通过一 定的联接方式互相关 联起来的完成预期动 作的一个整体
若干元件组成的承受 外部作用力并有特定 功能的整体
多
体 机
械
结 构
系
统
统 系
由构件和运动副组成 的机械系统
机
构
系
统
图1 机械系统的分类
东北大学机械设计及理论研究所
2 机械可靠性
机械可靠性模型 载荷 分布参数
设置抽样样本数N， 初始值m=0, j=0 j=j+1
由概率密度函数fX(x)参数随机样本点xj
将样本xj 代人功能函数，计算功能函数值g(xj)
IF(xj)=0
否 g(xj)≤0 ?
是 IF(xj)=1
N 1 ˆ Var Var P I ( x ) F j f N j 1
xj相互独立
N 1 ˆ Var P Var I F ( x j ) f N2 j 1
ˆ P ˆ2) N 1 N N (P 2 f f ˆ I F ( x j ) Pf N 1 N j 1 N 1
东北大学机械设计及理论研究所
开始
依据随机变量的分布形式和 参数，由随机样本的产生方法 产生 N 组随机样本向量的样本 xj=(xj1, xj2, …, xjn) (j=1, 2, …, N)。
将随机向量样本 xj 代人极限 状态方程，并根据状态指示函 数IF(xj)进行累加。
ˆ 求得失效概率估计值 P f ˆ 的方 估计失效概率估计值 P f 差及变异系数。
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
2. Monte Carlo 失效概率估计值的方差分析 由式(1)可以看出，失效概率估计值的随机样本xj (j=1, 2, …, N) 的函数，因此也是一个随机变量。为了对计算出的的收敛性 有一个清楚的认识，有必要对的方差进行分析。 对式(1)两边求数学期望，可得失效概率估计值 P ˆ 的期望 E P ˆ f f 如下所示：
2 机械可靠性
假定机械系统的强度随机变量 xR，应力随机变量为 xS，机械 系统极限状态函数（功能函数）为
Z g ( x R , x S ) x R xS
极限状态的定义为：整个机械系统的一部分超过某一特定状 态就不能满足设计规定的某一功能要求，此特定状态为该功 能的极限状态。
0 失效状态 Z g ( xR , xS ) 0 极限状态 0 可靠状态
东北大学机械设计及理论研究所
2 机械可靠性
f(x)
fS(xS) fR(xR) 应力xS落在附近dxS小区间内的概率为
P xS dxS f S ( xS )dxS
在应力 xS 落在小区间dxS 内的条件下，
强度xR小于应力xS的概率为
xS dxS
x
P xR xS xS dxS
3 Monte Carlo 可靠性分析
N 1 ˆ I (x ) P f F j N j 1
IF(xj)方差的无偏估计
1 N 2 2 S I ( x ) NI F j F N 1 j 1
2 2 N 1 N 2 1 N I F ( x j ) I F ( xk ) N 1 N k 1 N j 1
3 Monte Carlo 可靠性分析
3.1 Monte Carlo 模拟的理论基础 根据概率论的大数定律，若有来自同一母体且有相同分布 的n个相互独立的随机样本 x1, x2, …, xn, 他们具有相同的有 限均值 μ 和方差 σ2，则对于任意的ε >0, 有
1 n lim P xi 源自文库 1 n n i 1
东北大学机械设计及理论研究所
2 机械可靠性
Pf P ( Z 0) f R (r ) f S ( s )drds
r s xS f S ( xS ) f R ( xR )dxR dxS
FR ( xS ) f S ( xS )dxS
f (x )
机械系统 构件分布参数 fS(xS)
fR(xR)
机械系统 性能分布参数
应力分析 应力分布参数
强度分析
强度分布参数
x
机械系统 构件分布参数 应力是指对产品功能有影响的各种因素，比如机械应力、变形、磨损等。 强度是指产品承受应力的能力，比如机械强度、刚度、抗裂度等。
东北大学机械设计及理论研究所
N 1 ˆ E Pf E I F ( x j ) N j 1 样本xj与母体 N x独立同分布 1 ˆ E P E I F ( x j ) E I F ( x j ) Pf f N j 1
上式表明事件发生的频率是依概率收敛于事件发生的概率的。
东北大学机械设计及理论研究所
3 Monte Carlo 可靠性分析
3.2 Monte Carlo 可靠性分析的原理与计算公式 1. Monte Carlo 方法求解失效概率估计值的计算公式 失效概率可以改写成下式所示的失效域指示函数 IF(x)的数学期 望形式：
同理，也可以求得失效概率的另一种表达式
Pf P(Z 0)
f R ( xR ) f S ( xS )dxS dxR 1 FS ( xR ) f R ( xR )dxR xR
具有n设计变量的机械系统的功能函数可以表示为
xj与母体x独立同分布
1 ˆ Var Pf Var I F ( x) N
(2)
由于样本方差依概率收敛于母体的方差，所以可以用IF(.)的样
本方差
N 1 2 2 2 S I F ( x j ) NI F N 1 j 1
东北大学机械设计及理论研究所
Pf
Rn
g ( x ) 0
f X ( x1 , x2 ,
, xn )dx1dx2 , xn )dx1dx2
dxn dxn E[ I F ( x )]
I F ( x ) f X ( x1 , x2 ,
1, x F 式中， I F ( x) 为失效域指示函数；Rn为n维变量空间； 0, x F