CN111177960A

CN111177960A - 一种基于ancf的薄壳碰撞接触计算方法

Info

Publication number: CN111177960A
Application number: CN201911316609.9A
Authority: CN
Inventors: 刘吉凡; 王庆涛
Original assignee: Nanjing University of Science and Technology
Current assignee: Nanjing University of Science and Technology
Priority date: 2019-12-19
Filing date: 2019-12-19
Publication date: 2020-05-19
Anticipated expiration: 2039-12-19
Also published as: CN111177960B

Abstract

本发明公开了一种基于ANCF的薄壳碰撞接触计算方法，基于有限元理论和连续介质力学理论计算薄壳单元应变能，求得单元的弹性力和弹性力的雅克比矩阵，采用主从面算法作为接触检测算法，罚函数法计算法向接触力和摩擦力，最后用广义‑

算法计算柔性多体系统动力学，求解薄壳碰撞接触问题。本发明与现有技术相比，其显著优点在于：（1）与KED法与浮动坐标法相比，本方法质量矩阵为常数阵，不存在科氏力和离心力，解决柔性多体系统动力学中的动力刚化问题，很好完成刚柔耦合，使得结果更加精确。（2）采用主从面算法作为接触检测算法，适用于正常的接触问题。（3）采用罚函数法计算法向接触力和摩擦力，将接触问题变为无约束问题，计算效率较高。

Description

一种基于ANCF的薄壳碰撞接触计算方法

技术领域

本发明属于研究柔性多体系统动力学问题，具体涉及一种基于绝对节点坐标(ANCF)的薄壳碰撞接触计算方法。

背景技术

由于薄壳等可展开结构的轻质化等优点，薄膜结构如太阳帆板等被广泛应用于现代航空航天领域中。而对于柔性系统的大范围位移和大变形之间的耦合，是学术界和工程界的难点，其数值模拟困难且迫切，随着我国航天事业发展，太空中可展开结构如薄膜结构的研究愈发迫切，基于绝对节点坐标法的薄壳结构碰撞接触计算方法有着相当的重要性。本发明能为航空航天、车辆、机械等领域的多体动力学仿真提供一定的计算帮助

王庆涛在2016年公开发表了《经历大范围运动和大变形的细梁接触动力学》，论文分析了空间结构中可展开柔细绳索的复杂接触动力学问题，但对于薄壳结构的接触碰撞问题并未做进一步的分析。

发明内容

本发明的目的在于提供一种基于ANCF的薄壳碰撞接触计算方法，对薄壳的接触碰撞动力学问题进行较好的动力学分析。

实现本发明目的技术解决方案为：一种基于ANCF的薄壳碰撞接触计算方法，包括以下步骤：

步骤1、任意选取两个由具有梯度缺陷的不同薄壳单元分别构成薄壳，基于绝对节点坐标法确定两个薄壳单元参数，只考虑中面拉压变形和剪切变形，基于Kirchhoff板壳理论，计算出每个薄壳单元的应变能，并求得薄壳单元弹性力和弹性力的雅克比矩阵，转入步骤2；

步骤2、对上述两个薄壳进行接触检测，对薄壳单元采用主从面接触算法作为接触检测算法，选取接触区域大的薄壳单元接触面为从面，另一接触区域小的薄壳单元接触面为主面，获得接触区域和最小嵌入深度，转入步骤3；

步骤3、对接触薄壳进行动力学分析，考虑切向接触力即摩擦力，由步骤2所得最小嵌入深度使用罚函数法计算法向接触力和摩擦力，根据步骤2所得接触区域求得广义接触力以及广义接触力的雅克比矩阵，转入步骤4；

步骤4、使用广义-α算法作为柔性多体系统动力学隐式求解方法，将步骤1所得薄壳单元弹性力和弹性力的雅克比矩阵，和步骤3所得的广义接触力以及广义接触力的雅克比矩阵代入多体系统动力学方程，进行数值计算，获得动力学方程的雅克比矩阵和每一时间步长下柔性多体系统的节点坐标向量。

本发明与现有技术相比，其显著优点在于：

