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

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

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

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

薄膜結構自由振動時,僅考慮面外振動,振動構形函數為wv,考慮褶皺的存在,引入褶皺變形函數wz,則振動過程中薄膜彎曲形變勢能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是切應變。它們與振動構形函數和褶皺變形函數的關系可以表示為:

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

薄膜自由振動時的動能為:

$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),得到振動過程中總勢能為:

$\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),得到動能的表達式:

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

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

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

如圖6所示,在薄膜結構頂點附近安裝壓電作動器,通過垂直膜面的運動有效抑制薄膜振動。由于薄膜結構振動模型是高度耦合的非線性方程,很難基于模型設計控制器,因此在ANSYS中建立有限元模型進行瞬態分析作為被控對象,在MATLAB軟件中設計不基于模型的PID控制器,如圖7所示,yr(t)為參考輸入信號,y(t)為ANSYS模型后處理提取的P點處振動量。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)為輸出誤差,壓電作動器輸入電壓ua與壓電作動器輸出位移u之間近似成線性關系:

u=Kaua (11)

式中:Ka為壓電常數。

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

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