火箭发动机喷流的“监察队长”：基于飞桨探索火箭发动机真空羽流流场的快速计算

真空环境中，火箭发动机喷流向外部环境自由膨胀形成羽毛状流场，称为真空羽流[1]。真空羽流对航天器产气动力、气动热、污染、电磁干扰和视场干扰等效应统称为羽流效应。羽流效应会干扰航天器正常工作状态，甚至影响航天器寿命和任务成败。因此，真空羽流及其效应评估和防护是航天领域的重要科学和工程问题。

本文作者：张百一，北京航空航天大学博士研究生、飞桨开发者技术专家（PPDE）

项目意义

当前，我国正在开展载人登月关键技术攻关，真空羽流与月尘相互作用评估及控制是其中一项关键技术。真空羽流与月面作用会激起月尘，引发仪器读数错误和视野遮挡等问题，因此实时预测羽流场及月尘可为保障航天员/航天器安全及任务成功提供关键数据支撑[2]。但受限于计算效率，当前的主流数值模拟方法——直接模拟蒙特卡洛（Direct Simulation Monte Carlo, DSMC）方法无法实现实时预测。由此可见，大幅提高真空羽流数值模拟效率是当前亟待解决的问题。‍

图1 嫦娥五号落月过程影像

图2 火箭发动机羽流和月尘的相互作用

为了解决这一问题，本项目基于飞桨框架提出了一种基于卷积神经网络的直接模拟蒙特卡洛方法[3]（Convolutional Neural Networks-based Direct Simulation Monte Carlo Method, CNN-DSMC）。该方法使用卷积神经网络从DSMC数值模拟结果中提取特征，以此训练真空羽流智能化计算模型，然后基于DSMC数值模拟中不同的几何拓扑和边界条件，输入计算模型，可以实现火箭发动机真空羽流的快速计算。

技术方案

本文针对月面探测器月面着陆过程，分别通过CNN-DSMC和DSMC方法实现月面探测器在不同悬停高度时的真空羽流流场计算。

问题描述

本文研究的月面着陆过程羽流仿真所采用的计算域如下图所示。除了月面，计算域的其余五个边界均设置为开放边界。月面和月球探测器表面设置为固体边界，热适应系数设置为1.0。所有边界的温度设置为300K，粒子与边界相互作用模型使用Maxwell模型。DSMC数值模拟的时间步长设置为10-7s 。

图3 月面着陆过程羽流仿真计算域示意图

计算流程

在CNN-DSMC方法中，计算分为数据预处理和模型训练两个过程。在数据预处理中，真空羽流仿真模型中的几何拓扑信息被抽象为符号距离函数（Signed Distance Function, SDF），边界条件信息被抽象为标识符矩阵（Identifier Matrix, IM）。SDF和IM共同作为训练集的输入，由DSMC数值模拟得到的真空羽流速度场（三个方向）和密度场作为训练集的输出，测试集为未参与训练的DSMC数值模拟算例，用于验证CNN-DSMC方法的准确性。在完成训练之后，就得到了真空羽流智能计算模型。

图4 CNN-DSMC方法的计算总体流程

下面是基于飞桨构建的训练函数：

• 数据集训练与测试
  

def epoch(loader, training=False):
    total_loss = 0
    if training:
        model.train()
    else:
        model.eval()

•  使用GPU

    with paddle.static.device_guard('gpu'):

• 遍历数据集

 for tensors in loader:
            loss, output = loss_func(model, tensors)
            if training:
                optimizer.clear_grad()
                loss.backward()
                optimizer.step()
            total_loss += loss.item()

    return total_loss

• 训练主循环

def train():
    for epoch_id in range(1, epochs + 1):
        begin = time.time()
        print("Epoch #" + str(epoch_id))

•  训练过程

 train_loss = epoch(train_loader, training=True)
        print("\tTrain Loss = " + str(train_loss))

•  测试过程

with paddle.no_grad():
            val_loss = epoch(test_loader, training=False)
        print("\tValidation Loss = " + str(val_loss))

• 输出一个epoch的时间

   if (epoch_id != 0):
            print("运行1个epochs的时间为{:.2f} s".format(time.time() - begin))

• 保存模型

   if (epoch_id % 2 == 0):
            paddle.save(model.state_dict(), os.path.join(save_path, "CNN-DSMC" + str(epoch_id) + ".pdparams"))
            print("Model saved!")

下面是使用飞桨框架训练的主代码。其中，优化器采用了飞桨提供的AdamW优化器，这个优化器具有权重衰减功能，相比Adam可以有效地提高模型的泛化性能：

import os
import pickle
from utils import *
from paddle.io import TensorDataset, DataLoader
from net import CNN_DSMC
import time

• 创建保存路径

save_path = r"./Run"
if not os.path.exists(save_path):
    os.makedirs(save_path)

• 加载原始数据

x = pickle.load(open("./data/dataX.pkl", "rb"))
y = pickle.load(open("./data/dataY.pkl", "rb"))
x = paddle.to_tensor(x, dtype="float32")
y = paddle.to_tensor(y, dtype="float32")

• 分割数据集

