多孔介质是存在于自然界的重要介质之一, 在各类工程实际中均有广泛的应用。在土木工程领域, 大量多孔介质结构以承担横向荷载为主, 其动力响应以横向弯曲变形为主, 例如混凝土梁、木质结构梁等。近些年来, 学术界已对多孔介质梁的各类动力响应问题做了大量的研究^[1-4]。Li等^[1]通过变分原理, 给出流体饱和多孔弹性梁轴向扩散的混合有限元格式, 并数值分析其位移和孔隙压力随时间的响应。Wang等^[3]给出了饱和多孔弹性梁纯弯曲变形的三维解析解。周凤玺等^[4]基于不可压多孔介质理论和弹性地基模型理论, 利用Fourier级数展开法研究了饱和多孔弹性梁的自由振动特性。

然而, 对于短粗梁, 由于剪切变形效应显著, 应当采用考虑剪切效应的梁模型。Timoshenko梁模型理论是被广泛应用于各种工程结构的分析和研究中的梁模型^[5-9]。吴峰等^[7]研究了Timoshenko梁的波传播问题, 分析了Timoshenko梁的能带结构和波散射特性。Kiani等^[8]基于Timoshenko梁理论和Biot固结模型, 研究了横向移动荷载作用下的饱和多孔弹性梁结构的横向振动, 结果揭示了剪切变形效应对某些多孔介质梁结构的重要性。Chouiyakh等^[9]研究了移动质量荷载作用下Timoshenko梁的振动和多裂隙检测问题。

对于多孔介质梁结构的动力响应问题, 目前的研究方法多数为积分变换、有限元法等, 在解决实际问题的过程中, 也会出现诸如求解过程复杂繁长、数值结果的稳定性较差等情况, 特别是已有算法对多孔介质流固耦合系统的固有几何性质鲜有涉及。基于这一点, 本文采用其他数值方法研究饱和多孔弹性短粗梁的动力问题。

对于保守哈密顿动力学系统, 辛算法和多辛算法因其良好的精确性、有效性和保结构性, 受到众多研究者的广泛关注^[10-14]。高强等^[10]采用保辛数值积分方法, 研究了拉压刚度不同的桁架非线性动力问题。Mcdonald等^[11]采用多辛离散格式研究了半线性波动方程的行波解。彭海军等^[12]采用标准辛格式解决了基于线性时变系数的滚动时域控制问题。胡伟鹏等^[13-14]采用多辛算法, 分析和模拟了Ⅱ类超导体混合态的电磁特性和拟Degasperis-Procesi方程的非线性特性。对于大量存在耗散效应的无限维非保守哈密顿动力学系统的求解问题, 近年来, 广义多辛算法应运而生, 其相关研究成果已经被广泛报道^[15-18]。胡伟鹏^[16]基于多辛积分理论, 将多辛数学方法推广到可应用于无限维非保守动力学系统的广义多辛数学方法。张宇等^[18]对弦的有阻尼受迫振动过程作了广义多辛数值算法研究。

对于饱和多孔介质短粗梁结构, 由于梁的可渗透性, 随着时间的推移, 固相骨架将承担所有外力, 故固相骨架可采用与单相弹性梁相同的梁模型^{[2, 8, 19]}。因此, 本文将Timoshenko梁模型应用于饱和多孔弹性短粗梁, 并采用广义多辛算法研究梁的流固耦合动力响应特性。本文构造了梁动力响应控制方程的广义多辛离散格式和各种局部守恒律离散格式, 数值模拟了两端可渗透饱和多孔弹性Timoshenko悬臂梁在自由端受集中荷载作用的动力响应过程, 并考察了固相挠度和固相有效应力以及孔隙流体压力等效力偶的响应特性, 同时数值分析了系统的广义多辛守恒律误差、局部动量误差和局部能量误差。

1 不可压饱和多孔弹性Timoshenko梁的横向动力响应模型

考虑承受横向分布载荷q(x, t)的饱和均匀多孔线弹性Timoshenko梁, 如图 1所示, 梁长为L，横截面积为A，横截面绕oy轴的惯性矩为I。梁由互不相溶的微观不可压流体和不可压横观各向同性固相骨架组成, 其侧表面不透水, 并忽略流相和固相的体积力、两相间的质量交换及能量交换。

对于自然界和实际工程中的一些多孔介质梁结构, 例如植物根茎、动物骨骼和传热管道等, 由于梁骨架结构的特殊性, 使得孔隙流体以轴向运动为主, 可忽略其横向的运动, 为此本文首先假定孔隙流体仅沿梁的轴向运动。参阅文献[19], 假定固相小变形、梁横截面中性轴无纵向位移, 忽略孔隙流体的转动惯性效应, 利用不可压流体饱和多孔介质的混合物动量方程、孔隙流体动量方程及体积分数方程, 并根据Timoshenko梁模型假设^[8]

