基于STL几何的粒子模拟

电热输运

基于STL几何的粒子模拟

确定性方法需要网格文件输入，比如下面这张图是鲍华老师GiftBTE求解器的Workflow，前处理通过GMSH或者COMSOL生成网格，输入到GiftBTE中进行求解，然后再用开源的Paraview做后处理可视化。对于粒子的蒙特卡洛模拟，求解本身并不需要网格，而关注的是粒子和边界之间的相互作用，处理粒子-边界碰撞后的透射或反射过程，此时需要的是表面的信息，比如法向量、组成面的各个点的位置、以及如何计算线-面交点和判断粒子是否在结构内部的问题。这时候如果我们想要使用前处理CAD软件进行建模，就需要借助STL文件格式。

img

STL文件

正十二面体

STL文件可以以两种格式存储，二进制格式和ASCII（文本）格式。二进制格式更为紧凑，读写速度更快。但无论哪种格式，STL文件都通过一些列三角面片（facet）描述了一个3D对象的表面几何形状。每个三角面片由12个浮点数来表示，包括三个顶点 (vertex)和一个法线向量 (normal)，法线向量表示该面片的朝向。法线方向指向结构的外侧，顶点的顺序遵循右手法则。四个手指按顺序握住三个顶点，拇指指向的方向就是法线方向。比如下面是一个在COMSOL中搭建的1 x 2 x 3的长方体结构导出的STL文件的一部分，一个长方体包含6个面，因此被拆分成了12个小三角片面，每个长方形面有两个小三角片面组成。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
solid COMSOL rendering object fin
facet normal -1 0 0
outer loop
vertex 0 0 3
vertex 0 2 3
vertex 0 2 0
endloop
endfacet
facet normal -1 0 0
outer loop
vertex 0 2 0
vertex 0 0 0
vertex 0 0 3
endloop
endfacet
......
endsolid COMSOL rendering object fin

ASCII STL有以下几个关键字：

  • solid：表示一个STL文件的开始。在某些情况下，后面可能跟随一个名称，但这个名称通常被视为可选的。
  • facet normal：后面跟着三个数值，表示接下来定义的面片的法线向量。法线向量是一个从面片指向外部的向量，用于指示面片的朝向。
  • outer loop：标志着一个面片的顶点定义开始。每个面片由三个顶点定义，这三个顶点按照逆时针顺序排列（从法线向量的方向看去）。
  • vertex：后面跟着三个数值，表示一个顶点的x、y、z坐标。在一个facet定义中，会有三个这样的vertex行，分别对应面片的三个顶点。
  • endloop：表示一个面片的顶点定义结束。
  • endfacet：标志着一个面片的定义结束。
  • endsolid：表示STL文件的结束。在某些情况下，后面可能跟随一个名称，但这通常被视为可选的。

在STL文件标准中，理论上可以包含多个solid和对应的endsolid声明，每一对solid...endsolid表示一个单独的三维模型或对象。但是大多数的3D打印软件和CAD软件都假设一个STL文件代表一个单一的连续的几何体，并可能只识别文件中的第一个solid声明，因此一般一个STL就代表一个几何体了。当然，这和蒙特卡罗模拟的需求也是一样的，因为一个结构体是一种材料，我们每次只跟踪粒子在同一个结构内部的运动。此外，对于solid和对应的endsolid后面的名称，理论上应该匹配以表示它们之间包含的是同一个对象的定义，但是在许多处理STL文件的软件中，也并不强制要求这些名称必须匹配，甚至根本不使用这些名称。

Python的STL库

Python有一些库用来解析和处理STL格式的几何。比如如果我们只想读取STL几何并且做一些基本的计算，一个轻量级的库是numpy-stl，可以通过pip安装：

1
pip install numpy-stl

用numpy-stl加载STL文件很容易

1
2
3
4
5
6
7
8
9
10
import numpy as np
from stl import mesh

# 加载STL文件
your_mesh = mesh.Mesh.from_file('path_to_your_file.stl')

# 查看模型顶点
print(your_mesh.vectors)
# 查看模型法向量
print(your_mesh.normals)

提取的结果都是NumPy数组。如果希望进行一些更复杂的几何运算，比如判断一个坐标是否在STL结构内部，或者判断一条射线是否与STL几何相交并求出交点，或者判断两个STL几何是否碰撞，可以用更复杂一些的trimesh库

1
pip install trimesh

调用方法和numpy-stl也基本上是一致的

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
import numpy as np
import trimesh

# 加载模型
mesh = trimesh.load_mesh('path_to_your_model.stl')

# 定义一组要测试的点
points = np.array([
    [0, 0, 0],    # 示例点1
    [10, 10, 10]  # 示例点2
])

# 使用 contains 方法判断点是否在网格内部
inside = mesh.contains(points)

# 打印结果
for i, point in enumerate(points):
    print(f"Point {i} {point} is inside the mesh: {inside[i]}")

# 定义射线的起点和方向
ray_origins = np.array([[0, 0, 0]]) # 射线起点
ray_directions = np.array([[0, 0, 1]]) # 射线方向

# 检测射线与网格的交点
locations, index_ray, index_tri = mesh.ray.intersects_location(
    ray_origins=ray_origins,
    ray_directions=ray_directions
)

# 打印交点位置
print(locations)

蒙特卡罗

