[发明专利]CFD/CSD耦合求解非线性气动弹性仿真方法

权利要求书

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

K‾=a0Mn+a1Cn+KT,n]]>

ΔF=Fsn+1+Mn(a2u·n+a3u··n)+Cn(a4u·n+a5u··n)-Fi,n]]>

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

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

K‾iΔui=ΔFi]]>

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

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

un+1i=un+1i-1+Δuiu·n+1i=a1Δui-a4Δuni-a5u··niu··n+1i=a0Δui-a2Δuni-a3u··ni]]>

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

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

ψi=Fsn+1-(Mu··n+1i+Fi,n+1i);]]>

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

K‾δui=ψ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技术将计算模型表面上的气动网格点的位移插值计算整个空间流场的网格点的位移；

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