微分求积法在结构动力分析中的应用1

任争争 梅雨辰 李鸿晶

（南京工业大学，土木工程学院，南京 211816）

引言

地震等灾害环境的作用使结构产生复杂的动态响应，可能导致结构失效、破坏甚至倒塌，从而形成灾害、造成损失。对新建结构进行抗灾设计及对已建结构进行抗灾加固是防御和减轻工程灾害及其损失的有效途径，这就要求认识结构在这些灾害作用下的性态和响应行为，而结构动力分析则是实现这一目标的基本手段。

结构动力分析的本质是实现对动力荷载激励下的结构运动微分方程（组）的求解，属于常微分方程的初值问题。目前用于分析结构动力的方法大致可以分为2类：一类是变换方法，如振型叠加法利用振型的正交性和完备性将结构动态响应向各阶振型分解，再通过叠加各阶振型响应以获得结构动态响应的结果，这类方法采用叠加原理，一般只用于线弹性结构动力分析，且要求结构具有经典阻尼特性；另一类为直接方法，即直接对结构运动微分方程进行求解而不必引入任何假定，这类方法既可用于线弹性结构，也可用于非线性结构的分析，以Newmark-β法（Newmark，1959）等逐步积分法为代表，一般采用数值求解手段。

由于现代结构不断向大型化、复杂化发展，加之结构精细化模型的采用，导致结构动力分析的计算需求呈爆发式增长，因而需要寻求高效率的分析方法。微分求积法（DifferentialQuadrature Method，DQM）是由Bellman等（1971，1972）发展起来的1种求解微分方程的数值方法，可通过较小的计算工作量获得较高的计算精度。DQM要求求解域比较规则，在工程分析中一般用于时间无关问题的求解。如在结构力学领域，多用来求解静力问题和固有振动问题等（Bert等，1988，1993，1996，1997；Jang等，1989；Kang等，1995，1996；Liew等，1996a，1996b；Sherbourne等，1991；Striz等，1988；Wang等，1993，1994；Zeng等，2001）。这类问题都属于边值问题，即在DQM中计算的是解函数对空间坐标的导数。而结构动力分析属于初值问题，使用DQM对其求解的研究工作相对较少，Fung（2001a，2001b）、Liu等（2008）、李鸿晶等（2011a，2011b）和廖旭等（2013）开展过相关的研究工作。本文在此基础上发展1种基于DQM的结构动力分析的高精度方法。不同于静力边值问题，实际结构的动力响应问题有其特殊性，许多情况下时间跨度很长，像边值问题一样一次性对所有时间区域进行离散求解，将出现病态问题而产生错误。本文借鉴单元法的思想，以期提高结构动力分析的计算效率。

1 微分求积法基本原理

微分求积法是1种用于求解微分方程的数值方法。它的实质是将函数在某一离散节点处的各阶导数值，近似表示成计算域内所有节点处离散函数值的线性加权和，从而将复杂的微分方程化为关于离散点的线性方程（组）。由于本文讨论的是结构动力反应的计算，求解的运动微分方程仅是关于时间t的常微分方程，因此仅介绍一维区域内的微分求积原理。

设函数f(x)为在区间［a，b］上k阶连续可微，将区间［a，b］划分为m段，共（m＋1）个互不相同的节点，分别记为x0，x1，……，xm-1，xm，其中x0＝a，xm＝b。

根据计算数学的函数逼近理论，函数f（x）可做如下逼近：

其中，qj(x)为函数空间中各线性无关的基函数。

对式（1）求k阶导数，然后将所有节点代入，得到：

式（3）即为一维区域微分求积的基本公式。

常用的函数空间是（m＋1）维多项式空间，选择该空间中的所有基函数都能得到相同的权系数。最常见的多项式基函数是幂指数插值函数和拉格朗日插值函数2种，分别如式（4）、（5）所示：

选用式（4）的幂指数函数计算权系数，得到权系数的隐式表达式，需要求解范德蒙矩阵的逆矩阵，但当m较大时，不但计算量大，而且矩阵将容易出现病态。为了解决这些问题，目前大多采用式（5）给出的拉格朗日插值基函数，因为拉格朗日插值的各项系数为各节点的函数值，形式与式（1）完全相同，直接求k阶导数，便可得到相应的k阶权系数的显示表达式，计算效率大幅度提高。

实际计算权系数时往往只需计算1阶导数的权系数，其它各阶导数的权系数可由1阶权系数以矩阵的形式方便地表示：

