CFD/CSD耦合求解非线性气动弹性仿真方法

技术领域

本发明属于力学数值仿真。CFD/CSD耦合涉及计算流体力学、计算结构动力学和计算机程序设计技术领域。主要应用于飞行器气动弹性分析。

背景技术

气动弹性力学主要研究流场中弹性结构变形与气动力作用之间的耦合作用，这种耦合导致的飞行器的气动弹性问题，例如操纵面反转、翼面发散、抖振、颤振等动响应，轻则缩短疲劳寿命，重则毁坏飞行器。目前采用地面实验的方法解决这一类问题还不具备技术条件，因此通过数值仿真是一条节省物力劳力和快捷的实用途径，具有重要的实际应用价值。

非线性气动弹性问题中，结构场和气动场两个不同性质的物理场在耦合界面上相互作用，彼此影响。一方面，流体问题本身涉及到大量的非线性现象，如复杂湍流运动，高速运动产生的激波、激波导致的边界层分离、非定常涡脱落、运动及演化、结构变形或振动导致的非定常流体运动等；另一方面，结构问题中涉及到非线性几何大变形、弹塑性材料非线性和接触面不定的界面非线性问题等。另外，即使各自的物理场是线性的，CFD和CSD的共同界面上的不确定耦合也会引起新的非线性问题。因此，研究CFD和CSD的耦合不是这两个问题的简单叠加。通过界面上的连续性相容条件进行耦合的CFD/CSD系统是一个高度的非线性问题，处理不当会使整个系统的计算失败。

目前国内在工程领域中的气动弹性预测，主要使用现有的线性化商业软件，对于非线性气动弹性问题，只是做一些小的局部修正，没有从根本机理上加以研究。

气动弹性分析中通常有如下近似假设：

(1)机翼在其本身平面内是绝对刚硬的，因而其本身平面内的变形很小，结构小变形对气动力影响不大。在静气动弹性分析时，把非定常气动力作为不随时间变化的量来处理。气动力计算采用线性理论处理，如面元法、涡格法、活塞理论等。

(2)飞行器机翼翼根处的机身静变形非常小，因此可以认为机翼的翼根是固定支撑的。机翼各表面网格的切向平面上的面积可以用弦平面面积来表示，机翼法线方向上的变形位移也可以用竖直方向上的变形位移表示，可以应用结构线性模态的方法来处理这类问题。

显然，这些假设是典型的线性弹性系统。目前技术不能用于跨音速、超音速、高超声速柔性物体的气动弹性分析。对于具有柔性效应、大变形的结构非线性问题也不再成立。无法处理流场的非线性和结构的非线性问题。

本发明通过计算信息传输、耦合推进技术等实现了流体与固体之间的时域耦合，采用直接对计算流体动力学方程和计算结构动力学方程耦合(CFD/CSD)求解的气动弹性仿真系统是获得非线性气动弹性分析的途径。

发明内容

为了克服现有技术无法处理流场的非线性和结构的非线性问题的不足，本发明提供一种CFD/CSD耦合时域仿真方法，采用CFD/CSD松耦合的思想，通过计算信息传输和动网格技术实现了流体与固体之间的时域耦合，解决了气动非线性、结构非线性耦合的问题和几何大变形动网格问题，可应用于航空航天类飞行器非线性气动弹性分析和民用高层建筑和桥梁非线性气动弹性分析。

本发明解决其技术问题所采用的技术方案包括以下步骤：

步骤一：计算初始化

建立计算对象的实体模型，划分该实体模型的结构有限元网格，定义结构有限元网格的单元属性和材料特性，设置有限单元的边界条件；划分该实体模型表面气动网格和流场网格，设置计算状态(马赫数、密度、温度、飞行高度、雷诺数、姿态角、速度、时间步长和残差收敛标准)。

步骤二：非定常气动力的CFD求解

求解N-S积分方程，获取非线性非定常气动载荷，具体步骤如下：

(1)计算流场网格的面积向量和体积；

(2)确定流场计算初始状态：初始计算(n＝1)时以设置计算状态的定常流场计算的收敛值为初值，耦合求解时以n-1/2时刻的流场计算收敛值作为n+1/2时刻的计算初值；

(3)对流场网格的边界面进行边界条件(物面、对称和远场边界)赋值；

(4)选择湍流模式(Baldwin-Lomax、Spalart-Allmaras或K-ωSST)，求解粘性系数；

(5)根据CFL条件确定流场计算的当地时间步长；

(6)求解N-S积分方程的通量项。对于无粘通量项，根据计算对象、计算状态、计算效率和精度选择Jameson中心离散格式、Van Leer离散格式、通量差分分裂Roe格式或对流迎风分裂AUSMpw+格式，对于有粘通量项，利用中心差分格式即可；

