【《基于蒙特卡罗的某隧道工程风险因子条件概率分析计算过程案例》6700字（论文）】

基于蒙特卡罗的某隧道工程风险因子条件概率分析计算过程案例目录TOC\o"1-3"\h\u151751蒙特卡罗基本概念 1128382风险因子之间理论关系 2241813风险因子条件概率分析 61蒙特卡罗基本概念蒙特卡罗方法，也称为统计模拟方法，是一种基于电子计算的以数学概率理论为指导的重要研究方法。该方法起源于十八世纪法国数学家提出的投针实验，通过统计针与平行线相交的情况，推算圆周率的估计值。二十世纪以来，随着计算机的投入使用，蒙特卡罗方法不仅在数学、物理等领域得以应用，而且在还广泛用于金融工程学、宏观经济学及计算物理学等领域。蒙特卡罗方法主要针对于实际中难以直接处理的问题，可以通过实验模拟的方式，基于实验结果发生特征作为该问题的解。该方法解答问题的过程主要分为四个主要步骤，①基于问题本身性质，构造或转化为随机性质的问题，建立对应的概率模型；②基于概率模型，产生符合概率分布的随机变量，即随机抽样；③将随机抽样获取的值代入推算流程，得出一个无偏估计，作为问题的一次解；④重复以上操作，获取大量估计量，采用统计的方法反映该问题可能发生的情况。蒙特卡罗方法中最关键步骤，即建立适合问题且便于实现的概率模型，常见的连续性概率模型有三类，分别是均匀分布式、正态分布式和指数分布式。（1）均匀分布式均匀分布式中，其产生的随机变量落在规定区间或子区间内的概率只和区间长度有关，与区间位置无关，即均匀分布式产生的随机变量具有均匀性。其概率密度函数如式4-1所示。（4-1）（2）正态分布式正态分布式中，其随机变量的产生概率服从正态分布函数，主要由期望值和标准差控制，其概率密度函数如式4-2所示。（4-2）（3）三角分布式三角分布式中，其随机变量的产生概率服从三角分布函数，主要由最小值a、中间值b和最大值c控制，其概率密度函数如式4-3所示。（4-3）（4）指数分布式指数分布式是一种常见的分布模型，符合人类寿命、设备寿命规律特征，又如粒子运输过程中路径轨迹的分布也符合指数分布[58]。该抽样方法中，产生的随机变量概率服从指数分布函数，如式4-4所示。（4-4）2风险因子之间理论关系（1）风管阻力与其父节点风险因子理论关系风管阻力主要由风管沿程阻力和局部阻力组成，其中沿程阻力是由于空气自身粘滞性及其与风管管壁之间摩擦产生的沿程风压损失，局部阻力是由于风管变向或变径导致边界层分离产生旋涡而造成的风压损失。风管阻力计算公式如下：（4-5）式中：Pf——风管阻力；λ——摩擦阻力系数，λ=0.02；L——风管长度，m；D——风管直径，m；v——风管进风口风速，m/s；ξ——局部阻力系数，风管进口处ξ=0.6，风管出口处ξ=1，风管弯曲处（2）风管漏风率与其父节点风险因子理论关系通风风管由于材料或制作工艺原因，其密闭性能欠佳，在隧道施工过程中难免会出现破损情况，进而造成新风泄漏[56]。风管内百米漏风率计算公式如4-6所示。（4-6）式中：β——风管百米漏风率；Qm——风管单位破损孔口新风漏风量；N——风管破损数目；Q0——风管进风口风量，m3/s；i——风管更换程度，完全未更换情况下i=0%，完全更换情况下i=100%，其他情况介于两者之间。其中，单位孔口新风漏风量Qm与风管内风速vm、孔口面积f、风管孔口平面与漏风风流之间的夹角φ相关。单位孔口漏风量计算如公式5-4所示。（4-7）式中：μ——风管破损孔口流量系数，μ=0.63；vm——风管孔口处新风流速，m/s；f——风管孔口面积，m2；φ——风管孔口平面与漏风风流之间的夹角，°。（3）风管末端出风量与其父节点风险因子理论关系风管末端出风量主要受风阻、百米漏风率和通风长度的控制，其中风阻通过对流体的能量影响进而影响到风量。风管内风压主要分为静压与动压两个部分，静压是用于克服风阻的一部分能量，另一部分能量表现动能形式。风管内风压计算如公式4-8所示。（4-8）式中：Pt——通风全压，Pa；Ps——通风静压，Pa；Pv——通风动压，Pa。假定风管不存在漏风情况，那么风管内存在能量守恒关系，如公式4-9所示。（4-9）式中：Ps1、Pv1——断面1位置下的静压与动压，Pa；Ps2、Pv2——断面2位置下的静压与动压，Pa；Pf——断面1至断面2之间的风阻，Pa。其中动压计算公式如下所示。（4-10）在此情况下，结合5-3与5-4公式，可以推导出在无漏风情况下风管末端出风量计算公式，如公式4-11所示。（4-11）式中：Qs——风管末端出风量，m3/s；Q0——风管进风量，m3/s；A——风管断面面积，m2；Pf——风管风阻，Pa；ρ——风流密度，kg/m3。当风管沿程存在漏风情况，那么风管末端出风量应该如公式4-12表示。（4-12）（4）工作面需风量与其父节点风险因子理论关系工作面需风量主要考虑施工人员生理需求、掘进爆破工序反应需求、内燃机设备运转需求与隧道施工通风标准要求四个方面，计算如公式4-13所示。（4-13）式中：Qn——隧道工作面需风量，m3/s；Kr——重率高程校正系数，即海拔高度H处的空气重率与海平面处的空气重率之比，其确定如公式4-14所示；k——隧道工作面单位施工人员需风量，k=0.05m3/(人.s)；n——隧道工作面施工人员最大数目，人；G——单次爆破用药量，kg；A——隧道开挖面积，m2；L0——爆破排烟长度，其确定如公式4-15所示；H——内燃机械设备功率，kW；q——内燃机单位功率需风量，平原隧道环境下q=0.05m3/(kW.s)，高原隧道环境下q=0.075m3/(kW.s)；vmin——隧道洞内最低风速，正洞内要求不低于0.15m/s。（4-14）式中：γz——海拔高度为H处的空气重率，N/m3；γ0——海拔高度为0处的空气重率，N/m3；H——海拔高度，秦岭天台山隧道项目平均海拔1500m。（4-15）（5）新风量供需关系与其父节点风险因子理论关系新风量供需关系可以理解为工作面处新风量与需风量之间关系，然而风管末端出风量并非工作面新风量，但两者之间处在一定函数关系，而且该关系受风管末端至工作面距离控制。为探究工作面新风量与风管末端出风量之间关系，论文拟采用正交试验结合数值仿真模拟的方式实现。秦岭天台山隧道施工断面为马蹄形，断面面积98.3m，通风风管采用1800mm，论文将其简化为相等面积的圆形，风管布设于中心位置。正交试验的控制变量为风管末端至工作面距离l，距离在30~70m之间，每5m进行一次数值仿真模拟，分别进行了9组模拟。模拟过程中，其余参数保持不变，其中风管末端出风速10m/s（25.45m3/s），模拟结果如图4-1所示。图4-1(a)l=30m风速云图图4-1(b)l=35m风速云图图4-1(c)l=40m风速云图图4-1(d)l=45m风速云图图4-1(e)l=50m风速云图图4-1(f)l=55m风速云图图4-1(g)l=60m风速云图图4-1(h)l=65m风速云图图4-1(i)l=70m风速云图图4-1风管末端至工作面不同距离下隧道水平断面风速云图结合4-1速度云图分析可知，随着风管末端至工作面距离的增大，新风的扩散逐渐紊乱，冲击能力相对降低。另一方面，工作面位置处风量也随着末端至工作面距离有所不同，具体情况如图4-2所示。图4-2隧道工作面位置风速与风管末端至工作面距离数据模拟结果指出，随着风管末端至工作面距离的增大，工作面风量先呈现增大趋势再呈现降低趋势，且两个趋势均是线性关系。为拟合出两者之间关系，将风管末端至工作面为50m作为标准位置，取风管末端至工作面距离与最优距离的风量差∆q作为因变量，取两距离差∆l作为自变量，求解函数关系。拟合结果如公式4-16所示，展开结果如公式4-17所示。（4-16）式中：∆q——风管末端至工作面距离与最优距离间的工作面风量差，m3/s；qs——风管末端至工作面最优距离下工作面风量，m3/s；∆l——风管末端至工作面距离与最优距离间的距离差，m；ls——风管末端至工作面位置的最优距离，m。（4-17）式中：ql——风管末端至工作面距离为l下工作面风量，m3/s；l——风管末端至工作面位置的距离，m。根据相似准则，相同断面下断面新风量与整体风量之间存在公式4-18关系。（4-18）式中：Ql——风管末端至工作面距离为l下工作面新风风量，m3/s；Qs——风管末端至工作面最优距离下工作面新风风量，m3/s。因而可以推导出隧道工作面位置新风供需风量关系，如公式4-19所示。（4-19）式中：∆Q——隧道工作面位置通风供需风量差，m3/s。3风险因子条件概率分析基于蒙特卡罗的风险因子概率分析，主要是通过编程软件JAVA实现，通过上万次仿真实验计算出不符合标准要求的频次，进而获得风险因子的条件概率。仿真实验流程如图4-3所示。模拟实验的关键在于计算模型的形成、父节点风险因子状态范围的概率化与标准状态的识别，其中计算模型即是风险因子与其父节点风险因子之间的理论关系。父节点风险因子状态范围采用三角分布式概率模型模拟，三角分布式适用于样本数据较少且已知最小状态、平均状态和最大状态的情况下，适用于隧道施工通风风险概率分析。风险因子的标准状态，即风险因子发生与不发生两种状态的边界状态，在隧道施工通风风险研究中该状态定义为实际状态。图4-3基于蒙特卡罗的仿真实验流程（1）风管阻力条件概率分析基于风管阻力与其父节点风险因子理论关系，结合仿真实验流程，对其条件概率进行了分析，分析按照正交试验进行。根据理论关系，控制条件分为三类，通过不同组合形式形成八种，如表4-1所示。其中标准状态：风管直径1.8m，风管弯曲3处，每一处弯曲角度30°，风管长度2000m，标准状态下风管阻力为588Pa。表4-1风管阻力仿真正交试验控制条件表组合编号控制条件状态（Y/N）控制条件（a,b,c）IF105=Y,F207=Y,F209=Y风管直径(1.2m,1.5m,1.8m)；风管弯曲(3处,6处,12处)，每一处弯曲角度(15°,30°,45°)；风管长度(2000m,2760m,3520m)。IIF105=Y,F207=Y,F209=N风管直径(1.2m,1.5m,1.8m)；风管弯曲(3处,6处,12处)，每一处弯曲角度(15°,30°,45°)；风管长度2000m。IIIF105=Y,F207=N,F209=Y风管直径(1.2m,1.5m,1.8m)；风管弯曲3处，每一处弯曲角度(15°,30°,45°)；风管长度(2000m,2760m,3520m)。IVF105=Y,F207=N,F209=N风管直径(1.2m,1.5m,1.8m)；风管弯曲3处，每一处弯曲角度(15°,30°,45°)；风管长度2000m。VF105=N,F207=Y,F209=Y风管直径1.8m；风管弯曲(3处,6处,12处)，每一处弯曲角度(15°,30°,45°)；风管长度(2000m,2760m,3520m)。VIF105=N,F207=Y,F209=N风管直径1.8m；风管弯曲(3处,6处,12处)，每一处弯曲角度(15°,30°,45°)；风管长度2000m。VIIF105=N,F207=N,F209=Y风管直径1.8m；风管弯曲3处，每一处面积每一处弯曲角度(15°,30°,45°)；风管长度(2000m,2760m,3520m)。VIIIF105=N,F207=N,F209=N风管直径1.8m；风管弯曲3处，每一处面积每一处弯曲角度(15°,30°,45°)；风管长度2000m。基于以上组合，分别进行了蒙特卡罗仿真实验，实验结果如图4-4所示。图4-4(a)组合I下风管阻力频次图图4-4(b)组合II下风管阻力频次图图4-4(c)组合III下风管阻力频次图图4-4(d)组合IV下风管阻力频次图图4-4(e)组合V下风管阻力频次图图4-4(f)组合VI下风管阻力频次图图4-4(g)组合VII下风管阻力频次图图4-4(h)组合VIII下风管阻力频次图图4-4各控制条件组合下风管阻力频次图上图中红线表示风管阻力在标准状态下的数值，大于该红线部分面积与整体面积之比即是该控制条件组合下的条件概率，整理结果如表4-2所示。表4-2各控制条件组合下的风管阻力条件概率组合编号父节点风险因子状态P(F404=Y)/%P(F404=N)/%IF105=Y,F207=Y,F209=Y1000IIF105=Y,F207=Y,F209=N1000IIIF105=Y,F207=N,F209=Y1000IVF105=Y,F207=N,F209=N99.990.01VF105=N,F207=Y,F209=Y1000VIF105=N,F207=Y,F209=N97.302.70VIIF105=N,F207=N,F209=Y1000VIIIF105=N,F207=N,F209=N39.9360.07（2）风管百米漏风率条件概率分析同理，对风管百米漏风率条件概率进行分析，其控制条件组合形式分为八种，如表4-3所示。其中风管百米漏风率标准状态为2%。表4-3风管百米漏风率仿真正交试验控制条件表组合编号控制条件状态（Y/N）控制条件（a,b,c）IF209=Y,F410=Y,F504=Y风管长度(2000m,2760m,3520m)；风管破损(5处,22处,44处)，每处破损面积(1cm2,5cm2,10cm2)；风管更换程度0%。IIF209=Y,F410=Y,F504=N风管长度(2000m,2760m,3520m)；风管破损(5处,22处,44处)，每处破损面积(1cm2,5cm2,10cm2)；风管更换程度(0%,5%,10%)。IIIF209=Y,F410=N,F504=Y风管长度(2000m,2760m,3520m)；风管破损5处，每处破损面积(1cm2,5cm2,10cm2)；风管更换程度0%。IVF209=Y,F410=N,F504=N风管长度(2000m,2760m,3520m)；风管破损5处，每处破损面积(1cm2,5cm2,10cm2)；风管更换程度(0%,5%,10%)。VF209=N,F410=Y,F504=Y风管长度2000m；风管破损(5处,22处,44处)，每处破损面积(1cm2,5cm2,10cm2)；风管更换程度0%。VIF209=N,F410=Y,F504=N风管长度2000m；风管破损(5处,22处,44处)，每处破损面积(1cm2,5cm2,10cm2)；风管更换程度(0%,5%,10%)。VIIF209=N,F410=N,F504=Y风管长度2000m；风管破损5处，每处破损面积(1cm2,5cm2,10cm2)；风管更换程度0%。VIIIF209=N,F410=N,F504=N风管长度2000m；风管破损5处，每处破损面积(1cm2,5cm2,10cm2)；风管更换程度(0%,5%,10%)。基于以上组合，分别进行了蒙特卡罗仿真实验，实验结果如图4-5所示。图4-5(a)组合I下风管百米漏风率频次图图4-5(b)组合II下风管百米漏风率频次图图4-5(c)组合III下风管百米漏风率频次图图4-5(d)组合IV下风管百米漏风率频次图图4-5(c)组合V下风管百米漏风率频次图图4-5(d)组合VI下风管百米漏风率频次图图4-5(g)组合VII下风管百米漏风率频次图图4-5(h)组合VIII下风管百米漏风率频次图图4-5各控制条件组合下风管百米漏风率频次图上图中红线表示风管阻力在标准状态下的数值，大于该红线部分面积与整体面积之比即是该控制条件组合下的条件概率，整理结果如表4-4所示。表4-4各控制条件组合下的风管百米漏风率条件概率组合编号父节点风险因子状态P(F403=Y)/%P(F403=N)/%IF209=Y,F410=Y,F504=Y99.410.59IIF209=Y,F410=Y,F504=N39.4960.51IIIF209=Y,F410=N,F504=Y31.5265.48IVF209=Y,F410=N,F504=N5.7891.22VF209=N,F410=Y,F504=Y99.930.07VIF209=N,F410=Y,F504=N71.2025.80VIIF209=N,F410=N,F504=Y70.8229.18VIIIF209=N,F410=N,F504=N25.5271.48（3）风管末端出风量条件概率分析同理，对风管末端出风量条件概率进行分析，其控制条件组合形式分为八种，如表4-5所示。根据通风设计计算，风管末端出风量标准状态为31.2m3/s。表4-5风管末端出风量仿真正交试验控制条件表组合编号控制条件状态（Y/N）控制条件（a,b,c）IF209=Y,F403=Y,F404=Y风管长度(2000m,2760m,3520m)；风管百米漏风率(0%,2%,10%)；风管阻力(500Pa,2000Pa,4000Pa)。IIF209=Y,F403=Y,F404=N风管长度(2000m,2760m,3520m)；风管百米漏风率(0%,2%,10%)；风管阻力588Pa。IIIF209=Y,F403=N,F404=Y风管长度(2000m,2760m,3520m)；风管百米漏风率2%；风管阻力(500Pa,2000Pa,4000Pa)。IVF209=Y,F403=N,F404=N风管长度(2000m,2760m,3520m)；风管百米漏风率2%；风管阻力588Pa。VF209=N,F403=Y,F404=Y风管长度2000m；风管百米漏风率(0%,2%,10%)；风管阻力(500Pa,2000Pa,4000Pa)。VIF209=N,F403=Y,F404=N风管长度2000m；风管百米漏风率(0%,2%,10%)；风管阻力588Pa。VIIF209=N,F403=N,F404=Y风管长度2000m；风管百米漏风率2%；风管阻力(500Pa,2000Pa,4000Pa)。VIIIF209=N,F403=N,F404=N风管长度2000m；风管百米漏风率2%；风管阻力588Pa。基于以上组合，分别进行了蒙特卡罗仿真实验，实验结果如图4-6所示。图4-6(a)组合I下风管末端出风量频次图图4-6(b)组合II下风管末端出风量频次图图4-6(c)组合III下风管末端出风量频次图图4-6(d)组合IV下风管末端出风量频次图图4-6(e)组合V下风管末端出风量频次图图4-6(f)组合VI下风管末端出风量频次图图4-6(g)组合VII下风管末端出风量频次图图4-6(h)组合VIII下风管末端出风量频次图图4-6各控制条件组合下风管末端出风量频次图上图中红线表示风管末端出风量在标准状态下的数值，大于该红线部分面积与整体面积之比即是该控制条件组合下的条件概率，整理结果如表4-6所示。表4-6各控制条件组合下的风管末端出风量条件概率组合编号父节点风险因子状态P(F405=Y)/%P(F405=N)/%IF209=Y,F403=Y,F404=Y85.4611.54IIF209=Y,F403=Y,F404=N81.9015.10IIIF209=Y,F403=N,F404=Y89.9310.07IVF209=Y,F403=N,F404=N89.3610.64VF209=N,F403=Y,F404=Y71.1325.87VIF209=N,F403=Y,F404=N72.9427.06VIIF209=N,F403=N,F404=Y0100VIIIF209=N,F403=N,F404=N0100（4）隧道工作面需风量条件概率分析同理，对隧道工作面需风量条件概率进行分析，其控制条件组合形式分为四种，如表4-7所示。其中标准状态：爆破用药量375kg，内燃机设备2台，总功率为322kW，标准状态下隧道工作面需风量为26.21m3/s。表4-7隧道工作面需风量仿真正交试验控制条件表组合编号控制条件状态（Y/N）控制条件（a,b,c）IF411=Y,F413=Y爆破用药量(350kg,375kg,400kg)；内燃机设备数目(1台,2台,4台)，每台功率(107kW,161kW,215kW)。IIF411=Y,F413=N爆破用药量(375kg,380kg,390kg)；内燃机设备数目0台。IIIF411=N,F413=Y爆破用药量375kg；内燃机设备数目(1台,2台,4台)，每台功率(107kW,161kW,215kW)。。IVF411=N,F413=N爆破用药量375kg；内燃机设备数目0台。基于以上组合，分别进行了蒙特卡罗仿真实验，实验结果如图4-7所示。图4-7(a)组合I下隧道工作面需风量频次图图4-7(b)组合II下隧道工作面需风量频次图图4-7(c)组合III下隧道工作面需风量频次图图4-7(d)组合IV下隧道工作面需风量频次图图4-7各控制条件组合下隧道工作面需风量频次图上图中红线表示隧道工作面需风量在标准状态下的数值，大于该红线部分面积与整体面积之比即是该控制条件组合下的条件概率，整理结果如表4-8所示。表4-8各控制条件组合下的隧道工作面需风量条件概率父节点风险因子状态P(F416=Y)/%P(F416=N)/%F411=Y,F413=Y58.4741.53F411=Y,F413=N49.8550.15F411=N,F413=Y16.7383.27F411=N,F413=N0100（5）隧道工作面供需风量关系条件概率分析同理，对隧道工作面供需风量关系条件概率进行分析，其控制条件组合形式分为八种，如表4-9所示。其中隧道工作面供需风量关系标准状态为0m3/s。表4-9隧道工作面供需风量关系仿真正交试验控制条件表组合编号控制条件状态（Y/N）控制条件（a,b,c）IF210=Y,F405=Y,F416=Y风管末端至工作面距离(30m,50m,70m)；风管末端出风量(10m3/s,30m3/s,50m3/s)；隧道工作面需风量(24m3/s,26m3/s,34m3/s)。IIF210=Y,F405=Y,F416=N风管末端至工作面距离(30m,50m,70m)；风管末端出风量(10m3/s,30m3/s,50m3/s)；隧道工作面需风量26m3/s。IIIF210=Y,F405=N,F416=Y风管末端至工作面距离(30m,50m,70m)；风管末端出风量30m3/s；隧道工