朱媛媛， 王笑梅

(上海師范大學(xué) 信息與機(jī)電工程學(xué)院，上海 200234)

0 引言

飽和多孔介質(zhì)熱-流-固耦合問(wèn)題不僅在傳統(tǒng)的土力學(xué)、水文學(xué)等經(jīng)典領(lǐng)域得到廣泛應(yīng)用，而且也廣泛應(yīng)用于核廢料污染物處置、地?zé)豳Y源、熱能貯存、人體關(guān)節(jié)軟骨組織分析等工程領(lǐng)域。因此，研究飽和多孔介質(zhì)的熱力學(xué)性能的理論、數(shù)值方法及其應(yīng)用有著重要的意義和廣泛的背景。早在1955年，Biot[1]建立了有關(guān)飽和多孔介質(zhì)的熱彈性理論，已在眾多工程領(lǐng)域成功地取得了許多研究成果[2-4]。在連續(xù)介質(zhì)混和物公理體系和體積分?jǐn)?shù)概念的基礎(chǔ)上，De Boer[5]于2005年建立了更加精確完整的多孔介質(zhì)熱彈性理論，一些微觀性質(zhì)可以直接用來(lái)描述宏觀性質(zhì)，容易反映非線(xiàn)性效應(yīng)。利用混合物理論，De Boer等[6]為一維流體飽和不可壓縮多孔介質(zhì)的動(dòng)力學(xué)固結(jié)問(wèn)題提供了解析解。Heider等[7]研究了波在飽和多孔半空間中的傳播。在De Boer等[8]提出的多孔介質(zhì)熱力學(xué)本構(gòu)關(guān)系的基礎(chǔ)上，He等[9]給出了一種在熱局部非平衡條件下的流體多孔介質(zhì)的數(shù)學(xué)模型，Yang[10]建立了該模型的相應(yīng)的Gurtin型變分原理，朱媛媛等[11]對(duì)空間軸對(duì)稱(chēng)流體飽和多孔熱彈性柱體動(dòng)力學(xué)特性進(jìn)行了分析。但目前有關(guān)飽和多孔介質(zhì)熱-流-固耦合系統(tǒng)非線(xiàn)性理論和應(yīng)用以及數(shù)值模擬的報(bào)道很少見(jiàn)。

為探索所提問(wèn)題理論和數(shù)值方法分析及其可行性，筆者基于多孔介質(zhì)混合理論(porous media theory， PMT)[5,8]，研究了在有限變形和熱局部非平衡條件下，不可壓流體飽和多孔熱彈性半平面在表面溫度載荷作用下的動(dòng)力學(xué)特性。在建立該問(wèn)題非線(xiàn)性數(shù)學(xué)模型的基礎(chǔ)上，通過(guò)采用微分求積法(differential quadrature method, DQM)在空間域內(nèi)建立離散問(wèn)題的非線(xiàn)性數(shù)學(xué)模型，用二階后向差分格式來(lái)處理時(shí)間導(dǎo)數(shù)，用Newton-Raphson法求解非線(xiàn)性代數(shù)方程組，從而可得到問(wèn)題的數(shù)值結(jié)果，分析半平面的動(dòng)力學(xué)特性。通過(guò)和De Boer等[6]的解析解的對(duì)比研究，表明本文方法是有效可靠的，具有計(jì)算量小、精度高等優(yōu)點(diǎn)。

1 數(shù)學(xué)模型

考察圖1所示平面直角坐標(biāo)系xoz中飽和不可壓多孔熱彈性半平面的動(dòng)力學(xué)問(wèn)題。假設(shè)作用于半平面表面的載荷或者溫度具有某種對(duì)稱(chēng)性，取oz軸為對(duì)稱(chēng)軸，由表面指向半平面內(nèi)部的方向?yàn)檎较?圖1)。

圖1 半平面問(wèn)題的數(shù)學(xué)模型Figure 1 Mathematical model of half-plane problem

設(shè)平面處于熱局部非平衡狀態(tài)(固相介質(zhì)和液相介質(zhì)的變溫不同)，并假設(shè)平面介質(zhì)為各向同性線(xiàn)性熱彈性多孔固相骨架和理想流體，根據(jù)PMT，在不可壓和小流速的假設(shè)下，半平面的控制微分方程如下：

