悬臂梁结构的应变模态积分法与灵敏度分析

邢建伟， 郑钢铁

(清华大学航天航空学院， 北京 100084)

悬臂梁结构的应变模态积分法与灵敏度分析

邢建伟， 郑钢铁

(清华大学航天航空学院， 北京 100084)

应变模态振型决定了结构动应变分布规律。基于梁应变与曲率的正比关系，提出一种直接求解悬臂梁结构曲率模态的积分方法，适用于任意形式的刚度和质量函数。通过与数值和解析解比较模态频率与振型，验证了该方法的准确性。此外，该方法将连续体模态求解问题转化为有限自由度的广义特征值形式，进一步证明可以应用Nelson法求解该形式的特征灵敏度，获得了曲率模态振型对于结构设计参数的灵敏度显式表达。为了表现动应变变化，提出了曲率梯度模态的概念，其参数灵敏度可由曲率模态灵敏度获得。算例结果表明：该方法可以有效地优化结构参数、避免动应变集中。该方法同样也适用于固支-固支和固支-简支这两种边界条件。

应变模态； 曲率模态； 积分法； 灵敏度分析； 曲率梯度模态

引 言

分析激励作用下的动应力状态及其分布规律是工程结构动强度设计的关键，而应力又需要由应变通过本构关系间接得到。结构动位移响应可以采用位移模态综合法获得，同理也能够采用应变模态综合[1-2]直接获得动应变响应。应变模态振型是结构本身的固有特性，与激励载荷无关，决定了空间动应变的分布。工程师可以通过修改结构参数来优化应变模态振型，以改变动应变分布、避免动应变集中，从而提高结构的可靠性和使用寿命。

梁是许多工程结构的简化模型，尤其是悬臂梁模型，如机翼、高层建筑等。欧拉梁截面的应变与中性面变形曲率成正比，所以梁中性面曲率模态可以等效应变模态，下文所述的曲率与应变二者等同。曲率模态在损伤识别[3-8]领域有广泛的应用，但是这些文献一般都是先通过有限元法得到位移模态，再采用中心差分法近似、间接地获得曲率模态。据目前所调研的公开发表文献所示，还没有一种直接求解任意梁曲率模态的解析法或半解析法。

对于任意变截面非均质欧拉梁，其运动方程是以横向位移为变量的变系数四阶偏微分方程，如果能够得到位移模态的解析解，也可对该解析解求两次导数得到曲率模态。但是，文献[9-10]中的位移模态解析解仅适用于等截面均质梁。Naguleswaran[11]采用Bessel函数求解了楔形和锥形这两种特殊截面变化形式的欧拉梁位移模态解析解。Melnikov[12]在其著作里采用Green函数法求解任意欧拉梁位移模态，但是不同边值条件下Green函数的求取非常困难，这就限制了该方法的应用。Li[13]提出了一种求解任意梁剪切变形的精确方法，但是刚度函数与质量函数之中只能有一个是任意的。Elishakoff和Candan[14]给出了非均匀梁自振频率的封闭形式表达，但其要求密度和弹性模量都表示成多项式函数的形式。Huang等[15]提出一种求解轴向参数变化梁模态频率的积分法，相比Green函数法更加简便，且适用于多种刚度函数与质量函数。但是Huang在文中仅通过比较模态频率验证了方法的正确性，而没有分析模态振型，且据所调研公开发表文献，还没有学者将这种积分方法推广到求解曲率模态。

除了准确求解任意梁的应变模态，在结构设计阶段，工程师们更关心结构设计参数对应变分布的影响。模态振型灵敏度分析是连接结构设计参数与动应变分布的桥梁，也是前面所讲的优化动应变分布的关键课题，学术界在这方面的研究也很多。Sutler等[16]比较了4种典型模态振型灵敏度的计算方法：有限差分法、模态法[17]、改进模态法[18]和Nelson法[19]，结果表明Nelson方法最为有效、准确。Yu等[20]将Nelson方法得到的结果作为比较另外5种近似灵敏度求解方法的标准解。但是Nelson方法处理的是矩阵形式的特征值和特征向量灵敏度问题，如果利用该方法求解欧拉梁这种连续结构的应变模态灵敏度，还需要将问题转化为包含有限自由度的矩阵形式。

