葉東生 朱媛媛 王笑梅 王玉善

摘 要： 研究了分層不可壓流體飽和熱彈性多孔介質(zhì)軸對(duì)稱問(wèn)題的動(dòng)力響應(yīng)問(wèn)題，基于多孔介質(zhì)理論（PMT），給出了該問(wèn)題的數(shù)學(xué)模型.在空間域內(nèi)采用微分求積單元法（DQEM）設(shè)置離散的控制微分方程、邊界條件和連接條件，在時(shí)間域內(nèi)采用二階向后差分格式處理時(shí)間導(dǎo)數(shù).在離散化的初始條件下，運(yùn)用Newton-Raphson法進(jìn)行迭代求解，得到各離散點(diǎn)處未知物理量的數(shù)值結(jié)果.研究表明：該方法有效、可靠，且具備精度較高、計(jì)算量較小、數(shù)值穩(wěn)定等優(yōu)點(diǎn).

關(guān)鍵詞： 流體飽和熱彈性多孔介質(zhì); 多孔介質(zhì)理論（PMT）; 微分求積單元法（DQEM）; 動(dòng)力學(xué)響應(yīng)

中圖分類號(hào)： TU 311文獻(xiàn)標(biāo)志碼： A文章編號(hào)： 1000-5137（2019）02-0141-10

0 引 言

飽和多孔介質(zhì)熱-流-固耦合系統(tǒng)的研究不僅在土力學(xué)、水文學(xué)等經(jīng)典應(yīng)用領(lǐng)域發(fā)揮重要作用，也已成為許多新興學(xué)科和應(yīng)用技術(shù)發(fā)展的關(guān)鍵，其相關(guān)理論和數(shù)值方法的研究工作具有重要的理論意義和廣泛的應(yīng)用背景.

1955年，BIOT[1]建立了飽和多孔介質(zhì)熱彈性和熱動(dòng)力理論.隨后，CARTER等[2]研究了飽和土球形熱源附近的固結(jié)問(wèn)題.CUI等[3]引入熱孔隙度狀態(tài)表面的概念，提出了飽和土熱-水-力耦合分析的理論模型.劉干斌等[4]通過(guò)對(duì)Biot波動(dòng)方程的修正，研究了簡(jiǎn)諧均布荷載作用下，多孔彈性地基土體的熱-水-力耦合動(dòng)力響應(yīng)問(wèn)題.白冰[5]對(duì)循環(huán)溫度荷載下飽和多孔介質(zhì)熱-水-力耦合響應(yīng)的一維情形進(jìn)行了研究.戴清晨等[6]研究了熱局部非平衡條件下，橫觀各向同性飽和多孔介質(zhì)中柱形空洞的熱應(yīng)力.de BOER等[7-8]基于連續(xù)介質(zhì)混和物公理體系和體積分?jǐn)?shù)概念，建立了較完整的多孔介質(zhì)理論.de BOER等[9]利用拉普拉斯變換，給出了多孔介質(zhì)一維動(dòng)力學(xué)響應(yīng)的解析解.劉占芳等[10]等利用拉普拉斯變換和卷積定理，得到了邊界自由排水時(shí)，任意應(yīng)力和位移邊界條件下，瞬態(tài)波動(dòng)過(guò)程的解析表達(dá).李向約等[11]推導(dǎo)了描述飽和多孔介質(zhì)在熱固結(jié)耦合作用下的數(shù)學(xué)模型.HE等[12]建立了熱局部非平衡條件下飽和不可壓多孔彈性介質(zhì)熱-流-固耦合模型.YANG[13]建立了相應(yīng)的Gurtin型廣義變分原理.朱媛媛等[14]利用微分求積法（DQM）分析了流體飽和熱彈性多孔介質(zhì)圓柱體的動(dòng)力響應(yīng).嚴(yán)俊等[15]利用連續(xù)介質(zhì)理論，提出了非飽和多孔介質(zhì)的熱本構(gòu)關(guān)系以及孔隙流體的熱運(yùn)動(dòng)規(guī)律.

