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

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

引用

CSDN

1.

https://blog.csdn.net/FastCAE/article/details/142333770

计算流体力学（CFD）通过对Navier-Stokes方程的精确求解，能够精准获取流体在不同状态下的物理量分布详情。然而，在某些复杂的应用场景中，如机翼优化和流体与结构相互作用方面，需要使用千万级甚至上亿的网格对问题进行建模，导致CFD的计算量非常巨大。因此，使用数据驱动的机器学习方法，使用一小部分资源来生成这些模拟的精确近似是非常有吸引力的。本文将介绍基于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）参数说明

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 周围流动的绝对误差。