本文提出了基于悬臂欧拉梁结构的曲率模态积分法，将以曲率为变量的运动微分方程转化为积分方程形式，从而可以直接求解曲率模态。该方法适用于任意刚度与质量函数，通过与解析解及精细有限元结果的对比，验证了方法的准确性。此外，这种积分法能够将连续结构的模态求解转化为矩阵形式的广义特征值问题，同时基于矩阵分析，证明了采用Nelson法求解参数灵敏度的可行性。最后，计算了曲率梯度模态及其参数灵敏度，曲率梯度能够更加明显地反映应变集中。通过典型算例展示了曲率模态振型设计流程并验证了方法的有效性。

此外，通过公式推导与算例验证，证明该方法也适用于固支-固支和固支-简支两种边界条件。

1 应变模态积分法

任意截面非均质欧拉梁的无阻尼运动微分方程可写成

(1)

式中x为轴向坐标，y(x,t)表示中性面的横向位移；梁单位长度的质量记为m(x)，且m(x)=ρ(x)·A(x)，其中ρ(x)为梁材料密度，A(x)为梁的横截面面积；梁横截面的抗弯刚度记为S(x)，且S(x)=E(x)I(x)，其中E(x)为梁材料的杨氏模量，I(x)为梁横截面绕中心轴z轴的转动惯量；F(x,t)表示梁单位长度上的横向载荷，梁长L。

采用分离变量法，令

y(x,t)=Y(x)T(t)

(2)

ϖm(ξ)Y(ξ)=0,

0≤ξ≤1

(3)

对各向同性材料，仅考虑梁截面轴向的正应力，则应力应变关系为

σ(ξ)=E(ξ)ε(ξ)

(4)

式中σ(ξ)和ε(ξ)分别为轴向ξ处的正应力与正应变，应力不可直接测量，需通过公式(4)由应变换算得到，二者相互对应。

(5)

可见截面上任意点的应变均正比于该截面处的中性面曲率，所以可以用欧拉梁的中性面曲率表征应变，即下文中所述的应变与曲率二者等同。

1.1 曲率微分方程

由于曲率是梁横向位移的二阶导数，则有

式中r和s为内层积分变量。

对于悬臂梁结构，边界条件为

(7)

考虑到边界条件(7)，将式(6)代入到方程(3)中得

由高等微积分知识可知

(9)

则式(8)可化简为

(10)

上式即为以曲率Yc(ξ)为变量的自由振动微分方程。

将方程(10)等号两端在[0,ξ]范围内依次积分两次得

r)Yc(r)drdu=C1

(11)

r)Yc(r)drdu=C2+C1ξ

(12)

式中C1和C2为待定系数，u和υ为积分变量，式(12)进一步化简得

r)Yc(r)drdu=C2+C1ξ

(13)

将悬臂梁边界条件(7)代入到式(11)和(13)得

(14)

求得C1和C2后代回至方程(13)，最终得到

r)Yc(r)drdu=0, 0≤ξ≤1

(15)

至此，对于悬臂梁结构最初的以位移为变量的运动微分方程，通过积分方法转化为以曲率Yc(ξ)为变量、适用于任意弯曲刚度S(ξ)和质量函数m(ξ)的积分方程(15)。求解该方程可直接获得梁中性面曲率，进而确定结构应变。

1.2 曲率模态方程

为求解积分方程(15)，将曲率展开成幂级数形式

(16)

(17)

用ξk(k=0,1,…,N-1)乘以上式等号两边，并对等号两边在[0,1]范围内积分

(18)