(1)采用绝对节点坐标法，KED法与浮动坐标法相比，质量矩阵为常数阵，且不存在科氏力和离心力，解决柔性多体系统动力学中的动力刚化问题，很好完成刚柔耦合，使得结果更加精确。

(2)采用主从面算法作为接触检测算法，将接触面分为变形较大的从面和变形较小的主面，适用于正常的接触问题。

(3)采用罚函数法计算法向接触力，将接触问题变为无约束问题，不增加未知数，计算效率较高。

附图说明

图1为本发明基于ANCF的薄壳碰撞接触计算方法的流程图。

图2为本发明薄壳单元中面示意图。

图3为本发明薄壳单元上任意一点示意图。

图4为本发明计算多体系统动力学方程的广义-α算法迭代流程图。

图5为本发明薄壳主从面接触检测方法的接触点示意图。

图6为本发明验证薄壳碰撞接触动力学算法的实施例示意图。

图7为使用ANSYS商业软件LS-DYNA程序验证实施例建模图。

图8为使用ANSYS商业软件LS-DYNA程序验证实施例，参数中滑移因子为10且划分单元较小。

图9为使用ANSYS商业软件LS-DYNA程序仿真结果图。

图10为本发明计算实施例后某一时刻数值位形图。

图11为使用ANSYS商业软件LS-DYNA程序计算实施例四个点沿重力方向位移时间图。

图12为本发明计算实施例在四个点沿重力方向位移时间图。

图13为本发明计算考虑摩擦力的实施例在0.4秒时的数值位形图。

具体实施方式

结合图1，本发明所述的一种基于ANCF的薄壳碰撞接触计算方法，步骤如下：

步骤1、任意选取两个由具有梯度缺陷的不同薄壳单元分别构成薄壳，基于绝对节点坐标法确定两个薄壳单元参数包括形函数S以及薄壳单元的四个顶点的坐标向量参数e，只考虑中面拉压变形和剪切变形，基于Kirchhoff板壳理论，计算出每个薄壳单元中面上的应变能U^mid、薄壳单元中面上的弹性力Q^mid、Q^mid的雅克比矩阵J^mid，弯曲应变能U^κ、U^κ对应的弹性力Q^κ，及Q^κ的雅克比矩阵J^κ，转入步骤2：

该单元自由度为36，形函数为3×36的矩阵S，如下所示：

[S₁I S₂I … S₁₂I] (1)

其中的I为3×3的单位阵，形函数中具体表达式如下：

其中0≤ξ≤1与0≤η≤1为薄壳单元中面的正则参数a为薄壳单元中面的长度，b为薄壳单元中面的宽度，a和b由选取动力学算例的物理模型具体参数和单元数的划分确定。

薄壳单元的四个顶点的坐标向量参数e为36×1的矩阵，e＝[e_B,e_C,e_D,e_E]^T，其中e_B,e_C,e_D,e_E分别为如图2所示B点、C点、D点、E点四个顶点的坐标向量，每个顶点有9个自由度。

将有限元理论和连续介质力学理论结合，应用微分几何和坐标变换相关理论，通过张量分析推导出单元应变能，弹性力及其雅克比；

其中图2所示中面的应变能为：

ε^mid为中面上的格林拉格朗日应变张量，表示为ε^mid＝[ε₁₁ ε₂₂ 2ε₁₂]^T，E为四阶弹性张量，V为薄壳单元的体积；

由平面应力关系可得

其中E为杨氏模量，μ为泊松比，ε₁₁，ε₂₂为纵向应变，和ε₁₂为切向应变；

图3所示薄壳单元弯曲应变能为：

上式(5)中的Ω为中间变量，Ω₀为初始时刻中间变量，Ω＝-z[κ₁₁ κ₂₂ 2κ₁₂]^T且Ω₀＝-z[(κ₀)₁₁ (κ₀)₂₂ 2(κ₀)₁₂]^T。ε^κ为薄壳单元弯曲应变张量，表示为ε^κ＝-zT^T(κ-κ₀)T，z为薄壳单元上任意面与中面的垂直距离，T为坐标转换矩阵；

