CN110717285A - 一种大幅六自由度运动的流固耦合模拟方法 - Google Patents

Abstract

本发明的目的在于提供一种大幅六自由度运动的流固耦合模拟方法，包括以下步骤：(1)划分结构物六自由度网格；(2)设定参数；(3)读入网格、设置边界条件和初场；(4)重构球面supermesh；(5)求解流场；(6)求解各物理量；(6)由流场求解结构物再地球坐标系下的位移；(7)更新球形区域网格节点坐标；(8)由弹簧法更新非球形区域网格节点坐标；返回步骤(4)，开始下一个时间步的计算，直至计算完成。本发明可以适用于结构物大幅度六自由度流固耦合模拟，并保持网格运动时的网格质量；并且，本发明构造的球面supermesh重构方法计算流程简单，容易实现，且计算速度快，可靠性高；突破了传统滑移网格只能处理单一自由度网格运动的局限性。

Description

一种大幅六自由度运动的流固耦合模拟方法

技术领域

本发明涉及的是一种计算机辅助工程(CAE)、计算流体力学(CFD)、流固耦合方法。

背景技术

随着计算机技术和数值方法的飞速发展，计算流体力学方法(CFD)开始活跃于空气动力学和船舶与海洋工程的各个领域当中，已有越来越多研究工作基于CFD方法开展研究。然而CFD的发展一直受困于多个瓶颈，其中一个关键问题就是动网格技术。在结构物流固耦合运动问题当中，需要处理其六个自由度的运动。

结构物在流体中具有六个自由度的运动，并且在流体的作用下，其运动规律极其复杂。因此，要准确地模拟计算结构物的运动，并考虑流固耦合作用，就必须要让网格“动”起来。

传统的动网格技术主要有两大类:

一种方案是滑移网格方法。滑移网格技术是常用的动网格方法，这种动网格方法能够处理沿直线平移或者绕着固定轴旋转等网格运动，但是只适用于简单的单个自由度运动，并且运动轨迹固定。在船舶与海洋工程中，滑移网格方法在实际应用中仅限于船舶固定、螺旋桨在船后旋转的情况，无法适应于船体的旋转运动，且一旦船舶拥有较大幅度的运动，并且同时考虑螺旋桨的旋转时，传统的滑移网格就难以适用了。

另一种方案是网格变形法。传统的网格变形有弹簧法、局部网格重构法等，也仅能处理结构物小幅度运动的情况，当发生大幅度的旋转时，网格会出现严重变形、扭曲，导致网格质量下降，甚至计算发散。

因此，传统动网格技术上的局限性使CFD在应用过程中，始终很难做到真正意义上的结构物大幅度六自由度流固耦合模拟计算。

发明内容

本发明的目的在于提供克服传统动网格技术不足的一种大幅六自由度运动的流固耦合模拟方法。

本发明的目的是这样实现的：

本发明一种大幅六自由度运动的流固耦合模拟方法，其特征是：

(1)划分结构物六自由度网格，以结构物的旋转中心为球心，将计算域划分为包含结构物的球形计算域Ω_S和球外计算域Ω_R，并分别生成计算网格τ_S和τ_R；定义τ_S和τ_R的球面交界面网格为Π_S和Π_R；

(2)设定参数：分别给定结构物的旋转中心和重心在地球坐标系下的初始坐标和

结构物的质量m、基于结构物重心的主转动惯量(I_xcg,I_ycg,I_zcg)；

(3)读入网格、设置边界条件和初场；

(4)重构球面supermesh；

(5)求解流场；

(6)求解各物理量；

(7)由当前流场和网格节点位置，求解结构物在地球坐标系下的位移 (x₁,x₂,x₃,φ,θ,ψ)，其中x₁,x₂,x₃和φ,θ,ψ分别为结构物在地球坐标系下的位移和旋转；

(8)由(x₁,x₂,x₃,φ,θ,ψ)更新结构物重心、旋转中心、网格τ_S所有节点的坐标；

(9)采用弹簧法，即将τ_R看作一张具有弹性的网，网格节点之间的联系看作是一根独立的弹簧，在球面交界面网格Π_R运动时，由刚度系数重新分布τ_R网格节点；返回步骤(4)，开始下一个时间步的计算，直至计算完成。

本发明还可以包括：

1、重构球面supermesh包括以下步骤：

