整流罩地面分离过程仿真

1. 整流罩地面分离过程仿真

1.1 仿真动机

对于整流罩地面分离过程仿真的必要性，从以下几个方面考虑：

1）为了准确模拟整流罩在高空的分离过程，需要在大型真空罐内进行分离试验。然而，由于整流罩真空分离试验所需成本和对设备要求均较高，目前国内尚不具备实验条件。

2）随着我国航天水平的进步，对大推力火箭的需求日益迫切。为能承载更多有效载荷，大推力火箭要求配备大型整流罩。整流罩的尺寸增大使其结构刚度相对减小，地面环境下气动载荷带来的影响相对增大，可能导致出现如弯曲、扭转、呼吸等复杂的变形情况，不能视其为刚体。

1.2 仿真难点

整流罩地面分离过程仿真中存在的难点包括：

1）为提升仿真结果的精确性，要求整流罩模型具有一定精细度。省略一些不必要的结构或作适当简化，加强关键结构的准确性和精确度，如何掌握这两者的平衡是难点之一。

2）为实现仿真过程的高效性，必须考虑已有结构的网格划分。如何选取网格类型和网格密度，使之在不影响结果精确性的前提下实现较高的计算速度，又是一大难点。

3）整流罩地面分离是一典型的流固耦合问题。在此过程中，罩外大气附加整流罩气动力，影响其运动和变形；整流罩的运动与变形反过来又会影响罩外大气的流动。如何选择具有解决流固耦合问题能力的大型有限元软件以及如何选取合理高效的算法，成为整流罩地面分离过程仿真的首要难点。

1.3 仿真路线

以上，在仿真动机驱动下，重点关注仿真难点，同时考虑甲方对于大型有限元软件的要求，选择基于Abaqus/Explicit求解器的耦合欧拉-拉格朗日（CEL，coupled Eulerian-Lagrangian）算法分析整流罩在气动力作用下的分离特性，从而实现地面分离过程仿真。

1.4 耦合欧拉-拉格朗日算法优势分析

1.4.1 拉格朗日算法和欧拉算法

拉格朗日算法常用于固体力学中的受力与变形分析。其重要特征是有限元网格固连于材料区域且两者共享边界，所形成的拉格朗日单元内充满材料。因此，结构变形一致反映于有限元网格的变化，可跟踪节点的运动，从而简化控制方程的求解过程。但在大变形情况下，网格发生严重畸变，此时拉格朗日算法丧失了其准确性。

与拉格朗日算法相比，欧拉算法可有效应用于大变形问题，如液体晃动、气体流动以及渗流等。有限元网格固定于空间，其形状、大小、位置不随结构变形而变化。一般地，欧拉算法所形成的欧拉单元难以被同种材料填满或者无任何材料。因此，难以准确描述结构的材料边界。图1对比了拉格朗日算法和欧拉算法的单元特性。

a）拉格朗日算法 b）欧拉算法

图1 拉格朗日算法和欧拉算法的单元特性

1.4.2 耦合欧拉-拉格朗日算法

1.4.2.1 概述

耦合欧拉-拉格朗日算法由学者Noh提出，最初应用于带有移动边界的二维流体动力学问题。CEL算法吸收了拉格朗日算法和欧拉算法的优点并克服了两者的缺点。对固体建立拉格朗日模型，划分拉格朗日网格；对流体建立欧拉模型，划分欧拉网格。两类网格重叠处是耦合区，能够高效传递计算中的信息。

1.4.2.2 理论基础