本文作者分析了分層流體飽和熱彈性多孔介質(zhì)的動(dòng)力學(xué)行為，基于多孔介質(zhì)理論，建立了一維流體飽和熱彈性多孔介質(zhì)軸對(duì)稱問(wèn)題的數(shù)學(xué)模型，采用微分求積單元法（DQEM）和微分-代數(shù)方法求解耦合系統(tǒng)的動(dòng)力學(xué)問(wèn)題，獲得耦合系統(tǒng)一些定性的結(jié)果.將DQEM用于求解分層流體飽和熱彈性問(wèn)題分析中，研究工作有一定的探索性，建模思路和研究方法可為工程實(shí)踐提供參考.

1 問(wèn)題的數(shù)學(xué)描述

1.1 控制微分方程

圖1為由n層介質(zhì)組成的分層軸對(duì)稱不可壓流體飽和多孔熱彈性體.其中，（x，y）表示直角坐標(biāo)系，（r，θ）為極坐標(biāo)系，n表示第n層介質(zhì)，q（t）表示外加垂直機(jī)械載荷，φ（t）表示外加溫度載荷.

式（1）中第1，2個(gè)方程分別為第i層固相和流相介質(zhì)的動(dòng)量守恒方程，第3個(gè)方程為質(zhì)量守恒方程，第4個(gè)方程為能量守恒方程;uir為第i層固相的位移分量，u·ir為其相應(yīng)的速度分量，u··ir為其相應(yīng)的加速度分量;wir為第i層流相相對(duì)速度，w·ir為其相應(yīng)的相對(duì)加速度;θi=Θi-Θi0為第i層系統(tǒng)的變溫，其中Θi為第i層系統(tǒng)的絕對(duì)溫度，Θi0為第i層系統(tǒng)的初始絕對(duì)溫度，θ·i表示第i層溫度的變化速度; σSEir，σSEiθ分別表示第i層系統(tǒng)極徑r上的固相有效應(yīng)力和極角θ的固相有效應(yīng)力，pi為第i層流體的有效孔隙壓力，nαi（α=S，F(xiàn)）分別代表固相和流相的體積分?jǐn)?shù)，且nSi+nFi=1，ραi （α=S，F(xiàn)）分別代表固相和流相的宏觀質(zhì)量密度，它與微觀質(zhì)量密度ραRi之間有關(guān)系ραi=nαiραRi，在兩相不可壓縮的情況下ραRi為常數(shù);KSi=λSi+2μSi/3為體積模量，μSi，λSi為固相材料的Lame系數(shù)，αSi為固相材料的熱膨脹系數(shù);Siv=（niF）2γFRiκFi為固相與流相間的耦合系數(shù)，其中，γFRi為流相的有效比重，κFi為達(dá)西滲透系數(shù);βi為與流速有關(guān)的附加熱交換因子;ki為物相之間的熱傳導(dǎo)系數(shù);ρic=ρSicSi+ρFicFi為第i層系統(tǒng)的總比熱系數(shù)，其中，cαi （α=S，F(xiàn)）分別代表固相和流相的比熱系數(shù);ρir=ρSirSi+ρFirFi為第i層系統(tǒng)的總熱源強(qiáng)度，其中，rαi（α=S，F(xiàn)）分別代表固相和流相的熱源強(qiáng)度.

1.3 界面之間的連接性條件

對(duì)于分層不可壓流體飽和多孔熱彈性體而言，必須滿足界面之間的連接性條件.本研究的連接性條件為：1） 界面之間的固相位移分量必須相等;2） 界面之間的總應(yīng)力分量必須相等;3） 界面之間的流體有效孔隙壓必須相等;4） 界面之間的流體流量必須相等;5） 界面之間的變溫必須相等;6） 界面之間的熱流強(qiáng)度必須相等.

1.4 初始條件

2 DQEM和控制方程的微分求積（DQ）離散化

DQM是BELLMAN等[16-17]在20世紀(jì)70年代初提出的一種求解偏微分方程的數(shù)值算法.該算法的基本思路是用解區(qū)域中所有離散點(diǎn)處沿某個(gè)方向函數(shù)值的線性加權(quán)和作為該未知函數(shù)和它的各階導(dǎo)數(shù)在某一離散點(diǎn)的近似值，權(quán)系數(shù)只與解區(qū)域中所選擇的離散點(diǎn)和試函數(shù)有關(guān).因此，任何一個(gè)微分方程都可以轉(zhuǎn)化成相應(yīng)的代數(shù)方程.