(4a)分别生成Π_R和Π_S的单元列表

和

n和m分别为Π_R和Π_S的单元

和

的总数；

(4b)遍历

的单元

对遍历的单元

由AABB-BOX算法在地球坐标系下判断

和

是否空间重叠；

(4c)若重叠，采用球心投影，将地球坐标系下

和的投射至承影平面得到

和

(4d)若不重叠，则返回步骤(4b)；

(4e)由逐边裁剪法对和

进行相交计算，得到相交部分

(4f)由逆球心投影将

转换为三维笛卡尔坐标系下的K^int；

(4g)在地球坐标系下根据K^int和

的面积比值判断K^int是否为有效相交；

(4h)若有效相交，则计算K^int的几何信息，并储拓扑关系；将

加入supermesh列表

(4i)若无效相交，则返回步骤(4b)。

2、步骤(4b)在地球坐标系下由AABB-BOX算法判断

和

是否重叠包括以下步骤：

(4ba)分别对和

建立直平行六面体

和

所述ω指对于由N个节点x_i＝(x_i,y_i,z_i),i＝1,2,...,N组成的空间平面K，分别取点 A＝(max{x_i},max{y_i},max{z_i})和点B＝(min{x_i},min{y_i},min{z_i})， i＝1,2,...,N，建立以线段AB为对角线的直平行六面体，且其面均垂直于坐标轴；

(4bb)若

和

有重叠部分，则

和

为空间重叠；反之，则无空间重叠。

3、步骤(4c)球心投影包括以下步骤：

(4ca)旋转球面，使球面任意位置的

和

投影至承影平面内；球面旋转为绕过球心的旋转轴的一次旋转，或绕过以球心为原点的笛卡尔坐标系的坐标轴的组合旋转；

(4cb)由球心投影将

和

投影至乘影平面，得到极坐标系下的多边形

和

(4cc)将极坐标系变换为笛卡尔坐标系，得到笛卡尔坐标系下的

和

4、步骤(7)中求解结构物在地球坐标系下的位移包括以下步骤：

(7a)规定结构物绕固定于结构物的结构物坐标系的坐标轴的转动顺序，并计算该顺序下的转换矩阵J₁和J₂；

(7b)计算结构物在地球坐标系下所受的力F_e和力矩M_e

(7c)将F_e和M_e转换至结构物坐标系下，得到F_s和M_s：

(7d)求解六自由度刚体运动方程:

得到结构物在结构物坐标系下的速度v＝(v₁,v₂)＝(u,v,w,p,q,r)；式中

为绕j结构物旋转中心的主转动惯量，x_g,y_g,z_g为结构物重心到旋转中心的向量在结构物坐标系下的三个分量；

(7e)将v＝(v₁,v₂)转换为地球坐标系下的速度

(7f)将加速度积分得到η＝(η₁,η₂)＝(x₁,x₂,x₃,φ,θ,ψ)。

5、步骤(8)由(x₁,x₂,x₃,φ,θ,ψ)更新结构物重心坐标、旋转中心坐标，网格τ_S所有节点坐标包括以下步骤：

(8a)根据(φ,θ,ψ)，计算由步骤(3a)规定的转动顺序所对应的向量旋转矩阵J；

(8b)根据(x₁,x₂,x₃)和旋转矩阵J，计算结构物重心、旋转中心、网格τ_S所有节点的坐标。

6、步骤(9)采用弹簧法，重新分布网格τ_R所有节点包括以下步骤：

(9a)计算网格τ_R任意节点i与j之间的刚度系数α_ij；

(9b)计算球面网格Π_R所有节点的位移

并给定网格τ_R所有边界节点位移

(9c)对于网格τ_R的任意节点i，其与相连的所有节点j之间的相对位移

满足胡克定律：

(9d)由胡克定律计算网格τ_R所有节点的位移，

(9e)更新所有节点的坐标。

本发明的优势在于：本发明可以适用于结构物大幅度六自由度流固耦合模拟，并保持网格运动时的网格质量；并且，本发明构造的球面supermesh重构方法计算流程简单，容易实现，且计算速度快，可靠性高；突破了传统滑移网格只能处理单一自由度网格运动的局限性。

附图说明

图1为本发明流程图；

图2为潜艇六自由度网格划分示意图；

图3飞行器六自由网格划分示意图；

图4为重构球面supermesh的计算流程图；

