第1章
索缆动力学中的动刚度法基础

本章为本书的理论部分，首先介绍了动刚度法和Wittrick-Williams算法（W-W法）。这两种算法是本书拉索动力分析理论的基础。同时，笔者在此基础上对原DSM进行了改进，使原DSM得以扩展至阻尼结构中，同时避免了复超越方程的求解。另外，通过提出“假象结构法”使W-W法得以应用于复杂梁式结构中，弥补了原W-W法只能应用于简单或特殊结构中的局限性。同时，整个求解过程仍是解析的。

1.1 动刚度法

动力学问题关乎工程结构的工作性能和制造成本，与静力问题具有同等重要的意义。无论是设计阶段、建造期间，还是结构服役阶段，分析结构的动力特性、预测特定激励下的动力响应，对工程结构的设计、施工和运营期安全保障都具有重要的作用 ^［33-38］ 。寻找精确和高效的动力分析方法一直是工程结构动力学问题研究的基础和核心内容。在长期探索中形成了两种不同的工作方向：一是将工程结构离散化，采用集中质量法、Ritz法或有限元法等 ^［39-41］ ，得到结构动力学问题的数值近似解；二是基于结构振动微分方程的分析思路。前者的分析精度依赖于离散化的方式和程度，效率则以相反的模式依赖于离散化程度。后者由于出发点和基本分析过程是闭式的，因此分析精度在推导过程中基本不损失，对一些简单结构而言，甚至会得到精度完全不损失的解析解。动刚度法（dynamic stiffness method，DSM）正是这种求解思路的产物。该方法是Koloušek于20世纪40年代提出的 ^［42］ ，是一种高效、精确和稳定的结构动力分析方法。动刚度是描述结构控制微分方程的另外一种表现形式，其解答在结构上任意点处均满足控制微分方程，且在结点处满足所有的平衡条件、位移协调或约束条件。

时至今日，动刚度法已经被广泛应用于各类工程学科中 ^［43］ ，用来解决结构动力特性分析、波传播（wave propagation）及屈曲等动力分析问题 ^［44］ 。一些典型的工程领域应用包括：①Langley ^［45］ 应用动刚度法研究了航天器板结构的振动特性及振动能量分布问题。② Wittrick W.H.等 ^［46］ 根据板结构动刚度法开发了VIPASA程序，该程序被集成于NASA的分析程序PASCO中，成为COSMIC ^［47］ 软件中的一部分，被广泛应用于航空航天工业。③周平等 ^［48］ 采用动刚度法计算了船体振动的固有频率，成为动刚度法在船舶结构动力分析的典型案例。另外，在土木工程结构领域，尤其是对动力荷载较为敏感的预拉结构的动力分析领域 ^［49-51］ ，动刚度法也得到了应用，其在模态分析中的精确性成为这类研究的突出亮点 ^［52-54］ 。

一般来说，动刚度法是基于边界条件约束下结构控制微分方程的“强”形式解而建立的。它从系统频域的控制微分方程出发，通过求其通解可以获得系统的精确形函数，进而获得系统的动刚度形式的平衡方程。其中，所建立的动刚度矩阵是频率的闭式函数，无须依赖结构的离散化，即可用于计算任意阶次的模态，其计算效率和精度均好于有限元法。动刚度法的这一性质使其能适用于任意边界条件和频率范围，尤其在高精度和高阶模态求解时是强有力的。它能准确获知结构在更宽频范围内的动力特性，而不以牺牲计算效率为代价。同时，借助有限元法划分网格和单元集组方法，动刚度法还可扩展应用于复杂工程结构。

动刚度法在实施中主要有两个关键问题：一是精确动刚度矩阵及频率方程的建立；二是频率方程的精确求解。

1.1.1 基本理论

动刚度法的关键之一是获得结构的动刚度矩阵。为此，首先需要建立结构的振动微分方程。这部分工作可以借鉴相应的文献或经典理论来实现。一般来说，无阻尼结构的自由振动可用以下控制微分方程加以描述 ^［55］ ：

式中 L——微分算子；

u——对应的位移向量。

通过引入简谐运动假设，则位移u可表示为：

式中 U——位移幅值；

ω——圆频率（rad/s）；

t——时间；

i＝ 。

将式（1.2）代入式（1.1），可消去与时间相关的项，得到频域控制微分方程：

式中 L ₁ ——微分算子。其解可通过以下形式获得

式中 C——常数向量；

A——与频率有关的方阵。

将式（1.4）代入位移边界条件，可得节点的位移向量δ：

式中 B——将式（1.4）代入位移边界条件后从A中获得的方阵。

另外，结合力的边界条件，同样可以得到节点力和向量C的关系：

式中 F——节点力向量；

D——与频率有关的方阵。

通过式（1.5）和式（1.6）消去常数向量C后，可得

式中 K——所求的动刚度矩阵，K＝DB ^-1 。

式（1.7）即以动刚度形式给出的结构动力平衡方程。结构自由振动时，动力平衡方程变为

式（1.8）是一个齐次方程组，通常情况下，为了得到非平凡解，要求：

式中 ｜·｜——行列式的值。

满足特征方程式（1.9）的频率即所求结构的自振频率，亦即丨K（ω）丨的零点频率为结构的频率。这就是动刚度法求解结构模态频率的基本思路。必须指出，尽管形式上与有限元模态参数求解方法类似，但由于动刚度矩阵的特点，如果不做特别的控制，该方法会有存在遗漏模态和增加虚假模态的可能。下面对动刚度矩阵特点做出详细说明。

1.1.2 动刚度矩阵特性

当结构由于几何或材料不连续需要特别划分单元时，各单元在局部坐标系下的动刚度矩阵K ^e 仍可按照式（1.1）～式（1.8）得出，而整体坐标系下的动刚度矩阵K ^e 则通过常规的坐标变换即可得到。在此基础上，仿照静力直接刚度法的集组流程（实质是按位移法叠加各单元的贡献），可形成任意复杂结构的整体动刚度矩阵K。

