飽和多孔彈性桿熱傳導(dǎo)的廣義多辛方法及其數(shù)值實(shí)現(xiàn)

?

飽和多孔彈性桿熱傳導(dǎo)的廣義多辛方法及其數(shù)值實(shí)現(xiàn)

劉雪梅1，2，鄧子辰1，胡偉鵬1

(1．西北工業(yè)大學(xué)力學(xué)與土木建筑學(xué)院，陜西西安710072; 2．長安大學(xué)工程力學(xué)系，陜西西安710064)

摘要:首先根據(jù)多孔介質(zhì)理論，利用飽和多孔介質(zhì)的能量方程和本構(gòu)關(guān)系，推導(dǎo)出飽和多孔彈性桿局部熱平衡的熱傳導(dǎo)方程;繼而引入正交變量，將熱傳導(dǎo)方程導(dǎo)入Hamilton系統(tǒng)，得到飽和多孔彈性桿熱傳導(dǎo)方程的廣義多辛形式和多種局部守恒律形式;接著采用中點(diǎn)離散方法對(duì)熱傳導(dǎo)方程的廣義多辛形式進(jìn)行數(shù)值離散;最后利用計(jì)算機(jī)數(shù)值實(shí)現(xiàn)了飽和多孔彈性桿的熱傳導(dǎo)過程，并且討論了參數(shù)取值的不同對(duì)熱傳導(dǎo)過程的影響，同時(shí)在數(shù)值模擬過程中記錄了廣義多辛格式的局部動(dòng)量誤差。研究結(jié)果表明，構(gòu)造的廣義多辛方法能夠很好地模擬系統(tǒng)的熱傳導(dǎo)過程和耗散效應(yīng)，同時(shí)也可長時(shí)間保持系統(tǒng)的固有幾何性質(zhì)。

關(guān)鍵詞:多孔介質(zhì);廣義多辛;耗散;熱傳導(dǎo)

飽和多孔介質(zhì)通常是指孔隙中充滿液體的多孔連通介質(zhì)，其相關(guān)理論在地震工程學(xué)、土動(dòng)力學(xué)、地球物理學(xué)、生物力學(xué)等領(lǐng)域有著廣泛應(yīng)用。目前，研究多孔介質(zhì)宏觀力學(xué)行為的主要理論有: Biot理論、多孔介質(zhì)理論和混合物理論。在Biot理論基礎(chǔ)上，徐曾和等研究了含單一圓井的飽和多孔地層中，以定流量方式開采時(shí)的流固耦合問題［1］; de Boer等基于多孔介質(zhì)理論，研究了流固兩相微觀不可壓飽和多孔介質(zhì)穩(wěn)態(tài)振動(dòng)的基本解［2］;秦冰等以混合物理論為基礎(chǔ)，將非飽和土視為由土骨架、液態(tài)水、水蒸氣、干燥氣體及溶解氣體共5種組分構(gòu)成的混合物，系統(tǒng)研究了非飽和土的熱-水-力多場(chǎng)耦合問題［3］。然而，有關(guān)各類飽和多孔結(jié)構(gòu)熱傳導(dǎo)過程的數(shù)值算法研究較少，尤其是熱傳導(dǎo)方程的求解過程中如何體現(xiàn)由于熱量傳遞引起的耗散效應(yīng)以及如何保持系統(tǒng)的局部能量守恒律和局部動(dòng)量守恒律等局部幾何性質(zhì)更是數(shù)值算法的難點(diǎn)。

1984年是應(yīng)用數(shù)學(xué)和計(jì)算力學(xué)研究中具有重要意義的一年，也是辛幾何算法開創(chuàng)性發(fā)展的一年，我國數(shù)學(xué)家馮康于這一年首次系統(tǒng)地提出了哈密頓方程和哈密頓算法，并指出了從辛幾何內(nèi)部系統(tǒng)構(gòu)成算法并研究其性質(zhì)的途徑，從而開創(chuàng)了辛幾何算法的新領(lǐng)域［4］。自此以后，辛幾何算法得到了長足的發(fā)展，學(xué)術(shù)界將純理論的辛幾何和現(xiàn)代科學(xué)工程計(jì)算有機(jī)地結(jié)合起來，使其廣泛地應(yīng)用于天體動(dòng)力學(xué)和分子動(dòng)力學(xué)等諸多領(lǐng)域［5］。在此基礎(chǔ)上，Bridge教授針對(duì)無窮維Hamilton系統(tǒng)提出了多辛算法［6］，用以在數(shù)值求解過程中保持系統(tǒng)的局部幾何性質(zhì)。然而，實(shí)際力學(xué)系統(tǒng)中耗散效應(yīng)的存在卻成了多辛算法的“瓶頸”，近年來，廣義多辛算法理論的提出［7-9］，有效地解決了這一"瓶頸"問題，將保結(jié)構(gòu)算法理論體系拓寬至耗散系統(tǒng)。