(1a)

(1b)

(1c)

(1d)

(1e)

Θ0βS(wx,x+wz,z)-eθθFS+

(1f)

(1g)

在有限變形的情況下，對(duì)于各向同性線(xiàn)性熱彈性的固相材料，幾何關(guān)系和本構(gòu)關(guān)系為：

(2)

(3)

同時(shí)，固相材料的有效應(yīng)力與總應(yīng)力之間的關(guān)系為：

(4)

假設(shè)流體飽和多孔彈性半平面的上表面處于理想排水狀態(tài)，同時(shí)承受變溫φ(x,t)和垂直載荷q(x,t)的作用；底部為不排水剛性絕熱底面；側(cè)面不排水絕熱且在x方向位移被約束 (圖1)。設(shè)φ(x,t)和q(x,t)關(guān)于oz軸對(duì)稱(chēng)，可以在區(qū)域0≤x≤L0,0≤z≤H0內(nèi)考慮問(wèn)題，這時(shí)在熱局部非平衡狀態(tài)下，有如下邊界條件和對(duì)稱(chēng)條件。

(5)

(6)

(7)

(8)

假設(shè)t≤0時(shí)半平面處于靜止和熱力平衡狀態(tài)，在熱局部非平衡狀態(tài)下初始條件為：

(9)

因此，式(1)～(8)構(gòu)成了在熱局部非平衡條件下流體飽和多孔熱彈性半平面動(dòng)力學(xué)問(wèn)題的一個(gè)數(shù)學(xué)模型。當(dāng)固相和液相介質(zhì)的變溫相同時(shí)，式(1)～(8)退化為熱局部平衡條件下的數(shù)學(xué)模型。

不難看到，很難獲得上述動(dòng)力學(xué)問(wèn)題的解析或半解析解。在本文中，對(duì)于熱局部非平衡條件下的非線(xiàn)性數(shù)學(xué)模型在空間域中使用DQM離散，時(shí)間導(dǎo)數(shù)采用二階后向差分格式處理，離散后的非線(xiàn)性代數(shù)方程組采用Newton-Raphson法迭代求解，可得到問(wèn)題的數(shù)值模擬結(jié)果。

2 求解方法

1970年，Bellman等[12-13]提出了DQM。DQM具有公式簡(jiǎn)單、精度高、計(jì)算量少等優(yōu)點(diǎn)，已有很多成功的應(yīng)用[14-15]。DQM的基本思想是將函數(shù)的偏導(dǎo)數(shù)近似為離散點(diǎn)處相應(yīng)函數(shù)值的加權(quán)和。由于DQM權(quán)函數(shù)的選取與特定問(wèn)題無(wú)關(guān)，微分方程能DQ離散為相應(yīng)的代數(shù)方程求解。

設(shè)半平面占有區(qū)域Ω={〈x,z〉|0≤x≤L0,0≤z≤H0}，分別在x、z軸方向布置Nx×Nz個(gè)節(jié)點(diǎn)，節(jié)點(diǎn)坐標(biāo)為Chebyshev-Lobatto多項(xiàng)式的零點(diǎn)。根據(jù)DQM，未知函數(shù)ψ(x,z)對(duì)x的n階偏導(dǎo)數(shù)可近似表示為：

(10)

利用式(10)，在空間域內(nèi)對(duì)數(shù)學(xué)模型(1)～(8)DQM離散化后，產(chǎn)生了關(guān)于時(shí)間t的微分代數(shù)方程組。使用二階后向差分格式來(lái)處理關(guān)于微分代數(shù)方程組中的時(shí)間導(dǎo)數(shù)：

(11)

式中：Δt是時(shí)間步長(zhǎng)；ψ(tn,x)是在時(shí)刻t=tn的函數(shù)值。

最后利用DQ離散的初始條件(9)，對(duì)離散化方程組進(jìn)行Newton-Raphson迭代求解，從而可得到問(wèn)題的數(shù)值結(jié)果。

3 分析與討論

3.1 驗(yàn)證方法