从动刚度矩阵的组建流程可以看出，动刚度矩阵是由控制微分方程和边界条件解的“强”形式建立的。不同于有限元法中的静力刚度矩阵，其特点如下：

（1）K和K ^e 同时包含质量和刚度的信息，矩阵中的各元素，即刚度系数通常不再是ω ² 的线性函数，而是ω ² 的超越函数。

（2）K ^e 是对称矩阵，但不再是奇异矩阵，且当ω→0时，动力问题的刚度矩阵K和K ^e 趋于静力问题的单元刚度矩阵。

（3）根据Leung ^［56］ 的研究，动刚度矩阵K（ω）与连续形式的质量矩阵M（ω）之间满足：

（4）若频率ω等于单元固端频率（称单元端点固定情况下单元的自振频率为固端频率，记为ω _F ）时，某些单元刚度系数就变得奇异，即 ＝∞。同理，当结构的频率与其某个单元的固端频率相同时，整体刚度矩阵K中会出现无穷大元素，即K _ij ＝∞，此时K是无意义的。

（5）丨K（ω）丨不仅在零点变号，在∞点也可能变号，因此根据正负号来判断零点和求零点的方法可能会漏根或增加虚假根。同时，当结构的某阶模态ω恰好等于某个单元的固端频率ω _F 时，特征函数丨K（ω）丨为∞。这意味着用式（1.9）不能求出所有的结构频率。

（6）由式（1.9）求出的频率是非平凡解，但结构可能存在δ＝0的振型。该振型没有结点位移，只有单元内部的“泡状”振动（图1.1）。这说明如果只求非平凡解δ≠0的频率可能会丢根。

由动刚度矩阵的以上特点可以看出，若要求得精确的频率和振型，就需要一个可靠的频率方程的解法，使其能够处理δ＝0和δ≠0两种情况，不遗漏地准确求出前n个根。下一节介绍的Wittrick-Williams算法就是这种算法。

1.1.3 存在的问题

动刚度法是基于控制微分方程和边界条件解的“强”形式建立的。得益于闭合形式的动刚度矩阵，动刚度法既保证了高精度，又节省了大量的计算时间 ^［57］ 。与之形成对照的是，传统的有限元分析法在模态分析时，其高精度是借助细化网格或采用高阶差值函数来取得的。因此，必然会以牺牲计算效率为代价。

Lee ^［58］ 在其论著中分析了有限元法、动刚度法、谱元法等的建模理论并指出：①有限元法是基于系统运动微分方程解的“弱”形式，通过近似的形函数表示结构的位移场，通过变分法推导出结构的单元刚度和质量矩阵。由于计算结果依赖于单元划分的质量和数量及形函数的近似程度，同时也只能得到数值近似解，因此有一定的误差。②由于有限元的形函数采用的是与频率无关的插值函数，未能考虑结构频率对结构变形的影响，从而使有限元法不利于计算高阶模态频率。③有限元法在计算高阶模态或提高计算精度时需要借助细化网格，这无疑会增加计算量并减小计算效率。而动刚度法除了结构的几何或材料不连续处需要特别划分单元外，无须对结构离散，因而能够保证在给出精确解的同时具有很高的计算效率 ^［59］ 。

虽然动刚度法是一种精确解法，但其高效性和准确性是通过增加形函数对实际结构的逼近程度和频率方程求解难度来实现的。结构准确的形函数可通过两个方面的努力来取得，一是需要确保控制微分方程（组）的通解能够解析给出；二是确保结构的频率方程的求解不遗漏根、不增加虚假根。这两个条件带来的数学上的困难，在一定程度上限制了动刚度法的应用范围。为了实现频率方程的精确求解，就需要引入W-W法。

1.2 W-W法

结构模态频率计算是结构动力分析的核心之一，它可以通过求解动刚度法建立的结构特征频率方程［式（1.9）］来解决。在早期的研究中，一些数学技巧被用来解决特定的简单结构的频率方程，如Timoshenko ^［60-61］ 和Huang ^［62］ 的工作。为了求解复杂系统频率方程，Cheng ^［63］ 、Li ^［64］ 和Henshell ^［65］ 等提出了“图解法”，实质上是通过零根搜索法来确定结构的各阶模态频率。这类方法不仅计算量大，而且也存在缺点，难以区分距离非常靠近的两阶模态频率，难以区分零点和正负无穷处的变号，可能导致漏根或增根 ^［66］ ，从而丢失模态或混入虚假模态。Muller搜索法 ^［67］ 就属于此类方法。

为此，Wittrick和Williams在20世纪70年代提出的W-W法 ^{［66，68-70］} ，用于分析结构自由振动和屈曲问题 ^［71］ 与波传导问题 ^［72］ 。由于W-W法适用范围广，因此此后针对该方法的应用和改进研究成为结构动刚度法的热点话题。Zhong等 ^［44］ 利用Rayleigh定理和Sturm序列的性质，对W-W法进行了推广，在理论上解决了频率方程求解的漏根问题，使得困扰动刚度法多年的瓶颈问题得到解决，推进了动刚度法理论的发展和工程应用。

实际上，W-W法并不直接用于频率计算，而是通过计数得到结构小于试探频率ω*的模态频率的个数，从而确定任意阶次频率的上下界，最后再结合二分法或牛顿法求出任意精度的频率解。由于该方法可使动刚度矩阵保持Sturm序列特性，能保证方程求解过程中不遗漏根，其稳定性和正确性也在理论上被证明。这是大多数解析法和近似解法所不具备的优点 ^［69］ 。

1.2.1 基本理论

先引入Rayleigh定理：