以上方程最终可简化为如下形式

(19)

式中K和M分别为曲率刚度和质量矩阵，令m,n=1,2,…,N，其中j=n-1,k=m-1，则刚度和质量矩阵的任一元素Kmn和Mmn可表示为

(20)

2 模态设计

结构的应变模态振型决定了激励载荷作用下的动应变分布，尤其是当激励频率处于模态频率附近的时候。通过修改结构设计参数，以优化应变模态振型，从而避免动应变集中、提高结构寿命。

在初步设计阶段，许多工程结构可以简化成梁模型，例如传动轴、输运管道、高层建筑等。这些结构的质量函数与刚度函数在轴向都是变化的，有限个设计参数决定了结构的刚度与质量函数形式，分析模态振型对这些参数的灵敏度是进行模态设计的关键。

2.1 参数灵敏度

由于前面已经将连续梁的曲率模态求解转化为矩阵形式的广义特征值问题，所以就可以采用Nelson法求解特征向量(曲率模态展开系数)参数灵敏度。虽然Nelson在其文献[19]中的分析对象是标准特征值问题(质量矩阵为单位矩阵)，但是下面将证明其方法也同样适用于本文情况。

(21)

(22)

由式(20)可知，刚度矩阵K是对称的而质量矩阵M不一定对称。将右特征向量对质量矩阵进行归一化

ηiTMηj=δij,i,j=1,2,…,N

(23)

(24)

其中

今年以来中央出台了一系列减税降费措施，10月中旬，财政部预计全年可减轻税费负担1.3万亿元以上。去年是一万亿元，今年年初计划是1.1万亿元。

(26)

式中Vi是方程(24)的一个特解，μi为待定常数。令i=j，再对方程(23)求λ的导数并联合式(26)可得

(27)

(28)

(29)

(30)

(31)

(32)

至此就确定了原方程组(21)中奇异系数矩阵的N-1阶非奇异子矩阵。不妨令Vi的第k个元素为零，消去方程组(24)中的第k个方程，求解方程组

(33)

2.2 曲率梯度

为了能够明显地表现应变的变化，引入曲率梯度模态的概念，其是曲率模态振型在梁轴向方向的导数，反映了应变在轴向的梯度变化。

(34)

曲率梯度模态对于设计参数的变化也更加敏感，其参数灵敏度为

(35)

可以利用之前已经求得的曲率模态参数灵敏度结果来计算曲率梯度模态的灵敏度，缩短了设计周期。

3 算 例

3.1 等截面均质梁

表1 等截面均质梁前4阶等效模态频率

Tab.1 First four equivalent natural frequencies for a uniform cantilever beam

nExactN=3N=5N=7N=911．87511．87541．87511．87511．875124．69414．71524．69424．69414．694137．854810．86947．95237．85597．8548410．995511．336611．005410．9955

由表1可知，当幂级数项增加到N=9时便能精确求得前4阶模态频率(小数点后四位)，此时的模态振型如图1所示，积分法求得的模态振型与解析解高度一致。增大N值可提高精度。

图1 等截面均质梁前4阶曲率模态振型Fig.1 First four curvature mode shapes of a uniform cantilever beam

3.2 变截面梁

变截面梁构件在工程中有广泛应用，例如大跨度桥梁、传动轴等，研究其应变模态对结构动强度设计有重要意义。

假设一种变截面悬臂梁结构，材料为钢，弹性模量为E0、密度为ρ0，截面形状为圆形，半径变化规律为

(36)

式中R0为基础半径，R1(0≤R1≤1)为缩放率，反映了梁末端与初始端的截面积缩减度。当R0=0.05 m时，不同R1取值的梁截面半径如图2所示。

图2 变截面梁截面半径Fig.2 The radius of a non-uniform cross-section beam

表2 变截面梁前4阶模态频率(Hz)

Tab.2 First four natural frequencies for a non-uniform cross-section cantilever beam(Hz)