DQM具有公式簡(jiǎn)單、使用方便、計(jì)算量少、精度高等優(yōu)點(diǎn).傳統(tǒng)的DQM對(duì)于求解具有非規(guī)則區(qū)域和間斷性條件的問(wèn)題，存在一些局限性.因此，研究人員構(gòu)建了微分求積單元法（DQEM），并取得了一系列的研究成果[18-19].DQEM基本步驟是：1） 將求解區(qū)域分割成若干個(gè)子區(qū)域或單元;2） 利用DQM，將各子區(qū)域的微分方程和邊界條件轉(zhuǎn)化為離散的代數(shù)方程組或者常微分方程組;3）將各單元的離散化方程連接起來(lái)，組成一個(gè)整體離散化的代數(shù)方程組或者常微分方程組;4） 采用適當(dāng)方法求解，從而得到各節(jié)點(diǎn)的未知量.

考慮在區(qū)域Ω={x0≤x≤a}內(nèi)的未知函數(shù)ψ（x），設(shè)沿x方向布置Nx個(gè)節(jié)點(diǎn)，根據(jù)DQM，函數(shù)ψ（x）在節(jié)點(diǎn)x=xξ（ζ=1，2，3，…，Nx）處對(duì)自變量x的n階導(dǎo)數(shù)可近似表示為：

其中，ψk=ψ（xk），為相應(yīng)節(jié)點(diǎn)的函數(shù)值，A（n）ζk為試函數(shù)對(duì)于x的n階偏導(dǎo)數(shù)的權(quán)系數(shù).本研究中，權(quán)系數(shù)由Lagrange插值多項(xiàng)式?jīng)Q定.

按照分層不可壓流體飽和多孔熱彈性體，耦合系統(tǒng)被劃分為n個(gè)單元，每個(gè)單元內(nèi)布置N個(gè)節(jié)點(diǎn)（圖2），節(jié)點(diǎn)坐標(biāo)由Chebyshev-Lobatto多項(xiàng)式的零點(diǎn)決定.

2.1 空間域內(nèi)控制方程的DQ離散化

2.2 邊界條件和對(duì)稱性條件的DQ離散化

2.4 對(duì)稱軸上奇異性條件的處理

3 分層不可壓軸對(duì)稱流體飽和多孔熱彈性體的動(dòng)力學(xué)特性

3.1 數(shù)值結(jié)果的驗(yàn)證

耦合系統(tǒng)被分別劃分為2，3，4個(gè)單元，每個(gè)單元內(nèi)布置Ni=11個(gè)節(jié)點(diǎn)（圖2），節(jié)點(diǎn)坐標(biāo)由Chebyshev-Lobatto多項(xiàng)式的零點(diǎn)來(lái)決定.

圖3給出了熱彈性體不同深度處的位移ur曲線，Δt=1 s.其中，ra表示飽和多孔介質(zhì)半徑.實(shí)線和圈劃線分別對(duì)應(yīng)n=2個(gè)單元和n=4個(gè)單元的情形下，利用DQEM得到的數(shù)值解;點(diǎn)線為利用文獻(xiàn)[14]中DQM模型得到的結(jié)果.

圖4給出了二階向后差分格式的步長(zhǎng)對(duì)熱彈性體不同深度處位移ur的影響，單元數(shù)n=3.實(shí)線和圈劃線分別對(duì)應(yīng)Δt=1 s和Δt=2 s的情形下，利用DQEM得到的數(shù)值解.點(diǎn)線為利用文獻(xiàn)[14]中DQM模型得到的結(jié)果.

通過(guò)計(jì)算發(fā)現(xiàn)：每個(gè)單元內(nèi)布置Ni=7個(gè)節(jié)點(diǎn)，能得到令人滿意的結(jié)果.為了節(jié)省篇幅，不在此給出示例.

從圖3和4中可以看到：采用兩種模型求得的解趨于一致，證明本方法具有較高的精度和收斂性.

3.2 分層不可壓軸對(duì)稱流體飽和多孔熱彈性體的動(dòng)力學(xué)特性