设某结构的自振频率按升序排列为ω _i （i＝1，2，3，…），对该结构施加一个约束，得到新的频率 （i＝1，2，3，…），则有ω _i ≤ ≤ 。

Rayleigh定理揭示了单自由度约束对结构频率的影响：当结构某自由度被约束后，结构低于某试探频率ω*的频率数保持不变或减少一个。这一推论是W-W法的重要理论基础。应用W-W法计算结构模态频率问题的关键步骤是确定结构的模态频率计数J，J表示结构小于ω*的模态频率的个数。J一般难以直接求解，但可通过另外两种计数之和来表示，即

式中 J _K ——与结构整体动刚度矩阵有关的计数，等于三角阵K ^Δ （ω*）主对角线上负元素的个数；

K ^Δ （ω*）——将结构的动刚度矩阵K（ω*）进行高斯消元后形成的上三角矩阵；

s｛ ｝——计数符号；

J ₀ ——结构的固端频率计数，它等于小于试探频率ω*的结构固端频率 个数。

当结构因材料或截面不连续需要被分为若干个单元e时，仿照结构固端频率计数J ₀ 的定义，若用J _m 表示单元的固端频率计数，则两者间的关系可表示为：

有关W-W法的详细证明过程可参阅 ^{［66，68-69］} 。在具体实施W-W法时，在频率ω*值给定时，动刚度矩阵K就随之决定，进而可以很容易地算出J _K 。而J ₀ 的计算则和具体的结构形式有关 ^［73］ ，是W-W法的难点和关键。能否成功求出J ₀ ，将决定动刚度法能否适用于分析对象。对于复杂的结构来说，由于其固端频率通常难以解析给出，因而J ₀ 的表达式也难以解析给出。这成为W-W法难以应用于复杂结构中的主要原因 ^{［74，75］} 。

按照求解方式的不同，可将J ₀ 的解法分为直接法、间接法和数值法。其中，直接法通过求解结构的固端频率方程，对固端结构小于试探频率ω*的模态频率进行计数得出J ₀ 。间接法避免了对J ₀ 的直接求解，且整个求解过程仍是解析的，因而适用范围更广。数值法通过常微分方程求解器求解结构的控制微分方程，得到单元动刚度矩阵的数值解，进而通过增加各单元的内部结点试算出结构的固端频率计数J ₀ 。

1）直接法

直接法求解有两个前提：一是动刚度矩阵中的各元素能够解析表达；二是固端结构的频率方程能够解析表达。这两个条件极大地限制了直接法只能在少数几种简单结构中使用；对于稍复杂一点的结构，就需要对结构进行简化处理。目前还没有直接应用动刚度法解决复杂结构的案例见诸报道。为了将W-W法推广至更为复杂的结构，其关键的瓶颈在于J ₀ 的求解方法问题，即需要研究一种针对复杂结构的J ₀ 求解方法。这通常需要借助间接的方法来实现。

2）间接法

该方法求解过程仍保持解析的特点，不直接对J ₀ 进行求解，而是采用迂回的方式间接计算得到J ₀ 。正因为如此，间接法的适用范围比直接法更广。其主要思路是：通过改变原结构的边界条件，构造一个非固结边界条件的中间结构，由于仅改变了边界条件，该“中间结构”与原结构具有相同的固端结构和单元刚度矩阵，从而具有相同的固端频率计数J ₀ 。在此条件下，先计算“中间结构”的整体动刚度矩阵 （ω*），对 （ω*）进行高斯消去成上三角矩阵后，得到其主对角线上负元素的个数 。然后求解“中间结构”频率方程，得到小于ω*的固有频率个数 。最后利用关系式J ₀ ＝ - 计算得到“中间结构”的固端频率J ₀ ，从而也就间接地求出原结构的固端频率计数 ^［76-77］ 。

3）数值法

严格来说，数值法实质上也是间接法的一种，该方法最早是由袁驷和叶康生等 ^［78］ 提出。它通过常微分方程求解器COLSYS求解结构的控制微分方程 ^［79］ ，可以得到任意复杂结构的动刚度矩阵K的数值精确解 ^［80-81］ 。由于COLSYS是一个基于高斯样条配点法且具有自适应网格求解功能的求解器，这使其在计算动刚度矩阵时具有数值精确解的精度，同时还能够将一组微分方程的求解结果保存起来供随时调用。COLSYS的这些特点和功能保证了其在计算结构的动刚度矩阵K和固端频率计数时的高效性和可靠性。

虽然数值法能够求解任意复杂结构的固端频率计数J ₀ ，但由于其刚度矩阵及固端频率计数的求解过程均是通过常微分方程求解器完成，因此冗余计算量大，计算效率不高。同时，工程结构几何和材料的不连续因素往往导致数值法的计算结果不稳定，甚至可能得到错误结果。

1.2.2 计算流程

1）模态频率的计算

当求出结构在试探频率ω*下的固端频率计数J ₀ （ω*）后，即可按照式（1.11）求出结构小于试探频率ω*的模态频率的个数J（ω*）。此后若想计算结构的第i阶模态频率ω _i ，就只需首先粗略地确定出该阶模态频率的上下界ω _l 和ω _u ，使满足

则结构的第i阶模态频率满足ω _i ∈（ω _l ，ω _u ）。此后可采用二分法、牛顿法等数值方法，通过不断调整上下界的值来逼近真实频率ω _i 。当满足ω _l -ω _u ≤Tol·（1＋ω _u ）时，可以得到允许误差范围Tol内的模态频率。