train_data, test_data = split_tensors(x, y, ratio=10 / 12)
train_dataset, test_dataset = TensorDataset([train_data[0], train_data[1]]), TensorDataset(
    [test_data[0], test_data[1]])

• 超参数设置

lr = 0.001
epochs = 20000
batch_size = 2

• 网络设置

kernel_size = 5
filters = [8, 16, 32, 32, 64, 64, 128]
model = CNN_DSMC(2, 4, filters=filters, kernel_size=kernel_size)
wd = 0.005

• 优化器设置

optimizer = paddle.optimizer.AdamW(learning_rate=lr, parameters=model.parameters(), weight_decay=wd)

• 加载数据集

train_loader = DataLoader(train_dataset, batch_size=batch_size, shuffle=True)
test_loader = DataLoader(test_dataset, batch_size=batch_size, shuffle=False)

• 开始训练

train()

模型介绍

CNN-DSMC方法中使用的整体网络结构如下图所示。该网络由1个编码器和2个解码器组成，其中每个编码器（解码器）都由7个（反）卷积块组成。每个（反）卷积块的结构组成包括3个（反）卷积层和1个最大（反）池化层。单个（反）卷积层包括（反）卷积、激活函数和批量正则化三个过程。

图5 CNN-DSMC使用的卷积神经网络结构

使用飞桨进行模型构建代码如下，其中，特别使用到了飞桨中的Conv3D、 Conv3DTranspose、max_pool3d和max_unpool3d用来构成模型中的Autoencoder结构：

class CNN_DSMC(nn.Layer):
    def __init__(self, in_channels, out_channels, kernel_size=3, filters=[16, 32, 64], layers=3, weight_norm=False, batch_norm=True, activation=nn.ReLU, final_activation=None):
        super().__init__()
        assert len(filters) > 0
        self.final_activation = final_activation
        self.encoder = create_encoder(in_channels, filters, kernel_size, weight_norm, batch_norm, activation, layers)
        decoders = []
        for i in range(out_channels):
            decoders.append(create_decoder(1, filters, kernel_size, weight_norm, batch_norm, activation, layers))
        self.decoders = nn.Sequential(*decoders)

•   定义编码器

def encode(self, x):
        tensors = []
        indices = []
        sizes = []
        for encoder in self.encoder:
            x = encoder(x)
            sizes.append(x.shape)
            tensors.append(x)
            x, ind = F.max_pool3d(x, 2, 2, return_mask=True)
            indices.append(ind)
        return x, tensors, indices, sizes

•  定义解码器

 def decode(self, _x, _tensors, _indices, _sizes):
        y = []
        for _decoder in self.decoders:
            x = _x
            tensors = _tensors[:]
            indices = _indices[:]
            sizes = _sizes[:]
            for decoder in _decoder:
                tensor = tensors.pop()
                size = sizes.pop()
                ind = indices.pop()
                x = F.max_unpool3d(x, ind, 2, 2, output_size=size)
                x = paddle.concat([tensor, x], axis=1)
                x = decoder(x)
            y.append(x)
        return paddle.concat(y, axis=1)

    def forward(self, x):
        x, tensors, indices, sizes = self.encode(x)
        x = self.decode(x, tensors, indices, sizes)
        if self.final_activation is not None:
            x = self.final_activation(x)
        return x

数据处理

在CNN-DSMC方法中，真空羽流仿真模型中的几何拓扑信息被抽象为SDF，定义Ω为度量空间X的一个子空间，SDF定义如下：

式中：d为欧几里得距离，∂Ω为Ω的边界，Ω^c为Ω的内部。

SDF的物理意义为空间中的某个点到边界的最小距离，其符号由该点是否在边界内决定。如下图所示，月球探测器内部区域的SDF为负值，月球探测器外部的SDF为正值，且值的大小随着接近月球探测器边界而减小。

图6 代表几何拓扑信息的符号距离函数（SDF）

在CNN-DSMC中，边界条件信息被抽象为IM。IM本质上是一个三维的矩阵，其矩阵元素作为区分三维空间不同区域的标识符。在本文中，共选取4种不同的标识符：开放边界、航天器边界、月面边界和真空羽流区域，这4种标识符也与实际DSMC数值模拟中使用的边界条件相对应。具体设置如下图所示。

图7 代表边界条件信息的标识符矩阵（IM）

结果展示

将DSMC计算得到的数据输入CNN-DSMC网络中进行训练，共训练40000步，利用如下代码进行测试：

•  模型路径

path = r"./Run/CNN-DSMC.pdparams"

• 初步设置网络

kernel_size = 5
filters = [8, 16, 32, 32, 64, 64, 128]
net = CNN_DSMC(2, 4, filters=filters, kernel_size=kernel_size)

• 加载模型参数

net.set_state_dict(paddle.load(path))

• 加载数据

x = pickle.load(open("./data/dataX.pkl", "rb"))
y = pickle.load(open("./data/dataY.pkl", "rb"))
x = paddle.to_tensor(x, dtype="float32")
y = paddle.to_tensor(y, dtype="float32")