忽略體積力的影響，De Boer等[6]分析了幾何小變形和常溫條件下不可壓飽和多孔彈性介質(zhì)的一維動(dòng)力固結(jié)問(wèn)題，并給出了在PMT下唯一的解析解。為了驗(yàn)證本文方法的正確性和有效性，本文在第1節(jié)給定的數(shù)學(xué)模型中設(shè)平面表面變溫為φ(x,t)=0，從而可以認(rèn)為在熱局部非平衡條件下的解近似等于常溫狀態(tài)時(shí)的解，然后用一個(gè)高H0=10 m，寬L0=10 m的流體飽和多孔彈性體來(lái)模擬De Boer等[6]所考慮的類(lèi)似問(wèn)題，并與解析解進(jìn)行比較。在驗(yàn)證中，分別考慮了階梯載荷q(x,t)=q0h(t)和周期載荷q(x,t)=q0[1-cos(ωt)]的影響，其中，q0=3 kN/m2，ω=0.4 s-1，h(t)為Heaviside函數(shù)。本文計(jì)算中選取的材料是廢料存儲(chǔ)中使用的黏土材料[2]，材料參數(shù)見(jiàn)表1。

圖2和圖3分別給出了在階梯載荷和周期載荷下，對(duì)稱(chēng)軸處不同深度的位移uz(沉降)，其中，實(shí)線(xiàn)和虛線(xiàn)分別對(duì)應(yīng)了在7×7和9×9布點(diǎn)下利用本文方法得到的數(shù)值解，坐標(biāo)點(diǎn)為文獻(xiàn)[6]中的解析解。時(shí)間步長(zhǎng)Δt=1 s。從圖中看到，本文提供的數(shù)值結(jié)果與解析結(jié)果吻合良好，當(dāng)布置結(jié)點(diǎn)數(shù)為7×7個(gè)時(shí)，便能得到令人滿(mǎn)意的結(jié)果。說(shuō)明用本文方法來(lái)計(jì)算流體飽和多孔彈性介質(zhì)的動(dòng)力學(xué)特性是有效的。

表1 材料參數(shù)Table 1 Parameters of material

圖2 位移時(shí)程曲線(xiàn)圖(階梯載荷)Figure 2 Time-history curves of displacement (step loading)

圖3 位移時(shí)程曲線(xiàn)圖(周期載荷)Figure 3 Time-history curves of displacement (cycle loading)

3.2 耦合系統(tǒng)的動(dòng)力學(xué)特性分析

在熱局部非平衡條件下，考察高(深)H0=10 m，寬L0=10 m的半平面域，邊界條件由式(5)給定?？疾毂砻骐A梯溫度載荷φ(x,t)=φ0h(t)的影響，其中φ0=10 ℃，h(t)是Heaviside函數(shù)。計(jì)算中布點(diǎn)數(shù)為7×7，時(shí)間步長(zhǎng)Δt=1 s，表面垂直外載q(x,t)=3 kN/m2，材料參數(shù)見(jiàn)表1。

3.2.1 達(dá)西滲透系數(shù)κF的影響

圖4顯示了在熱局部非平衡條件下，達(dá)西滲透系數(shù)κF對(duì)半平面對(duì)稱(chēng)軸(x=0)處未知量uz、wz、p及θS和θF的影響。圖4(a)表明隨著時(shí)間的逐漸增加，沉降uz最終達(dá)到相同的穩(wěn)定狀態(tài)。κF較大時(shí)，初始階段固結(jié)速度比熱傳導(dǎo)快，導(dǎo)致熱體積膨脹不明顯，所以半平面的沉降uz大。圖4(b)表明相對(duì)流速wz的變化趨勢(shì)。當(dāng)系統(tǒng)趨于穩(wěn)定時(shí)，相對(duì)流速wz趨于零。但在初始階段，平面內(nèi)部的流速小于上表面附近，并且κF較大時(shí)流速的峰值較大。圖4(c) 表明孔隙壓p隨著時(shí)間由初始孔壓逐漸至零，半平面上表面附近孔隙壓的面內(nèi)部快，同時(shí)孔隙壓的消散速度隨κF增大而增大。消散速度較平圖4(d) 表明固相和液相的溫度從上表面?zhèn)鞑サ桨肫矫嫔钐?，并且逐漸到達(dá)至給定的表面溫度，最終達(dá)到等溫狀態(tài)。由于在熱局部非平衡條件下，流固兩相之間的溫度是不同的，在給定的參數(shù)下，固相溫度高于流相溫度，并且在初始階段溫差明顯。同時(shí)，κF對(duì)變溫θS和θF的影響很小。

