基于本征方程求解復合材料梁幾何非線性靜力平衡

李園園， 何　歡， 陳國平， 劉龐輪(南京航空航天大學 機械結構力學及控制國家重點實驗室，南京　210016)

?

基于本征方程求解復合材料梁幾何非線性靜力平衡

李園園， 何歡， 陳國平， 劉龐輪(南京航空航天大學 機械結構力學及控制國家重點實驗室，南京210016)

摘要：對基于本征梁理論求解復合材料梁的幾何非線性大變形屈曲問題進行了研究，根據材料屬性利用漸近變分法確定復合材料梁的剛度矩陣，再根據本構方程和平衡方程求得其靜力學行為，結果表明：對單層鋪層的復合材料梁來說，剛度矩陣的耦合項可以忽略，其變形構型及梁末端的位移及轉角的變化趨勢與各項同性材料相同；對一個一般的復合材料梁來說，其剛度矩陣的耦合項不可忽略，耦合項對位移和轉角的影響與施加在梁上的載荷大小有關，在載荷小于30 N，以耦合項50%的變化量為界，當變化量小于50%時，位移和轉角的變化趨勢與初始時相同，當變化量大于50%時，位移和轉角的變化趨勢發生很大的改變，但與解耦后的變化趨勢相似。

關鍵詞：本征梁理論；復合材料梁；漸近變分法；幾何非線性；大變形；屈曲

在飛行器、運載火箭、人造衛星等的航空航天領域，以彈性體為主裝配部件的大范圍運動裝置比比皆是。隨著國民經濟的高速發展，各種機械系統在向著構型復雜化，結構輕巧化，運轉高速化方向發展，所有的這些要求導致系統部件要有較大的柔性，因此，空間大范圍運動柔性結構的研究成為一個引人注目的話題，其基礎結構包括梁，板，殼等。

梁作為一種常見的單元應用于不同的機械和結構單元中，它在航空工業中的典型應用是大展弦比機翼和直升機旋翼，在分析梁的大范圍運動時，由于變形過程的復雜性，單純的線性理論已不足以描述這種復雜的運動，而非線性變形的量化開始倍受關注，很多學者對其進行了研究。Christensen等[1-2]從應變與位移的非線性角度出發，Boutaghou等[3]則采用Hamilton原理，各自建立了柔性梁的非線性力學模型，Hodges[4]基于幾何精確的本征梁理論，它采用線應變和角應變代替位移和轉角作為獨立變量。Gatti-Bono等[5]在不考慮剪切變形的情況下得到了和文獻[4]類似的方程。馮志華[6]基于Kane方程的建模方法研究了大范圍運動柔性梁的非線性行為。齊朝暉[7]運用虛功原理推倒了大變形梁的動力學公式，并給出線彈性材料的本構關系。利用上述學者建立的模型求解梁的動力學或靜力學行為時，基本都是基于能量法和有限元法，這些方法一般都需要多次迭代，且隨著幾何的變形，有限元方法必須不斷修正載荷和剛度矩陣，增加了計算的費用和仿真時間。

由于梁可以被看成是維數降階的結構，以其軸線方向x1作為參考，其他所有量都可以表達為x1的函數，因此梁的非線性方程中只有一個獨立變量，而迭代打靶法就是針對只有一個獨立變量的邊界值問題的數值求解方法[8]，后來，Shvartsman[9-11]提出了一種求解受伴隨載荷的平面懸臂梁問題的直接打靶方法，這些方法求解的非線性方程組階數較低，但它仍需要多次迭代，計算量較大，且直接打靶法局限于利用Euler-Bernoulli梁理論，不能推廣到三維的情況。

為克服上述的計算困難，Karlson等[12]提出了基于本征梁方程的打靶法，用于計算三維空間內，具有預彎曲型梁的大變形平衡問題。對于受非保守載荷的懸臂梁邊界值問題，該方法不需要迭代就能得到一個精確解。同時，這些方程是以1階形式定義的，因此僅需要一個初始條件就可以計算。

本文基于文獻[12]的建模和求解方法，結合Hodges提出的漸近變分法[13-16](Variational-Asymptotic Beam Sectional Analysis，VABS)求解了復合材料梁的幾何非線性大變形問題，比較了復合材料梁與各項同性材料梁的變形構型及末端的位移和轉角的變化趨勢，同時探討了某種復合材料梁的耦合項對變形構型及末端位移和轉角的影響，為后續進一步研究提供新的思路和理論依據。

1本征梁方程

1.1運動學方程

