基于图像数据生成网格的方法

你有没有想过利用三维成像技术获得的数据生成仿真网格？这篇文章，我们将解释如何使用 COMSOL Multiphysics® 软件来实现这一目标。这个主题是对我们在以前的文章中讨论过的模拟不规则形状主题的扩展，今天我们将要讨论的工作流程，一部分已经在验证拓扑优化结果研究时使用过。这篇文章，我们将对这个工作流程进行概括，使其适合用于通过三维成像技术获得的数据。

COMSOL Multiphysics® 支持导入不同格式的数据：txt 文件、Excel 文件、CSV 文件、图像、DEM 文件、MATLAB® 函数，以及用 C 语言编写的外部函数。扫描的数据通常以一组截面图的形式出现，包括我们在这篇文章中介绍的例子：人体股骨的一部分(这个模型中的某些数据来自 Visible Human 项目）。通常，COMSOL Multiphysics 只支持导入一张图像用于生成二维数据。因此，在使用软件之前，通过以下步骤生成数据文件：手动挑选带有股骨上半部分的图像，人工识别包含其中一根股骨的像素，以减少数据量。随后，识别出的像素亮度，连同它在空间中的位置将一起被读取和写入一个txt 文件。图像中的像素位置用 x 坐标和 y 坐标表示，图像编号用 z 坐标表示。电子表格是一种适合使用的格式，因为每个像素都是从图像中读取并逐一写入 txt 文件的：

x1 y1 z1 data1
x2 y2 z2 data2
x3 y3 z3 data3
…

此外，还对数据进行了清理，以去除股骨以外的所有东西，并按比例调整以适应股骨的实际尺寸。为了方便使用，还对数据进行了二值化处理。

此模型的 MPH 文件可以点击文末“阅读原文”查看，所以这里的说明只是对需要采取的步骤做一个大致的描述。在这篇文章的最后，我们将给出两个包含周围域的人类股骨的网格，在以后的计算中可以使用。我们还将讨论如何添加其他计算域和划分表面网格，以便能够在网格上应用适当的边界条件。

首先，将数据导入软件中的一个插值函数中，我们把这个函数命名为 human_femur。我们将通过两种可能的工作流程来获得股骨的网格：

1. 一个快速通道，在这个流程中，我们可以得到一个非常粗糙的数据网格表示，并将使用重新划分的网格将粗糙的表面变平。
2. 一个包含更多步骤的工作流程，使用亥姆霍兹平滑来确保我们可以得到一个平滑的形状。

我们从更快的工作流程开始。在这个例子中，我们将使用栅格数据集来表示数据，并使用过滤器数据集滤除数据来表示股骨。这将产生一个比较大块的表面网格，因此我们使用自由三角形操作来平滑不平整的表面。

在这个工作流程中，我们不需要几何节点，甚至不需要模型开发器中的组件节点，尽管我们稍后将在添加额外的计算域和创建股骨的计算网格时添加一个组件。在全局定义节点下设置插值函数，然后直接进入结果，定义一个栅格数据集。

添加插值函数后，直接设置栅格数据集。

COMSOL Multiphysics UI 显示了模型开发器，突出显示了网格数据集功能，相应的设置窗口与数据，参数边界和网格部分展开。

将栅格数据集作为一个规则的网格使用，并在这个网格上应用插值函数。这里，大约一半的图像分辨率被用作数据集的分辨率。现在，添加一个过滤器 数据集，输入表达式 human_femur(x,y,z)，并将下限 设置为 0.99。这将创建一个值为 0.99 的等值面，并包括所有超过该值的数据。换句话说，就是人类股骨的等值面。如果你想在最后得到一个域，过滤器数据集很适用。在第二个工作流程中，我们将改用分割数据集，其中你可以定义多个值来创建多个等值面，产生多个计算域。

由过滤器数据集生成的等值面。请注意，股骨在 Z 方向上被压缩，因为坐标轴经过了缩放。

过滤器数据集所代表的网格可以导入到网格划分序列中，在这里我们将重新划分面网格以平滑粗糙度。我们可以通过右击数据集 > 过滤器节点，选择创建网格零件来完成这个工作流程。

在第二个工作流程中，我们将使用在新组件中创建网格选项。在当前流程中，我们选择创建一个网格零件，因为我们将使用这个网格作为构造几何，这一点你很快就会看到。这样操作后，我们将在全局定义下自动设置一个网格零件，并添加一个导入 节点，导入过滤器数据集的网格。过滤器数据集的域单元通常质量很差，所以我们只导入三角形网格。然后，当表面网格达到要求时，我们将重新生成四面体网格。

通常情况下，我们需要一个特定的边界划分来施加边界条件。使用“与平面相交”操作进行平面切割，也可以对其他形状进行分割操作。您可以阅读文章“在 COMSOL Multiphysics® 中编辑和修复导入的网格”，看看相交和分割网格的比较。这里，我们将使用用球分割操作来得到一个大致的圆形边界。

用球分割操作（粉红色边）分割出一个大致的圆形边界（用蓝色突出显示），稍后我们可以在上面添加一个边界条件。为了便于观察，关闭了网格的渲染。

由于单元大小不一，并且有许多条状的三角形，使网格表面不光滑，不适合模拟。因此，需要对网格进行重新划分后再使用。使用更大的网格尺寸重新划分面网格，也可以使形状更加平滑。

使用自由三角形操作重新划分面网格。为了避免求解网格中的粗糙度，对网格进行均匀化操作，可以设置一个比网格中的不均匀度大的单元尺寸。这可以通过设置最大单元尺寸和最小单元尺寸为 0.005m 来实现。如下图所示，现在的表面已经平滑了很多，但仍然有点不均匀。在工作流程的后期，将再次重新划分边界网格，所以这些不规则的地方将被进一步均匀化。

图片显示的是重新划分面网格之后的股骨顶部。

如果你尝试对几乎彼此相交的表面重新进行网格划分，或者正在对那些几乎与自身相交的表面重新划分网格，请试试亥姆霍兹平滑功能。

现在，表面网格已经足够好了，可以继续在这个网格上添加额外的几何形状。

股骨域可能是你唯一感兴趣的域，但假设你还想在股骨周围再建立一个模拟域或一些螺钉来固定股骨中的裂缝，你会怎么做呢？这可以通过几何序列来实现。首先，在股骨模型上添加一个 3D 组件。然后，添加几何 > 导入节点，在其设置窗口中选择网格零件 1。这个节点将使我们能够添加和定位新的几何体，并且应该被标记为构造几何，因为我们将只使用股骨的面来定位所述的几何。一旦确定最终几何，构造几何对象就会被移除。你可以看到它是一个构造几何，因为它的边是虚线而不是实线。

将股骨作为构造几何对象导入几何序列，以定位与股骨有关的额外几何形状。

我们要创建的几何形状稍后将在网格划分序列中将与网格结合。作为一种替代方案，也可以将股骨的表面网格和创建的几何形状与几何序列中的几何并集操作结合起来。然而，通常在网格序列中使用网格并集操作更为稳健，这也是这里要遵循的工作流程。

你可以画出你在仿真过程中所需的任何几何形状，或者导入 CAD 设计。有了可用的构造几何，很容易绘制出所需的几何形状。在这个例子中，我们将添加两个螺钉来加固股骨，但无论你是创建还是导入一个或多个几何对象，原理都相似。螺钉的几何形状和股骨的构造几何形状如下面的左图所示。一旦创建完几何结构，并且已经构建了形成联合体，股骨的构造几何就会从生成的几何结构中移除，只剩下螺钉，如下图中间所示。对于导入的 CAD 设计，也可以添加移除细节或虚拟操作对几何形状进行简化，然后再与网格结合。关于如何结合几何和网格的更详细说明，请至官网查看 STL 导入教程系列的第二个教程。

在同一个组件中，转到网格 序列。确保添加两个导入节点，一个导入螺钉的几何，另一个导入节点从网格部分导入股骨的网格（如下图右）。使用网格并集操作将两个表面网格结合起来，并使用自由三角形节点对块的平面进行重新网格化，其中最大单元尺寸和最小单元尺寸设置为 0.005m。

左图：螺钉的几何形状与股骨的构造几何形状。中间：构建形成联合体节点后的螺钉的几何形状。右图：股骨和螺钉的网格导入到网格划分序列中。

最后，添加一个自由四面体操作，用四面体网格填充域。如果需要的话，还可以选择添加一个扫掠网格或边界层网格，以满足当前的物理场的需求。这个网格所定义的域和边界将在设置和求解模拟时使用。

现在网格已经创建好了，接下来可以进行仿真了，如下面的左图所示。但是，如果我们加入高阶离散化的物理场，会发生什么？比方说，我们想使用离散化设置二次巧凑边点单元创建一个固体力学仿真。弯曲单元的二阶节点将被放置在网格表面的弯曲表示上，类似于具有基础几何形状时的放置方式。为了查看弯曲单元在网格中的样子，我们可以通过将网格数据集的几何形函数设置为二次拉格朗日，将网格绘图的节点设置为几何形函数。二次拉格朗日绘图的设置并不与模型中的任何物理离散化相耦合，它只用于绘图。

左图：股骨和螺钉的体网格图。白色单元代表螺钉的域，股骨的单元被表达式 x> -0.015 过滤掉了。右图：在股骨的弯曲边界上显示的具有二阶节点的最终网格。绿色和黄色代表偏度。

为了比较生成的网格如何跟随数据的一种方式，我们可以添加一个表达式 human_femur(x,y,z) 的等值面图，以删格数据集的网格分辨率绘制。从下面的图中可以看出，网格很好地遵循了等值面。

生成的表面网格（灰色）和插值数据 human_femur(x,y,z) 的等值面（蓝色）的比较。

对于这个工作流程，简单来说，我们将添加一个应用了股骨数据的规则块状域，然后通过求解亥姆霍兹方程对数据进行平滑处理，得到一个光滑的股骨表面网格。我们仍然需要对股骨的表面网格进行重新网格化，但这次只是为了提高网格单元的质量。

为了求解亥姆霍兹方程，需要定义一个域，可以是长方体、球或其他一些简单适用于股骨数据的形状，如这篇文章中所述。要确保这个域尽可能的小，这样就可以用一个精细的网格来求解这个域。网格需要足够精细，以解析你想在表面网格中看到的数据的任何细节。

一个定义了导入数据的规则域上的足够细的六面体网格。

接下来，设置方程，添加一个瞬态研究并进行求解，以获得我们可以使用的解数据集。

为了避免数值伪影，并使之有可能使用较粗的网格，插值函数使用亥姆霍兹滤波的偏微分方程（PDE）进行低通滤波：

其中，R 是滤波器半径，近似代表低通滤波器的宽度。添加一个系数形式偏微分方程接口，如下图所示。

系数形式偏微分方程接口的设置。

这里，我们使用表达式 max((h[1/m])^2,R^2) 来表示扩散系数。软件会对网格中的每个节点使用最大的网格大小 h 和参数 R。在这个例子中，局部的网格大小 h 是恒定的，但是如果使用非结构化的四面体网格，h 会变化。请注意，由于所有操作都是在一个由网格表示的几何结构上进行的，每个单元网格都是线性的，所有的单元离散化阶次被都降低到线性单元，然后使用一个稳态求解器计算求解。

在六面体网格上绘制的插值函数 human_femur 的切面图（左）和平滑处理后的解数据（右）。

我们在平滑的解数据中使用一个介于 0（蓝色）和 1（红色）之间的值来过滤股骨的表面网格。在这个例子中，假设我们也希望股骨外的区域成为我们模拟的一部分，因此应该选择添加一个分割数据集。如果只对模拟股骨域感兴趣，可以使用过滤数据集，就像第一个工作流程中做的那样。水平值 为 0.6 的分割数据集就已经足够满足我们的目标了，如下图所示。

分割数据集的图，水平值设置为 0.6。等值面将模拟域划分为两个域，一个是股骨域，一个是股骨周围的域。

当对水平值满意时，右击分割或过滤数据集节点，选择在新组件中创建网格。如果需要创建额外的模拟域，请使用创建网格组件选项，并按照第一个工作流程的步骤进行。

这将自动设置一个带有导入节点的新组件，导入分割数据集的网格。清除导入域单元复选框，并将边界分割 设置为最小值。这将只导入表面网格，并将形成尽可能少的边界；7 个平面边界和 1 个股骨边界。我们将对这个表面网格进行改进，然后重新创建一个体网格。

使用自由三角形操作对同一网格划分序列中的表面进行重新划分。为了得到大小基本相同的三角形，将最大单元大小和最小单元大小都设置为 0.005m。

重新划分网格后的股骨和长方体表面。大小已经被修改，以生成一个大小基本相同的三角形的表面网格。为了便于观察，一些边界被隐藏起来。

现在，我们可以用一个四面体网格来填充这些域。添加的物理场接口将被应用到这个网格所定义的域和边界上。这样，我们就有了高质量的体网格，可以用于仿真对于其他的形状和应用，还可以在域中添加扫掠和边界层网格，但在这个例子中，四面体网格已经满足我们的需求了。

左图：股骨和周围长方体的体网格。白色的单元代表股骨域，周围的单元使用表达式 x > 0 过滤。右图：在股骨的弯曲边界上显示二阶节点的最终的网格。颜色代表偏度。为了便于观察，隐藏了一些边界起来。

在这个工作流程中生成的网格可以与导入的数据的等值面 human_femur(x,y,z) 进行比较，该等值面与分割数据集使用的等值面相同，被绘制在原始长方体域的网格上。为了生成一个更接近导入数据的网格，可以在求解亥姆霍兹方程的原始长方体域中使用更细的网格，相应的也可以降低扩散系数。

至此，这篇由扫描的数据创建表面网格的文章就讲完了。我们对表面网格进行了分割，使用了两种不同的技术对数据进行了平滑处理，并了解了如何将网格和几何图形结合起来。

作者：Hanna Gothäll
来源：COMSOL
原文：https://mp.weixin.qq.com/s/Tb0rSNE7uO_MknbDoX7HGg