大型空间薄膜结构由于其质量轻、收藏体积小、成本低、易折叠等优点得到广泛关注,然而在张拉薄膜结构中通常存在褶皱[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内,验证了控制方法对抑制振动的有效性。