Numerical model for hydro-mechanical-damage coupling of rocks based on TOUGHREACT
-
摘要: 从细观力学角度出发,充分考虑水-力耦合条件下岩石细观特征及其演化,结合热力学理论,建立基于TOUGHREACT的岩石细观水力损伤耦合数值模型。模型可较好地考虑任意微裂纹滑移剪胀、损伤扩展和法向压缩闭合等细观力学行为对岩石宏观变形破坏、渗透性演化和水流运动过程的影响。采用室内煤岩注水破坏试验成果对数值模型的正确性和有效性进行验证,进而开展现场尺度下岩石注水响应的应用模拟研究。模拟结果表明,注水引起的岩石损伤与压力增高区的分布同时受注入流量、现场应力水平与初始微裂纹各向异性分布等因素的影响。岩石宏观水力耦合响应的模拟有赖于内部微裂纹结构的准确表征。研究成果对深化岩石水力耦合研究具有一定参考意义。Abstract: For better modelling the microscopic characteristics and evolution of materials under hydro-mechanical coupling conditions, a TOUGHREACT-based hydro-mechanical-damage coupled numerical model for saturated rocks is established by using the microscopic homogenization method and the thermodynamics theory. The proposed model well accounts for the influences of sliding dilatancy, damage propagation and normal compression of arbitrary microcracks on the macroscopic deformation and failure characteristics, permeability evolution and fluid flow process. The numerical method is successfully validated through the experimental data of water injection tests on coal sample at the laboratory scale and then used to carry out application simulations of water injection responses at the field scale. The numerical simulation results demonstrate that the distributions of injection-induced rock damage and elevated pressure are affected by the injection rate, in-situ stresses and anisotropic distribution of the initial microcracks, and they are more developed in the directions with larger in-situ stress and dominant development of microcracks. Better simulations of the macroscopic hydro-mechanical responses of rocks depend on the accurate characterization of the internal microscopic structures. The research may provide a useful reference for deepening the study on hydro-mechanical coupling of rocks.
-
0. 引言
渗流计算是评价水工建筑物渗流安全及进行渗控方案设计的依据。渗流计算的实践表明有限元方法是一种非常成熟、非常有效的数值计算方法。有限元方法分为协调有限元方法和非协调有限元方法,要求近似函数在剖分单元边界上连续的称为协调有限元方法,允许近似多项式函数在单元边界上不连续的称为非协调有限元方法,又称为间断有限元方法。通常讲的有限元方法是指协调有限元方法,称为一般有限元方法。一般有限元方法连续性的要求,使得剖分网格只能使用单一的单元类型,同一次数的近似多项式函数,因而一般有限元方法不适用于混合剖分网格。在渗流计算中,往往希望通过加密网格的方法来提高计算精度,而这有可能导致计算格式不稳定、计算结果不收敛。国内外众多学者一直在探索适用于混合剖分网格的间断有限元方法,主要成果有:内罚间断有限元法(IPDG)[1],局部间断有限元法(LDG)[2]和可杂交的间断有限元法(HDG)[3]。与一般有限元法比较,非协调有限元方法有两方面的优势:①近似函数灵活,可以在不同剖分单元采用不同次数的近似多项式函数;②网格生成灵活,适用混合剖分网格。但是,间断有限元方法的变分形式比协调有限元方法的变分形式要复杂得多,计算格式的稳定性往往需要靠调整参数来保证、而辅助变量的引入和变分的复杂化则使得单元矩阵非局部化,这些困难使得间断有限元方法难以在渗流分析中充分发挥其优势[4]。弱有限元法(WG)是Wang等[5]提出的,它属于非协调有限元方法,通过在剖分单元上定义弱函数和弱梯度,从而用弱梯度算子替代经典的梯度算子,并在计算格式中引入了一个不含有参数的稳定子或光滑子,保证了格式的绝对稳定性[5-7]。因此,弱有限元方法不但具有间断有限元方法的优点,还克服了间断有限元方法的缺点,目前,弱有限元方法可以用来求解Stokes方程、Helmholtz方程、Maxwell方程、线弹性方程等[8-12]。
本文主要研究弱有限元法针对混合剖分网格在渗流分析上的具体应用。安排如下:在第一节中将介绍弱有限元方法的原理,并通过一个简例介绍弱有限元方法处理混合网格的流程;在第二节中运用弱有限元方法,采用混合剖分网格求解均值土坝渗流自由面和闸基渗流场。
1. 弱有限元方法
1.1 变分形式
稳定渗流的数学模型为椭圆型方程边值问题:
{−∇(p∇u)=fin Ωu=con Γ1∂u∂n=gon Γ2, (1) 式中:Γ1∪Γ12=∂Ω,n为Γ2上的外法线单位向量;u为水头函数,p为渗透系数矩阵。此问题的一般有限元法的变分形式为:求u∈H1(Ω),使得在Γ1上u=c,且满足
a(u,v)=(f,v)+⟨g,v⟩Γ2 ∀v∈H10(Ω)。 式中,a(u,v)=(p∇u,∇v)。在弱有限元法中,允许近似函数在单元边界上不连续,为此引入弱函数和弱函数的梯度算子。假设Th是区域Ω的一个剖分,剖分单元可以是三角形、四边形或其它多边形。单元T上的弱函数v定义为
v={v0,vb}={v0in Tvbon ∂T。 即,弱函数v在单元T上包含内部函数v0和边界函数vb两部分。用Pm(T)表示定义在T上且次数不超过m的多项式函数,弱函数v的弱梯度算子∇wv定义为对∀ττ∈[Pm(T)]2=(Pm(T),Pm(T)),有
(∇wv,ττ)=−(v0,∇⋅ττ)T+⟨vb,ττ⋅n⟩∂T。 (2) 单元T上的弱有限元空间定义为
Wj,l(T)={v={v0,vb}|v0(T)∈Pj(T),vb(e)∈Pl(e),e∈∂T}, 弱梯度算子是由Wj,l(T)到[Pm(T)]2上的一个映射。剖分Th上的弱有限元法空间定义为
Vh={v|v(T)∈Wj,l(T),∀T∈Th}。 令
V0h={v={v0,vb}∈Vh,vb|Γ1=0,}, 问题(1)的弱有限元法的变分形式为,求uh∈Vh,使得在Γ1上uh=c,且满足
ah(uh,v)+s(uh,v)=(f,v)+⟨g,v⟩Γ2, ∀v∈V0h。 (3) 其中ah(uh,v)=(p∇wuh,∇wv),s(uh,v)是为了保证格式稳定的一个补充项,称为稳定子或光滑子,且
s(uh,v)=∑T∈Th(J∂(uh),J∂(v))∂T。 (4) 式中,J∂(v)=h−1/2T(Q(v0)−vb),hT表示单元T的直径,即T上任意两点间距离的最大值(编程计算时一般取为多边形单元任意两节点间距离的最大值),Q(v0)表示内部函数v0在边界上的投影。Θ(∇u0)表示∇u0在单元边界∂T上的投影。根据弱梯度的定义和分部积分公式计算得
(∇wuh,∇wv)T=(∇u0,∇v0)T−⟨Θ(∇u0)⋅n,Q(v0)−vb⟩∂T−⟨Q(u0)−ub,∇wv⋅n⟩∂T, 注意到
Θ(∇u0)⋅n≈(Q(u0)−ub)/hT, ∇wv⋅n≈(Q(v0)−vb)/hT, 故
(∇u0,∇v0)T≈(∇wuh,∇wv)T+2h−1T⟨Q(u0)−ub,Q(v0)−vb⟩∂T。 因此在变分式中用ah(uh,v)近似a(u,v),把s(uh,v)作为补偿项加入变分式中是合理的[7]。
1.2 计算格式
弱有限元空间Wj,l(T)表示在单元T内部近似函数采用j次多项式,在单元边界∂T上近似函数采用l次多项式。[Pm(T)]2表示弱函数的弱梯度是分量为m次多项式的向量。选取不同的j,l和m,会得到单元局部矩阵不同的计算格式,Wj,l(T)和[Pm(T)]2简记为Wj,l(T)/[Pm(T)]2,下面假设单元T是一个三角形单元(图 1)。
若选取弱有限元空间W1,0(T)/[P0(T)]2,取基函数为
φ1={1,in T,0,on ∂T.,φ2={x,in T,0,on ∂T.,φ3={y,in T,0,on ∂T., ϕ3+i={0,其它,1,on ei. (i=1,2,3) 。 则,
Wj,l(T)=span{φ1,φ2,⋯,φ6}, uh(T)={uh0(T),uhb(T)}=6∑i=1aTiφi。 变分式(3)右端项,
(f,v)T=(f,v0)T, (5) ⟨g,v⟩∂T∩Γ2=∑e∈∂T∩Γ2⟨g,vb⟩e, (6) 依次取v=φi,i=1,2,⋯,6,由式(5)和(6)可得单元常数项矩阵.单元局部系数矩阵AT包含两部分,
\boldsymbol{A}_T=\boldsymbol{M}_T+\boldsymbol{S}_T。 (7) 其中,
\boldsymbol{M}_T=\left(\nabla_{\mathrm{w}} \varphi_i, \nabla_{\mathrm{w}} \varphi_j\right)_{6 \times 6}, {\mathit{\boldsymbol{S}}_T} = {\left( {\sum\limits_{i = 1}^3 {{{\left\langle {{J_\partial }\left( {{\varphi _i}} \right),{J_\partial }\left( {{\varphi _j}} \right)} \right\rangle }_{{e_i}}}} } \right)_{6 \times 6}}{\rm{ }}。 根据弱梯度的定义(2),由 \boldsymbol{\tau} \in\left[P_0(T)\right]^2 ,有 {\nabla _{\text{w}}}{\varphi _i} = 0, \nabla_{\mathrm{w}} \varphi_{3+i}=\frac{\left|e_i\right|}{|T|} \boldsymbol{n}_i, i=1,2,3 ,其中, |{e_i}| 为单元边界 {e_i} 的长度, |T| 为单元T的面积,ni为单元边界 {e_i} 的外法线单位向量。令
\boldsymbol{M}_T=\left(\begin{array}{cc} O & O \\ O & \boldsymbol{M}_{b b} \end{array}\right) (8) 这里 M_{b b}=\left(\frac{\left|e_i \| e_j\right|}{|T|} \boldsymbol{n}_i \boldsymbol{n}_j\right)_{3 \times 3} 。ST通过数值积分计算,用 ({x_i}, {y_i}), i = 1, 2, 3 表示单元边界 {e_i} 的中点坐标,则
{\mathit{\boldsymbol{S}}_T} = \sum\limits_{i = 1}^3 {\left| {{e_i}} \right|\mathit{\boldsymbol{N}}_i^{\rm{T}}{\mathit{\boldsymbol{N}}_i}} (9) 其中
\boldsymbol{N}_1=h_T^{-1 / 2}\left(1, x_1, y_1,-1,0,0\right) \text {, } \boldsymbol{N}_2=h_T^{-1 / 2}\left(1, x_2, y_2, 0,-1,0\right) \text {, } \boldsymbol{N}_3=h_T^{-1 / 2}\left(1, x_3, y_3, 0,0,-1\right) 。 式(5)~(9)称为弱有限元法的计算格式。对于四边形单元的计算格式参考文献[6],其他多边形的计算格式按照这个过程也易于推得。
由于相邻单元内部函数之间没有联系,所以可以选取不同次数的近似多项式函数,因而弱有限元方法适用于混合剖分网格,对于边中点网格,如剖分图 2中第4号单元是个四边形,一条边中间有个节点,那么这条边上的近似函数可以理解为分段近似,这在局部矩阵计算格式上相当于两条边,这个点把该边一分为二,这个单元就可以被认定为五边形。类推,如果一条边中间有多个节点,就可把该条边分成多条边,该单元就被认定为相应类型的多边形单元。因此,对于弱有限元方法,网格剖分比较自由。
1.3 混合剖分网格应用简例
下面通过一个简单的例子来了解弱有限元法处理混合网格的计算过程和计算效果。
实例1:考虑
\left\{\begin{array}{rc} -\nabla \cdot(\nabla u)=0, & (x, y) \in[0,2] \times[0,2] \\ u=1, & x=0 \\ u=0, & x=2 \\ \frac{\partial u}{\partial \boldsymbol{n}}=0, & y=0 \text { 或 2。 } \end{array}\right. (10) 对正方形区域 \Omega = [0, 2] × [0, 2] 做剖分,如图 2所示。区域被分成4个单元,其中第1、2号单元为三角形,第3号单元为四边形,第4号单元有一条边中间有一个节点,因而该单元被认定为五边形,节点共8个,单元边界共11条。单元信息见表 1,各节点坐标信息见表 2。
表 1 剖分单元信息Table 1. Information of elements单元编号 单元类型 节点编号 1 3 3 2 4 2 3 3 1 2 3 4 3 4 5 6 4 5 3 6 7 8 1 表 2 节点坐标信息Table 2. Information of node coordinates节点 1 2 3 4 5 6 7 8 x 1 2 1 2 2 1 0 0 y 0 0 1 1 2 2 2 0 每个单元的弱有限元空间取 {W_{1, 0}}(T)/{[{P_0}(T)]^2} ,根据式(7)~(9)计算各单元的局部系数矩阵AT,用 \mathit{\boldsymbol{A}}_T^i 表示第 i 号单元的局部系数矩阵,则 \mathit{\boldsymbol{A}}_T^1,\mathit{\boldsymbol{A}}_T^2 为6阶方阵, \mathit{\boldsymbol{A}}_T^3,\mathit{\boldsymbol{A}}_T^4 分别为7阶、8阶方阵。为了把单元局部矩阵组合成整体矩阵,需要把所有自由度排序,每个单元的内部近似函数有3个自由度,每条单元边界对应一个自由度,为此先给所有单元边界编号如下(图 3)。
如果把所有单元内部近似函数的自由度按单元序号排前面,即,第 i 号单元内部近似函数的3个自由度序号依次为 3(i - 1) + 1 , 3(i - 1) + 2 , 3(i - 1) + 3 ,把边界自由度按序号排在所有单元内部近似函数自由度序号的后面,比如,第1号边界的自由度排在第3×4+1=13号,第(8)号边界的自由度排在第3×4+8=20号。对于这个剖分网格,一共有23个自由度。依据自由度序号把单元局部矩阵组合成整体系数矩阵A。例如,第2号单元的局部矩阵的元素 \mathit{\boldsymbol{A }}_T^2(3, 6) ,它是第2号单元上的基函数 {\varphi _3} 和 {\varphi _6} 作用的结果,它们对应的自由度序号分别为6和13,所以应该累加到总体系数矩阵的A(6,13),把所有局部单元矩阵组合成整体系数矩阵A,A是一个23阶的方阵,此问题的总体代数方程组 \mathit{\boldsymbol{A X}}{\rm{ = }}\mathit{\boldsymbol{b}} ,由于此例 f = 0 ,且在第二类边界外法线方向上的方向导数为零,所以此时常数项矩阵b是一个23×1阶的零矩阵,未知数X是由所有自由度构成的一个23×1阶矩阵。问题的第二类边界条件属于自然边界条件,已在变分式中体现了,接下来需要处理第一类边界条件。处理第一类边界条件有不同的方法,这里采用如下处理方法:由于第2号、第6号和第10号单元边界属于第一类边界,它们对应的自由度序号分别为14,18,22,第2号和第6号边界上函数值为0,第10号边界上函数值为1,故对总体系数矩阵和常数项矩阵做如下修正:
\begin{aligned} & \boldsymbol{A}(14, :)=0, \boldsymbol{A}(14, 14)=1, \boldsymbol{b}(14)=0, \\ & \boldsymbol{A}(18, :)=0, \boldsymbol{A}(18, 18)=1, \boldsymbol{b}(18)=0, \\ & \boldsymbol{A}(22, :)=0, \boldsymbol{A}(22, 22)=1, \boldsymbol{b}(22)=1 \text { 。 } \end{aligned} 其中, \mathit{\boldsymbol{A}}(i, :) 表示矩阵A的第i行的所有元素。处理边界条件后就得到该问题的弱有限元法的最终的总体线性方程组 \mathit{\boldsymbol{A X}}{\rm{ = }}\mathit{\boldsymbol{b}} 。解总体线性方程组,得每个单元的内部近似函数,结果都为
{u_{h0}}(T) = 1.002 - 0.5001x 。 而此问题的解析解为 u = 1 - 0.5x 。这说明弱有限元法处理混合剖分网格是可行,且有效的。
2. 渗流分析
2.1 自由面问题
在岩土边坡、土坝、地下洞室及地下水运动等渗流分析中,均存在有渗流自由面问题。渗流自由面的确定是一个重点,也是个难点。求解渗流自由面的方法有虚单元法[13]、截止负压法[14]、高斯点法[15]、变渗透系数法[16]等,这些方法都是以一般有限元方法为基础,因而精度都受制于一般有限元方法的不足。现在用弱有限元法来求解,看看计算结果的精度如何。
计算模型1:假设有一10 m×10 m的均质土坝,上游水位10 m,下游水位2 m,底部为不透水边界。现把分析区域进行剖分,剖分网格为网格Ⅰ(图 4):区域被均匀划分成200个直角边为1 m的等腰直角三角形单元,一共有320条单元边界。坝体上游和下游水位以下部分为第一类边界,初始计算时假定溢出段为: x=10,\text{ }2\le y\le 5 。弱有限元空间取 {W_{1, 0}}(T) / {\left[ {{P_0}(T)} \right]^2} ,由于溢出段也是第一类边界,溢出段上水头 u = y ,对于落在溢出段上的单元边界 {e_i} ,需做特殊处理,该单元边界上对应的基函数应取为
\varphi_{3+i}=\left\{0, \varphi_{b i}\right\}=\left\{\begin{array}{cc} y & (x, y) \in e_i \\ 0 & \text { 其他 } \end{array}\right. \text { 。 } 通过计算得
\nabla_{\mathrm{w}} \varphi_{3+i}=\frac{\left|e_i\right|}{|T|} y_i \boldsymbol{n}_i, \quad J_{\partial}\left(\varphi_{3+i}\right)=-h_T^{-1 / 2} y_i \text { 。 } 所以单元局部矩阵中的Mbb和ST部分的计算公式需做相应调整。比如,假设三角形单元(图 1)的边界 {e_2} 属于溢出段上, {e_1} 和 {e_3} 属于内部边界,则该单元局部矩阵
\boldsymbol{A}_T=\boldsymbol{M}_T+\boldsymbol{S}_T, 其中
\boldsymbol{M}_T=\left(\begin{array}{cc} O & O \\ O & \boldsymbol{M}_{b b} \end{array}\right), \begin{aligned} & \boldsymbol{M}_{b b}=\frac{1}{|T|}\left(\begin{array}{ccc} \left|e_1\right|^2 & \left|e_1 \| e_2\right| y_2 \boldsymbol{n}_2 \boldsymbol{n}_1 & \left|e_1 \| e_3\right| \boldsymbol{n}_3 \boldsymbol{n}_1 \\ \left|e_1 \| e_2\right| y_2 \boldsymbol{n}_2 \boldsymbol{n}_1 & \left|e_2\right|^2 y_2^2 & \left|e_2 \| e_3\right| y_2 \boldsymbol{n}_2 \boldsymbol{n}_3 \\ \left|e_1 \| e_3\right| \boldsymbol{n}_3 \boldsymbol{n}_1 & \left|e_2 \| e_3\right| y_2 \boldsymbol{n}_2 \boldsymbol{n}_3 & \left|e_3\right|^2 \end{array}\right) \\ & \boldsymbol{S}_T=\sum\limits_{i=1}^3\left|e_i\right| \boldsymbol{N}_i^T \boldsymbol{N}_i, \end{aligned} 其中,
\boldsymbol{N}_1=h_T^{-1 / 2}\left(1, x_1, y_1,-1,0,0\right) \text {, } \boldsymbol{N}_2= h_T^{ - 1/2}\left(1, {x_2}, {y_2},0,- {y_2},0\right) \text{,} \boldsymbol{N}_3 = h_T^{ - 1/2}\left(1, {x_3}, {y_3}, 0, 0, - 1\right) \text{。} 处理第一类边界条件时,该边界 {e_2} 上自由度的值为1。根据变单元渗透系数法[16]的原理,先假设整个坝体渗透系数为1,用弱有限元法初始计算后,根据计算结果把自由面以上单元的渗透系数赋值为1/1000,再进行计算,循环计算3次后得数值结果(表 3),表中数据精确到0.01。
表 3 计算结果对比ⅠTable 3. Comparison of calculated results Ⅰ横坐标 解析解 FEM 误差1 WG 误差2 0 10.00 10.00 0.00 10.00 0.00 1 9.59 9.73 0.14 9.62 0.03 2 9.17 9.39 0.22 9.26 0.09 3 8.72 8.99 0.27 8.79 0.07 4 8.25 8.53 0.28 8.30 0.05 5 7.75 8.03 0.28 7.78 0.03 6 7.21 7.46 0.25 7.23 0.02 7 6.63 6.83 0.20 6.63 0.00 8 6.00 6.11 0.11 6.04 0.04 9 5.29 5.21 -0.08 5.34 0.05 10 4.47 4.50 0.03 4.50 0.03 确定溢出点的方法:根据计算结果得到各节点的水头值后,搜索在坝体下游坝面( x = 10, y \geqslant 2 )上的单元边界 {e_i} ,计算 {e_i} 的两个节点处的纵坐标值与水头值的差值,如果两个差值异号,则表明出溢点在该单元边界内,再利用计算结果求该单元内部近似函数在边界 x = 10 上的投影,解方程 y = {Q_{\text{b}}}({u_{\text{h}}}(T)) = {a_{T1}} + {a_{T2}} \times 10 + {a_{T3}}y 即可求得溢出点的位置。自由面图见图 5。该问题的自由面有解析解: {y^2} = 100 - 8x ,参考文献[17]中有针对剖分网格Ⅰ,利用虚单元法求解的自由面数据结果,虚单元法是以有限元方法为基础,所以在表 3中对应FEM列。
结果对比显示,在相同剖分网格下,弱有限元法的最大绝对误差仅为0.09,最大相对误差为0.0098,而一般有限元法的最大绝对误差为0.28,最大相对误差为0.0361。这表明弱有限元法的计算结果比一般有限元法的精度明显要高。弱有限元法不足的地方在于:按照1.2节算例中第一类边界条件的处理方法,一般有限元法的自由度个数为121个,而弱有限元法的自由度个数为3×200+320=920个。为了发挥弱有限元法适用于混合剖分网格的优势,减小弱有限元法的自由度个数,根据自由面的预估位置,对区域进行有针对性的剖分(见图 6),在自由面附近区域采用边长为1 m的正方形网格,其它部位采用粗网格单元,一共剖分成57个单元,其中52个四边形,2个五边形和3个六边形,单元边界有133条,弱有限元空间仍然选取 {W_{1, 0}}(T)/{[{P_0}(T)]^2} 。自由度个数减少为3×57+133 =304个。仍然循环计算3次后得计算结果(见表 4)
表 4 计算结果对比ⅡTable 4. Comparison of calculated results Ⅱ横坐标 解析解 WG 误差 0 10.00 10.00 0.00 1 9.59 9.59 0.00 2 9.17 9.18 0.01 3 8.72 8.77 0.05 4 8.25 8.23 -0.02 5 7.75 7.74 -0.01 6 7.21 7.16 -0.05 7 6.63 6.61 -0.02 8 6.00 6.00 0.00 9 5.29 5.34 0.05 10 4.47 4.50 0.03 总体方程组的自由度大大地减少了,而数值解的精度却没有降低,最大绝对误差仅为0.05,最大相对误差为0.0095。这说明,弱有限元法能处理混合网格这一优势是非常显著的。
2.2 闸基渗流场
弱有限元法适用于混合剖分网格,那么在进行复杂的渗流场分析时,可以根据需要采用有针对性的剖分网格,这样可以用较少的自由度得到精度较高的数值计算结果。下面用弱有限元法来分析有防渗帷幕的闸基渗流场。
计算模型2:模型长和高分别为60,20 m,闸基底板长20 m,上游20 m,下游20 m,上游水位深4 m,下游水位默认为0 m,闸基底部默认为不透水边界,在闸基底部靠上游端下方有1 m厚防渗墙(见图 7)。
闸基渗流分析重点关心防渗墙下端处和闸基底板下游出口处的水力梯度,因此采用有针对性的混合剖分网格(图 8)。
防渗墙周围和闸基底板下游出口处两个部位采用细网格的四边形单元,其它部位采用大网格四边形单元,大网格与细网格相邻的大网格被认定为八边形单元,因此,单元类型共两种:四边形和八边形,剖分单元一共有192个单元,其中四边形180个,八边形12个,节点共241个,有单元边界432条。弱有限元空间取 {W_{10}}(T)/{[{P_0}(T)]^2} ,闸基渗流不含有需要特殊处理的边界,所以计算格式按1.2节中介绍的公式计算即可。自由度总数为3×192+432=1008个。
防渗墙和地基黏土的渗透系数比假定为1∶1000,对防渗墙深5,10,15 m 3种工况采用弱有限元法进行了渗流计算,计算结果的水头等值线分布如图 9~11(单位:m),防渗墙底端和闸基出口处的水力梯度值见表 5。
表 5 水力梯度值Table 5. Values of hydraulic gradient防渗墙
深/m防渗墙底端处水力梯度 闸基出口处
水力梯度5 0.323 0.332 10 0.333 0.262 15 0.358 0.190 计算结果表明随着防渗墙的加深下游水头值逐渐降低,防渗墙的防渗效果得到充分体现,也说明了采用有针对性的混合剖分网格,运用弱有限元法来处理类似闸基这种比较复杂的渗流场问题是非常有效的。
3. 结论
计算渗流自由面模型的数值结果表明弱有限元法的计算结果的精度比一般有限元法的要高,并且采用较稀疏的混合剖分网格也能得到精度很好的数值结果。利用弱有限元法适用混合剖分网格的特点,对较复杂的闸基渗流场进行分析,根据数值结果绘制的等水头线表明结果符合闸基渗流的特征。高精度和能适用于混合剖分网格这两个优势,使得弱有限元法在渗流分析中有很好的应用前景。
-
表 1 注水试验模拟的计算参数
Table 1 Parameters used in water injection experiment
Es/GPa V(dc) /MPa b0 ks /m2 dc β0 /m2 Δt/s 4.0 3×10-4 0.8 1.0×10-18 1.5 0.3 1.0 0.01 1.4×10-15 30 表 2 现场尺度注水响应模拟的计算参数
Table 2 Parameters used for water injection modeling in field scale
Es/GPa V(dc)/MPa b0 ks /m2 dc β0 /m2 T/s 10.0 3×10-4 0.8 1.0×10-16 1.9 0.3 0.5 0.01 1.4×10-13 480 -
[1] 陈益峰, 胡冉, 周创兵, 等. 热-水-力耦合作用下结晶岩渗透特性演化模型[J]. 岩石力学与工程学报, 2013, 32(11): 2185-2195. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201311003.htm CHEN Yi-feng, HU Ran, ZHOU Chuang-bing, et al. A permeability evolution model for crystalline rocks subjected to coupled thermo-hydro- mechanical loading[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 32(11): 2185-2195. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201311003.htm
[2] WU G J, CHEN W E, RONG C, et al. Elastoplastic damage evolution constitutive model of saturated rock with respect to volumetric strain in rock and its engineering application[J]. Tunnelling and Underground Space Technology, 2020, 97: 103284. doi: 10.1016/j.tust.2020.103284
[3] 姚池, 姜清辉, 位伟, 等. 复杂裂隙岩体水-力耦合模型及溶质运移模拟[J]. 岩石力学与工程学报, 2013, 32(8): 1656-1665. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201308020.htm YAO Chi, JIANG Qing-hui, WEI Wei, et al. Numerical simulation of hydro-mechanical coupling and solute transport in complex fractured rock masses[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 32(8): 1656-1665. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201308020.htm
[4] 胡亚元. 基于混合物理论的饱和岩石弹塑性模型[J]. 岩土工程学报, 2020, 42(12): 2161-2169. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202012002.htm HU Ya-yuan. Elastoplastic model for saturated rock based on mixture theory[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(12): 2161-2169. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202012002.htm
[5] 陈益峰, 李典庆, 荣冠, 等. 脆性岩石损伤与热传导特性的细观力学模型[J]. 岩石力学与工程学报, 2011, 30(10): 1959-1969. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201110003.htm CHEN Yi-feng, LI Dian-qing, RONG Guan, et al. A micromechanical model for damage and thermal conductivity of brittle rocks[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(10): 1959-1969. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201110003.htm
[6] JIANG T, SHAO J F, XU W Y, et al. Experimental investigation and micromechanical analysis of damage and permeability variation in brittle rocks[J]. International Journal of Rock Mechanics and Mining Sciences, 2010, 47(5): 703-713. doi: 10.1016/j.ijrmms.2010.05.003
[7] ZHU Q Z, SHAO J F. Micromechanics of rock damage: Advances in the quasi-brittle field[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2017, 9(1): 29-40. doi: 10.1016/j.jrmge.2016.11.003
[8] 刘武. 考虑多尺度结构的贯通节理岩体损伤摩擦耦合模型[J]. 岩土工程学报, 2018, 40(1): 147-154. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201801019.htm LIU Wu. Coupled damage and friction model for persistent fractured rocks considering multi-scale structures[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(1): 147-154. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201801019.htm
[9] CHEN Y F, HU S H, ZHOU C B, et al. Micromechanical modeling of anisotropic damage-induced permeability variation in crystalline rocks[J]. Rock Mechanics and Rock Engineering, 2014, 47(5): 1775-1791.
[10] 胡大伟, 朱其志, 周辉, 等. 脆性岩石各向异性损伤和渗透性演化规律研究[J]. 岩石力学与工程学报, 2008, 27(9): 1822-1827. doi: 10.3321/j.issn:1000-6915.2008.09.009 HU Da-wei, ZHU Qi-zhi, ZHOU Hui, et al. Research on anisotropic damage and permeability evolutionary law for brittle rocks[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(9): 1822-1827. (in Chinese) doi: 10.3321/j.issn:1000-6915.2008.09.009
[11] LIU W, CHEN Y F, HU R, et al. A two-step homogenization- based permeability model for deformable fractured rocks with consideration of coupled damage and friction effects[J]. International Journal of Rock Mechanics and Mining Sciences, 2016, 89: 212-226.
[12] DORMIEUX L, KONDO D. Micromechanics of damage propagation in fluid-saturated cracked media[J]. Revue Européenne De Génie Civil, 2007, 11(7/8): 945-962.
[13] XIE N, ZHU Q Z, SHAO J F, et al. Micromechanical analysis of damage in saturated quasi brittle materials[J]. International Journal of Solids and Structures, 2012, 49(6): 919-928.
[14] ZHU Q Z. Strength prediction of dry and saturated brittle rocks by unilateral damage-friction coupling analyses[J]. Computers and Geotechnics, 2016, 73: 16-23.
[15] 朱其志, 王岩岩, 仇晶晶, 等. 准脆性岩石水力耦合不排水多尺度本构模型[J]. 河海大学学报:自然科学版, 2018, 46(2): 165-170. https://www.cnki.com.cn/Article/CJFDTOTAL-HHDX201802014.htm ZHU Qi-zhi, WANG Yan-yan, QIU Jing-jing, et al. Multiscale hydro-mechanical constitutive model for qusi-brittle rocks under undrained condition[J]. Journal of Hohai University (Natural Science), 2018, 46(2): 165-170. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-HHDX201802014.htm
[16] XU T F, SPYCHER N, SONNENTHAL E, et al. TOUGHREACT Version 2.0: a simulator for subsurface reactive transport under non-isothermal multiphase flow conditions[J]. Computers & Geosciences, 2011, 37(6): 763-774.
[17] HU C, LEMARCHAND E, DORMIEUX L, et al. Quasi-isotropic Biot’s tensor for anisotropic porous rocks: experiments and micromechanical modelling[J]. Rock Mechanics and Rock Engineering, 2020, 53(19): 4031-4041.
[18] CHEN Y F, WEI K, LIU W, et al. Experimental characterization and micromechanical modelling of anisotropic slates[J]. Rock Mechanics and Rock Engineering, 2016, 49(9): 3541-3557.
[19] ZHU Q Z, SHAO J F. A refined micromechanical damage- friction model with strength prediction for rock-like materials under compression[J]. International Journal of Solids and Structures, 2015, 60/61: 75-83.
[20] BAZANT Z P, OH B H. Efficient numerical integration on the surface of a sphere[J]. ZAMM Journal of Applied Mathematics and Mechanics, 1986, 66(1): 37-49.
[21] WU C F, ZHANG X Y, WANG M, et al. Physical simulation study on the hydraulic fracture propagation of coalbed methane well[J]. Journal of Applied Geophysics, 2018, 150: 244-253.
[22] GAO Q, GHASSEMI A. Three dimensional finite element simulations of hydraulic fracture height growth in layered formations using a coupled hydro-mechanical model[J]. International Journal of Rock Mechanics and Mining Sciences, 2020, 125: 104137.
[23] YI L P, LI X G, YANG Z Z, et al. A fully coupled fluid flow and rock damage model for hydraulic fracture of porous media[J]. Journal of Petroleum Science and Engineering, 2019, 178, 814-828.
[24] 刘武, 张振华, 叶晓东, 等. 层状岩体渗透特性多尺度演化模型研究[J]. 岩土工程学报, 2018, 40(增刊2): 68-72. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC2018S2016.htm LIU Wu, ZHANG Zhen-hua, YE Xiao-dong, et al. Multi-scale permeability evolution model for layered rocks[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(S2): 68-72. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC2018S2016.htm
-
期刊类型引用(22)
1. 王智超,彭柱,彭峰,闫实. 聚氨酯固化钙质砂物理力学特性. 长江科学院院报. 2024(01): 107-113 . 百度学术
2. 韩庆华,王永超,刘铭劼,李浩斌. 振动台试验饱和机制砂模型土动力特性研究. 土木工程学报. 2024(03): 110-122 . 百度学术
3. 王家全,和玉,祝梦柯,钱弘毅. 相对密实度和固结应力比对北部湾海砂动力特性影响的试验研究. 安全与环境工程. 2024(04): 20-28 . 百度学术
4. 王钰轲,李俊豪,邵景干,余翔. 不同影响因素下路用黄河泥沙动剪切模量和阻尼比试验及理论模型研究. 工程科学学报. 2023(03): 509-519 . 百度学术
5. 胡健,肖杨,肖鹏,王林,丁选明,仉文岗,刘汉龙. 基于机器学习预测微生物加固钙质砂统一动强度. 中国公路学报. 2023(02): 80-88 . 百度学术
6. 刘鑫,李飒,尹福顺,姚婷. 基于动态图像技术的南海钙质土颗粒形态特征研究. 岩土工程学报. 2023(03): 590-598 . 本站查看
7. 邱筱童,尹训强,王桂萱. 考虑不同影响因素的珊瑚岛礁场地地震响应分析. 自然灾害学报. 2023(01): 131-138 . 百度学术
8. 尹福顺,李飒,刘鑫. 钙质粗粒料颗粒强度和压缩特性的试验研究. 岩土力学. 2023(04): 1120-1129+1152 . 百度学术
9. 李能,吴杨,周福霖,谭平. 岛礁吹填珊瑚砂不排水单调和循环剪切特性试验. 中国公路学报. 2023(08): 152-161 . 百度学术
10. 郭桢,蒲建,卢劲锴,黄雨. 南海西沙岛礁非饱和珊瑚砂共振柱-弯曲元试验研究. 工程地质学报. 2023(05): 1552-1562 . 百度学术
11. 吴杨,崔杰,李晨,温丽维,单振东,廖静容. 细粒含量对岛礁吹填珊瑚砂最大动剪切模量影响的试验研究. 岩石力学与工程学报. 2022(01): 205-216 . 百度学术
12. 赵云辉,孟凡超,郑志华. 相对密实度对结构性砂土动剪切模量和阻尼比影响的试验研究. 工程抗震与加固改造. 2022(01): 152-159 . 百度学术
13. 吴琪,杨铮涛,刘抗,陈国兴. 细粒含量对饱和珊瑚砂动力变形特性影响试验研究. 岩土工程学报. 2022(08): 1386-1396 . 本站查看
14. 陈龙珠,顾晓强. 共振柱试验确定土动剪切模量和阻尼比的理论辨析. 地基处理. 2022(05): 445-450 . 百度学术
15. 王伟,李犇,罗佳乐,胡俊,姜屏,李娜. 动荷载作用历史对水泥固化钙质砂三轴力学特性影响. 自然灾害学报. 2022(05): 158-167 . 百度学术
16. 周正龙,丁芷萱,刘杰,赵凯,梁珂,鹿庆蕊. 南海海域饱和粉土动剪切模量和阻尼比试验研究. 土木工程学报. 2022(S1): 227-233 . 百度学术
17. 宋前进,程磊,贺为民. 孔隙比对土体动力特性参数的影响——以豫东平原粉土为例. 科学技术与工程. 2021(07): 2830-2835 . 百度学术
18. 高盟,彭晓东,陈青生. 南海非饱和钙质砂动力特性三轴试验研究. 北京工业大学学报. 2021(06): 625-635 . 百度学术
19. ZHANG Yan-ling,DING Xuan-ming,CHEN Zhi-xiong,WU Qi,WANG Cheng-long. Seismic responses of slopes with different angles in coral sand. Journal of Mountain Science. 2021(09): 2475-2485 . 必应学术
20. 郭聚坤,王瑞,尹斌,卞贵建,魏道凯. 钢-钙质砂界面循环剪切特性试验研究. 建筑结构. 2021(S2): 1613-1617 . 百度学术
21. 宋前进,程磊,贺为民. 豫东平原粉质黏土动剪切模量与阻尼比试验研究. 地震工程学报. 2020(04): 1013-1018 . 百度学术
22. 信鹏飞,李飒,吴文娟. 冲击荷载作用下钙质砂破碎与变形特性研究. 水力发电学报. 2019(12): 102-111 . 百度学术
其他类型引用(16)