本發(fā)明屬于計算土力學(xué)和工程災(zāi)害防治
技術(shù)領(lǐng)域：
，涉及一種數(shù)值計算方法，尤其涉及一種膨脹土在降雨入滲條件下增濕膨脹數(shù)值模擬方法。
背景技術(shù)：
：全世界分布有膨脹土的國家有40多個，而相對其他國家，我國的膨脹土問題尤為突出。首先，膨脹土在我國分布范圍很廣，基本覆蓋了西南云貴高原到華北平原之間的大片土地。其次，膨脹土種類較多，不同地區(qū)膨脹土的形成原因及工程特性不盡相同，具有顯著的地域性。隨著我國基礎(chǔ)設(shè)施建設(shè)進(jìn)程加速，許多工程建設(shè)將不可避免地穿越膨脹土地區(qū)。由于我國現(xiàn)有膨脹土地區(qū)建筑技術(shù)規(guī)范還不成熟和完善，膨脹土地區(qū)工程設(shè)計施工多依據(jù)以往工程經(jīng)驗。然而膨脹土具有顯著的地域差異性，盲目套用其他地區(qū)膨脹土設(shè)計施工經(jīng)驗和參數(shù)將導(dǎo)致許多嚴(yán)重工程問題的產(chǎn)生，造成巨大的經(jīng)濟損失和人員傷亡。眾多膨脹土工程實例表明，降雨入滲是誘發(fā)膨脹土工程災(zāi)害事故發(fā)生的主要元兇。然而，由于其涉及復(fù)雜的非飽和滲流和膨脹變形過程，采用理論求解時，為計算方便，常對其實際情況進(jìn)行過度簡化，導(dǎo)致計算結(jié)果精度較低，且耗時耗力。雖然目前已經(jīng)探索了相關(guān)膨脹土增濕膨脹變形的數(shù)值模擬方法，但仍存在一些不足，例如，基于“濕度應(yīng)力場”理論，采用溫度場模擬降雨增濕時，其實際上只適用于飽和滲流分析，不能考慮降雨過程中雨水重力以及非飽和區(qū)基質(zhì)吸力和滲透系數(shù)變化的影響，因此不能真實反映膨脹土在降雨入滲條件下的非飽和滲流過程。同時在考慮膨脹土膨脹效應(yīng)時，通常是通過在單元上施加垂直于臨空面的外力來模擬膨脹力，但實際中膨脹變形是向四周發(fā)展的，膨脹力并不能單純以外力方式施加。綜上所述，目前針對膨脹土增濕膨脹數(shù)值模擬研究中，均存在一定的問題，尚不能合理且全面地模擬膨脹土在增濕條件下的非飽和滲流和膨脹變形過程，計算結(jié)果與實際情況存在較大誤差，不能很好地指導(dǎo)膨脹土地區(qū)工程的設(shè)計和施工。技術(shù)實現(xiàn)要素：本發(fā)明的目的在于解決現(xiàn)有數(shù)值計算中不能真實反映膨脹土在增濕條件下膨脹變形的技術(shù)問題，進(jìn)而提供了一種合理且全面考慮降雨入滲條件下非飽和滲流和膨脹變形過程的數(shù)值計算方法。為了實現(xiàn)上述目的，本發(fā)明采用如下技術(shù)方案：一種膨脹土在降雨入滲條件下增濕膨脹數(shù)值模擬方法，其特征在于：所述方法包括以下步驟：1)建立膨脹土工程力學(xué)計算模型；采用FLAC3D程序建立三維幾何模型并劃分計算網(wǎng)格，劃分有計算網(wǎng)格的三維幾何模型包括模型單元以及模型節(jié)點，所述膨脹土工程力學(xué)計算模型是根據(jù)具體膨脹土工程的實際情況建立的；2)對步驟1)所建立的膨脹土工程力學(xué)模型進(jìn)行初始平衡地應(yīng)力計算；3)結(jié)束步驟2)的計算后進(jìn)行非飽和滲流計算：4)結(jié)束步驟3)的計算后進(jìn)行熱傳導(dǎo)計算；5)對數(shù)值模擬結(jié)果進(jìn)行分析。作為優(yōu)選，本發(fā)明所采用的步驟3)的具體實現(xiàn)方式是：3.1)對模型單元設(shè)置相應(yīng)的滲流參數(shù)，并根據(jù)具體膨脹土工程的實際情況設(shè)置初始飽和度和孔隙水壓力，膨脹土工程力學(xué)模型的上邊界設(shè)置為流量邊界以模擬降雨入滲過程，并設(shè)置滲流計算時間；3.2)在每一滲流計算時步中，提取模型節(jié)點飽和度Sr，并判斷模型節(jié)點飽和度Sr是否等于1.0；若是，則直接進(jìn)行步驟3.3)；若否，則依據(jù)通過試驗得出的土-水特征曲線擬合公式計算該模型節(jié)點的基質(zhì)吸力值，并將模型節(jié)點的基質(zhì)吸力值賦值給節(jié)點的孔隙水壓力后進(jìn)行步驟3.3)；3.3)首先通過反距離加權(quán)插值法，將模型節(jié)點飽和度轉(zhuǎn)化為模型單元飽和度，然后根據(jù)滲透系數(shù)和強度參數(shù)與飽和度之間的關(guān)系計算相應(yīng)飽和度下模型單元的滲透系數(shù)以及強度參數(shù)值，并將相應(yīng)飽和度下模型單元的滲透系數(shù)以及強度參數(shù)值賦給模型單元；3.4)檢查地表節(jié)點的孔隙水壓力pˊ，若孔隙水壓力大于0，將該地表節(jié)點的流量邊界修改為壓力邊界，固定該地表節(jié)點的孔隙水壓力為0；若地表節(jié)點的孔隙水壓力小于或等于零，則直接進(jìn)行步驟3.5)，所述地表節(jié)點是膨脹土工程力學(xué)模型上邊界節(jié)點；3.5)判斷是否達(dá)到非飽和滲流計算的終止條件，若是，則結(jié)束計算過程并保存結(jié)果文件；若否，重復(fù)步驟3.2)-3.4)直至達(dá)到非飽和滲流計算的終止條件；所述非飽和滲流計算的終止條件是滲流計算時間。作為優(yōu)選，本發(fā)明所采用的步驟4)的具體實現(xiàn)方式是：4.1)基于濕度應(yīng)力場理論，根據(jù)滲流連續(xù)性微分方程和熱傳導(dǎo)微分方程的相似性，建立滲流參數(shù)與熱力學(xué)參數(shù)的等效對應(yīng)關(guān)系；其中，所述滲流連續(xù)性微分方程的表達(dá)形式是：∂∂x(kx∂hm∂x)+∂∂y(ky∂hm∂y)+∂∂z(kz∂hm∂z)=Cw∂hm∂t]]>所述熱傳導(dǎo)微分方程的表達(dá)形式是：∂∂x(λx∂T∂x)+∂∂y(λy∂T∂y)+∂∂z(λz∂T∂z)=ρCv∂T∂t]]>式中：kx、ky以及kz分別為x、y以及z這三個方向的滲透系數(shù)；hm是基質(zhì)吸力水頭；Cw是比水容重；λx、λy以及λz分別為x、y以及z這三個方向的熱傳導(dǎo)系數(shù)；T是溫度；ρ是介質(zhì)密度；Cv是介質(zhì)的比熱容；t是時間；比較滲流連續(xù)性微分方程以及熱傳導(dǎo)微分方程可知，滲流參數(shù)與熱力學(xué)參數(shù)存在等效對應(yīng)關(guān)系，其中，滲透系數(shù)ki對應(yīng)熱傳導(dǎo)系數(shù)λi；基質(zhì)吸力水頭hm對應(yīng)溫度T；比水容重Cw對應(yīng)比熱容ρCv；根據(jù)以上對應(yīng)關(guān)系，確定出熱傳導(dǎo)計算中模型單元和模型節(jié)點的熱力學(xué)參數(shù)；4.2)提取步驟3)計算結(jié)束后保存的結(jié)果文件，給模型單元賦予上述按等效原理計算出的熱力學(xué)參數(shù)；4.3)設(shè)定當(dāng)模型單元含水率從殘余含水率θr增加到飽和含水率θs時對應(yīng)的溫度為100℃，各模型節(jié)點對應(yīng)的等效溫度Tp采用線性插值法計算，計算公式如下：Tp=100×Δθθs-θr]]>4.4)固定降雨入滲邊界處模型節(jié)點的溫度，并進(jìn)行熱傳導(dǎo)計算，在每一熱力計算時步中，判斷各模型節(jié)點溫度是否全部達(dá)到其等效溫度Tp，若是，則結(jié)束計算過程；若否，則重復(fù)進(jìn)行步驟4.3)直至達(dá)到等效溫度Tp后結(jié)束計算。作為優(yōu)選，本發(fā)明所采用的步驟2)中是FLAC3D程序采用顯式有限差分計算方法進(jìn)行迭代結(jié)算，在迭代過程中通過監(jiān)控最大不平衡力比率，檢測計算是否達(dá)到收斂，一旦最大不平衡力比率小于預(yù)先設(shè)置的值，則計算過程終止。作為優(yōu)選，本發(fā)明所采用的步驟5)中的計算結(jié)果分析包括但不限于不同降雨歷時下圍巖含水率分布、孔隙水壓力分布、圍巖應(yīng)力以及工程結(jié)構(gòu)的受力和變形。與現(xiàn)有技術(shù)相比，本發(fā)明的優(yōu)點是：(1)通過本發(fā)明的數(shù)值計算方法，能夠全面模擬膨脹土在降雨入滲條件下的增濕膨脹變形全過程，克服了現(xiàn)有數(shù)值計算方法的不足和缺陷，從而能夠精確計算膨脹土在增濕膨脹變形過程中對工程結(jié)構(gòu)受力和變形的影響，揭示膨脹土工程災(zāi)害孕育機理和演化規(guī)律，為膨脹土地區(qū)相關(guān)工程的設(shè)計和施工提供參考和依據(jù)。(2)本發(fā)明的數(shù)值計算方法具有簡單、靈活以及通用性強的特點，方法便于實現(xiàn)，既能保證計算精度，又能節(jié)約時間成本。附圖說明圖1是本發(fā)明膨脹土增濕膨脹數(shù)值模擬方法流程圖；圖2是隧道平面力學(xué)模型網(wǎng)格劃分圖；圖3是采用流-固耦合模塊模擬非飽和滲流過程的數(shù)值計算流程圖；圖4是采用熱力學(xué)模塊模擬膨脹變形過程的數(shù)值計算流程圖；圖5是隧道初期支護(hù)監(jiān)測點布置示意圖；圖6是初期支護(hù)彎矩隨降雨歷時變化曲線；圖7是降雨過程中初期支護(hù)變形轉(zhuǎn)化示意圖。具體實施方式下面結(jié)合附圖和具體的實施例對本發(fā)明的膨脹土在降雨入滲條件下增濕膨脹數(shù)值模擬方法做進(jìn)一步的詳細(xì)說明：膨脹土在降雨入滲條件下增濕膨脹數(shù)值模擬方法流程圖見附圖1，本發(fā)明以某淺埋膨脹土隧道工程為例，對其在降雨入滲條件下支護(hù)結(jié)構(gòu)的動力響應(yīng)進(jìn)行模擬分析，具體按照以下步驟實施：1、建立淺埋膨脹土隧道力學(xué)計算模型；某山嶺隧道為淺埋膨脹土隧道，單洞凈寬10.5m，凈高5.0m，洞口段埋深8.4m。為簡化計算，考慮圍巖處于平面應(yīng)變狀態(tài)，取模型厚度1m。為減小邊界效應(yīng)的影響，模型寬度取100m，隧道底部至模型下邊界35m，埋深8.4m，如附圖2所示。對淺埋膨脹土隧道力學(xué)計算模型四周施加法向約束，并將底部固定，淺埋膨脹土隧道力學(xué)計算模型四周及底部均設(shè)置為不透水邊界，頂部設(shè)置為降雨邊界。圍巖采用Mohr-Coulomb模型，初期支護(hù)采用Shell單元，厚度為26cm，其彈性模量由鋼拱架、鋼筋網(wǎng)片和噴射混凝土按剛度等效原理進(jìn)行折算，圍巖及支護(hù)結(jié)構(gòu)的力學(xué)參數(shù)列于表1，其中土體飽和滲透系數(shù)ksat取為10-4cm/s，非飽和區(qū)滲透系數(shù)是關(guān)于飽和度的函數(shù)，可采用以下經(jīng)驗公式進(jìn)行計算。k(s)＝ksats2(3-2s)(1)表1圍巖及支護(hù)結(jié)構(gòu)的力學(xué)參數(shù)密度/(kg/m3)彈性模量/GPa泊松比圍巖20000.080.35錨桿—2100.25初支2500230.20采用濾紙法測定現(xiàn)場膨脹土樣的基質(zhì)吸力，并以VanGenuchten(V-G)提出的土-水特征曲線模型進(jìn)行參數(shù)擬合，V-G模型表達(dá)式及相關(guān)擬合參數(shù)如下：θ=θr+θs-θr(1+(aψ)n)m---(2)]]>式中：θ為含水率，θs為飽和含水率，θr為殘余含水率，ψ為基質(zhì)吸力。擬合參數(shù)a＝0.0116，n＝2.2610，m＝0.2841，θr＝9.02％，θs＝31.07％。通過室內(nèi)直剪試驗測得現(xiàn)場膨脹土樣的黏聚力c和內(nèi)摩擦角隨含水率θ變化的擬合關(guān)系式為：2、進(jìn)行淺埋膨脹土隧道力學(xué)計算模型初始平衡狀態(tài)求解計算；FLAC3D程序采用顯示有限差分計算方法進(jìn)行迭代結(jié)算，在迭代過程中通過監(jiān)控最大不平衡力比率，檢測計算是否達(dá)到收斂，一旦最大不平衡力比率小于預(yù)先設(shè)置的值，則計算過程終止。3、進(jìn)行非飽和滲流計算；采用FLAC3D內(nèi)置的流-固耦合模塊進(jìn)行非飽和滲流計算的流程圖如附圖3所示，具體過程表述如下：①對模型單元設(shè)置相應(yīng)的滲流參數(shù)，并根據(jù)實際情況設(shè)置初始飽和度和孔隙水壓力，并將模型上邊界設(shè)置為降雨入滲邊界，本實例中假設(shè)降雨強度較大，地表在短時間內(nèi)便達(dá)到飽和狀態(tài)，進(jìn)入壓力入滲階段，因此設(shè)置地表孔隙水壓力為零，并設(shè)置滲流計算時間為72h。②在每一滲流計算時步中，提取模型節(jié)點飽和度Sr，若模型節(jié)點飽和度不等于1.0，則依據(jù)土-水特征曲線擬合公式求得該模型節(jié)點的基質(zhì)吸力值，并將該模型節(jié)點的基質(zhì)吸力值賦值給相應(yīng)節(jié)點的孔隙水壓力(負(fù)孔隙水壓)，若飽和度等于1.0，則直接進(jìn)行下一步計算；③由于強度參數(shù)以及滲透系數(shù)是單元變量，而飽和度是節(jié)點變量，因此首先通過反距離加權(quán)插值法，將模型節(jié)點飽和度轉(zhuǎn)化為模型單元飽和度，然后計算相應(yīng)飽和度下模型單元的滲透系數(shù)以及強度參數(shù)值，并將其值賦給模型單元。④檢查地表節(jié)點的孔隙水壓力pˊ，若孔隙水壓力大于0，將該地表節(jié)點的流量邊界修改為壓力邊界，固定其孔隙水壓力為0。(由于本實例中將地表空隙水壓力設(shè)置為零來模擬強降雨情況，因此可跳過此步驟)⑤檢查是否達(dá)到滲流計算終止條件(預(yù)先設(shè)定的滲流計算時間72h)，若未達(dá)到則重復(fù)步驟②～⑤，否則結(jié)束計算。4、進(jìn)行熱傳導(dǎo)計算；采用FLAC3D內(nèi)置的熱力學(xué)模塊進(jìn)行膨脹變形計算的流程圖如附圖4所示，具體過程表述如下：①基于“濕度應(yīng)力場理論”，根據(jù)滲流連續(xù)性微分方程和熱傳導(dǎo)微分方程的相似性，建立滲流參數(shù)與熱力學(xué)參數(shù)的等效轉(zhuǎn)換關(guān)系。其中，滲流連續(xù)性微分方程和熱傳導(dǎo)微分方程的表達(dá)形式如下：∂∂x(kx∂hm∂x)+∂∂y(ky∂hm∂y)+∂∂z(kz∂hm∂z)=Cw∂hm∂t---(4)]]>∂∂x(λx∂T∂x)+∂∂y(λy∂T∂y)+∂∂z(λz∂T∂z)=ρCv∂T∂t---(5)]]>式中：kx、ky、kz分別為三個主方向的滲透系數(shù)，hm為基質(zhì)吸力水頭，Cw為比水容重，λx、λy、λz分別為三個主方向的熱傳導(dǎo)系數(shù)，T為溫度，ρ為介質(zhì)密度，Cv為介質(zhì)的比熱容。比較式(4)和式(5)可知，滲流參數(shù)與熱力學(xué)參數(shù)存在等效關(guān)系，其中，滲透系數(shù)ki對應(yīng)熱傳導(dǎo)系數(shù)λi；基質(zhì)吸力水頭hm對應(yīng)溫度T；比水容重Cw對應(yīng)比熱容ρCv。根據(jù)以上對應(yīng)關(guān)系，可確定出熱傳導(dǎo)計算中，模型單元和節(jié)點對應(yīng)的熱力學(xué)參數(shù)。②提取非飽和滲流計算結(jié)果文件，給模型單元賦予上述按等效原理計算出的熱力學(xué)參數(shù)(λi、Cv、α)。③設(shè)定當(dāng)模型單元含水率從殘余含水率θr增加到飽和含水率θs時對應(yīng)的溫度為100℃，各模型節(jié)點對應(yīng)的等效溫度Tp采用線性插值法計算，計算公式如下：Tp=100×Δθθs-θr---(6)]]>④固定降雨入滲邊界處模型節(jié)點的溫度，并進(jìn)行熱傳導(dǎo)計算，檢查各模型節(jié)點溫度是否全部達(dá)到其等效溫度，如未達(dá)到則繼續(xù)熱傳導(dǎo)計算，否則，結(jié)束計算。5、對數(shù)值模擬結(jié)果進(jìn)行分析。本實施步驟主要包括：降雨入滲過程中圍巖孔隙水壓力和圍巖應(yīng)力分布變化規(guī)律分析以及支護(hù)結(jié)構(gòu)受力和變形情況分析。(1)降雨入滲過程中圍巖孔隙水壓力和圍巖應(yīng)力分布變化規(guī)律分析依據(jù)上一步計算所得的不同降雨歷時下圍巖孔隙水壓力以及圍巖應(yīng)力分布圖，分析降雨入滲過程中圍巖孔隙水壓力和圍巖應(yīng)力分布變化規(guī)律?？梢钥闯鲭S降雨歷時增加，浸潤鋒線向下部推進(jìn)，降雨影響范圍內(nèi)圍巖孔隙水壓力升高(基質(zhì)吸力減小)。降雨過程中圍巖垂直應(yīng)力變化不大，但水平應(yīng)力顯著增大，降雨前拱頂處水平應(yīng)力為0.06MPa，降雨72h后拱頂水平應(yīng)力增加至0.33MPa，且在隧道拱頂上方產(chǎn)生應(yīng)力集中。(2)支護(hù)結(jié)構(gòu)受力和變形情況分析取襯砌結(jié)構(gòu)上5處進(jìn)行監(jiān)測，監(jiān)測點布置如附圖5所示，監(jiān)測點處彎矩隨降雨歷時變化曲線如附圖6所示，其中正值表示襯砌臨空一側(cè)受拉，負(fù)值表示襯砌圍巖一側(cè)受拉。從圖中可看出，降雨增濕過程中，襯砌在監(jiān)測點4(即拱腳)處彎矩絕對值最大。拱頂處彎矩隨降雨歷時增加逐漸減小至0，后反向增大；而拱腰及邊墻處逐漸由負(fù)彎矩轉(zhuǎn)變?yōu)檎龔澗?。這是由于在降雨初期，垂直荷載較大，支護(hù)結(jié)構(gòu)在上覆荷載作用下發(fā)生豎向擠壓變形，此時拱頂處支護(hù)臨空一側(cè)受拉，拱腰和邊墻處支護(hù)圍巖一側(cè)受拉；降雨持續(xù)一段時間后，圍巖增濕膨脹效應(yīng)引起水平壓力逐漸增大，支護(hù)結(jié)構(gòu)逐漸由豎向擠壓變形轉(zhuǎn)化為水平向擠壓變形，這使得拱頂、拱腰及邊墻處襯砌彎矩發(fā)生正負(fù)值變化，降雨入滲過程中初期支護(hù)結(jié)構(gòu)變形轉(zhuǎn)化示意圖如附圖7所示。最后說明的是，以上實施例僅用以說明本發(fā)明的技術(shù)方案而非限制，盡管參照實施例對本發(fā)明進(jìn)行了詳細(xì)說明，凡對本發(fā)明的技術(shù)方案進(jìn)行修改或者等同替換，而不脫離本發(fā)明技術(shù)方案的宗旨和范圍，其均應(yīng)涵蓋在本發(fā)明的權(quán)利要求范圍當(dāng)中。當(dāng)前第1頁1 2 3