如圖1描述了一個長度為l的復合材料梁的變形過程，可用三個構型來描述：參考構型Ωref，初始構型Ω0和變形構型Ωf，其中，參考構型表示沒有預彎曲的直梁，而初始構型存在預彎曲γ0和κ0，變形構型是用內力表示的應力狀態下的構型，此時的線應變和角應變相對于參考構型分別表示為γf和κf。在上述三種構型中分別定義三組基：[I1,I2,I3]T，[b1,b2,b3]T和[B1,B2,B3]T，這三組基在各自的構型中是互相正交的，且[B1,B2,B3]T隨著參考軸線坐標的變化而變化，b1=b2×b3，B1=B2×B3，其中，b1和梁的參考軸線x1相切，而由于剪切變形的存在，B1不一定相切于參考軸線s，參考軸線由x1變為s的位移向量表示為u(x1,t)。變形構型的基與參考構型的基(笛卡爾坐標系)之間有關系：Βi=BixI1+BiyI2+BizI3，其中，i=1,2,3，Bix,Biy,Biz表示Bi在笛卡爾坐標系下的坐標分量，初始構型的基與參考構型的基之間也存在類似的關系。由圖1可以看出，梁從構型Ω0變化到構型Ωf后橫截面發生翹曲，其中的虛線表示初始構型的橫截面形狀，陰影表示翹曲了的橫截面形狀，橫截面翹曲后存在翹曲位移。

根據文獻[14]，在初始構型Ω0中，沿參考軸線x1的任意一點處，梁中心線的位置向量r與橫截面上任意一點的位置向量r*以及基坐標之間存在如下的關系：

r*(x1,x2,x3)=r(x1)+xαbα

(1)

式中,α=2,3，x2,x3是在橫截面上與x1正交的單位矢量且初始時分別相切于b2,b3。

在Ωf處，梁中心線的位置向量R與橫截面上任意一點的位置向量R*以及基坐標和翹曲位移之間存在如下的關系：

R*(x1,x2,x3)=R(x1)+xαBα(x1)+wiBi(x1)=

r(x1)+u(x1,t)+xαBα(x1)+wiBi(x1)

(2)

式中，wi是與x1,x2,x3有關的翹曲矢量。

根據文獻[12]，在變形構型Ωf中，中心線位置向量R和基矢量以及線應變和彎曲變形之間存在如下的關系式：

R′j=(1+γ11)B1j+2γ12B2j+2γ13B3j

(3)

B′ij=κ×Bij

(4)

(5)

圖1　梁的變形簡圖Fig.1 Schematic of beam deformation

1.2平衡方程

文獻[17]提出的矩陣形式的運動方程定義了梁的線速度、角速度、彎曲應變及線應變的空間和時間變化，如式(6)和(7)：

(6)

M′+κf×M+(e1+γf)×F+m=

(7)

式中，上標“·”表示關于時間t的偏導數，內力F=[F1,F2,F3]T；彎矩M=[M1,M2,M3]T；單位長度上的外力和外力偶f=[f1,f2,f3]T和m=[m1,m2,m3]T；對應于線速度V=[V1,V2,V3]T和角速度Ω=[Ω1,Ω2,Ω3]T的單位長度線動量和角動量P=[p1,p2,p3]T和H=[H1,H2,H3]T；e1表示B1方向的單位矢量[1,0,0]T，γf=[γ11,2γ12,2γ13]T。

略去運動方程式(6)和(7)中的時間項與時間導數項，得到平衡方程式(8)和(9)：

F′+κf×F+f=0

(8)

M′+κf×M+(e1+γf)×F+m=0

(9)

內力和彎矩與彎曲變形和應變的本構方程式如下[18]：

(10)

其中，

(11)

為剛度矩陣，Sij(i=1,2,3,4,5,6，j=1,2,3,4,5,6)表示橫截面的剛度系數，且SijT=Sji(i≠j)，下標1代表拉伸，2、3代表剪切， 4代表扭轉，5、6代表彎曲。

對應變很小的各項同性材料梁，本構方程式(10)是線性的，此時的剛度矩陣D只存在對角線項，耦合項忽略不計，該剛度矩陣利用材料屬性通過維數降階即可獲得[12]。然而，對初始彎曲變形不足夠小的各項同性材料或對各項異性材料來說，剛度矩陣中的某些非對角線元素不能省略，這些非對角元素為耦合項，此時的剛度矩陣已不能通過簡單的維數降階獲得。

注意到，在忽略剪切變形的情況下，式(10)得到彎曲變形和彎矩的關系：

(12)

其中，

(13)

2剛度矩陣的求解-漸近變分法(VABS)

