饱和多孔弹性杆热传导的广义多辛方法及其数值实现



饱和多孔弹性杆热传导的广义多辛方法及其数值实现

刘雪梅1，2，邓子辰1，胡伟鹏1

(1．西北工业大学力学与土木建筑学院，陕西西安710072; 2．长安大学工程力学系，陕西西安710064)

摘要:首先根据多孔介质理论，利用饱和多孔介质的能量方程和本构关系，推导出饱和多孔弹性杆局部热平衡的热传导方程;继而引入正交变量，将热传导方程导入Hamilton系统，得到饱和多孔弹性杆热传导方程的广义多辛形式和多种局部守恒律形式;接着采用中点离散方法对热传导方程的广义多辛形式进行数值离散;最后利用计算机数值实现了饱和多孔弹性杆的热传导过程，并且讨论了参数取值的不同对热传导过程的影响，同时在数值模拟过程中记录了广义多辛格式的局部动量误差。研究结果表明，构造的广义多辛方法能够很好地模拟系统的热传导过程和耗散效应，同时也可长时间保持系统的固有几何性质。

关键词:多孔介质;广义多辛;耗散;热传导

饱和多孔介质通常是指孔隙中充满液体的多孔连通介质，其相关理论在地震工程学、土动力学、地球物理学、生物力学等领域有着广泛应用。目前，研究多孔介质宏观力学行为的主要理论有: Biot理论、多孔介质理论和混合物理论。在Biot理论基础上，徐曾和等研究了含单一圆井的饱和多孔地层中，以定流量方式开采时的流固耦合问题［1］; de Boer等基于多孔介质理论，研究了流固两相微观不可压饱和多孔介质稳态振动的基本解［2］;秦冰等以混合物理论为基础，将非饱和土视为由土骨架、液态水、水蒸气、干燥气体及溶解气体共5种组分构成的混合物，系统研究了非饱和土的热-水-力多场耦合问题［3］。然而，有关各类饱和多孔结构热传导过程的数值算法研究较少，尤其是热传导方程的求解过程中如何体现由于热量传递引起的耗散效应以及如何保持系统的局部能量守恒律和局部动量守恒律等局部几何性质更是数值算法的难点。

1984年是应用数学和计算力学研究中具有重要意义的一年，也是辛几何算法开创性发展的一年，我国数学家冯康于这一年首次系统地提出了哈密顿方程和哈密顿算法，并指出了从辛几何内部系统构成算法并研究其性质的途径，从而开创了辛几何算法的新领域［4］。自此以后，辛几何算法得到了长足的发展，学术界将纯理论的辛几何和现代科学工程计算有机地结合起来，使其广泛地应用于天体动力学和分子动力学等诸多领域［5］。在此基础上，Bridge教授针对无穷维Hamilton系统提出了多辛算法［6］，用以在数值求解过程中保持系统的局部几何性质。然而，实际力学系统中耗散效应的存在却成了多辛算法的“瓶颈”，近年来，广义多辛算法理论的提出［7-9］，有效地解决了这一"瓶颈"问题，将保结构算法理论体系拓宽至耗散系统。

多孔介质传热现象遍及于工农业生产的各个领域，例如:埋地电缆和直流接地极的热耗散;地源热泵和地冷空调中换热器埋管;高温原件的多孔冷却;生物多孔介质的传热问题等等。因此，多孔介质传热理论和应用研究是一个具有重要学术和应用价值的研究方向，而作为多孔介质一维传热问题中最基础的饱和多孔弹性杆的热传导就更加重要，它的研究将为多孔介质杆系结构以及植物根茎等多种多孔介质结构的传热传质问题提供有力的依据，尤其是其数值求解的应用将更加广泛。正是基于这一点，本文试图通过广义多辛保结构算法来寻求新的研究饱和多孔弹性杆热传导过程的数值方法，同时也为了揭示其热传导过程中的耗散效应以及系统的广义多辛守恒律等多种局部几何性质，从而为多孔介质传热问题的数值求解提供新的途径。

本文基于饱和多孔介质理论，首先建立饱和多孔弹性杆热传导的一维数学模型。其次，通过正则变换，构造饱和多孔弹性杆热传导方程的广义多辛形式，并研究了热量传递的耗散效应在广义多辛形式中的数学表述。随后，采用中点离散方法离散广义多辛形式，得到其中点Box广义多辛离散格式。最后，利用中点Box广义多辛离散格式数值模拟饱和多孔弹性杆的热传导过程，并将数值解与精确解进行比对，同时研究了热传导过程中的耗散效应。本文的广义多辛分析将完善保结构算法理论，也为多孔介质耗散系统的数值求解提供新的途径。

1　饱和多孔弹性杆的热传导方程

