CN113505514B - 一种复杂扰动条件下岩体弹塑性损伤-渗流耦合计算方法 - Google Patents
一种复杂扰动条件下岩体弹塑性损伤-渗流耦合计算方法 Download PDFInfo
- Publication number
- CN113505514B CN113505514B CN202110891535.2A CN202110891535A CN113505514B CN 113505514 B CN113505514 B CN 113505514B CN 202110891535 A CN202110891535 A CN 202110891535A CN 113505514 B CN113505514 B CN 113505514B
- Authority
- CN
- China
- Prior art keywords
- rock mass
- damage
- elastoplastic
- disturbance
- under
- 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
Links
- 239000011435 rock Substances 0.000 title claims abstract description 147
- 230000008878 coupling Effects 0.000 title claims abstract description 63
- 238000010168 coupling process Methods 0.000 title claims abstract description 63
- 238000005859 coupling reaction Methods 0.000 title claims abstract description 63
- 238000004364 calculation method Methods 0.000 title claims abstract description 43
- 230000006378 damage Effects 0.000 claims abstract description 71
- 238000009826 distribution Methods 0.000 claims abstract description 20
- 230000001550 time effect Effects 0.000 claims abstract description 15
- 238000010276 construction Methods 0.000 claims abstract description 11
- 238000011156 evaluation Methods 0.000 claims abstract description 9
- 238000005422 blasting Methods 0.000 claims description 47
- 230000035699 permeability Effects 0.000 claims description 38
- 238000000034 method Methods 0.000 claims description 35
- 230000006870 function Effects 0.000 claims description 22
- 230000009471 action Effects 0.000 claims description 15
- 239000000463 material Substances 0.000 claims description 14
- 230000008859 change Effects 0.000 claims description 13
- 230000000694 effects Effects 0.000 claims description 10
- 238000004458 analytical method Methods 0.000 claims description 9
- 238000010257 thawing Methods 0.000 claims description 9
- 125000004122 cyclic group Chemical group 0.000 claims description 5
- 238000005315 distribution function Methods 0.000 claims description 5
- 230000007774 longterm Effects 0.000 claims description 5
- 238000009827 uniform distribution Methods 0.000 claims description 5
- 230000001186 cumulative effect Effects 0.000 claims description 4
- 238000013507 mapping Methods 0.000 claims description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 4
- 238000013459 approach Methods 0.000 claims description 3
- 238000003491 array Methods 0.000 claims description 3
- 230000008014 freezing Effects 0.000 claims description 3
- 238000007710 freezing Methods 0.000 claims description 3
- 239000011148 porous material Substances 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 238000011426 transformation method Methods 0.000 claims description 3
- 238000000518 rheometry Methods 0.000 claims 1
- 238000013461 design Methods 0.000 abstract description 4
- 230000009977 dual effect Effects 0.000 abstract description 3
- 238000006073 displacement reaction Methods 0.000 description 22
- 238000012544 monitoring process Methods 0.000 description 14
- 238000009412 basement excavation Methods 0.000 description 10
- 230000008569 process Effects 0.000 description 10
- 238000005553 drilling Methods 0.000 description 3
- 238000004062 sedimentation Methods 0.000 description 3
- 238000012360 testing method Methods 0.000 description 3
- 102100021503 ATP-binding cassette sub-family B member 6 Human genes 0.000 description 2
- 101100000375 Homo sapiens ABCB6 gene Proteins 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 2
- 238000007906 compression Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000003902 lesion Effects 0.000 description 2
- 230000003068 static effect Effects 0.000 description 2
- 208000027418 Wounds and injury Diseases 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000005094 computer simulation Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 239000003673 groundwater Substances 0.000 description 1
- 230000002706 hydrostatic effect Effects 0.000 description 1
- 208000014674 injury Diseases 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003204 osmotic effect Effects 0.000 description 1
- 238000005325 percolation Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000003313 weakening effect Effects 0.000 description 1
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
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- 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/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Mathematical Optimization (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Computer Hardware Design (AREA)
- Operations Research (AREA)
- Evolutionary Computation (AREA)
- Algebra (AREA)
- Geometry (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Investigation Of Foundation Soil And Reinforcement Of Foundation Soil By Compacting Or Drainage (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种复杂扰动条件下岩体弹塑性损伤‑渗流耦合计算方法，包括：S1、建立复杂扰动下岩体弹塑性损伤耦合模型；S2、建立复杂扰动下岩体弹塑性损伤‑渗流耦合模型；S3、将岩体弹塑性损伤模型修改为服从Weibull分布的概率表达形式；S4、建立考虑时间效应的岩体弹塑性损伤‑渗流耦合模型；S5、进行复杂扰动下考虑时间效应的岩体弹塑性损伤‑渗流耦合场进行分析，获取对应的安全性评价数据。本发明建立了同时考虑扰动因素与应力重分布双重作用下的岩体破坏准则，同时考虑了岩体的流变特性；将复杂扰动下考虑时间效应的岩体弹塑性损伤‑渗流耦合模型应用于实际工程的稳定性评价当中，为类似工程的安全施工提供了一定的设计依据。
Description
技术领域
本发明涉及隧道稳定性分析技术领域，尤其涉及一种复杂扰动条件下的岩体弹塑性损伤-渗流耦合计算方法。
背景技术
隧道开挖过程中，周围岩体不可避免的会受到扰动影响。钻爆施工过程中的爆破荷载，寒区岩体的冻融作用等均会对掌子面附近岩体的力学参数造成不同程度的影响。同时岩石的流变特性也是影响岩土工程长期稳定性的重要因素之一。大量力学试验及现场监测数据表明，岩石的强度和变形与时间有着密切的关系。此外，岩石是一种非均匀的地质材料，其物理力学性质在空间上的分布有着明显的非连续性。岩石的非均匀特性并不是静态不变的，随着时间和加载历史的改变，岩石内部原有的非均匀性和内部缺陷在外荷载作用下不断发展，进一步地增加了自身的非均匀特性。已有岩石弹塑性损伤模型大多采用Mohr-Coulomb(M-C)或Druker-Prager(D-P)等线性强度准则，广义Hoek-Brown(H-B)屈服准则更能够反映出岩体的非线性特征及结构面、开挖扰动等因素对岩体强度的影响。由于H-B准则在棱线和尖点处存在不连续性，导致其有限元数值求解过程变得十分困难。X.D.PAN、R.G.Wan和R.S.Merifield等通过将H-B准则的角点进行圆滑处理或修改屈服函数的方法，避免了数值求解中的奇异点问题。然而，等价参数法存在一定的使用范围，角点圆滑法本质上修改了屈服函数形式，使其在求解一些经典问题时会产生偏差。
除开挖造成的塑性损伤以外，岩体还会受到施工方法及环境扰动等因素的影响。一方面，岩体中的水冰相变导致岩体冻胀力，岩体微裂纹扩展；另一方面，外部荷载作用改变了岩体的应力，也会引起岩体形变和破坏的产生。在钻爆法施工过程中，频繁的爆破作业会对围岩造成较大扰动,威胁施工安全,因此,爆破振动的预测与控制也成为影响工程顺利进行的关键因素之一。此外，随着损伤变量的演化，岩体内部裂隙不断发育贯通，岩体渗透率也随之发生重大变化，应力-渗流耦合相应愈发明显。地下水充斥于边坡裂隙当中，大大降低了其抗剪强度，产生的静水压力和动水压力还会对潜在崩塌体产生托浮作用，导致岩体的整体稳定性降低。
发明内容
本发明提供一种复杂扰动下考虑时间效应的岩体弹塑性损伤-渗流计算方法，以克服上述技术问题。
本发明包括以下步骤：
S1、建立复杂扰动下岩体弹塑性耦合损伤模型，所述复杂扰动包括：冻融循环扰动、循环爆破扰动；
S2、引入渗透系数演化方程，建立复杂扰动下岩体弹塑性损伤-渗流耦合模型；
S3、将岩体弹塑性损伤模型参数进行服从Weibull分布的概率表达；
S4、建立考虑时间效应的岩体弹塑性损伤-渗流耦合模型；
S5、给定工程条件，所述工程条件即施加边界条件并输入对应的材料参数，基于步骤S1-S4进行复杂扰动作用下考虑时间效应的岩体弹塑性损伤-渗流场耦合分析，获取当前工程条件对应的安全性评价数据。
进一步地，S1中所述复杂扰动下岩体弹塑性耦合损伤模型为：
f＝σ1-σ3-σci[(mbσ3/σci+s)]a (1)
其中，σ1和σ3分别为岩体的最大主应力和最小主应力；σci为完整岩石的单轴抗压强度；mb、a为针对不同岩体的量纲为一的强度参数；s为反映岩体破碎程度的强度参数；三个强度参数的计算公式为：
其中，GSI为地质强度指标；mi为力学参数，用于反映岩石的软硬程度；D为扰动系数，反映了施工作业对岩体的扰动程度，是一个损伤结果量；
复杂扰动下岩体的耦合损伤Dc的计算公式为：
其中，Dc为复杂扰动下岩体的耦合损伤，α、β为弹塑性损伤演化方程中的常参数；N为冻融循环或循环爆破扰动次数；Nt为当前的扰动次数；Dt为当前扰动次数下的岩体损伤值；U、W为与岩体材料性质相关的常参数；为等效塑性应变；
等效塑性应变的计算公式为：
其中，εp1、εp2、εp3为笛卡尔坐标系中x、y、z三个方向上的主塑性应变。
进一步地，S2引入岩体渗透系数演化方程，建立复杂扰动下岩体弹塑性损伤-渗流耦合模型，包括：
S21、建立岩体渗透系数演化方程；
当岩体处于弹性状态时，渗透系数的演化方程为：
其中，k0为初始渗透系数；εv为体积应变；φ0为初始孔隙率。
当岩石处于塑性状态时，渗透系数的演化方程为：
其中，ξ为跳突系数，用于描述泥岩破裂时渗透性突增；A'和B'分别为常参数，其表达式为A'＝1/(e-1/θ-1)，B'＝-1/(e-1/θ-1)；θ为经验参数。
S22、建立岩体孔隙度演化方程如下：
其中，b和R为经验参数，用于描述损伤变量对孔隙度的影响，且R>1。
进一步地，S3将岩体弹塑性损伤模型进行服从Weibull分布的概率化表达：
S31、第i个RVE力学参数Pi的计算公式为，所述RVE为模拟代表性体积元：
Pi＝P0·xi (10)
其中，P0为岩石式样的宏观性质；xi为相应的Weibull随机数。
双参数Weibull分布下，x的累积分布函数W(x)为：
x的概率密度函数w(x)为：
x的期望值为：
x的方差为：
其中，x是大于零的服从Weibull分布的随机数，用于表征岩体的非均质性，x0为随机数的平均值；xu为x的阈值；ζ为形状参数，即均质性指数；x0为尺寸参数，与数学期望有关；Г为Gamma函数；E(x)为数学期望函数；Var(x)为方差函数。
S32、根据Monte-Carlo随机抽样方法，对各个单元进行初始化赋值；岩石的力学参数mb符合Weibull分布规律，其余参数均为统一值，建立累积分布函数：
在(0,1)范围内生成n个服从均匀分布的随机数u，所述n为计算时模型划分的单元数，用逆变法将这组随机数生成为服从Weibull分布的随机数，映射函数为：
mb＝mb0(-ln(u))1/ζ (16)
其中，mb0为强度参数mb的平均值，u为在(0,1)范围内生成n个服从均匀分布的随机数，n为计算时模型划分的单元数；
将生成的n个服从Weibull分布的随机数，通过PYTHON语言中的循环语句生成n个材料数组，对n个单元进行遍历赋值。
进一步地，S4中建立考虑时间效应的岩体弹塑性损伤-渗流耦合模型：
对于RVE单元而言，力学参数P随时间的变化规律为：
Pt＝P∞+(P0-P∞)e-At (17)
其中，P为P0为力学参数的初始值；P∞为强度参数的长期强度值，当t趋近于无限大时，P∞＝P0；A为强度衰减因子，表示强度衰减的速率。
进一步地，S5中给定工程条件，所述工程条件即施加边界条件并输入对应的材料参数，基于步骤S1-S4进行复杂扰动作用下考虑时间效应的岩体弹塑性损伤-渗流场耦合分析，包括：
S51、根据耦合损伤Dc和流变时间t，对初始强度参数进行修正；
S52、对所述岩体弹塑性损伤模型进行弹性预测，求解对应的预测应力；求解对应的预测应力的公式为：
其中，tn为计算时步，σn为tn时刻的应力；为tn时刻的内变量；Δεn+1为tn+1时刻的应变增量；
S53、将预测应力代入S1中所述复杂扰动下岩体弹塑性耦合损伤模型，判断所述损伤模型是否大于零；若大于零，对预测应力进行塑性修正，对预测应力进行塑性修正的计算公式为：
其中，Δλ为塑性因子增量；h为塑性势函数对应力的导数，g为塑性势函数；Δσp为塑性应力增量，Δσp＝Δλh；
S54、根据损伤变量，对应力再次进行修正：
其中，σ'n+1为计算所得的最终应力；
S55、根据应力场计算所得的损伤值和塑性应变值、初始渗透系数和孔隙率、当前荷载步下的渗透系数和孔隙率、以及通过渗流有限元求解器对孔隙水压p进行求解，对最终应力再次进行修正，对最终应力进行修正的公式为：
其中，为有效应力；α为Biot系数；δij为克罗内克尔符号。
S56、重复步骤S51至S55，直到两时步之间的计算差值小于预设收敛值为止，所述收敛值的取值范围为1e10-5～1e10-3。
本发明从Hoek-Brown准则中的扰动系数出发，将循环爆破和冻融循环造成的扰动视为疲劳损伤，建立了同时考虑扰动因素与应力重分布双重作用下的岩体破坏准则；进一步的，引入渗透系数演化方程，建立岩体弹塑性损伤-渗流耦合模型；将模型中的力学参数进行；将复杂扰动下岩体应力-损伤-渗流耦合模型应用于实际工程的稳定性评价当中，为类似工程的安全施工提供了一定的设计依据。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明的整体流程图；
图2为本发明复杂扰动下岩体弹塑性损伤-渗流耦合计算流程图；
图3为本发明单轴压缩有限元计算模型图；
图4为本发明损伤单元加载时步-应变-渗透系数曲线图；
图5为本发明渗透系数与突跳系数ξ的关系图；
图6为本发明渗透系数与参数θ的关系图；
图7为本发明孔隙度与参数b的关系图；
图8为本发明孔隙度与参数R的关系图；
图9为本发明扰动作用下岩石耦合损伤模型；
图10为本发明用于某爆破开挖地铁车站稳定性评价中的计算模型图；
图11为本发明不同爆破次数下高边墙水平位移；
图12为本发明不同爆破次数下地表沉降值；
图13为本发明不同爆破次数下导洞1拱顶沉降值；
图14为本发明不同爆破次数下导洞2拱顶沉降值；
图15为本发明不同爆破次数下大拱顶沉降值；
图16为本发明不同爆破次数下拱底隆起值；
图17本发明中衰减因子A＝0.01下高边墙水平位移随时间变化图；
图18本发明中衰减因子A＝0.1下高边墙水平位移随时间变化图；
图19本发明中衰减因子A＝0.2下高边墙水平位移随时间变化图；
图20本发明中衰减因子A＝0.5下高边墙水平位移随时间变化图。
附图标号说明：
1、地表沉降线；2、大拱顶监测点；3、导洞1拱顶监测点；4、导洞2拱顶监测点；5、爆破影响区；6、高边墙水平位移线；7、拱底监测点。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，本实施例的方法可以包括以下步骤：
S1、建立复杂扰动下岩体弹塑性耦合损伤模型；
具体而言，如图9所示，为扰动作用下的岩石耦合损伤模型图；
S2、引入渗透系数演化方程，建立复杂扰动下岩体弹塑性损伤-渗流耦合模型；
S3、将岩体弹塑性损伤模型参数进行服从Weibull分布的概率表达；
S4、建立考虑时间效应的岩体弹塑性损伤-渗流耦合模型；
S5、给定工程条件即施加边界条件并输入对应的材料参数，基于步骤S1-S4进行复杂扰动作用下考虑时间效应的岩体弹塑性损伤-渗流场耦合分析，获取当前工程条件对应的安全性评价数据。
可选的，在其中一个实施例中，所述S1中各向异性岩体弹塑性损伤模型对应的模型公式为：
f＝σ1-σ3-σci[(mbσ3/σci+s)]a (1)
其中，σ1和σ3分别为岩体的最大主应力和最小主应力；σci为完整岩石的单轴抗压强度；mb、a为针对不同岩体的量纲为一的强度参数；s为反映岩体破碎程度的强度参数。
使用过程中参数mb，s和a的取值方法为：
其中，GSI为地质强度指标(Geological strength index)；mi反映了岩石的软硬程度；D为扰动系数，反映了施工作业对岩体的扰动程度，是一个损伤结果量。本发明中将D视为冻融循环或爆破扰动下岩体的耦合损伤Dc，其表达式为：
其中，α、β为弹塑性损伤演化方程中的常参数；N为冻融循环或循环爆破扰动次数；Nt为当前的扰动次数；Dt为当前扰动次数下的岩体损伤值；U、W为与岩体材料性质相关的常参数；为等效塑性应变，其具体表达式为：
为等效塑性应变，其表达式为：
其中，εp1、εp2、εp3分别为三个方向上的主塑性应变。
引入岩体渗透系数演化方程，建立复杂扰动下岩体弹塑性损伤-渗流耦合模型。当岩体处于弹性状态时，其渗透率总体变化不大，根据Kozeny-Carman渗透系数演化方程可以表示为：
其中，k0为初始渗透系数；εv为体积应变；φ0为初始孔隙率。
当岩石处于塑性状态时，其内部会发生损伤，渗透性也会发生突变，此时，渗透系数的演化方程可以表示为：
其中，ξ为跳突系数，用于描述泥岩破裂时渗透性突增；A'＝1/(e-1/θ-1)，B'＝-1/(e-1/θ-1)；θ为经验参数。
建立岩体孔隙度演化方程如下：
其中，b和R为经验参数，用于描述损伤变量对孔隙度的影响，且R>1。
将岩体弹塑性损伤模型参数进行服从Weibull分布的概率表达。假设第i个RVE，模拟代表性体积元，中的某个力学参数Pi可以用Weibull统计分布函数进行描述，即
Pi＝P0·xi (10)
其中，P0为岩石式样的宏观性质；xi为相应的Weibull随机数。
计算时应根据Monte-Carlo随机抽样方法，对各个单元进行初始化赋值。在前述模型基础上，假设所研究的岩石试样中，其力学参数mb符合Weibull分布规律，其余参数均为统一值，即：
其中，mb0为强度参数mb的平均值。
首先在(0,1)范围内生成n个(n为总的单元个数)服从均匀分布的随机数u，再利用逆变法将这组随机数映射为服从Weibull分布的随机数，映射函数为：
mb＝mb0(-ln(u))1/ζ (12)
对各划分单元材料进行随机属性赋值，首先根据上述方法通过编写程序生成n个服从Weibull分布的随机数，然后通过PYTHON语言中的循环语句生成n个材料数组，最后对n个单元进行遍历赋值。
设其某一力学参数P随时间的变化规律可以用公式(17)进行描述，具体如下：
Pt＝P∞+(P0-P∞)e-At (13)
其中，P0为力学参数的初始值；P∞为强度参数的长期强度值，当t趋近于无限大时，P∞＝P0；A为强度衰减因子，表征了强度衰减的速率。
进一步的，给定工程条件即施加边界条件并输入对应的材料参数，对复杂扰动作用下考虑时间效应的岩体弹塑性损伤-渗流场耦合计算，获得工程稳定性评价参数，具体过程包括：
(1)根据根据耦合损伤Dc和流变时间t，对初始强度参数进行修正；
(2)对所述岩体弹塑性损伤模型进行弹性预测，求解对应的预测应力；
求解对应的预测应力的公式为：
其中，σn为tn时刻的应力；为tn时刻的内变量；Δεn+1为tn+1时刻的应变增量。
(3)将预测应力代入屈服函数，判断屈服函数是否大于零。若大于零，利用公式(19)对预测应力进行塑性修正：
其中，Δλ为塑性因子增量；h为塑性势函数对应力的导数，g为塑性势函数；Δσp为塑性应力增量，Δσp＝Δλh。
(4)根据公式(5)计算出新的损伤变量，并利用公式(20)对应力再次进行修正：
其中，σ'n+1为最终计算所得的应力。
(5)根据应力场计算所得的损伤值和塑性应变值，以及初始渗透系数和孔隙率，根据式(7)、(8)和(9)计算当前荷载步下的渗透系数和孔隙度，并在渗流场中对孔隙水压p进行求解，并利用有效应力公式(17)对应力再次进行修正：
其中，为有效应力；α为Biot系数；δij为克罗内克尔符号。
(6)重复步骤S51和S55，直到两时步之间的计算差值小于预设收敛值为止。
本发明中的应力-损伤-渗流完全耦合计算过程通过ABAQUS子程序接口UMAT和USDFLD的联合使用进行实现。计算过程中，首先利用UMAT中独立编写的复杂扰动下考虑时间效应的岩体弹塑性损伤模型及ABAQUS中自带的Soil单元进行耦合计算，同时利用PYTHON语言实现模型力学参数在单元上的随机赋值。在同一时步内保持积分单元的渗透系数和孔隙率不变；然后通过USDFLD子程序接口提取当前时步计算所得的体积应变εv和损伤值Dc，利用式公式(7)、公式(8)和公式(9)中，对渗透系数和孔隙度进行计算，并将计算结果以场变量Field Variables的形式进行输入，作为下一时步的计算参数，耦合计算流程图如图2所示。
为验证应力-渗流耦合场的合理性，建立单轴压缩数值计算模型如图2所示。模型高0.01m，宽0.005m，顶部为位移加载边界条件，底部为固定边界条件。模型共分为5151个节点和5000个单元。计算模型的计算参数如表1所示，初始孔隙率为0.3。本节只对单元破坏过程中渗透系数的变化规律及程序的准确性进行验证，因此整个计算过程只在静力学条件下进行。
表1计算参数
设置突跳系数ξ＝2，θ＝0.3，提取压缩过程中损伤结点的应力-应变曲线及渗透系数扩大倍数，绘值曲线如图4所示。
由图4可以看出，当节点处于弹性状态时，渗透系数演化规律由式公式(7)控制，此时其随加载时步的变化并不明显；当岩石内部出现塑性损伤时，渗透系数由公式(8)控制，随着损伤值的累积发生一定突跳，这与相关文献的试验结果相对应。计算结果验证了本文所建弹塑性-损伤-渗流耦合模型及其求解程序的正确性。
图5和图6分析了不同参数对计算结果的影响。突跳系数ξ和经验参数θ越大，渗透系数随损伤的变化速率也越大。
图7和图8为不同参数b和R对孔隙度的影响规律。由图可以看出，b值越大，随着损伤地累积，孔隙度增长速率越快；R值越大，孔隙度变化速率越慢。
某地铁站主要采用钻爆法进行开挖，频繁的爆破作业会降低围岩质量，因此有必要研究爆破损伤作用下车站的开挖稳定性。由于缺乏车站现场的声波监测数据，计算过程中以声波监测数据为基础，设置爆破区域如图10所示，包括：地表沉降线1、大拱顶监测点2、导洞1拱顶监测点3、导洞2拱顶监测点4、爆破影响区5、高边墙水平位移线6、拱底监测点7。对爆破影响区域分别设置8次爆破作用下的爆破损伤值，数值分别为0.032，0.088，0.13，0.15，0.17，0.18，0.19和0.21。
提取不同爆破次数下，车站关键测点的位移值如图11至图16所示。由计算结果可以看出，爆破作用降低了围岩的力学参数，在开挖卸荷后，各测点的位移均随着爆破损伤的累积逐渐增大。高边墙水平位移最大值在初次爆破作用下为7.73mm，第八次爆破后增加为9.55mm。地表沉降最大值由初次爆破作用下的-24.16mm增加为-29.49mm。不同开挖步序下，导洞拱顶监测点的位移值也均有所增大。考虑一次爆破的情况下，不同步序下导洞1拱顶沉降值分别为-0.78mm、-0.92mm、-2.18mm、-3.6mm、-3.87mm、-4.16mm、-4.39mm和-5.28mm。在第八次爆破冲击作用后，不同步序下导洞1拱顶沉降值分别为-0.95mm、-1.16mm、-2.67mm、-4.34mm、-4.63mm、-5.03mm、-5.31mm和-6.44mm。导洞2、大拱顶以及拱底的位移值变化规律均与导洞1拱顶沉降规律相同。1次爆破冲击作用下，车站开挖完成后导洞2拱顶沉降值为-5.33mm，经历8次爆破冲击作用后，该值变为-6.52mm。1次爆破冲击作用下，大拱顶最大沉降值为-7.06mm，经历8次爆破冲击作用后，其沉降值变为-8.61mm。而高边墙拱底最大隆起值在1次爆破冲击作用下为3.61mm，在第8次爆破冲击作用增长为4.37mm。计算结果表明，在爆破荷载的冲击作用下，围岩的力学性质发生弱化，导致在车站开挖卸荷后，各关键测点的位移均有所增大。因此，在爆破次数频繁的断面，设计时应该充分考虑爆破荷载作用对围岩的弱化作用，以免造成工程事故。而本文所建立的模型和方法能够很好的反映爆破扰动、弹塑性损伤以及渗流作用对围岩稳定性的影响，但在使用时需要通过反分析或室内试验的手段对参数进行准确标定，从而更为合理的对工程稳定性进行评价。
图17至图20为不同衰减因子下高边墙水平位移变化曲线。由计算结果可以看出，高边墙开挖完毕初期，位移值随时间有明显地增加，并逐渐趋于稳定。衰减因子决定了位移值的变化速率，衰减因子越大，位移的增长速率越快。由于长期强度不变，最终位移的稳定值也基本相同。根据前文分析，位移变化速率受支护措施和地质环境的影响，过快的变化速率会导致结构失稳，发生工程事故。因此在高边墙开挖完毕初期应密切关注监测数据，并适时增强支护，待监测数据稳定后方可进行下一步的施工。
整体有益效果：
本发明从02版的Hoek-Brown准则中的扰动系数出发，将循环爆破和冻融循环造成的扰动视为疲劳损伤，建立了同时考虑扰动因素与应力重分布双重作用下的岩体破坏准则；给出了体积应变-损伤-渗流耦合演化方程；基于有限元方法给出了复杂扰动下岩体应力-损伤-渗流耦合模型的数值求解算法；将复杂扰动下岩体应力-损伤-渗流耦合模型应用于实际工程的稳定性评价当中，为类似工程的安全施工提供了一定的设计依据。
最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。
Claims (6)
1.一种复杂扰动条件下岩体弹塑性损伤-渗流耦合计算方法，其特征在于，包括：
S1、建立复杂扰动条件下岩体弹塑性耦合损伤模型，所述复杂扰动条件包括：冻融循环扰动、循环爆破扰动；
S2、引入渗透系数演化方程，建立复杂扰动条件下岩体弹塑性损伤-渗流耦合模型；
S3、将岩体弹塑性损伤模型参数，转换为服从Weibull分布的概率表达形式；
S4、建立考虑时间效应的岩体弹塑性损伤-渗流耦合模型；
S5、给定工程条件，所述工程条件即施加边界条件并输入对应的材料参数，基于步骤S1-S4进行复杂扰动条件作用下考虑时间效应的岩体弹塑性损伤-渗流场耦合分析，获取当前工程条件对应的安全性评价数据。
2.根据权利要求1所述的一种复杂扰动条件下的岩体弹塑性损伤-渗流耦合计算方法，其特征在于，S1中所述复杂扰动下岩体弹塑性耦合损伤模型为：
f＝σ1-σ3-σci[(mbσ3/σci+s)]a (1)
其中，σ1和σ3分别为岩体的最大主应力和最小主应力；σci为完整岩石的单轴抗压强度；mb、a为针对不同岩体的量纲为一的强度参数；s为反映岩体破碎程度的强度参数；三个强度参数的计算公式为：
其中，GSI为地质强度指标；mi为力学参数，用于反映岩石的软硬程度；D为扰动系数，反映了施工作业对岩体的扰动程度，是一个损伤结果量；
复杂扰动下岩体的耦合损伤Dc的计算公式为：
其中，Dc为复杂扰动下岩体的耦合损伤，α、β为弹塑性损伤演化方程中的常参数；N为冻融循环或循环爆破扰动次数；Nt为当前的扰动次数；Dt为当前扰动次数下的岩体损伤值；U、W为与岩体材料性质相关的常参数；为等效塑性应变；
等效塑性应变的计算公式为：
其中，εp1、εp2、εp3为笛卡尔坐标系中x、y、z三个方向上的主塑性应变。
3.根据权利要求1所述的一种复杂扰动条件下的岩体弹塑性损伤-渗流耦合计算方法，其特征在于，所述S2引入岩体渗透系数演化方程，建立复杂扰动条件下的岩体弹塑性损伤-渗流耦合模型，包括：
S21、建立岩体渗透系数演化方程；
当岩体处于弹性状态时，渗透系数的演化方程为：
其中，k0为初始渗透系数；εv为体积应变；φ0为初始孔隙率；
当岩石处于塑性状态时，渗透系数的演化方程为：
其中，ξ为跳突系数，用于描述泥岩破裂时渗透性突增；A'和B'分别为常参数，其表达式为A'＝1/(e-1/θ-1)，B'＝-1/(e-1/θ-1)；θ为经验参数；
S22、建立复杂扰动条件下的岩体弹塑性损伤-渗流耦合模型：
其中，b和R为经验参数，用于描述损伤变量对孔隙度的影响，且R>1。
4.根据权利要求1所述的一种复杂扰动条件下的岩体弹塑性损伤-渗流耦合计算方法，其特征在于，所述S3将岩体弹塑性损伤模型进行服从Weibull分布的概率化表达：
S31、RVE为模拟代表性体积元，第i个RVE单元的力学参数Pi的计算公式为：
Pi＝P0·xi (10)
其中，P0为岩石式样的宏观性质；xi为相应的Weibull随机数；
双参数Weibull分布下，x的累积分布函数W(x)为：
x的概率密度函数w(x)为：
x的期望值为：
x的方差为：
其中，x是大于零的服从Weibull分布的随机数，用于表征岩体的非均质性，x0为随机数的平均值；xu为x的阈值；ζ为形状参数，即均质性指数；x0为尺寸参数，与数学期望有关；Г为Gamma函数；E(x)为数学期望函数；Var(x)为方差函数；
S32、根据Monte-Carlo随机抽样方法，对各个RVE单元进行初始化赋值；岩石的力学参数mb符合Weibull分布规律，其余参数均为统一值，建立累积分布函数：
在(0,1)范围内生成n个服从均匀分布的随机数u，所述n为RVE单元的单元数，用逆变法将这组随机数生成为服从Weibull分布的随机数，映射函数为：
mb＝mb0(-ln(u))1/ζ (16)
其中，mb0为强度参数mb的平均值，u为在(0,1)范围内生成n个服从均匀分布的随机数，n为RVE单元的单元数；
将生成的n个服从Weibull分布的随机数，通过PYTHON语言中的循环语句生成n个材料数组，对n个RVE单元进行遍历赋值。
5.根据权利要求1所述的一种复杂扰动条件下的岩体弹塑性损伤-渗流耦合计算方法，其特征在于，S4中所述建立考虑时间效应的岩体弹塑性损伤-渗流耦合模型：
对于RVE单元而言，力学参数P随时间的变化规律为：
Pt＝P∞+(P0-P∞)e-At (17)
其中，P为P0为力学参数的初始值；P∞为强度参数的长期强度值，t为流变时间，当t趋近于无限大时，P∞＝P0；A为强度衰减因子，表示强度衰减的速率。
6.根据权利要求1所述的一种复杂扰动条件下的岩体弹塑性损伤-渗流耦合计算方法，其特征在于，S5中给定工程条件，所述工程条件即施加边界条件并输入对应的材料参数，基于步骤S1-S4进行复杂扰动作用下考虑时间效应的岩体弹塑性损伤-渗流场耦合分析，包括：
S51、根据耦合损伤Dc和流变时间t，对初始强度参数进行修正；
S52、对所述岩体弹塑性损伤模型进行弹性预测，求解对应的预测应力；
求解对应的预测应力的公式为：
其中，tn为计算时步，σn为tn时刻的应力；为tn时刻的内变量；Δεn+1为tn+1时刻的应变增量；
S53、将预测应力代入S1中所述复杂扰动下岩体弹塑性耦合损伤模型，判断所述损伤模型是否大于零；若大于零，对预测应力进行塑性修正，对预测应力进行塑性修正的计算公式为：
其中，Δλ为塑性因子增量；h为塑性势函数对应力的导数，g为塑性势函数；Δσp为塑性应力增量，Δσp＝Δλh；
S54、根据损伤变量，对应力再次进行修正：
其中，σ'n+1为计算所得的最终应力；
S55、根据应力场计算所得的损伤值和塑性应变值、初始渗透系数和孔隙率、当前荷载步下的渗透系数和孔隙率、以及通过渗流有限元求解器对孔隙水压p进行求解，对最终应力再次进行修正，对最终应力进行修正的公式为：
其中，为有效应力；α为Biot系数；δij为克罗内克尔符号；
S56、重复步骤S51至S55，直到两时步之间的计算差值小于预设收敛值为止，所述收敛值的取值范围为1e10-5～1e10-3。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110891535.2A CN113505514B (zh) | 2021-08-04 | 2021-08-04 | 一种复杂扰动条件下岩体弹塑性损伤-渗流耦合计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110891535.2A CN113505514B (zh) | 2021-08-04 | 2021-08-04 | 一种复杂扰动条件下岩体弹塑性损伤-渗流耦合计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113505514A CN113505514A (zh) | 2021-10-15 |
| CN113505514B true CN113505514B (zh) | 2024-01-05 |
Family
ID=78015012
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110891535.2A Active CN113505514B (zh) | 2021-08-04 | 2021-08-04 | 一种复杂扰动条件下岩体弹塑性损伤-渗流耦合计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113505514B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115270523B (zh) * | 2022-09-27 | 2022-12-20 | 成都理工大学 | 一种岩体损伤影响深度预测方法 |
| CN115600398B (zh) * | 2022-10-10 | 2023-06-16 | 长安大学 | 一种基于蒙特卡洛模拟的大型洞室岩体参数概率估计方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105401939A (zh) * | 2015-11-30 | 2016-03-16 | 中国石油大学(北京) | 一种多因素耦合作用下的煤层井壁稳定性分析方法 |
| CN110261573A (zh) * | 2019-05-16 | 2019-09-20 | 同济大学 | 一种高位岩质滑坡稳定性动态数值评价方法 |
| CN111444641A (zh) * | 2020-02-10 | 2020-07-24 | 大连海事大学 | 一种考虑冻融环境下的岩体工程稳定性分析方法 |
| CN111695285A (zh) * | 2020-06-17 | 2020-09-22 | 大连海事大学 | 一种各向异性岩体应力-损伤-渗流耦合数值模拟方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110222439B (zh) * | 2019-06-12 | 2020-01-07 | 四川大学 | 基于Abaqus平台疲劳损伤与寿命评估方法 |
-
2021
- 2021-08-04 CN CN202110891535.2A patent/CN113505514B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105401939A (zh) * | 2015-11-30 | 2016-03-16 | 中国石油大学(北京) | 一种多因素耦合作用下的煤层井壁稳定性分析方法 |
| CN110261573A (zh) * | 2019-05-16 | 2019-09-20 | 同济大学 | 一种高位岩质滑坡稳定性动态数值评价方法 |
| CN111444641A (zh) * | 2020-02-10 | 2020-07-24 | 大连海事大学 | 一种考虑冻融环境下的岩体工程稳定性分析方法 |
| CN111695285A (zh) * | 2020-06-17 | 2020-09-22 | 大连海事大学 | 一种各向异性岩体应力-损伤-渗流耦合数值模拟方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113505514A (zh) | 2021-10-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113505514B (zh) | 一种复杂扰动条件下岩体弹塑性损伤-渗流耦合计算方法 | |
| Ghorbani et al. | A comprehensive solution for the calculation of ground reaction curve in the crown and sidewalls of circular tunnels in the elastic-plastic-EDZ rock mass considering strain softening | |
| CN111705808A (zh) | 一种适用于悬挂式深基坑工程的土体参数动态反演分析方法 | |
| CN111444641A (zh) | 一种考虑冻融环境下的岩体工程稳定性分析方法 | |
| Kostina et al. | An Applicability of Vyalov’s equations to ice wall strength estimation | |
| Parsa-Pajouh et al. | Analyzing consolidation data to predict smear zone characteristics induced by vertical drain installation for soft soil improvement | |
| CN110348098B (zh) | 一种黏质黄土隧道开挖模拟参数赋值方法 | |
| CN112329287A (zh) | 一种基于试桩监测数据的p-y曲线贝叶斯学习方法 | |
| Consoli et al. | Numerical analysis of pressuremeter tests and its application to the design of shallow foundations | |
| Guilbaud | Damage plastic model for concrete failure under impulsive loadings | |
| Adel et al. | Evaluation of the bearing capacity of a single pile by numerical analysis and various methods | |
| CN115221727B (zh) | 一种基于含水率的岩体的数值仿真模型参数确定方法 | |
| Ju et al. | Implicit integration of an anisotropic egg-shaped elastoplastic model for saturated soft clay | |
| Gusev et al. | Research on regularities of deformation behavior of building structures in the areas of technogenic impact caused by mining | |
| Dyakonov et al. | Reverse analysis of geotechnical monitoring results for the estimation of the diaphragm walls stress-strain | |
| CN115186513B (zh) | 一种块石填料长期剪应变的预测方法 | |
| CN109783947B (zh) | 一种采水型地裂缝数值模拟方法 | |
| Bakeer et al. | Determination the capacity reduction factor of masonry walls under buckling-a numerical procedure based on the transfer-matrix method | |
| CN116484568A (zh) | 一种考虑疲劳和时间效应的岩体稳定性计算方法 | |
| Tao et al. | A comparison between EnKF and MCMC-based Bayesian updating for consolidation settlement prediction | |
| CN115438729A (zh) | 一种基于单元线性回归的隧道隆起预测方法 | |
| Raines et al. | Numerical modelling of membrane penetration effects on undrained triaxial tests | |
| CN116561983A (zh) | 建立多轴约束下uhpc修正连续损伤塑性本构模型的方法及系统 | |
| Sharafutdinov et al. | A study of the ground volume loss modeling technique influence the soil displacement in course of shield tunneling | |
| Liew et al. | Bending moment interpretation of structural element with measured deflection profile |
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 |