前面提到，對初始彎曲變形不足夠小的各項同性材料或對各項異性材料來說，剛度矩陣D不能通過簡單的維數降階獲得，Hodges等開始將漸近變分法應用于求解復合材料梁的橫截面彈性常數和翹曲場。對于梁這種維數降階的結構，VABS是強有力的工具，它可以將一個一般的三維非線性彈性問題的求解轉化為二維的線性橫截面分析和一維的非線性分析兩部分。經過不斷發展，該方法已形成單獨的工程軟件包來進行具有任意橫截面形狀的各向異性材料梁的2D橫截面分析。

在用漸近變分法求解橫截面彈性常數時，所依據的基本理論是能量法，使容易求解的一維應變能近似代替原三維梁的應變能。式(14)給出了Timoshenko梁(1D)單位長度的能量與彈性常數和應變之間的關系[13]：

2U=εTXε+2εTYγs+γsTGγs

(14)

式中，ε=[γ11,κ1,κ2,κ3]T，分別代表廣義應變(1D)中的拉伸，扭轉和兩個彎曲，γs=[2γ12,2γ13]T，X,Y,G即為組成6×6的剛度矩陣D的子矩陣。

式(14)可以轉化為本構方程的如下形式：

{FM}=[Ds]{s}

(15)

式中，{FM}={F1M1M2M3F2F3}T，[Ds]=[XY;YTG]，{s}={εγs}T。

為求得X,Y,G，需將基于二階漸進理論的梁橫截面上的能量轉化為1D的Timoshenko梁形式，轉化過程涉及1D的動力學方程、本構方程和靜力學平衡方程。基于二階漸進理論的梁橫截面上的應變能如下[15]：

(16)

將式(8)和(9)中所有的外力置為零，整理得到梁的靜力學平衡方程變為如下形式：

{FM}′+[DM]{FM}=0

(18)

(19)

式中,R,S,T為橫截面柔度矩陣的子矩陣。由式(19)遞推得到：

(20)

將式(19)和(20)代入式(16)中，并使整理后的方程與式(14)對應，根據等號兩邊對應項相等即可求得X,Y,G的表達式，即確定了剛度矩陣。

3打靶法

對于受伴隨力作用的懸臂梁來說，“打靶”從梁的自由端開始，一次“打靶”就能確定梁的非線性變形，不需要迭代，打靶過程自動實施了梁末端的固定邊界。基于本征方程的打靶法就是將式(3)，(4)，(5)，(8)和(9)利用式(10)和(11)展開為具有初值問題的如下形式：

y′(x1)=f(x1,y(x1))

(21)

[uvw]T=R

(23)

式中，u，v，w分別對應全局坐標系下的X，Y，Z方向的位移。

基于上述理論，圖2給出了本文基于本征方程求解復合材料梁的靜力學問題的流程圖。

圖2　復合材料梁靜力學分析的流程圖Fig.2 Flow chart of the static anlysis for composite beams

4算例

4.1算例1

算例1在忽略剪切變形，沒有預彎曲的情況下，給出了末端分別承受不同大小單點力和分布力的各向異性材料直梁的變形構型、梁末端位移及轉角的變化，并與文獻[12]中承受單點力和分布力的各項同性材料直梁的變形構型、梁末端的位移及轉角變化進行了對比。

表1給出了鋪層只有一層的正交各向異性材料梁的材料屬性、尺寸以及VABS輸出的彈性常數，需要指出的是，表1給出的數值是在文獻[20]的基礎上進行單位換算得到的。其中，梁的長度是任意的，只要比橫截面尺寸大的多即可。

其中，寬度方向沿橫截面的x2，厚度方向沿橫截面的x3，梁的長度L在本算中例取0.53 m，μ為梁單位長度的質量，i2,i3為橫截面質量慣性矩，由表1發現，單層鋪層的各向異性材料的剛度矩陣為對角陣，耦合項很小忽略不計。

表1　材料屬性、橫截面維數以及彈性常數

圖3和圖5分別為末端受不同大小單點力和分布力的各向異性材料直梁的變形構型，圖4和圖6分別為對應的梁末端位移及梁末端關于Y軸轉角的變化趨勢。

對比圖3～圖6以及文獻[12]發現：對鋪層只有一層的正交各向異性材料梁，其剛度矩陣耦合項很小以致忽略不計，只存在對角線元素，其變形構型、梁末端位移及轉角的變化趨勢與各項同性材料梁相同。

圖3　受垂直梁末端單點伴隨力作用的變形構型Fig.3 Deformation configuration of the beam loaded with a perpendicular, point, follower force