需要指出的是，W-W法最初是采用二分法求解模态频率 ^［66］ ，然而二分法求解时收敛速度较慢。Simpson ^［82］ 采用牛顿法改进了其收敛速度，数值算例表明相同条件下，牛顿法的求解时间仅为二分法的一半左右。但该方法首先需要将动刚度矩阵对频率进行求导，这就增加了计算的复杂性。为了避免对刚度矩阵求导，Williams和Kennedy ^［83］ 提出了一种可靠的多重行列式抛物线插值法，该方法比Simpson提出的牛顿法节省10％的计算时间。齐朝辉和Williams ^［84］ 还提出了以能量范数为标准的四阶龙格库塔法，从而很大程度上提高了收敛速度，且精度可以用绝对误差来衡量，但该法同样需要对动刚度矩阵进行求导。袁驷等 ^［85］ 在二分法的基础上提出了一种具有二阶精度的牛顿法，该方法通过求解一个广义特征值问题

然后，利用每次求出的特征值μ外插出一个更接近真实值的结果ω _μ 。上式中，ω _a 表示由二分法计算出的具有一阶精度的近似频率值，K _a ＝K（ω _a ）， ＝dK（ω _a ）/dω，w _μ ＝ω _a -μ。

2）模态振型的计算

振型的求解过程实际上是问题式K _a Δ＝ 的广义特征向量的求解问题，求解该问题的一个自然的选择是逆幂迭代法。考虑如下的标准广义特征值问题：

逆幂迭代法的步骤是：

（1）给定初始向量Δ ^（0） 。

（2）求解新的特征向量： ＝K ^-1 MΔ ^（k） ，k＝0，1，…。

（3）找绝对值最大元素：设第i个元素的绝对值最大，即 。

（4）估计特征值： 。

（5）特征向量归一化： 。

（6）检验是否满足用户的精度要求或超过了允许的迭代次数k≥k _max ，如“是”则停止。

（7）否则，回到第（2）步。

对于第（6）步的精度要求，可以用特征向量来控制，例如，

式中 Tol——用户预先指定的误差限。

依据逆幂迭代法，袁驷等利用W-W法，在牛顿法的基础上对迭代过程加以有效地保护和引导，以确保步长方向的正确和大小的适当，最终将牛顿法转变为一个全局大范围收敛的方法。文献［86］详细介绍了对牛顿法的保护和引导过程，并给出了完整的导护性牛顿法的算法步骤。该算法可以高效、精确地求解结构的振型，与W-W法已形成了“标配”。

从上述求解流程可以看出，W-W法能够很好地解决频率方程的求解问题。同时相较于传统的零根求解法，该方法在求解频率方程时主要有以下优势：①理论完善、方法简单，同时结果可以达到计算机允许的任意精度；②不丢根，单根、重根均可处理；③可以直接求解第k阶频率，而无须预先解出前k-1阶频率值。

1.2.3 存在的问题

从上文可以看出：W-W法多用于规则或简单结构，对于复杂结构由于固端频率计数J ₀ 求解上的困难，其固端频率通常难以得到，从而导致了难以得到J ₀ 及其解析表达式。这使得W-W法在复杂结构中难以推广应用。为了继续在复杂梁式结构中应用W-W法来得到高精度、不漏阶的模态频率，就必须解决J ₀ 的求解问题。有鉴于此，1.3节将在总结归纳已有的J ₀ 求解方法的基础上，提出一种用于复杂梁式结构的固端频率计数J ₀ 求解的改进W-W法。并结合算例验证该方法的准确性和有效性，使W-W法应用于复杂索缆结构的动力分析中成为可能。

1.3 适用于复杂梁式结构的改进W-W法

1.3.1 基本理论

多数复杂工程结构均可看作由若干种基本结构组合而成。这些基本结构的动刚度通常可以解析表示，如能利用这一特点求解J ₀ ，将能很好地克服数值法的上述缺陷。复杂梁式结构是一种常见的工程结构形式，针对该类结构的动力特性分析问题，笔者提出了一种改进W-W法——“假想结构法”。该方法通过构造一个与原结构边界条件不同的基本结构，称为“假想结构”，用以间接求解原结构的频率方程。求得“假想结构”的动刚度矩阵后，通过解其频率方程可以得到“假想结构”的固有频率计数 和刚度矩阵特性计数 ，从而间接计算出原结构的固端频率计数J ₀ 。根据J ₀ 可以确定出原结构各阶模态频率的上下界，为二分法或牛顿法求解频率方程建立基础，最终解决复杂梁式结构的频率和阵型计算问题。

图1.2给出了直接法和假想结构法的求解过程比较。由图可知，由于原结构和假想结构共享同一个固端结构，因此两者具有相同的固端频率计数J ₀ 。但由于直接法在面对复杂结构时难以奏效，因此采用“假想结构法”，根据“假想结构”的固有频率计数 和动刚度矩阵特性计数 之差来计算原结构的J ₀ 。其中， 可通过对假想结构的动刚度矩阵进行上三角变换后，对主对角线上为负的元素计数得出，其求解过程与J _K 完全一致。而假想结构的固有频率计数 的求解问题相对困难，因此是研究关键之一。

直接法和假想结构法固端频率计数求解过程的区别可分别由式（1.16）和式（1.17）表示：

式中 p ₁ ——原结构的材料和几何特性参数，如抗弯刚度EI、质量m、长度l等；

p ₂ ——固端结构的边界条件，全部为固结；

——假想结构的边界条件，一般是简支边界条件；

ω——固端频率；

——试探频率，需要预先指定；

——“假想结构”的固有频率计数，指小于ω*的“假想结构”模态频率的个数；

——“假想结构”的动刚度矩阵，亦即对原结构应用边界条件后得到的动刚度矩阵；

——代表用高斯消元法将K消成上三角阵K ^Δ 后，其主对角线上元素为负的个数；

count｛ ｝——计数符号，表示满足条件的自变量的个数。

表1.1是对J ₀ 的三种不同求解方法的对比，从中可以看出，假想结构法比其他已有方法优势明显。

1.3.2 计算流程

假想结构法的具体求解步骤如下：

（1）构造假想结构。假想结构通常是将原结构的边界铰接后得到的。