图5为AABB-BOX算法流程图；

图6为球心投影处理球面相交单元示意图；

图7为二维平面相交计算流程图；

图8为二维平面Sutherland-Hodgma算法示意图；

图9为逆球心投影算法流程图。

具体实施方式

下面结合附图举例对本发明做更详细地描述：

本发明公开了一种大幅六自由度运动的流固耦合模拟方法，包括以下步骤：

(1)划分结构物六自由度网格，以结构物的旋转中心为球心，将计算域划分为包含结构物的球形计算域Ω_S和球外计算域Ω_R，并分别生成计算网格τ_S和τ_R；定义τ_S和τ_R的球面交界面网格为Π_S和Π_R；

(2)设定参数；分别给定结构物的旋转中心和重心在地球坐标系下的初始坐标

和

结构物的质量m，基于结构物重心的主转动惯量(I_xcg,I_ycg,I_zcg)；

(3)读入网格、设置边界条件和初场；

(4)重构球面supermesh；

(5)求解流场；

(6)求解各物理量；

(7)由当前流场和网格节点位置，求解结构物在地球坐标系下的位移 (x₁,x₂,x₃,φ,θ,ψ)，其中x₁,x₂,x₃和φ,θ,ψ分别为结构物在地球坐标系下的位移和旋转；

(8)由(x₁,x₂,x₃,φ,θ,ψ)更新结构物重心、旋转中心、网格τ_S所有节点的坐标；

(9)采用弹簧法，即将τ_R看作一张具有弹性的网，网格节点之间的联系看作是一根独立的具有一定刚度系数的弹簧，在球面交界面网格Π_R运动时，由刚度系数重新分布τ_R网格节点；返回步骤(4)，开始下一个时间步的计算，直至计算完成；

步骤(4)重构球面supermesh包括以下步骤：

(4a)分别生成Π_R和Π_S的单元列表

和

n和m分别为Π_R和Π_S的单元

和

的总数；

(4b)遍历

的单元

对

遍历

的单元由AABB-BOX算法在三维笛卡尔坐标系(即地球坐标系)下判断和是否空间重叠；

(4c)若重叠，采用球心投影(又称日晷投影或大圆投影或心射切方位投影)，将三维笛卡尔坐标系下

和

的投影至投影平面得到

和

(4d)若不重叠，则返回(b)；

(4e)由逐边裁剪法(Sutherland-Hodgma算法)对

和

进行相交计算，得到相交部分

(4f)由逆球心投影将

转换为三维笛卡尔坐标系下的K^int；

(4g)在三维笛卡尔坐标系下根据K^int和

的面积比值判断K^int是否为有效相交；

(4h)若有效相交，则计算K^int的几何信息，并储拓扑关系；将

加入 supermesh列表

(4i)若无效相交，则返回(b)；

步骤(4b)在三维笛卡尔坐标系下由AABB-BOX算法判断

和

是否重叠包括以下步骤：

(4ba)分别对和

建立直平行六面体和

所述ω指对于由N 个节点x_i＝(x_i,y_i,z_i),i＝1,2,...,N组成的空间平面K，分别取点 A＝(max{x_i},max{y_i},max{z_i})和点B＝(min{x_i},min{y_i},min{z_i})， i＝1,2,...,N，建立以线段AB为对角线的直平行六面体，且其面均垂直于坐标轴；

(4bb)若

和

有重叠部分，则

和

为空间重叠；反之，则无空间重叠；

步骤(4c)所述的球心投影包括以下步骤：

(4ca)旋转球面，使和

在投影平面的投影

和易于相交计算；由于球心投影方法只能投影半球以内的单元，且方位角在接近π/2处的单元节点会投影至无穷远处，无法进行相交计算，故通过旋球面或改变投影平面的位置(旋转球面较改变投影平面位置简单且易于编程实现)，可使球面任意位置的

和投影至投影平面内，且便于相交计算；球面旋转可以为绕过球心的旋转轴的一次旋转，或绕过以球心为原点的笛卡尔坐标系的坐标轴的组合旋转；

(4cb)由球心投影将

和投影至乘影平面，得到极坐标系下的多边形

和

由于球面投影可将任意大圆弧投影为投影平面内的线段，故对于分别由

和

的节点构成的球面多边形

和

其在投影平面内的投影为多边形；

