懸臂梁結構的應變模態積分法與靈敏度分析

邢建偉， 鄭鋼鐵

(清華大學航天航空學院， 北京 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