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

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

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

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

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

薄膜的應(yīng)變勢(shì)能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方向的應(yīng)變,γxy是切應(yīng)變。它們與振動(dòng)構(gòu)形函數(shù)和褶皺變形函數(shù)的關(guān)系可以表示為:

將式(3)代入到式(2),得到薄膜應(yīng)變勢(shì)能的表達(dá)式為:

薄膜自由振動(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均表示成級(jí)數(shù)形式[11]:

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

$\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)能的表達(dá)式:

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

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

$\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)代入拉格朗日方程,最終得到褶皺薄膜結(jié)構(gòu)的非線性振動(dòng)微分方程:

當(dāng)膜面存在褶皺時(shí),膜面上各點(diǎn)的振動(dòng)平衡位置不再是0,根據(jù)上一節(jié)建立的有限元模型,當(dāng)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)施加節(jié)點(diǎn)號(hào)和褶皺峰值P點(diǎn)節(jié)點(diǎn)號(hào)對(duì)應(yīng)的振動(dòng)值,通過(guò)/OUTPUT指令將時(shí)間與振動(dòng)值保存到文本文檔中,繪制得到圖4所示的擾動(dòng)點(diǎn)處振動(dòng)曲線圖,可以看到,在1.3 s之前該點(diǎn)處于靜止?fàn)顟B(tài),1.3 s時(shí)出現(xiàn)擾動(dòng)信號(hào),之后膜面開始振動(dòng),振動(dòng)頻率約為28.5 Hz,與圖2中的一階固有頻率值相吻合。

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

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

$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為比例系數(shù),Ki為積分系數(shù),Kd為微分系數(shù),e(t)=y(t)-yr(t)為輸出誤差,壓電作動(dòng)器輸入電壓ua與壓電作動(dòng)器輸出位移u之間近似成線性關(guān)系:

u=Kaua (11)

式中:Ka為壓電常數(shù)。

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

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