2. 中国工程物理研究院 核物理与化学研究所, 四川 绵阳 621000;
3. 空间物理重点实验室, 北京 100071
2. Institute of Nuclear Physics and Chemistry, Chinese Academy of Engineering Physics, Mianyang Sichuan 621000, China;
3. Key Laboratory of Space Physics, Beijing 100076, China
飞机在含有过冷水滴的云层中飞行时,飞机迎风面处会有结冰产生。结冰问题会直接威胁到飞机的飞行安全[1, 2, 3, 4, 5, 6],因此研究可行的飞机防/除冰系统是必要的。防/除冰系统可以分为热力防/除冰系统、机械防/除冰系统和化学防/除冰系统三种类型。热力防/除冰系统被用于大多数的防/除冰系统,但是其能耗较大(电热式),或者对发动机较大的引气量会影响发动机的性能(热气式)。化学防/除冰系统只能用于地面防除冰过程,并且大量的化学除冰物质会对环境有影响,另外其防除冰效率不高,不适用于飞行状态下的防除冰过程。机械防/除冰系统主要包含电脉冲除冰方法[7]、电斥除冰方法[8]、气囊除冰方法[9]和压电除冰方法等。其中电脉冲和电斥的除冰方法在能耗上较大,同时瞬时较大的冲击力会对机翼蒙皮造成损害;气囊除冰方法利用机翼蒙皮的不断变形进行除冰,这也会对机翼蒙皮的寿命产生影响;因此近年来在防/除冰系统“效率高、能耗低、重量轻、对气动外形影响小”的要求下,压电除冰方法得到了重视和不断研究。
目前国内外针对压电除冰方法的研究主要分为两个方面:一个是利用超声频率范围内的超声除冰方法[10, 11, 12, 13],另外一种是在低频范围内进行的压电振动除冰方法。针对低频下的压电振动除冰方法研究而言,印度科学理工学院的 Kandagal 等[14]详细研究了平板振动的振型特性,并且通过具体实验观测了除冰效果;艾克朗大学的 Venna 等[15]通过对机翼前缘数值模拟和相关节点数据的研究,得到特定压电元件尺寸和排布下的除冰最佳频率,然后实验验证了压电振动除冰方法的可行性;国内在压电振动除冰方面的研究处于起步阶段,南京航空航天大学的严春霞[16]研究了基于弹性波的机翼振动除冰机理;南京航空航天大学的姚赛金[17]也在基于压电驱动器的基础上研究了机翼振动除冰方法,利用有限元分析手段得到平板和翼型前缘结构上的压电振动特性,然后进行了平板和翼型的压电振动除冰实验,得到较好的除冰效果;南京航空航天大学的白天等[18]对平板上的压电振动除冰方法进行了细致研究,通过针对探索压电元件的尺寸与所激发振动模态的关系,以及振动模态与除冰剪应力之间的关系来进行压电除冰的设计和研究,取得了不错成果。
目前的压电振动除冰方法的研究过程可以认为有正向和反向两种研究模式[18]。正向的研究模式为:把压电元件的分布位置作为确定量,研究在该分布位置下的最佳振动频率。该研究模式是现今绝大多数学者采用的方法。反向的研究模式为:把需要进行除冰的位置作为已知量,研究针对特定除冰位置下的最佳模态振型,然后依据最佳振型研究合适的压电元件尺寸和布局方式。这种研究模式更加贴近于实际应用。本文的研究即采用反向的研究模式。由于白天等的研究主要针对于平板结构,和实际机翼相比有很大的差距,因此本文主要针对缩比后的机翼前缘翼型结构进行压电振动除冰方法的研究,对翼型结构的相关几何处理、压电元件的贴片位置和翼型上的压电元件排布规律等内容进行了研究,最后通过具体地面实验进行除冰效果的验证。
1 压电元件的可贴范围由于翼型结构为曲面结构,而压电元件为平直的薄板结构,因此在一定的曲率下,压电元件和翼型结构之间会存在一定的间隙。较小的间隙可以通过胶层的填充而实现压电元件的粘接过程,较大的间隙下通过胶层的填充来实现压电元件的粘接过程,会使压电元件的振动激励效果受到影响,因此需要对压电元件在曲面翼型上的可贴位置和范围进行确定。
对于研究用翼型结构,其展向各个位置的截面均一致,因此研究某一截面下的翼型曲线上的贴片范围,即可确定整体翼型曲面上的贴片范围。而对于NACA0012翼型而言,其翼型的描述函数同实际数据点差距较大,通过解析的方式求取相应的弧长是不可行的,因此需要借助有限元的离散数据点来求解压电元件的相关贴片位置。
1.1 粘接层厚度的相关处理由于间隙中需要填充胶层,因此需要得到满足振动要求下的间隙厚度,保证压电元件在此间隙范围内能够达到较好的激励效果。Overmeyer的研究结果表明,粘接层厚度要大于0.381mm的情况下,才能满足较好的振动效果[12]。因此借鉴其研究过程,设定压电元件和结构之间的间隙要小于0.381mm,这样才能保证有更好的贴合效果。
由于椭圆曲线方程形式简洁,并且为连续曲线,因此小尺寸的计算过程中,在确定的间隙和压电元件边长的前提下,近似采用椭圆形弧长来代替描述机翼曲线形状,可以进行最大边长与弧长差值的计算,并且该结果可以作为编程计算的参考结果。依据椭圆公式,确定的压电元件边长pl和确定的间隙s可以得到如图 1所示的椭圆的方程
由误差为十万分之三的椭圆圆周公式,可以得到确定压电元件边长下的椭圆的弧长:
由于选取的两种方形压电元件的边长分别为10mm和15mm,因此得到对应下的弧长和边长的差值分别为1.005×10-4m和1.3949×10-3m。计算时设置的弧长和边长的差值为1×10-5m,满足了研究的精度要求。
1.2 具体算法过程翼型的有限元离散点通过有限元软件导出,网格单元长度为0.001m。提取出的离散点为翼型中间截面的一系列坐标点,这些坐标点可以描述该截面翼型曲线。求解算法以这些坐标点(按照曲线顺序进行了排列)为基础进行开展。
关于位置的计算过程中,涉及间距和弧长的求解。间距的求解采用简单的两点距离求解方法,而弧长采用近似求解方法,即在组成一段曲线的所有坐标点中,通过邻近两点直线距离的叠加得到整个弧长,如式(3)所示:
大致的算法过程如图 2所示。由于翼型曲线具有分布规律的坐标点,因此可以认为曲线是连续的,因此计算压电元件可贴区域时采用弧长和弦长对比的方式来实现。而压电元件的某一可贴位置可以认为是满足压电元件粘接条件的两个端点坐标。
对于压电元件而言,设定压电元件边长所代表的二维线段端点分别为贴片起始点A和终止点B。从压电元件的贴片起始点A开始,参照已知的按顺序排列的节点坐标,分别求解A点和沿曲线顺序增大方向的节点之间的距离,记为弦长chord,当某一弦长同压电元件长度非常接近时,满足式(4),即可认定该弦长近似代表压电元件的长度,而此时弦长对于的节点视为压电元件的终止点B。由于压电元件的规格和可贴位置的定位尺寸单位为毫米,并且实际压电元件贴片时的贴片位置也是允许有1~2mm的误差,所以当压电元件的终止点B不为机翼曲面上某节点时,可以近似选取最接近的节点作为压电元件的终止点。
此时压电元件的计算边长为chord,而可贴位置对应的曲线上弧长(A点和B点之间的曲线长度)为arc。
对比A、B两点之间计算弧长arc和压电元件计算边长chord的差距,设定一个参考值suit,当满足公式(6)时,即可认为该位置为可贴压电元件位置。
从曲线起始位置开始向曲线终点方向逐点进行判定,即可得到不同尺寸的压电元件在翼型曲面上的可贴区域。
1.3 具体算例对于生成的有限元分析模型,提取了翼型展向起始截面的节点坐标信息,并按照曲线生成方向排列了数据点。利用算法计算了边长为10mm和15mm下压电元件的压电元件可贴范围,结果如图所示。计算结果表明,边长为10的压电元件在上翼面实际的可贴压电元件范围为(0.0162,0.08),边长为15的压电元件在上翼面实际的可贴压电元件范围为(0.02696,0.08)。该结果为实际条件下压电元件的排布有工程上的限定作用。
2 翼型相关几何关系处理翼型曲面不同于一般规则曲面,其无法用准确的方程来描述,所以曲面某一点处的切向向量和法向向量无法直接得到,并且由于有限元分析软件中的坐标系一般为直角坐标系、柱坐标系和球坐标系,任意一种坐标系都不能够直接得到翼型曲面切法向的相关分量。另外,有限元分析软件中建立的局部坐标系仅对某点有作用,要得到整体翼型曲面上各点的切法向分量而建立大量的局部坐标系是不可行的。因此,借鉴第1部分中相关内容,利用有限元模型导出的离散数据点为基础,采取算法来求解每一点对应的切法向的方向矢量,进而通过每点对应切法向的方向矢量来求解翼型曲面上的相关量。
2.1 理论基础Tookey和Ball[19]研究发现,对于组成样条曲线的离散点而言,利用三点共圆法来求解曲线上某点的切法向方向矢量的误差是比较大的,他们引入了五点共圆法来提高求解精度。
如图 4所示,五点共圆法针对连续分布的五个离散点,分别选择两组点(r0,r2,r4和r1,r2,r3)确定两个法向矢量ka,kb,然后两个法向矢量通过Richardson外推法求得目标点的法向矢量。
两组点确定的法向矢量分别为:
其中: 若求ka,则A=(r1-r2),B=(r3-r2),C=(r1-r3);若求kb,则A=(r0-r2),B=(r4-r2),C=(r1-r4)
对于目标点的法向矢量,则有:
五点共圆法适用于整个曲线上除去端点之外的内部离散点。由于对端点的处理比较欠缺,因此Ma和Cripps[20]提出了利用四点共圆法求解曲线端点处法向矢量,如图 5所示。
方法与五点共圆法一致,法向矢量的公式如式(7)所示。其中若求ka,则A=(r1-r0),B=(r2-r0),C=(r1-r2);若求kb,则A=(r1-r0),B=(r3-r0),C=(r1-r2)。
对于目标点的法向矢量,则有:
由于翼型展向无形状变化,仍然提取某一截面上的数据点来进行计算。针对翼型数据点组成的翼型截面曲线,利用五点共圆法处理中间点和四点共圆法处理端点的组合形式来求解切法向矢量。
对于翼型中间节点:
对于翼型曲线起始端点,:
对于翼型曲线末端端点:
采用该算法计算翼型曲线的各点切法向分量,并借助切向分量求得距离原始曲线一定切向距离下的新翼型曲线,如图 6所示,整体结果比较相似,证明该算法对于处理翼型切法向几何关系是可行的。另外,三点法与五点/四点共圆法的计算结果比较接近,但是在实际算法编程计算中,由于处理的数据数量级较小,三点法对应的相关计算表达式容易出现分母为零的现象,从而造成数据错误,需要进行大量的数据类型之间的变换;相比之下,五点/四点共圆法的表达式采用基本向量运算的方式,更适用于进行编程计算,因此在求解算法中具有更好的应用优势。
3 翼型上压电元件布局规律研究依据压电元件的贴片位置和贴片范围,以及翼型上各点的切法向几何关系,创建“翼型-压电元件-冰层”的压电耦合分析模型。利用CATIA和Hypermesh组合创建有限元分析模型,ANSYS软件进行模型的结构动力学分析。翼型尺寸为展向200mm,翼型前缘弦长为80mm,材料为铝,材料参数为ρ=2780kg·m-3,E=7.05×1010Pa,μ=0.33;压电元件材料选用PZT-8;冰的材料参数为ρ=919.7kg·m-3,E=9.33×109Pa,μ=0.325。
对初始翼型进行模态分析,翼型固支方式为四边固支,得到翼型前八阶模态频率。由于设定目标除冰区域为翼型前缘中间区域,因此综合考虑各阶模态振型特点和压电元件可贴位置,选择第一阶模态振型为除冰振型,由此确定压电元件在翼型曲面上的贴片位置中心点坐标为(0.03215,-0.0149,0.1),该中心点为翼型振型的一个波峰极值点,压电元件的排布规律则依此中心点为基准来研究。
由于结构的模态是结构自身的一种属性,不会随着尺寸的缩放进行改变,并且布局规律研究是以翼型的模态振型为基准的,因此以按照尺寸缩放的机翼前缘为研究对象进行的压电原件布局规律研究,对于缩放前后的机翼而言,在规律上是通用的。所以研究是具有普遍性的。
3.1 压电元件在翼型展向上的排布规律为了保证振动和结构排布的对称性,上下翼面分别设置两个长条形压电元件,压电元件尺寸如图 7所示,计算了间距从10mm到190mm情况下19种排布形式的压电振动效果,对比了翼型前缘冰与翼型交界面处的位移情况,结果如图 8所示。
结果表明压电元件在翼型展向上的贴片对振动效果的影响,随着压电元件间距的增大,大致呈下降趋势,因此压电元件的布置在接近中心点区域时能够保证较好的激励效果。比较间距10mm和20mm下的振动效果,发现压电元件全部分布在振型波峰位置处不一定能产生最大的振动效果,因为压电元件附加刚度的影响,压电元件布置在波峰处会影响激励效果[18],所以排布在稍微偏离波峰位置处有更佳的激励效果。
3.2 压电元件相对贴片数量的影响文章提出了压电元件的相对贴片数量的定义:在相同的压电元件与基底结构接触面积下,单个压电元件分为若干个较小尺寸的压电元件,小尺寸压电元件的数量即为相对贴片数量。因此,选定初始单个压电元件的尺寸为20mm×20mm×1mm,即保证整个压电元件与基底结构的接触面积为400mm2,压电元件的厚度均为1mm。研究的相对贴片数量和规格如下表 1所示。
结果如图 9所示,该图表示冰与翼型交界面处展向上节点法向位移的分布情况。结果表明,在相同的接触面积下,单个压电元件具有最佳的激励效果,其次为单侧2片的布置方式。随着相对贴片数量的增加,压电元件的激励效果也逐渐减弱,由于压电元件激励点的分布从集中到分散,对于整个翼型曲面而言,较分散的激励作用对减弱对整体结构的激励效果。因此在相同的贴片接触面积下,适当减少压电元件的贴片数量会增大压电元件的激励效果。
3.3 压电元件贴片集中度的影响对于无法粘贴单个尺寸较大的压电片的曲面基底结构,只能采取多个稍小尺寸的压电元件的布置方式,因此3.2节的研究结果是有实际应用价值的。对于单个压电元件附加刚度的影响,若要把压电元件布置在振型波峰位置处,只能选择分为小尺寸的压电元件,为了确保较好的激励效果而又不分散布置压电元件,需要研究在固定贴片位置处的贴片集中度的影响,文章提出的贴片集中度的定义是:在固定的贴片位置附近,把单个大尺寸压电元件分为若干小的压电元件,但压电元件之间的间距很小,基本保证整体压电片贴片位置不发生明显变化,在此情况下分出的压电元件数量越多,认为贴片集中度越小。因此,仍采用3.2节中单个压电元件的规格,具体的分块情况如下图 10所示
不同集中度下的激励效果如图 11所示,结果表明贴片集中度越高,压电元件的激励效果越好。对比单侧均为2片的情况,纵向2片比横向2片有更好的激励效果,原因在于在翼型展向上,压电片集中度的影响较小,而在翼型截面曲线方向上,压电片集中度的影响较大,因此集中度较高的纵向压电元件比横向压电元件有更大的激励效果。因此,在固定贴片位置区域的情况下,贴片集中度越高,激励效果越好,同时翼型截面曲线方向上的贴片集中度会对整体结构激励效果有影响。
3.4 最佳排布方式设计及除冰效果模拟依据上述的压电元件布局规律,设计了两种压电元件布局方式,如下图 12所示。布局方式1中依据相对贴片数量的规律,在实验用压电元件尺寸为10mm×10mm×1mm和15mm×15mm×1mm两种规格的情况下选取了较大尺寸的压电元件,依据相对贴片数量设置了单侧4片的布局方式;布局方式2中依据压电片集中度和压电片在展向上的排布规律来设计,在振型波峰位置附近,依据压电元件贴片集中度规律,设置了稍小尺寸的压电元件,另外设置的两个大尺寸压电元件的布置位置在展向允许的布置位置范围内。
采用与姚塞金[17]等研究中相同的仿真手段,通过对整体模型进行结构动力学分析的方式,对两种布局方式进行除冰效果模拟,得到了两种情况下的冰与翼型交界面处的剪切应力,对于两种布局方式,其最大剪切应力值分别为0.634MPa和0.53MPa。相关学者文献中关于冰的剪切粘附强度的实验数据范围大致在0.14MPa~1.52MPa内,该数据范围是所有实验中冰层脱落时对应的粘附强度测量值的整体统计结果,本文的计算结果处于该范围内,因此在理论上具有冰脱落的可能性。另外,文章研究中加载在压电元件上的负载电压为100V,整体压电元件对应的负载功率较小,而在增大负载电压的过程中,相关剪切应力值也会增大,因此该研究可以满足除冰需要。
4 除冰实验依据布局方式1进行了实际地面除冰实验研究。实验包含了实验用机翼缩比翼型结构、信号发生器、功率放大器和示波器等,大致结构如图 13所示。对于实验中用到的机翼前缘压电元件的排布方式,其三维结构如图 14所示。
地面除冰实验在冰箱的冷冻环境中进行,在效果较好的实验中,压电片上施加的电压激励一般在320Vp-p,电流最大为0.06Ap-p。由文献[18]可得除冰功率为:
本文使用的机翼前缘结构尺寸为200mm×172mm×1mm,总面积为0.0344m2,则单位面积的功耗为69.77W·m-2。该功耗值同文献[17]和文献[18]的结果相比稍大,但小于一般的电热除冰系统的功率损耗,原因在于上述文献的除冰对象为平板,其尺寸较大,振动比较强烈。本文的除冰对象为实际的缩比翼型前缘结构,曲面结构的振动同平板相比较困难,而且使用的压电元件尺寸稍小,因此在除冰功耗上会有所增加。
5 结论(1)依据有限元模型导出的离散数据点,通过相关算法可以求解得到不同尺寸下压电元件在翼型曲面上的可贴区域,为实际压电元件的布局位置提供了参考。计算得到,边长为10mm的压电元件在上翼面实际的可贴压电元件范围为(0.0162,0.08),边长为15mm的压电元件在上翼面实际的可贴压电元件范围为(0.02696,0.08)。
(2)本文为研究翼型曲面的相关变量,通过五点共圆法和四点共圆法的组合,来求解翼型曲面上每个点对应的切法向矢量。通过对原始机翼的法向扩展机翼的计算,验证了切法向矢量求解的正确性,对翼型相关分量的求解提供了帮助。
(3)本文通过仿真的方法研究翼型上压电元件排布的相关规律。对于压电元件在翼型展向上的排布而言,随着两个对称压电元件间距的增大,压电元件对整体结构的激励效果呈下降趋势,而考虑压电元件附加刚度的影响,布局在振型波峰附近的压电元件可以产生最佳的激励效果,但在一定间距范围内进行压电元件的布置是可以满足实验要求的。对于压电元件的相对贴片数量而言,在接触面积不变的前提下,单个压电元件具有最强的激励效果,随着相对贴片数量的增加,压电元件对于整体结构的激励效果会减弱,因此适当减少相对贴片数量可以提高压电元件的激励效果。对于压电元件的贴片集中度而言,贴片集中度越高,激励效果越好,并且在翼型截面曲线方向上具有较高的贴片集中度,也可以具有很好的激励效果。
(4)依据压电元件的排布规律设计了两种布局方式,通过除冰数值模拟计算得到相应的最大剪切应力值分别为0.634MPa和0.53MPa,满足除冰需要。针对布局方式1进行了实际除冰实验,在不考虑电子设备的功耗等因素,除冰功率最大为69.77W·m-2。该功率同相关平板压电除冰实验相比稍大一些,但远小于电热除冰系统的功耗,因此对于翼型压电除冰方法而言,具有良好的研究价值。
[1] | 肖春华, 桂业伟, 杜雁霞. 电热除冰传热特性的结冰风洞实验研究[J]. 实验流体力学, 2010, 24(04): 21-24. Xiao C H, Gui Y W, Du Y X. Experimental study on heat transfer characteristics of aircraft electro-thermal deicing in icing wind tunnel[J]. Journal of Experiments in Fluid Mechanics, 2010, 24(04): 21-24. |
[2] | 易贤, 赵萍, 陈坤, 等. 水平轴风力机结冰探测器设计[J]. 空气动力学学报, 2013, 31(2): 260-265. Yi X, Zhao P, Chen K, et al. Method of designing icing prober for an horizontal axis wing turbine[J]. Acta Aerodynamica Sinica, 2013, 31(2): 260-265. |
[3] | Gent R W, Dart N P, Cansdale J T. Aircraft icing[J]. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2000, 358(1776): 2873-2911. |
[4] | Lynch F T, Khodadoust A. Effects of ice accretions on aircraft aerodynamics[J]. Progress in Aerospace Sciences, 2001, 37(8): 669-767. |
[5] | Cancilla D A, Martinez J, Van Aggelen G C. Detection of aircraft deicing/antiicing fluid additives in a perched water monitoring well at an international airport[J]. Environmental Science & Technology, 1998, 32(23): 3834-3835. |
[6] | Myers T G, Hammond D W. Ice and water film growth from incoming supercooled droplets[J]. International Journal of Heat and Mass Transfer, 1999, 42(12): 2233-2242. |
[7] | Haslim L A, Lee R D. Electro-expulsive separation system[P]. U.S. Patent 4, 690, 353, 1987-9-1. |
[8] | Reinmann J J, Shaw R J, Ranaudo R J. NASA's program on icing research and technology[R]. Cleveland: Lewis Research Center, 1989. |
[9] | Martin C A, Putt J C. Advanced pneumatic impulse ice protection system (PIIP) for aircraft[J]. Journal of Aircraft, 1992, 29(4): 714-716. |
[10] | Palacios J L, Smith E C. Dynamic analysis and experimental testing of thin-walled structures driven by shear tube actuators[D]. Pennsylvania State University, 2004. |
[11] | Palacios J L, Zhu Y, Smith E C, et al. Ultrasonic shear and lamb wave interface stress for helicopter rotor de-icing purposes[J]. 47th AIAA Structural Dynamics & Materials, Newport, Rhode Island, 2006. |
[12] | Overmeyer A, Palacios J L, Smith E C. Actuator bonding optimization and system control of a rotor blade ultrasonic deicing system[C]//53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference 20th AIAA/ASME/AHS Adaptive Structures Conference 14th AIAA. 2012: 1476. |
[13] | Palacios J L, Smith E C, Gao H, et al. Ultrasonic shear wave anti-icing system for helicopter rotor blades[C]//Annual Forum Proceedings-American Helicopter Society. American Helicopter Society, Inc, 2006, 62(3): 1492. |
[14] | Ven Katraman S B K K. Piezo-actuated vibratory deicing of a flat plate[J]. 2005. |
[15] | Venna S, Lin Y J, Botura G. Piezoelectric transducer actuated leading edge deicing with simultaneous shear and impulse forces[J]. Journal of Aircraft, 2007, 44(2): 509-515. |
[16] | 严春霞. 基于弹性波的机翼除冰机理研究[D]. 南京: 南京航空航天大学, 2009. Yan C X. Study on the mechanism of wing deicing method based on elastic wave[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2009. |
[17] | 姚赛金. 基于压电驱动器的机翼除冰方法研究[D]. 南京: 南京航空航天大学, 2010. Yao S J. Deicing method based on the piezoelectric actuators for aircraft wings[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2010. |
[18] | Bai T, Zhu C, Miao B, et al. Vibration deicing method with piezoelectric actuators[J]. Journal of Vibroengineering, 2015, 17(1): 61-73. |
[19] | Tookey R M, Ball A A. Estimation of curvatures from planar point data[J]. The Mathematics of Surfaces VII, Information Geometers, Winchester, 1997: 131-141. |
[20] | Ma X, Cripps R J. Estimation of end curvatures from planar point data[M]//Mathematics of Surfaces XII. Berlin Heidelberg: Springer, 2007: 307-319. |