多孔介質(zhì)傳熱現(xiàn)象遍及于工農(nóng)業(yè)生產(chǎn)的各個(gè)領(lǐng)域，例如:埋地電纜和直流接地極的熱耗散;地源熱泵和地冷空調(diào)中換熱器埋管;高溫原件的多孔冷卻;生物多孔介質(zhì)的傳熱問題等等。因此，多孔介質(zhì)傳熱理論和應(yīng)用研究是一個(gè)具有重要學(xué)術(shù)和應(yīng)用價(jià)值的研究方向，而作為多孔介質(zhì)一維傳熱問題中最基礎(chǔ)的飽和多孔彈性桿的熱傳導(dǎo)就更加重要，它的研究將為多孔介質(zhì)桿系結(jié)構(gòu)以及植物根莖等多種多孔介質(zhì)結(jié)構(gòu)的傳熱傳質(zhì)問題提供有力的依據(jù)，尤其是其數(shù)值求解的應(yīng)用將更加廣泛。正是基于這一點(diǎn)，本文試圖通過廣義多辛保結(jié)構(gòu)算法來尋求新的研究飽和多孔彈性桿熱傳導(dǎo)過程的數(shù)值方法，同時(shí)也為了揭示其熱傳導(dǎo)過程中的耗散效應(yīng)以及系統(tǒng)的廣義多辛守恒律等多種局部幾何性質(zhì)，從而為多孔介質(zhì)傳熱問題的數(shù)值求解提供新的途徑。

本文基于飽和多孔介質(zhì)理論，首先建立飽和多孔彈性桿熱傳導(dǎo)的一維數(shù)學(xué)模型。其次，通過正則變換，構(gòu)造飽和多孔彈性桿熱傳導(dǎo)方程的廣義多辛形式，并研究了熱量傳遞的耗散效應(yīng)在廣義多辛形式中的數(shù)學(xué)表述。隨后，采用中點(diǎn)離散方法離散廣義多辛形式，得到其中點(diǎn)Box廣義多辛離散格式。最后，利用中點(diǎn)Box廣義多辛離散格式數(shù)值模擬飽和多孔彈性桿的熱傳導(dǎo)過程，并將數(shù)值解與精確解進(jìn)行比對(duì)，同時(shí)研究了熱傳導(dǎo)過程中的耗散效應(yīng)。本文的廣義多辛分析將完善保結(jié)構(gòu)算法理論，也為多孔介質(zhì)耗散系統(tǒng)的數(shù)值求解提供新的途徑。

1　飽和多孔彈性桿的熱傳導(dǎo)方程

為了研究飽和多孔彈性桿的傳熱問題，本文建立如圖1的模型:設(shè)長為L、橫截面高度為H的流體飽和多孔彈性桿由不相溶的微觀不可壓流體和微觀不可壓彈性多孔固相骨架組成，其側(cè)表面不透水。根據(jù)多孔介質(zhì)理論，每一相物質(zhì)被理想化地分布在整個(gè)區(qū)域中，并擁有各自獨(dú)立的運(yùn)動(dòng)，記物質(zhì)粒子的位移為:

圖1　LH的飽和多孔彈性桿

式中，i =S、F分別表示固相和流相，V為多孔介質(zhì)初始時(shí)的空間區(qū)域。忽略兩相間的質(zhì)量交換，則流固兩相的能量方程可表示為［10］:

式中，ρi為兩相的宏觀質(zhì)量密度;εi、σi、qi和ri分別為兩相的內(nèi)能、宏觀應(yīng)力張量、熱流矢量和內(nèi)部熱源;i和i分別為兩相間的相互作用力和能量交換量，滿足+= 0+= 0; ()i=+ grad(…)·i表示隨物相的物質(zhì)時(shí)間導(dǎo)數(shù)。