上式(4)中的E^κ为系数矩阵，且

H为中间变量，且

其中θ为局部坐标两坐标基向量的夹角。

其中κ为曲率向量，κ₀为初始时刻曲率向量，κ和κ₀表示为：

上式中，r为中面上任意一点的全局位置，r₀初始时刻中面上任意一点的全局位置，n为垂直于中面的单位向量，n₀为初始时刻垂直中面的单位向量，κ₁₁和κ₂₂为横向弯曲曲率分量，κ₁₂为扭转曲率分量，(κ₀)₁₁和(κ₀)₂₂为初始时刻横向弯曲曲率分量，(κ₀)₁₂为初始时刻扭转曲率分量；

中面上的弹性力为：

根据式(3)，可算得弹性力Q^mid为：

根据式(4)可得到弯曲弹性力Q^κ为：

中面上的弹性力的雅克比矩阵为：

可计算得到J^mid：

其中有：

弯曲变形弹性力为：

弹性力对应雅克比矩阵为：

由上述步骤1可得到每个薄壳单元中面上的应变能U^mid、薄壳单元中面上的弹性力Q^mid、Q^mid的雅克比矩阵J^mid，弯曲应变能U^κ、U^κ对应的弹性力Q^κ，及Q^κ的雅克比矩阵J^κ。

步骤2、对上述两个薄壳进行接触检测，对薄壳单元采用主从面接触算法作为接触检测算法，选取接触区域大的薄壳单元的为从面，另一个接触区域小的薄壳单元为主面，获得接触区域X和最小嵌入深度g，转入步骤3：

如图5所示接触面上的接触点示意图，选取主面上P点对应主薄壳的中面上i点的位置矢量r(ξ_i,η_i)，从面上Q点对应从薄壳的中面上j点的位置矢量r(ξ_j,η_j)，h_i为主面厚度，h_j为从面厚度；设两壳之间的最近点为主面上P点和从面上Q点，则PQ之间的距离称作最小嵌入深度，由g表示如下:

为求得最小嵌入深度g，要满足i点和j点的连线垂直于P、Q两点的切平面，可得以下四个方程：

采用迭代法确定最小嵌入深度g对应i点局部坐标(ξ_i,η_i)、j点的局部坐标(ξ_j,η_j)，从而得出P点的全局坐标r_P和Q点的全局坐标r_Q：

n是中面上点i、j两点的法向接触向量；

由此得到最小嵌入深度g，当最小嵌入深度g小于0时，判断发生接触，且要求计算接触区域大小：

通过接触区域边界点来确定接触区域，显然从面Q点在接触区域内，取Q点对应中面j点附近一点j'点，j'点局部坐标(ξ_j′，η_j′)，其中ξ_j′＝ξ_j+Δξ，η_j′＝η_j+Δη；通过如下式正交条件得主面的中面上对应i'点的局部坐标ξ_i′和η_i'：

其中，r(ξ_i',η_i')为i'点的位置矢量；

通过牛顿-拉弗森迭代法解出上式，得到点i'的局部坐标(ξ_i′,η_i')，根据i'、j'两点的局部坐标，得到i'点对应主面上P'点的全局坐标r_P'和j'点对应从面上Q'点的全局坐标r_Q'：

当主从面上P'、Q'两点相互嵌入深度小于容许误差tol时，则认为P'、Q'两点为接触区域的边界点，从而确定接触区域X，其条件方程为：

由上述步骤2可获得接触区域X和最小嵌入深度g。

步骤3、对接触薄壳进行动力学分析，考虑切向接触力即摩擦力，由步骤2所得最小嵌入深度g使用罚函数法计算法向接触力F_n和摩擦力F_t，根据步骤2所得接触区域X求得广义接触力P以及广义接触力的雅克比矩阵J，转入步骤4；

采用罚函数法计算步骤2接触面上的主薄壳单元中面上i点和从薄壳单元中面上j点的法向接触力如下：

F_n＝-pgn (21)

