基于DeepCFD模型和CNN/U-Net模型的流场预测

1.遇到问题

计算流体力学（Computational fluid dynamics, CFD）通过对Navier-Stokes方程（简称N-S方程）的精确求解，能够精准获取流体在不同状态下的物理量分布详情，这些物理量包括但不限于密度、压力及速度等关键参数。以下图为例，展示了圆柱体周边二维非均匀稳定层流通道流动的模拟结果。

在某些复杂的应用场景中，如机翼优化和流体与结构相互作用方面，需要使用千万级甚至上亿的网格对问题进行建模，导致CFD的计算量非常巨大。因此，目前亟需发展出一种相比于传统CFD方法更高效，且可以保持计算精度的方法。

因此，使用数据驱动的机器学习方法，使用一小部分资源来生成这些模拟的精确近似是非常有吸引的。这些近似解决方案以低错误率为代价加速结果的潜力可以更有效地开发依赖 CFD 研究的产品。问题来了，提出DeepCFD网络模型，该如何设计网络模型来解决流体力学中流场的问题呢？数据集应该如何构造呢？

2.前情提要

1）流场问题

二维不可压缩流体的 Navier-Stokes方程：

假设非均匀稳态流动条件：

输入：障碍物、边界条件

输出：流体的 x 方向速度、y 方向速度和流体压强 p

障碍物：

基本形状(圆形、正方形、前向三角形、后向三角形和菱形)，如 primitive。

生成器生成障碍物形状为sample1 ~ 3的 2D 河道流量数据集的障碍物。

边界条件：

如下图所示：

顺流方向的域尺寸为260 mm，垂直于流动方向的域尺寸为120 mm。

网格元素的数量根据所使用的形状而变化，但是对于大约 1 mm 的基本单元大小，平均单元数大约为30000。

边界条件保持固定，入口(左壁)的恒定径向速度为0.1 米/秒，出口(右壁)的梯度为零，顶部/底部和障碍物壁的边界条件无滑 ，层流动态粘度设置为 1 * 10^-4 m^2/s，中心差分方案(CDS)用于动量方程对流项和扩散项的离散化。

2）数据集

数据集使用原作者利用OpenFOAM计算的CFD算例，共981组，分为两个文件（dataX.pkl, dataY.pkl），两个文件大小都是152 MB，形状均为[981, 3, 172, 79]。dataX.pkl包括三种输入：障碍物的SDF、计算域边界的SDF和流动区域的标签；dataY.pkl包括三种输出：流体的x方向速度、y方向速度和流体压强。数据获取使用的计算网格为172×79。

数据集地址：Data_DeepCFD_数据集-飞桨AI Studio星河社区

或https://www.dropbox.com/s/kg0uxjnbhv390jv/Data_DeepCFD.7z?dl=0

数据效果图如下：

3)深度学习模型

Net网络：CV领域中图像分割的算法，其主要特点是输入和输出的尺寸是一样的。

典型的U-Net网络如下图所示：

由于其网络结构“U”型而得名，是由卷积、下采样、上采样和拼接操作组成的编码器—解码器对称网络。

网络左半部分是输入路径，右半部分是输出路径。

3.解决方案

1)网络结构

网络结构：编译器-解码器，如下图：

a图是具有SDF和流动区域多类别标记的输入通道。b图是下采样卷积操作从输入中创建流几何的潜在表示。c图是上采样反卷积将LGR显示出特定的变量。

2)DeepCFD模型

DeepCFD U-Net网络：模仿U-Net网络，由卷积、下采样、上采样和拼接操作组成。

其中，实验中4种变形模型如下：

3)损失函数

其中，速度v是L2，压强是L1。

4)训练细节

最后选择是：kernel size = 5。

4.代码结构及参数说明

1）自定义代码结构

loss_func()：定义了计算一个模型（model）和一组批次数据（batch）的损失的损失函数。

    def loss_func(model, batch):
        x, y = batch
        output = model(x)
        lossu = ((output[:,0,:,:] - y[:,0,:,:]) ** 
2).reshape((output.shape[0],1,output.shape[2],output.shape[3])) 
        lossv = ((output[:,1,:,:] - y[:,1,:,:]) ** 
2).reshape((output.shape[0],1,output.shape[2],output.shape[3])) 
        lossp = torch.abs((output[:,2,:,:] - 
y[:,2,:,:])).reshape((output.shape[0],1,output.shape[2],output.shape[3])) 
        loss = (lossu + lossv + lossp)/channels_weights
        return torch.sum(loss), output

