劉昊 瞿葉高 孟光

(上海交通大學機械系統(tǒng)與振動國家重點實驗室,上海 200240)

引言

軸向流動下彈性梁或板結構的穩(wěn)定性問題是流固耦合動力學中經典的問題之一,廣泛存在于航空工程、海洋工程、核能工程[1-3]等領域.置于流場中前緣固支的柔性梁或板狀結構在來流速度超過臨界流速時,會產生周期性的大撓度自激振動現(xiàn)象,這種振動現(xiàn)象可以歸因于運動誘發(fā)激勵機制(movementinduced excitation,MIE)[4].

國內外針對均勻軸向流中柔性梁或板結構振動問題開展了一些研究.早期Taneda[5]基于實驗方法研究了旗幟的各種振動模式.Eloy 等[6]采用穩(wěn)定性分析方法分析了軸向流中懸臂柔性板的周期性極限環(huán)振動和混沌振動機制.Zhu 和Peskin[7]采用浸入邊界法研究了細絲在流動的肥皂膜中的振動現(xiàn)象,揭示了更長的細絲向不穩(wěn)定振動模式的過渡和雙穩(wěn)態(tài)區(qū)域的存在.文獻[8]采用基于有限差分格式的直接模擬(fluid-structure direct simulation,FSDS)方法,以質量比為變量,分析了均勻流作用下柔性板的振動特性.Lui 等[9]提出CEFI (combined field with explicit interface)公式,采用有限元方法分析質了量比和雷諾數的改變對柔性板振動特性的影響.Zhang 等[10]采用一種松耦合格式研究了強附加質量效應下二維柔性板的流致振動現(xiàn)象.以往的研究[11-12]揭示了均勻流作用下柔性板具有三種不同的變形模式:定點穩(wěn)定模式、周期性極限環(huán)振動模式以及混沌振動模式.最近,Saravanakumar 等[13]開展了復合材料層合梁在均勻流中的振動特性研究,分析了鋪層剛度和密度對層合梁振幅、頻率和渦脫落模式的影響.

實際工程中存在很多剪切流流場,而剪切流動作用下彈性結構的流固耦合動力學行為更加復雜.Cheng 等[14]采用LBM (lattice boltzmann method)方法分析了作用在圓柱和方柱上的二維不可壓縮線性剪切流.結果表明,圓柱后的渦結構與剪切速率有很大關系,作用在圓柱上的升力和阻力一般隨剪切速率的增大而減小,通過方形圓柱的渦流脫落頻率隨剪切速率的增大而減小.Yu 等[15]采用高精度譜方法研究了剪切流對NACA0012 翼面上渦結構的演變及相應的氣動性能的影響.Liu 等[16]研究了水翼在不同剪切速率的剪切流作用下振動的能量收集效率問題.文獻[17-18]對二維線性剪切來流作用下彈性支撐圓柱體結構物雙自由度流致振動問題進行了數值計算,發(fā)現(xiàn)圓柱在剪切流中的運動軌跡呈液滴形狀,這與均勻流的“8”形軌跡明顯不同.

復合材料彈性結構已廣泛應用于船舶艦艇、海洋工程等領域結構設計[19-20].研究剪切流動作用下復合材料層合梁振動機制對于上述領域結構設計具有重要的意義.關于單層梁在均勻流作用下的振動已經開展了大量研究,但對于復合材料層合梁在剪切流作用下的振動研究較少,剪切流分布、層合材料的差異對于振動特性的影響和機理尚不清楚.本文針對水下剪切流動下復合材料層合梁振動問題,基于復合材料層合梁高階剪切鋸齒理論[21-23],建立了不可壓縮黏性剪切流中層合梁高精度的非線性流固耦合動力學數值模型,研究了不同剪切流動下不同剛度比、不同鋪層厚度、不同鋪層角度層合梁的大變形非線性振動行為與機理.本文為水下剪切流動作用下復合材料梁的動力學響應預測和復合材料結構設計提供參考依據.

1 模型描述

1.1 計算模型

本文考慮一個沉浸在水中的復合材料層合梁,層合梁的長寬尺寸為L×H,H=0.01L,其幾何中心距流體域左邊界距離為2L,到右邊界距離為10L,到流體域上下邊界的距離均為2L,如圖1 所示.不考慮流體和結構所受的重力.將水介質假設為不可壓縮黏性流體,流體密度為 ρf,動力黏度為 μf.左端入口處的剪切流速度分布函數為U(Y)=U0(Y/2L)α,層合梁中軸線處的流速為U0,對應的雷諾數為Re=ρf U0L/μf.

圖1 計算模型Fig.1 Computational model

