技术领域本发明涉及一种适用于CFD(计算流体力学)数值模拟的周期非定常流场的预测方法，该方法能够高效高精度地获得对飞行器周期非定常流场分布，并根据获得的流场分布能够得到飞行器的周期非定常气动力，然后根据获得的周期非定常气动力能够指导飞行器的动稳定性设计。
背景技术：
周期性非定常流动具有广泛的工程应用背景，比如直升机的旋翼绕流，包含定子与转子的涡轮机的内流流动，风机叶片的绕流，以及扑翼动力装置的研究。此外数值强迫振动法计算动导数，让物体作给定频率的强迫振动，关键也是获得周期运动条件下与振动频率直接相关的迟滞气动力。上述所有流动都有一个明显的特征频率，宏观上这类流动的周期性都是物体周期运动的结果。目前，通常采用双时间步法数值模拟这种非定常流动，但双时间步法并未考虑流动的周期特性，因此计算效率不高，尤其对于高超声速流动由于特征时间急剧减少，计算成本会更高。为减少计算量，一些学者提出了一系列的简化算法，然而非定常流动计算的关键是在计算效率和计算精度之间取得较好的平衡，高精度高分辨率的方法通常都被当前计算机的计算能力所限制，而低精度算法通常又不能反应一些重要的物理信息。因此，开发一种能快速计算周期性非定常流动的算法，要综合考虑其效率和精度两方面的要求。
技术实现要素：
本发明的技术解决问题是：克服现有技术的不足，提供一种适用于CFD数值模拟的周期非定常流场的预测方法。本发明的技术解决方案是：一种适用于CFD数值模拟的周期非定常流场的预测方法，步骤如下：(1)根据需要在计算域内生成计算网格；(2)将一个计算周期分成等时间间距的N个采样点，此时，N-S控制方程中第n个采样点的守恒变量Wn对时间的导数项离散为：∂∂tWn=Σj=0N-1D[n,j]Wj]]>其中n＝0,1,2,...,N-1，j＝0,1,2,...,N-1，N为自然数，D[n,j]表示维度为N×N矩阵D的第n行、第j列的元素；(3)根据步骤(1)获得的计算网格采用刚性动网格法获得各采样点的计算网格；(4)在流场解算器中设置来流初始条件，将各个采样点设置步骤(3)中相应时刻的计算网格，将流场解算器中的时间导数项改写为步骤(2)中的形式，计算得到N个采样时刻收敛流场与气动力；(5)将步骤(4)中得到的N个采样时刻的气动力，得到傅里叶级数展开系数：C^k≈1NΣn=0N-1Cne-ikωnΔt]]>其中Δt＝T/N，ω＝2π/T，T为一个周期的时间，k为谐波数，k＝0,+1,-1,+2,-2,...，Cn为第n个采样点的气动力，为傅里叶级数第k阶系数，重构周期内时刻t时的气动力Ct为：Ct≈Σk=-(N-1)/2(N-1)/2C^keikωt.]]>步骤(2)中，Wn＝(ρ,ρu,ρv,ρw,ρE)nT，其中ρ为流体密度，(u,v,w)为直角坐标系下的速度分量，E为单位质量气体的总能量。步骤(2)中，D[n,j]的表达式为：N为奇数N为偶数。步骤(3)中刚性动网格法为对初始网格进行整体旋转或平移。步骤(4)中来流初始条件包括雷诺数Re、马赫数Ma及来流温度T∞。步骤(4)中的流场解算器为三维可压缩雷诺时均Navier-Stokes(RANS)控制方程求解器，解算器中使用SA、MenterSST两种湍流模型模拟湍流流动。有益效果(1)本发明的时间谱方法将时域内随时间变化的周期非定常问题变换为N个耦合的定常方程，只需耦合计算N个定常时刻，即可重构整个周期，从而较大幅的提高计算效率。对高超声速流动，本发明的时间谱方法在不降低计算精度的条件下，计算效率能提高一个量级以上；(2)本发明的方法中刚性动网格法获得各采样点的计算网格，无需改变网格拓扑结构。(3)本发明的时间谱方法对现有解算器有较好的继承性，其中解算器中通量计算模块无需修改，只需要修改时间导数项。(4)本发明的时间谱方法对周期非定常混合雷诺时均Navier-Stokes方程(RANS)方程中湍流模型方程的时间导数项的离散采用统一的格式。附图说明图1为实施例1中的初始计算网格；图2为初始0.00s时刻的网格旋转得到0.02s时刻的计算网格；图3为初始0.00s时刻的网格旋转得到0.04s时刻的计算网格；图4为初始0.00s时刻的网格旋转得到0.06s时刻的计算网格；图5为初始0.00s时刻的网格旋转得到0.08s时刻的计算网格；图6为NACA0015翼型周期的俯仰力矩系数随攻角变化的迟滞曲线，其中，横坐标为攻角，纵坐标为俯仰力矩系数，TSM为采用本发明的方法当N＝5时得到的迟滞曲线，Dual-Time-Step为传统的双时间步法得到的迟滞曲线，Exp为迟滞曲线的试验结果曲线；图7为本发明的方法与传统的双时间步法CPU计算时间的比较；其中，TSM为采用本发明的方法CPU计算的时间，Dual-Time-Step为传统的双时间步法CPU计算的时间。具体实施方式建立计算非定常流动的有效方法，应该充分考虑其流动特征。周期性非定常流动其特征为经过若干个时间间隔后，流场的流动特性会重复出现。传统上，这种周期非定常与一般非定常问题同样对待。过去开发的算法大部分都是采用时间推进。这是由于在时间方向，任何时刻的解都只影响未来时刻的解。这种时间方向的抛物化处理方法，对大多数非定常流动都是可以接受的，比如飞行器收到脉冲扰动和作机动飞行，这时确定的初始条件必须准确，过渡过程也很重要。然而在周期性非定常问题中，这种时间方向只向后影响的假设不再完全合理，因为若只考虑一个收敛的周期，而不考虑不同周期的相互影响，这时可以认为周期内任何时刻的解都影响这个周期内其他时刻的解。对这类流动，通常由初始条件到稳定状态解之间过渡过程是非物理的，这种过渡过程很有可能收敛的非常慢，而大量的CPU资源都耗费于此。对这种周期非定常问题，时间推进方法都太过于昂贵，显然开发一种计算方法直接求解周期状态，而避免求解非物理的过渡过程，将较大的提高计算效率。对时间谱方法，当前所有通量计算都在时域内进行，对现有解算器有较好的继承性。1D矩阵计算采用傅里叶变化，n时刻的时间导数可以表示为D矩阵与周期内采样序列的乘积。2采样点网格计算要耦合计算周期内各采样点的流场，先要获得各采样点的计算网格。采用刚性动网格法获得各采样点的计算网格。采样点N＝5时，各采样点的计算网格示意图如图所示。3初始化流场对于流场的初始化有两种方式:(1)各采样点的初始数据都给定为来流条件；(2)计算一个初始定常流场，各采样点的初始数据都给定为该初始定常流场。经比较研究，两种初始化方式对收敛性的影响不大。4通量计算无粘通量和粘性通量计算都可继承已有解算器，但通量计算中要加入包含D矩阵的时间导数。5时间推进采用定常的虚拟时间步方法进行时间推进。采用全隐式的时间谱方法进行时间推进，可改善收敛稳定性。第4步至第5步可重复多次，直到获得满意的结果为止。6重构周期流场对计算得到的各采样时刻的流场和气动力，利用傅里叶反变换可整个周期流场。一种高效周期非定常流动预测方法，采用傅里叶变换将时域内随时间变化的周期非定常问题变换为N个耦合的定常方程，只需耦合计算N个定常时刻，即可重构整个周期，从而较大幅的提高计算效率。对SA湍流模型、MenterSST湍流模型方程也采用了时间谱方法离散。采用全隐式的时间谱方法进行时间推进。采用MPI框架，实现了并行计算，提高了计算规模。本发明介绍了一种周期非定常流动的高效预测方法。采用傅里叶变换将时域内随时间变化的周期非定常问题变换为N个耦合的定常方程，只需耦合计算N个定常时刻，即可重构整个周期，在不降低计算精度的条件下，可较大幅的提高计算效率。本发明所介绍的方法与已有的方法相比，可应用于周期非定常雷诺时均Navier-Stokes方程(RANS)方程，对SA(Spalart-Allmaras)湍流模型、MenterSST(Shear-StressTransport)湍流模型也采用了时间谱方法离散，提高了数值模拟精度和复杂流动适应性。所有通量计算都在时域内进行，对现有解算器有较好的继承性。采用时间谱方法的全隐式格式，改善了随着采样点数增大而出现的稳定性问题。可实现并行计算，适用于大规模的工程应用问题。一种适用于CFD数值模拟的周期非定常流场的预测方法，步骤如下：(1)根据需要在计算域内生成计算网格；(2)将一个计算周期分成等时间间距的N个采样点，此时，N-S控制方程中守恒变量Wn对时间的导数项离散为：(其中N为自然数，n＝0，1，2…N-1，j＝0，1，2…N-1)∂∂tWn=Σj=0N-1D[n,j]Wj]]>其中D[n,j]为维度为N×N矩阵D的第n行，第j列的元素，D[n,j]的表达式为：N为奇数N为偶数W＝(ρ,ρu,ρv,ρw,ρE)T为守恒变量，其中ρ为流体密度，(u,v,w)为直角坐标系下的速度分量，E为单位质量气体的总能量；则Wn＝(ρ,ρu,ρv,ρw,ρE)nT表示第n个采样点的守恒变量；(3)根据步骤(1)获得的计算网格采用刚性动网格法获得各采样点的计算网格；(4)在流场解算器中设置来流初始条件，包括雷诺数Re，马赫数Ma，及来流温度T∞，将各个采样点设置步骤(3)中相应时刻的计算网格，将流场解算器中的时间导数项改写为步骤(2)中形式，计算得到收敛流场与气动力；(5)将步骤(4)中得到的N个采样时刻的气动力，得到傅里叶级数展开系数：C^k≈1NΣn=0N-1Cne-ikωnΔt]]>其中Δt＝T/N，ω＝2π/T，T为一个周期的时间，k为谐波数，k＝0，+1，-1，+2，-2…；Cn为第n个采样点的气动力，为傅里叶级数第k阶系数，最终，可重构周期内任意时刻的气动力：Ct≈Σk=-(N-1)/2(N-1)/2C^keikωt.]]>实施例(1)计算模型采用NACA0015翼型，翼型弦长为c＝0.3048m，在十倍弦长范围内即计算域内生成计算网格，如图1所示，翼型表面有260个点，壁面第一层网格间距为3×10-6m，NACA0015翼型作正弦俯仰振动，即α＝4.0+4.2sin(20πt)其中，α为NACA0015翼型的攻角，t为振动时间；(2)将一个计算周期分成等时间间距的5个采样点，分别为(0.00s，0.02s,0.04s,0.06s,0.08s)，此时，D矩阵为：D=053.45-33.0333.03-53.45-53.45053.45-33.0333.0333.03-53.45053.45-33.03-33.0333.03-53.45053.4553.45-33.0333.03-53.450]]>(3)根据初始0.00s时刻的网格采用刚性动网格法旋转分别得到0.02s,0.04s,0.06s,0.08s时刻的计算网格，依次如图2-5所示；(4)在流场解算器中设置来流初始条件，其中雷诺数Re＝1.95×106，马赫数Ma＝0.29，及来流温度T∞＝288.16K；(5)将流场解算器中的时间导数项采用步骤(2)中D矩阵进行离散，即∂W0/∂t∂W1/∂t∂W2/∂t∂W3/∂t∂W4/∂t=053.45-33.0333.03-53.45-53.45053.45-33.0333.0333.03-53.45053.45-33.03-33.0333.03-53.45053.4553.45-33.0333.03-53.450W0W1W2W3W4]]>(6)流场解算器根据步骤(3)得到的五个采样时刻的计算网格、步骤(4)设置的来流初始条件以及步骤(5)离散后的时间导数项计算得到各采样点的俯仰力矩系数，为Cmz0Cmz1Cmz2Cmz3Cmz4=6.096E-04-1.192E-02-1.919E-02-9.032E-032.250E-03]]>(7)将步骤(6)中得到的N个采样时刻的俯仰力矩系数，得到傅里叶级数展开系数：C^k≈1NΣn=0N-1Cne-ikωnΔt]]>C^-2C^-1C^0C^1C^2=-5.77E-05+i2.66E-044.09E-03-i3.89E-03-7.464.09E-03+i3.89E-03-5.77E-05-i2.66E-04]]>(8)根据步骤(7)的结果并结合公式(2)得到周期内任意时刻的俯仰力矩系数Ct，并结合攻角随时间的变化规律公式(3)，得到NACA0015翼型周期的俯仰力矩系数随攻角变化的迟滞曲线，如图6所示；Ct≈Σk=-(N-1)/2(N-1)/2C^keikωt---(2)]]>α＝4.0+4.2sin(20πt)(3)根据图6可知，采用本发明的方法得到的结果与试验值吻合较好。根据图7可知，采用本发明的方法的CPU计算的时间计算约为传统双时间步方法CPU计算的时间的1/10，计算效率能提高一个量级以上；本发明未详细描述内容为本领域技术人员公知技术。