split_tensors()：将数据集（x，y）拆分成一定百分比的数据集和测试集

def split_tensors(*tensors, ratio):
    assert len(tensors) > 0
    split1, split2 = [], []
    count = len(tensors[0])。
    for tensor in tensors:
        assert len(tensor) == count
        split1.append(tensor[:int(len(tensor) * ratio)])
        split2.append(tensor[int(len(tensor) * ratio):])
    if len(tensors) == 1:        split1, split2 = split1[0], split2[0]
    return split1, split2

下方代码是一个用于在给定数据加载器上执行一个 epoch 的训练或者是验证的函数。

def epoch(scope, loader, on_batch=None, training=False):
    model = scope["model"]
    optimizer = scope["optimizer"]
    loss_func = scope["loss_func"]
    metrics_def = scope["metrics_def"]
    scope = copy.copy(scope)
    scope["loader"] = loader
    metrics_list = generate_metrics_list(metrics_def)
    total_loss = 0
    if training:
        model.train()
    else:
        model.eval()
    for tensors in loader:
        if "process_batch" in scope and scope["process_batch"] is not None:
            tensors = scope["process_batch"](tensors)
        if "device" in scope and scope["device"] is not None:
            tensors = [tensor.to(scope["device"]) for tensor in tensors]
        loss, output = loss_func(model, tensors)
        if training:
            optimizer.zero_grad()
            loss.backward()
            optimizer.step()
        total_loss += loss.item()
        scope["batch"] = tensors
        scope["loss"] = loss
        scope["output"] = output
        scope["batch_metrics"] = {}
        for name, metric in metrics_def.items():
            value = metric["on_batch"](scope)
            scope["batch_metrics"][name] = value
            metrics_list[name].append(value)
        if on_batch is not None:
            on_batch(scope)
    scope["metrics_list"] = metrics_list
    metrics = {}
    for name in metrics_def.keys():
        scope["list"] = scope["metrics_list"][name]
        metrics[name] = metrics_def[name]["on_epoch"](scope)
    return total_loss, metrics

2）参数说明

参数	推荐值	额外说明
batch_size	64	批次大小，默认64
train_test_ratio	0.7	训练集占数据集的比例，0.7即训练集70%测试集30%
learning_rate	0.001	学习率
weight_decay	0.005	AdamW专用，若修改优化算法需要修改train.py
epochs	1000	训练轮数
kernel_size	5	卷积核大小
filters	8, 16, 32, 32	卷积层channel数目
batch_norm	0	批量正则化，0为False，1为True
weight_norm	0	权重正则化，0为False，1为True
data_path	./data	数据集路径，视具体情况设置
save_path	./result	模型和训练记录的保存路径，视具体情况设置
model_name	DeepCFD_965.pdparams	具体加载的模型名称，后缀不能省略

3）部分训练日志如下：

Epoch #1

        Train Loss = 884808909.0

        Train Total MSE = 10197.3000353043

        Train Ux MSE = 3405.3426083044824

        Train Uy MSE = 4334.0962839376825

        Train p MSE = 2457.8616943359375

        Validation Loss = 53205074.5

        Validation Total MSE = 1027.7523040254237

        Validation Ux MSE = 419.7688029661017

        Validation Uy MSE = 543.9674920550848

        Validation p MSE = 64.01604872881356

Epoch #2

        Train Loss = 75408434.25

        Train Total MSE = 603.198411591199

        Train Ux MSE = 277.9321616481414

        Train Uy MSE = 303.4222437021684

        Train p MSE = 21.843986488987337

        Validation Loss = 17892356.5

        Validation Total MSE = 312.7194186970339

        Validation Ux MSE = 169.64230501853814

        Validation Uy MSE = 140.46789757680085

        Validation p MSE = 2.6092084981627384

5.实验对比

1）误差对比

2）计算效果对比

3）效果展示

1）ground-truth CFD (simpleFOAM) 和 DeepCFD 预测之间的比较，显示了速度分量和压力场，以及基于圆的形状 1 周围流动的绝对误差。

2）ground-truth CFD (simpleFOAM) 和 DeepCFD 预测的比较，显示了速度分量和压力场，以及围绕基于方形的形状2的流动绝对误差。

3）ground-truth CFD (simpleFOAM) 和 DeepCFD 预测的比较，显示了速度分量和压力场，以及基于菱形的形状 1 周围流动的绝对误差。