为了研究饱和多孔弹性杆的传热问题，本文建立如图1的模型:设长为L、横截面高度为H的流体饱和多孔弹性杆由不相溶的微观不可压流体和微观不可压弹性多孔固相骨架组成，其侧表面不透水。根据多孔介质理论，每一相物质被理想化地分布在整个区域中，并拥有各自独立的运动，记物质粒子的位移为:

图1　LH的饱和多孔弹性杆

式中，i =S、F分别表示固相和流相，V为多孔介质初始时的空间区域。忽略两相间的质量交换，则流固两相的能量方程可表示为［10］:

式中，ρi为两相的宏观质量密度;εi、σi、qi和ri分别为两相的内能、宏观应力张量、热流矢量和内部热源;i和i分别为两相间的相互作用力和能量交换量，满足+= 0+= 0; ()i=+ grad(…)·i表示随物相的物质时间导数。

对于各向同性线弹性多孔固相骨架和无粘性流体，在两相不可压和小变形的假定下，并忽略由于流固两相的相对运动引起的附加热对流和附加热传导，可得流固两相的宏观应变张量和宏观应变率张量分别为:

本构方程可表示为［10］:

式中，上标T表示转置，p和Sν分别表示有效孔隙水压力和与流体渗流系数有关的两相相互作用耦合系数，eΘ为非局部热平衡引起的两相间能量交换系数，μS和λS为固相宏观拉梅常数，kii和θi(i = S，F)分别为热传导系数和流固两相的温度，wF=F－S为孔隙流体相对固相骨架的速度。

以下考虑局部热平衡情况，即在每一个空间点上，固相骨架和孔隙流体具有相同的温度，其温度可表示为θS=θF=θ，将(3)式和(4)式代入能量方程(2)中，并采用如图1所示的一维空间坐标轴，同时忽略非线性项对热传导的影响［10］，可得饱和多孔弹性杆局部热平衡时的热传导方程:

式中，θ和ci(i = S，F)分别表示局部热平衡温度和流固两相的比热容，再引入如下无量纲量和常量:

2　饱和多孔弹性杆热传导方程的广义多辛形式

对于饱和多孔弹性杆热传导方程(6)，引入正交变量φ=əxθ，上述二阶偏微分方程(6)可转化为以下一阶耦合Hamilton方程组形式:

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

从(8)式容易看到，系数矩阵M并不是反对称的，因此，按照广义多辛算法的概念［8］，(7)式为广义多辛形式。这是由于本文涉及的饱和多孔弹性杆的热传导问题中存在着热量传递，从而使系数矩阵M并不是反对称的，也正是由于非完全反对称矩阵的存在，上述广义多辛形式(7)式就不能满足多辛算法的多辛守恒律和局部动量守恒律。

基于多辛理论，由多辛守恒律概念可给出广义多辛形式(7)式的广义多辛守恒律为［8］:

同时，由局部动量守恒律误差的概念，亦可给出广义多辛形式(7)式的局部动量守恒律误差为:

这里需要指出，由于系数矩阵K是严格反对称的，广义多辛形式(7)式严格满足局部能量守恒律。

3　饱和多孔弹性杆热传导方程的广义多辛离散

本文选用中点差分离散方法分别在空间方向和时间方向对广义多辛形式(7)式进行差分离散，并分别选取时间步长Δt和空间步长Δx对计算区域均匀划分，用zji表示状态变量z在(xi，tj)网格点上的近似值，则可得时间和空间2个方向上的中点Box离散格式:

式中

将系数矩阵(8)和状态变量z =［θ，φ］T以及Hamilton函数代入上述离散格式(11)，便可得到广义多辛偏微分方程(7)的等价形式:

对于上述已得到的广义多辛守恒律(9)式，亦可通过有限差分运算得到其离散格式。这里采用向前差分运算，得到广义多辛局部动量误差的离散形式［7］:

为了得到局部动量误差随时间变化的平面图，本文取广义多辛局部动量误差在第j时间步的绝对值［7］:

4　数值实验

这一节中，本文将利用以下饱和多孔弹性悬臂杆的数值算例来探讨和验证前述已构造出的广义多辛离散格式的精确性、有效性和数值稳定性等性能。

在广义多辛离散格式(12)中，考虑如下边界条件和初始条件:

选取计算空间步长Δx = 0. 01，时间步长Δt =0. 5，λ2= 0. 000 1，分别取参数C = 0. 043和C = 0. 02，在t∈［0，200］内，得到局部热平衡时饱和多孔弹性杆不同时刻的温度随空间坐标的变化图(见图2和图3)。为了验证广义多辛离散格式(12)的精确性，采用分离变量法，可得无量纲偏微分方程(6)的定解问题(15)～(16)的精确解:

式中

因此，亦可得局部热平衡时饱和多孔弹性杆不同时刻温度场的精确解随空间坐标的变化图(图2和图3)。