nIntegralMethodFineFEMError157．537357．53680．0009%2281．4566281．45400．0009%3737．0197737．01190．0011%41418．43321418．4815-0．0034%

图3 变截面梁前4阶曲率模态振型Fig.3 First four curvature mode shapes of a non-uniform cross-section beam

从表2和图3可以看出，积分法在计算梁刚度和质量函数为三角函数类时，同样拥有很高的精确度。

3.3 部段连接

在很多工程结构中，整体结构都是由多个部段通过连接结构组装在一起的，连接段导致了局部刚度和质量变化。在梁模型中，连接段的局部动力学特性可以表示为变化的刚度和质量函数，有限个设计参数决定了函数的具体形式。

(37)

式中α1和α2分别表示杨氏模量在ξ1和ξ2两处的变化率，H1(ξ)至H4(ξ)是Hermite插值多项式

(38)

图4 不同参数下的杨氏模量曲线Fig.4 The Young’s modulus curves with different parameters

图5 α1=-2时二阶曲率模态振型及灵敏度Fig.5 The second curvature mode shape and its sensitivity with respect to α1 when α1=-2

由于参数变化仅发生在中间段，对模态振型的影响也相应主要体现在中间段局部，故下文将以该局部为分析对象。

采用本文中的积分法，计算α1=-2时的第2阶曲率模态振型(中间段处于振型波谷位置，能够明显地反映参数变化)，如图5所示，同时计算该阶振型关于设计参数α1的灵敏度。图中曲率模态振型已经向起始点归一化，中间段是振型波谷且为负值，参数灵敏度在该区域为正值，表明在α1=-2处增大α1可以提高该部位的振型值、减小其绝对值，从而降低该区域的应变。

α1取不同值的二阶曲率模态振型如图6(a)所示，可见正如灵敏度分析预测：α1从-4增加到0的过程中，模态振型也相应地上移，中间段区域的应变值也随之减小。图6(b)中所示的曲率梯度模态振型更加直观地显示了参数变化对模态振型的影响，α1的增大使得曲率梯度的变化越来越平缓，应变集中程度越来越小。

由此可知，连接段刚度与部段刚度过渡越平滑越能够减小应变集中。这就要求在结构设计阶段，连接结构形式和连接材料(可参考功能梯度材料[21])的选择与设计应充分考虑刚度平滑过渡的要求，避免动应变集中、防止动载荷作用下的结构破坏。

图6 不同参数下的二阶曲率及曲率梯度模态振型Fig.6 The second curvature and curvature gradient mode shapes with different parameters

4 边界条件适用性

本章将分析上述积分方法对不同边界条件的适用性。

Yc(r)dr+Yr(0)ξ+Y(0)=0

(39)

边界条件就是对端点的位移Y(ξ)、转角θ(ξ)、弯矩M(ξ)以及剪力Q(ξ)的约束，后面3个物理量可分别表示为

(40)

将方程(39)等号两端在[0,ξ]范围内依次积分两次，并最终化简得

S(1)Yc(1)+Q(1)L3(1-ξ)

(41)

由于将曲率展开成式(16)的形式，所以有

(42)

ξ=0端为固支边界，ξ=1端的边界条件一般有3种形式：自由、简支和固支。ξ=1端为自由边界时，上文已经详细分析过，在此主要讨论另外两种边界形式。

4.1 固支-简支

在起始端固支、末端简支这种边界条件下，有

(43)

将上式分别代入到式(16)和(42)中，可得

(44)

则有

(45)

(46)

这时末端的剪力可以写成

(47)

将式(45)，(46)和(47)代入到方程(41)中可得

(48)

用ξk(k=0,1,…,N-3)乘以上式等号两边，并对等号两边在[0,1]范围内积分，最终可化简为式(19)的形式。令m,n=1,2,…,N-2，其中n=j-1,m=k+1，则刚度和质量矩阵的任一元素Kmn和Mmn可表示为