對(duì)于各向同性線彈性多孔固相骨架和無粘性流體，在兩相不可壓和小變形的假定下，并忽略由于流固兩相的相對(duì)運(yùn)動(dòng)引起的附加熱對(duì)流和附加熱傳導(dǎo)，可得流固兩相的宏觀應(yīng)變張量和宏觀應(yīng)變率張量分別為:

本構(gòu)方程可表示為［10］:

式中，上標(biāo)T表示轉(zhuǎn)置，p和Sν分別表示有效孔隙水壓力和與流體滲流系數(shù)有關(guān)的兩相相互作用耦合系數(shù)，eΘ為非局部熱平衡引起的兩相間能量交換系數(shù)，μS和λS為固相宏觀拉梅常數(shù)，kii和θi(i = S，F(xiàn))分別為熱傳導(dǎo)系數(shù)和流固兩相的溫度，wF=F－S為孔隙流體相對(duì)固相骨架的速度。

以下考慮局部熱平衡情況，即在每一個(gè)空間點(diǎn)上，固相骨架和孔隙流體具有相同的溫度，其溫度可表示為θS=θF=θ，將(3)式和(4)式代入能量方程(2)中，并采用如圖1所示的一維空間坐標(biāo)軸，同時(shí)忽略非線性項(xiàng)對(duì)熱傳導(dǎo)的影響［10］，可得飽和多孔彈性桿局部熱平衡時(shí)的熱傳導(dǎo)方程:

式中，θ和ci(i = S，F(xiàn))分別表示局部熱平衡溫度和流固兩相的比熱容，再引入如下無量綱量和常量:

2　飽和多孔彈性桿熱傳導(dǎo)方程的廣義多辛形式

對(duì)于飽和多孔彈性桿熱傳導(dǎo)方程(6)，引入正交變量φ=?xθ，上述二階偏微分方程(6)可轉(zhuǎn)化為以下一階耦合Hamilton方程組形式:

M?tz + K?xz =▽zS(z)，z =［θ，φ］T∈R2(7)

從(8)式容易看到，系數(shù)矩陣M并不是反對(duì)稱的，因此，按照廣義多辛算法的概念［8］，(7)式為廣義多辛形式。這是由于本文涉及的飽和多孔彈性桿的熱傳導(dǎo)問題中存在著熱量傳遞，從而使系數(shù)矩陣M并不是反對(duì)稱的，也正是由于非完全反對(duì)稱矩陣的存在，上述廣義多辛形式(7)式就不能滿足多辛算法的多辛守恒律和局部動(dòng)量守恒律。

基于多辛理論，由多辛守恒律概念可給出廣義多辛形式(7)式的廣義多辛守恒律為［8］:

同時(shí)，由局部動(dòng)量守恒律誤差的概念，亦可給出廣義多辛形式(7)式的局部動(dòng)量守恒律誤差為:

這里需要指出，由于系數(shù)矩陣K是嚴(yán)格反對(duì)稱的，廣義多辛形式(7)式嚴(yán)格滿足局部能量守恒律。

3　飽和多孔彈性桿熱傳導(dǎo)方程的廣義多辛離散

本文選用中點(diǎn)差分離散方法分別在空間方向和時(shí)間方向?qū)V義多辛形式(7)式進(jìn)行差分離散，并分別選取時(shí)間步長Δt和空間步長Δx對(duì)計(jì)算區(qū)域均勻劃分，用zji表示狀態(tài)變量z在(xi，tj)網(wǎng)格點(diǎn)上的近似值，則可得時(shí)間和空間2個(gè)方向上的中點(diǎn)Box離散格式:

式中

將系數(shù)矩陣(8)和狀態(tài)變量z =［θ，φ］T以及Hamilton函數(shù)代入上述離散格式(11)，便可得到廣義多辛偏微分方程(7)的等價(jià)形式:

對(duì)于上述已得到的廣義多辛守恒律(9)式，亦可通過有限差分運(yùn)算得到其離散格式。這里采用向前差分運(yùn)算，得到廣義多辛局部動(dòng)量誤差的離散形式［7］:

為了得到局部動(dòng)量誤差隨時(shí)間變化的平面圖，本文取廣義多辛局部動(dòng)量誤差在第j時(shí)間步的絕對(duì)值［7］:

4　數(shù)值實(shí)驗(yàn)

這一節(jié)中，本文將利用以下飽和多孔彈性懸臂桿的數(shù)值算例來探討和驗(yàn)證前述已構(gòu)造出的廣義多辛離散格式的精確性、有效性和數(shù)值穩(wěn)定性等性能。

在廣義多辛離散格式(12)中，考慮如下邊界條件和初始條件:

選取計(jì)算空間步長Δx = 0. 01，時(shí)間步長Δt =0. 5，λ2= 0. 000 1，分別取參數(shù)C = 0. 043和C = 0. 02，在t∈［0，200］內(nèi)，得到局部熱平衡時(shí)飽和多孔彈性桿不同時(shí)刻的溫度隨空間坐標(biāo)的變化圖(見圖2和圖3)。為了驗(yàn)證廣義多辛離散格式(12)的精確性，采用分離變量法，可得無量綱偏微分方程(6)的定解問題(15)～(16)的精確解:

式中

因此，亦可得局部熱平衡時(shí)飽和多孔彈性桿不同時(shí)刻溫度場(chǎng)的精確解隨空間坐標(biāo)的變化圖(圖2和圖3)。

圖2　C=0. 043時(shí)不同時(shí)刻溫度的數(shù)值解(點(diǎn))與精確解(曲線)的比較

圖3　C=0. 02時(shí)不同時(shí)刻溫度的數(shù)值解(點(diǎn))與精確解(曲線)的比較

圖2反映了參數(shù)C = 0. 043時(shí)飽和多孔彈性桿局部熱平衡的熱傳導(dǎo)過程，這一過程中，由于桿x = 1端與外界絕熱，初始溫度從桿x = 1端向桿x = 0端逐漸降低，因此隨著時(shí)間的推移，熱量將從桿x = 1端向桿x = 0端廣泛傳遞，從初始時(shí)刻到t = 200的過程中，桿x = 1端的溫度從0. 5降低至0. 16;熱傳導(dǎo)過程中，由于桿內(nèi)部無熱源，桿x = 0端溫度始終為零，桿將由x = 0端持續(xù)向外界傳熱，桿內(nèi)各點(diǎn)的溫度均逐漸下降。與圖2類似，圖3反映了參數(shù)C = 0. 02時(shí)飽和多孔彈性桿局部熱平衡的熱傳導(dǎo)過程，這一過程中明顯可以看到，相同時(shí)刻桿內(nèi)各點(diǎn)的溫度均比C = 0. 043時(shí)低，t = 200時(shí)，桿內(nèi)各點(diǎn)溫度均已接近零;由此可以說明，參數(shù)C的取值越小，熱傳導(dǎo)過程將越快;這種現(xiàn)象是因?yàn)閰?shù)C取值的減小，相當(dāng)于飽和多孔彈性桿熱傳導(dǎo)方程的導(dǎo)熱系數(shù)增大，則必然導(dǎo)致熱傳導(dǎo)過程加快。此外，由圖2和圖3均可看到，熱傳導(dǎo)過程中溫度場(chǎng)分布的數(shù)值解與精確解非常吻合。以上這些結(jié)果說明本文采用的廣義多辛方法具有很好的有效性和精確性。

為了進(jìn)一步驗(yàn)證廣義多辛格式(12)的局部保結(jié)構(gòu)性，本文將無量綱模擬時(shí)間延長至500，利用(14)式，將第j時(shí)間步中誤差絕對(duì)值最大的空間點(diǎn)的數(shù)值誤差作為該時(shí)間步的局部動(dòng)量數(shù)值誤差，并取參數(shù)C足夠小(其中C分別取0. 01和0. 043) ;繼而在計(jì)算機(jī)模擬過程中記錄每一時(shí)間步的局部動(dòng)量誤差，從而得到廣義多辛局部動(dòng)量數(shù)值誤差隨時(shí)間的變化曲線圖(見圖4)。

