大型空間薄膜結構由于其質(zhì)量輕、收藏體積小、成本低、易折疊等優(yōu)點(diǎn)得到廣泛關(guān)注,然而在張拉薄膜結構中通常存在褶皺[1],褶皺的存在不僅會(huì )降低薄膜結構的表面精度,還會(huì )改變薄膜中的應力分布,影響薄膜結構的動(dòng)態(tài)性能[2,3]。已有文獻中,KUKATHASAN et al[4]通過(guò)實(shí)驗分析了不同褶皺形態(tài)下薄膜振動(dòng)特性的變化,驗證了褶皺的存在可改變薄膜結構的振型和固有頻率;HOSSAIN et al[5]在A(yíng)BAQUS軟件中引入罰參數修正材料模型,說(shuō)明褶皺對薄膜結構的振動(dòng)特性有影響。因此,對薄膜結構進(jìn)行動(dòng)態(tài)特性分析時(shí),應當考慮褶皺引起的質(zhì)量矩陣和剛度矩陣的變化。

本文在A(yíng)NSYS有限元軟件中采用殼單元建立正方形薄膜結構,首先利用非線(xiàn)性有限元屈曲分析得到薄膜結構的褶皺形態(tài)[6],引入褶皺信息更新有限元模型進(jìn)行模態(tài)分析,對比有無(wú)褶皺2種情況下薄膜結構各階固有頻率和模態(tài)振型。以初始褶皺形態(tài)作為振動(dòng)平衡位置,在薄膜表面一點(diǎn)施加垂直于膜面方向的沖擊擾動(dòng)進(jìn)行瞬態(tài)分析[7,8],得到褶皺薄膜的振動(dòng)曲線(xiàn)。在此基礎上,設計PID控制系統,通過(guò)MATLAB/ANSYS聯(lián)合仿真驗證控制方法對抑制振動(dòng)的有效性,從而提高薄膜結構的在軌穩定性[9,10]。

建立圖1所示正方形(邊長(cháng)L=500 mm)薄膜結構的有限元模型,薄膜楊氏模量E=2.5 GPa,泊松比ν=0.34,密度ρ=1 400 kg/m3,厚度h=25 μm。為避免應力集中、求解不收斂的情況,將4個(gè)角切除,L1=25 mm。邊界條件為四角固支,在4個(gè)頂點(diǎn)處的直線(xiàn)上施加沿法線(xiàn)方向的位移載荷VX=1.0 mm,VY=0.25 mm。

根據已有研究成果,當X方向和Y方向拉力比值大于等于4時(shí)形成貫穿中心的橫向褶皺[6],更新模型并進(jìn)行模態(tài)分析,去除噪聲振型后提取前4階固有頻率和模態(tài)振型,如圖2所示。從振型圖上即可看到褶皺形態(tài),說(shuō)明褶皺對薄膜結構的動(dòng)態(tài)特性產(chǎn)生了影響。

薄膜結構自由振動(dòng)時(shí),僅考慮面外振動(dòng),振動(dòng)構形函數為wv,考慮褶皺的存在,引入褶皺變形函數wz,則振動(dòng)過(guò)程中薄膜彎曲形變勢能Ep1為:

薄膜的應變勢能Ep2可以表示為:

$\begin{array}{l}E{}_{\text{p}2}=\frac{Eh{}_m}{2(1-\nu {}^2)}{\displaystyle \iint (}\epsilon {}_x^2+\epsilon {}_y^2+\\ 2\nu \epsilon {}_x\epsilon {}_y+\frac{1-\nu }{2}\gamma {}_{xy}^2)\text{d}x\text{d}y(2)\end{array}$

式中:εx,εy分別是X方向和Y方向的應變,γxy是切應變。它們與振動(dòng)構形函數和褶皺變形函數的關(guān)系可以表示為:

將式(3)代入到式(2),得到薄膜應變勢能的表達式為:

薄膜自由振動(dòng)時(shí)的動(dòng)能為:

$E{}_{\text{v}}=\frac{1}{2}\rho h{}_m{\displaystyle \iint (}\frac{\partial w{}_{\text{v}}}{\partial t}){}^2\text{d}x\text{d}y$ (5)

將wv和wz均表示成級數形式[11]:

將式(6)分別代入式(1)和(4),得到振動(dòng)過(guò)程中總勢能為:

$\begin{array}{l}E{}_{\text{p}}=E{}_{\text{p}1}+E{}_{\text{p}2}=\\ \frac{1}{2}(k{}_{ijkl}q{}_{ijkl}+k{}_{ij}q{}_{ij}+k^{\prime }{}_{ijkl})(7)\end{array}$

將式(6)代入式(5),得到動(dòng)能的表達式:

