Establishment and application of two-dimensional axisymmetric temperature field model for high-level radioactive waste disposal repository
-
摘要: 高放废物处置库温度场长期演化特征是处置库缓冲层安全性能评估和设计、处置容器空间优化布局的重要依据之一。建立了处置单元二维轴对称双层热分析模型,以揭示处置库多重屏障系统近场温度演化规律。对传热控制方程应用拉普拉斯和有限傅里叶正弦变换,得到缓冲层和围岩层温度的拉普拉斯域解,借助Crump方法对拉普拉斯域解进行数值反演获得对应的时域解;通过与线热源解和数值解的对比验证了模型的正确性,分析了处置单元近场多屏障系统温度场的时空分布特征,探讨了不同参数对缓冲层峰值温度的影响;借助模型半解析解确定最小处置容器间距以及对原位试验结果进行预测。结果表明:①与线热源解相比,模型半解析解可以更准确地模拟缓冲层温度变化;②当处置隧道间距超过40 m时,其对缓冲层峰值温度影响较小;③在处置隧道间距取40 m时,最小处置容器间距为7.7 m;④模型半解析解可以较好地预测原位试验结果。Abstract: The long-term evolution of temperature field in high-level radioactive waste disposal repository is one of the important bases for the safety evaluation and design of buffer layer and the spatial optimization of disposal container. A two-dimensional axisymmetric thermal analysis model with two layers of disposal unit is established to reveal the near-field temperature evolution of multi-barrier system. The Laplace and finite Fourier sine transforms are applied to the governing equations for heat transfer to obtain the temperature solutions to the buffer layer and the surrounding rock layer in the Laplace domain. The corresponding solutions in time domain are obtained by numerical inversion of the Laplace domain solutions with the help of the Crump method. The correctness of the model is verified by comparing with the linear heat source solution and the numerical solution. The temporal and spatial distribution characteristics of the temperature field in the near-field multi-barrier system of the disposal unit are analyzed, and the effects of different parameters on the peak temperature of the buffer layer are discussed. The semi-analytical solution of the model is used to determine the minimum disposal container spacing and predict the results of the in-situ tests. The results show that the semi-analytical solution of the model can simulate the temperature change of buffer layer more accurately by comparing with the linear heat source solution. When the tunnel spacing exceeds 40 m, it has small effects on the peak temperature of buffer layer. When the tunnel spacing is 40 m, the container spacing is 7.7 m. The semi-analytical solution of the model can well predict the results of the in-situ test.
-
0. 引言
土体参数空间变异性对土坡稳定性的影响已受到岩土工程界的广泛关注,繁冗的计算量是空间变异土坡可靠度分析面临的瓶颈问题。为解决该问题,诸多学者在蒙特卡洛模拟(MCS)分析中采用高效的代理模型取代耗时的边坡稳定分析模型,从而减少计算量、提高分析效率。在处理空间变异土坡可靠度问题时,常用的随机场模拟方法需要将土体离散为成百上千个随机变量,带来的“维度灾难”问题,导致许多代理模型难以适用。尽管目前存在一些降维方法,如采用Karhunen-Loève(KL)级数展开模拟土坡随机场,并在随机变量数量较少的低维标准正态空间构建代理模型。然而,当土性参数自相关距离较小、场地区域较大或采用指数型自相关函数时,为保证随机场模拟具有足够精度,KL级数展开所需的随机变量数量仍然较多[1-4]。开发不受随机变量维度影响,并在高维物理参数空间(X空间)中适用的代理模型,进而进行空间变异土坡可靠度分析的研究依然有限[5]。
边坡代理模型的构建需要足够的训练样本点,一般而言,训练样本点越多,模型越精确、预测误差越小。然而,过多的训练样本点意味着需要调用更多次数的边坡稳定分析模型,从而大幅度增加计算量,影响计算效率。如何采用较少的训练样本点获得精度较高的边坡代理模型仍是一个关键技术难题。当前基于代理模型的可靠度方法,其训练样本点的选取主要分为全局取样法和局部取样法两大类。全局取样法试图构建一个在不确定参数的全域内均具有较高预测精度的代理模型,因此往往需要较多的训练样本点[1-3, 5]。需要指出的是,基于MCS法的边坡可靠度分析属于分类问题,即失效概率的求解只需准确判断安全系数Fs是否小于1.0,而无需得出其精确值。因此,全局取样法会产生大量对求解边坡失效概率贡献不大的训练样本点,从而降低计算效率。局部取样法旨在参数空间的关键区域内(如极限状态面LSS附近)选取训练样本点,进而构建边坡的局部代理模型,提高其在LSS附近区域的分类预测精度[6]。同时,局部取样法能够较好地与主动学习策略相结合,通过借助不同的主动学习函数,逐个识别LSS附近的最优训练样本点,从而减少训练样本点的数量[6-8]。然而,现有的局部取样法大多只适用于随机变量较少(低维)的情况,在考虑参数空间变异性的随机场工况时仍然面临着高维变量拟合困难的问题。
强度折减法是求解边坡安全系数最普遍的方法之一,该方法将边坡安全系数定义为将土体抗剪强度折减以使得边坡接近失稳状态的折减系数[9]。换而言之,若将土体抗剪强度除以Fs,则折减后的土体强度参数可恰好使边坡处于极限平衡状态。基于强度折减理论,可以直接识别出LSS上的临界样本点,进而有助于局部代理模型的构建。目前,关于这方面的尝试仍十分有限。
针对上述问题,从强度折减理论和局部取样法角度出发,提出一种基于强度折减采样(strength reduction sampling, SRS)与高斯过程回归(Gaussian process regression, GPR)的空间变异土坡自适应可靠度分析方法(SRS-GPR)。该方法基于SRS生成位于LSS上的临界训练样本点,通过GPR模型构建土体高维随机变量与边坡安全系数之间的函数关系。在此基础上,提出一种主动学习策略来有效识别靠近LSS区域附近的最优训练样本点,通过迭代更新逐步提高GPR模型在失效区域附近的预测精度,从而降低所需训练样本点数量,提高计算效率。最后,通过两个空间变异土坡算例验证所提方法的有效性。
1. 土体参数空间变异性模拟
土体性质的空间变异性表现为土体参数在空间不同位置间存在相关性。为合理描述土体性质的空间变异性,随机场模型得到众多学者的青睐。在随机场模型中,由于现场勘察数据有限,通常采用理论自相关函数来模拟土体参数间的自相关性,采用的二维高斯型自相关函数表达式为[10]
ρ(Δh,Δ)=exp[−(Δ2hlh2+Δ2vl2v)]。 (1) 式中:Δh,Δv分别为水平和垂直方向上任意两点间的相对距离;lh,lv分别表示水平和垂直方向的自相关距离。此外,高斯型自相关距离l为波动范围δ的1/√π倍。
采用计算精度和效率较高的KL级数展开法离散土体强度参数非高斯随机场。基于KL展开的独立高斯随机场H(x, y)表达式为
H(x,y)≈ˆH(x,y)=μ′+σ′ns∑i=1√λifi(x,y)ξi。 (2) 式中:(x, y)为二维计算区域Ω = {(x, y)| x∈[0, L1], y∈[0, L2]}内任意点的坐标,其中L1,L2分别为随机场计算区域的水平和垂直尺寸;μ′,σ′分别为随机场的均值和标准差;λi,fi(·)分别为自相关函数的第i个特征值和特征函数;ξi为第i组独立标准正态随机变量;Ĥ(x, y)为考虑截断误差的近似随机场;ns为KL展开截断项数,其数值大小决定随机场模型精度,一般采用期望能比率因子ε来确定:
ε=∑nsi=1λi∑∞i=1λi=∑nsi=1λiL1L2。 (3) 根据文献[3]建议,可采用ε ≥ 95%作为确定截断项数ns的依据,以确保离散后的随机场具有较高精度。当土体参数服从对数正态分布时,便可将式(2)的独立高斯随机场通过等概率转换方法得到独立非高斯随机场如下:
ˆH(x,y)=exp[μ′lnx+σ′lnxns∑i=1√λifi(x,y)ξi]。 (4) 式中:μlnx′,σlnx′分别为参数对数化后的均值和标准差。最后,考虑参数之间的互相关性(如黏聚力和内摩擦角之间的负相关性),借助乔列斯基法分解参数间的互相关系数矩阵,便可得到土体多参数的相关非高斯随机场。
2. 空间变异土坡自适应可靠度分析的SRS-GPR方法
2.1 失效概率评估
在边坡可靠度分析中,可将代理模型与MCS方法结合进行失效概率Pf计算[11]:
Pf≈1NMCSNMCS ∑i=1I[Fs(xi)<1]≈1NMCSNMCS∑i=1I[R(xi)<1]。 (5) 式中:xi为土体随机场参数的输入向量;Fs(xi)为利用边坡稳定分析模型计算的安全系数;R(xi)为代理模型预测的边坡安全系数;I[·]为指示性函数,当Fs(xi) < 1(即边坡失效),I[·] = 1,否则I[·] = 0;NMCS为模拟次数。当代理模型构建完成后,Pf的计算仅需基于该模型求解边坡安全系数,而无需重复调用边坡稳定分析模型,因此可以显著提高边坡可靠度分析效率[12]。
2.2 高斯过程回归
由于高斯过程回归在处理复杂非线性问题时具有较强的拟合优势[13-16],采用GPR构建边坡Fs的代理模型,以下简要介绍GPR原理。高斯过程是随机变量的集合,且假定任意有限个随机变量均服从联合高斯分布[17]。高斯过程f(x)可通过均值函数和协方差函数来表征:
m(x)=E(f(x)), (6) k(x,x′)=E[(f(x)−m(x))(f(x′)−m(x′))]。 (7) 式中:m(x)为均值函数,通常取零;E(·)表示期望算子;k(x, x′)为输入样本x和x′的协方差函数,这里采用平方指数型[16]。当给定一组训练样本Xt = [x1, x2, …, xn]T时,其相应的响应向量为Yt = [y1, y2, …, yn]T,则可构建输入变量与响应量之间的关系如下:
vi=f(xi)+ε。 (8) 式中:ε为噪声值,服从均值为零、方差为σ2的独立高斯分布。以f值为条件的Yt的分布可由各向同性的高斯分布表示:
p(Yt∣f,Xt)=N(f,σ2I)。 (9) 式中:I为单位矩阵。由于GPR是一种贝叶斯方法,它假设回归函数的均值先验为零,因而Yt的先验分布可表示为
p(Yt∣Xt)=∫p(Yt∣f,Xt)p(f∣Xt)df=N(0,K(Xt,Xt)+σ2I)。 (10) 式中:K(·)表示分量为Kij = k(xi, xj)的协方差矩阵。因此,在给定训练输入样本的情况下,可以通过最大化输出的对数似然来求解式(10)中的超参数θ = [λ, σf, σ]。在确定最优超参数θ后,先验条件下新样本X*处的响应值Yt与预测值Y*的联合分布可表示为
[YtY∗]∼N(0,[K(Xt,Xt)+σ2IK(Xt,X∗)K(X∗,Xt)K(X∗,X∗)])。 (11) 进而可得Y*的后验分布:
p(Y∗∣X∗,Xt,Yt)=N(μ∗,σ∗2), (12) μ∗=K(X∗,Xt)[K(X∗,Xt)+σ2I]−1Yt, (13) σ∗2=K(X∗,X∗)−K(X∗,Xt)[K(Xt,Xt)+ σ2I]−1K(Xt,X∗)。 (14) 式中:μ*,σ*2分别为预测值Y*在新样本X*处的均值和方差向量。
2.3 强度折减采样法
强度折减法的原理是将边坡土体抗剪强度参数c和ϕ同时除以折减系数Fr,得到一组新的cr和ϕr值,然后作为一组新的边坡输入参数,再进行试算。土体的抗剪强度折减公式为
τr=(c+σntanϕ)/Fr=cr+σntanϕr, (15) cr=c/Fr, (16) tanϕr=tanϕ/Fr。 (17) 式中:c,ϕ分别为土体的黏聚力和内摩擦角;σn,τr分别为作用于土体的正应力和剪应力;cr,ϕr分别为折减后的黏聚力和内摩擦角。在强度折减法中,折减系数Fr通过迭代计算不断调整,直至边坡恰好达到失稳破坏(极限平衡状态),则此时的折减系数等于边坡的安全系数,即Fr = Fs。
根据强度折减理论,本文提出强度折减采样法,该方法可直接基于已知的边坡安全系数,得到相应的临界抗剪强度参数为
ce=c/Fs, (18) tanϕe=tanϕ/Fs。 (19) 式中:ce,ϕe为使边坡处于极限平衡状态的临界抗剪强度参数,即Fs(ce, ϕe) = 1.0,因此在随机参数空间由临界抗剪强度参数构成的点恰好位于边坡功能函数的极限状态面上。结合式(18),(19)可得
cetanϕe=ctanϕ。 (20) 由式(20)可知,在c-tanϕ参数空间,坐标原点O、任意一组强度参数构成的随机点(cj, tanϕj)和其对应的临界点(cej, tanϕej) 三者共线,j = 1, 2, …, k,如图 1所示。
当考虑土体强度参数空间变异性时,令x0 = {c1, …, cM, ϕ1, …, ϕM}表示参数随机场离散后的c-ϕ向量,其中ci,ϕi分别为黏聚力和内摩擦角的随机场实现值,i = 1, 2, …, M。基于强度折减理论,可得临界样本点xe0 = {ce1, …, ceM, ϕe1, …, ϕeM},其中
cei=ci/Fs,ϕei=arctan(tanϕi/Fs)。} (21) 综上,当存在一组随机场实现值x0时,可基于式(21)获取其相应的临界样本点xe0,此时的xe0恰好位于边坡极限状态面上,即Fs(xe0) = 1.0。值得说明的是,基于SRS法生成的临界样本点xe0不仅可以极大地提升边坡代理模型在极限状态面附近的预测精度,而且xe0的获取无须额外调用边坡稳定分析模型,可大幅提高边坡可靠度分析效率。
2.4 主动学习策略
选取合适的训练样本点不仅可以提高代理模型预测精度,而且能够减少调用边坡稳定分析模型次数。本文提出一种主动学习策略,该策略可在蒙特卡洛样本池中选取对失效概率评估贡献最大的样本点。该样本点应当满足:①尽可能靠近预测的LSS;②远离已有的训练样本点,以避免样本点过度聚集,造成信息冗余。首先从蒙特卡洛(MC)样本池中提取靠近预测的LSS的候选样本点,该候选样本点满足:
|R(xi)−1|max[|R(xi)−1|]<10−3。 (22) 式中:R(xi)为给定输入xi时的GPR模型预测的安全系数值。然后比较候选样本点与已有训练样本点之间的欧氏距离,并根据如下公式选择最优样本点xa作为新的训练样本点:
dmin,i=min(||xi−xj||), (23) xa=arg[max(dmin,i)]。 (24) 式中:xi为由式(22)提取的候选样本点,i = 1, 2, …, Nc;xj为已有的训练样本点,j = 1, 2, …, n。
式(22)~(24)可保证识别的最优训练样本点分布在边坡功能函数的极限状态边界附近。
2.5 计算流程
本文融合SRS、GPR模型和主动学习策略,提出空间变异性边坡自适应可靠度分析方法(SRS-GPR)。其中,SRS用以生成临界样本点,GPR模型作为边坡安全系数的代理模型,而主动学习策略用以识别新的最优训练样本点来更新GPR模型。综上,所提SRS-GPR方法计算流程如图 2所示,其具体实施步骤如下:
(1)确定边坡随机土体参数的统计特征,如概率分布、均值、标准差、波动范围(自相关距离)、自相关函数及互相关系数等,建立边坡稳定性分析模型。
(2)根据随机参数的概率分布,在X空间中生成NMCS个随机场蒙特卡洛样本点,形成样本池。
(3)生成一组随机场初始样本点S0,并调用边坡稳定模型求解相应Fs。在主动学习策略的协助下,初始样本数对边坡可靠度分析影响较小,少量的初始样本点(本文取20)即可满足计算精度要求[7]。为保证样本点空间分布的均匀性,通常可采用拉丁超立方抽样(LHS)、均匀设计和中心复合设计等方法生成初始样本点[14, 18-19]。
(4)基于SRS法生成与初始样本点对应的临界样本点Se0。随后,将初始样本点S0、临界样本点Se0及相应的安全系数组成训练样本集。
(5)在X空间内,基于步骤(4)中的训练样本集构建边坡安全系数的GPR代理模型。
(6)基于步骤(5)得到的GPR模型和步骤(2)生成的MC样本点进行边坡可靠度分析,采用式(5)计算边坡失效概率Pf。
(7)判断迭代计算是否终止。若σ(Pf,iter(5))/μ(Pf,iter(5)) < 0.3%[6],计算转至步骤(9);否则,转至步骤(8)。其中,Pf,iter(5)为GPR模型在当前连续5次迭代计算的失效概率[20],μ(·),σ(·)分别为相应的均值和标准差。
(8)采用2.4节所提主动学习策略确定新的最优训练样本xa,调用边坡稳定模型计算点xa处的安全系数,并基于式(21)生成相应的临界样本点xea。将xa,xea及其安全系数补充至训练样本集中,随后转至步骤(5)。
(9)终止迭代计算,输出步骤(6)最终所得边坡失效概率Pf。
3. 算例分析
本节将通过两个土坡算例来验证所提SRS-GPR方法的准确性、高效性、鲁棒性和适用性。其中,边坡安全系数采用简化Bishop法(BSM)计算,将基于直接MCS法求解的Pf作为精确解来验证SRS-GPR方法的准确性。采用边坡稳定性模型的平均调用次数Ncall量化计算效率,并与文献中其它方法比较来验证SRS-GPR方法的高效性。进行20次独立的边坡可靠度分析,采用95%置信区间验证SRS-GPR方法的鲁棒性。同时,设计多种计算工况,进一步验证SRS-GPR方法的适用性。
3.1 不排水饱和黏性土坡
以不排水饱和黏性土坡为例进行分析,二维边坡几何尺寸如图 3所示,坡高5 m,坡角为26.6°。边坡黏土重度为20 kN/m3,视为常量。采用高斯型自相关函数模拟土体不排水抗剪强度cu的空间变异性,根据文献[3,21~23],假设cu服从均值为23 kPa,变异系数为0.3的对数正态分布随机场,水平和垂直自相关距离分别取20,2 m。基于图 3几何尺寸建立边坡稳定性分析模型,采用KL展开方法离散cu随机场,边坡区域共划分为910个边长为0.5 m的随机场网格单元。利用MATLAB软件编写BSM程序求解边坡Fs。基于cu参数均值计算的Fs为1.356,与文献中采用BSM方法求解的1.356[3, 21]和1.355[22-23]一致,其相应最危险滑动面位置如图 3所示。
图 4给出cu随机场的一次典型实现,颜色越深表示cu值越大,颜色越浅则表示cu值越小。从图 4中不难发现参数cu值随着空间位置的变化而变化,且相邻位置间的cu值也存在着一定的相关性。
根据2.5节的计算流程,首先生成20万组cu随机场的MC样本点,从中选取20个初始样本点S0,求解边坡Fs,其相应的临界样本点Se0可通过式(21)得到,并将初始样本点S0和临界样本点Se0及其安全系数组合构成训练样本集。随后,利用训练数据集构建边坡Fs的GPR模型,并基于20万组MC样本点求解边坡Pf。当迭代计算不满足终止条件时,基于式(22)~(24)识别最优样本点xa,将点xa和相应临界样本点xea及其安全系数补充至训练数据集中,接着采用新的训练样本集重新训练GPR模型,并求解对应的边坡Pf。当终止条件满足时,迭代计算结束,输出最终的边坡Pf。图 5给出一次随机分析下,采用SRS-GPR方法进行边坡可靠度计算的迭代过程。由图 5(a)可见,随着迭代计算的进行,本文提出方法计算的Pf逐渐收敛到MCS精确解,计算终止时得到的边坡Pf(0.0883)与精确解(0.0882)吻合较好,此时的边坡稳定模型调用次数为78,GPR模型的超参数取值为θ = [466.14, 0.52, 0.0017]。由图 5(b)可见,大部分训练样本点的Fs值都接近1.0(即靠近LSS),说明本文提出的主动学习策略可有效识别边坡极限状态面附近的样本点。此外,由于临界样本点可直接利用式(21)获得,无需进行额外的边坡确定性分析,因此在迭代计算过程中,边坡确定性分析次数仅为用于构建GPR模型的训练样本数量的1/2。
为说明所提方法建立局部代理模型的特点,图 6给出所提方法与MCS方法计算的边坡安全系数累积分布函数(CDF)曲线。由图 6可知,当Fs取值在1.0附近时,所提方法得到的CDF曲线与MCS方法吻合较好;而当Fs取其它值时(如1.5附近),所提方法的CDF曲线精度较差。这是由于所提方法针对的是影响边坡稳定的关键区域,建立的是局部代理模型,故而所选取的训练样本点大多集中在LSS附近(如图 5(b)所示),因此所提方法得到的CDF曲线在Fs = 1.0附近相较于其它位置的预测精度更高。考虑到本文旨在准确计算边坡失效概率,故而所提方法通过聚焦于边坡极限状态区域,能够有效选取对失效概率评估贡献较大的关键样本点。
为降低随机模拟过程中的不确定性,采用SRS-GPR方法对该算例进行20次重复独立可靠度分析,并对结果取平均。表 1给出所提方法与其它方法的计算结果比较。由表 1可知,相较于其它可靠度分析方法,SRS-GPR方法平均仅需边坡确定性分析73.3次,便可得出与MCS精确解精度吻合的结果,显著提高边坡可靠度分析效率,且20次独立分析所得Pf的95%置信区间较窄,表明所提方法具有较好的鲁棒性。同时,相较于CNN-LHS、EQP-MRSM和RSSs-RSM方法,不难发现本文方法通过有效结合强度折减采样、主动学习策略和GPR模型,在X空间处理空间变异性土坡可靠度问题上具有显著优势。为了验证所提方法分析低概率问题的有效性,根据文献[3]将cu的变异系数取为0.15,所提方法计算的边坡失效概率为5.16×10-4,边坡确定性分析次数为100,与500万次直接MCS所得结果(5.19×10-4)的一致性较好。为进一步验证所提方法在考虑不同自相关距离时的适用性,设置水平和垂直自相关距离取值范围分别为[3]lv = [0.5 m, 3.0 m]和lh = [10 m, 40 m]。图 7给出边坡Pf随自相关距离的变化关系曲线。由图 7可知,对于不同工况下自相关距离值,本文方法均可给出较为准确的边坡失效概率。
表 1 不排水黏性土坡可靠度分析结果Table 1. Reliability results of undrained cohesive slope方法类别 Ncall Pf Fs计算方法 MCS 100000 0.0882 BSM LHS[3] 1000 0.0830 BSM 2阶多重SRSM[3] 1000 0.0790 BSM 3阶多重SRSM[24] 1000 0.0810 BSM GPR-RSM[4] 150 0.0740 BSM CNN-LHS[5] 100 0.1070 FDM EQP-MRSM[22] 3021 0.0740 BSM RSSs-RSM[23] 1821 0.0755 BSM SRS-GPRa 73.3 0.0881 BSM SRS-GPRb [64.3,84.2] [0.0854,0.0908] BSM 注:上标a为20次独立分析的平均结果;上标b为95% 置信区间;FDM为有限差分法;SRSM为随机响应面法;GPR-RSM为基于高斯过程回归的响应面法;CNN-LHS为基于卷积神经网络的拉丁超立抽样法;EQP-MRSM为基于等效参数的多重响应面法;RSSs-RSM为基于代表性滑面的响应面法。 3.2 黏性-摩擦土坡
以文献[10,21]中的黏性-摩擦(c-ϕ)土坡为例,对SRS-GPR方法进行验证。该边坡模型尺寸如图 8所示,坡高为10 m,坡角为45°,土体重度为20 kN/m3。边坡可靠度分析的土体统计参数如表 2所示,采用相关对数正态随机场模拟土体c和ϕ的空间变异性。基于表 2中参数均值,采用BSM方法求解的Fs为1.204,与文献中计算结果1.204[21],1.205[1, 23],1.206[3],1.208[10]相一致,其相应的临界滑动面见图 8所示。随机场离散后的边坡区域共剖分1190个边长为0.5 m的四边形单元和坡面附近过渡区20个三角形单元,如图 8所示。图 9给出c-ϕ随机场的一次典型实现,可见参数c与ϕ之间存在明显的负相关性,与实际情况吻合。
表 2 黏性-摩擦土坡可靠度分析的土体参数Table 2. Parameters for reliability analysis of cohesive-frictional slope土体参数 均值 变异系数 分布类型 波动范围 互相关系数 c/kPa 10 0.3 对数
正态δv = 4 m
δh = 40 mρc,ϕ = −0.5 ϕ/(°) 30 0.2 对数
正态δv = 4 m
δh = 40 m基于表 2中的参数统计特征,在一次随机模拟下采用本文方法计算的边坡Pf为0.0263,与10万次MCS计算的结果(Pf = 0.0269)吻合较好,边坡确定性分析次数为86,GPR模型的超参数取值为θ = [1456.51,0.82,0.0011],对应的可靠度分析迭代计算过程如图 10所示,其中归一化Pf定义为所提方法计算的失效概率与MCS精确解的比值。为保持与文献中的研究工况一致,图 10中还给出另外两种不同波动范围时的计算结果。由图 10(a)可见,在迭代计算初期,由于训练样本点较少,本文方法在LSS附近的预测精度较差,因此归一化Pf与精确解(即1.0)的差别较大。随着迭代计算的持续进行,用于训练GPR模型的有效随机场样本点逐渐被准确识别(如图 10(b)中大多数训练样本点的Fs值在1.0附近),本文方法的预测精度不断提升并最终收敛至精确解。值得注意的是,尽管该边坡所含随机变量的数量高达2420(包括1210个c和1210个ϕ),但本文方法仍然具有较高的计算精度,呈现出良好的适用性。
与3.1节类似,为降低随机模拟过程中的不确定性,对该黏性-摩擦土坡进行20次独立可靠度分析,并对结果取平均。将文献中的不同方法,如LHS,MARS,SIR-MARS等结果也列于表 3进行比较。由表 3可知,相较于其它可靠度分析方法,提出的SRS-GPR方法仅需调用边坡稳定模型小于100次,便可获得与精确解吻合的结果,同时,其较窄的置信区间也表明本文方法鲁棒性较好。
表 3 不同工况下的黏性-摩擦土坡可靠度分析结果Table 3. Reliability results of cohesive-frictional slope in different cases计算工况 方法类别 Ncall Pf Fs计算方法 δv = 4 m
δh = 40 mMCS 100000 0.0269 BSM LHS[1] 10000 0.0272 BSM MARS[1] 280 0.0292 BSM SIR-MARS[25] 120 0.0261 FDM SRS-GPRa 87.9 0.0267 BSM SRS-GPRb [78.4, 97.4] [0.0252, 0.0282] BSM δv = 8 m
δh = 40 mMCS 100, 000 0.0582 BSM LHS[1] 10000 0.0570 BSM MARS[1] 280 0.0546 BSM SIR-MARS[25] 120 0.0538 FDM SRS-GPRa 67.8 0.0593 BSM SRS-GPRb [60.4, 75.2] [0.0569, 0.0617] BSM δv = 4 m
δh = 80 mMCS 100, 000 0.0277 BSM LHS[1] 10000 0.0286 BSM MARS[1] 280 0.0305 BSM SIR-MARS[25] 120 0.0295 FDM SRS-GPRa 88.4 0.0283 BSM SRS-GPRb [81.8, 94.9] [0.0271, 0.0294] BSM 注:上标a为20次独立分析的平均结果;上标b为95%置信区间;FDM为有限差分法;MARS为多元自适应回归样条;SIR-MARS为切片逆回归与MARS组合方法。 为进一步系统研究本文提出方法在不同工况下的适用性,根据文献[10],选取参数ρc, ϕ,COVc,COVϕ,δv的取值范围为ρc, ϕ∈ [−0.7, 0.5],COVc ∈ [0.1, 0.7],COVϕ ∈ [0.05, 0.20],δv ∈ [1 m, 6 m]。当某一参数变化时,其余参数取定值:ρc, ϕ = 0,COVc = 0.3,COVϕ = 0.2,δv = 4 m,δh = 40 m。图 11~13分别给出SRS-GPR方法与其它方法的边坡可靠度分析结果比较。由图 11~13可见,当参数ρc, ϕ,COVc,COVϕ和δv在不同范围内取值时,本文方法计算的Pf与精确解基本保持一致。因此,所提的SRS-GPR方法能够高效、准确地评估不同土体空间变异程度的高维边坡可靠度问题。
4. 结论
将强度折减采样、高斯过程回归和主动学习策略有机结合,提出一种可在原始物理参数空间处理高维变量的空间变异边坡自适应可靠度分析方法。通过不排水饱和黏性土坡和黏性-摩擦土坡两个算例验证所提方法的有效性,得到4点结论。
(1)基于强度折减理论,推导边坡极限平衡状态下的强度参数公式,提出强度折减采样法,该方法无需进行边坡确定性分析,便可有效生成极限状态面上的临界样本点。
(2)融合强度折减采样和主动学习策略,可自适应搜索最优样本点,有效提高GPR模型在极限状态面附近的预测精度,从而避免因样本点过度集中造成信息冗余,显著提高可靠度分析效率。
(3)算例结果表明,相较于其它方法,所提方法不受随机场模拟方法的影响,可直接在原始高维物理参数空间应用,且边坡确定性分析次数少,可靠度计算效率高,Pf置信区间较窄。同时,所提方法在不同工况下,仍能给出与直接MCS法相吻合的Pf结果。上述结果均表明所提方法具备较好的准确性、高效性、鲁棒性和适用性,可为解决考虑参数空间变异性的高维边坡可靠度问题提供一条有效途径。
(4)在验证方法有效性时,采用极限平衡法(简化Bishop法)求解边坡安全系数。需要指出的是,所提方法中的边坡确定性分析并不限于极限平衡法,边坡稳定性分析的有限元或有限差分法同样适用。
-
图 12 ATLAS Ⅲ现场测试水平示意图[24]
Figure 12. Horizontal view of ATLAS Ⅲ heating tests
图 13 AT89E中央加热钻孔现场试验图[25]
Figure 13. Field photo of central heating borehole AT89E
表 1 参数取值
Table 1 Values of parameters
参数 量值 参数 量值 λB/(W·(K·m)-1) 0.68 z2/m 502.625 λR/(W·(K·m)-1) 2.50 t0/℃ 10 cB/(J·(K·kg)-1) 1014 t1/℃ 34 cR/(J·(K·kg)-1) 832 P(0)/W 1545 ρB/(kg·m-3) 1 770 ϕ 0.91 ρR/(kg·m-3) 2 630 a1 0.049 a/m 0.525 a2 0.696 b/m 0.825 a3 -0.059 h/m 5.25 a4 0.271 L/m 1000 a5 0.027 z0/m 0 a6 -0.010 z1/m 487.375 a7 0.026 表 2 ATLAS Ⅲ试验加热和冷却信息
Table 2 Heating and cooling information of ATLAS Ⅲ heating tests
步骤 阶段 热功率/W 日期 天数序号 持续/d 1 加热 0→400 2007-04-02→2007-04-05 0→4 4 2 稳定 400 2007-04-06→2007-05-20 4→49 45 3 加热 400→900 2007-05-21→2007-05-25 49→54 4 4 稳定 900 2007-05-26→2007-07-29 54→120 66 5 加热 900→1400 2007-07-30→2007-08-03 120→125 5 6 稳定 1400 2007-08-04→2008-04-16 125→381 256 7 冷却 0 2008-04-17→2009-11-02 381→945 564 表 3 参数取值
Table 3 Values of parameters
参数 量值 参数 量值 λG/(W·(K·m)-1) 16.3 ρBoom/(kg·m-3) 2000 λBoom/(W·(K·m)-1) 1.53 a/m 0.16 cG/(J·(K·kg)-1) 500 b/m 0.19 cBoom/(J·(K·kg)-1) 1434 t0/℃ 16 ρG/(kg·m-3) 7850 注:表中下标G和Boom分别表示不锈钢钢管和Boom黏土。 -
[1] 王驹, 陈伟明, 苏锐, 等. 高放废物地质处置及其若干关键科学问题[J]. 岩石力学与工程学报, 2006, 25(4): 801-812. doi: 10.3321/j.issn:1000-6915.2006.04.015 WANG Ju, CHEN Weiming, SU Rui, et al. Geological disposal of high-level radioactive waste and its key scientific issues[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(4): 801-812. (in Chinese) doi: 10.3321/j.issn:1000-6915.2006.04.015
[2] 叶为民, 王琼, 潘虹, 等. 高压实高庙子膨润土的热传导性能[J]. 岩土工程学报, 2010, 32(6): 821-826. http://cge.nhri.cn/cn/article/id/13418 YE Weimin, WANG Qiong, PAN Hong, et al. Thermal conductivity of compacted GMZ01 bentonite[J]. Chinese Journal of Geotechnical Engineering, 2010, 32(6): 821-826. (in Chinese) http://cge.nhri.cn/cn/article/id/13418
[3] 陈航, 张虎元, 郭永强, 等. 混合型缓冲回填材料导热性能测试与预测研究[J]. 岩石力学与工程学报, 2014, 33(增刊2): 4312-4320. CHEN Hang, ZHANG Huyuan, GUO Yongqiang, et al. Measurement and prediction of thermal properties of bentonite-sand mixtures as buffer backfilling materials for high level radioactive waste[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(S2): 4312-4320. (in Chinese)
[4] 谈云志, 李辉, 王培荣, 等. 膨润土受热作用后的水-力性能研究[J]. 岩土力学, 2019, 40(2): 489-496. TAN Yunzhi, LI Hui, WANG Peirong, et al. Hydro-mechanical performances of bentonite respond to heat-treated history[J]. Rock and Soil Mechanics, 2019, 40(2): 489-496. (in Chinese)
[5] ZHANG J R, SUN D A, YU H H, et al. Swelling of unsaturated GMZ07 bentonite at different temperatures[J]. Bulletin of Engineering Geology and the Environment, 2020, 79(2): 959-969. doi: 10.1007/s10064-019-01595-y
[6] 陈正汉, 秦冰. 缓冲/回填材料的热-水-力耦合特性及其应用[M]. 北京: 科学出版社, 2017. CHEN Zhenghan, QIN Bing. Thermal-Hydraulic-Mechanical Coupling Characteristics of Buffer/Backfill Materials and its Application[M]. Beijing: Science Press, 2017. (in Chinese)
[7] 陈正汉, 郭楠. 非饱和土与特殊土力学及工程应用研究的新进展[J]. 岩土力学, 2019, 40(1): 1-54. CHEN Zhenghan, GUO Nan. New developments of mechanics and application for unsaturated soils and special soils[J]. Rock and Soil Mechanics, 2019, 40(1): 1-54. (in Chinese)
[8] 陈正汉. 非饱和土与特殊土力学[M]. 北京: 中国建筑工业出版社, 2022. CHEN Zhenghan. Mechanics for Unsaturated and Special Soils[M]. Beijing: China Architecture & Building Press, 2022. (in Chinese)
[9] CHO W J, KIM G Y. Reconsideration of thermal criteria for Korean spent fuel repository[J]. Annals of Nuclear Energy, 2016, 88: 73-82. doi: 10.1016/j.anucene.2015.09.012
[10] 刘月妙, 王驹, 蔡美峰, 等. 热-力耦合条件下高放废物处置室间距研究[J]. 铀矿地质, 2009, 25(6): 373-379. doi: 10.3969/j.issn.1000-0658.2009.06.009 LIU Yuemiao, WANG Ju, CAI Meifeng, et al. Study on disposal pit space for high-level radioactive waste in thermal-mechanical coupling conditiion[J]. Uranium Geology, 2009, 25(6): 373-379. (in Chinese) doi: 10.3969/j.issn.1000-0658.2009.06.009
[11] CARSLAW H S, JAEGER J C. Conduction of Heat in Solids[M]. 2d ed. Oxford: Clarendon Press, 1959.
[12] CLAESSON J, PROBERT T. Thermoelastic stress due to a rectangular heat source in a semi-infinite medium. Presentation of an analytical solution[J]. Engineering Geology, 1998, 49(3/4): 223-229.
[13] HÖKMARK H, FÄLTH B. Thermal dimensioning of the deep repository[R]. Stockholm: Svensk Kärnbränslehantering AB, 2003.
[14] HÖKMARK H, LÖNNQVIST M, KRISTENSSON O. Strategy for thermal dimensioning of the final repository for spent nuclear fuel[R]. Stockholm: Svensk Kärnbränslehantering AB, 2009.
[15] 刘东东, 项彦勇. 高放射核废处置库温度场的分布线热源解析模型[J]. 岩石力学与工程学报, 2019, 38(增刊1): 2816-2822. LIU Dongdong, XIANG Yanyong. A distributed line heat-source analytical model for the temperature field of a high level nuclear waste repository[J]. Chinese Journal of Rock Mechanics and Engineering, 2019, 38(S1): 2816-2822. (in Chinese)
[16] IKONEN K. Thermal analysis of repository for spent EPR-type fuel[R]. Olkiluoto: Posiva Oy, 2005.
[17] HARTLEY L, HOCH A, JACKSON P, et al. Groundwater flow and transport modeling during the temperate period for the SR-Can assessment. Forsmark area-version 1.2[R]. Stockholm: Svensk Kärnbränslehantering AB, 2006.
[18] 陈永贵, 贾灵艳, 叶为民, 等. 施工接缝对缓冲材料水-力特性影响研究进展[J]. 岩土工程学报, 2017, 39(1): 138-147. doi: 10.11779/CJGE201701012 CHEN Yonggui, JIA Lingyan, YE Weimin, et al. Advances in hydro-mechanical behaviors of buffer materials under effect of technological gaps[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(1): 138-147. (in Chinese) doi: 10.11779/CJGE201701012
[19] XU X, HE L Q, SUN D A, et al. Validation of the fully-analytical solution with temperature superposition for the nuclear waste repository[J]. Nuclear Engineering and Design, 2023, 409: 112367. doi: 10.1016/j.nucengdes.2023.112367
[20] 陶文铨. 传热学[M]. 5版. 北京: 高等教育出版社, 2019. TAO Wenquan. Heat Transfer[M]. 5th ed. Beijing: Higher Education Press, 2019. (in Chinese)
[21] HE L Q, ZHOU X Y, SUN D A. Fully analytical solution in time and space domains on temperature in multi-barrier nuclear waste repository[J]. Computers and Geotechnics, 2023, 154: 105164. doi: 10.1016/j.compgeo.2022.105164
[22] CRUMP K S. Numerical inversion of Laplace transforms using a Fourier series approximation[J]. Journal of the ACM, 1976, 23(1): 89-96. doi: 10.1145/321921.321931
[23] 徐国庆. 国际高放废物处置研发工作在花岗岩地区的进展[J]. 世界核地质科学, 2016, 33(2): 119-124. doi: 10.3969/j.issn.1672-0636.2016.02.010 XU Guoqing. Abroad progress in R & D work on high-level radioactive waste disposal in granite areas[J]. World Nuclear Geoscience, 2016, 33(2): 119-124. (in Chinese) doi: 10.3969/j.issn.1672-0636.2016.02.010
[24] CHEN G J, SILLEN X, VERSTRICHT J, et al. ATLAS Ⅲ in situ heating test in boom clay: field data, observation and interpretation[J]. Computers and Geotechnics, 2011, 38(5): 683-696. doi: 10.1016/j.compgeo.2011.04.001
[25] DE BRUYN D, LABAT S. The second phase of ATLAS: the continuation of a running THM test in the HADES underground research facility at Mol[J]. Engineering Geology, 2002, 64(2/3): 309-316.
-
期刊类型引用(0)
其他类型引用(1)
-
其他相关附件