(7)将N-S积分方程的半离散形式转化为关于时间的常微分方程，并引入伪时间项后设置成子迭代格式，子迭代推进中可以选取Runge-Kutta显式推进或LU-SGS隐式推进方法，Runge-Kutta显式方法计算量小，易于实现，但时间步长受限，LU-SGS隐式方法无条件稳定，但计算量较大，需要综合考虑加以选择；

(8)若子迭代计算的流场残差值满足ε_n+1/2≤(10^-9～10^-6)的收敛标准，则退出子迭代，若不满足，则迭代计算，直至满足收敛标准；

(9)由模型表面的压强积分得到n+1/2时刻计算模型表面上的所有气动网格点上的非定常气动载荷Fa_n+1/2；

步骤三：边界载荷信息转换

将步骤二计算得到的非定常气动载荷转换到结构有限元网格点上，具体步骤如下：

(1)为计算模型表面上的每一个气动网格点j基于能量最小找出其相邻的一个结构有限元网格点的三角单元Δs₁s₂s₃；

(2)求解气动点j在该三角单元Δs₁s₂s₃中的投影的面积坐标α、β和γ，构成子矩阵

S^j＝[αβγ]^T    (1)

其中S^j为第j个气动网格点所对应的结构有限元网格三角单元Δs₁s₂s₃的面积坐标构成的子矩阵；

(3)对于计算模型表面上的所有气动网格点和结构有限元网格点，用矩阵形式来表示这种转换

Fs_n+1/2＝SFa_n+1/2    (2)

其中Fs_n+1/2表示n+1/2时刻计算模型表面上的所有结构有限元网格点上的等效载荷，S为整体转换矩阵，由子矩阵S^j投放构成。

步骤四：非线性结构动态响应的CSD求解

求解非线性结构动力学方程，获取结构非线性动态响应，具体步骤如下：

(1)确定结构有限元网格点的位移、速度及加速度初值：初始计算(n＝1)时以结构的平衡状态为初值，耦合求解时以n时刻的有限元网格点的位移、速度及加速度的收敛值作为n+1时刻的计算初值；

(2)确定结构n+1时刻等效外载荷

Fs_n+1＝2Fs_n+1/2-Fs_n    (3)

其中Fs_n+1、Fs_n+1/2和Fs_n分别表示第n+1、n+1/2和n时刻的等效外载荷；

(3)基于Co-rotational理论求解结构在当前平衡状态下的总体刚度矩阵和内力阵。其实施步骤为：①建立结构有限元的局部坐标；②求解局部坐标系下的刚度矩阵和内力；③将局部坐标下的刚阵和内力向总体坐标转换得到第n时刻的总体切线刚阵K_T，n和内力矩阵F_i，n；

(4)由第n时刻的总体切线刚阵K_T，n、质量阵M_n、阻尼阵C_n、速度、加速度ü_n和n+1时刻的等效外载荷Fs_n+1等形成等效刚度矩阵和等效载荷列阵ΔF

$\bar{K}=a_0{M}_n+a_1{C}_n+K_{T,n}$(4)

$\mathrm{\Delta F}={\mathrm{Fs}}_{n+1}+M_n(a_2{\stackrel{·}{u}}_n+a_3{\stackrel{··}{u}}_n)+C_n(a_4{\stackrel{·}{u}}_n+a_5{\stackrel{··}{u}}_n)-F_{i,n}$

系数a₁＝γ/(βΔt)，a₂＝1/(βΔt)，a₃＝1/(2β)-1，其中β的取值为0.25-0.5；

(5)求解以下子迭代格式的平衡关系

${\bar{K}}^i{\mathrm{\Delta u}}^i={\mathrm{\Delta F}}^i---\left(5\right)$

其中上标i表示第i次的迭代过程，得到第i次迭代位移增量初值Δuⁱ；

(6)由Δuⁱ求得下次迭代的位移、速度及加速度

$\left(\begin{array}{c}{u}_{n+1}^i=u_{n+1}^{i-1}+{\mathrm{\Delta u}}^i\\ {\stackrel{·}{u}}_{n+1}^i=a_1{\mathrm{\Delta u}}^i-a_4{\mathrm{\Delta u}}_n^i-a_5{\stackrel{··}{u}}_n^i\\ {\stackrel{··}{u}}_{n+1}^i=a_0{\mathrm{\Delta u}}^i-a_2{\mathrm{\Delta u}}_n^i-a_3{\stackrel{··}{u}}_n^i\end{array}\right)---\left(6\right)$

系数a₀＝1/(βΔt²)，a₄＝γ/β-1，a₅＝Δt/2·(γ/β-2)，其中β的取值同上(4)；

(7)由步骤(3)计算n+1时刻第i次迭代时的内力矩阵，求解非平衡力

${\psi }^i={\mathrm{Fs}}_{n+1}-(M{\stackrel{··}{u}}_{n+1}^i+F_{i,n+1}^i)---\left(7\right)$