(4cc)将极坐标系变换为笛卡尔坐标系，得到笛卡尔坐标系下的

和

步骤(7)中求解结构物在地球坐标系下的位移包括以下步骤：

(7a)规定结构物绕固定于结构物的结构物坐标系的坐标轴的转动顺序，并计算该顺序下的转换矩阵J₁和J₂；

(7b)计算结构物在地球坐标系下所受的力F_e和力矩M_e

(7c)将F_e和M_e转换至结构物坐标系下，得到F_s和M_s：

(7d)求解六自由度刚体运动方程:

得到结构物在结构物坐标系下的速度v＝(v₁,v₂)＝(u,v,w,p,q,r)；式中

为绕j结构物旋转中心的主转动惯量，x_g,y_g,z_g为结构物重心到旋转中心的向量在结构物坐标系下的三个分量；

(7e)将v＝(v₁,v₂)转换为地球坐标系下的速度η＝(η₁,η₂)＝(x₁,x₂,x₃,φ,θ,ψ)：

(7f)将加速度积分得到η＝(η₁,η₂)＝(x₁,x₂,x₃,φ,θ,ψ)；

步骤(8)由(x₁,x₂,x₃,φ,θ,ψ)更新结构物重心坐标、旋转中心坐标，网格τ_S所有节点坐标包括以下步骤：

(8a)根据(φ,θ,ψ)，计算由步骤(3a)规定的转动顺序所对应的向量旋转矩阵J；

(8b)根据(x₁,x₂,x₃)和旋转矩阵J，计算结构物重心、旋转中心、网格τ_S所有节点的坐标；

步骤(9)采用弹簧法，重新分布网格τ_R所有节点包括以下步骤：

(9a)计算网格τ_R任意节点i与j之间的刚度系数α_ij；

(9b)计算球面网格Π_R所有节点的位移

并给定网格τ_R所有边界节点位移

(9c)对于网格τ_R的任意节点i，其与相连的所有节点j之间的相对位移

满足胡克定律：

(9d)由胡克定律计算网格τ_R所有节点的位移，

(9e)更新所有节点的坐标。

结合图1-9，本具体实施方式包括如下步骤：

(1)划分结构物六自由度网格；潜艇和飞行器的六自由度网格划分分别如图2和图3所示，计算中使用两个坐标系：地球坐标系oxzy和结构物坐标系 OXYZ，地球坐标系静止不动，结构物坐标系固定于结构物，其原点O在结构物的旋转中心；结构物坐标系相对于地球坐标系的位置为η＝(η₁,η₂)＝(x₁,x₂,x₃,φ,θ,ψ)，其中x₁,x₂,x₃和φ,θ,ψ分别为结构物绕地球坐标系x,y,z轴的位移和旋转角度；当结构物正浮时，结构物坐标系的坐标轴的方向与地球坐标系相同；以结构物的旋转中心为圆心，取半径为R的球，使结构物置于球形区域Ω_S内，并在结构物表面和球面围成的区域生成网格τ_S；取足够大的长方体区域Ω_R包含区域Ω_S，在Ω_R和Ω_S的边界形成区域内生成网格τ_R；定义Ω_R和Ω_S的一对球面交界面网格为Π_R和Π_S；

(2)给定结构物的旋转中心和重心在地球坐标系下的初始坐标

和

结构物的质量m，基于结构物重心的主转动惯量(I_xcg,I_ycg,I_zcg)；

(3)读入网格、设置边界条件和初场；

(4)重构球面supermesh；

(5)求解流场；

(6)求解各物理量；

(7)由当前流场和网格节点位置，求解结构物在地球坐标系下的位移 (x₁,x₂,x₃,φ,θ,ψ)，其中x₁,x₂,x₃和φ,θ,ψ分别为结构物在地球坐标系下的位移和旋转；

(8)采用弹簧法，即将τ_R看作一张具有弹性的网，网格节点之间的联系看作是一根独立的具有一定刚度系数的弹簧，在球面交界面网格Π_R运动时，由刚度系数重新分布τ_R网格节点；返回步骤(4)，开始下一个时间步的计算，直至计算完成；

(9)求解流场，回到步骤(3)，开始下一个时间步的计算；

步骤(4)中：重构球面supermesh。

(4a)如图4所示，分别生成Π_R和Π_S的单元列表

和

n和m分别为Π_R和Π_S的单元

和