圖4　單點力作用下的位移及轉角Fig.4 Normalized displacements and rotations under perpendicular, point, follower force

圖5　受垂直梁均布伴隨載荷作用的變形構型Fig.5 Deformation configuration of a straight beam loaded with a perpendicular, distributed, follower force

圖6　均布載荷作用下的位移及轉角Fig.6 Normalized displacements and rotations under perpendicular, distributed, follower force

4.2算例2

該算例同樣在沒有剪切變形及預彎曲的情況下，針對剛度矩陣耦合項不可忽略的各項異性材料梁，當末端承受不同大小的單點力時，觀察梁末端位移及轉角受耦合項的影響。表2給出了該各向異性材料梁的材料屬性，橫截面維數以及彈性常數。需要指出的是，表2給出的數值是在文獻[21]的基礎上進行單位換算得到的。

表2　材料屬性、橫截面維數以及彈性常數

為便于計算，該梁的長度L仍取0.53 m,由表2可以看出，該材料的拉-剪、彎-扭耦合S12,S45不可忽略。圖7～圖9分別給出了復合材料梁在單點力作用下，彎-扭耦合對梁末端位移及梁末端關于Y軸轉角的影響。

圖7　單點力作用下彎扭耦合對Y方向位移的影響Fig.7 The normalized Y displacements impacted by bending and torsion coupling under perpendicular, point, follower force

圖8單點力作用下彎扭耦合對X方向位移的影響Fig.8 The normalized X displacements impacted by bending and torsion coupling under perpendicular, point, follower force

圖9　單點力作用下彎扭耦合對轉角的影響Fig.9 The normalized rotations impacted by bending and torsion coupling under perpendicular, point, follower force

圖7～圖9是彎扭耦合由0.008 kg·m2按圖示的箭頭分別減小10%，20%，35%，40%，50%，60%，70%，75%，80%，85%，90%，93%，95%和98%得到的彎扭耦合對位移和轉角的影響曲線，其中，細實線表示位移和轉角隨彎扭耦合的變化規律，粗實線表示解耦后位移和轉角的變化趨勢。由這三幅圖可以看出，當載荷約小于10 N時，彎扭耦合對位移和轉角的影響較小，當載荷約大于30 N時，彎扭耦合對位移和轉角的影響較復雜，當載荷在10～30 N之間時，彎扭耦合的影響開始呈現一定的規律：當彎扭耦合的變化量在0～50%之間時，隨著彎扭耦合的變化，梁末端位移及轉角的變化趨勢與彎扭耦合為0.008 kg·m2(初始時)時的變化趨勢相同，當彎扭耦合的變化量大于50%時，位移和轉角按另一種變化趨勢有規律的變化，但與解耦后相近，當完全解耦后，梁末端的位移及轉角的變化趨勢與各項同性材料一致，且載荷在10～20 N之間變化時，位移和轉角隨耦合項成比例的變化。

5結論

(1) 對單層鋪層的各向異性材料梁來說，由于其剛度矩陣的耦合項忽略不計，此時的本構方程形式類似于各項同性材料，最終求得的變形構型、梁末端的位移及轉角的變化趨勢與各項同性材料的一致。

(2) 對剛度矩陣的耦合項不可忽略的各向異性材料來說，剛度矩陣沒有解耦之前，耦合項對梁末端位移和轉角的影響與載荷有關，載荷較小時耦合項的影響也小，幾乎成比例變化，載荷較大時影響復雜，當載荷在10～30 N之間變化時，彎扭耦合以50%的變化量為分界點，分別呈現不同的變化規律，當完全解耦后，位移和轉角的變化趨勢與各項同性材料一致。

參 考 文 獻

[ 1 ] Christensen E R, Lee S W. Nonlinear finite element modeling of the dynamics of unrestrained flexible structures[J]. Computers and Structures, 1986, 23(4): 819-829.

[ 2 ] Belytschko T, Hsieh B. Nonlinear transient finite element analysis with convected coordinates[J]. International Journal for Numerical Methods in Engineering,1973,26(2):255-271.

[ 3 ] Boutaghou Z E, Erdman A G, Stolarski H K. Dynamics of flexible beams and plates in large overall motion[J]. ASME Journal of Applied Mechanics, 1992, 59(4):991-999.

[ 4 ] Hodges D. Geometrically exact, intrinsic theory for dynamics of curved and twisted anisotropic beams[J]. AIAA Journal, 2003, 41 (6): 1131-1137.