（2）构造位移函数。位移函数是由满足端部边界条件的谐振函数形式构造的，可自动满足假想结构的边界条件。

（3）求解假想结构的频率方程和固有频率计数 。将构造出的位移函数代入原结构的控制微分方程中，得到假想结构的动刚度 和频率方程。该方程一般是关于频率的一元高次方程，运用相应的求根公式即可求出假想结构的各阶模态频率，最后计数得到小于试探频率ω*的固有频率个数 。

（4）求解假想结构的 。对假想结构的动刚度矩阵 进行高斯消元，得到试探频率下的上三角阵 ，对其主对角线上为负的元素进行计数求出。

（5）求解J ₀ 。根据J ₀ ＝ 即可求得原结构的固端频率计数J ₀ 。

（6）计算原结构的J。它可由J＝J ₀ ＋J _K 求出。

（7）确定第k阶模态频率的上下界ω _u 和ω _l ，使J（ω _u ）≥k和J（ω _l ）＜k。

（8）以二分法、牛顿法等数值求根方法求解原结构的频率方程，得到上下界ω _u 和ω _l 内的第k阶模态频率。

1.3.3 准确性验证

以Bernoulli-Euler梁的完曲振动为例，对作者提出的改进W-W法进行验证 。

1）J的计算

对于轴向力作用下的Bernoulli-Euler梁，其横向自由振动微分方程为

式中 EI——梁的弯曲刚度；

m——梁的单位长度线质量；

φ——梁的位移函数的幅值函数；

P——梁的轴向力（以梁受拉为正）。

以两端简支的欧拉梁单元作为假想结构，如图1.3所示，其边界条件有

故边界条件也应满足以下通解形式：

将式（1.20）代入式（1.18）后可得

上式具有非平凡解的条件为

式（1.22）即假想结构的频率方程，其中 ， ， 。该方程是一元二次方程，其频率解为

由上式可以看出，一个n值对应一个频率解，因此通过改变n的值即可找出假想结构小于ω*的模态频率个数 。

2） 和J ₀ 的计算

文献［88］给出了简支边界条件下的动刚度矩阵 为：

其中

其中

根据假想结构对应的刚度矩阵 ，即可求出假想结构的 。最后通过J ₀ ＝ 即可得出轴力作用下Bernoulli-Euler梁的J ₀ 。

文献［88］给出了轴力作用下Bernoulli-Euler梁J ₀ 的计算公式

式中 floor（x）——求小于x最大整数的函数；

sgn（x）——获取x符号的函数，当x＞0、x＜0和x＝0时，sgn（x）分别为＋1、-1和0。

下面对本书提出的改进W-W法进行验证。设m＝EI＝1，P＝2，梁长l＝1，频率ω*取不同值时，两种算法所得结果的对比见表1.2。

由表1.2可知，对试探频率取不同值时，本书的计算结果均与文献给出的结果吻合。这在一定程度上说明了本书提出的假想结构法对Bernoulli-Euler梁的适用性，也证明了其在以梁模型为基础的拉索动力分析中的有效性。文献［87］以考虑弯扭耦合的Timoshenko梁及曲梁等更复杂的梁式结构的固端频率计数J ₀ 进行了计算，进一步说明了该方法的准确性和普适性。

1.3.4 改进W-W法优势

（1）可以应用于动刚度矩阵能够解析表示的复杂结构，且求解过程中的中间变量均可表示为初等函数形式，因此得到的J ₀ 是一种闭式解。

（2）给出了假想结构频率方程的一般化建立方法，保证了对复杂梁式结构的适用性。

（3）除了频率计算外，其余计算过程均以解析的方式完成，因此比数值法具有更高的计算效率高。

1.4 W-W法在弱几何非线性结构中的扩展

由于动刚度法适用于线性系统，当考虑拉索垂度效应后系统将具有非线性特征，此时系统动刚度矩阵不再严格对称（参见第3章），也就无法直接应用W-W法进行频率计算。然而，对于小垂度拉索系统，由于其弱非线性特点，垂度对各阶模态频率的影响有限，尤其对于高阶模态的影响几乎可以忽略 ^［31］ 。

基于这一特点，可以做出如下假设：相比于无垂度系统（即h＝0），小垂度索缆系统的各阶模态频率ω _sag 与无垂度系统的模态频率ω ₀ 具有一一对应的关系，即第n阶模态频率 将出现在 的邻域内，这一假设也是符合实际工程情况的 ^［89-90］ 。

1.4.1 计算步骤

基于上述假设，本书结合改进W-W法给出了一种适用于小垂度拉索体系的频率方程解法，具体计算步骤如下：

（1）分别计算无垂度系统和有垂度系统的系统动刚度矩阵K（ω）和K _sag （ω）。

（2）应用改进W-W法计算无垂度系统的第n阶模态频率 ，从而“锁定”小垂度系统的第n阶模态频率 可能落在的区间，并以 为初值在其邻域开始搜索 。

（3）将 代入动刚度阵K _sag 并计算丨K _sag （ ）丨。若丨K _sag （ ）丨≤Tol，则认为 ＝ 。这说明该垂度对系统该阶模态频率的影响很小。

（4）若丨K _sag （ ）丨＞Tol，则需要进一步判断搜索方向，即判断abs（丨K _sag （ ）丨）和abs（丨K _sag （ ＋Tol）丨）两者的大小。其中，abs（·）为绝对值。若abs（丨K _sag （ ）丨）＞abs（丨K _sag （ ＋Tol）丨）则不断增加初值直到相邻两个搜索步内有丨K _sag （ ＋（m-1）Tol）丨×丨K _sag （ ＋mTol）丨＜0。此时停止搜索并输出 ＝ ＋（m-0.5）Tol，其中m为计算次数。