圖4 達(dá)西滲透系數(shù)對(duì)半平面時(shí)-程曲線(xiàn)的影響Figure 4 Effect of Darcy permeability coefficient on time-history curves of half-plane

3.2.2 熱能交換系數(shù)eθ的影響

圖5給出了在熱局部非平衡條件下，不同物相之間熱能交換系數(shù)eθ對(duì)半平面對(duì)稱(chēng)軸(x=0)處未知量uz及θS和θF的影響。同時(shí)，圖5也給出了在熱局部平衡條件下系統(tǒng)的時(shí)程曲線(xiàn)。在熱局部非平衡條件下，隨著eθ增大，固相溫度降低并逐漸接近熱局部平衡條件下系統(tǒng)的溫度，而流相溫度升高并逐漸接近熱局部平衡條件下系統(tǒng)的溫度，因此流固兩相溫差逐漸變小。另外，隨著eθ增大，固相溫度降低，在初始階段系統(tǒng)熱傳導(dǎo)過(guò)程小于固結(jié)作用，沉降變大。

與此同時(shí)，筆者也考慮了熱交換系數(shù)βS、βF和熱傳導(dǎo)系數(shù)kSS、kFF、kSF對(duì)半平面的影響。限于篇幅，僅列出結(jié)論如下：當(dāng)熱交換系數(shù)βS、βF較小時(shí)，流相和固相之間的耦合作用減弱，初始階段半平面的沉降uz較大。當(dāng)熱傳導(dǎo)系數(shù)kSS、kFF、kSF較大時(shí)，初始階段時(shí)熱傳導(dǎo)比固結(jié)過(guò)程快，半平面的沉降uz較小，溫度場(chǎng)達(dá)到等溫狀態(tài)所需的時(shí)間更短。

3.2.3 非線(xiàn)性的影響

圖6給出了在熱局部非平衡條件下，對(duì)稱(chēng)軸不同深度處沉降uz的幾何非線(xiàn)性結(jié)果與相應(yīng)的線(xiàn)性解的比較。可以看到，當(dāng)載荷較小時(shí)(q(x,t)=3 kN/m2)，非線(xiàn)性解和線(xiàn)性解趨于一致。當(dāng)載荷較大時(shí)(q(x,t)=30 kN/m2)，非線(xiàn)性解略大于線(xiàn)性解。

圖5 熱能交換系數(shù)對(duì)半平面時(shí)-程曲線(xiàn)的影響Figure 5 Effect of energy exchange coefficient on time-history curves of half-plane

圖6 非線(xiàn)性解與線(xiàn)性解的對(duì)比Figure 6 Comparison between nonlinear solutions and linear solutions

4 結(jié)論

在幾何非線(xiàn)性條件和熱局部非平衡條件下，對(duì)流體飽和不可壓多孔熱彈性半平面在表面溫度載荷作用下的動(dòng)力學(xué)特性進(jìn)行研究。首先基于PMT并考慮幾何非線(xiàn)性和熱局部非平衡的影響，給出了問(wèn)題的數(shù)學(xué)模型；其次提出了一個(gè)系統(tǒng)的綜合數(shù)值計(jì)算方法來(lái)模擬問(wèn)題的數(shù)值結(jié)果，通過(guò)DQM和二階后向差分格式分別在空間域和時(shí)間域離散數(shù)學(xué)模型，利用Newton-Raphson法求解非線(xiàn)性代數(shù)方程組，從而可得到問(wèn)題的數(shù)值結(jié)果，分析半平面的動(dòng)力學(xué)特性。

和現(xiàn)有解析解的對(duì)比研究驗(yàn)證了本文方法是正確可靠的，具有精度高、計(jì)算量小、數(shù)值穩(wěn)定等優(yōu)點(diǎn)。在此基礎(chǔ)上，研究了在表面溫度載荷作用下半平面的熱力學(xué)特性，考察了系統(tǒng)的材料參數(shù)以及幾何非線(xiàn)性對(duì)半平面動(dòng)力學(xué)特性的影響，獲得了一些有益的結(jié)論。