(49)

4.2 固支-固支

在两端都固支的这种边界条件下，有

(50)

此时有

(51)

(52)

这时末端的剪力可以写成

(53)

与上一节固支-简支边界条件的处理方式相同，令m,n=1,2,…,N-2，其中n=j-1,m=k+1，则该边界条件下，刚度和质量矩阵的任一元素Kmn，Mmn可表示为

(54)

4.3 算 例

以等截面均质梁为例，分别利用上两节所示的积分法求得(N=11)固支-简支、固支-固支两种边界条件的前3阶模态频率和曲率模态振型，与解析解[10]相比如表3所示。表中的“曲率模态振型平均误差(Curvature Modal Shapes Average Error)”指的是，在ξ∈[0,1]区间等距取100个数据点，所有数据点对应积分法与解析解振型数据误差的代数平均值。

表3 等截面均质梁两种边界条件前3阶曲率模态

Tab.3 First three curvature modal frequencies and shapes for a uniform beam with two boundary conditions

BoundaryconditionnFrequencieserrorCurvatureshapesaverageerrorClamped⁃simplysupported1-0．0001%-0．0077%2-0．0015%0．2558%30．0121%0．2999%Clamped⁃clamped10．0001%-0．0168%20．0192%-1．9071%30．0063%-0．5837%

从上表可以看出，积分法同样适用于这两种边界条件形式，且求解曲率模态频率和模态振型的精度都很高。

5 结 论

本文通过直接求解曲率模态及其参数灵敏度，将模态振型分析与动应变分布优化联系起来。提出曲率梯度模态的概念，更加直观、灵敏地展现动应变的空间变化及结构设计参数的影响。

曲率和曲率梯度的模态振型展开成幂级数，通过积分以曲率为变量的变系数运动微分方程，将连续体系统的模态求解问题转化为包含刚度和质量矩阵的广义特征值形式。幂级数系数为特征向量，矩阵分析证明了采用Nelson法求解该形式特征灵敏度的可行性，得到曲率和曲率梯度模态关于结构设计参数的精确灵敏度。3个典型算例证明，采用这种积分方法可准确地求解曲率模态、参数灵敏度分析能有效地优化动应变分布以避免动应变集中。

文中方法适用于任意刚度与质量函数的悬臂梁结构，同时也适用于固支-固支和固支-简支这两种边界条件。

[1] Li D B, Zhuge H C, Wang B. The principle and techniques of experimental strain modal analysis [A]. The 7th International Modal Analysis Conference [C]. Las Vegas, NV, USA, 1989: 1 285—1 289.

[2] Yam L H, Leung T P, Li D B, et al. Theoretical and experimental study of modal strain analysis [J]. Journal of Sound and Vibration, 1996, 192 (2): 251—260.

[3] Pandey A K, Biswas M, Samman M M. Damage detection from changes in curvature mode shapes [J]. Journal of Sound and Vibration, 1991, 145 (2): 321—332.

[4] 李德葆, 陆秋海, 秦权. 承弯结构的曲率模态分析[J]. 清华大学学报(自然科学版), 2002, 42 (2): 224—227.

Li Debao, Lu Qiuhai, Qin Quan. Curvature modal analysis of bending structures [J]. J. Tsinghua Univ. (Sci. & Tech.), 2002, 42 (2): 224—227.

[5] Sazonov E, Klinkhachorn P. Optimal spatial sampling interval for damage detection by curvature or strain energy mode shapes [J]. Journal of Sound and Vibration, 2005, 285 (4): 783—801.

[6] Qiao P Z, Lu K, Lestari W, et al. Curvature mode shape-based detection in composite laminated plates [J]. Composite Structures, 2007, 80 (3): 409—428.

[7] Abdo M A B. Damage detection in plate-like structures using high-order mode shape derivatives [J]. International Journal of Civil and Structure Engineering, 2012, 2 (3): 801—816.