的总数；

(4b)遍历

的单元

对

遍历

的单元

由AABB-BOX算法在三维笛卡尔坐标系下判断

和

是否空间重叠；

(4c)若重叠，采用球心投影，将三维笛卡尔坐标系下

和

的投影至投影平面得到

和

(4d)若不重叠，则返回(b)；

(4e)由逐边裁剪法(Sutherland-Hodgma算法)对

和

进行相交计算，得到相交部分

(4f)由逆球心投影将

转换为三维笛卡尔坐标系下的K^int；

(4g)在三维笛卡尔坐标系下计算K^int和

的面积比值，若则为K^int无效相交；反之，则有效相交；V₁为给定的小量；

(4h)若有效相交，则计算

的面心坐标、面积矢量，储存

和

的面编号和其所属的单元编号；将

加入supermesh列表

(4i)若无效相交，则返回(b)；

步骤(4b)中：在三维笛卡尔坐标系下由AABB-BOX判断和

是否重叠。：

(4ba)如图5所示，由

和

节点坐标计算点A^R、B^R、A^S、B^S：

(4bb)若满足下列之一：

即判定

和

不重叠；反之，则判定

和重叠；V₂为给定的小量；

步骤(4c)中：由球心投影将三维笛卡尔坐标系下和的投影至投影平面得到

和

(4ca)如图6所示，以结构物旋转中心

为原点，建立局部三维笛卡尔坐标系x′y′z′，其坐标轴与地球坐标系平行；将

和

的节点

和

转换为局部坐标系x′y′z′下的坐标

和取的第一个点P₁ ^R′为旋转点，计算旋转轴r：

其中z′＝(0,0,1)；

(4ca)计算P₁ ^R′在以结构物旋转中心

为原点的球坐标系下的极角θ；

(4ca)绕旋转轴r按右手旋转至切点O所转的角度Δθ：

Δθ＝π-θ

θ为旋转点P在以结构物旋转中心

为原点的球坐标系下的极角；

(4cb)计算绕旋转轴r按右手旋转Δθ的旋转矩阵J_r：

其中(x,y,z)＝(x′_r,y′_r,z′_r),θ＝Δθ；

(4cc)将

和绕r按右手旋转Δθ得到和

P″＝J_r·P′

(4cd)将

和转为球坐标系下的

和

(4ce)由球心投影，将

和

投影至投影平面，并通过坐标系转换得到二维笛卡尔坐标系下的和

θ和

分别为P_i ^R″′、

在球下坐标系下的极角和方位角；

和

分别为

和

节点；

步骤(4e)中：由逐边裁剪法(Sutherland-Hodgma算法)对

和进行相交计算。

(4ea)如图7所示，顺时针分别定义

和

的节点为P₁,P₂,...,P_N和 S₁,S₂,...,S_M，产生线段

和

其中

(4eb)依次遍历

计算垂直于

且指向

外侧的矢量n：

点C为

的几何中心；

(4ec)对于依次遍历

分别判断P_n和P_n+1是否在

内侧；

(4ed)如图8所示的四种情况：若P_n和P_n+1均在内侧，则输出P_n+1；若 P_n和P_n+1均不在内侧，不输出点；若P_n在

内侧，P_n+1不在内侧，则输出交点；若P_n不在

内侧，P_n+1在

内侧，则输出P_n+1和交点；返回(c)；

(4ee)将输出的点记为P₁,P₂,...,P_N，返回(b)

(4ef)遍历完后，输出的点即为和

的相交部分

的节点，记为Q_i,i＝1,2,...,N；

步骤(4f)中：逆球心投影将

转换为三维笛卡尔坐标系下的K^int。

(4fa)如图9所示，将

的节点Q_i,i＝1,2,...,N转换为球坐标系下的Q_i″′：

(4fb)将Q_i″′转为以球心为原点的笛卡尔坐标系下的Q_i″；

(4fc)将Q_i″绕r旋转按左手旋转Δθ得到Q_i′：

(4fc)将Q_i′转换为统一笛卡尔坐标系下的P_i；

(4fd)输出由P_i组成的和

相交部分K^int；

步骤(7)中求解结构物在地球坐标系下的位移包括以下步骤：

(7a)上指标n表示时刻，规定结构物依次绕结构物坐标系Z轴、Y轴、X 轴旋转，在该顺序下，由上一时间刻(φ,θ,ψ)^n-1计算转换矩阵J₁和J₂：