图2　C=0. 043时不同时刻温度的数值解(点)与精确解(曲线)的比较

图3　C=0. 02时不同时刻温度的数值解(点)与精确解(曲线)的比较

图2反映了参数C = 0. 043时饱和多孔弹性杆局部热平衡的热传导过程，这一过程中，由于杆x = 1端与外界绝热，初始温度从杆x = 1端向杆x = 0端逐渐降低，因此随着时间的推移，热量将从杆x = 1端向杆x = 0端广泛传递，从初始时刻到t = 200的过程中，杆x = 1端的温度从0. 5降低至0. 16;热传导过程中，由于杆内部无热源，杆x = 0端温度始终为零，杆将由x = 0端持续向外界传热，杆内各点的温度均逐渐下降。与图2类似，图3反映了参数C = 0. 02时饱和多孔弹性杆局部热平衡的热传导过程，这一过程中明显可以看到，相同时刻杆内各点的温度均比C = 0. 043时低，t = 200时，杆内各点温度均已接近零;由此可以说明，参数C的取值越小，热传导过程将越快;这种现象是因为参数C取值的减小，相当于饱和多孔弹性杆热传导方程的导热系数增大，则必然导致热传导过程加快。此外，由图2和图3均可看到，热传导过程中温度场分布的数值解与精确解非常吻合。以上这些结果说明本文采用的广义多辛方法具有很好的有效性和精确性。

为了进一步验证广义多辛格式(12)的局部保结构性，本文将无量纲模拟时间延长至500，利用(14)式，将第j时间步中误差绝对值最大的空间点的数值误差作为该时间步的局部动量数值误差，并取参数C足够小(其中C分别取0. 01和0. 043) ;继而在计算机模拟过程中记录每一时间步的局部动量误差，从而得到广义多辛局部动量数值误差随时间的变化曲线图(见图4)。

从局部动量数值误差的结果可以看出:整个模拟过程中，局部动量误差均在10－5数量级以下，这说明本文构造的广义多辛算法能够很好地模拟系统的热量传递耗散效应，并且保持了系统的固有几何性质。从图4可以看到，参数C取0. 043时，无量纲时间从t = 450之后，局部动量误差为零;然而参数C取0. 01时，无量纲时间从t = 100之后，局部动量误差就已经为零。这就是说，参数C的取值越小，局部动量误差将越早等于零，这也就表示，参数C的取值越小，越能体现本文构造的广义多辛形式的局部保结构性。

图4　不同参数C下的局部动量数值误差

5　结论

本文首先基于多孔介质理论，通过多孔介质能量方程和本构关系，首先推导出饱和多孔弹性杆局部热平衡时的热传导方程;继而在多辛算法理论的基础上，采用广义多辛算法研究了这一热传导问题，构造出饱和多孔弹性杆局部热平衡时热传导方程的广义多辛形式和广义多辛守恒律以及广义多辛局部动量守恒律误差，并且得到其相应的中点Box广义多辛离散形式;最后，在已得的离散格式的基础上，数值考察了饱和多孔弹性杆局部热平衡时的热传导过程，将温度场的数值结果与精确解进行比较，并且数值考察了一段长时间内的局部动量误差。得到如下结论:

1)本文构造的广义多辛格式能够很好地模拟饱和多孔弹性杆的热传导过程，数值实现了热传导过程中温度场随时间的变化，并且数值解与精确解非常吻合，这充分显示出本文采用的广义多辛方法具有很好的有效性和精确性;

2)在局部动量数值误差的研究中，本文选取了不同的参数，均得到数量级在10－5以下的误差结果，这表明本文构造的广义多辛格式能够很好地保持系统的固有几何性质，并且具有良好的长时间数值稳定性。

参考文献:

［1］徐曾和，徐小荷．饱和多孔地层中定量抽放的流固耦合问题［J］．岩土工程学报，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］．应用数学与力学，2010，31: 1476-1488 Qin Bing，Chen Zhenghan，Fang Zhendong，Sun Shuguo，Fang 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］赵长印，廖新浩，刘林．辛算法在动力天文中的应用［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，Fan 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］胡伟鹏，邓子辰．大阻尼杆振动的广义多辛算法［J］．动力学与控制学报，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］杨骁，程昌钧．流体饱和多孔介质的动力学Gurtin型变分原理和有限元模拟［J］．固体力学学报，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

作者简介:刘雪梅(1980—)，女，西北工业大学博士研究生，主要从事多孔介质问题的保结构算法研究。

收稿日期:2014-09-10基金项目:国家自然科学基金(11372252、11172239和11372253)与中央高校基金(2014G1121096)资助

文章编号:1000-2758(2015) 02-0265-06

文献标志码:A

中图分类号:O241．82