摘要： 凍土在我國(guó)分布十分廣泛，凍脹現(xiàn)象是凍土區(qū)經(jīng)常出現(xiàn)的問(wèn)題，其導(dǎo)致的工程病害屢見(jiàn)不鮮。凍脹主要是由土體內(nèi)部溫度變化及水分遷移造成的，是一類(lèi)極其復(fù)雜的溫度、滲流及應(yīng)力多場(chǎng)耦合問(wèn)題。通過(guò)對(duì)凍土中的溫度場(chǎng)進(jìn)行分析，求解含冰量及孔隙率隨溫度變化條件下的凍土一般性瞬態(tài)熱傳導(dǎo)方程，獲得其近似解析解并應(yīng)用于凍土熱傳導(dǎo)過(guò)程的預(yù)測(cè)。通過(guò)與有限元計(jì)算所得精確解對(duì)比，驗(yàn)證其有效性。研究成果對(duì)凍土凍脹問(wèn)題的深入研究具有一定的理論和實(shí)踐意義。

關(guān)鍵詞： 凍土; 溫度場(chǎng)方程; 含冰量; 孔隙率; 溫度變化; 解析解

中圖分類(lèi)號(hào)： TU44 文獻(xiàn)標(biāo)志碼：A 文章編號(hào)： 1000-0844（2024）06-1373-07

DOI：10.20000/j.1000-0844.20230106001

Analytical solution for the temperature

field equation of frozen soil

WANG Dongyuan¹， ZHANG Ke¹， CHEN Xi²， WU Xiaopeng²

（1. China Petroleum Pipeline Engineering Corporation， Langfang 065000， Hebei， China;

2. College of Civil Engineering and Architecture， Jiaxing University， Jiaxing 314001， Jiangsu， China）

Abstract： Frozen soil is widely distributed in China. As a common problem in frozen soil areas， frost heave frequently causes engineering diseases. Frost heave， which is mainly formed by the temperature change and water migration in the soil during freezing and thawing， is an extremely complex multifield coupling problem of temperature， seepage， and stress. In this paper， the general transient heat conduction equation of frozen soil under the changes of ice content and porosity with temperature was solved by analyzing the temperature field in frozen soil， and its approximate analytical solution was obtained and applied to predict the heat conduction process of frozen soil. Compared with the exact solution obtained from finite element calculation， the effectiveness of the analytical solution was verified. The research results show certain theoretical and practical significance for the intensive study of the problem of frost heave in frozen soil.

Keywords： frozen soil; temperature field equation; ice content; porosity; temperature change; analytical solution

0 引言

我國(guó)是一個(gè)凍土大國(guó)，多年凍土主要分布在高緯度的東北大小興安嶺、長(zhǎng)白山，以及西部的青藏高原、祁連山、天山和阿爾泰山等高山、高原地區(qū)，多年凍土區(qū)約占中國(guó)國(guó)土面積的21.5%；季節(jié)凍土主要分布在黑龍江、吉林、遼寧、內(nèi)蒙古、甘肅、寧夏、青海、新疆北部等地，面積約占中國(guó)國(guó)土面積的53.5%。多年凍土和季節(jié)性?xún)鐾撩娣e總和超過(guò)國(guó)土總面積的2/3^［1］。在凍土區(qū)，土體會(huì)發(fā)生凍脹現(xiàn)象^［2］，導(dǎo)致橋梁、道路、管線、水渠等基礎(chǔ)設(shè)施頻受凍害困擾。

土體凍脹通常會(huì)使得橋梁樁基、路基邊坡，油氣管道及水渠發(fā)生較大的凍脹變形，導(dǎo)致各種開(kāi)裂和破壞。此外，凍脹還會(huì)引起凍脹丘、冰錐^［3］，并造成凍土區(qū)鐵軌大變形、基坑支護(hù)結(jié)構(gòu)破壞、隧道襯砌掛冰、排水系統(tǒng)癱瘓等病害^{［2，4-5］}。大量事實(shí)表明^［2］，凍土區(qū)發(fā)生的土體凍脹現(xiàn)象會(huì)造成各類(lèi)公用及民用設(shè)施破壞，影響其正常運(yùn)營(yíng)，嚴(yán)重時(shí)甚至?xí)绊懙饺嗣竦墓藏?cái)產(chǎn)及生命安全。針對(duì)凍土凍脹問(wèn)題，國(guó)內(nèi)外學(xué)者開(kāi)展了大量的研究工作，相關(guān)資料顯示^{［2，6-7］}，土體的凍脹是指土壤中溫度降低，孔隙水凝結(jié)成冰引起的土體體積膨脹現(xiàn)象。受控于溫度狀況，土體中的水分會(huì)向凍結(jié)鋒面遷移、聚集并結(jié)晶形成冰透鏡，是土體凍脹量的主要來(lái)源。可見(jiàn)，土壤中的熱狀況、水分狀況及其變化規(guī)律是引起凍害嚴(yán)重與否的主要因素。土體的凍脹問(wèn)題歸根結(jié)底是熱質(zhì)遷移問(wèn)題。因此，本文主要對(duì)凍土中的溫度場(chǎng)進(jìn)行分析，探討含冰量及孔隙率隨溫度變化條件下的凍土一般性瞬變熱傳導(dǎo)方程的求解問(wèn)題，研究成果有利于推動(dòng)凍土凍脹問(wèn)題的深入研究。