從局部動(dòng)量數(shù)值誤差的結(jié)果可以看出:整個(gè)模擬過程中，局部動(dòng)量誤差均在10－5數(shù)量級(jí)以下，這說明本文構(gòu)造的廣義多辛算法能夠很好地模擬系統(tǒng)的熱量傳遞耗散效應(yīng)，并且保持了系統(tǒng)的固有幾何性質(zhì)。從圖4可以看到，參數(shù)C取0. 043時(shí)，無量綱時(shí)間從t = 450之后，局部動(dòng)量誤差為零;然而參數(shù)C取0. 01時(shí)，無量綱時(shí)間從t = 100之后，局部動(dòng)量誤差就已經(jīng)為零。這就是說，參數(shù)C的取值越小，局部動(dòng)量誤差將越早等于零，這也就表示，參數(shù)C的取值越小，越能體現(xiàn)本文構(gòu)造的廣義多辛形式的局部保結(jié)構(gòu)性。

圖4　不同參數(shù)C下的局部動(dòng)量數(shù)值誤差

5　結(jié)論

本文首先基于多孔介質(zhì)理論，通過多孔介質(zhì)能量方程和本構(gòu)關(guān)系，首先推導(dǎo)出飽和多孔彈性桿局部熱平衡時(shí)的熱傳導(dǎo)方程;繼而在多辛算法理論的基礎(chǔ)上，采用廣義多辛算法研究了這一熱傳導(dǎo)問題，構(gòu)造出飽和多孔彈性桿局部熱平衡時(shí)熱傳導(dǎo)方程的廣義多辛形式和廣義多辛守恒律以及廣義多辛局部動(dòng)量守恒律誤差，并且得到其相應(yīng)的中點(diǎn)Box廣義多辛離散形式;最后，在已得的離散格式的基礎(chǔ)上，數(shù)值考察了飽和多孔彈性桿局部熱平衡時(shí)的熱傳導(dǎo)過程，將溫度場(chǎng)的數(shù)值結(jié)果與精確解進(jìn)行比較，并且數(shù)值考察了一段長時(shí)間內(nèi)的局部動(dòng)量誤差。得到如下結(jié)論:

1)本文構(gòu)造的廣義多辛格式能夠很好地模擬飽和多孔彈性桿的熱傳導(dǎo)過程，數(shù)值實(shí)現(xiàn)了熱傳導(dǎo)過程中溫度場(chǎng)隨時(shí)間的變化，并且數(shù)值解與精確解非常吻合，這充分顯示出本文采用的廣義多辛方法具有很好的有效性和精確性;

2)在局部動(dòng)量數(shù)值誤差的研究中，本文選取了不同的參數(shù)，均得到數(shù)量級(jí)在10－5以下的誤差結(jié)果，這表明本文構(gòu)造的廣義多辛格式能夠很好地保持系統(tǒng)的固有幾何性質(zhì)，并且具有良好的長時(shí)間數(shù)值穩(wěn)定性。

參考文獻(xiàn):

［1］徐曾和，徐小荷．飽和多孔地層中定量抽放的流固耦合問題［J］．巖土工程學(xué)報(bào)，1999，21: 737-741 Xu Zenghe，Xu Xiaohe．Fluid-Solid Coupling Problem in the Liquid Extraction at Fixed Flux［J］．Chinese Journal of Geotechnical Engineering，1999，21: 737-741 (in Chinese)

［2］De Boer R，Svanadze M．Fundamental Solution of the System of Equations of Steady Oscillations in the Theory of Fluid-Saturated Porous Media［J］．Transport in Porous Media，2004，56: 39-50

［3］秦冰，陳正漢，方振東，孫樹國，方祥位，王駒．基于混合物理論的非飽和土的熱-水-力耦合分析模型Ⅰ［J］．應(yīng)用數(shù)學(xué)與力學(xué)，2010，31: 1476-1488 Qin Bing，Chen Zhenghan，F(xiàn)ang Zhendong，Sun Shuguo，F(xiàn)ang Xiangwei，Wang Ju．Analysis of Coupled Thermo-Hydro-Mechanical Behavior of Unsaturated Soils Based on Theory of MixturesⅠ［J］．Applied Mathematics and Mechanics，2010，31: 1476-1488 (in Chinese)

［4］Feng K．On Difference Schemes and Symplectic Geometry［C］∥Proceeding of the 1984 Beijing Symposium on D D，1984: 42-58