3.2.1 熱交換系數(shù)βi對(duì)熱彈性體動(dòng)力學(xué)特性的影響

從圖5中可以看到：ur隨著時(shí)間的增加趨于相同穩(wěn)定值;wr隨時(shí)間增加逐漸趨于0;p由初始值逐漸消散至0，且表面附近孔隙壓的消散速度快于內(nèi)部的消散速度;θ隨時(shí)間的增加逐漸上升并由表面向縱深處傳導(dǎo)和擴(kuò)散，最后達(dá)到等溫狀態(tài).同時(shí)，當(dāng)βi較小時(shí)，流、固兩相之間相互作用力較小，初始階段固相的熱體積膨脹效應(yīng)被抑制，熱彈性體沉降大.當(dāng)βi較大時(shí)，熱傳導(dǎo)過(guò)程快于機(jī)械載荷下的固結(jié)過(guò)程，初始階段表現(xiàn)為固相的熱體積膨脹效應(yīng)，而后固結(jié)作用才逐漸顯現(xiàn).

從圖6中可以看到，在分層不可壓軸對(duì)稱流體飽和多孔熱彈性體中，由于各個(gè)層中βi不同，p和θ在界面處不連續(xù).

3.2.2 體積分?jǐn)?shù)對(duì)熱彈性體動(dòng)力學(xué)特性的影響

圖7為不同的體積分?jǐn)?shù)ci對(duì)分層熱彈性體動(dòng)力學(xué)特性的影響.實(shí)線為分層（n=3）多孔熱彈性體的實(shí)驗(yàn)結(jié)果，此時(shí)c1=0.6，c2=0.8，c3=0.6;虛線為均勻多孔熱彈性體的實(shí)驗(yàn)結(jié)果，即c1=c2=c3=0.6.從圖7中可以看到，在分層不可壓軸對(duì)稱流體飽和多孔熱彈性體中，由于各個(gè)層中體積分?jǐn)?shù)不同，wr在界面處明顯不連續(xù).

4 結(jié) 論

在熱局部平衡條件下，基于PMT，研究了分層軸對(duì)稱流體飽和多孔熱彈性體在表面溫度載荷作用下的動(dòng)力學(xué)特性，提出了問(wèn)題的數(shù)學(xué)模型，采用DQEM、二階向后差分法及Newton-Raphson迭代法模擬問(wèn)題的數(shù)值結(jié)果.為了驗(yàn)證本方法的正確性，研究了不可壓流體飽和多孔彈性體的動(dòng)力固結(jié)問(wèn)題，并與

現(xiàn)有結(jié)果進(jìn)行比較，二者能良好地吻合，證明DQEM具備精度高、計(jì)算量小、數(shù)值穩(wěn)定等優(yōu)點(diǎn).研究和比較了一維分層軸對(duì)稱流體飽和多孔熱彈性體在表面受到溫度載荷時(shí)的動(dòng)力學(xué)特性，考察了材料參數(shù)對(duì)熱彈性體動(dòng)力學(xué)特性的影響.

參考文獻(xiàn)：

[1] BIOT M A.Theory of elasticity and consolidation for a porous anisotropic solid [J].Journal of Applied Physics，1955，26：182-185.

[2] CARTER J R，SAVVIDOU C.Consolidation around a spherical heat source [J].International Journal of Solids and Structures，1984，20：1079-1090.

[3] CUI Y J，SULTAN N，DELAGE P.A thermomechanical model for saturated clays [J].Canadian Geotechnical Journal，2000，37：607-620.

[4] 劉干斌，姚海林，楊洋，等.考慮熱-水-力耦合效應(yīng)多孔彈性地基的動(dòng)力響應(yīng) [J].巖土力學(xué)，2007，28（9）：1784-1795.

LIU G B，YAO H L，YANG Y，et al.Coupling thermo-hydro-mechanical dynamic response of a porous elastic medium [J].Rock and Soil Mechanics，2007，28（9）：1784-1795.

[5] 白冰.循環(huán)溫度荷載作用下飽和多孔介質(zhì)熱-水-力耦合響應(yīng) [J].工程力學(xué)，2007，24（5）：87-92.