在凍脹過(guò)程中，土體內(nèi)溫度分布通常是不均勻的，會(huì)形成溫度梯度，如油氣管道溫度往往高于附近土體的溫度，在管壁附近的土體中水分會(huì)朝著溫度降低的方向遷移，引起對(duì)流放熱。此外，水分在凍結(jié)鋒面處由于溫度降低會(huì)發(fā)生相變，凝結(jié)成冰放出熱量，影響著凍結(jié)鋒的移動(dòng)速度。這些遷移的水分會(huì)聚集在凍結(jié)鋒面附近，形成冰透鏡體^{［2，8］}，使管道附近的凍結(jié)區(qū)土體發(fā)生膨脹，而未凍結(jié)區(qū)水分的減少又使土體發(fā)生固結(jié)，改變了土體原有的孔隙結(jié)構(gòu)，進(jìn)而對(duì)土體的滲透及變形過(guò)程產(chǎn)生影響?？梢?jiàn)，土體的凍脹過(guò)程跟溫度變化及水分遷移密切相關(guān)，而水分的遷移又受到溫度變化的影響。土體凍脹主要是由溫度的變化引起的。根據(jù)Li等^［9］、Nassar等^［10］及Hansson等^［11］的研究成果，土體凍脹的溫度場(chǎng)方程可表示為：

式中：ρ為土體密度;C為比熱容;T為土體溫度;λ為熱傳導(dǎo)系數(shù);L為土體潛熱;ρ_i為土體中凍結(jié)冰的密度;θ_i為土體孔隙中冰的含量;

式（1）是一個(gè)帶有源項(xiàng)的拋物型偏微分方程，可以定量地描述凍土中熱傳導(dǎo)過(guò)程以及土體孔隙中含冰量的變化規(guī)律。近年來(lái)，大量學(xué)者對(duì)其進(jìn)行研究，獲得了許多近似解：如Lunardini^［12］給出了同時(shí)考慮土體凍融表層熱傳導(dǎo)和相變的解析式;Lunardini等^［13］利用熱平衡積分法提出寒帶土壤凍融的Neumann問(wèn)題的精確解;Klement'Ev等^［14］對(duì)固結(jié)凍土融化過(guò)程進(jìn)行求解，獲得一維土體相變問(wèn)題的自相似解析表達(dá)式;Shao等^［15］利用傅里葉變換，得到一維土體熱傳導(dǎo)-對(duì)流方程的解析解;李方政等^［16］利用指數(shù)積分函數(shù)，獲得凍土壁擴(kuò)展深度和速度與凍結(jié)時(shí)間的理論函數(shù)表達(dá)式;Kurylyk等^［17］基于Lunardini解，對(duì)一維擴(kuò)散-對(duì)流凍土溶解溫度方程進(jìn)行求解。最近，Hu等^［18］開(kāi)始推導(dǎo)凍結(jié)土壤雙排及三排管道附近溫度場(chǎng)分布公式。

上述求解大多數(shù)是基于線性化熱傳導(dǎo)方程，將含冰量θ_i或者凍結(jié)鋒面視為時(shí)間變量的函數(shù)，采用積分、分離變量、函數(shù)展開(kāi)或傅里葉變換等方法進(jìn)行的。然而在實(shí)際工程中，含冰量通常為溫度的函數(shù)^［9］：

式中：n為土體孔隙率;S_i為孔隙中含冰量，且有：

式中：a為經(jīng)驗(yàn)系數(shù);T_f為冰凍溫度。