（5）若abs（丨K _sag （ ）丨）＜abs（丨K _sag （ ＋Tol）丨），则不断减小初值直到相邻两个搜索步内有丨K _sag （ -（m-1）Tol）丨×丨K _sag （ -mTol）丨＜0。此时停止搜索并输出 ＝ -（m＋0.5）Tol。

在上述求解步骤中，首先通过W-W法对小垂度系统的各阶模态频率进行了“定位”，大大缩小了搜索范围，从而提高了计算效率且避免了漏根的可能。

1.4.2 准确性验证

需要说明的是，本书提出的方法及研究对象是垂跨比e＜1/8的小垂度拉索，不考虑大垂度及大变形的情况。为了说明本书所提出的改进W-W法在计算小垂度拉索时的准确性，以文献［90］中给出的两个拉索为例，其参数见表1.3，文献［90］、原W-W法、改进W-W法的计算结果见表1.4。

从表1.4可以看出：由于原W-W法无法计入垂度影响，其在计算低阶模态时容易产生较大误差，而改进W-W法所得结果与文献结果更加接近，从而说明了改进W-W法的计算精度。需要说明的是，对于偶数阶模态而言，由于拉索振型具有反对称特点，其附加索力h＝0 ^［90］ 。因此，垂度将不会引起模态频率的改变，此时三种方法的计算结果一致。对于奇数阶模态，垂度对低阶模态的影响更为显著。

为了进一步说明改进W-W法在频率方程求解时的准确性，以及以文献［90］为代表的一类解法在频率方程求解时存在的漏根可能性，仍以上述拉索为例，表1.5给出了文献［90］、有限元解 ^［91］ 及张紧弦解。

从表1.5可以看出：相比文献方法，本书提出的方法与文献结果及有限元解吻合得很好，说明本书提出的方法在进行频率方程求解时具有更高的精度。此外，对于拉索2的第三阶频率，本书方法与有限元解一致，同时参照弦理论可以判断出文献算法遗失了第三阶模态，所得频率实为该索的第四阶模态。

从上述分析可以看出，改进W-W法较文献不仅具有更高的计算精度，还能够有效地避免在求解超越频率方程时的漏根可能性。此外，采用有限元法求解结构的高阶模态特性时往往需要加密单元划分，这将导致计算效率的降低且计算精度难以保证。与有限元法不同的是，改进W-W法在对结构进行动力分析时，仅需在约束位置及几何不连续处增设节点。以本书讨论的拉索为例，全索只需建立一个单元即可求解出各阶模态频率，因此改进W-W法比有限元法等数值计算方法具有更高的计算效率。

1.5 扩展动刚度法

当考虑结构的阻尼特性时，如裸索系统的内阻尼或复合拉索系统填充层材料的阻尼特性时，系统的频率方程将变为复超越方程，从而使W-W法不再适用。为了继续使用W-W法求解频率，就需要避免复超越方程的产生，因此需要对动刚度法进行改进。鉴于此，本书提出了一种适用于线性阻尼系统的扩展动刚度法（EDSM） ^［92］ 。

扩展动刚度法通过在拉氏域内对结构位移函数进行变量分离，建立了阻尼和无阻尼系统统一形式的动刚度矩阵和频率方程；并在此基础上，给出了复杂结构阻尼比的确定方法，从而将难以求解的阻尼频率转换为求解系统阻尼比及对应的无阻尼频率。在应用扩展动刚度法求解时，其计算流程与动刚度法基本一致，只需要在分离变量时将结构的位移向量和荷载向量从时域转换到拉氏域，进而得到系统在拉氏域的控制微分方程及频率方程丨K（λ）丨＝0，其中λ＝α＋iβ，i＝ 。该处理方式不仅能有效地避免复动刚度矩阵的产生，还不会额外增加频率方程的求解困难。此外，在分析阻尼系统时，可以非常方便地考虑时滞阻尼及黏滞阻尼的影响。特别是，对于无阻尼系统，有α＝0，λ＝iω。

为了介绍扩展动刚度法在分析不同阻尼系统时的计算步骤，以欧拉梁为例，分别针对时滞阻尼和黏滞阻尼两种类型进行分析。

1.5.1 基本理论

考虑一个匀质等截面欧拉梁的弯曲振动，其无阻尼自由振动方程可写为 ^［93］

式中 EI——梁的抗弯刚度；

m——梁的线质量；

u（x，t）——梁的横向位移函数。

假设解的形式为u（x，t）＝ ，则式（1.27）可化为

令a ₀ ＝1，b ₀ ＝ ，则上式可改写为更一般的形式，如下

设上式的通解形式为U＝Ae ^{κ x} ，则微分方程式（1.29）的特征方程为：

其中， ＝- /EI；p ₀ 与波数κ有关，本书称其为波数因子。记该特征方程的四个根为κ ₁ ～κ ₄ ，则式（1.29）的通解可表示为

其中A ₁ ～A ₄ 为待定系数。由此可确定位移函数的具体形式。

然后，结合位移边界条件

和结点力边界条件

可以直接写出以动刚度形式表达的结点力平衡方程

其中，q＝｛q ₁ ，q ₂ ，q ₃ ，q ₄ ｝，Q＝｛Q ₁ ，Q ₂ ，Q ₃ ，Q ₄ ｝。上式中动刚度阵K是λ ₀ 的函数，而非传统意义上的频率ω的函数。

故可将上述动刚度阵表示为如下简洁形式。该形式是以波数因子p ₀ 表示的函数形式，与前述的函数形式K（·）有区别。为了方便，仍记作K（·），应该有K（λ ₀ ）＝K（p ₀ ）＝K（p ₀ （λ ₀ ））。

矩阵中各系数为

1）黏滞阻尼系统

若考虑材料的黏滞性阻尼时，上述欧拉梁的控制微分方程变为：

式中 γ和c——与应变率和速度相关的阻尼系数。

仍假设u（x，t）＝ ，则式（1.34）化为