(8)修正位移增量，求解第i次迭代位移增量的修正量δuⁱ，得到第i+1次迭代的位移增量Δuⁱ⁺¹

$\bar{K}{\mathrm{\delta u}}^i={\psi }^i$(8)

Δuⁱ⁺¹＝Δuⁱ+δuⁱ

(9)判断位移增量是否满足收敛准则(Δuⁱ⁺¹≤10^-6～10^-9)，否则返回(3)；

(10)由以上计算的收敛值得到n+1时刻的结构有限单元网格点上的位移u_n+1、速度及加速度ü_n+1。

步骡五：判断收敛

以结构有限元网格点位移u_n+1的响应为是否退出计算的判断准则：若u_n+1的变化趋向于收敛情况，则当|u_n+1|≤(10^-8～10^-6)时退出计算；若u_n+1的变化趋于等幅振荡情况，则当u_n+1响应进入稳定等幅振荡周期后退出计算；若u_n+1的变化趋于发散情况，则当|u_n+1/u_n|＞10²时退出计算。

步骤六：边界位移信息转换

若不满足步骤五的准则，进行计算模型表面位移信息的转换，步骤如下：

(1)若已知n+1时刻结构有限单元网格点上的位移u_n+1、速度，预测n+3/2时刻的结构有限单元网格点上的位移u_n+3/2

$u_{n+3/2}=u_{n+1}+\mathrm{\Delta t}/2{\stackrel{·}{u}}_{n+1}---\left(9\right)$

(2)将结构点上n+3/2时刻的位移转换到计算模型表面上的气动网格点上，位移转换公式为

x_n+3/2＝S^Tu_n+3/2+AN_s    (10)

其中x_n+3/2和u_n+3/2分别为计算模型表面气动网格和结构有限单元网格点上n+3/2时刻的位移，整体转换矩阵S与步骤三相同，矩阵A为计算模型表面上的每个气动网格点到所对应结构有限元三角单元投影的距离构成的对角矩阵，N_s为计算模型表面上的所有结构有限元网格点三角形单元的法向向量组成的向量；

步骤七：流场动态网格变形

利用TFI技术将计算模型表面上的气动网格点的位移插值计算整个空间流场的网格点的位移；

步骤八：返回到步骤二计算，直到满足步骤五所定义的收敛准则，计算结束。

本发明的有益效果是：

(1)针对现有技术无法处理既有流场非线性、又有结构非线性问题的能力。本发明步骤二和步骤四的非线性求解器的引入，具备处理弹性飞行器结构域和气动域流场非线性耦合的能力。

(2)CFD/CSD耦合方法可用于分析亚、跨、超和高超声速飞行器气动弹性的绝大部分问题：弹性变形、载荷分布求解、静稳定性分析、气动弹性动稳定性分析、动态响应模拟、极限环模拟等。

(3)CFD/CSD耦合时域计算可获得结构响应的所有信息(如：位移响应，加速度响应，速度响应，惯性载荷响应)和流场响应的所有信息(如：非定常气动载荷，压强系数，密度，温度，马赫数)。

下面结合附图和实施例对本发明进一步说明。

附图说明

图1是本发明CFD/CSD耦合求解非线性气动弹性的基本流程的示意图；

图2是本发明步骤二CFD/CSD耦合求解非线性气动弹性的非定常气动力求解的基本流程的示意图。

图3是本发明步骤四CFD/CSD耦合求解非线性气动弹性的非线性结构响应求解的基本流程的示意图。

图4是方法实施例的机翼蒙皮结构网格和机翼梁-肋板结构图。

图5是方法实施例的机翼表面气动网格和机翼C-H型流场网格图。

图6是方法实施例的初始时刻气动网格点载荷和进行载荷转换后结构有限元网格点上的等效载荷分布图。

图7是方法实施例的基于CR理论的结构单元局部坐标系定义。

图8是方法实施例收敛时刻结构有限元网格点变形和进行位移转换后的和气动网格点上的变形图。

图9是方法实施例的初始及收敛时刻的流场网格剖面图。

图10是方法实施例的Patran、刚性机翼与CFD/CSD耦合计算得到的机翼后缘沿Z向变形比较。

图11是方法实施例的初始及收敛时刻表面马赫数云图。

图12是方法实施例的初始及收敛时刻表面压力系数云图。

图13是方法实施例的初始及收敛时刻展向升力系数曲线。

图14是方法实施例的初始及收敛时刻翼梢剖面位移曲线。

图15是方法实施例的升力系数的变化历程。

图16是方法实施例的阻力系数的变化历程。

图17是方法实施例的力矩系数的变化历程。

图18是方法实施例的特征点处X向位移变化历程。

图19是方法实施例的特征点处Y向位移变化历程。

图20是方法实施例的特征点处Z向位移变化历程。