1）控制方程由质量守恒、动量守恒、能量守恒及连续性方程组成。本次仿真的欧拉材料(即罩外大气）由本构方程与连续性方程描述，分别为：

(1)

其中：—应变张量

—大气压强

—剪切粘度

—应变变化率

(2)

其中：—大气压强

—标准大气压

—大气密度

—气体常数

—大气温度

—绝对零度

2）有限差分形式的时间积分。当前解由前一步获得而不需迭代。其数值稳定性由临界时间步判断

(3)

其中：—单元特征尺寸

—波速

3）罚函数耦合系数追踪拉格朗日节点和欧拉材料界面上锚定点的相对位移，检查每一个节点对欧拉物质界面的贯穿，如果节点不出现贯穿，则无需进行处理；如果发生节点对物质界面的贯穿，界面力就会分布到欧拉材料的锚定点上，如图2所示。界面力的大小与发生的贯穿相对位移成正比，即

(4)

其中：—罚刚度系数，仅取决于材料特性。

图2 罚函数耦合算法

1.4.2.3 研究现状

CEL算法在解决流固耦合问题时具有强大的优势，目前已有广泛的应用。左辉等人对比了简单气体压缩过程的仿真结果和解析结果，成功证明了Abaqus中理想气体状态方程的可靠性与正确性，并应用于产品分析。姚小虎等人采用CEL算法对水陆两栖飞机水上降落的流固耦合问题进行了结构响应分析，得到的高精度结果为结构设计提供了十分有效的参考。徐文杰采用CEL算法模拟分析了滑坡涌浪这一复杂的流固耦合问题，模拟结果与试验结果的高吻合度证明了CEL算法的巨大优越性。宋庆华利用ABAQUS的CEL功能实现燃油箱晃动的流固耦合分析，为产品设计提供了宝贵的参考。

1.4.2.4 Abaqus应用

Abaqus/Explicit求解器包含了CEL算法，且具有如下特性与优势：

1）高速计算。Abaqus/Explicit求解器对CEL算法进行优化，即使计算资源有限，仍能以较快速度运行，保证了计算的稳定性。

2）高度集成。将流体前处理与典型的固体前处理一同集成于Abaqus/CAE模块中，消除调用其它模块甚至其他软件的麻烦。

3）运动流体网格。于Abaqus 6.10引入的运动流体网格仅需用户在初始分析步内包含固体结构，而非整个运动轨迹，极大减少了计算成本。

4）欧拉体积分数（EVF，Eulerian volume fraction）。基于流体体积法，用户可定义每个欧拉单元内的欧拉体积分数，若材料充满某一单元，则体积分数为1；若单元内无材料，则体积分数为0。欧拉体积分数的定义极大方便了复杂欧拉材料的建模。

5）通用接触设置。用户只需在通用接触中选择接触类型，即可定义流体材料与结构几何边界的欧拉-拉格朗日接触。

1.4.3 光滑粒子流体动力学（SPH，smoothed particle hydrodynamic）算法

光滑粒子流体动力学是一种无网格数值方法，用大量离散的光滑粒子的集合定义流体材料，从而通过插值离散化连续方程组。其理论核心为核函数，实质是一定光滑长度范围内其他临近粒子对研究粒子影响程度的权函数，如图3所示。

图3 核函数

SPH可有效解决流体流动和自由表面等问题，但存在以下不足：

1）属于传统的拉格朗日方法，在处理大变形问题时结果准确性不如CEL。

2）其计算效果依赖于对于光滑粒子的建模，若粒子布置不均匀，可能会造成流体材料质量分布不均。

3）适合单线程计算，当粒子数量较多时对多线程计算有限制。

4）Abaqus使用SPH的步骤较为繁琐，首先在Abaqus/CAE中创建质量单元，再写入输入文件。

1.5 Abaqus/Explicit优势分析

有限元仿真领域对流固耦合问题一般采用计算流体动力学（CFD，Computational Fluid Dynamics），目前较为主流的求解器包括CFX、Fluent、COMSOL。这些传统求解器相对于Abaqus/Explicit的劣势包括：

1）以流体为核心，流体材料产生运动或变形后再将信息传递给固体结构，与整流罩地面分离过程截然相反，难以满足初始时刻的情况。

2）不同程度简化了气动载荷，对结果精度产生较大影响。落䶮寿曾利用Fluent进行整流罩地面分离试验仿真；张小伟等分析了整流罩分离过程的模态。两学者基于相似的简化原则，未考虑气动力的时变影响和空间分布不均匀性。

3）对固体材料仿真的求解能力较差。由于流固耦合分析对边界条件、模型尺寸、单元大小等因素极为敏感，要求同其他求解器联合求解，不利于仿真的二次开发。

1.6 技术方案

1.6.1 概述

图4为整流罩地面分离仿真的基本步骤。其中，固体结构的建模如前所述，在此不赘述。下文将重点阐述Abaqus CAE中的流体材料建模、装配部件、设置相互作用及施加载荷环节，为二次开发作铺垫。

图4 整流罩地面分离仿真的基本步骤

1.6.2 关键步骤

1） 流体材料建模。

流体材料为罩外大气，在创建部件时命名为“atmosphere”，类型为三维空间下的欧拉区域。由于整流罩半罩为半圆柱，宜采用旋转方式建模（见图5）。考虑到二次开发的实现，应简化建模操作，不妨设为半圆柱并进行partition操作（见图6）。为创建重叠网格区域并设置耦合欧拉-拉格朗日接触，两模型应有重合区域，故取分割面位置由内靠近整流罩蒙皮内表面。

图5 图6

        采用理想气体状态方程（2）式描述流体材料的属性

(2)

选择Abaqus材料库EOS模型中的Ideal Gas类型（见图7）。其中：Gas Constant为气体常数，取为；Ambient Pressure为环境压强，取为。此外还需设置气体密度为（见图8），动力粘度系数为（见图9），和绝对零度为-273（见图10）。

图7

图8

图9

图10

    选择EC3D8R型欧拉单元为其划分网格（见图11）。

图11

2） 装配部件

如前所述，调整部件相对位置，使固体结构与流体材料有重合部分，重合部分为整流罩蒙皮。

3） 设置相互作用

对接触性质进行编辑（见图12），选择罚函数定义切向行为，其中Friction Coeff为罚刚度系数，仅取决于材料特性。在通用接触中定义相互作用，赋予接触性质，关键是定义相互接触面（见图13）。对于固体结构，接触面为整流罩半罩蒙皮外表面，类型为几何；对于流体材料，接触面为整体，类型为网格（见图14）。

图12

图13

图14

4） 施加载荷

仅为整流罩固体结构施加重力场。创建预定义场（见图15），为已进行partition操作的罩外大气进行材料设置，罩外部分设为1，罩内部分使用缺省值0（见图16）。

图15 图16