温度冲击下粘弹性装药结构完整性快速计算方法①

邓奕铭，王彦涛，惠卫华

(1.北京机电工程研究所，北京 100074；2.郑州机电工程研究所，郑州 450015；3.西北工业大学 燃烧、热结构与内流场重点实验室，西安 710072)

0 引言

固体推进剂作为固体火箭发动机的重要构成部件之一，其力学性能将直接关系到发动机能否正常工作[1]。固体火箭发动机装药在温度冲击载荷下容易出现裂纹和脱粘，使发动机结构完整性遭到破坏，从而影响发动机的安全性和使用可靠性[2]。但是长期以来，温度冲击下的装药力学响应仿真分析过程复杂，三维建模与有限元计算软件前处理之间的数据传递过程费时费力、效率低下，给发动机装药设计工作带来了极大的不便[3]。而且温度冲击实验是每个型号的发动机在列装前必须进行的，如果能找到一种快速有效的方法来分析温度冲击装药结构响应，在发动机设计完成前验证装药的温度冲击性能，可有效提高发动机设计效率，避免时间浪费。在药柱结构完整性方面国内外学者有很多研究，H C Yildirim和S Ozupek以及Heller等[4-6]研究了低温贮存温度下固体推进剂的应力应变和累积损伤；邓康清等[7]通过热-机耦合仿真分析，分别模拟了固体火箭发动机在降至常温的慢速固化降温过程和降至低温的快速降温过程。在参数化建模方面，李磊等[8]利用MSC.Patran软件的二次开发工具PCL，实现了圆管伞盘型药柱、星型药柱、车轮型药柱的二维参数化建模分析。

上述固体火箭发动机参数化建模研究中对推进剂药柱进行参数化建模时可能会过多地简化药柱模型或忽略药柱的一些微小几何特征，无法全面地反映药柱的几何形状，而对装药结构响应的研究中对不均匀温度场对装药结构响应的影响也分析甚少。因此，研究温度冲击下的装药结构响应快速分析方法具有重要意义。其能简化繁琐的发动机装药设计仿真流程，提高设计、分析效率，降低试验、计算成本，而且能够实现特定工况下的分析计算功能。

1 装药参数化建模

1.1 装药几何分析及参数化

固体火箭发动机从初样设计到最终定型是一个不断重复的过程，要历经多轮设计，但每次修改仅是小范围的尺寸修改，如壳体厚度，药柱肉厚等，如果每次修改都要求设计人员重新开始，将极大地增加设计人员的负担，而固体发动机的参数化建模能很好地解决这一问题。

参数化建模[9-10]有基于特征、尺寸约束、关系驱动等特点，因此在参数化建模的过程中要注意以下要点：

(1)参数化几何模型要力求精确还原实际物体基本特征，如星型药柱要准确反映星根倒角、星尖弧度等基础特征；

(2)尺寸链完备且闭合，要有完备的尺寸链约束且不能过约束，过约束会导致参数传递错误，造成尺寸链无法闭合，从而引起参数无法修改的情况；

(3)选择适当的关系驱动，尽量选取关键参数舍弃次要参数。

参数化设计起始于草图，进而形成实体模型。利用草图进行相关参数化尺寸驱动，然后对草图上的对象进行位置关系约束。需要注意的是，建模时有且只能有一个主特征，其他的特征都利用它的基准点等进行定位，从而确定位置关系[11-13]。

药柱通常是由药柱基体和药柱芯模通过布尔差构造的。在对药柱型面进行修改时，往往药柱外轮廓不动而对内型面进行微小改动，这样就只需修改药柱芯模，再运用一次布尔差即可，可以避免药柱外型面和绝热层之间的联动而造成的错误。

1.1.1 药柱基体

药柱基体即为实心药柱。它的主要参数包括长度L、椭球封头的长轴L3、短轴L2、药柱外半径R等参数。典型的星型药柱外形轮廓如图1所示。药柱基体是典型的回转体，在得到药柱外轮廓后，将此外轮廓沿药柱中轴旋转360°后可得到三维药柱基体。

图1 药柱基体外轮廓尺寸

1.1.2 药柱芯模