其中a _vd ＝1＋λ _vd γ，b _vd ＝ ＋λ _vd c。假设U＝ ，可得上式的特征方程为

这里波数因子 ＝- ，其中 的下标代表viscous damping。故而微分方程式（1.35）通解可表示为

其中 ～ 是式（1.36）的四个根。

对比式（1.30）和式（1.36）可以看出，当考虑材料的黏滞性阻尼特性后，系统与对应无阻尼系统的特征方程和通解具有完全相同的形式，唯一的区别在于指数κ和 不同。由此可见，在黏滞阻尼系统的动刚度矩阵 中，各元素的表达式应与无阻尼系统的动刚度矩阵 ，因此只需用 ＝- 替代 即可。

2）时滞阻尼

当结构阻尼类型为时滞阻尼时，微分方程式（1.29）化为

式中 ＝m（1-iω ₀ η _m ）， ＝E（1＋iω ₀ ηE）；η _m 和ηE是由结构阻尼引起的损失因子；a _hd ＝1＋iω ₀ ηE，b _hd ＝ （1-iω ₀ η _m ），其中ω ₀ 为无阻尼欧拉梁的频率。故式（1.38）的特征方程为：

式中波数因子 ＝- 。因此，式（1.38）的通解可表示为

由于通解的形式与无阻尼系统一致，因而时滞阻尼系统动刚度矩阵的形式与无阻尼系统K ₀ （p ₀ ）完全相同，只需将中间变量 改为 。其中p _hd 的下标表示hysteretic damping。

上述三种情况中，p ₀ 、p _vd 及p _hd 均为系统的波数因子，统一记作p；λ ₀ 、λ _vd 及λ _hd 统一记作λ。

由特征方程式（1.39）、式（1.36）及式（1.30）可看出，基于扩展动刚度法得出的阻尼系统和无阻尼系统的频率方程实质上是一个关于波数因子的特征值求解问题，其中p是λ的函数。因为p和λ具有一一对应的关系，当λ确定后，p也将唯一确定，因此满足丨K（λ）丨＝0的λ值也必然使得丨K［p（λ）］丨＝0。因而系统频率方程丨K（λ）丨＝0可用丨K（p）丨＝0来代替。

当通过丨K（p）丨＝0确定出p的值后，即可分别由p ₀ 、p _vd 及p _hd 的表达式确定出各系统的频率特征值如下

由上述分析可知，对于线性阻尼系统，扩展动刚度法首先通过对位移函数进行拉氏变换得到了与无阻尼系统形式完全一致的动刚度矩阵和频率方程。随后对统一频率方程丨K（p）丨＝0进行求解，得到了与频率相关的波数因子p（λ）。最后根据p与λ的关系即可确定出系统的频率特征值λ。直观起见，将扩展动刚度法的计算原理用图1.4表示。

可见，当采用传统动刚度法分析有阻尼系统时，由于阻尼项的存在使得动刚度矩阵的表达式中将不可避免地出现虚数单位i，因此频率方程通常难以求解。而扩展动刚度法在进行变量分离时由于采用了更为一般的复频率参数，使得动刚度矩阵的表达式中不再包含虚数单位，从而有效地避免了对复超越方程的求解。

3）与有限元中Rayleigh阻尼的等效

式（1.34）和式（1.38）中分别考虑了系统阻尼的影响。在传统的有限元法中，有多种阻尼模型。在通用有限元软件ANSYS中可以定义五种形式的阻尼，包括结构的Rayleigh阻尼、材料相关的阻尼、恒定阻尼比、振型阻尼和单元阻尼。其中Rayleigh阻尼是结构有限元分析最常用的阻尼模型。根据Rayleigh阻尼的定义及计算方法，Rayleigh阻尼模型考虑了两种形式的阻尼：

一是与质点运动速度成正比的阻尼力，即式（1.34）中的阻尼系数c。由有限元定义，单位体积上作用的与质点运动速度成正比的阻尼力可表示为：

式中 α ₀ ——比例常数；

ρ——材料常数；

｛f｝——位移函数向量；

［N］——形函数矩阵；

——单元结点位移向量。

则单元结点上的阻尼力向量可以表示为

式中 V——结构体积。

式（1.42）同样可写成

式中 ［c］——阻尼矩阵，可以表示为

式中 ［m］——单元质量矩阵。

二是与应变速率成正比的阻尼力，即式（1.34）中的γ。根据有限元法的定义，与应变速度成正比的阻尼力可写为

则单元的阻尼力向量可以表示为

上式亦可表示为

式中 ［c］——阻尼矩阵，可以表示为

式中 ［k］——有限元法中的单元刚度矩阵。

根据有限元法，阻尼矩阵可看作质量矩阵和刚度矩阵的线性组合，即

式中 α ₀ 和β ₀ ——结构Rayleigh阻尼系数。

比较式（1.44）、式（1.48）及式（1.36）中的阻尼系数可知，扩展动刚度法中结构与速度相关的阻尼系数c和与应变率相关的阻尼系数γ与Rayleigh阻尼系数α ₀ 、β ₀ 的关系为

由上式可知，通过扩展动刚度法中的两种阻尼系数与结构Rayleigh阻尼系数的关系，可以十分方便地考虑结构的Rayleigh阻尼。

1.5.2 系统阻尼比的计算

由图1.4可以看出，EDSM的优势在于能够使阻尼系统和无阻尼系统在波数因子p层面上的频率方程丨K（p）丨＝0具有统一性。这一特性使得阻尼系统的模态频率ω _D 可以通过无阻尼频率ω ₀ 来计算，而无须重新推导系统在拉氏域下的动刚度矩阵。其中，ω ₀ 可直接由无阻尼系统的频率方程丨K（ω ₀ ）丨＝0求出，而ω _D 与ω ₀ 之间的关系可通过系统阻尼比ξ确定。为此，下面对阻尼比ξ的计算方法进行介绍。以黏性阻尼系数为c的欧拉梁为例，其控制微分方程为