[ 5 ] Gatti-Bono C, Perkins N. Physical and numerical modeling of the dynamic behavior of a fly line[J]. Journal of Sound and Vibration, 2002, 255 (3): 555-577.

[ 6 ] 馮志華. 大范圍運動柔性梁非線性動力學[D]. 南京：南京航空航天大學, 2002.

[ 7 ] 齊朝暉. 大變形梁的動力學公式及本構關系[J]. 吉林工業大學學報, 1994, 24(73): 24-29.

QI Zhao-hui.Kinetic equations and constitutive relations of large deformation beam[J]. Journal of Jilin University of Technology, 1994, 24(73): 24-29.

[ 8 ] Roberts S, Shipman J. Two-point boundary value problems: shooting methods[M].New York:Elsevier,1972.

[ 9 ] Shvartsman B. Direct method for analysis of flexible beam under a follower load[J]. Proceedings of Computational Mechanics for the Next Millennium, 1999:155-160.

[10] Shvartsman B. Large deflections of a cantilever beam subjected to a follower force[J]. Journal of Sound and Vibration, 2007, 304 (3): 969-973.

[11] Shvartsman B. Direct method for analysis of flexible cantilever beam subjected to two follower forces[J]. International Journal of Non-Linear Mechanics,2009,44(2):249-252.

[12] Karlson K N, Leamy M J. Three-dimensional equilibria of nonlinear pre-curved beams using an intrinsic formulation and shooting[J]. International Journal of Solids and Structures, 2013, 50: 3491-3504.

[13] Hodges D H. Nonlinear composite beam theory[M]. Reston, VA:AIAA, 2006.

[14] Yu W, Hodges D H, Volovoi V V, et al. On timoshenko-like modeling of initially curved and twisted composite beams[J]. International Journal of Solids and Structures,2002,39(19):5101-5121.

[15] Yu W, Hodges D H. Generalized timoshenko theory of the variational asymptotic beam sectional analysis[J]. Journal of the American Helicopter Society, 2005, 50 (1): 46-55.

[16] Yu Wen-bin, Hodges D H, Volovoi V V, et al. On Timoshenko-like modeling of initially curved and twisted composite beams[J]. International Journal of Solids and Structures, 2002, 39: 5101-5121.

[17] Hodges D H. A mixed variational formulation based on exact intrinsic equations for dynamics of moving beams[J]. International Journal of Solids and Structures,1990,26(11):1253-1273.

[18] Leamy M. Intrinsic finite element modeling of nonlinear dynamic response in helical springs[J]. Journal of Computational and Nonlinear Dynamics,2012,7(3):031007.

[19] Cesnik C E S, Hodges D H. VABS: a new concept for composite rotor blade cross-sectional modeling[J]. Journal of the American Helicopter Society, 1997, 42(1): 27-38.

[20] Kovvali R, Hodges D H. Verification of variational asymptotic sectional analysis for initially curved and twisted beams [J]. Journal of Aircraft, 2012, 49(3): 861-869.

[21] Popescu B, Hodges D H. On asymptotically correct Timoshenko-like anisotropic beam theory [J]. International Journal of Solids and Structures, 2000, 37: 535-558.

Geometric nonlinear static equilibria of composite beams using intrinsic formulation

LIYuan-yuan,HEHuan,CHENGuo-ping,LIUPang-lun(The State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China)

Abstract:Based on intrinsic beam theory, this paper solved the large deformation buckling problem of geometric nonlinear composite beams.Using the asymptotic variational method, we can get the stiffness matrix of the composite beam in light of material properties.Then, the static behavior of composite beams could be obtained through the balance equation and constitutive equation.The results show that if the composite beam has a single layer, the coupling terms of the stiffness matrix can be ignored and the trends of the deformation configuration, normalized displacements and rotations of the beam end are the same as the isotropic material beams.However, for a general composite beam, the coupling terms of its stiffness matrix cannot be ignored, and the impact of coupling term on displacement and rotation change regularly according to load.When the load is less than 10lb.if the amount of the coupling term change is less than 50%, the displacement and rotation are the same as the initial trends; however, if the change is greater than 50%, they will undergo great changes.

Key words:intrinsic beam theory; composite beam; variational asymptotic method; geometric nonlinearity; large deformation; buckling

中圖分類號:V214.3

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.08.010

通信作者陳國平 男，博士，教授，1956年7月生

收稿日期：2014-12-03修改稿收到日期：2015-04-20

基金項目：國家自然科學基金資助項目(11472132)；中央高校基本科研業務費資助項目(NS2014002)；江蘇高校優勢學科建設工程資助項目

第一作者 李園園 女，博士生，1986年3月生