p是罚参数，g是嵌入深度。在已知法向接触力的F_n情况下，罚参数p和嵌入深度g为反比关系，罚参数p其取值要适当，需要根据经验选取；

n是主薄壳单元中面上点i和从薄壳单元中面上点j间的法向量，如下所示：

其中，Δr＝r(ξ_i,η_i)-r(ξ_j,η_j)，Δr为从i点指向j点的向量；

由虚功原理得法向接触力的广义力形式：

其中δW_n是法向接触力F_n的虚功，

其中e_i、为i接触点所在薄壳单元的节点坐标，e_i为j接触点所在薄壳单元的节点坐标；P_n是两接触点的广义法向接触力，满足：

式子中的中间变量即广义形函数

广义法向接触力P_n的雅可比矩阵J_n为：

对于切向接触力即库伦摩擦力F_t而言，摩擦力为：

其中μ_m为摩擦系数，v_t为两薄壳单元上P接触点和Q接触点相对速度；

切向接触力F_t的广义力可由虚功原理得到，则两接触点之间的广义切向接触力为：

广义切向接触力P_t的雅克比矩阵J_t为：

将广义法向力接触力P_n和广义切向接触力P_t相加，并在步骤2所得接触区域X中积分，可得广义接触力P：

广义接触力P的雅克比矩阵J为：

由步骤3可得到广义接触力P和P的雅克比矩阵J。步骤4、使用广义-α算法作为柔性多体系统动力学隐式求解方法，将步骤1所得薄壳单元弹性力和弹性力的雅克比矩阵，和步骤3所得的广义接触力以及广义接触力的雅克比矩阵代入多体系统动力学方程，使用C++编程进行数值计算，获得动力学方程的雅克比矩阵和每一时间步长下柔性多体系统的节点坐标向量，将所得坐标向量数据在中MATLAB进行分析，并和ANSYS商业软件仿真结果进行对比，以验证该方法的有效性优越性。

基于绝对节点坐标，对含约束的柔性多体系统动力学而言，其动力学方程为如下微分代数形式方程组：

其中M是整个柔性多体系统的质量矩阵，q为系统广义坐标列阵，F(q)为系统的弹性力列阵，Φ(q,t)为系统的约束方程，Φ_q为系统的约束方程对广义坐标列阵的偏导，λ为拉格朗日乘子，

为柔性多体系统的广义外力列阵。

采用隐式的广义-α算法进行求解，和牛顿迭代法类似，将上式柔性多体系统动力学方程组离散为以下代数形式的方程，具体的迭代步骤为：，

其中，t为时间；

迭代步骤中第(k+1)步的启动值如下：

上式中h是积分步长，a是牛顿迭代法的基础上新引入的算法参数阵，满足如下方程组：

β、γ、α_m和α_f为参数，取值方法如下：

其中

称为谱半径，对于算法的能量耗散而言，谱半径和柔性多体系统的能量损失成反比，当谱半径最大时系统无能量的耗散，反之系统有最大的能量耗散。具体的迭代的步骤如图4所示，图中

和

满足：

图4中

和

使得公式

和

成立。其中I为单位矩阵，tol是收敛误差；T_final为整个设定的数值计算时间。J为动力学方程的雅克比矩阵，表示如下：

其中λ为拉格朗日乘子矩阵，G为多体系统动力学方程；

对于多体系统动力学方程，由上述步骤4可得到多体系统动力学方程的雅克比矩阵J，和每一时间步长下多体系统的节点坐标。结合图6所示实施例建立数学模型并编写程序进行数值仿真计算，对所得数据进行分析，将所得结果和商业软件仿真结果进行对比，以验证该方法的有效性优越性。

实施例

步骤1、任意选取两个由具有梯度缺陷的不同薄壳单元分别构成薄壳，基于绝对节点坐标法确定两个薄壳单元参数，只考虑中面拉压变形和剪切变形，基于Kirchhoff板壳理论，计算出每个薄壳单元的应变能，并求得薄壳单元弹性力和弹性力的雅克比矩阵，其中图2为薄壳单元中面示意图，图3为薄壳单元上任意一点位置示意图，转入步骤2；