其中，A表示1阶权系数矩阵，A（k）表示k阶权系数矩阵。1阶权系数可通过选定拉格朗日基函数，直接得到如下的表达式：

应用微分求积法时还需确定采样网格节点的位置，大致分为均匀网格点和非均匀网格点2大类。虽然均匀网格点的精度总体上没有非均匀点高，但是其在处理离散荷载，如地震荷载或风荷载时，可直接使用原始采样点作为节点，不需要对荷载进行额外的插值，计算效率比不均匀网格点更高。

2 结构动力反应微分求积分析方法

用上述微分求积法求解结构的动力反应，并以线弹性单自由度体系为例进行讨论。虽然实际结构大都为多自由度体系，且为非线性体系，但其动力反应都可以通过线性迭代和振型分解转化为线弹性单自由度体系。因而，讨论线弹性单自由度体系的微分求积法更具普遍意义，且简单直观、易于理解。

线弹性单自由度体系的动力反应运动方程为：

其中，ω、ξ分别表示体系的自振频率和阻尼比；u表示体系的位移；分别表示体系的速度和加速度；p（t）为结构所受到的随时间变化的荷载。

在实际工程中，动荷载随类型的不同，作用时间差异很大，短则瞬间（如冲击荷载），长则几分钟甚至几十分钟（如风荷载）。对于一些长时间作用的荷载，在整个荷载作用时域内实施微分求积法，要保证计算结果的精确性，需要成千上万的时间节点，这将导致求解时系数矩阵的阶数过于庞大，引起矩阵的条件数大，病态效果严重，无法得到可靠的结果。为了解决此问题，可以借鉴有限单元法的思想，将整个荷载持时划分为多个等间距的时间单元，称为1个时步，在每个时步内，使用微分求积法求解。

考虑动力荷载持时内长度Δt的时步［tj，tk］内的反应，tj、tk一般与荷载p（t）的采样时刻点重合。在时步内定义一局部坐标tt∈ [0，Δt]，坐标起点为tj，方向同时间坐标。为了方便处理，正则化局部坐标，作则定义域被正则化为［0，1］，时步内关于τ的运动方程可写为：

将时步离散为m段，记节点为τ0，τ1，……，τm，将分别简记为则在各离散节点处的方程为：

对这些节点使用微分求积法：

其中，aij是微分求积的权系数，将式（11）、（12）写成矩阵的形式：

已知时步的初始位移和速度分别为u0和0v，且将式（11）、（12）改写为：

将式（15）、（16）代入式（10），并将已知量与未知量分离在等式的两侧，整理后得到如下方程：

式（17）为各时步内微分求积的基本方程，求解该方程可得到时步内各点的位移反应，再运用微分求积原理，可进一步求得各节点速度反应：

将本时步末的位移um和速度作为下一时步的初始位移和速度，从初始0时刻开始逐时步进行求解，可得到荷载作用的所有时刻的位移和速度反应。除了位移和速度外，工程上关心的加速度可通过运动方程（8）变形得到：

3 算例

通过具体的算例验证上述微分求积法求解结构动力反应计算结果的精确性和可靠性。选择3种不同自振周期的单自由度体系，其频率范围大致覆盖低频、中频和高频，阻尼比都取为工程中常见的0.05，在上面施加不同频率的正弦荷载，体系的基本参数如表1所示，荷载信息如表2所示。

表1 体系的基本特性Table 1 Basic characteristics of systems

表2 简谐荷载的信息Table 2 Information of simple harmonic load

用微分求积法求解3种体系在以上简谐荷载下的动力反应。为了方便计算，时步内采用均匀网格离散方案，时步的长度Δt取为简谐荷载的周期。由于1个周期的长度是简谐荷载的最小重复单元，这样取值不仅能方便地分析时步分段数m对数值稳定性和精度的影响，更能清楚地发现荷载周期与步长Δt取值的规律。

首先，计算该时步长度下采用不同m所得到的体系的位移和速度反应，然后采用体系在简谐荷载下的解析解作为精确解进行校核。表3—5列出了分段数m为20内的偶数时，简谐荷载激励下3种体系动力反应DQM解与精确解的平均相对误差。

表3 荷载周期1s的平均相对误差Table 3 Average relative error with load period of 1s

表4 荷载周期0.2s的平均相对误差Table 4 Average relative error with load period of 0.2s

