构造法在机翼气动载荷转换中的应用

构造法在机翼气动载荷转换中的应用

作者：胡亮文梁勇

来源：《计算机辅助工程》2014年第06期

摘要：针对机翼气动压力分布数据的离散性，构造任意点附近区域3种形式的压力场曲面函数，以有限元单元为积分区域积分得到该单元内的集中气动载荷，并将其按照最小变形能原理分配到有限元的节点之上.结合某型飞机的测压试验数据比较3种构造函数法求得的气动累积内力，结果表明构造连续的压力曲面函数求得的气动累积内力符合实际气动力的分布特性，满足实际工程需求.

关键词：构造函数法；曲面函数；载荷转换；最小变形能；内力；测压试验；气动力

中图分类号： V211.412文献标志码： B

0引言

在飞机结构有限元分析中，通常需要将气动载荷转换为结构有限元的节点载荷.基于这种需求，出现各种各样的气动载荷转换算法.在20世纪90年代之前，载荷转换的算法按照“三点排”分布方案，但是局部受力分配变化很大，可能出现负值以及零值分配点很多的情况.[1]近年，国内外的专家[28]和学者相继研究开发出不同的载荷转换算法.王仁宏等[9]提出基于最小变形能的载荷分布计算算法，即通过建立极值函数，将气动节点力转换到有限元节点上.2005年，腾春明等[10]依据此算法基于Nastran开发机翼的有限元载荷加载模块.王专利[11]通过工程算例证明此算法的计算结果符合气动载荷分布规律.于哲峰等[12]提出基于距离加权插值法和薄壁样条法的载荷转换算法，在Patran中实现三维流体向有限元节点载荷的转换.戴愚志等[13]提出离散化思想的载荷分布算法；林小厦等[14]提出基于特征函数分布的曲面有限元加载方法；尹晶等[15]提出椭圆和抛物线2种分布形式下的近似解析方法.

本文从积分的角度出发，构造3种不同形式的压力场曲面函数，通过积分求得任意微元面积内的气动力及其作用点，然后按照最小变形能原理将积分得到的气动载荷分配到有限元节点上，并比较3种算法求得的累计气动剪力、弯矩和扭矩.

1压力分布函数的构造法

气动网格点通常以离散格式，即气动网格点的压力离散场的形式出现.以有限元网格为出发点，求出有限单元网格内的集中气动力，首先假定作用于翼面上的气动力是垂直于弦平面即机翼曲面的垂直投影，假设气动外形和有限元模型的外形基本一致.将气动外形内任意一点（x，y）和压力场函数设为Cp=f（x，y），构造3种压力分布函数.

1.1构造函数1

将有限元网格内的任意一点设为（x0，y0），首先找到该点附近不在同一直线上的3点，分别设为（x1，y1），（x2，y2），（x3，y3），构造压力曲面函数为Cp=ax+by+c （1）由过不在同一平面的3点可得出Cp1=ax1+by1+c

Cp2=ax2+by2+c

Cp3=ax3+by3+c （2）求出该平面的系数a，b和c，将（x0，y0）代入式（1）即可求出任意一点的压力值.

1.2构造函数2

将有限元网格内的任意一点设为（x0，y0），找到该点附近不在同一直线上的4个点设为（x1，y1），（x2，y2），（x3，y3），（x4，y4），构造4点的曲面方程为

Cp=axy+bx+cy+d （3）由于该曲面过上述4个点，分别将4个点的坐标代入即可求出曲面方程的系数a，b，c和d.

1.3构造函数3

首先找到不在同一平面的4个点设为A，B，C和D，4个点的坐标依次为（x1，y1），（x2，y2），（x3，y3）和（x4，y4）.该4点构成2个连续且光滑的曲面ABC和曲面ACD，见图1.

图 1压力分布曲面

Fig.1Curved surface of press distribution

构造压力曲面函数为Cp=a1x2y+a2xy2+a3xy+a4x+a5y （4）设ABC构成的曲面为

Cp=a11x2y+a12xy2+a13xy+a14x+a15y （5）ACD构成的曲面为

Cp=a21x2y+a22xy2+a23xy+a24x+a25y （6）将4个点的坐标代入得

Cp1=a11x21y1+a12x1y21+a13x1y1+a14x1+a15y1

Cp2=a11x22y2+a12x2y22+a13x2y2+a14x2+a15y2

Cp3=a11x23y3+a12x3y23+a13x3y3+a14x3+a15y3

Cp1=a21x21y1+a22x1y21+a23x1y1+a24x1+a25y1

Cp3=a21x23y3+a22x3y23+a23x3y3+a24x3+a25y3

Cp4=a21x24y4+a22x4y24+a23x4y4+a24x4+a25y4 （7）由于真实气动压力场函数在任何一点均满足1阶连续性，A点和C点均在两曲面之上，则A点和C点满足Cp1x=Cp3x，

Cp1y=Cp3y （8）即

2a11x1y1+a12y21+a13y1+a14=

2a21x1y1+a22y21+a23y1+a24

a11x21+2a12x1y1+a13x1+a15=

a21x21+2a22x1y1+a23x1+a25

2a11x3y3+a12y23+a13y3+a14=

2a21x3y3+a22y23+a23y3+a24

a11x23+2a12x3y3+a13x3+a15=

a21x23+2a22x3y3+a23x3+a25 （9）

联立上述9个方程即可求得函数关系式系数aij（1≤i≤2，1≤j≤5）.判断（x0，y0）在曲面轮廓的区域位置：若该点在ABC的平面投影内则代入曲面ABC函数；若在ACD平面投影内则代入曲面ACD函数.对于特殊点不在曲面ABCD平面的投影内的，分别代入2个曲面函数后取平均值.

1.4气动力的求解和转换

在求出任意一点的压力函数库后，任意面元ds内的气动力F以及作用点（xf，yf）为

F=Cpidxdy

xf=xCpidxdyF

yf=yCpidxdyF（10）积分法示意见图2.

图 2积分法示意

Fig.2Schematic of integral method

由积分得到气动力和作用点的坐标，根据最小势能原理[9]可将某一单元的集中气动力转换到该单元的节点上.

2算例

根据某型机翼的测压试验数据，取一展长为18 m，弦长为2 m的矩形机翼为考核算例.气动网格在展向共划分为20个气动剖面，弦向共设置31个测压点，设沿展向为y轴、弦向为x 轴作为分析坐标系.某典型工况的上、下翼面三维压力分布见图3和4.