步骤2、对上述两个薄壳进行接触检测，对薄壳单元采用主从面接触算法作为接触检测算法，选取接触区域大的薄壳单元的为从面，另一个接触区域小的薄壳单元为主面，获得接触区域和最小嵌入深度，其中如图5所示为接触点示意图，实施例如图6所示，选取上方薄壁圆筒接触面所在薄壳单元为从面，下方圆筒接触面所在薄壳单元为主面，获得接触区域和最小嵌入深度，转入步骤3；

步骤4、使用如图4所示广义-α算法作为柔性多体系统动力学隐式求解方法，将步骤1所得薄壳单元弹性力和弹性力的雅克比矩阵，和步骤3所得的广义接触力以及广义接触力的雅克比矩阵代入多体系统动力学方程，使用C++编程进行数值计算，获得动力学方程的雅克比矩阵和每一时间步长下柔性多体系统的节点坐标向量，并用MATLAB对数据进行分析，和ANSYS商业软件仿真结果进行对比，以验证该方法的有效性优越性。

本发明采用C++编程计算了如图6所示两薄壁圆筒接触动力学算例，并和商业软件对如图6所示实施例仿真结果进行对比。

本实施例如图6所示，圆筒具体参数为母线长2米，周长1米1圆筒厚度0.01米，材料的杨氏模量为E＝1e7帕，泊松比μ＝0.3，材料密度为5000千克每立方米，两圆筒材料形状等参数都完全相同，下方圆筒两侧有固定约束，上方圆筒水平横放。初始时刻两圆柱筒由一微小间隙，给本实施所示多体系统有一个9.8牛顿每千克的重力场，在重力作用下两圆筒发生碰撞接触，先不考虑摩擦力。

使用ANSYS的LS-DYNA程序对实施例进行如图7和如图8所示建模，并进行动力学仿真。发现当滑移界面因子为10时，运算结果显示两圆筒发生穿透效应。如图9所示，为加大滑移界面因子为100时，且此时单元数约划分有10000个，计算运行后发现只能计算到0.57秒，此最终时刻运算结果如图10所示。可以看到传统商业软件不能很好地描述薄壳的碰撞接触动力学行为。

使用本发明所述的一种基于ANCF的薄壳碰撞接触计算方法，实施例中广义-α算法中谱半径为0.2，将每一个圆筒划分为32个单元，沿母线方向单元数为8个，沿圆周方向单元数为4，选取步长为2×10^-4秒，整个系统设定运行时间为0.8秒，系统罚参数选取为1×10⁷。若先不考虑摩擦力，使用本发明所述计算方法得到每一步长时间的两圆筒上每个单元上均匀分布的36个点的位置信息，将结果输出到MATLAB进行图像绘制以观察圆筒位移与变形情况，图10为其中某一时刻实施例所示多体系统位形图。可以看到采用本发明所述一种基于ANCF的薄壳碰撞接触计算方法能较好地描述薄壳接触碰撞的运动情况。

选取如图6所示A、B、C、D四点，其中D点位于下方圆筒最下方母线上的中点。

图11为使用ANSYS的LS-DYNA程序后，实施例A、B、C、D四点沿重力方向位移数值结果图。由图可以看到，计算结果图并未反映系统存在耗散和波动频率，其中A、B、C、D四点与实际物理情况不吻合。这也进一步说明了商业软件并不能很好地描述薄壳结构的碰撞接触动力学问题。

图12为本发明所述的一种基于ANCF的薄壳碰撞接触计算方法计算A、B、C、D四点四个点随重力方向位移的数值计算结果图。图12中可以看出，系统随时间变化逐渐趋于平衡，说明系统耗散存在作用。A、B、C、D四个点沿重力方向位移均随时间上下波动，且四个点波动频率与时间都一样，其中A点D点和C点位移较大，B点波动范围较小，其中D点位移较A点和C点又略微偏小，这与真实物理情况吻合，即两根薄壁圆筒在有规律的上下弹动，中间变形较小，两边变形较大。这进一步说明了采用本发明所述一种基于ANCF的薄壳碰撞接触计算方法能较好地描述薄壳接触碰撞的动力学行为。

