基于DQEM的分層流體飽和熱彈性多孔介質軸對稱問題的動力響應分析

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

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

關鍵詞： 流體飽和熱彈性多孔介質; 多孔介質理論（PMT）; 微分求積單元法（DQEM）; 動力學響應

中圖分類號： TU 311文獻標志碼： A文章編號： 1000-5137（2019）02-0141-10

0 引 言

飽和多孔介質熱-流-固耦合系統的研究不僅在土力學、水文學等經典應用領域發揮重要作用，也已成為許多新興學科和應用技術發展的關鍵，其相關理論和數值方法的研究工作具有重要的理論意義和廣泛的應用背景.

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

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

1 問題的數學描述

1.1 控制微分方程

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

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

1.3 界面之間的連接性條件

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

1.4 初始條件

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

DQM是BELLMAN等[16-17]在20世紀70年代初提出的一種求解偏微分方程的數值算法.該算法的基本思路是用解區域中所有離散點處沿某個方向函數值的線性加權和作為該未知函數和它的各階導數在某一離散點的近似值，權系數只與解區域中所選擇的離散點和試函數有關.因此，任何一個微分方程都可以轉化成相應的代數方程.

DQM具有公式簡單、使用方便、計算量少、精度高等優點.傳統的DQM對于求解具有非規則區域和間斷性條件的問題，存在一些局限性.因此，研究人員構建了微分求積單元法（DQEM），并取得了一系列的研究成果[18-19].DQEM基本步驟是：1） 將求解區域分割成若干個子區域或單元;2） 利用DQM，將各子區域的微分方程和邊界條件轉化為離散的代數方程組或者常微分方程組;3）將各單元的離散化方程連接起來，組成一個整體離散化的代數方程組或者常微分方程組;4） 采用適當方法求解，從而得到各節點的未知量.

考慮在區域Ω={x0≤x≤a}內的未知函數ψ（x），設沿x方向布置Nx個節點，根據DQM，函數ψ（x）在節點x=xξ（ζ=1，2，3，…，Nx）處對自變量x的n階導數可近似表示為：

其中，ψk=ψ（xk），為相應節點的函數值，A（n）ζk為試函數對于x的n階偏導數的權系數.本研究中，權系數由Lagrange插值多項式決定.

按照分層不可壓流體飽和多孔熱彈性體，耦合系統被劃分為n個單元，每個單元內布置N個節點（圖2），節點坐標由Chebyshev-Lobatto多項式的零點決定.

2.1 空間域內控制方程的DQ離散化

2.2 邊界條件和對稱性條件的DQ離散化

2.4 對稱軸上奇異性條件的處理

3 分層不可壓軸對稱流體飽和多孔熱彈性體的動力學特性

3.1 數值結果的驗證

耦合系統被分別劃分為2，3，4個單元，每個單元內布置Ni=11個節點（圖2），節點坐標由Chebyshev-Lobatto多項式的零點來決定.

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

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

通過計算發現：每個單元內布置Ni=7個節點，能得到令人滿意的結果.為了節省篇幅，不在此給出示例.

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

3.2 分層不可壓軸對稱流體飽和多孔熱彈性體的動力學特性

3.2.1 熱交換系數βi對熱彈性體動力學特性的影響

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

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

3.2.2 體積分數對熱彈性體動力學特性的影響

圖7為不同的體積分數ci對分層熱彈性體動力學特性的影響.實線為分層（n=3）多孔熱彈性體的實驗結果，此時c1=0.6，c2=0.8，c3=0.6;虛線為均勻多孔熱彈性體的實驗結果，即c1=c2=c3=0.6.從圖7中可以看到，在分層不可壓軸對稱流體飽和多孔熱彈性體中，由于各個層中體積分數不同，wr在界面處明顯不連續.

4 結 論

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

現有結果進行比較，二者能良好地吻合，證明DQEM具備精度高、計算量小、數值穩定等優點.研究和比較了一維分層軸對稱流體飽和多孔熱彈性體在表面受到溫度載荷時的動力學特性，考察了材料參數對熱彈性體動力學特性的影響.

參考文獻：

[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] 劉干斌，姚海林，楊洋，等.考慮熱-水-力耦合效應多孔彈性地基的動力響應 [J].巖土力學，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] 白冰.循環溫度荷載作用下飽和多孔介質熱-水-力耦合響應 [J].工程力學，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] 戴清晨，何錄武.熱局部非平衡條件下含柱形空洞橫觀各向同性飽和多孔介質的熱應力分析 [J].力學季刊，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] 劉占芳，姜乃斌，李思平.飽和多孔介質一維瞬態波動問題的解析分析 [J].工程力學，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] 李向約，李向維.飽和多孔介質的熱固結理論 [J].固體力學學報，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的空間軸對稱流體飽和多孔熱彈性柱體動力學特性研究 [J].振動與沖擊，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] 嚴俊，魏迎奇，蔡紅，等.非飽和多孔介質水-熱-力耦合數學模型研究 [J].水利學報，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] 聶國雋，仲政.用微分求積法求解梁的彈塑性問題 [J].工程力學，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.

（責任編輯：包震宇）