在大部分情况下，声子蒙特卡罗模拟的对象都是一系列长方体，比如纳米薄膜。即使对于实际晶体管HEMT或者FinFET（SISPAD. IEEE, 2014: 17-20.）结构，也就是一系列的方块堆起来。但是有CAD软件接口和没有CAD软件接口，听起来软件的完善程度是完全不一样的。为了实现CAD建模+蒙特卡罗模拟流程，需要将建立的几何结构导出为STL格式文件。

image-20240304155405202

东京大学Shiomi Lab的邵成博士（现山东高等技术研究院副研究员）前两年用Fortran编写了P-TRANS，用于模拟纳米结构材料的等效热导率，程序也已经开源在了Github上，可以作为好的学习和参考的资料。这个程序支持从CAD导出STL文件驱动蒙特卡罗模拟，这里简要介绍一下P-TRANS模拟的逻辑。

RTMC_illustration-17103122764571

首先处理STL文件中的一系列三角面片，将它们封装成一个个surface结构体，一系列surface结构体构成了grain结构体，作为一块结构。假如现在声子处于该grain中，所处的位置为O，群速度方向矢量为\(\vec{k}\)，现在想确定声子下一步将要运动到的位置，我们需要判断声子是否与结构体的边界发生了碰撞。在P-TRANS中，是通过比较距离来实现碰撞检测的。首先，抽样声子在没有障碍时将要运动的距离\(l_p\)， \[ l_p = - \Lambda \cdot \ln(R) \] 其中\(R\)是处于\((0, 1)\)间均匀的随机数，具体抽样原理可以参考之前的博客 蒙特卡洛中的频率抽样。假如O和\(\vec{k}\)表示的射线和结构体的交点为H，现在只需要判断OH的距离和\(l_p\)的关系就可以判断粒子是否发生碰撞了。如果\(\text{OH}>l_p\)则说明没有发生碰撞，声子下一步所处的位置还在结构内部，继续抽样运动即可；如果\(\text{OH}<l_p\)​则说明发生了碰撞，此时需要将声子的位置更新到交点处，按照边界条件对声子进行相应的处理。

image-20240313205128010

关于交点H的确定，P-TRANS所采用算法的大致思路是找到O和\(\vec{k}\)所组成的射线与结构体所有多边形面（对于STL格式几何就都是三角面）的交点，其中距离O最近的点，就是实际相交的点。那么如何找到与某一个多边形面，比如与V1-V4构成的长方形的交点呢？这里所采用的思路是，先找到这条射线与这个多边形面所在的无限大平面的交点，然后再判断交点是否处于V1-V4围成的多边形内部，具体过程如下：

  • 确定交点H

    要找到与无限大平面的交点H，由于已经知道了\(\vec{k}\)，因此只需要知道OH的长度就可以了。 而\(\vec{k}\)和法向量\(\vec{n}\)的夹角我们是知道的，因此只要求出O点到这个平面的垂直距离，除以\(\cos\theta\)，就是OH的长度。要求垂直距离，只需要随便连接O和平面上一点，然后计算这个向量在平面法向量上的投影就可以了。因此OH长度的表达式可以整理为 \[ O H=\frac{O H_{\perp}}{\cos \theta}=\frac{\overrightarrow{O V}_1 \cdot \vec{n}}{\vec{k} \cdot \vec{n}} \]

  • 确定交点H是否在多边形区域内

    P-TRANS采用这样的逻辑判断H是否在区域内部，做所有\(\overrightarrow{V_nV_{n+1}}\)和\(\overrightarrow{V_nH}\)的叉积，比如对于V1-V2-V3-V4组成的长方体，这里就有4个叉积。如果H在多边形内部，那么此时叉积向量都指向同一个方向。假如此时H在多边形外部，比如到了红色的位置，那么\(\overrightarrow{V_4V_{1}}\)和\(\overrightarrow{V_4H}\)的叉积向量和其他几个叉积向量的方向将相反。如果H在内部，则保留OH的长度。

  • 遍历所有的外表面，执行上述过程，找到最短的OH长度。最短的OH长度所对应的H点，则是实际交点。

多态设计

P-TRANS为为科研工作者进一步理解微纳米尺度声子传输提供很好的工具，大大促进了传热社区的繁荣与发展，但其实也有些设计的地方是可以修改和提高的。P-TRANS是用Fortran编写的，对表面对象和区域对象做了一些基本的封装，但基本上还是面向过程的。比如为了实现对不同类型几何体的支持，采用了一系列if-else语句来执行类型检查以进行不同的操作，当需要添加新的行为或条件时，这需要直接修改的is-else链，这是不满足开闭原则的（OCP）。另一方面，随着需求的增多，if-else 语句可能会变得越来越复杂，耦合度会越来越上升，维护的难度也会逐渐增加。为什么会这样呢，主要是我们做科学计算的逻辑大部分情况下还是面向过程的，直接从实现的角度去考虑。实际上当遵循依赖倒置原则（DIP），面向接口进行编程，让实现的细节依赖于抽象，整个耦合度就会降下来了。

另一方面，为了实现对各种复杂结构的支持，P-TRANS采用了相当通用的算法去做碰撞检测。但实际上，可能90%的情况下，我们都在和最简单的长方体打交道。判断一个点是否在一个长方体内部，可以用更加简单快速、容易维护的方法去处理。如何实现通用性支持的同时，又能针对具体的场景采用最合适的方式呢？同样是面向接口编程，通过多态实现这种灵活性。比如定义一个抽象基类把接口抽象出来，子类用策略模式，为不同的结构开发具体的算法。这样整个程序也更加灵活，易于维护和扩展。