Loading [MathJax]/extensions/TeX/boldsymbol.js

    热传导CPT探头的研发与应用

    刘松玉, 郭易木, 张国柱, 周遊

    刘松玉, 郭易木, 张国柱, 周遊. 热传导CPT探头的研发与应用[J]. 岩土工程学报, 2020, 42(2): 354-361. DOI: 10.11779/CJGE202002017
    引用本文: 刘松玉, 郭易木, 张国柱, 周遊. 热传导CPT探头的研发与应用[J]. 岩土工程学报, 2020, 42(2): 354-361. DOI: 10.11779/CJGE202002017
    LIU Song-yu, GUO Yi-mu, ZHANG Guo-zhu, ZHOU You. Development and application of heat conduction CPT probe[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(2): 354-361. DOI: 10.11779/CJGE202002017
    Citation: LIU Song-yu, GUO Yi-mu, ZHANG Guo-zhu, ZHOU You. Development and application of heat conduction CPT probe[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(2): 354-361. DOI: 10.11779/CJGE202002017

    热传导CPT探头的研发与应用  English Version

    基金项目: 

    国家自然科学基金项目 51578146

    国家自然科学基金项目 51778138

    江苏省研究生科研与实践创新计划项目 KYCX18_0106

    详细信息
      作者简介:

      刘松玉(1963— ),男,博士,教授,博士生导师,主要从事现代原位测试技术和特殊地基处理技术等方面的研究工作。E-mail:liusy@seu.edu.cn

    • 中图分类号: TU473.1

    Development and application of heat conduction CPT probe

    • 摘要: 土体热导率是能源地下工程、高压电线埋设、冻土路基等工程设计中的重要参数,是评估土体换热性能及地下温度场计算的关键,但目前缺乏有效的原位测试手段。基于瞬时线热源温度消散理论,研发了可测试原位土体热导率的静力触探(CPT)探头。根据理论假设与现有CPT系统尺寸,确定了探头具体长度、直径、内部构造、温度采集点位,并提出相应测试步骤与热导率计算方法。利用COMSOL有限元分析软件对测试过程进行模拟验证,结果表明探头实际传热符合线热源假定,且计算方法适用于一般热导率土体。对于热导率较小(小于0.6 W/(m·K))土体,需适当延长测试时间。现场应用表明,土体原位热导率略高于取样土室内测试结果,表明取样扰动可能降低土体导热性能,最后对取样测试及工程设计提出改进建议。
      Abstract: The thermal conductivity is the key parameter to the design of many projects, such as energy structures, high-voltage buried power cables and permafrost embankment, related to estimating the heat transfer capability and temperature field in the soil. However, at present there is no effective in-situ testing method. Based on the theory of instantaneous heat release along a line source, a heat conduction cone penetration test (CPT) probe for thermal conductivity evaluation of in-situ soil is developed. According to the theoretical assumptions and the sizes of CPT system, the length, diameter, internal structure and positions of the temperature sensors are introduced. Then, the corresponding test procedure and the method for thermal conductivity are proposed. The test process is simulated in COMSOL to verify the method, and the results validate that the actual heat transfer conforms to the line source theory. The interpretation method yields reasonable values within a general range of conductivities. For less conductive soil (<0.6 W/ (m·K)), longer duration of heat dissipation may be required. The field test results show that the in-situ soil conductivity is higher than that from laboratory tests on undisturbed samples, indicating the sampling disturbance may be responsible for this reduction. Finally, some suggestions on laboratory thermal conductivity tests and engineering designs are given.
    • 渗流计算是评价水工建筑物渗流安全及进行渗控方案设计的依据。渗流计算的实践表明有限元方法是一种非常成熟、非常有效的数值计算方法。有限元方法分为协调有限元方法和非协调有限元方法,要求近似函数在剖分单元边界上连续的称为协调有限元方法,允许近似多项式函数在单元边界上不连续的称为非协调有限元方法,又称为间断有限元方法。通常讲的有限元方法是指协调有限元方法,称为一般有限元方法。一般有限元方法连续性的要求,使得剖分网格只能使用单一的单元类型,同一次数的近似多项式函数,因而一般有限元方法不适用于混合剖分网格。在渗流计算中,往往希望通过加密网格的方法来提高计算精度,而这有可能导致计算格式不稳定、计算结果不收敛。国内外众多学者一直在探索适用于混合剖分网格的间断有限元方法,主要成果有:内罚间断有限元法(IPDG)[1],局部间断有限元法(LDG)[2]和可杂交的间断有限元法(HDG)[3]。与一般有限元法比较,非协调有限元方法有两方面的优势:①近似函数灵活,可以在不同剖分单元采用不同次数的近似多项式函数;②网格生成灵活,适用混合剖分网格。但是,间断有限元方法的变分形式比协调有限元方法的变分形式要复杂得多,计算格式的稳定性往往需要靠调整参数来保证、而辅助变量的引入和变分的复杂化则使得单元矩阵非局部化,这些困难使得间断有限元方法难以在渗流分析中充分发挥其优势[4]。弱有限元法(WG)是Wang等[5]提出的,它属于非协调有限元方法,通过在剖分单元上定义弱函数和弱梯度,从而用弱梯度算子替代经典的梯度算子,并在计算格式中引入了一个不含有参数的稳定子或光滑子,保证了格式的绝对稳定性[5-7]。因此,弱有限元方法不但具有间断有限元方法的优点,还克服了间断有限元方法的缺点,目前,弱有限元方法可以用来求解Stokes方程、Helmholtz方程、Maxwell方程、线弹性方程等[8-12]

      本文主要研究弱有限元法针对混合剖分网格在渗流分析上的具体应用。安排如下:在第一节中将介绍弱有限元方法的原理,并通过一个简例介绍弱有限元方法处理混合网格的流程;在第二节中运用弱有限元方法,采用混合剖分网格求解均值土坝渗流自由面和闸基渗流场。

      稳定渗流的数学模型为椭圆型方程边值问题:

      {(pu)=fin Ωu=con Γ1un=gon Γ2 (1)

      式中:Γ1Γ12=ΩnΓ2上的外法线单位向量;u为水头函数,p为渗透系数矩阵。此问题的一般有限元法的变分形式为:求uH1(Ω),使得在Γ1u=c,且满足

      a(u,v)=(f,v)+g,vΓ2 vH10(Ω)

      式中,a(u,v)=(pu,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,ττnT (2)

      单元T上的弱有限元空间定义为

      Wj,l(T)={v={v0,vb}|v0(T)Pj(T),vb(e)Pl(e)eT}

      弱梯度算子是由Wj,l(T)[Pm(T)]2上的一个映射。剖分Th上的弱有限元法空间定义为

      Vh={v|v(T)Wj,l(T),TTh}

      V0h={v={v0,vb}Vh,vb|Γ1=0,}

      问题(1)的弱有限元法的变分形式为,求uhVh,使得在Γ1uh=c,且满足

      ah(uh,v)+s(uh,v)=(f,v)+g,vΓ2, vV0h (3)

      其中ah(uh,v)=(pwuh,wv)s(uh,v)是为了保证格式稳定的一个补充项,称为稳定子或光滑子,且

      s(uh,v)=TTh(J(uh),J(v))T (4)

      式中,J(v)=h1/2T(Q(v0)vb)hT表示单元T的直径,即T上任意两点间距离的最大值(编程计算时一般取为多边形单元任意两节点间距离的最大值),Q(v0)表示内部函数v0在边界上的投影。Θ(u0)表示u0在单元边界T上的投影。根据弱梯度的定义和分部积分公式计算得

      (wuh,wv)T=(u0,v0)TΘ(u0)n,Q(v0)vbTQ(u0)ub,wvnT

      注意到

      Θ(u0)n(Q(u0)ub)/hT, wvn(Q(v0)vb)/hT

      (u0,v0)T(wuh,wv)T+2h1TQ(u0)ub,Q(v0)vbT

      因此在变分式中用ah(uh,v)近似a(u,v),把s(uh,v)作为补偿项加入变分式中是合理的[7]

      弱有限元空间Wj,l(T)表示在单元T内部近似函数采用j次多项式,在单元边界T上近似函数采用l次多项式。[Pm(T)]2表示弱函数的弱梯度是分量为m次多项式的向量。选取不同的j,lm,会得到单元局部矩阵不同的计算格式,Wj,l(T)[Pm(T)]2简记为Wj,l(T)/[Pm(T)]2,下面假设单元T是一个三角形单元(图 1)。

      图  1  三角形单元图
      Figure  1.  Triangular element

      若选取弱有限元空间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)}=6i=1aTiφi

      变分式(3)右端项,

      (f,v)T=(f,v0)T (5)
      g,vTΓ2=eTΓ2g,vbe (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号单元是个四边形,一条边中间有个节点,那么这条边上的近似函数可以理解为分段近似,这在局部矩阵计算格式上相当于两条边,这个点把该边一分为二,这个单元就可以被认定为五边形。类推,如果一条边中间有多个节点,就可把该条边分成多条边,该单元就被认定为相应类型的多边形单元。因此,对于弱有限元方法,网格剖分比较自由。

      图  2  剖分网格图
      Figure  2.  Diagram of partial grids

      下面通过一个简单的例子来了解弱有限元法处理混合网格的计算过程和计算效果。

      实例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
      下载: 导出CSV 
      | 显示表格
      表  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
      下载: 导出CSV 
      | 显示表格

      每个单元的弱有限元空间取 {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)。

      图  3  单元边界编号
      Figure  3.  Number of element boundary

      如果把所有单元内部近似函数的自由度按单元序号排前面,即,第 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),把所有局部单元矩阵组合成整体系数矩阵AA是一个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 。这说明弱有限元法处理混合剖分网格是可行,且有效的。

      在岩土边坡、土坝、地下洞室及地下水运动等渗流分析中,均存在有渗流自由面问题。渗流自由面的确定是一个重点,也是个难点。求解渗流自由面的方法有虚单元法[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} ,需做特殊处理,该单元边界上对应的基函数应取为

      图  4  坝体剖分网格Ⅰ
      Figure  4.  Diagram of dam body meshing Ⅰ
      \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 { 。 }

      所以单元局部矩阵中的MbbST部分的计算公式需做相应调整。比如,假设三角形单元(图 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
      下载: 导出CSV 
      | 显示表格

      确定溢出点的方法:根据计算结果得到各节点的水头值后,搜索在坝体下游坝面( 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列。

      图  5  自由面
      Figure  5.  Free surface

      结果对比显示,在相同剖分网格下,弱有限元法的最大绝对误差仅为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

      图  6  坝体剖分网格Ⅱ
      Figure  6.  Diagram of dam body meshing Ⅱ
      表  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
      下载: 导出CSV 
      | 显示表格

      总体方程组的自由度大大地减少了,而数值解的精度却没有降低,最大绝对误差仅为0.05,最大相对误差为0.0095。这说明,弱有限元法能处理混合网格这一优势是非常显著的。

      弱有限元法适用于混合剖分网格,那么在进行复杂的渗流场分析时,可以根据需要采用有针对性的剖分网格,这样可以用较少的自由度得到精度较高的数值计算结果。下面用弱有限元法来分析有防渗帷幕的闸基渗流场。

      计算模型2:模型长和高分别为60,20 m,闸基底板长20 m,上游20 m,下游20 m,上游水位深4 m,下游水位默认为0 m,闸基底部默认为不透水边界,在闸基底部靠上游端下方有1 m厚防渗墙(见图 7)。

      图  7  闸基模型
      Figure  7.  Model for brake base

      闸基渗流分析重点关心防渗墙下端处和闸基底板下游出口处的水力梯度,因此采用有针对性的混合剖分网格(图 8)。

      图  8  闸基剖分网格图
      Figure  8.  Meshes of gate base

      防渗墙周围和闸基底板下游出口处两个部位采用细网格的四边形单元,其它部位采用大网格四边形单元,大网格与细网格相邻的大网格被认定为八边形单元,因此,单元类型共两种:四边形和八边形,剖分单元一共有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

      图  9  防渗墙深5 m时水头等值线分布图
      Figure  9.  Distribution of water head contour line at depth of anti-seepage wall of 5 m
      图  10  防渗墙深10 m时水头等值线分布图
      Figure  10.  Distribution of water head contour line at depth of anti-seepage wall of 10 m
      图  11  防渗墙深15 m时水头等值线分布图
      Figure  11.  Distribution of water head contour line at depth of anti-seepage wall of 15 m
      表  5  水力梯度值
      Table  5.  Values of hydraulic gradient
      防渗墙
      深/m
      防渗墙底端处水力梯度 闸基出口处
      水力梯度
      5 0.323 0.332
      10 0.333 0.262
      15 0.358 0.190
      下载: 导出CSV 
      | 显示表格

      计算结果表明随着防渗墙的加深下游水头值逐渐降低,防渗墙的防渗效果得到充分体现,也说明了采用有针对性的混合剖分网格,运用弱有限元法来处理类似闸基这种比较复杂的渗流场问题是非常有效的。

      计算渗流自由面模型的数值结果表明弱有限元法的计算结果的精度比一般有限元法的要高,并且采用较稀疏的混合剖分网格也能得到精度很好的数值结果。利用弱有限元法适用混合剖分网格的特点,对较复杂的闸基渗流场进行分析,根据数值结果绘制的等水头线表明结果符合闸基渗流的特征。高精度和能适用于混合剖分网格这两个优势,使得弱有限元法在渗流分析中有很好的应用前景。

    • 图  1   热传导CPT探头基本结构及测试原理

      Figure  1.   Basic structure of heat conduction CPT probe and principle of tests

      图  2   原位热物性测试数值模型

      Figure  2.   Numerical model for in-situ thermal conductivity tests

      图  3   典型测试温度场与热流场分布云图(土体热导率1.8 W/(m·K))

      Figure  3.   Distribution of temperature fields and heat flows of a typical test in 1.8 W/(m·K) soil

      图  4   探头各边界热流分布变化

      Figure  4.   Variation of heat flows though each boundary of probe

      图  5   拟合温度消散与初始地温

      Figure  5.   Fitted temperature dissipation and initial ground temperature

      图  6   土体热导率计算

      Figure  6.   Calculated results of thermal conductivity of soil

      图  7   温度消散2000 s曲线与热导率计算(真实值0.6 W/(m·K))

      Figure  7.   Temperature dissipation for 2000 s and calculated thermal conductivity of soil (true value 0.6 W/(m·K))

      图  8   热传导CPT测试系统与现场测试

      Figure  8.   Heat conduction CPT probe system and field application

      图  9   现场测试与拟合温度消散曲线

      Figure  9.   In-situ test and fitted temperature dissipation

      图  10   现场土体热导率计算

      Figure  10.   Calculation of thermal conductivity of in-situ soil

      表  1   模型热物性参数表

      Table  1   Summary of thermal properties of numerical model

      材料热导率/(W·m-1·K-1)比热容/(J·kg-1·K)密度/(kg·m-3)
      空气0.02610151.16
      加热片24.28692548
      回填材料0.2310502200
      金属外壳469002700
      土体9002000
      注:研究对土体进行参数化扫描研究,每组热导率间隔为0.3 W/(m·K),范围为0.6~2.7 W/(m·K)。
      下载: 导出CSV

      表  2   参数化扫描结果

      Table  2   Summary of parametic sweep results

      土体热导率 /(W·m -1·K -1)拟合消散热导率/(W·m -1·K -1)误差/%原消散热导率/(W·m -1·K-1)误差/%
      0.60.85420.8948
      0.91.00111.0213
      1.21.2111.200
      1.51.46-21.37-9
      1.81.73-41.53-15
      2.12.05-21.70-19
      2.42.37-11.85-23
      2.72.7932.00-26
      注:误差为(拟合值-真实值)/真实值,+表示高估,-表示低估。
      下载: 导出CSV
    • [1]

      BRANDON T, MITCHELL J, CAMERON J. Thermal instability in buried cable backfills[J]. Journal of Geotechnical Engineering, 1989, 115(1): 38-55. doi: 10.1061/(ASCE)0733-9410(1989)115:1(38)

      [2]

      SLEGEL D L, DAVIS L. Transient heat and mass transfer in soils in the vicinity of heated porous pipes[J]. Journal of Heat Transfer, 1977, 99(4): 541-546. doi: 10.1115/1.3450739

      [3] 白冰, 赵成刚. 温度对黏性土介质力学特性的影响[J]. 岩土力学, 2003, 24(4): 533-537. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200304012.htm

      BAI Bing, ZHAO Cheng-gang. Temperature effects on mechanical characteristics of clay soils[J]. Rock and Soil Mechanics, 2003, 24(4): 533-537. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200304012.htm

      [4]

      LACKNER R, AMON A, LAGGER H. Artificial ground freezing of fully saturated soil: thermal problem[J]. Journal of Engineering Mechanics, 2005, 131(2): 211-220. doi: 10.1061/(ASCE)0733-9399(2005)131:2(211)

      [5] 商允虎, 牛富俊, 刘明浩, 等. 多年冻土区桥梁工程桩基础服役期温度场研究[J]. 岩石力学与工程学报, 2017, 36(9): 2313-2323. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201709026.htm

      SHANG Yun-hui, NIU Fu-jun, LIU Ming-hao, et al. Long-term effect of a pile foundation on ground temperatures in permafrost regions[J]. Chinese Journal of Rock Mechanics and Engineering, 2017, 36(9): 2313-2323. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201709026.htm

      [6]

      BRANDL H. Energy foundations and other thermo-active ground structures[J]. Géotechnique, 2006, 56(2): 81-122. doi: 10.1680/geot.2006.56.2.81

      [7]

      FAROUKI O T. Thermal Properties of Soils, CRREL Monograph 81-1[R]. New Hampshire: U. S. Army Cold Regions Research and Engineering Laboratory Hanover, 1981.

      [8]

      LOW J E, LOVERIDGE F A, POWRIE W, et al. A comparison of laboratory and in situ methods to determine soil thermal conductivity for energy foundations and other ground heat exchanger applications[J]. Acta Geotechnica, 2015, 10(2): 209-218. doi: 10.1007/s11440-014-0333-0

      [9] 郭志光, 白冰. 描述饱和土热固结过程的一个非线性模型及数值分析[J]. 岩土工程学报, 2018, 40(11): 2061-2067. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201811014.htm

      GUO Zhi-guang, BAI Bing. Nonlinear model and numerical simulation of thermal consolidation process of saturated soils[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(11): 2061-2067. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201811014.htm

      [10]

      BRANDON T, MITCHELL J. Factors influencing thermal resistivity of sands[J]. Journal of Geotechnical Engineering, 1989, 115(12): 1683-1698. doi: 10.1061/(ASCE)0733-9410(1989)115:12(1683)

      [11]

      ASTM. D5334-14:Standard Test Method for Determination of Thermal Conductivity of Soil and Soft Rock by Thermal Needle Probe Procedure[S]. West Conshohocken, PA, USA; ASTM International. 2014.

      [12]

      IEEE Std 442-1981 IEEE Guide for Soil Thermal Resistivity Measurements[S]. New York; Institute of Electrical and Electronics Engineers, 1981.

      [13]

      GANGADHARA RAO M, SINGH D. A generalized relationship to estimate thermal resistivity of soils[J]. Canadian Geotechnical Journal, 1999, 36(4): 767-773. doi: 10.1139/t99-037

      [14]

      ABU-HAMDEH N H, REEDER R C. Soil thermal conductivity effects of density, moisture, salt concentration, and organic matter[J]. Soil Science Society of America Journal, 2000, 64: 1285-1290. doi: 10.2136/sssaj2000.6441285x

      [15]

      SINGH D N, DEVID K. Generalized relationships for estimating soil thermal resistivity[J]. Experimental Thermal and Fluid Science, 2000, 22(3): 133-143.

      [16] 肖琳, 李晓昭, 赵晓豹, 等. 含水量与孔隙率对土体热导率影响的室内实验[J]. 解放军理工大学学报(自然科学版), 2008, 9(3): 241-247. https://www.cnki.com.cn/Article/CJFDTOTAL-JFJL200803009.htm

      XIAO Lin, LI Xiao-zhao, ZHAO Xiao-bao, et al. Laroratory on influences of moisture content and porosity on thermal conductivity of soils[J]. Journal of PLA University of Science and Technology, 2008, 9(3): 241-247. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-JFJL200803009.htm

      [17]

      BARRY-MACAULAY D, BOUAZZA A, SINGH R M, et al. Thermal conductivity of soils and rocks from the Melbourne (Australia) region[J]. Engineering Geology, 2013, 164: 131-138. doi: 10.1016/j.enggeo.2013.06.014

      [18]

      ABUEL-NAGA H, BERGADO D, BOUAZZA A. Thermal conductivity evolution of saturated clay under consolidation process[J]. International Journal of Geomechanics, 2008, 8(2): 114-122. doi: 10.1061/(ASCE)1532-3641(2008)8:2(114)

      [19]

      ASHRAE, Methods for determining soil and rock formation thermal properties from short-term field tests, ASHRAE Research Summary 1118-TRP[R]. American Society of Heating, Refrigerating and Air-Conditioning Engineers, 2002.

      [20]

      SANNER B, HELLSTR M G, SPITLER J, et al. Thermal response test-current status and world-wide application[C]//Proceedings of the World Geothermal Congress. Antalya, Turkey, 2005: 1-9.

      [21]

      ZHANG C, GUO Z, LIU Y, et al. A review on thermal response test of ground-coupled heat pump systems[J]. Renewable and Sustainable Energy Reviews, 2014, 40: 851-867. doi: 10.1016/j.rser.2014.08.018

      [22]

      DENG Y, FEDLER C. Multi-layered soil effects on vertical ground-coupled heat pump design[J]. Transactions of the ASAE, 1992, 35(2): 687-694.

      [23]

      FUJII H, OKUBO H, NISHI K, et al. An improved thermal response test for U-tube ground heat exchanger based on optical fiber thermometers[J]. Geothermics, 2009, 38(4): 399-406.

      [24]

      EWEN J, THOMAS H R. The thermal probe—measurement of the thermal conductivity and drying rate of soil in the field[J]. Geotechnical Testing Journal, 1992, 15(3): 256-263.

      [25]

      AKROUCH G A, BRIAUD J-L, SANCHEZ M, et al. Thermal cone test to determine soil thermal properties[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2015, 142(3): 04015085.

      [26]

      BLACKWELL J. The axial-flow error in the thermal- conductivity probe[J]. Canadian Journal of Physics, 1956, 34(4): 412-417.

      [27]

      Hukseflux Inc. Hukseflux Thermal Sensors. MTN01 Manual[M]. Cersion 1008. The Netherlands, Delft: Hukseflux Inc., 2003.

      [28]

      Decagon Devices, Inc. KD2 Pro Thermal Properties Analyzer Operator's Manual[M]. Version 5. Pullman, Washington: Decagon Devices, Inc., 2008.

      [29]

      CARSLAW H S, JAEGER J C. Conduction of Heat in Solids[M]. 2nd ed. Oxford: Clarendon Press, 1959.

    • 期刊类型引用(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)

    图(10)  /  表(2)
    计量
    • 文章访问数:  370
    • HTML全文浏览量:  29
    • PDF下载量:  257
    • 被引次数: 38
    出版历程
    • 收稿日期:  2019-03-10
    • 网络出版日期:  2022-12-07
    • 刊出日期:  2020-01-31

    目录

    /

    返回文章
    返回