BAI B.Thermo-hydro-mechanical response of saturated porous media under cycle thermal loading [J].Engineering Mechanics，2007，24（5）：87-92.

[6] 戴清晨，何錄武.熱局部非平衡條件下含柱形空洞橫觀各向同性飽和多孔介質(zhì)的熱應(yīng)力分析 [J].力學(xué)季刊，2014，35（1）：1-9.

DAI Q C，HE L W.Thermal stresses around a cylindrical hole in a transversely isotropic poroelastic medium considering local thermal non-equilibrium [J].Chinese Quarterly of Mechanics，2014，35（1）：1-9.

[7] de BOER R.Theoretical poroelasticity：a new approach [J].Chaos Solitons & Fractals，2005，25（4）：861-878.

[8] de BOER R，KOWALSKI S J.Thermodynamics of fluid-saturated porous media with a phase change [J].Acta Mechanica，1995，109（1/2/3/4）：167-189.

[9] de BOER R，EHLERS W，LIU Z.One-dimensional transient wave propagation in fluid-saturated incompressible porous media [J].Archive of Applied Mechanics，1993，63（1）：59-72.

[10] 劉占芳，姜乃斌，李思平.飽和多孔介質(zhì)一維瞬態(tài)波動(dòng)問(wèn)題的解析分析 [J].工程力學(xué)，2006，23（7）：19-24.

LIU Z F，JIANG N B，LI S P.An analysis on one-dimensional transient wave motion in saturated porous media [J].Engineering Mechanics，2006，23（7）：19-24.

[11] 李向約，李向維.飽和多孔介質(zhì)的熱固結(jié)理論 [J].固體力學(xué)學(xué)報(bào)，1990，11（4）：330-338.

LI X Y，LI X W.Theory of thermo-consolidation for saturated porous elastic media [J].Acta Mechanica Solida Sinica，1990，11（4）：330-338.

[12] HE L W，JIN Z H.A local thermal non-equilibrium poroelastic theory for fluid saturated [J].Journal of Thermal Stresses，2010，33：799-813.

[13] YANG X.Gurtin-type variational principles for dynamics of a non-local thermal equilibrium saturated porous medium [J].Acta Mechanica Solida Sinica，2005，18（1）：37-45.

[14] 朱媛媛，胡育佳，程昌鈞，等.基于DQM的空間軸對(duì)稱流體飽和多孔熱彈性柱體動(dòng)力學(xué)特性研究 [J].振動(dòng)與沖擊，2017，36（23）：83-91.

ZHU Y Y，HU Y J，CHENG C J，et al.The study on dynamic characteristics for a space-axisymmetrical fluid-saturated porous thermo-elastic cylinder based on DQM [J].Journal of Vibration and Shock，2017，36（23）：83-91.

[15] 嚴(yán)俊，魏迎奇，蔡紅，等.非飽和多孔介質(zhì)水-熱-力耦合數(shù)學(xué)模型研究 [J].水利學(xué)報(bào)，2014（增刊2）：152-160.

YAN J，WEI Y Q，CAI H，et al.A mathematical thermal hydraulic-mechanical coupling model for unsaturated porous media [J].Journal of Hydraulic Engineering，2014（Suppl.2）：152-160.

[16] BELLMAN R E，CASTI J.Differential quadrature and long term integration [J].Journal of Mathematical Analysis & Applications，1970，34（2）：235-238.

[17] BELLMAM R E，KASHEF B G，CASTI J.Differential quadrature：a technique for the rapid solution of nonlinear partial differential equations [J].Journal of Computational Physics，1972，10（1）：40-52.

[18] ZHU Y Y，HU Y J，CHENG C J.DQEM for analyzing dynamic characteristics of layered fluid-saturated porous elastic media [J].Acta Mechanica，2013，224（9）：1977-1998.

[19] 聶國(guó)雋，仲政.用微分求積法求解梁的彈塑性問(wèn)題 [J].工程力學(xué)，2005，22（1）：59-62.

NIE G J，ZHONG Z.Elasto-plastic analysis of beams by differential quadrature method [J].Engineering Mechanics，2005，22（1）：59-62.

（責(zé)任編輯：包震宇）