星孔药柱的芯模可以通过星形草图拉伸得到，星孔参数较为固定，如表1所示。星孔的几何参数化轮廓及得到的实体芯模如图2所示。

1.1.3 药柱装配

在完成药柱基体和芯模的参数化建模后，将芯模端面和药柱基体端面进行匹配定位，最后采用布尔差运算完成星型药柱的建模。图3为运用该方法得到的星型药柱。

表1 星型内腔参数化方案

图2 二维星孔建模参数及星孔芯模

图3 星型药柱实体模型

1.1.4 装配体建模

在药柱、壳体、绝热层模型参数化完成后，为方便模型的修改，建立装配体参数，然后将各个参数关联到零件的参数上。这样就可以实现只需修改装配参数，零件参数就可自动修改，进而整个模型(包括装配体中所有零件)实现参数联动，即实现重生所有零件，并且还能保证所有零件配合依然正确。

经过参数与关系的设置，即可最终完成整个模型的参数联动，完整的装配体模型如图4所示。

图4 星型药参数化完成后的装配体

1.2 装药建模自动化实现

在完成装药的参数化建模后，运用Pro/E二次开发来实现装药模型的修改、再生、保存等，二次开发涉及到调用函数对应关系如图5所示。可以看出，在二次开发中比较重要功能有加载参数、修改参数、模型再生、导出模型等几个关键步骤，下面对着几个功能进行介绍[14-15]。

图5 Pro/E与J-Link函数对应使用流程

(1)加载参数

设计模板涉及多个装药参数模型，每个装药模型又涉及多个参数，在单个设计模板中Pro/E随机返回参数，存在随机性，这样程序不够稳健。为加强稳健性，采用“参数名-取值”的结构，返回列表中的数据也采用这种结构，可有效联动图形用户界面，方便用户操作，实现模型参数界面修改

(2)修改参数

在参数加载完成后，同一模板各参数不可能一直保持不变，否则就失去了参数化建模的意义，如果想对特定尺寸进行修改，这时就要用到Pro/E二次开发中的接口Parameter.SetValue()函数，该函数可将参数修改值进行传递。

(3)执行再生

在数据参数修改完成后，执行Pro/E的模型再生功能，读入修改后尺寸数据，进行重新造型。再生步骤利用Pro/E提供的Solid. Regenerate()函数来完成，该函数可对Solid类型的模型进行操作。

(4)模型导出

Pro/E生成的部件和装配体的默认格式为*.prt和*.asm，但本文对装药温度冲击的仿真计算时需要模型往往是*.x_t或*.step格式文件，这就需要模型导出时进行文件格式的转化，借助J-Link中的Model.Export()函数，即可将零件或装配体转化为所需文件类型。

1.3 温度冲击过程

根据GJB 5021—2001有关温度冲击实验的规定，温度冲击循环次数一般规定为3次，如图6所示。发动机温度冲击保温时间与保温厚度之间的关系如表2所示。保温时间指药柱内腔温度达到保温温度后，还需要继续保温的时间。

图6 温度冲击实验流程

本次仿真所用星形装药发动机的装药肉厚为60 mm，故根据温冲实验保温要求，选择保温时间为t≥24 h。

表2 保温时间与厚度关系

本文采用第三类热传导边界条件进行温度冲击仿真，由流体横掠圆管表面平均表面传热系数的关联式、相关物性参数和发动机几何参数可得：调温期间对流换热系数h=8.0 W/(m2·K)；高温冲击期间对流换热系数h=13.0 W/(m2·K)；低温冲击期间对流换热系数h=9.0 W/(m2·K)。

对应图6实验流程设定温度冲击装药的仿真流程如下，装药模型初始温度为20 ℃，然后对发动机壳体外表面施加8.0 W/(m2·K)的对流边界条件，持续48 h；紧接着高温冲击对发动机壳体外表面施加13 W/(m2·K)的对流边界条件，也保持48 h；最后进行低温冲击时对壳体外表面施加9.0 W/(m2·K)的对流边界条件，同样保持48 h，第一次高低温仿真冲击完成，然后按照此流程重复3次即可，这时温度冲击装药结构完整性仿真分析完成。