由式（3）～（4）可見(jiàn)，θ_i不僅與時(shí)間相關(guān)，還與空間位置相關(guān)。式（1）是一個(gè)分段非線性偏微分方程，這給理論解的尋求造成較大的困難。本文基于玻爾茲曼變換，對(duì)一類(lèi)工程中常見(jiàn)的Dirichlet邊界凍土溫度場(chǎng)公式（1）～（4）進(jìn)行求解，通過(guò)修正Chen等^［19］的顯式級(jí)數(shù)解，構(gòu)造了在定常溫Dirichlet定解條件下式（1）～（4）的高精度近似解析解，并與精確解及有限元高精度數(shù)值解對(duì)比，驗(yàn)證其正確性。

1 凍土初-邊界條件

本文主要研究的凍土溫度場(chǎng)初始-邊界條件如下：

式中：T₁gt;T₀，均為已知定值。式（5）～（6）為一類(lèi)定常溫Dirichlet定解條件，該初-邊界在凍土溫度場(chǎng)問(wèn)題中較為常見(jiàn)，如季節(jié)性?xún)鐾镣翆訜醾鲗?dǎo)問(wèn)題、凍結(jié)鑿井問(wèn)題、凍土區(qū)隧道襯砌隔熱層邊界土體的熱傳導(dǎo)問(wèn)題等^［20-22］。

2 凍土孔隙率的簡(jiǎn)化

當(dāng)只考慮溫度對(duì)凍土孔隙率的影響時(shí)，根據(jù)Lewis等^［23］的研究成果，有：

式中：K_T為土體骨架的體積模量;K_s為土顆粒的體積模量;β_s為溫度耦合系數(shù)。對(duì)式（7）關(guān)于時(shí)間t進(jìn)行積分，有：

當(dāng)初始孔隙率n₀為定值且土體顆粒不可壓縮時(shí)，有：

將式（2）～（4）及式（9）代入式（1），有：

其中，A為：

上述即為凍土孔隙率的簡(jiǎn)化情況。

3 凍土溫度場(chǎng)方程的求解

在式（10）中，令：

其中：

將式（12）代入式（10），則有：

其中：

則初-邊界條件式（5）～（6）簡(jiǎn)化為：

其中：

上述式（14）～（17）為關(guān)于U的非線性偏微分方程。近年來(lái)，學(xué)者們對(duì)其進(jìn)行了大量研究，例如，Parlance等^［24］給出其積分形式為：

式中：?=x/t^1/2為玻爾茲曼變量。Parlance等^［24］基于玻爾茲曼變換，提出了擴(kuò)散系數(shù)D（U）為冪律和指數(shù)律形式的半無(wú)限邊界的二階隱式解。Witelski^［25］基于攝動(dòng)法，對(duì)Newman邊界給出了攝動(dòng)解。Omidvar等^［26］采用微分變換及同倫攝動(dòng)法，給出無(wú)限邊界的解析解。最近，Chen等^［19］通過(guò)玻爾茲曼和級(jí)數(shù)展開(kāi)技術(shù)，給出初-邊界條件式（16）～（17）的顯式級(jí)數(shù)解：

式中：u_i（i=1，2，3，…，n）為計(jì)算參數(shù)。

將式（20）代入式（19），在U=U₁處逐項(xiàng)求導(dǎo)，有：

式中：j=0，1，2，…，n-1。

聯(lián)立式（21）即可求解得u_i。該近似解對(duì)高非線性問(wèn)題具有較好的收斂性^［27］，然而對(duì)于線性或近似線性問(wèn)題，其收斂情況并不理想，例如式（22）中無(wú)量綱化之后的熱傳導(dǎo)問(wèn)題：

此時(shí)，U=T，利用級(jí)數(shù)展開(kāi)技術(shù)^［27］進(jìn)行計(jì)算，并與精確解對(duì)比（圖1）。

當(dāng)i為偶數(shù)階時(shí)，u_i和f值均為復(fù)數(shù)，上述計(jì)算失效。由圖1可見(jiàn)，伴隨近似階數(shù)n的增加，其相對(duì)誤差逐步增加，該級(jí)數(shù)解是發(fā)散的。

本文對(duì)顯式解式（20）進(jìn)行修正，構(gòu)造玻爾茲曼變量的顯式解為：

式中：u_i（i=1，2，3，…，n）及c為計(jì)算參數(shù)。

將式（19）代入式（18），在求解域內(nèi)進(jìn)行配點(diǎn)得：