1.2 邊界條件

浸沒在流體中的層合梁前端(x=0)固支,右端(x=L)自由.流場區(qū)域的左端入口采取流速的邊界條件:u=U(Y),v=0 ;右端出口設置壓力出口邊界條件:p=0 ;層合梁結構表面為流固耦合界面,設置無滑移邊界條件;流域上下邊界設置為滑移固壁邊界條件:?u/?y=0,v=0 .

2 流固耦合建模與計算方法

2.1 流體數值計算模型

不可壓縮黏性流體的控制方程為Navier-Stokes方程,本文采用ALE (arbitrary Lagrangian Eulerian)方法描述流體網格的變形.在ALE 框架下,流體控制方程為[24]

其中,uf是流體速度矢量，ρf表示流體密度,c是ALE 速度矢量c=uf?um,um表示流體網格變形速度.將水介質假設為不可壓縮牛頓流體,應力張量σf定義為

采用有限體積法對流體控制方程進行離散[25],任意一個控制體單元內的離散方程可以表示為

其中,V代表控制體體積,S代表控制體表面積,nj是控制體表面的單位法向向量,um,j表示控制體j方向的網格變形速度.時間離散采用二階隱式歐拉方法,對流項離散采用線性迎風差分方法,擴散項離散采用中心差分方法,速度和壓力耦合項采用PIMPLE(PISO-SIMPLE)算法[26]求解.

2.2 層合梁結構動力學數值模型

基于考慮層合梁面內位移鋸齒效應的高階剪切梁鋸齒理論來描述層合梁的振動變形.層合梁內任意一點P(x,z)在x和z方向的位移可以表示為[27]

其中,u,w,? 和 η 表示層合梁中性軸上的廣義位移,僅與坐標x和時間t有關.f(z)和g(z) 是廣義位移分布形函數,反映了位移沿厚度方向的分布情況,其表達式為[28]

式(5)中 φ (z,k) 是鋸齒函數(zig-zag 函數),與坐標z及鋪層數目k有關,定義為[29]

基于von Kármán 位移?應變關系來考慮層合梁的幾何非線性大變形效應.層合梁的應變分量表達式為

忽略層合梁厚度方向的正應力和正應變,正交各向異性材料鋪層的應力應變關系如下[30]

采用有限元方法對層合梁的運動方程進行離散,采用2 節(jié)點10 自由度梁單元對層合梁進行離散,對軸向位移分量u和廣義位移變量?,η 采用線性插值,對橫向位移分量w采用Hermite 函數進行插值.層合梁非線性振動有限元方程可寫為[30]

式中,M表示層合梁的質量矩陣,KL和KNL分別表示層合梁的線性和非線性剛度矩陣,C是阻尼矩陣,動力學分析中一般采用瑞利阻尼描述.Ff為外部流場載荷向量.

采用直接積分Newmark-β方法結合Newton-Raphson 迭代方法對層合梁結構的非線性振動有限元方程進行求解.在每個時間步長內,需要根據層合梁與流體的耦合作用來迭代計算Ff.

2.3 雙向流固耦合算法

考慮流體與層合梁的強耦合作用.在流體與結構耦合界面上,需要滿足力平衡條件和速度協(xié)調條件

式中,n為流固耦合界面的單位法向量,為流固耦合界面處結構的速度.

本文采用分區(qū)流固耦合計算方法來迭代求解流動響應和結構振動響應,如圖3 所示.由于流體和結構網格是非匹配的,需要對流體域計算后獲得的流固耦合界面力ff si=FΓ進行插值計算,將其作為層合梁的外載荷Ff.對于不可壓縮牛頓流體,邊界上的力Ff可由Cauchy 應力張量與表面法向向量的內積,再乘以表面面積得到[26]

圖2 復合材料層合梁的示意圖Fig.2 Schematic diagram of composite laminated beam

圖3 雙向流固耦合計算流程圖Fig.3 Flowchart of bidirectional fluid-structure interaction

結構在流場載荷作用下產生振動變形,更新結構位移xf si=us,將其插值后傳遞給流場作為下一個耦合迭代步流場的運動邊界條件.在一個時間步內進行多次耦合迭代,直到流場、結構以及流固耦合界面都達到設定的力收斂和位移收斂準則,程序再推進下一時間步計算.本文采用徑向基函數(radial basis functions,RBF)[31]方法對流固耦合界面上流場和結構非匹配節(jié)點的物理量數據進行插值.

3 數值計算和結果分析

3.1 算例驗證