在上述实施例中考虑摩擦力，给上方圆筒加一沿下方圆筒母线向右且大小为1米每秒的初速度，其他条件都不变。图13为本发明所述的一种基于ANCF的薄壳碰撞接触计算方法所得摩擦系数μ_m为0.1且时间为0.4秒时刻的薄壳碰撞接触仿真结果位形图。可以看到采用本发明所述的一种基于ANCF的薄壳碰撞接触计算方法能较好地描述考虑摩擦的薄壳接触碰撞动力学问题。

结合图11、图12和图13，对比使用商业软件ANSYS的LS-DYNA程序进行薄壳接触碰撞动力学仿真，可以看到对于柔性薄壳间的接触碰撞的动力学计算商业软件存在较大缺陷，而本发明所述一种基于ANCF的薄壳碰撞接触计算方法可以准确描述薄壳接触动力学。

Claims

1.一种基于ANCF的薄壳碰撞接触计算方法，其特征在于，包括以下步骤：

2.根据权利要求1所述的基于ANCF的薄壳碰撞接触计算方法，其特征在于：步骤1中，任意选取两个由具有梯度缺陷的不同薄壳单元分别构成薄壳，基于绝对节点坐标法确定两个薄壳单元参数包括形函数S以及薄壳单元的四个顶点的坐标向量参数e，只考虑中面拉压变形和剪切变形，基于Kirchhoff板壳理论，计算出每个薄壳单元中面上的应变能U^mid、薄壳单元中面上的弹性力Q^mid、Q^mid的雅克比矩阵J^mid，弯曲应变能U^κ、U^κ对应的弹性力Q^κ，及Q^κ的雅克比矩阵J^κ，具体如下：

形函数S表示：

[S₁I S₂I…S₁₂I]

其中，I为3×3的单位矩阵；

具体形函数参数如下：

S₁＝-(ξ-1)(η-1)(2η²-η+2ξ²-ξ-1)，S₂＝-aξ(ξ-1)²(η-1)

S₃＝-bη(ξ-1)(η-1)²，S₄＝ξ(2η²-η+2ξ²-3ξ)(η-1)

S₅＝-aξ²(ξ-1)(η-1)，S₆＝bξη(η-1)²

S₇＝-ξη(2η²-3η+2ξ²-3ξ+1)，S₈＝aξ²η(ξ-1)

S₉＝bξη²(η-1)，S₁₀＝η(ξ-1)(2η²-3η+2ξ²-ξ)

S₁₁＝aξη(ξ-1)²，S₁₂＝-bη²(ξ-1)(η-1)

其中，ξ和η为标准插值参数，a为薄壳单元中面的长度，b为薄壳单元中面的宽度；

薄壳单元的四个顶点的坐标向量参数e为：

e＝[e_B,e_C,e_D,e_E]^T

其中e_B,e_C,e_D,e_E分别为B点、C点、D点、E点四个顶点的坐标向量，每个顶点有9个自由度，e是36×1的矩阵；

中面应变能U^mid为：

式子中E为弹性张量，ε^mid为格林拉格朗日应变张量；V为薄壳单元的体积；

中面弹性力Q^mid为：

Q^mid的雅克比矩阵为：

弯曲应变能U^κ如下：

其中ε^κ为薄壳单元的弯曲应变张量，E^κ为系数矩阵，Ω、Ω₀为中间变量；

U^κ对应的弹性力Q^κ：

Q^κ对应的雅克比矩阵J^κ：

由步骤1可得到薄壳单元中面上的应变能U^mid、薄壳单元中面上的弹性力Q^mid、Q^mid的雅克比矩阵J^mid，弯曲应变能U^κ、U^κ对应的弹性力Q^κ，及Q^κ的雅克比矩阵J^κ。

3.根据权利要求1所述的基于ANCF的薄壳碰撞接触计算方法，其特征在于：步骤2中，对上述两个薄壳进行接触检测，对薄壳单元采用主从面接触算法作为接触检测算法，选取接触区域大的薄壳单元的接触面为从面，另一个接触区域小的薄壳单元接触面为主面，获得接触区域X和最小嵌入深度g，具体如下；