式中：j=0，1，2，3，…，n-1。

此外，令L_f為：

此為積分方程式（19）在全域內(nèi)帶權(quán)?的積分，當(dāng)?為精確解時(shí)，有：

當(dāng)?為近似解時(shí)，總存在一個(gè)c=c^*值，使得L_f取最小，即：

式中：min（X）為對(duì)X的最小值。

聯(lián)立式（20）及式（22），即可求解u_i及c。具體求解步驟為：

（1） 預(yù)先設(shè)定若干c值，分別代入式（24），求解對(duì)應(yīng)的u_i;

（2） 將步驟（1）中的c及u_i代入式（25），計(jì)算對(duì)應(yīng)的L;

（3） 比較步驟（2）中的L_f值，選擇最小的L_f，其對(duì)應(yīng)的c及u_i即可滿足式（27）中的L_f|_c=c^*。

4 凍土溫度場(chǎng)數(shù)值模擬

對(duì)一維凍土溫度場(chǎng)Dirichlet定解條件進(jìn)行分析，采用第3節(jié)提出的計(jì)算方法分別對(duì)其進(jìn)行計(jì)算，藉此驗(yàn)證本文方法的正確性。

4.1 凍土溫度場(chǎng)線性熱傳導(dǎo)問(wèn)題模擬

當(dāng)Tgt;T_f時(shí)，凍土溫度場(chǎng)方程退化為普通線性熱傳導(dǎo)方程，將其無(wú)量綱化之后即為式（22），采用式（23）～（27）計(jì)算1～3階近似解u_i和c值，如表1所列。

將表1中計(jì)算參數(shù)代入式（23），可得玻爾茲曼變量1～3階近似解，并與精確解對(duì)比（圖2）。

由圖2可見(jiàn)，隨著玻爾茲曼變量的增加，溫度在降低，當(dāng)?=x/t^1/2在4.65附近時(shí)，溫度趨近于0。本節(jié)計(jì)算了1～3階近似解，伴隨著階數(shù)增加，其計(jì)算結(jié)果是逐步收斂于精確解的。表2給出了3階近似解與精確解計(jì)算結(jié)果及其相對(duì)誤差。

在表2中，3階近似解與精確解的相對(duì)誤差從T=1 ℃開(kāi)始逐步增大，在T=0.7 ℃處達(dá)到最大，為4.072 9%，而后隨著溫度的降低逐漸減小，在靠近0處回落到3.159 13%。由表2計(jì)算結(jié)果可知，目前的結(jié)果與精確解較為吻合。對(duì)比圖1可知，引入?yún)?shù)c之后，對(duì)于線性問(wèn)題可以改善玻爾茲曼變量顯式解式（20）的收斂性，加快其收斂速度。

4.2 凍土分段非線性熱傳導(dǎo)問(wèn)題模擬

對(duì)于凍土熱傳導(dǎo)問(wèn)題，根據(jù)Li等^{［2，9］}的研究成果，選擇土體參數(shù)（表3）。

選擇邊界溫度T₁為20 ℃，初始溫度T₀為-10 ℃，初始孔隙率為0.33，則式（10）可變?yōu)椋?/p>

其中，A為：

則式（28）的近似解析解為：

式中：U及U₁詳見(jiàn)式（12）～（13）和式（18）。采用式（23）～（27）計(jì)算1～3階近似解，u_i和c值如表4所列。

將表4中計(jì)算參數(shù)代入式（23），可得玻爾茲曼變量1～3階近似解，并與有限元（FEM）解（劃分100 000個(gè)單元，可視為精確解）對(duì)比，如圖3所示。

由圖3可見(jiàn)，對(duì)于分段非線性的凍土熱傳導(dǎo)問(wèn)題，本文提出的玻爾茲曼變量?=x/t^1/2解析解是隨著階數(shù)n的增加逐步收斂于FEM解的，當(dāng)n=3時(shí)，已經(jīng)與FEM精確解較為吻合。表5給出了3階近似解的相對(duì)誤差。

由表5中可知，3階近似解析解的計(jì)算結(jié)果與FEM解較為接近，其最大相對(duì)誤差為4.776 9%，位于T=12 ℃位置附近，在冰凍溫度分界點(diǎn)-0.35 ℃處的誤差為1.993 6%。此外，計(jì)算10 d，20 d，40 d，60 d，80 d及100 d的溫度分布如圖4所示。