(7b)将结构物所有表面单元受到的流体动力dF_e积分，得到结构物所受的力和力矩F_e和力矩M_e：

F_e＝mg+∫_shipdF_e

M_e＝r_cg×mg+∫_shipr×dF_e

dF_e＝τ·dS_f+pdS_f

τ＝-ρν_eff(▽U+▽U^T)

p＝p_d+ρg·x

τ为剪切应力，p为总压，U为结构物壁面流体相对于结构物壁面的速度， g为重力加速度，x为水深，ν_eff为粘度，ρ为水密度；

(7c)将F_e和M_e转换至结构物坐标系下，得到F_s和M_s：

(7d)将六自由度刚体运动方程线性化处理为：

x＝A^-1b

采用四阶龙格库塔法由上一时刻(u,v,w,p,q,r)^n-1求解得到当前时刻速度 (u,v,w,p,q,r)ⁿ；

(7e)将vⁿ＝(v₁,v₂)ⁿ＝(u,v,w,p,q,r)ⁿ转换为地球坐标系下的速度ηⁿ＝(η₁,η₂)ⁿ＝(x₁,x₂,x₃,φ,θ,ψ)ⁿ：

η₁＝J₁·v₁ η₂＝J₂·v₂

(7f)将加速度积分得到ηⁿ＝(η₁,η₂)ⁿ＝(x₁,x₂,x₃,φ,θ,ψ)ⁿ：

ηⁿ＝η^n-1+ηⁿdt

dt为时间步长；

步骤(8)由(x₁,x₂,x₃,φ,θ,ψ)ⁿ更新结构物重心坐标、旋转中心坐标，网格τ_S所有节点坐标包括以下步骤：

(8a)根据(φ,θ,ψ)ⁿ，计算由步骤(7a)规定的转动顺序所对应的旋转矩阵 J:

(8b)计算当前时刻结构物旋转中心

重心

网格τ_S所有节点坐标

和

分别为网格τ_S所有节点坐标初始坐标和当前时间步的坐标；

步骤(9)采用弹簧法，重新分布网格τ_R所有节点包括以下步骤：

(9a)计算网格τ_R所有相连节点i和j之间的刚度系数α_ij：

α_ij＝|x_j-x_i|^-N

一般取N＝2；

(9b)计算球面网格Π_R所有节点的位移

(9b)给定网格τ_R所有边界节点位移

(9c)根据胡克定律，由Jacobi迭代求解计算网格τ_R所有内部节点i的位移

为与内部节点i相连的所有节点的位移，上标k表示迭代次数，k可直接给定或通过残差限制；

(9d)更新网格τ_R所有节点的坐标：

xⁿ＝x^n-1+δ。

Claims (7)

1.一种大幅六自由度运动的流固耦合模拟方法，其特征是：

(1)划分结构物六自由度网格，以结构物的旋转中心为球心，将计算域划分为包含结构物的球形计算域Ω_S和球外计算域Ω_R，并分别生成计算网格τ_S和τ_R；定义τ_S和τ_R的球面交界面网格为Π_S和Π_R；

(2)设定参数：分别给定结构物的旋转中心和重心在地球坐标系下的初始坐标

和

结构物的质量m、基于结构物重心的主转动惯量(I_xcg,I_ycg,I_zcg)；

(3)读入网格、设置边界条件和初场；

(4)重构球面supermesh；

(5)求解流场；

(6)求解各物理量；

(7)由当前流场和网格节点位置，求解结构物在地球坐标系下的位移(x₁,x₂,x₃,φ,θ,ψ)，其中x₁,x₂,x₃和φ,θ,ψ分别为结构物在地球坐标系下的位移和旋转；

(8)由(x₁,x₂,x₃,φ,θ,ψ)更新结构物重心、旋转中心、网格τ_S所有节点的坐标；

(9)采用弹簧法，即将τ_R看作一张具有弹性的网，网格节点之间的联系看作是一根独立的弹簧，在球面交界面网格Π_R运动时，由刚度系数重新分布τ_R网格节点；返回步骤(4)，开始下一个时间步的计算，直至计算完成。

2.根据权利要求1所述的一种大幅六自由度运动的流固耦合模拟方法，其特征是：重构球面supermesh包括以下步骤：