采用不同數目網格來分析層合梁(?=0)在均勻流作用下的振動響應收斂特性.采用3 節(jié)點三角形非結構單元對流體區(qū)域進行離散,采用2 節(jié)點梁單元對層合梁進行空間離散.考慮了三種計算網格:Mesh A (NE=11860,NP=6109,NM=50)、Mesh B(NE=19426,NP=10220,NM=60)和Mesh C (NE=41624,NP=21184,NM=80),其中NE表示流場單元數,NP表示流場節(jié)點數,NM表示梁單元數目.圖4給出了網格C 對應的流場非結構網格,層合梁的表面布置三層邊界層網格,其中第一層網格高度根據無量綱壁面距離y+=1 計算得到.三種網格得到的層合梁無量綱振幅uy/L和無量綱振動頻率f L/U0如表1所示,表中還給出了文獻[13]數值計算得到的結果.結果表明,由Mesh C 網格得到的層合梁振動幅值與文獻[13]結果之間的最大偏差為2.4%,后續(xù)的數值計算采用網格C.

圖4 流場非結構網格C Fig.4 Overview of the Mesh C

表1 網格收斂性分析Table 1 Grid independence test

考慮層合梁不同鋪層彎曲剛度之差 ? 情況下,分析層合梁上下兩層彈性模量的差異對層合梁振動響應的影響.圖5 和圖6 給出了10 個無量綱時長內,層合梁右端(x=L)節(jié)點振動最大幅值和主導頻率隨鋪層彎曲剛度之差的變化.圖中還將本文數值計算結果與文獻[13]的結果進行對比.結果表明,層合材料的差異對梁的振動響應具有顯著影響,隨著 ? 的增加,層合梁振動幅值增大,主導頻率有下降的趨勢.需要指出,文獻[13]基于CEFI 公式[9],采用有限元方法統(tǒng)一求解流體和結構,層合梁的非線性建模采用了Saint Venant Kirchhoff線性超彈性材料和Green-Lagrangin 位移?應變關系.而本文采用分區(qū)強流固耦合方法,流體部分采用有限體積法求解,結構部分采用有限元方法求解,層合梁的幾何非線性建模采用von Kármán 位移?應變關系,導致本文計算得到的層合梁振幅和頻率結果與文獻解略有差別.

圖5 不同 ? 情況下,層合梁右端(x=L)節(jié)點橫向位移的最大振幅Fig.5 Maximum amplitude of transverse displacement at laminated beam tip as a function of ?

圖6 不同 ? 情況下,層合梁右端(x=L)節(jié)點橫向位移的主導頻率Fig.6 Dominant frequency of transverse displacement at laminated beam tip as a function of ?

3.2 剪切流速度分布對單層梁振動的影響

本節(jié)研究剪切流的速度分布對單層梁振動特性的影響.考慮剪切流不同的來流速度分布 α ∈[0,2],采用如下無量綱參數:Re=1000,αavg=10,βavg=6000,h=0.01.來流速度分布與系數 α 的關系如圖7 所示.

圖7 不同 α 情況下的來流速度分布Fig.7 Velocity profile as a function of α

圖8 給出了剪切流中單層梁右端節(jié)點10 個運動周期的橫向位移標準差和平均值隨 α 變化的曲線.圖9 給出了振動頻率隨 α 變化的曲線.結果表明,隨著剪切流速度分布系數 α 的增加,單層梁振動的幅度和主導頻率(渦脫落頻率)均是先減小后增加,單層梁的振動平均值隨著 α 的增加逐漸增大,振動的向下偏斜更加明顯.

圖8 單層梁右端(x=L)節(jié)點的橫向位移標準差曲線和平均值曲線Fig.8 Standard deviation and mean value of transverse displacement at single isotropic beam tip as a function of α

圖9 單層梁右端(x=L)節(jié)點的橫向振動頻率曲線Fig.9 Dominant frequency of transverse displacement at single isotropic beam tip as a function of α

圖10 給出不同 α 情況下單層梁在一個完整運動周期內的變形包絡圖,紅色虛線為中性軸,圖10(c)標記了單層梁拍動過程中彎曲變形的三個拐點(x4=0.175L,x4=0.475L,x3=0.775L) 和右端節(jié)點x4=L.圖11 給出了 α=1 情況下上述四個節(jié)點在5 個振動周期內的橫向位移時域曲線,位移曲線峰值的依次出現(xiàn)表明單層梁的振動模式區(qū)別于類模態(tài)振動模式,是一種行波形式的振動模式[32].圖12 給出不同 α 情況下,20 個振動周期內,單層梁右端節(jié)點位移的Lissajous 曲線,紅色虛線表示橫向位移的平均值.結果進一步表明,隨著剪切流速度分布系數 α的增加,單層梁的振動平均偏斜量增加,右端節(jié)點的運動軌跡發(fā)生改變,當 α=2 時(圖12),Lissajous 曲線的不對稱性變得明顯.