由圖4可見(jiàn)，加入了計(jì)算參數(shù)c的3階近似解析解可以用于預(yù)測(cè)不同時(shí)間的溫度分布。以表3給出的計(jì)算參數(shù)可知，在100 d時(shí)，在x=0.22 m附近溫度開(kāi)始升高，相對(duì)于初始時(shí)刻，冰凍臨界點(diǎn)T=-0.35 ℃向右推進(jìn)了約0.055 m。目前的計(jì)算結(jié)果與FEM高精度數(shù)值解吻合較好。

5 結(jié)語(yǔ)

本文對(duì)含冰量及孔隙率同時(shí)變化時(shí)凍土熱傳導(dǎo)方程進(jìn)行求解，主要結(jié)論如下：

（1） 基于玻爾茲曼變換，對(duì)Chen等^［19］的顯式解析解進(jìn)行修正，構(gòu)造出收斂性較好的凍土熱傳導(dǎo)方程級(jí)數(shù)解析解;

（2） 將上述解應(yīng)用于線性無(wú)量綱化的熱傳導(dǎo)方程及含冰量為分段非線性函數(shù)的凍土熱傳導(dǎo)方程的求解，獲得1～3階近似解析解，并與精確解及高精度FEM解進(jìn)行對(duì)比，對(duì)比結(jié)果證實(shí)，本文的求解方法可以同時(shí)求解線性及非線性熱傳導(dǎo)方程。

參考文獻(xiàn)（References）

［1］王澄海，靳雙龍，施紅霞.未來(lái)50 a中國(guó)地區(qū)凍土面積分布變化［J］.冰川凍土，2014，36（1）：1-8.

WANG Chenghai，JIN Shuanglong，SHI Hongxia.Area change of the frozen ground in China in the next 50 years［J］.Journal of Glaciology and Geocryology，2014，36（1）：1-8.

［2］李智明.基于復(fù)合混合物理論的凍土多場(chǎng)耦合研究［D］.哈爾濱：哈爾濱工業(yè)大學(xué)，2021.

LI Zhiming.Research on coupling of multifield in frozen soil based on hybrid mixture theory［D］.Harbin：Harbin Institute of Technology，2021.

［3］王小軍，曾輝輝，賈海鋒，等.青藏鐵路多年凍土區(qū)次生不良凍土現(xiàn)象的調(diào)查分析［J］.鐵道工程學(xué)報(bào)，2003，20（增刊1）：1-4.

WANG Xiaojun，ZENG Huihui，JIA Haifeng，et al.Investigation and analysis of seocndary unfavourable phenomenon of permafrost in permafrost area on Qinghai-Tibet railway［J］.Journal of Railway Engineering Society，2003，20（Suppl01）：1-4.

［4］沈宇鵬，王篤禮，林園榕，等.越冬基坑水平凍脹的防治措施效果分析［J］.巖土力學(xué)，2021，42（5）：1434-1442.

SHEN Yupeng，WANG Duli，LIN Yuanrong，et al.On the effect of prevention measures against horizontal frost heave of foundation pits over winter［J］.Rock and Soil Mechanics，2021，42（5）：1434-1442.

［5］歐爾峰，黃志軍.寒冷地區(qū)水渠的熱力耦合數(shù)值分析［J］.蘭州交通大學(xué)學(xué)報(bào)，2010，29（4）：66-71.

OU Erfeng，HUANG Zhijun.Numerical analysis of the temperature and stress fields of ditch in cold regions［J］.Journal of Lanzhou Jiaotong University，2010，29（4）：66-71.

［6］何平，程國(guó)棟，朱元林.土體凍結(jié)過(guò)程中的熱質(zhì)遷移研究進(jìn)展［J］.冰川凍土，2001，23（1）：92-98.

HE Ping，CHENG Guodong，ZHU Yuanlin.The progress of study on heat and mass transfer in freezing soils［J］.Journal of Glaciology and Geocryology，2001，23（1）：92-98.

［7］程國(guó)棟，何平.多年凍土地區(qū)線性工程建設(shè)［J］.冰川凍土，2001，23（3）：213-217.

CHENG Guodong，HE Ping.Linearity engineering in permafrost areas［J］.Journal of Glaciolgy and Geocryology，2001，23（3）：213-217.

［8］BESKOW G.Soil freezing and frost heaving with special application to roads and railroads［J］.Cold Regions Research and Engineering，1947，26（3）：410-416.

［9］LI Z M，CHEN J，SUN K，et al.Numerical simulation and experimental validation of moisture-heat coupling for saturated frozen soils［J］.Sciences in Cold and Arid Regions，2017，9（3）：250-257.