$E{}_{\text{v}}=\frac{1}{2}m{}_{ij}\dot{q}{}_i\dot{q}{}_j$ (8)

式中:kijkl,kij,k′ijkl,mij的具體表達式如下:

$\begin{array}{l}k{}_{ijkl}=\frac{Eh{}_m}{4(1-\nu {}^2)}{\displaystyle \iint (}\frac{\partial w{}_{\text{v}i}}{\partial x}\frac{\partial w{}_{\text{v}j}}{\partial x}\frac{\partial w{}_{\text{v}k}}{\partial x}\frac{\partial w{}_{\text{v}l}}{\partial x}+\\ \frac{\partial w{}_{\text{v}i}}{\partial y}\frac{\partial w{}_{\text{v}j}}{\partial y}\frac{\partial w{}_{\text{v}k}}{\partial y}\frac{\partial w{}_{\text{v}l}}{\partial y}+2\frac{\partial w{}_{\text{v}i}}{\partial x}\frac{\partial w{}_{\text{v}j}}{\partial x}\frac{\partial w{}_{\text{v}k}}{\partial y}\frac{\partial w{}_{\text{v}l}}{\partial y})\text{d}x\text{d}y\\ k{}_{ij}=\text{Δ}{\displaystyle \iint (}\frac{\partial {}^{2}w{}_{\text{v}i}}{\partial x{}^2}\frac{\partial {}^{2}w{}_{\text{v}j}}{\partial x{}^2}+\frac{\partial {}^{2}w{}_{\text{v}i}}{\partial y{}^2}\frac{\partial {}^{2}w{}_{\text{v}j}}{\partial y{}^2}+\\ 2\nu \frac{\partial {}^{2}w{}_{\text{v}i}}{\partial x{}^2}\frac{\partial {}^{2}w{}_{\text{v}j}}{\partial y{}^2}+2(1-\nu )\frac{\partial {}^{2}w{}_{\text{v}i}}{\partial x\partial y}\frac{\partial {}^{2}w{}_{\text{v}j}}{\partial x\partial y})\text{d}x\text{d}y-\\ \frac{Eh{}_m}{2(1-\nu {}^2)}{\displaystyle \iint (}\frac{\partial w{}_{\text{v}i}}{\partial x}\frac{\partial w{}_{\text{v}j}}{\partial x}\frac{\partial w{}_{\text{z}i}}{\partial x}\frac{\partial w{}_{\text{z}j}}{\partial x}+\\ \frac{\partial w{}_{\text{v}i}}{\partial y}\frac{\partial w{}_{\text{v}j}}{\partial y}\frac{\partial w{}_{\text{z}i}}{\partial y}\frac{\partial w{}_{\text{z}j}}{\partial y}+\nu \frac{\partial w{}_{\text{v}i}}{\partial x}\frac{\partial w{}_{\text{v}j}}{\partial x}\frac{\partial w{}_{\text{z}i}}{\partial y}\frac{\partial w{}_{\text{z}j}}{\partial y}+\\ \nu \frac{\partial w{}_{\text{v}i}}{\partial y}\frac{\partial w{}_{\text{v}j}}{\partial y}\frac{\partial w{}_{\text{z}i}}{\partial x}\frac{\partial w{}_{\text{z}j}}{\partial x}+\\ 2(1-\nu )\frac{\partial w{}_{\text{v}i}}{\partial x}\frac{\partial w{}_{\text{z}i}}{\partial x}\frac{\partial w{}_{\text{v}j}}{\partial y}\frac{\partial w{}_{\text{z}j}}{\partial y})\text{d}x\text{d}y\\ k^{\prime }{}_{ijkl}=\frac{Eh{}_m}{4(1-\nu {}^2)}{\displaystyle \iint (}\frac{\partial w{}_{\text{z}i}}{\partial x}\frac{\partial w{}_{\text{z}j}}{\partial x}\frac{\partial w{}_{\text{z}k}}{\partial x}\frac{\partial w{}_{\text{z}l}}{\partial x}+\\ \frac{\partial w{}_{\text{z}i}}{\partial y}\frac{\partial w{}_{\text{z}j}}{\partial y}\frac{\partial w{}_{\text{z}k}}{\partial y}\frac{\partial w{}_{\text{z}l}}{\partial y}+2\frac{\partial w{}_{\text{z}i}}{\partial x}\frac{\partial w{}_{\text{z}j}}{\partial x}\frac{\partial w{}_{\text{z}k}}{\partial y}\frac{\partial w{}_{\text{z}l}}{\partial y})\text{d}x\text{d}y\\ m{}_{ij}=\rho h{}_m{\displaystyle \iint w}{}_{\text{v}i}w{}_{\text{v}j}\text{d}x\text{d}y\end{array}$