圖10 單層梁在一個完整振動周期內的變形包絡圖Fig.10 Deformation envelope of single isotropic beam in a complete motion period

圖11 單層梁不同位置節(jié)點在5 個振動周期內的橫向位移時域曲線(α=1)Fig.11 Time domain responses of transverse displacement of beam at different positions in 5 flapping motion periods (α=1)

圖12 不同 α 情況下,單層梁在20 個完整振動周期內的右端節(jié)點位移Lissajous 曲線Fig.12 Lissajous curves of single isotropic beam in 20 flapping motion periods at beam tip with different α

圖13～圖15 分別給出了剪切速度分布 α=0,α=0.5,α=1.8 情況下,單層梁一個完整振動周期內的流場渦量圖.在α=0 情況下觀察到典型的2 S 渦模態(tài)(每個振動周期脫落兩個單獨的渦,兩個渦的強度幾乎相等,旋轉方向相反);隨著剪切速度分布系數 α 增加,單層梁后方流場中脫落的漩渦不再呈現(xiàn)對稱趨勢,在α=0.5 和 α=1.8 情況下,上方的漩渦尺寸較大且形狀較圓,下方的漩渦尺寸較小且形狀較細長.

圖13 一個完整振動周期內的渦量圖(α=0)Fig.13 Vorticity diagrams in a complete flapping motion period at α=0

圖14 一個完整振動周期內的渦量圖(α=0.5)Fig.14 Vorticity diagrams in a complete flapping motion period at α=0.5

圖15 一個完整振動周期內的渦量圖(α=1.8)Fig.15 Vorticity diagrams in a complete flapping motion period at α=1.8

3.3 剛度比對剪切流中層合梁振動的影響

本節(jié)研究剪切流中層合梁上下兩層彈性模量的差異對振動響應的影響,考慮不同剛度比 ?=2.4 ×10?3,4.8 × 10?3,7.2 × 10?3,9.6 × 10?3.圖16 給出了10 個無量綱時長區(qū)間內,不同剛度比的層合梁右端(x=L)節(jié)點的橫向位移時域曲線,可以發(fā)現(xiàn)隨著剛度比 ? 的增加,雙層層合梁振動的振幅增大,主導頻率下降,這一規(guī)律與均勻流作用下層合梁的振動特性變化規(guī)律一致.

圖16 層合梁右端(x=L)節(jié)點的橫向位移時域響應Fig.16 Time domain response of transverse displacement at beam tip

進一步對比上下兩層彈性模量的差異對振動響應運動軌跡的影響.不同剛度比情況下,20 個振動周期內,層合梁右端節(jié)點位移的Lissajous 曲線如圖17 所示,紅色虛線表示橫向位移的平均值.當剛度比較小時,Lissajous 曲線上下對稱,運動軌跡是‘8’字形;當剛度比 ? 較大時可以發(fā)現(xiàn)Lissajous 曲線不對稱的現(xiàn)象,層合梁向下振動的彎曲程度大于向上振動的彎曲程度.

圖17 層合梁右端(x=L)節(jié)點位移的Lissajous 曲線Fig.17 Lissajous curves of displacements at beam tip

3.4 厚度占比對剪切流中層合梁振動的影響

本節(jié)研究剪切流中層合梁上下兩層的不同厚度占比對振動特性的影響,定義上鋪層厚度占比 γ=Ht/(Hb+Ht).考慮剪切速度分布 α=1,分析了剛度比 ?=0.0048 和 ?=0.0096 兩種情況下,上鋪層厚度 γ ∈[0,1] 的層合梁振動特性,其他無量綱參數與3.2 節(jié)相同.

圖18 和圖19 分別給出了兩種剛度比情況下,10 個振動周期內,層合梁右端節(jié)點橫向位移標準差和平均值曲線隨厚度占比 γ 變化的曲線.圖18 表明,在γ <0.8 的區(qū)間內,隨著厚度比的增加,層合梁的振動幅度逐漸降低;在0.8 <γ <0.9 的區(qū)間內,層合梁的振動從大幅度周期極限環(huán)振動模式(P)轉變到靜變形為主的定點穩(wěn)定模式(S),兩種模式下的變形包絡圖在圖18 中給出.圖19 表明,隨著剛度較大的上鋪層厚度占比的增加,剪切流作用下層合梁的平均偏斜量減少.結合圖18 和圖19 中的結果分析,厚度占比 γ=0 和 γ=0.1 時,剛度比差異大(? ＝ 0.0096)的層合梁平均偏移量更大,但振幅更小;在0.2 <γ <0.8 的區(qū)間內,剛度比差異大(? ＝ 0.0096)的層合梁平均偏移量和振幅均大于剛度比差異小的層合梁;當 γ >0.8 時,層合梁處于定點穩(wěn)定變形模式,剛度比大的層合梁的平均偏斜量小.

