弹性三角翼的静气动弹性方法研究

弹性三角翼的静气动弹性方法研究 弹性三角翼的静气动弹性方法研究 第24卷第4期 2006年l2月 气动研究与实验 AERODYNAMICRESEARCH&EXPERIMENT Vo1.24,No.4 Dec.,2006 弹性三角翼的静气动弹性方法研究 郭承鹏董军杨庆华李俊甫张铁军 (沈阳110034) 摘要静气动弹性问题的研究关键是要解决流固数据耦舍的问题,即把气动力从 气动网格点插值到结构点,由柔度影响系数法求得结构变形后,再把结构变形插值到气 动网格上.本文在O—H网格的基础上,采用Jameson的中心差分有限体积法求解Euler/ N—S方程,采用结构影响系数法计算结构的弹性变形,用三角元面积加权法和常体积转 换法(CVT)实现流固力的耦合,采用无限平板样条法(IPS)实现位移的耦合. 关键词静气动弹性流固耦合三角元面积加权法常体积转换法无限平板 样条法动网格技术 引言 为了满足现代飞机设计的要求,如第四代战斗机超音速巡航,无人机的高空长航时, 运输机减轻结构重量等,都使得结构弹性变形对飞机飞行性能的影响越来越严重.而过 去的以刚性飞机为设计对象,并加以弹性修正的方法已经过时.需要以真实的弹性机体 作为设计对象. 静气动弹性问题的研究关键是要解决流固数据耦合的问题,即把气动力从气动网格 点插值到结构点,由柔度影响系数法求得结构变形后,再把结构变形插值到气动网格 一 t-.本文在O—H网格的基础上,采用Jameson的中心差分有限体积法…求解Euler/N— s方程,采用结构影响系数法…计算结构的弹性变形,用三角元面积加权法…和常体积转 换法(CV'I')…实现流固力的耦合,采用无限平板样条法(IPS)HI实现位移的耦合. 1数值方法,边界条件和计算网格 1.1控制方程 本文采用的控制方程为直角坐标系下的可压缩流Euler/Navier—Stokes方程,其积分 形式为: 旦 Ot+一惫? 其中: 第4期弹性三角翼的静气动弹性方法研究9 O=,pu,pv,pw,pE)r I= I= pu +P puv puw (pE+p)u l00 r"rJ v tlt啦 r工rr:v 0 rn r r:: —— R为控制体,S为控制体面,dS为S控制体表面微元的外法向面积向量,that,为 正.若L=0则方程(1)退化为Euler方程. 1.2空1"oq离散 流场划分网格,假定流体守恒量Q定义在每个格子中心,且在每个格子中是常数.将 积分方程用于每个小格子,可得如下离散方程: 砌k? 为网格单元的体积,为流体流出网格单元的无粘净通量,占f为流体流出 网格单元的有粘净通量. 在求解Navier—Stokes方程组时,无粘项和粘性项的离散方法是不同的.粘性通量的 计算通常采用二阶的中心差分格式;关于无粘净通量的汁算,不同的计算方法可以引出 不同的计算格式,本文中采用的是二阶的中心差分格式. 1.3边界条件 对于Euler方程物面上采用无穿透速度边界条件,而对于N—S方程物面上采用无滑 移速度边界条件.Euler方程和N—S方程都采川绝热壁边界条件,对称面边界条件和远 场无反射边界条件. 2气动一结构耦合数据交换方法 2.1力的耦合一三角元面积加权法 由于气动计算的网格和结构计算的网格不同,一般气动网格较密,结构的网格较稀,需 要将气动网格上的压力信息传递到结构网格点上,将变形后的结构网格点的位移信息传递 到气动网格点上.本文采用三角元面积加权法.对每个气动点(此处为网格节点的体心点), 帅 , 10弹性三角翼的静气动弹性方法研究第4期 寻找它附近最近的三个结构点点1,点2,点3(空心点),具体寻找方法参照了文献,并做 这个气动点在三个结构点构成平面内的投影点4(实心点),这样由三结构点和一个气 动点构成了4个三角形,其面积分别为S.,S:,S,S,其中面积的计算 Y2 3 4 lYzl : I,zI,Il lX1YlzlI .S4如果 气动投影点落在结构三角形的内部,则S.,S,S,S都为正,反之某些值则会出现负值, 如图1(b)4点落在三角形外部靠23边处,则S(0,则气动点4上的离散压力FiJ分解到 三个结构点上,其离散压力分别为/1)=F-/S,/2)=F$2/s,/3)=fi3/&, 式中S./S,5/S,53/S就是所说的权值. (b) 图1三角元面积加权法示意图 2.2位移的耦合 在气动弹性的问题中,需要把结构在气动力载荷的作用下得到的变形传递到气动网 格点上,以实现气动计算的需要. 2.2.1常体积转换法 常体积转换法(CVT),首先是由Gour~q[I等人提出的,该方法是一种与结构模态无关 的局部插值方法,其基本原理是,由三个结构点和一个气动网格点构成了一个空间四面 体,求出这时四面体的体积和气动点在结构点平面的投影点相对位置(由三角元面积权 值来确定),当结构点发生位移,到达一个新的位置时,这时的气动点和结构点构成了一 ZZZ ; = = = 绎?期弹性三角翼的静气动弹性方法研究11 个新的四面体,这时的四面体体积和变形前四面体体积相等,并且满足气动点在结构点 平面投影相对位置相同.该方法的优点是它只与局部信息有关(也就足说,不需要计算和 储存矩阵),计算量比较小. 2 图2CVT方法原理示意图 (b) 用(Y?,?),(y2,:z),(y3,)来表示变形前结构三角元I,2,3的坐标,(毛 y,)表示变形前气动点在123平面的投影坐标,(,),4,)表示变形前气动点坐标.由上 文中的三角元面积加权法求出投影点在123平面内的权值,,,则: 其中++=1,123平面的法向量(1,m,n)为: } Y2——Yl y一yl兰f=,+,巧+气动网格点及其投影点的关系为: +v, J,-4-啪 z-4-1-' (4) 由气动网格点与三角单元组成的四面体的有向体积为: 3 .. t , 2-+m2+nF). f.1(x4-x,)+—m(y~-y,)—+n(z-z,)?,++2 = 吉)+)+(z.(6) ++卜 ++} {{== Z ,?????J,????,??L 1.一一 12弹性三翼的静气动弹性方法研究第4期 变形后结构三角元1,2,3的坐标,气动网格点及其投影点的坐标分0为【.,,Y"' -),(,Y',),(',y,:),(,',),(',Y',:'),变形前后气动投影点在123内的 相对位置没有改变,就是说仪值,i3,没有改变,则变形后气动网格点标勾: I:=Xt+V','=Q+++V',{:=-)If.}?川?:+十刀'.vtt'7' 【z:.z+'=02;十十;+V' 变形后四体体积为: er = 言(,'(',t)4-m(.~]t4--t)(8) ?'【z'4,z'l 将(7)式代人(8)式中,并H.有变形前后四面体体积相等.=,可以得出未知数 V'.将解出的vlf?入(7)式中可街变肜后气动网格点坐标(:山Y,:).此就完成了结构 点位移往气动网格点位移的粥台.. 2.2.2置跟平板样条法 我们假定…"1--17341F;.IH/,j分布载荷下产生的垂直J扳囊仃玎的位移满足下式: (jf,)=,!/,J!).:,(:1,2…,)(9) 其中,=(.一)!+(Y一YJ),这样构成了有N+3个未知数的N个方程细成的方程组, 还缺少三个J,由下式补充: ?=?-=?=0?i=lJ!i=I 把它写成阵形式: ,,J, 1 In2 】 , In I /川 In rl}1rN N{N /.2 lN J.I 2 3 y 解此线性方稗组可得出多项式的系数-,…F,‰n.若要求N个结构点之 外的任意点的垂直于平板表面的位移则可用下式求得: ,):?N厂zln(r)+..+口.+口其中为N个结构点之外的任意点的坐标, 一 ' 1I1??,l,l ,|..鸵,,,,,;J,,r,,r,,,,frf22,^2r? ,0OO ,OO0 ,OO0 一 , . ,,,,, ? 一 ,l-? 一,? i『.三『一O00 第4期弹性三角冀的静气动弹性方法研究13 ?f,),J为该点在垂直于平板表面的位移,r2=(一.)+(Y—y). 由于无限平板样条法式一种二维插值方法,它不能够求解非给定的位移方向的其他 方向的位移.IPS和CVq,插值方法的区别见下图,可见IPS是二维插值方法,CVT是三维 插值方法.然而cV.r方法插值结果网格不光滑,而IPS插值方法得到的网格比较光滑,在 计算姬t-行,1tS方法算量也比较小,由于变形量是个小量,三维办法CVT干?二维方法 IPS方法区别也不明显.故采用tPS方法较之CVT方法要好. 图3三维方法和二维方法的区别 图4IPS方法插值出的表面网格图5CVrr方法插值出的表面网格3 3动网格技术 本算例采用的网格是0一H型的结构网格,是一个沿着流向的二维网格的组合,因此 动网格技术也可以采用准二维动网格技术.其基本原理是,找这样一个函数,满足在{边 界上是1,在外边界t足0,函数值从内边界往外边界上逐渐变小,这样越靠近内边界的 气动网格变形越接近刚性变形,越靠近外边界变形越小.余弦函数"就满足这样的条件. 取口=cos()for:1, U^一l 则当=1(内边界)时=1,当J.=(外边界)时cd=0,由于可以将计算网格当作准 二维网格,而且气动弹性变形仅仅在法向变形量值又很小,故气动网格点法向坐标变化 由下式确定,另外两个方向的坐标不变: Znew . =z+j" "j='其[}1"=z一z篙,表示内边界弹性变形量,为了控制网 l4弹性三角翼的静气动弹性方法研究第4期 格质量,避免网格重叠,引入c=-(血'.一血,),?I1,2I 如果zurig+lIid^一z<,Ij~iJ二:1. 由此在离内边界较近的地方很密的网格做近乎刚性的移动,保证了网格在内边界附 近不重叠. 4气动弹性迭代过程 由于机翼是弹性的,在气动力的作用下将发生弹性变形,机翼的弹性变形将造成机 翼气动状态的改变,新的机翼状态将产生新的气动力,新的气动力又将使机翼发生新的 变形,如此往返循环.所以,求解弹性机翼平衡时的变形和载荷就将是一个迭代过程.此 静气动弹性迭代过程实现如下: 第一步:把机翼作为刚性的,求解Euler/N—S方程得到机翼承受的气动力载荷; 第二步:将第一步中的气动力载荷耦合到结构点上,并由结构影响系数法得到弹性 机翼变形后的新位置; 第三步:判断机翼是否达到弹性平衡位置; 是:迭代结束,得到机冀的弹性平衡位置和气动力载荷; 否:返回第一步,继续进行迭代. 第四步:将机翼结构点的变形量耦合到气动网格点上,并生成新的气动网格. 5计算结果与分析 本文采用的计算模型是45.后掠角的菱形三角翼,计算网格为HO型网格,是一个顺 流向的准二维网格,数量为61×51×41,其中顺流向网格数为5l,在三角翼面网格数为 l9,法向网格数为51,绕机翼周向网格数为4l.采用Jameson的有限体积法求解Euler控 制方程,用上述流固耦合方法及动网格技术等完成静气动弹性问题的计算.图6为 本文 图6三角翼计算网格 ?坩蛹心档?期曲霄帅} 图7弹性:三角翼曲线 撇,懈眦姗懈一 第4期弹性三角翼的静气动弹性方法研究15 计算的弹性三角翼c一q曲线和实验数据以及文献中计算结果的比较曲线. 由图7的曲线可知本文算得的在同一速压q一F比实验数据大,同文献6中算的 相比,除去速压很小和很大时稍有偏差外,在实验速压范围内吻合很好,且计算数据 连线几乎同实验数据连线平行.由此,验证了本文所述的流固耦合方法是_?J行的,静气 动弹性程序是正确的,可以为研究翼身组合体,甚至是全机静气动弹性问题做好方法,程 序铺垫. 参考文献 [1JJamesonN.SchmitW.TurkelE.NumericalsolutionoftheEulerequationsby finitevolumemethodsusingRunge—KuttalimesteppingschemeIRI.AIAAPaper1981一 I259.198I. I2lManojK.Bhardwaj.ACFD/CSDINTERACTIONMETHODOLOGYFORAIR. CRAFTWINCS. I3lG..S.L.Goura,K.J.Badcoeketc.Adataexchangemethodforfluid— structureinter—actionproblemsIJI,JournalofAircraft,Vo1.38,No.2,March—April 2001.PI282,342. I4I Transferin I5I l6l M.Sadeghi,F.Liu.ApplicationofFhree—DimensionalInterfacesforData Aero—elasticComputations.AIAA2004—5376. 李俊甫.弹性三角翼结构柔度矩阵原始数据. 李俊甫.简化全机静弹性气动力计算报告.2001年.