考虑孔隙流体压缩性的土体平面应变固结半解析数值分析.pdf

第 34 卷 第 1 期 岩 土 工 程 学 报 Vol.34 No.1 2012 年 .1 月 Chinese Journal of Geotechnical Engineering Jan. 2012 考虑孔隙流体压缩性的土体平面应变固结 半解析数值分析 徐 进 1，2，蔡正银1，王旭东2 1. 南京水利科学研究院岩土工程研究所，江苏 南京 210024；2. 南京工业大学岩土工程研究所，江苏 南京 210009 摘 要针对 Biot 固结基本方程的特点，在构建满足边界条件的位移和超静孔压试探函数的基础上，由 Galerkin 法建 立了考虑孔隙流体可压缩性渗透各向异性层状土体中平面应变固结问题的半解析数值求解格式，并利用三角函数系的 正交性实现了加权余量方程按不同级数项的解耦。在此基础上，编制相应计算程序实现了半解析数值方程的求解。通 过算例对比分析验证了半解析数值方法的正确性，并说明了其处理渗透各向同性、孔隙流体可压缩性和土体层状的能 力。 关键词半解析数值法；平面应变固结；可压缩性；层状地基 中图分类号TU433 文献标识码A 文章编号1000–4548201201–0089–05 作者简介 徐 进1982– ， 男， 江苏盐城人， 博士研究生， 主要从事流固耦合和环境岩土工程方面的研究工作。 E-mail jinxu1031。 Semi-analytical numerical analysis for plane strain consolidation of anisotropic soil with compressible constituents XU Jin1, 2, CAI Zheng-yin1, WANG Xu-dong2 1. Department of Geotechnical Engineering, Nanjing Hydraulic Research Institute, Nanjing 210024, China; 2. Institute of Geotechnical Engineering, Nanjing University of Technology, Nanjing 210009, China Abstract Based on the characteristics of the governing equations of the Biot’s consolidation, the suitable trial functions of displacement and excess pore water pressure satisfying the boundary conditions are obtained, and a semi-analytical numerical scheme for solving the plane strain consolidation of multi-layered and anisotropic soil with compressible constituents is presented by using the Galerkin’s . Moreover, by means of the orthogonality of the trigonometric functions series, the weighted residual equations are decoupled to several independence ones associated with different series modes. A computer program is developed for solving the semi-analytical numerical equations. Several numerical examples are presented, and the numerical results are compared with other available solutions to verify the validity of the present . The applicabilities of the present in dealing with the anisotropy of permeability, the compressibility of pore fluid and the layered characteristics of soil are demonstrated. Key words semi-analytical numerical ; plane strain consolidation; compressibility; layered soil 0 引 言 Biot 固结理论[1]反映了超静孔压消散和土骨架变 形间的耦合作用，从而能够更加准确地描述多孔介质 的力学行为。由于数学上的困难，Biot 固结方程只能 在少数简单情况下得出解析解，因此，实现该问题的 有效求解近几十年来得到了较多关注。如 Mcnamee 等[2]、Booker[3]、Ai 等[4]分别对平面应变及轴对称、 单层土体和柱坐标系下有限厚土体等多种情形下 Biot 固结问题的求解作了有意义的探索。 天然地基土大多具有层状非均质和各向异性的特 点。 艾智勇等[5-6]利用传递矩阵法对层状横观各向异性 地基 Biot 固结轴对称问题进行了研究， 运用积分变换 法分析了渗透各向异性地基的平面应变固结问题。宰 金珉等[7]提出了层状地基土的比奥固结问题的有限层 分析方法，考虑了土体的横观各向同性，实现了有限 层法在固结问题上的应用。 ─────── 基金项目江苏省六大人才高峰基金资助项目（09-F-004） 收稿日期2010–11–26 90 岩 土 工 程 学 报 2012 年 即使实际土体的饱和度高于 95，其孔隙水中往 往由于有封闭气泡的存在，而表现出一定的可压缩性 [8-10]，且研究表明孔隙流体的压缩性对土体固结行 为有较大影响[8-12]。因此，在这种情形下孔隙流体的 可压缩性不可忽略。但是，目前针对孔隙流体可压缩 性渗透各向异性层状地基固结平面应变问题的研究仍 很少,艾智勇等利用传递矩阵法和Fourier-Laplace积分 变换技术对该问题进行了研究，研究结果表明渗透 各向异性、孔隙流体可压缩性和层状特性对土体固结 行为均有较大影响[12]。 本文在构建满足边界条件的位移和超静孔压试探 函数的基础上，利用 Galerkin 法和三角函数系的正交 性建立考虑孔隙流体可压缩性和渗透各向异性的土体 Biot 平面应变固结问题的半解析数值求解格式，并根 据求解格式编制相应计算程序。通过算例对比分析验 证半解析数值法的正确性，探讨本文方法考虑孔隙流 体可压缩性、渗透各向异性和土体层状分布等因素的 计算能力。 1 基本方程和半解析数值分析模型 假设各向异性层状土体中孔隙流体可压缩，对于 平面应变问题，规定应力以压为正，则考虑孔隙流体 可压缩性的 Biot 固结方程为 2v 2v 22 v hv 22 w 0 , 12 0 , 12 1 0 0 Gp Gux zD xx Gp Gwx zD zz ppp nKKt ttxz                             ， ， 。 1 式中 u，w分别为x和z方向的位移；p为超静孔压； v 表示体积应变；G，分别为土体的剪切模量和泊 松比；为孔隙流体的压缩系数，对应于不同的饱和 度95～100，的大小在4.510 -1～4.510-4 MPa-1 间变化[10]；n为孔隙率； w 为水的重度；Kh，Kv分别 为水平向和竖向渗透系数；D表示计算区域，t为时间。 微分方程组（1）对应的上、下边界条件可根据不 同实际情形确定，一般取地表完全透水，且位移自由； 下边界取为位移固定不透水边界或位移固定透水边 界； 侧边界取为竖向约束， 水平向位移自由且透水[7, 13]。 对于图 1 所示情形，假设地表作用有任意分布的 荷载qx，对应的边界条件可表示为 0 0 0 0 0 zz xzz z q x p              ， ， ； 2 0 0 z cz c z c uw p z          ， ； 3 0 0 0 0 0 x x Ax x A x x A u w x p               ， 。 4 将计算区域D沿z方向离散成L层， 如图 1 所示， 层元l的上下节面的编号记为l和l1，且层元l可具 有独立的物理力学参数 hv ,,,, lllll GKK 和几何参数 l c，从而得到问题的半解析数值分析模型。 图 1 各向异性层状土体及边界条件示意图 Fig. 1 Schematic diagram of layered and anisotropic soil with compressible constituents 2 半解析数值求解格式 2.1 位移和超静孔压试探函数 根据计算区域的离散形式，水平节面上采用满足 边界条件（4）的解析函数表示，并利用z方向上的标 准线性形函数可得到全域上位移和超静孔压的试探函 数   1 11 LM mjmj jm fN       。 5 式中   T fu w p    ， 其中 u  ， w  和 p  分别为位移和 超静孔压试探函数； T uwp mjmjmjmj ttt  ，其 中 u mj t， w mj t和 p mj t为待求位移和孔压系数；M 为 级 数 计 算 项 数 ，L为 层 元 数 ； 形 函 数 cos 0 0 0 sin 0 0 0 sin mj mjmj mj k x Nz Nk x Nz k x Nz          ， 其中π/ m kmA， j Nz为z方向上标准线性形函数。 2.2 Galerkin法 根据Galerkin法， 取形函数作为权函数， 将式 （5） 代入基本方程（1） ，利用分部积分和边界条件（2）～ （4） ，由三角函数系的正交性[14]可实现各级数项的解 第 1 期 徐 进，等. 考虑孔隙流体压缩性的土体平面应变固结半解析数值分析 91 耦， 从而得到对应于第m项的孔隙流体可压缩土体固 结整体半解析数值方程组   Tdd 0 dd mmmmm mmmmmm KKpF KSpKp tt             ， 。 6 式（6）是关于t的常微分方程组，采用全隐式差 分格式，取时间步为t，可得 T mmmtmt mtmt mmm KKF pR KStK          。 7 式中 m K 为整体劲度矩阵，  m K 为整体耦合矩阵，  m S为整体压缩矩阵， m K   为整体渗透矩阵； mt F，  mt R分别为整体荷载向量和流量向量； mt ， mt p 分别为整体位移系数向量和孔压系数向量。 2.3 层元分析 对于某一个特定层l l1,2,,L，式（7）中相关 项只在该层域内积分计算，可得到对应于l的层元矩 阵及层元向量，分别记为 T l mm mmm KK KStK      ， l mt mt p    ， l mt mt F R    。其中， l m K  ， l m K ，和 l m K   的表达式与 文献[7]中给出的有限层求解格式具有相同形式，只需 将原三维问题退化到平面应变情形；层元位移和孔压 系数分别为 T 1 1 , ll uwuw mtmlmlm lm lmt p  T 1 pp mlm l     ； 层元荷载向量 l mt F 11 xzxz mlmlm lm l FFFF   ； 二阶的层元压缩矩阵 l m S  T d lmlml D NnN，其中 ml N sin ml k x N z  1 sin ml k x Nz  ； 二维的层元流量向量 T ll mtm RK  ll l m ttmm tt Sp    。 为了方便地通过层元矩阵叠加得到整体矩阵，将 式 （7） 中同一节面对应的线性方程调整在一起重新排 序[7]，则层元矩阵和层元向量也作相应调整，如调整 后的待求系数向量为  T 111 l uwpuwp mjmlmlmlm lm lm l     。8 调整后的层元矩阵仍为对称矩阵，其与荷载–流 量向量可一起展开为 11 2122 313233 141424344 5152535455 1 616263636566 x ml z ml ml x m l z m l F k F kk R kkk Fkkkk kkkkk F kkkkkk             ， 1m l R             。 式中， 22 114414 /3/ 2 l ml Ac kkd kd c ， 2124 /2 2 l ml Ac kddkc ， 2 22554 /3 2 l m Ac kkd k 2 3/ l d c ， 31 /3 2 l m Ac kk ， 32 4 A k ， 3366 kk 22 hmv w 3/ 36 ll l ActAc K kK cn    ， 2 411m /6 2 l Ac kd k 2 4/ l d c ， 4224m /2 2 l l Ac kddkc ， 43 k m 12 l Ac k； 5124m /2 2 l l Ac kdd kc  ， 53 4 A k  ， 52 2 l Ac k 22 4m3 /6/ l d kd c ， 61m 12 l Ac kk ， 62 4 A k ， 64m 6 l Ac kk ， 65 4 A k  ， 5424m /2 2 l l Ac kdd kc  ， 63 w 6 l Act k     22 hmv 6/ 12 l l Ac K kK cn， uw mmm m 232 lll l k cA R       uwp p mm1m1m1 m1 62236 llll lml l tt k c Ac nR           ， uw uwp mm1m1 mmmm 2623226 lll lllll k c k cAcA n          p m1 3 l tt       。其中， 2 1 1d， 2 1d， 2 3 1d， 4 dG为弹性参数，/ 1E 2 12  ，2 1EG。 根据不同荷载情形，可相应地积分计算得到层元 荷载向量中的元素， 如对于作用于地表宽为a1的均布 荷 载q的 情 形 可 得 1 m1m m cos 22 z aqA Fk k      1 m cos 22 aA k       ，其余项均为零。 2.4 整体分析 在得到各层元矩阵和层元列阵后，按常规方式可 叠加得到整体矩阵[13]，建立对应于第m项的3L1 阶线性方程组，且其半带宽仅为6。 在式（5）中已经考虑了侧边界上边界条件，对于 式（2）～（4）中上、下边界处位移和孔压为已知的 边界条件可参考文献[15]进行处理，特别地，对于边 界值均为零的情形，只需对与已知边界值无关的位移 和孔压系数建立代数方程式即可。 3 算例分析 为了验证半解析数值方法和计算程序的正确性， 说明本文方法考虑渗透各向同性、孔隙流体可压缩性 和土体层状分布等因素的计算能力，进行了算例对比 分析。 92 岩 土 工 程 学 报 2012 年 3.1 均质各向同性 利用本文方法对条形荷载作用下有限厚度均质各 向同性土体固结问题进行了计算， 并与Booker的经典 解答进行了对比[3, 12]，计算参数和计算结果见图2。 由图可见， 半解析数值解与Booker的解答之间吻合良 好，验证了本文理论和计算程序的正确性。 3.2 渗透各向异性 对3种不同渗透系数各向异性参数 hv KK 情形 下固结问题进行了计算分析。图3给出了荷载作用中 心处地表沉降随时间的计算结果，半解析数值结果与 Ai等[6]的解答吻合良好。同时计算结果也表明了渗透 各向异性对土体固结行为的较大影响竖向渗透系数 一定时，渗透各向异性越明显，土体的固结越快，而 地表的最终沉降保持一致。 图 2 各向同性均质土体固结度 U 随时间变化曲线 Fig. 2 Variation of consolidation degree U of isotropic and ..homogenous soil 图 3 不同渗透各向异性下地表沉降曲线 Fig. 3 Curves of surface settlement under different anisotropic permeabilities 3.3 孔隙流体可压缩性 为了验证半解析数值方法考虑孔隙流体可压缩性 的能力，引入参数1/ 1n MnMG，其 中2 1/ 12MG， 随增大， 孔隙流体压缩性 越大。对1.0，1.05，1.15和1.25四种 不同压缩性参数取值情形进行了分析，其中泊松比 取0.35， hv KK 1。计算结果见图4，5，由图可见， 在其他条件相同时，压缩性参数越大，地表初始沉 降越大，孔压消散越快，但是地表最终沉降相同[12]。 3.4 多层各向异性孔隙流体可压缩土体 最后选取了一个5层非均质地基算例[12]来说明半 解析数值方法求解层状渗透各向异性地基中孔隙流体 可压缩土体固结问题的能力，计算参数和分析结果见 图6。本文方法的解耦性使得编制程序中需求解的线 性方程组规模较小， 同时即使在计算项数M仅为8的 情形下，其计算结果已经与传递矩阵解保持了较好的 一致性，体现了解析函数快速收敛的特点，这些都说 明了半解析数值方法较高的计算效率。 图 4 不同压缩性参数下地表沉降曲线 Fig. 4 Curves of surface settlement under different compressibility parameters 图 5 不同压缩性参数下孔压随深度分布 Fig. 5 Distribution of excess pore water pressure under different .compressibility parameters 图 6 多层土体地表沉降曲线 Fig. 6 Curves of surface settlement of layered soil 第 1 期 徐 进，等. 考虑孔隙流体压缩性的土体平面应变固结半解析数值分析 93 4 结 论 （1）基于Galerkin法和三角函数系正交性建立 了考虑孔隙流体可压缩性各向异性层状土体Biot平面 应变固结问题的半解析数值求解格式。由于层元矩阵 间不相耦合， 减少了整体矩阵带宽， 提高了计算效率， 便于编程实现和实际应用，体现了半解析数值方法求 解层状介质问题的优势。同时应该指出的是，在处理 具有不规则边界以及复杂非均质问题时，本文方法无 法取代有限元等传统数值方法的地位。 （2）通过算例对比分析验证了半解析数值法的 正确性，说明了其处理渗透各向异性、孔隙流体可压 缩性和土体层状的能力，计算结果很好地反映了渗透 各向异性和孔隙流体可压缩性对土体固结行为的影 响。 参考文献 [1] BIOT M A. General theory of three-dimensional consolidation[J]. Journal of Applied Physics, 1941, 122 155–164. [2] MCNAMEE J, GIBSON R E. Plane strain and axially symmetric problem of the consolidation of a semi-infinite clay stratum[J]. Quarterly Journal of Mechanics and Applied Mathematics, 1960, 132 210–227. [3] BOOKER J R. The consolidation of a finite layer subject to surface loading[J]. International Journal of Solids and Structure, 1974, 109 1053–1065. [4] AI Z Y, WANG Q S, WU C. A new for solving Biot’s consolidation of a finite soil layer in the cylindrical coordinate system[J]. Acta Mechanica Sinica, 2008, 246 691–697. [5] 艾智勇, 王全胜. 多层横观各向同性地基轴对称固结的传 递 矩 阵 解 [J]. 岩 土 力 学 , 2009, 304921 – 925. AI Zhi-yong, WANG Quan-sheng. Transfer matrix solutions for axisymmetric consolidation of multilayered transversely isotropic soils[J]. Rock and Soil Mechanics, 2009, 304921 –925. in Chinese [6] AI Z Y, WU C. Plane strain consolidation of soil layer with anisotropic permeability[J]. Applied Mathematics and Mechanics, 2009, 3011 1437–1444. [7] 宰金珉, 梅国雄. 有限层法求解三维比奥固结问题[J]. 岩 土工程学报, 2002, 241 31–33. ZAI Jin-min, MEI Guo-xiong. Finite layer analysis of three dimensional Biot consolidation[J]. Chinese Journal of Geotechnical Engineering, 2002, 241 31–33. in Chinese [8] YUE Z Q, SELVADURAI A P S, LAW K T. Excess pore water pressure in a porelastic seabed saturated with a compressible fluid[J]. Canadian Geotechnical Journal, 1994, 31989– 1003. [9] CHEN G J. Axi-symmetric solutions for layered and poro-elastic media with compressible constituents[C]// Proceeding of the 3rd International Conference on Soft Soil Engineering, Hong Kong, 2001 251–256. [10] CHEN G J. Consolidation of multilayered half space with anisotropic permeability and compressible constituents[J]. International Journal of Solids and Structures, 2004, 4116-17 4567–4586. [11] SENJUNTICHAI T, RAJAPAKSE R K N D. Exact stiffness for quasi-statics of a multi-layered poroelastic medium[J]. International Journal of Solids and Structures, 1995, 32111535–1553. [12] 艾智勇, 吴 超. 渗透各向异性可压缩地基固结的平面应 变分析[J]. 力学学报, 2009, 415 801–806. AI Zhi-yong, WU Chao. Analysis on plane strain consolidation of a multi-layered soil with anisotropic permeability and compressibility constituents[J]. Chines Journal of Theoretical and Applied Mechanics, 2009, 415 801 – 806. in Chinese [13] 张问清, 赵锡宏, 宰金珉. 任意力系作用下的层状弹性半 空间的有限层分析方法[J]. 岩土工程学报, 1983, 32 27 –42. ZHANG Wen-qing, ZHAO Xi-hong, ZAI Jin-min. Finite layer analysis for layered and elastic half-space under arbitrary force systems[J]. Chinese Journal of Geotechnical Engineering, 1983, 32 27–42. in Chinese [14] HARRY F. Davis. Fourier series and orthogonal functions[M]. Boston Allyn and Bacon INC, 1963. [15] 朱百里, 沈珠江, 等. 计算土力学[M]. 上海 上海科学技 术出版社, 1990. ZHU Bai-li, SHEN Zhu-jiang, et al. Computational soil mechanics[M]. Shanghai Shanghai Science and Technology Press, 1990. in Chinese