其对应的无阻尼系统的控制微分方程为

由1.5.1分析可知，阻尼系统的频率特征值λ _D 和无阻尼系统的频率特征值λ ₀ 均可由波数因子p确定，即对于同一阶模态应有

因而有

上式可重新改写为如下标准形式

其中ξ＝ 。式（1.55）可看作关于 的一元二次方程，其根为

其中，λ _D 的实部反映了系统运动的衰减幅度，而虚部 即有阻尼系统的固有频率ω _D 。有关无阻尼模态频率ω ₀ 的计算可按照改进W-W法求得。

由此可将有阻尼系统模态分析的流程归纳如下：

（1）对系统位移函数在拉氏域进行变量分离，分别计算阻尼系统和无阻尼系统控制微分方程的特征方程。

（2）由特征方程计算系统阻尼系统阻尼比ξ的表达式。

（3）计算无阻尼系统的动刚度矩阵。

（4）通过W-W法求解丨K（λ ₀ ）丨＝0，确定无阻尼系统的模态频率ω ₀ 。

（5）根据ω _D ＝ 计算阻尼系统的模态频率。

技术流程图如图1.5所示。

1.5.3 方法验证

为了验证扩展动刚度法的正确性，下面结合文献［94］中的案例进行说明。文献［94］研究了一类由黏弹性阻尼连接的简支双梁系统，该系统中两个梁的抗弯刚度、长度及质量完全相同，且忽略了双梁内阻尼。可见，该双梁系统是上文所述复杂双梁系统的一个特例，即E ₁ I ₁ ＝E ₂ I ₂ ＝EI，m ₁ ＝m ₂ ＝m，c ₁ ＝c ₂ ＝0时的情况，其发生自由振动时的控制微分方程为：

对该系统应用扩展动刚度法，为了获得该双梁系统的动刚度矩阵，首先需要将位移函数u ₁ 和u ₂ 进行分离变量。令u ₁ （x，t）＝φ ₁ （x）e ^{λ t} ，u ₂ （x，t）＝φ ₂ （x）e ^{λ t} 。其中，φ ₁ （x）和φ ₂ （x）分别为1号和2号梁的振型函数。假设两者的形式为

其中，A和B为与边界条件有关的待定参数。将式（1.58）代入式（1.57）中可以得该系统的特征方程组为

根据有阻尼系统和无阻尼系统具有统一形式频率方程的特点，若计无阻尼双梁系统的动刚度矩阵为K ₀ （κ ₀ ），则对于同一阶模态应有

由式（1.59）可得无阻尼双梁系统的控制微分方程为（c＝0）

再根据κ ₀ ＝κ即可构造出黏弹性双梁系统的阻尼比。将式（1.61）和式（1.59）中的二式分别作差后得

再用式（1.63）减去式（1.62）可得

上式可改写为以下形式

其中λ ₀ ＝iω ₀ ，ω ₀ 为无阻尼双梁系统的模态频率。

定义ξ＝ 为系统的阻尼比，则式（1.65）可进一步改写为

上式的解为

由式（1.67）可以看出，通过构造阻尼比，有阻尼双梁系统的模态频率可表示为与多自由度系统相似的形式。可见，阻尼双梁系统的模态频率ω _D 依然可通过无阻尼频模态频率ω ₀ 和阻尼比ξ来确定。

文献［94］中给出了该双梁系统第n阶模态频率的表达式

其中Ω _n 为连接层阻尼系数c＝0时对应模态频率， 为阻尼比，两者表达式为

不难看出，ξ与 的形式完全相同，从而验证了扩展动刚度法的正确性。值得说明的是，扩展动刚度法不仅可以用于处理两个梁的质量、刚度及阻尼系数不同的情况，还能很方便地解决除简支边界外的其他边界条件，且计算结果是精确的，因此能够很好地用于解决各类复杂工程结构的动力问题。

1.5.4 优势

（1）时滞性和黏滞性阻尼均可应用。

（2）基于该理论原W-W法可被扩展应用于复杂阻尼结构。

（3）该方法易于实施，且很好地保留了原方法的优点。

1.5.5 复合索缆体系

考虑连接层阻尼特性的复合拉索体系控制微分方程为［95］

式中 u _i （x，t）——系统第i片梁的横向位移函数；

E _i I _i 、 P _i 及m _i ——第i片梁的抗弯刚度、轴向力及单位长度线质量；

t和x——时间和空间坐标；

材料常数k、c和m ₃ ——黏弹性层的刚度、阻尼系数和单位长度线质量。

其中i＝1，2分别代表上梁和下梁。

仿照同样的流程，根据扩展动刚度法令u ₁ （x，t）＝Ae ^{κ x} e ^{α t} ，u ₂ （x，t）＝Be ^{κ x} e ^{α t} ，则可以得到阻尼系统的频率因子 和无阻尼系统（即c＝0时）的频率因子λ间的关系

定义系统阻尼比ξ＝ ，则可得如下标准形式

其解为 ＝-ξλ± ，从而可得阻尼系统的模态频率 ＝ 。其中ω为无阻尼系统的频率。

综上所述，当采用传统动刚度法分析有阻尼系统时，由于阻尼项的存在使得动刚度矩阵的表达式中将不可避免地出现虚数单位，因此频率方程通常难以求解。而扩展动刚度法在进行变量分离时由于采用了更为一般的复频率参数，使得动刚度矩阵的表达式中不再包含虚数单位，从而有效地避免了对复超越方程的求解。此外，该方法通过定义系统的模态阻尼比，无须再次推导阻尼系统的动刚度矩阵即可实现模态阻尼频率 的计算。 oBCaUOQsdoFQLtHghvwdCvvXDeRGigjF7IZKZ9GCzRq9oLkfwvZVgBjVoH4pA64a