圖18 層合梁右端(x=L)節(jié)點的橫向位移標準差曲線以及 γ=0.33,γ=1 時的變形包絡圖Fig.18 Standard deviation of transverse displacement at beam tip as a function of γ and deformation envelope of beam at γ=0.33,γ=1,respectively

圖19 層合梁右端(x=L)節(jié)點的橫向位移平均值曲線Fig.19 Mean value of transverse displacement at beam tip as a function of γ

結果表明,層合材料對復合材料梁的動力學行為有顯著影響,相同的流場條件下層合特性的輕微改變會導致動力學行為的突變,這可以歸因于厚度比的變化影響了層合梁的等效剛度,進而影響了剪切流?層合梁系統(tǒng)的穩(wěn)定性.

3.5 鋪層角度對剪切流中層合梁振動的影響

本節(jié)研究剪切流中層合梁上下兩層厚度相同,材料相同,不同的鋪層角度對振動特性的影響,采用與3.2 節(jié)相同的無量綱參數,各向異性復合材料參數如下:E2=E1/10,G23=G13=G12=E1/20 .定義鋪層角度 θ,層合梁的鋪層方式表示為 [ θ,?θ] .考慮剪切速度分布 α=1,分析了鋪層夾角 θ ∈[0?,90?] 的復合材料層合梁振動特性.

圖20 給出了20 個振動時長內層合梁右端節(jié)點橫向位移標準差值曲線隨鋪層角度 θ 變化的曲線,其中藍色散點表示不同角度 θ 情況下20 個振動時長內的位移曲線幅值極大值.結果表明,隨著鋪層角度的增加,層合梁的振動模式從周期極限環(huán)振動模式(P)轉變到非周期振動模式(NP).非周期振動模式在大多數情況下無法識別極限環(huán),通常是概周期或混沌運動,在這種非周期狀態(tài)下,由于運動的隨機性導致了位移曲線的幅值極大值具有不穩(wěn)定性;另一方面,位移的標準差曲線保持了連續(xù)性,非周期振動模式與周期極限環(huán)振動模式下的標準差具有相同的數量級.圖21 給出了20 個振動時長內,層合梁右端節(jié)點橫向位移平均值曲線隨鋪層角度 θ 變化的曲線,θ=35°和 θ=40°情況下右端節(jié)點的Lissajous 曲線也在圖21 中給出.隨著鋪層角度的增加,層合梁的等效剛度減少,導致剪切流作用下層合梁的平均偏斜量增大.

圖20 層合梁右端(x=L)節(jié)點的橫向位移標準差曲線以及 θ=35°,θ=40°時的時域曲線Fig.20 Standard deviation of transverse displacement at beam tip as a function of θ and time domain response at θ=35°,θ=40°,respectively

圖21 層合梁右端(x=L)節(jié)點的橫向位移平均值曲線以及 θ=35°,θ=40°時的Lissajous 曲線Fig.21 Mean value of transverse displacement at beam tip as a function of θ and Lissajous curves at θ=35°,θ=40°,respectively

4 結論

本文對剪切流作用下復合材料層合梁的振動特性進行數值建模和計算,得到了以下結論.

(1)剪切流作用下單層梁的振動特性與均勻流作用下不同,梁的振動受剪切流影響向下偏斜,隨著速度分布系數增加,單層梁振動的不對稱性逐漸明顯,振動的幅度和主導頻率(渦脫落頻率)均是先減小后增加,尾部流場中的渦結構發(fā)生改變.

(2)剛度比對剪切流作用下層合梁的振動特性有顯著影響.隨著剛度比的增加,層合梁振動的振幅增大,主導頻率下降,運動軌跡由“8”字形逐漸變得不對稱,層合梁向下振動的彎曲程度大于向上振動的彎曲程度.

(3)厚度占比對剪切流作用下層合梁的振動特性有顯著影響.隨著厚度占比的增加,觀察到兩種不同的響應狀態(tài):大幅度周期極限環(huán)振動模式(P)轉變到靜變形為主的定點穩(wěn)定模式(S).

(4)復合材料鋪層角度對剪切流作用下層合梁的振動特性有顯著影響.隨著鋪層角度的增加,觀察到振動模式由周期極限環(huán)振動模式(P)向非周期振動模式(NP)轉變.