2 温度冲击装药结构响应计算模板

在进行典型装药模板的设计时，为保证模板的通用性和精准性，要建立一套标准的分析流程，用通用的方法完成网格的剖分和离散、边界条件的设置、载荷工况的加载等，在保证精确度的情况下，尽可能的减少人工干预和工作量，保持模板的合理性。

典型的温度冲击装药模板的建立流程如下：

(1)对目标装药的温度冲击过程进行详细分析，确定所建立的有限元模型的理想假设情况；(2)确定有限元模型的几何参数、材料参数、边界条件、载荷工况等，根据装药特征确定重点关注部位，如星型药柱的星尖位置；(3)确定装药的网格离散方法及通用建模方法；(4)在通用建模方法的基础上，提取*.rpy日志文件，并在提取的*.rpy日志文件的基础上添加相关命令流，实现参数化的定义；(5)编写GUI，使Pro/E的参数化模型和Python脚本得以联动，实现数据的内部传递及计算参数的灵活输入。

2.1 控制方程和物性参数

推进剂材料的热粘弹性本构方程为[15]：

(1)

式中t为具体某时刻；τ为时间变量；T0为参考温度；εT(τ)为热应变，可通过热膨胀系数α0和温度变化来表示；ζ(t)为时间折减因子；ζ(t)=t/αT，αT为温度-时间移位因子。

其他装药具体材料参数如表3所示。

该药柱为典型的粘弹性材料，采用大应变模型来描述，用Prony级数的形式来表征其松弛模量，如式(2)所示。

(2)

其前6阶参数如表4所示。

药柱材料的时间-温度等效因子的WLF方程为

(3)

式中C1，C2为材料常数，由材料自身特性决定，C1=17，C2=52；T0为基准温度，T0=20 ℃。

表3 装药材料参数

表4 药柱粘弹性材料参数

确定完装药各部件的材料参数后，以药柱材料的传导率为例，定义材料参数代码如下：

mdb.models[′Model-1′].materials[′grain′].Conductivity(table=daorexishu_yao)

对实体部件赋予不同材料属性：

p=mdb.models[′Model1′].parts[′GRAIN_XING′]

c=p.cells

region=p.Set(cells=cells, name=′Set-1′)

c1=mdb.models[′Model-1′].parts[′grain′] .cells

p.SectionAssignment(region=region, sectionName=′grain′,offset=0.0,

offsetType=MIDDLE_SURFACE, offsetField=′′,

thicknessAssignment=FROM_SECTION)

2.2 网格划分与边界条件

对于简单的几何体，可以采用扫掠法。但是设计药柱时，为内弹道的平稳性通常会有含有翼槽等几何结构，就不能直接生成全六面体网格，需要做适当的剖分，分别对剖分后的每一部分几何元素采用扫掠法生成六面体网格。

如翼柱形装药的网格划分的难点集中在翼槽的位置，这里几何形状较为复杂，有众多线段、倒角、圆弧等几何，需要将这一段单独剖分出来，与中间的圆柱段分离。如图7所示，利用Face1可将翼槽单独剖分出来，之后利用Face3沿着轴向剖分翼槽。需要说明的是，Face3与翼槽的拉伸面平行，且距离翼槽中心的面距离为Thick/2-3，Thick为翼宽，如图8所示。

图7 翼柱形药柱(1/2n)的几何剖分

图8 Face3示意图

其中Face1和Face2可以通过创建模型时由Pro/E输出面的x点坐标确定，Face3由偏移翼槽剖面得到。三个面的创建相关代码如下：

Face1=mdb.models[′Model-1′].parts[′grain′].DatumPlaneByPrincipalPlane(offset=offset1,principalPlane=YZPLANE)

Face2=mdb.models[′Model-1′].parts[′grain′].DatumPlaneByPrincipalPlane(offset=offset2,principalPlane=YZPLANE)

Face3=mdb.models[′Model-1′].parts[′grain′].DatumPlaneByOffset(flip=SIDE2,offset=(Wing_Thick*0.5-3),plane=plane1)

