CN102012953A - Cfd/csd耦合求解非线性气动弹性仿真方法 - Google Patents

Abstract

本发明公开了一种CFD/CSD耦合求解非线性气动弹性仿真方法，首先进行计算初始化，然后求解非定常气动力的CFD，然后进行边界载荷信息转换，求解非线性结构动态响应的CSD，以结构有限元网格点位移的响应为收敛准则判断是否退出计算，若不退出，则进行计算模型表面位移信息的转换后进行流场动态网格变形，然后继续计算，直到满足收敛准则，计算结束。本发明解决了气动非线性、结构非线性耦合的问题和几何大变形动网格问题，可应用于航空航天类飞行器非线性气动弹性分析和民用高层建筑和桥梁非线性气动弹性分析。

Description

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

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

$\mathrm{\Δ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)求解以下子迭代格式的平衡关系

${\stackrel{\‾}{K}}^i{\mathrm{\Δu}}^i={\mathrm{\Δ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{\Δu}}^i\\ {\stackrel{\·}{u}}_{n+1}^i=a_1{\mathrm{\Δu}}^i-a_4{\mathrm{\Δu}}_n^i-a_5{\stackrel{\·\·}{u}}_n^i\\ {\stackrel{\·\·}{u}}_{n+1}^i=a_0{\mathrm{\Δu}}^i-a_2{\mathrm{\Δ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次迭代时的内力矩阵

，求解非平衡力

${\mathrm{\ψ}}^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ⁱ⁺¹

$\stackrel{\‾}{K}{\mathrm{\δu}}^i={\mathrm{\ψ}}^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{\Δ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向位移变化历程。

具体实施方式

方法实施例：以某大展弦比机翼的气动弹性一体化的分析为例，并针对附图流程，说明利用CFD/CSD耦合求解非线性气动弹性问题的具体实施方法。

步骤一：计算初始化

(1)坐标原点取在机翼前缘顶点，X轴沿弦长指向后方，Y轴垂直于X轴，Z轴沿展长方向。利用CAD软件建立三维实体模型，该机翼是由两条翼梁、翼肋、桁条、蒙皮等构件铆接而成，翼肋铆在翼梁的腹板上，桁条铆在蒙皮上，蒙皮则铆在翼梁、翼肋等构件上；

(2)导入CAD实体模型，划分该机翼的结果有限元网格，共有972个结构点，包含2731个三点壳单元，材料弹性模量为7.0e+10，泊松比为0.3，材料密度为2820kg/m³。翼梁、翼肋、蒙皮的单元之间用公共点连接起来并对翼根部的点加六自由度的约束，即将翼根部固支，整个模型共有4440个自由度。附图4为该机翼的蒙皮结构有限元网格和机翼的梁和肋板的结构有限元网格单元；

(3)导入CAD实体模型划分CFD计算的气动网格体系，附图5为机翼的表面气动网格(点数为165×30)和机翼的空间C-H型流场网格(网格数为165×50×50)；

(4)设置计算状态，计算条件取6km高度、10°攻角下以巡航速度150m/s飞行，其马赫数为0.47，雷诺数10361947，温度为249.15K，密度为0.65959Kg/m³，声速为316.43m/s。CFD无量纲时间步长取0.001，残差收敛标准取1.0e-6；

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

(1)计算该模型流场六面体网格的面积及体积；

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

(3)对该机翼模型流场网格的内边界面计算模型表面施加物面边界条件，内边界的其他边界面和对称面施加对称边界条件，对外边界面施加远场边界条件；

(4)本实例中选择了Baldwin-Lomax湍流模式；

(5)给定CFL＝2.5，求解流场计算的当地时间步长；

(6)对于无粘项计算，由于计算状态为亚音速流场，本实例中采用了Jameson中心离散格式对N-S方程进行空间离散；

(7)考虑到流场为亚音速流场，本实例采用Runge-Kutta的子迭代格式推进求解；

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

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

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

将机翼模型表面气动网格点上的非定常气动载荷转换到结构有限元网格的点上，形成等效节点载荷。附图6显示了初始时刻气动表面网格上的载荷图及载荷信息转换后有限元网格上点的等效载荷分布图。

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

(1)本实例初始计算(n＝1)时结构点的位移、速度及加速度赋0。耦合求解时以n时刻的有限元网格节点的位移、速度及加速度的收敛值作为n+1时刻的计算初值；

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

(3)基于CR理论求解结构在当前平衡状态下的总体刚度矩阵、内力列阵，结构质量取集中质量阵，局部坐标定义见附图7；

(4)由第n时刻的总体切线刚阵K_T，n、质量阵M_n、阻尼阵C_n、速度

、加速度ü_n和n+1时刻的等效外载荷Fs_n+1等形成等效刚度矩阵和等效载荷列阵ΔF；

(5)求解子迭代格式的载荷和位移平衡关系，得到第i次迭代位移增量初值Δuⁱ；

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

，系数β＝0.25；

(7)计算n+1时刻第i次迭代时的内力矩阵，并求解非平衡力ψⁱ；

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

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

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

及加速度ü_n+1。

步骤五：判断收敛

本实例计算中，以机翼翼尖特征点上的位移为|u_n+1|≤10^-6作为收敛标准。

步骤六：若不满足收敛标准，边界位移信息转换

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

(2)将结构节点上n+3/2时刻的位移转换到计算模型表面上的气动网格点上，附图8为收敛时刻结构有限元网格点上的变形和进行位移信息转换后表面气动网格点上的变形图。

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

利用TFI技术将变形后的气动表面网格运动插值到整个流场域，计算整个域的动态网格变形。附图9显示了初始及收敛时刻的流场网格剖面图。

步骤八：返回到步骤二计算，直到满足收敛标准。

本实例计算中迭代计算n＝812步后满足收敛标准(步骤五)，计算结束。

利用计算结果的数据绘制流场云图和计算曲线，并对结果进行分析。

(1)表1比较了机翼翼梢处最大变形，与MSC.Patran的结果做了对比，可以看出Patran软件计算得到的结构变形最大，这是由于其气动力模型是平板模型，从而计算得到结构点上的载荷值较大，因此得到的结构变形较大。刚性机翼虽然建立了较为精确的气动模型，但没有考虑弹性效应，变形也较大。而用此耦合方法计算得到的弹性机翼变形考虑了弹性效应，即加入惯性力和阻尼力的作用，变形相对较小，最大变形为2.4668m。附图10显示了Patran、刚性机翼与CFD/CSD耦合计算得到的机翼后缘沿Z向变形比较。

表1机翼翼梢处最大变形比较(m)

刚性机翼变形	Patran软件计算的机翼变形	CFD/CSD耦合求解弹性机翼变形
			2.6226	2.6745	2.4668

(2)图11、12分别为该机翼初始及收敛时刻表面的马赫数云图和压力云图。从这些云图的变化看出，在亚声速阶段，机翼弹性变形并没有引起这些云图较大的变化，其分布大致相似，在变形较大的翼梢位置，云图分布有所差别。图13为沿着展向的机翼升力系数的变化，即展向的载荷分布问题。可以看出在靠近翼根的地方，收敛时刻的升力系数稍大于初始值，而靠近翼梢的地方，小于初始值。这种弹性效应使得机翼的总升力的作用线向机翼根部移动，从结构强度来看，这种大展弦比机翼的弹性变形对载荷分布是有利的。图14记录了翼梢剖面初始时刻与收敛时候的变形情况，在Z方向有明显的变形，而X方向变形相对小，这是由于在Z向有较大载荷的缘故。翼梢变形过程中的没有明显地转动。

图15、16、17分别是升力、阻力、力矩系数变化图，可以看出，升力系数和阻力系数都稍有增加，初始时刻的升阻比约为19.1058，收敛时刻的升阻比约为19.3405，考虑弹性效应后机翼的升阻比变大。图18、19、20分别记录了翼梢处特征点上X、Y和Z方向的位移变化情况。由图可以看出，这些曲线都有振荡过程，而且在刚开始时振荡比较大，原因是在刚开始的时间步长里位移变化大，速度变化也大，从而导致阻尼力和惯性力有较大的变化。经过2s左右，这些振荡逐步减小，最终达到收敛。对于这个大展弦比机翼来说，其Z向位移最大，这也是其在Z向的气动力大的缘故。

对于该大展弦比机翼，在空气动力载荷作用下结构响应较大，弹性效应明显。特别是在翼梢处有很大的变形，通过分析发现这种弹性效应并没有对机翼的结构和气动特性产生不良的影响，相反，弹性效应作用下的载荷分布对结构强度来说是有利的，而且提高了升阻比。

Claims (1)

1.CFD/CSD耦合求解非线性气动弹性仿真方法，其特征在于包括下述步骤：

步骤一：建立计算对象的实体模型，划分该实体模型的结构有限元网格，定义结构有限元网格的单元属性和材料特性，设置结构有限单元的边界条件；划分该实体模型表面气动网格和流场网格，设置计算状态；

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

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

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

(3)对流场网格的边界面进行边界条件赋值；

(4)选择湍流模型，求解粘性系数；

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

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

(7)将N-S积分方程的半离散形式转化为关于时间的常微分方程，并引入伪时间项后设置成子迭代格式，子迭代推进中可以选取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，S^j为第j个气动网格点所对应的结构有限元网格三角单元Δs₁s₂s₃的面积坐标构成的子矩阵；

(3)对于计算模型表面上的所有气动网格点和结构有限元网格点，用矩阵Fs_n+1/2＝SFa_n+1/2来表示这种转换，Fs_n+1/2表示n+1/2时刻计算模型表面上的所有结构有限元网格点上的等效载荷，S为整体转换矩阵，由子矩阵S^j投放构成；

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

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

(2)确定结构n+1时刻等效外载荷Fs_n+1＝2F_sn+1/2-Fs_n，其中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

$\stackrel{\‾}{K}=a_0{M}_n+a_1{C}_n+K_{T,n}$

$\mathrm{\Δ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)求解以下子迭代格式的平衡关系

${\stackrel{\‾}{K}}^i{\mathrm{\Δu}}^i={\mathrm{\ΔF}}^i$

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

(6)由Δuⁱ求得下次迭代的位移

、速度

及加速度

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

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

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

，求解非平衡力

${\mathrm{\ψ}}^i={\mathrm{Fs}}_{n+1}-(M{\stackrel{\·\·}{u}}_{n+1}^i+F_{i,n+1}^i);$

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

$\stackrel{\‾}{K}{\mathrm{\δu}}^i={\mathrm{\ψ}}^i$

Δ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时刻的结构有限单元网格点上的位移

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

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

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

Images