將式(7)-(8)代入拉格朗日方程,最終得到褶皺薄膜結構的非線(xiàn)性振動(dòng)微分方程:

當膜面存在褶皺時(shí),膜面上各點(diǎn)的振動(dòng)平衡位置不再是0,根據上一節建立的有限元模型,當VX=1.0 mm,VY=0.25 mm時(shí)得到圖3所示褶皺形態(tài),各點(diǎn)褶皺幅值即為無(wú)振動(dòng)時(shí)的平衡位。進(jìn)行瞬態(tài)分析,在1.3 s時(shí)施加沖擊擾動(dòng),擾動(dòng)幅值為4 mm,作用時(shí)間為0.02 s,仿真總時(shí)間為4 s。在后處理/POST26中提取擾動(dòng)施加節點(diǎn)號和褶皺峰值P點(diǎn)節點(diǎn)號對應的振動(dòng)值,通過(guò)/OUTPUT指令將時(shí)間與振動(dòng)值保存到文本文檔中,繪制得到圖4所示的擾動(dòng)點(diǎn)處振動(dòng)曲線(xiàn)圖,可以看到,在1.3 s之前該點(diǎn)處于靜止狀態(tài),1.3 s時(shí)出現擾動(dòng)信號,之后膜面開(kāi)始振動(dòng),振動(dòng)頻率約為28.5 Hz,與圖2中的一階固有頻率值相吻合。

圖5為褶皺峰值P點(diǎn)處的振動(dòng)曲線(xiàn),1.3 s前該點(diǎn)處初始平衡位置為褶皺變形量0.80 mm。

如圖6所示,在薄膜結構頂點(diǎn)附近安裝壓電作動(dòng)器,通過(guò)垂直膜面的運動(dòng)有效抑制薄膜振動(dòng)。由于薄膜結構振動(dòng)模型是高度耦合的非線(xiàn)性方程,很難基于模型設計控制器,因此在A(yíng)NSYS中建立有限元模型進(jìn)行瞬態(tài)分析作為被控對象,在MATLAB軟件中設計不基于模型的PID控制器,如圖7所示,yr(t)為參考輸入信號,y(t)為ANSYS模型后處理提取的P點(diǎn)處振動(dòng)量。PID控制器輸出表達式為:

$u{}_a={\rm K}{}_{\text{p}}e(t)+{\rm K}{}_{\text{i}}{\displaystyle {\int }_0^{t}e}(t)\text{d}t+{\rm K}{}_{\text{d}}\frac{\text{d}}{\text{d}t}e(t)$ (10)

式中:Kp為比例系數,Ki為積分系數,Kd為微分系數,e(t)=y(t)-yr(t)為輸出誤差,壓電作動(dòng)器輸入電壓ua與壓電作動(dòng)器輸出位移u之間近似成線(xiàn)性關(guān)系:

u=Kaua (11)

式中:Ka為壓電常數。

在MATLAB軟件中編寫(xiě)m文件,設置Kp,Ki,Kd以及Ka的值。進(jìn)行MATLAB與ANSYS聯(lián)合仿真時(shí),首先在A(yíng)NSYS軟件中進(jìn)行褶皺薄膜瞬態(tài)分析,利用*GET命令提取P點(diǎn)振動(dòng)值,然后用*VWRITE命令把數據寫(xiě)入由*CFOPEN打開(kāi)的中轉文本文件vibration.txt中,該文本文檔中的數據在每次仿真后會(huì )被覆蓋,在數據變化之前,由MATLAB通過(guò)textread命令讀取結果并進(jìn)行運算,運用save命令將計算得到的壓電作動(dòng)器輸出位移保存到文本文檔u.txt中,之后ANSYS采用*VREAD命令讀取u.txt文檔中的數據,引入有限元模型中,并進(jìn)入下一次瞬態(tài)分析。程序自動(dòng)按照先后順序進(jìn)行運算,設定運行時(shí)間為4 s,仿真結果如圖8所示,認為輸出誤差在±0.01 mm內即滿(mǎn)足要求,可以得到振動(dòng)時(shí)間約為3.17 s,說(shuō)明控制系統可有效抑制薄膜振動(dòng)。

通過(guò)仿真結果對比得到如下結論:引入褶皺后,在振型圖上可以清晰地看到褶皺形態(tài);從褶皺峰值點(diǎn)處振動(dòng)曲線(xiàn)可以看到,其平衡位置點(diǎn)即為褶皺變形值,說(shuō)明褶皺的存在改變了初始平衡狀態(tài);通過(guò)控制后的振動(dòng)曲線(xiàn)可以看到,振動(dòng)時(shí)間為3.17 s時(shí)將振幅控制在±0.01 mm內,驗證了控制方法對抑制振動(dòng)的有效性。