［5］趙長印，廖新浩，劉林．辛算法在動(dòng)力天文中的應(yīng)用［J］．計(jì)算物理，1992，4: 812-812 Zhao Changyin，Liao Xinhao，Liu Ling．Application of Symplectic Algorithm in Dynamic Astronomy［J］．Chinese Journal of Computational Physics，1992，9(4) : 812-812 (in Chinese)

［6］Bridge T J．Multi-Symplectic Structures and Wave Propagation［J］．Mathematical Proceedings of the Cambridge Philosophical Society，1997，121: 147-190

［7］Hu W P，Han S M，Deng Z C，F(xiàn)an W．Analyzing Dynamic Response of Non-Homogeneous String Fixed at Both Ends［J］．International Journal of Non-Linear Mechanics，2012，47: 1111-1115

［8］Hu W P，Deng Z C，Han S M，Zhang W R．Generalized Multi-Sympletic Integrators for a Class of Hamiltonian Nonlinear Wave PDEs［J］．Journal of Computational Physics，2013，235: 394-406

［9］胡偉鵬，鄧子辰．大阻尼桿振動(dòng)的廣義多辛算法［J］．動(dòng)力學(xué)與控制學(xué)報(bào)，2013，11: 1-4 Hu Weipeng，Deng Zichen．Generalized Multi-Sympletic Method for Vibration of Big Damping Rod［J］．Journal of Dynamics and Control，2013，11: 1-4 (in Chinese)

［10］Yang X．Gurtin-Type Variational Principles for Dynamics of a Non-Local Thermal Equilibrium Saturated Porous Medium［J］．Acta Mechanica Solida Sinica，2005，18: 37-45

［11］楊驍，程昌鈞．流體飽和多孔介質(zhì)的動(dòng)力學(xué)Gurtin型變分原理和有限元模擬［J］．固體力學(xué)學(xué)報(bào)，2003，24: 267-276 Yang Xiao，Cheng Changjun．Gurtin Variational Principle and Finite Element Simulation for Dynamical Problems of Fluid-Saturated Porous Media［J］．Acta Mechanica Solida Sinica，2003，24: 267-276 (in Chinese)

Generalized Multi-Symplectic Method and Numerical Experiment for Thermal Conduction of Saturated Poroelastic Rod

Liu Xuemei1，2，Deng Zichen1，Hu Weipeng1

(1．Department of Engineering Mechanics，Northwestern Polytechnical University，Xi'an 710072，China 2．Department of Engineering Mechanics，Chang'an University，Xi'an 710064，China)

Abstract:Based on the theory of porous media，the thermal conduction equation of fluid saturated poroelastic rod is established by using the energy equation and constitutive relations of the two constitutes firstly in this paper．Then introducing orthogonal variables，we use the generalized multi-symplectic method to derive a first-order generalized multi-symplectic form for thermal conduction equation and several errors of conservation laws illustrating the local properties of the system．Thirdly，a midpoint box generalized multi-symplectic scheme is constructed; furthermore，discrete errors of generalized multi-symplectic conservation law and generalized local momentum conservation law are also obtained．Finally，the dissipation effect in thermal conduction process of saturated poroelastic rod and generalized local momentum conversation law are investigated numerically; moreover，the influence of parameter values for thermal conduction process is established later．From results of the numerical experiments，it can be preliminarily concluded that the generalized multi-symplectic scheme constructed in this paper has excellent accuracy，longtime numerical behavior and good conservation properties．

Key words:boundary conditions，constitutive equations，errors，matrix model，momentum，numerical methods，porosity，temperature distribution，thermal conductivity; dissipation，generalized multi-symplectic，saturated poroelastic rod，thermal conduction

作者簡(jiǎn)介:劉雪梅(1980—)，女，西北工業(yè)大學(xué)博士研究生，主要從事多孔介質(zhì)問題的保結(jié)構(gòu)算法研究。

收稿日期:2014-09-10基金項(xiàng)目:國家自然科學(xué)基金(11372252、11172239和11372253)與中央高?；?2014G1121096)資助

文章編號(hào):1000-2758(2015) 02-0265-06

文獻(xiàn)標(biāo)志碼:A

中圖分類號(hào):O241．82