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

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

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

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

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

图 1 L>>H的饱和多孔弹性杆

图选项

点击浏览公式

点击浏览公式

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

点击浏览公式

本构方程可表示为^[10]:

点击浏览公式

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

点击浏览公式

点击浏览公式

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

点击浏览公式

点击浏览公式

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

基于多辛理论,由多辛守恒律概念可给出广义多辛形式(7)式的广义多辛守恒律为^[8]:

点击浏览公式

点击浏览公式

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

点击浏览公式

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

点击浏览公式

点击浏览公式

点击浏览公式

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

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

点击浏览公式

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

点击浏览公式

点击浏览公式

图 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)。

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

图选项

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

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

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

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