选取主面上P点对应主薄壳的中面上i点的位置矢量r(ξ_i,η_i)，从面上Q点对应从薄壳的中面上j点的位置矢量r(ξ_j,η_j)，h_i为主面厚度，h_j为从面厚度；设两个薄壳之间的最近点为主面上的P点和从面上的Q点，则PQ之间的距离称作最小嵌入深度，由g表示如下：

最小嵌入深度g满足i点和j点的连线垂直于P、Q两点的切平面，得以下四个方程：

采用迭代法确定最小嵌入深度g对应i点局部坐标(ξ_i,η_i)、j点的局部坐标(ξ_j,η_j)，从而得出P点的全局坐标r_P，Q点的全局坐标r_Q：

n是中面上点i、j两点的法向接触向量；

由此得到最小嵌入深度g，当最小嵌入深度g小于0时，判断发生接触，且要求计算接触区域大小；

通过接触区域边界点来确定接触区域，显然从面Q点在接触区域内，取Q点对应中面j点附近一点j'点，j'点局部坐标(ξ_j’，η_j’)，其中ξ_j’＝ξ_j+Δξ，η_j’＝η_j+Δη；通过如下式的正交条件得到在i点所在中面上且由j'点所对应的i'点局部坐标ξ_i’和η_i'：

其中，r(ξ_i',η_i')为i'点的全局位置矢量；

通过迭代法解出上式，得到点i'的局部坐标(ξ_i’,η_i')，根据i'、j'两点的局部坐标，得到i'点对应主面上P'点的全局坐标r_P'，j'点对应从面上Q'点的全局坐标r_Q'：

4.根据权利要求1所述的基于ANCF的薄壳碰撞接触计算方法，其特征在于：步骤3中，对接触薄壳进行动力学分析，考虑切向接触力即摩擦力，由步骤2所得最小嵌入深度g使用罚函数法计算法向接触力和摩擦力，根据步骤2所得接触区域X求得广义接触力以及广义接触力的雅克比矩阵，具体如下：

接触面上主薄壳单元中面上i点和从薄壳单元中面上j点的法向接触力如下：

F_n＝-pgn

式中p是罚参数，g是嵌入深度；在已知法向接触力的F_n情况下，n是主薄壳单元中面上点i和从薄壳单元中面上点j间的法向量，如下所示：

其中，Δr＝r(ξ_i,η_i)-r(ξ_j,η_j)，Δr为从i点指向j点的向量；

由虚功原理得到法向接触力的广义力形式：

其中δW_n是法向接触力F_n的虚功，

其中e_i、e_i分别为i接触点所在薄壳单元的节点坐标和j接触点所在薄壳单元的节点坐标；P_n是两薄壳单元上P接触点和Q接触点的广义法向接触力，满足：

其中，中间变量

广义法向接触力P_n的雅可比矩阵J_n为：

对于切向接触力即库伦摩擦力F_t而言，摩擦力为：

其中μ_m为摩擦系数，v_t为两薄壳单元上P接触点和Q接触点之间的相对速度；

切向接触力F_t的广义力P_t可由虚功原理得到，即两接触点之间的广义切向接触力P_t为：

广义切向接触力P_t的雅克比矩阵J_t为：

广义法向接触力P的雅克比矩阵J为：

5.根据权利要求1所述的基于ANCF的薄壳碰撞接触计算方法，其特征在于：步骤4中，使用广义-α算法作为柔性多体系统动力学隐式求解方法，将步骤1所得薄壳单元弹性力和弹性力的雅克比矩阵，和步骤3所得的广义接触力以及广义接触力的雅克比矩阵代入多体系统动力学方程，进行数值计算，获得动力学方程的雅克比矩阵和每一时间步长下柔性多体系统的节点坐标向量，具体如下：

所述柔性多体系统由两个薄壳构成，基于绝对节点坐标方法对含约束的柔性多体系统动力学而言，其动力学方程如下：