Offset1和Offset2由Pro/E输出面的x坐标来定位，plane1为翼槽的剖面。

通过上述方法的剖分，翼槽的两个几何元素就可采用扫掠法生成六面体网格，扫掠路径见图9。

图9 翼槽扫掠示意图

本文采用的通星型装药网格划分较为容易，从前端到后端即是扫掠路径，药柱的网格单元选择杂交公式的热力耦合C3D8HT单元。生成的网格模型如图10所示。

图10 药柱翼槽有限元网格

其实现命令如下：

p=mdb.models[′Model-1′].parts[′GRAIN_XING′]

p.seedPart(size=seedSize,deviationFactor=0.1,minSizeFactor=0.1)

elemType1 = mesh.ElemType(elemCode=C3D8HT,elemLibrary=STANDARD)

c = p.cells

pickedRegions =(c,)

p.setElementType(regions=pickedRegions,elemTypes=(elemType1)

生成的星型药柱网格如图11所示。

图11 星型装药网格

由于装药载荷和几何均具有对称性，计算通常只选取1/2n模型，为限制刚性位移，通常还需设置固支边界条件，为识别这些表面，可采用RGB编码来标识。其实现如下：

p1.SetFromColor("shell_duichen_set",(255,255,0))

p2.SetFromColor("jrc_duichen_set",(255,255,0))

p3.SetFromColor("grain_duichen_set",(255,255,0))

a=mdb.models[′Model-1′].rootAssembly

p1_in=a.instances[′XING_SHELL′]

p2_in=a.instances[′XING_JRC′]

p3_in=a.instances[′GRAIN_XING′]

a.SetByBoolean("duichen_set",(p1_in.sets[′shell_duichen_set′],

p2_in.sets[′jrc_duichen_set′],p3_in.sets[′grain_duichen_set′]))

由以上语句可知，对不同的表面进行颜色标识时，不同颜色的表面建立不同的集合，然后进行标识储存。

3 快速计算平台实现

3.1 界面设计思想

固体火箭发动机装药结构完整性快速计算分析软件包括了模型生成、材料参数输入、仿真计算及后处理等模块，如图12所示。

图12 温度冲击装药仿真分析平台

(1)图形用户界面：用于完成用户与程序间的数据传递，显示仿真结果。

(2)模型生成：用于输入模型几何参数，生成发动机模型。

(3)材料参数输入：用于输入装药各部件材料参数。

(4)仿真及后处理：用于提交inp文件给ABAQUS以及统计结构完整性仿真结果，包括应力、应变、位移。

软件采用模块化的方式建立了界面，软件分为4个模块，分别负责三维视图显示、连接Pro/E、处理Python脚本文件与界面可视化。采用模块化设计的优势不言而喻，模块之间相互独立，一个模块内部发生变化时，其他模块几乎不受影响。

3.2 GUI界面调用

利用JAVA GUI接口建立软件界面，实现软件可视化。这部分主要负责连接Pro/E与ABAQUS，并且提供与用户的交互功能。其工作流程大致如下：在后台连接Pro/E，打开模型模板，读取模板参数至界面供用户修改，在用户修改完毕后，重新生成模型，并转换为ABAQUS支持的格式；读取相应计算模块的ABAQUS脚本参数，用户修改后生成新的脚本；后台执行脚本，分析开始。

软件共有3个大的功能界面。每个大的功能界面底下有5个独立的小的功能，分别为模型几何参数修改、材料参数设置、网格与计算参数设置、结果云图查看、安全系数查看。

下文以材料参数设置和计算及后处理为例进行展示。

3.3 材料参数

设置材料参数是温度冲击仿真分析的关键步骤，精确完备的材料参数是完成高精度仿真分析的基础，如图13所示，该界面实现了药柱、绝热层、壳体的密度、比热、热传导率等材料参数的设置，其中绝热层与药柱的粘弹特性用Prony级数与时温等效因子来表征，简明清晰，可以实现材料参数的快速设定和修改。

图13 材料参数设置界面

3.4 计算及后处理

图14为求解计算界面，该界面完成温度冲击装药仿真的计算及后处理部分，在计算完成后可以对温度冲击分析的固化降温阶段、高温冲击及低温冲击阶段的温度和Mises应力及应变，无需进行繁琐的转换，用户可以自行选择所关注部件进行2D或者3D结果显示。

图14 计算与后处理

3.5 模型优化及收敛性分析

在建立通用建模方法时，应根据装药的形状和尺寸，以从底向上的方法，从低维图元向高维图元建立药柱几何模型。保证药柱模型中图元的数量和相对位置的固定，保证建模过程中图元的生成顺序固定，确保图元之间具有正确的拓扑约束关系。同时还需要对建模过程中用到的关键参数进行标记，固定其编号，避免因图元编号混乱而导致建模意外终止。为后续通过图元编号进行实体剖分和网格加密等步骤提供基础。

装药模型优化方法[16]框架中主要包含优化算法更新参数、模型重生及导出、温度冲击装药结构完整性计算。优化算法和模型重生及导出模块之间通过预留的参数配置文件作为数据交换的接口；模型重生与温度冲击装药计算模块之间采用颜色标识的方法传递初始信息，事先要给表面标记上特定的颜色，装药结构完整性计算运用模型中的颜色标识进行识别。

计算结果的收敛性主要受模型构造复杂度、网格划分精度、分析步类型、工况、边界条件设置五个方面，提高仿真分析的收敛性主要从以下几个方面入手：

(1)在进行几何建模时，可忽略药柱外表面的凸台等对结构完整性影响不大的细小几何特征，将药柱外表面看成光滑、连续的曲面，这样以来就无需考虑和凸台有关的几何模型参数，简化了建模过程。

(2)划分成精度更高的全六面体网格；采用映射法或扫掠法等高质量网格生成算法进行网格划分，对于复杂造型的药柱，需先分割成简单的子实体；对重点关注区域进行适当的网格加密操作。

(3)依据实际工作状态选择恰当的分析步类型，设置适当大小的分析步长，根据理论预估情况确定几何非线性的开关，可有效提高仿真分析的收敛性。

(4)定义工况时，确保工况参数的合理正确性，如温冲过程要准确设定高低温冲击时间、装药转移时间。初始温度、对流换热系数等，特别注意的是要保证各参数单位的一致性。

(5)边界条件设置既要合理又要简化，在进行装药温度冲击分析时，装药的边界条件是不变的，即在装药柱外表面上施加固支条件，约束外表面节点前三个方向的自由度，在药柱对称面上施加对称约束条件，约束对称面节点第二个方向的自由度。另外需要根据具体工况添加必要的边界条件，如温度冲击时需要添加温度边界条件。

通过以上五点可有效解决温度冲击装药仿真分析的收敛性问题。

3.6 快速性分析

在保证装药模型、材料参数、边界条件、网格单元等一致的情况下，在已验证手动计算结果相对正确的条件下，以该特定星型装药为例进行温度冲击装药结构完整性的快速性对比分析，各阶段用时对比结果如表5所示。

表5 原始手段与快速计算对比

在计算结果基本正确的前提下，从表5总的时间与手动计算的比对中可以看到，在计算精度一致的情况下，计算效率提升了54.03%。有效提高了温度冲击装药结构完整性的仿真分析效率。

4 结论

(1)通过分析装药参数化建模的核心思想和实现途径，基于J-Link的Pro/E二次开发技术，实现了药柱的修改、再生、保存、导出等一系列功能，建立起设计变量与药柱模型尺寸约束之间的映射关系，从而实现变量驱动药柱更改的功能，与计算模板做到了有效的数据传递。

(2)本章利用ABAQUS二次开发技术，与参数化建模模块形成了接口。建立了几种典型药型的结构完整性计算模板，通过修改模板内的参数，即可自动进行模型导入、网格离散、载荷及边界条件添加、有限元计算等功能。

(3)完成了温度冲击装药结构完整性的快速计算平台设计，整合了Pro/E的参数化建模、温度冲击计算模板，采用图形用户界面实现简明高效的温度冲击数值仿真，在保证计算精度情况下，最终实现温度冲击装药结构完整性的快速计算，能极大提高工作效率。