表5 荷载周期为0.1s的平均相对误差Table 5 Average relative error with load period of 0.1 seconds

从表3—5可以看出，当时步分段数m很小时，DQM计算结果严重偏离精确解，但随着m的增大，除少部分点外，误差大致呈迅速减小的趋势，个别数据（如m＝14时）反应误差巨大是由于算法在该时步分段下，时步长度取为荷载周期时出现了数值不稳定现象，初始误差随着逐时步推进不断放大，最后完全湮灭真实的结果。当m增大到10时，3种体系在不同周期的简谐荷载下的位移和速度反应都减小到了5%以内，对于实际工程已足够精确。为了更形象地描述这种精确程度，图1—3给出了0—20s内荷载周期为1s下，m＝10时DQM的位移计算结果与精确解，为了更清晰地进行对比，对每张图进行了局部的放大。

图1 简谐荷载激励下体系1的位移反应Fig.1 Displacement response of system 1 under simple harmonic load

图2 简谐荷载激励下体系2的位移反应Fig.2 Displacement response of system 2 under simple harmonic load

图3 简谐荷载激励下体系3的位移反应Fig.3 Displacement response of system 3 under simple harmonic load

从图中可以看出，DQM的计算结果与精确解几乎完全重合，充分说明了分段数m＝10时，DQM计算动力反应足够精确，且数值稳定性能得到很好的保证。以地震荷载为例，中低频地震荷载的大部分周期一般在0.2s以上，偏于保守取为0.2s，分10段后节点间距为0.02s。而目前一般的地震地面加速度采样周期仅为0.005s或0.1s，节点间距取采样周期的几倍以上，都能得到精确的结果，计算精度较高。

对比表3、4、5可以发现，当m由10继续向上增大时（排除m＝14时因失稳误差剧增的情况），误差基本上呈继续减小的趋势，很多情况下误差甚至降至1%或1‰以内。但这对工程实际已无太大的意义，反而会随着m的增大，不断地增大计算量。因此，对于结构在简谐荷载作用下的反应，采用均匀节点方案，即在1个荷载周期内取m＝10，时步内插入9个内节点，既能达到土木工程领域内精度的要求，又能最大限度地减少计算工作量，是相对较优的时步节点数量。而且，以往DQM求解工程结构静力问题的大量经验表明，m取8—10可得到满意的计算结果（王鑫伟，1995），这与本文得出的m＝10的取值较优也相符，这虽然是DQM解决静力问题的经验，但在1个周期内的动力问题与静力问题存在不少的联系，从而也可间接说明本文将m＝10作为较优参数的合理性。

由表3—5还可以发现，对于简谐荷载激励下的反应，当m＝10时，取时步长度Δt与荷载周期相等时，计算误差已在工程可接受的范围内。因此，计算体系在简谐荷载激励下的反应，在选定m＝10的前提下，将时步长度和简谐荷载的周期保持一致即可。

虽然简谐荷载是1种理想的荷载，但是对计算实际荷载激励下的反应具有重要的意义，这是因为工程中实际的动荷载都是一定持时的暂态荷载，都可以通过周期延拓表示成一系列简谐荷载的叠加。实际计算时，可以首先采用谐波分析，检测出其包含的不同周期的简谐波的幅值，得到该动力荷载的频谱；然后以最大幅值所对应的周期，即卓越周期为基准，确定等效周期；最后，取时步的长度为荷载的等效周期，将时步等分成10段进行计算，即可得到较精确的计算结果。

4 结论

本文将微分求积法引入结构动力反应的分析与计算，并通过数值算例得到如下的结论：

（1）采用微分求积法求解结构在动力荷载激励下的反应合理可行。在较大的节点距离下依然能够得到精确的结果，计算精度高，且具有普适性，对不同自振周期的结构、不同频率的动力荷载都适用。

（2）用DQM进行动力分析时，能一次性求得多个时刻的反应。相比传统的每次只能求1个时刻的逐步积分法，计算效率得到提高，计算成本也得到了降低。

（3）在时步长度Δt一定时，DQM求解动力反应的计算精度和数值稳定性与时步分段数m有关。排除一些失稳飘移的情况，一般m越大，计算精度越高，但计算量也越大。综合考虑，对于均匀网格离散方案，实际计算时取m＝10相对较优。

（4）使用DQM进行实际结构动力反应分析时，时间步长Δt可选为动力荷载的等效周期，然后将各时步等分成10段来计算，这样可获得较满意的计算结果。