西北工业大学学报基于欧拉方程的种机翼气动弹性计算方法叶正寅王刚杨永年杨炳渊21.西北工业大学飞机系，陕西西安710072；2，上海航天局，上海200233格，考虑了机翼变形时的网格生成问，通过与气动力方程的联立求解，在时间域内用阶龙格库塔方法求解机翼弹性运动方程。计算结果明，本文计算方法具有较高的计算效率，所计算的颤振临界速度与风洞实验致。
随着计算机和计算流体力学的发展，通过求解欧拉方程计算气动弹性问已成为可能，但由于国内计算机能力的限制，国内此方面的工作落后国外较多。本文研究种效功能强大的气动弹性计算方法。
1网格生成无限插值理论首先是根据求解欧拉方程的需要而发展起来的，是种较好的代数方法，其原理是从内边界即机翼面按定要求如内边界处网格与物面正交网格在物面处密离物面越远网格越稀等等，使网格逐渐光滑地连接到外边界，其公式可简单写为3尸3切2物体面对应点处法向导数值；2外边界上对应点坐标；此为收稿日期20000912基金项目航空科学基金0053001从机翼面到外边界之间的个变化参数，变化范围，切，为人工设计的融合函数，具体文献1.
2欧拉方程的求解所用的气动力基本方程为可压缩流非定常欧拉方程，其积分形式为户，分别为空气密度工，2方向的速度分量和单位体积的总内能为面积分的法向单位向量少为体积分域为包含体积分域的边界为通量项。
用中心格式的有限体积法，欧拉方程写成离散形式后，可为物理量平均值为该网格单元的通量项。为避免奇偶单元的耦合和激波附近的振荡，需要加上人工粘性项，即这里为人工粘性项，其构成取自文献2的方法。
与定常问求解不同，对于非定常流动，4式中的时间有物理意义，为了仍采用原来的显式龙格库塔推进方法，将4式的时间项用阶精度假设伪时间5式又可写为在伪时间名义下，6式与原来的定常龙格库塔方法求解定常问具有相同的形式，只是残值为5式中的形式，将6式推进到收敛满足某收敛指标后，便满足了6式中非定常问的离散方程。本文采用步龙格库塔方法对6式作伪时间的推进，具体过程参文献3.
3结构运动方程的求解以采用传统振动模态的概念来描述机翼的变形；阶振型，9，为与第阶振型对应的广义坐标。
应用拉格朗日方程，机翼的运动方程可以写为矩阵形式广义气动力。
为了便于用龙格库塔方法求解方程8，引入状态变量，即则，方程8变成方程9可以用阶精度的龙格库塔方法作时间推进，在每推进个时间步长时，由欧拉方程计算的压力分布提供所需的广义气动力矩阵。
对于有迎角的情况，机翼除动态变形外，定常气动力也会在弹性机翼上产生变形，这种静变形在线性结构的情况下对结构振动结果没有影响，但由于气动力已经是非线性问，其静变形对空气动力的影响是存在的。为了计入静变形的影响，只需将方程8的动态项去掉，即可计算出静变形计算气动弹性问中，其具体实现步骤为让机翼以定迎角此时迎角为基本迎角以希望的从数从静止状态启动，当绕流流场达到稳定后，便得到此迎角时的定常气动力，在此阶段，让弹性结构不产生变形。如果已有此迎角定常流场的话，也可以直接读人此已知流场。由于本文求解程序。
在得到基本迎角的定常流场后，计算翼面上的气动力载荷，由弹性结构与气动力的方程得到结构的静变形量，但此时的结构变形量不是此基本迎角的静变形，因为有了此变形量后，结构外形的变化会带来流场和气动力的变形，那么在基本迎角的基础上加上此静变形量，再求解空气动力方程计算新的流场。有了新的流场气动力又可计算出新轮同达到定指标，便认为此变形量为此基本迎角得到基本迎角的静变形和相应的流场后，此时的结构弹性力与气动力之间是平衡的，如加上定的初始扰动量，弹性结构将进人个动态响应过程。由于本项目采用双时间求解欧拉方程的方。法，在真实的物理时间域内的每时间步长中，由基本迎角加上此时的结构变形量产生相应的网格，求解欧拉方程，得到此时的气动力载荷，又由此气动力载荷计算出结构运动方程中的广义气动力项，通过级龙格库塔方法，求出新的结构变形量。此变形量将作为下时间步长求解空气动力方程的物面边界条件，反复如此进行时间推进，便得到此财数基本迎角和初始扰动量情况下的动态响应过程，由此响应过程的特征可判断出气动弹性的特性。通过改变从数基本迎角和初始扰动量，便可得到所有的气动弹性信息。
4算例为验证方法，第l个算例为以NACA0012翼型为剖面的直机翼，目的是检验非定常气动力计算的精度和收敛性，其机翼绕14弦作俯仰简谐振动，迎角的变化规律为≈=am+aA，sinco，t式中，为基本迎角，心为迎角变化的幅值为振动的频率，如频率变换成无因次形式，则各参数的计算的翼根处剖面升力系数和绕0.273倍弦线矩系数对迎角的变化曲线与有关文献3是致的。
在1中，给出了两种时间步长的计箅结果，种是每周期64步，另种是每周期只有8步近似8步。从计算结果看，即使是每周期只有8步，计算结果与64步的结果也是接近的，而且经历了十几个周期后，其结果仍没有偏移现象还是在它应该所处的曲线上，当然，由于个周期中只有8个点，且8个点并不是正好重复到原迎角位置，所以曲线看上去，似乎不同周期没有重复，可本文的方法具有很好00升力随迎角变化1升力随迎角变化放大幻力矩随迎角变化角抬头瞬时的剖面压力分布，可即使只有8个步长个周期的近似结果，压力分布与64个步长个周期的压力分布吻合很好，明本方法具有很好的精度保证。
第2个算例为气动弹性算例，3给出了机翼阵为两阶振型的自然振动频率分别为218.340和度的动态响应过程，从的结果看，速度为9，1 5基本上处在颤振临界状态而且从4和41定的，它不依赖于初始扰动的大小，用偶极子网格法计算的颤振临界速度为7mS，颤振风洞实验洞实验致。
3速度为40，8时的响应结果扰动速度为5，0.51速度为68时的响应结果扰动速度为1.，1.，速度为9，1！5时的响应结果扰动速度为5，0.51速度为908时的应结果扰动速度为1.，1.，e速度为100ms时的响应结果扰动速度为1.，1.，5结论它能保证在真实物理时间的推进过程中流场计算的收敛精度达到使用者所期望的精度。
它是种全隐式方法，物理时间步长可以按使用者的愿望足够大，从而保证了方法具有效重要。
此计算方案具有显式方法的特点，降低了对计箅机的要求，编程也方便，对于能向量运箅的计算机，这种方法具有明显的优势。从计算的箅例情况看，本报告中提出的计算方案具有良好的精度保证和较的计算效率。
叶正寅，赵令诚，杨永年。大迎角机翼的气动弹性计算方法。航空学报，1990，1336
QQ交流群