［10］NASSAR I N，HORTON R.Simultaneous transfer of heat，water，and solute in porous media：I.theoretical development［J］.Soil Science Society of America Journal，1992，56（5）：1350.

［11］HANSSON K，SIMUNEK J，MIZOGUCHI M，et al.Water flow and heat transport in frozen soil：numerical solution and freeze-thaw applications［J］.Vadose Zone Journal，2004，3（2）：693-704.

［12］LUNARDINI V J.Heat transfer with freezing and thawing［M］.New York：Elsevier Science Publishing Company，1991.

［13］LUNARDINI V J，VAROTTA R.Approximate solution to Neumann problem for soil systems［J］.Journal of Energy Resources Technology，1981，103（1）：76-81.

［14］KLEMENT'EV A F，KLEMENT'EVA E A.Self-similar solution of the problem of consolidation and thawing of frozen soil［J］.Journal of Engineering Physics，1988，54（4）：463-466.

［15］SHAO M G，HORTON R，JAYNES D B.Analytical solution for one-dimensional heat conduction-convection equation［J］.Soil Science Society of America Journal，1998，62（1）：123-128.

［16］李方政，夏明萍.基于指數(shù)積分函數(shù)的人工凍土溫度場(chǎng)解析研究［J］.東南大學(xué)學(xué)報(bào)（自然科學(xué)版），2004，34（4）：469-473.

LI Fangzheng，XIA Mingping.Study on analytical solution of temperature field of artificial frozen soil by exponent-integral function［J］.Journal of Southeast University （Natural Science Edition），2004，34（4）：469-473.

［17］KURYLYK B L，MCKENZIE J M，MACQUARRIE K T B，et al.Analytical solutions for benchmarking cold regions subsurface water flow and energy transport models：one-dimensional soil thaw with conduction and advection［J］.Advances in Water Resources，2014，70：172-184.

［18］HU X D，HAN L，HAN Y G.Analytical solution to temperature distribution of frozen soil wall by multi-row-piped freezing with the boundary separation method［J］.Applied Thermal Engineering，2019，149：702-711.

［19］CHEN X，HE D D，YANG G Y，et al.Approximate analytical solution for Richards' equation with finite constant water head Dirichlet boundary conditions［J］.Computational and Applied Mathematics，2021，40（7）：236.

［20］羅凱喬.季節(jié)性?xún)鐾翢醾鲗?dǎo)性質(zhì)研究［D］.哈爾濱：東北林業(yè)大學(xué)，2016.

LUO Kaiqiao.Study on thermal conductivity property of seasonal frozen soil［D］.Harbin：Northeast Forestry University，2016.

［21］吳嘉林，辛德林，張建平.井工開(kāi)采技術(shù)的創(chuàng)新與發(fā)展［J］.煤炭工程，2014，46（1）：4-8.

WU Jialin，XIN Delin，ZHANG Jianping.Innovation and development of underground coal mining technology［J］.Coal Engineering，2014，46（1）：4-8.

［22］夏才初，范東方，李志厚，等.隧道多年凍土段隔熱層厚度解析計(jì)算結(jié)果的探討［J］.土木工程學(xué)報(bào)，2015，48（2）：118-124.

XIA Caichu，F(xiàn)AN Dongfang，LI Zhihou，et al.Discussion on analytical calculation for thermal-insulation layer thickness of tunnel in permafrost area［J］.China Civil Engineering Journal，2015，48（2）：118-124.

［23］LEWIS R W，SHREFLER B A.The finite element method in the static and dynamic deformation and consolidation of porous media［M］.2nd ed.New York：John Wiley，1998.

［24］PARLANCE M B，PRASAD S N，PARLANGE J Y，et al.Extension of the Heaslet-Alksne technique to arbitrary soil water diffusivities［J］.Water Resources Research，1992，28（10）：2793-2797.

［25］WITELSKI T P.Motion of wetting fronts moving into partially pre-wet soil［J］.Advances in Water Resources，2005，28（10）：1133-1141.

［26］OMIDVAR M，BARARI A，MOMENI M，et al.New class of solutions for water infiltration problems in unsaturated soils［J］.Geomechanics and Geoengineering，2010，5（2）：127-135.

［27］CHEN X，DAI Y.An approximate analytical solution of Richards equation with finite boundary［J］.Boundary Value Problems，2017，2017（1）：167.

（本文編輯：張向紅）