(1)

得如下饱和多孔弹性Timoshenko梁横向动力响应的控制方程组：

(2)

(3)

(4)

其中固相挠度w^s(x, t)、梁横截面转角ϕ(x, t)及孔隙流体压力等效力偶M_p(x, t)为基本未知量, 分别记为^[19]

(5)

以上各式中, E^s, G^s和K^s分别为固相骨架的杨氏模量、剪切模量和弹性模量, k为Timoshenko梁剪切修正因子, n^f为流相的体积分数, ρ^s和ρ^f分别为固相和流相的宏观质量密度, σ^se和ε^s分别为固相弹性骨架的有效应力和应变, p为孔隙流体压力, 流固两相耦合作用系数S_v与Darcy渗透系数K′的关系为S_v=(n^f)²γ^fr/K′, 其中γ^fr为流体的真实比重。

2 饱和多孔弹性Timoshenko梁横向动力响应的广义多辛算法

为了计算方便, 引入无量纲量和变量

(6)

式中，s和K₁分别表示梁的长细比和剪切效应, 可得无量纲形式的饱和多孔弹性Timoshenko梁横向动力响应的控制方程组：

(7)

(8)

(9)

为了书写方便, 上述方程及后文中仍用量纲形式的符号表示无量纲量。

引入中间变量∂_xϕ=ξ, ∂_x=ζ, ∂_xw=χ, ∂_tϕ=α, ∂_tw=β, 并定义z=[ϕ, M_p, w, ξ, ζ, χ, α, β] ∈R⁸作为状态变量, 可将控制方程组(7)~(9)式写成一阶Hamilton偏微分方程组

(10)

式中，Hamilton函数

系数矩阵分别为

(11)

(12)

由于(11)式和(12)式存在对称矩阵, 则一阶偏微分方程组(10)式不是严格多辛形式, 且不满足多辛局部动量守恒律和局部能量守恒律。根据广义多辛理论^[15-16], 可对这一存在耗散效应的动力响应系统做广义多辛近似。(10)式的广义多辛条件为^[16]

(13)

式中, Δ和o(Δt², Δt, Δx², Δx)分别为动力响应系统的广义多辛守恒律误差和差分截断误差。一阶Hamilton偏微分方程组(10)式的广义多辛局部动量误差和局部能量误差分别为^[16]

(14)

(15)

为了考察广义多辛算法的数值特性, 采用Preissmann Box离散格式^[18]对耗散动力学系统(10)式进行如下数值离散

(16)

采用有限差分近似, 分别给出第j时间层广义多辛局部动量数值误差和局部能量数值误差^[15]

(17)

(18)

另外, 广义多辛守恒律误差的离散格式也可通过差分运算给出^[16]。

3 饱和多孔弹性Timoshenko悬臂梁的广义多辛动力响应特性

学术界已对多孔介质简支梁模型做了一定的研究^{[6, 8, 19]}, 因此本文进一步针对工程中存在的悬臂梁模型, 采用广义多辛算法研究饱和多孔弹性Timoshenko悬臂梁的动力响应特性。设梁长为L、高为H、横截面积为A, 两端可自由渗透(P-P), 自由端受集中荷载Q(t)的作用, 且初始未变形, 如图 2所示。

由于仅在自由端受集中荷载, 则横向振动控制方程组(7)~(9)式中, q(x, t)=0, 且该悬臂梁的边界条件和初始条件分别为^{[2, 8]}

(19)

(20)

引入无量纲量, 并仍用Q表示无量纲量Q, 则有如下无量纲边界条件和初始条件

(21)

(22)

分别选取Δt=0.1和Δx=0.01为时间步长和空间步长, 取η=0.9及γ=3。不失一般性, 取Q=10^-3, 采用离散格式(16)式模拟多孔弹性Timoshenko悬臂梁的动力响应过程。

通过图 3发现, 随着梁长细比增大, 梁各横截面固相骨架的横向衰减振动周期和振动幅值均增大, 且达到各自稳定值的时间加长。其稳定值为相应的经典单相弹性Timoshenko悬臂梁的静挠度值, 这是因为梁两端可渗透, 则随着时间的推移, 多孔弹性Timoshenko梁的固相骨架将承担所有外荷载, 从而其固相挠度将最终达到单相弹性梁的静响应值。

从图 3也可得到, 随着梁长细比增大, 梁各横截面固相骨架挠度的稳定值越接近于相应经典单相弹性Euler-Bernoulli梁的静挠度值(图中虚线所示), 这说明, 梁长细比越大, 造成梁的剪切变形效应越小, 则固相挠度的稳定值越接近于Euler-Bernoulli梁的静挠度值, 静挠度解析解如下所示