[8] Whalen T M. The behavior of higher mode shape derivatives in damaged beam-like structures [J]. Journal of Sound and Vibration, 2008, 309 (3): 426—464.

[9] Thomson W T, Dehleh M D. Theory of Vibration with Applications [M]. New Jersey: Prentice-Hall, 2005.

[10]Inman D J. Engineering Vibration [M]. New Jersey: Prentice-Hall, 2007.

[11]Naguleswaran S. A direct solution for the transverse vibration of Euler-Bernoulli wedge and cone beams [J]. Journal of Sound and Vibration, 1994, 172 (3): 289—304.

[12]Melnikov Y A. Influence Functions and Matrices [M]. New York: Marcel Dekker, 1999.

[13]Li Q S. A new exact approach for determining natural frequencies and mode shapes of non-uniform shear beams with arbitrary distribution of mass or stiffness [J]. International Journal of Solids and Structures, 2000, 37 (37): 5 123—5 141.

[14]Elishakoff I, Candan S. Apparently first closed-form solutions for vibrating inhomogeneous beams [J]. International Journal of Solids and Structures, 2001, 38 (19): 3 411—3 441.

[15]Huang Y, Li X F. A new approach for free vibration of axially functionally graded beams with non-uniform cross-section [J]. Journal of Sound and Vibration, 2010, 329 (11): 2 291—2 303.

[16]Sulter T R, Gamarda C J, Walsh J L, et al. Comparison of several methods for calculating vibration mode shape derivatives [J]. AIAA Journal, 1988, 26 (12): 1 506—1 511.

[17]Fox R L, Kapoor M P. Rates of change of eigenvalues and eigenvectors [J]. AIAA Journal, 1968, 6 (12): 2 426—2 429.

[18]Wang B P. An improved approximate method for computing eigenvector derivatives [J]. Finite Element in Analysis and Design, 1993, 14 (4): 381—392.

[19]Nelson R B. Simplified calculation of eigenvector derivatives [J]. AIAA Journal, 1976, 14 (9): 1 201—1 205.

[20]Yu M, Liu Z S, Wang D J. Comparison of several approximate modal methods for computing mode shape derivatives [J]. Computers & Structures, 1996, 62 (2):301—393.

[21]Birman V, Byrd L W. Modeling and analysis of functionally graded materials and structures [J]. Applied Mechanics Reviews, 2007, 60(5): 195—216.

The integral method and sensitivity analysis for strain modes of cantilever beam-like structures

XINGJian-wei,ZHENGGang-tie

(School of Aerospace, Tsinghua University, Beijing 100084, China)

The spatial distribution of dynamic strain is determined by the modal shape and hence the strain at a certain area can be reduced by modal shape design. Based on the proportional relation between strain and curvature, this paper proposes an integral method to solve curvature modal frequencies and shapes directly of cantilever beam-like structures with variable stiffness and mass functions. The method is validated by comparing the computation results with both numerical and analytical solutions. Furthermore, the presented method transforms the problem solving of continuous systems modes into a generalized eigenvalue form of finite degrees of freedom. The explicit expressions of curvature eigen-sensitivity have been established with Nelson's method, and the curvature gradient modes are also calculated to show the change of dynamic strain more visibly. An example of parameter design is presented with the proposed sensitivity analysis method. Moreover, this method is proved capable of treating the other two boundary conditions of clamped-simply supported and clamped-clamped.

strain mode； curvature mode； integral method； sensitivity analysis; curvature gradient mode

2014-04-09；

2014-09-25

V414.1

A

1004-4523(2015)06-0855-10

10.16385/j.cnki.issn.1004-4523.2015.06.001

邢建伟(1986—)，男，博士。电话：13811950902；E-mail: xingjianwei1019@126.com

郑钢铁(1960—)，男，教授。电话：(010)62783235；E-mail: gtzheng@tsinghua.edu.cn