train_data, test_data = split_tensors(x, y, ratio=10 / 12)
train_dataset, test_dataset = TensorDataset([train_data[0], train_data[1]]), TensorDataset(
    [test_data[0], test_data[1]])
test_x, test_y = test_dataset[:]

• 模型运行

out = net(test_x)

• 读取流场坐标
  

coord = pickle.load(open("./coord.pkl", "rb"))
XX = np.unique(coord["x"])
YY = np.unique(coord["y"])
ZZ = np.unique(coord["z"])

• 导出流场预测数据，并转换为国际制单位

s = 0
rho = 10 ** (-10 * out[s, 3, :, :, :]) - 1e-10
u = out[s, 0, :, :, :] * (-5010)
v = out[s, 1, :, :, :] * (-5010)
w = out[s, 2, :, :, :] * (-5010)
ux_pre = u.detach().numpy()
uy_pre = v.detach().numpy()
uz_pre = w.detach().numpy()
rho_pre = rho.detach().numpy()
ux_pre = np.flip(ux_pre)
uy_pre = np.flip(uy_pre)
uz_pre = np.flip(uz_pre)
rho_pre = np.flip(rho_pre)

• 写数据为Tecplot格式

with open("result_pre.dat", 'w') as f:
    f.writelines("variables=\"X\",\"Y\",\"Z\",\"ux\",\"uy\",\"uz\",\"rho\"\n")
    f.writelines("zone i=300 j=100 k=100\n")
    for idz, Z in enumerate(ZZ):
        for idy, Y in enumerate(YY):
            for idx, X in enumerate(XX):
                f.writelines(    "{:.6e}\t{:.6e}\t{:.6e}\t{:.6e}\t{:.6e}\t{:.6e}\t{:.6e}\n".format(X, Y, Z, ux_pre[idx, idy, idz], uy_pre[idx, idy, idz], uz_pre[idx, idy, idz], rho_pre[idx, idy, idz]))

首先得到了h=8.2m情况下的真空羽流速度场的DSMC数值模拟结果和CNN-DSMC计算结果，如下图所示。可以看出，两者得到的速度场几乎完全一致，CNN-DSMC计算的激波形状也与DSMC数值模拟结果一致。图中也给出了流线的分布，结果表明DSMC数值模拟结果和CNN-DSMC计算的粒子运动轨迹也基本相同。

 图8 真空羽流速度场的DSMC数值模拟结果（左）

CNN-DSMC计算结果（右）

下图分别给出了月球探测器轴线处速度和密度的变化曲线，范围为-9m到-1 m，其中，-9m对应于月面位置，-1m对应于月球探测器正下方与足垫同一高度的位置。结果表明，CNN-DSMC计算得到的速度和密度与DSMC数值模拟结果基本一致，其平均相对误差分别为6.0%和8.8%。

图9 月球探测器轴线处的速度（左）和密度（右）分布

然后对CNN-DSMC的计算效率进行评估。下表给出了三维情况下的DSMC和CNN-DSMC计算时间对比。其中，DSMC数值模拟的计算总时间取的是训练集中11个算例运行的平均时间。结果显示，随着算例数目的增加，平均单个算例时间从3.1906s降低到了0.2141s，相比传统DSMC方法的加速比可达4-6个量级，极大地提高了真空羽流的计算效率。卓越的加速性能和较高的流场计算精度表明，CNN-DSMC是一个非常有应用潜力的真空羽流智能化计算方法。

表1 CNN-DSMC和DSMC方法计算时间对比

最后，对利用CNN-DSMC从训练到推理整个的步骤进行简要介绍。考虑到大部分人可能没有充足的计算资源，这里可以使用AI Studio完成整个过程。AI Studio可以提供算力支持，可免费使用V100、A100等显卡进行训练，满足绝大多数网络的训练要求。

• 模型训练

python
git clone https://github.com/X4Science/CNN_DSMC.git # clone the repo
cd CNN_DSMC
chmod 777 main.py
python main.py # 开始训练，训练前请先下载数据

• 模型推理

python
python ./database.py # 先在文件中设置好模型路径

推理完成后，得到*.dat文件，可使用流体可视化软件如Tecplot查看流场计算结果。

参考文献

[1] HE B, ZHANG J, CAI G. Research on vacuum plume and its effects[J]. Chinese Journal of Aeronautics, 2013, 26(1): 27-36. 

[2] HE B, HE X, ZHANG M, et al. Plume aerodynamic effects of cushion engine in lunar landing[J]. Chinese Journal of Aeronautics, 2013, 26(2): 269-278.

[3] 蔡国飙, 张百一, 贺碧蛟, 等. 真空羽流智能化计算研究[J]. 航空学报, 2022, 43(10): 527352.

拓展阅读

相关地址

• AI Studio链接
  

https://aistudio.baidu.com/aistudio/projectdetail/4486133

• Github repo

https://github.com/X4Science/CNN_DSMC/tree/main/CNN_DSMC

• 论文链接

https://hkxb.buaa.edu.cn/CN/10.7527/S1000-6893.2022.27352

关注【飞桨PaddlePaddle】公众号

获取更多技术内容~