(23)

图 4反映出随着流固两相耦合作用系数S_v增大, 梁固相各截面的衰减振动过程缩短, 这是由于该系数具有阻尼效应, 则随该系数值增大, 系统能量衰减加快, 从而引起梁各截面固相骨架的衰减振动过程缩短。但两相耦合作用系数并不影响和改变梁的挠度稳定值。

通过图 5发现，x=0.2, 0.5, 0.8各横截面的z=0.5位置, 固相有效正应力均呈振荡减小的过程, 且随流固两相耦合作用系数S_v的增大, 有效正应力的振荡减小过程均缩短。最终, 应力均达到各自相应的稳定值, 其稳定值就是经典单相Timoshenko悬臂梁的弯曲应力值(虚线所示), 且有效正应力值均大于相应单相Euler-Bernoulli梁的稳态正应力值。应力振荡过程中, 耦合作用系数增大, 即流相的Darcy渗透系数减小, 固相应力响应的幅值减小。

通过图 6可以得出:随着梁长细比增大, 梁中点横截面孔隙流体压力等效力偶的衰减振荡过程越长, 振荡周期越大; 反之, 随着两相耦合作用系数增大, 孔隙流体压力等效力偶的衰减振荡过程明显缩短。另外, 由于梁固相骨架最终将承担所有外载荷, 则梁中点截面最终将达到完全渗透, 孔隙流体压力等效力偶为零。同理可以推知, 梁各横截面最终均将达到完全渗透, 孔隙流体压力等效力偶都将振荡减小至零。

将本文与文献[6]对比, 发现由于流固两相的耦合作用, 饱和多孔弹性Timoshenko悬臂梁的固相骨架均作衰减振动, 且耦合作用效应越强, 衰减振动过程越短, 本文结果与文献[6]的研究结果一致。对比本文数值结果和文献[8]中受移动荷载作用的饱和多孔弹性Timoshenko简支梁的有限元数值结果, 发现梁长细比越大, 多孔弹性Timoshenko梁的固相挠度、孔隙流体等效力偶随时间的响应越接近于Euler-Bernoulli梁的响应过程, 两者结论一致。这些对比结果进一步说明本文所得数值解的有效性和采用广义多辛数值算法的合理性。

为了考察本文研究的耗散动力学系统的局部保结构性能, 在表征耗散效应的无量纲参数满足广义多辛条件(13)式的前提下, 在无量纲时间段t∈[0, 100]内, 分别数值考察了饱和多孔弹性Timoshenko悬臂梁动力响应过程的广义多辛守恒律误差、广义多辛局部动量误差和广义多辛局部能量误差, 如图 7至9所示。3种数值误差的数量级均在10^-8或10^-8以下, 且随着造成出现对称矩阵的流固两相耦合作用系数S_v和参数1/2K₁的减小(1/2K₁的减小对应梁长细比s的增大), 广义多辛守恒律误差及各种局部守恒律误差的数量级也减小, 这说明非保守系统的耗散效应越小, 系统越接近于保守系统。另外, 数值误差在考察时间段内均稳定分布。

4 结论

本文通过引入中间变量和辛降阶, 将饱和多孔弹性Timoshenko梁动力响应控制方程转化为一阶偏微分方程组, 并构造广义多辛离散格式, 从而简化了计算。研究了两端可渗透和自由端受集中荷载的多孔弹性Timoshenko悬臂梁的流固耦合响应特性, 结果发现:

1) 梁各截面固相骨架均作衰减振动, 达到各自稳态值所需的时间随梁长细比的增大而加长, 随流固耦合作用系数的增大而缩短。固相挠度稳态值不受耦合系数的影响, 且梁长细比越大, 该稳态值越小, 越接近于相应经典单相弹性Euler-Bernoulli梁的静挠度值, 这说明梁剪切变形效应越小, 梁动力响应过程越接近多孔弹性Euler-Bernoulli梁的响应过程。

2) 梁各截面的固相有效应力均呈振荡减小状态, 随耦合作用系数增大, 固相应力幅值减小, 振荡过程缩短, 并最终达到相应的单相Timoshenko悬臂梁的应力稳定值, 且固相有效正应力值大于单相Euler-Bernoulli梁的稳态正应力值。

3) 孔隙流体压力等效力偶振荡衰减, 衰减过程随梁长细比的减小而缩短, 随流固耦合作用系数的增大而缩短, 并最终等于零。

4) 长时间内数值考察广义多辛离散格式的局部保结构性, 发现数值误差的数量级均非常小, 且表征耗散效应的参数越小, 数值误差的数量级越小。说明广义多辛离散格式具有长时间的数值稳定性, 且非保守系统的耗散效应越小, 系统越接近保守系统。