(4a)分别生成Π_R和Π_S的单元列表

和

n和m分别为Π_R和Π_S的单元

和

的总数；

(4b)遍历

的单元

对

遍历

的单元

由AABB-BOX算法在地球坐标系下判断

和

是否空间重叠；

(4c)若重叠，采用球心投影，将地球坐标系下

和

的投射至承影平面得到

和

(4d)若不重叠，则返回步骤(4b)；

(4e)由逐边裁剪法对

和进行相交计算，得到相交部分

(4f)由逆球心投影将

转换为三维笛卡尔坐标系下的K^int；

(4g)在地球坐标系下根据K^int和

的面积比值判断K^int是否为有效相交；

(4h)若有效相交，则计算K^int的几何信息，并储拓扑关系；将

加入supermesh列表

(4i)若无效相交，则返回步骤(4b)。

3.根据权利要求2所述的一种大幅六自由度运动的流固耦合模拟方法，其特征是：步骤(4b)在地球坐标系下由AABB-BOX算法判断

和是否重叠包括以下步骤：

(4ba)分别对

和

建立直平行六面体

和

所述ω指对于由N个节点x_i＝(x_i,y_i,z_i),i＝1,2,...,N组成的空间平面K，分别取点A＝(max{x_i},max{y_i},max{z_i})和点B＝(min{x_i},min{y_i},min{z_i})，i＝1,2,...,N，建立以线段AB为对角线的直平行六面体，且其面均垂直于坐标轴；

(4bb)若

和

有重叠部分，则

和

为空间重叠；反之，则无空间重叠。

4.根据权利要求3所述的一种大幅六自由度运动的流固耦合模拟方法，其特征是：步骤(4c)球心投影包括以下步骤：

(4ca)旋转球面，使球面任意位置的

和

投影至承影平面内；球面旋转为绕过球心的旋转轴的一次旋转，或绕过以球心为原点的笛卡尔坐标系的坐标轴的组合旋转；

(4cb)由球心投影将和投影至乘影平面，得到极坐标系下的多边形

和

(4cc)将极坐标系变换为笛卡尔坐标系，得到笛卡尔坐标系下的

和

5.根据权利要求4所述的一种大幅六自由度运动的流固耦合模拟方法，其特征是：步骤(7)中求解结构物在地球坐标系下的位移包括以下步骤：

(7a)规定结构物绕固定于结构物的结构物坐标系的坐标轴的转动顺序，并计算该顺序下的转换矩阵J₁和J₂；

(7b)计算结构物在地球坐标系下所受的力F_e和力矩M_e

(7c)将F_e和M_e转换至结构物坐标系下，得到F_s和M_s：

(7d)求解六自由度刚体运动方程:

得到结构物在结构物坐标系下的速度v＝(v₁,v₂)＝(u,v,w,p,q,r)；式中

为绕j结构物旋转中心的主转动惯量，x_g,y_g,z_g为结构物重心到旋转中心的向量在结构物坐标系下的三个分量；

(7e)将v＝(v₁,v₂)转换为地球坐标系下的速度

(7f)将加速度积分得到η＝(η₁,η₂)＝(x₁,x₂,x₃,φ,θ,ψ)。

6.根据权利要求5所述的一种大幅六自由度运动的流固耦合模拟方法，其特征是：步骤(8)由(x₁,x₂,x₃,φ,θ,ψ)更新结构物重心坐标、旋转中心坐标，网格τ_S所有节点坐标包括以下步骤：

(8a)根据(φ,θ,ψ)，计算由步骤(3a)规定的转动顺序所对应的向量旋转矩阵J；

(8b)根据(x₁,x₂,x₃)和旋转矩阵J，计算结构物重心、旋转中心、网格τ_S所有节点的坐标。

7.根据权利要求6所述的一种大幅六自由度运动的流固耦合模拟方法，其特征是：步骤(9)采用弹簧法，重新分布网格τ_R所有节点包括以下步骤：

(9a)计算网格τ_R任意节点i与j之间的刚度系数α_ij；

(9b)计算球面网格Π_R所有节点的位移

并给定网格τ_R所有边界节点位移

(9c)对于网格τ_R的任意节点i，其与相连的所有节点j之间的相对位移

满足胡克定律：

(9d)由胡克定律计算网格τ_R所有节点的位移，

(9e)更新所有节点的坐标。

Images