三周期极小曲面（TPMS）是一类具有零平均曲率的周期性曲面，在自然界中广泛存在，如蝴蝶翅膀的鳞片、骨骼的微结构等。基于TPMS设计的晶格结构因其优异的力学性能、高比表面积和可控的孔隙率，成为增材制造（3D打印）领域的研究热点。本文提供一个完全自包含的Python脚本，能够根据用户指定的TPMS类型、结构形式（杆状网络型或片状壳体型）和目标体积分数，自动生成封闭的STL网格文件，可直接用于3D打印或有限元分析。
代码的核心思想是将TPMS的隐函数方程离散化到规则网格上，通过Marching Cubes算法提取等值面，并对网格边界添加盖面以保证模型的封闭性。为了精确控制体积分数，程序采用先粗调后二分法的方式自动校准水平集常数r。整个过程无需任何外部几何建模软件，仅依赖NumPy、SciPy和scikit-image三个科学计算库，即可在数秒内生成高精度网格。
import numpy as np
from skimage.measure import marching_cubes, find_contours
from scipy.spatial import Delaunay
# ==================== TPMS 隐函数定义 ====================
deftpms_function(x, y, z, tpms_type):
"""
计算指定类型TPMS的隐函数值 F(x,y,z)
支持类型: 'P', 'G', 'D', 'I', 'S', 'F', 'N'
"""
pi2 = 2 * np.pi
pi4 = 4 * np.pi
if tpms_type == 'P':
return - (np.cos(pi2 * x) + np.cos(pi2 * y) + np.cos(pi2 * z))
elif tpms_type == 'G':
return (np.cos(pi2 * x) * np.sin(pi2 * y) +
np.cos(pi2 * y) * np.sin(pi2 * z) +
np.cos(pi2 * z) * np.sin(pi2 * x))
elif tpms_type == 'D':
return (np.sin(pi2 * x) * np.sin(pi2 * y) * np.sin(pi2 * z) +
np.sin(pi2 * x) * np.cos(pi2 * y) * np.cos(pi2 * z) +
np.cos(pi2 * x) * np.sin(pi2 * y) * np.cos(pi2 * z) +
np.cos(pi2 * x) * np.cos(pi2 * y) * np.sin(pi2 * z))
elif tpms_type == 'I':
return (2 * (np.cos(pi2 * x) * np.cos(pi2 * y) +
np.cos(pi2 * y) * np.cos(pi2 * z) +
np.cos(pi2 * z) * np.cos(pi2 * x)) -
(np.cos(pi4 * x) + np.cos(pi4 * y) + np.cos(pi4 * z)))
elif tpms_type == 'S':
return (np.cos(pi4 * x) * np.sin(pi2 * y) * np.cos(pi2 * z) +
np.cos(pi2 * x) * np.cos(pi4 * y) * np.sin(pi2 * z) +
np.sin(pi2 * x) * np.cos(pi2 * y) * np.cos(pi4 * z))
elif tpms_type == 'F':
return - (4 * (np.cos(pi2 * x) * np.cos(pi2 * y) * np.cos(pi2 * z)) -
(np.cos(pi4 * x) * np.cos(pi4 * y) +
np.cos(pi4 * y) * np.cos(pi4 * z) +
np.cos(pi4 * z) * np.cos(pi4 * x)))
elif tpms_type == 'N':
return3 * (np.cos(pi2 * x) + np.cos(pi2 * y) + np.cos(pi2 * z)) + \
4 * np.cos(pi2 * x) * np.cos(pi2 * y) * np.cos(pi2 * z)
else:
raise ValueError(f"不支持的TPMS类型: {tpms_type}")
# ==================== 网格生成辅助函数 ====================
defgenerate_grid(x_range, y_range, z_range, n):
"""
在指定长方体区域内生成规则网格
返回: 网格点数组, 各轴坐标序列, 网格间距
"""
x = np.linspace(x_range[0], x_range[1], n + 1)
y = np.linspace(y_range[0], y_range[1], n + 1)
z = np.linspace(z_range[0], z_range[1], n + 1)
xx, yy, zz = np.meshgrid(x, y, z, indexing='ij')
spacing = ((x_range[1] - x_range[0]) / n,
(y_range[1] - y_range[0]) / n,
(z_range[1] - z_range[0]) / n)
return xx, yy, zz, x, y, z, spacing
# ==================== 边界盖面生成 ====================
defadd_isocaps(x_coords, y_coords, z_coords, phi, level=0):
"""
为三维标量场phi的边界添加等值面盖面（三角剖分）
返回盖面的顶点和三角形索引（局部坐标）
"""
nx, ny, nz = phi.shape
all_verts = []
all_faces = []
offset = 0
# 六个边界平面的定义
planes = [
(0, 0, x_coords[0], lambda u, v: (x_coords[0], u, v)),
(0, -1, x_coords[-1], lambda u, v: (x_coords[-1], u, v)),
(1, 0, y_coords[0], lambda u, v: (u, y_coords[0], v)),
(1, -1, y_coords[-1], lambda u, v: (u, y_coords[-1], v)),
(2, 0, z_coords[0], lambda u, v: (u, v, z_coords[0])),
(2, -1, z_coords[-1], lambda u, v: (u, v, z_coords[-1]))
]
for axis, idx, fixed_val, to3d in planes:
# 提取二维切片
if axis == 0:
slice_phi = phi[idx, :, :] # shape (ny, nz)
u_coords = y_coords
v_coords = z_coords
elif axis == 1:
slice_phi = phi[:, idx, :] # shape (nx, nz)
u_coords = x_coords
v_coords = z_coords
else:
slice_phi = phi[:, :, idx] # shape (nx, ny)
u_coords = x_coords
v_coords = y_coords
# 提取等值线轮廓
contours = find_contours(slice_phi, level)
for contour in contours:
# 将轮廓点从像素坐标映射到实际坐标
u = u_coords[0] + contour[:, 1] * (u_coords[-1] - u_coords[0]) / (len(u_coords)-1)
v = v_coords[0] + contour[:, 0] * (v_coords[-1] - v_coords[0]) / (len(v_coords)-1)
points_2d = np.column_stack((u, v))
if len(points_2d) >= 3:
tri = Delaunay(points_2d) # 三角剖分
verts_3d = np.array([to3d(ui, vi) for ui, vi in points_2d])
faces_3d = tri.simplices
all_verts.append(verts_3d)
all_faces.append(faces_3d + offset)
offset += len(verts_3d)
if all_verts:
return np.vstack(all_verts), np.vstack(all_faces)
return np.empty((0, 3)), np.empty((0, 3), dtype=int)
# ==================== 体积与表面积计算 ====================
defstl_volume(vertices, faces):
"""
计算封闭三角网格的体积（有向）和表面积
"""
v0 = vertices[faces[:, 0]]
v1 = vertices[faces[:, 1]]
v2 = vertices[faces[:, 2]]
# 体积：散度定理
vol = np.sum((v0[:, 0] * (v1[:, 1] * v2[:, 2] - v2[:, 1] * v1[:, 2]) -
v1[:, 0] * (v0[:, 1] * v2[:, 2] - v2[:, 1] * v0[:, 2]) +
v2[:, 0] * (v0[:, 1] * v1[:, 2] - v1[:, 1] * v0[:, 2])) / 6.0)
# 表面积
normals = np.cross(v1 - v0, v2 - v0, axis=1)
area = 0.5 * np.linalg.norm(normals, axis=1)
return vol, np.sum(area)
# ==================== 二进制STL写入 ====================
defwrite_binary_stl(filename, vertices, faces, title="TPMS Lattice"):
"""
将三角网格写入二进制STL文件
"""
vertices = vertices.astype(np.float32)
v0 = vertices[faces[:, 0]]
v1 = vertices[faces[:, 1]]
v2 = vertices[faces[:, 2]]
normals = np.cross(v1 - v0, v2 - v0, axis=1)
norms = np.linalg.norm(normals, axis=1, keepdims=True)
norms[norms == 0] = 1
normals /= norms
facets = np.hstack([normals, v0, v1, v2]).astype(np.float32)
with open(filename, 'wb') as f:
f.write(title.encode('ascii')[:80].ljust(80, b'\x00'))
f.write(np.array(len(faces), dtype=np.uint32).tobytes())
for facet in facets:
f.write(facet.tobytes())
f.write(b'\x00\x00')
# ==================== 体积分数校准 ====================
defcalibrate_r(xx, yy, zz, x_coords, y_coords, z_coords, spacing, tpms_type, struct_type,
vol_frac_target, r_start=-0.07, r_step=0.1, tol=1e-4, max_iter=50):
"""
通过粗调+二分法自动寻找合适的水平集常数r，使体积分数接近目标值
struct_type: 'net' (杆状) 或 'sheet' (片状)
"""
F_val = tpms_function(xx, yy, zz, tpms_type)
# 粗调
r = r_start
vol_frac = 0
for _ in range(max_iter):
if struct_type == 'net':
phi = F_val + r
else:
phi = -(F_val + r) * (F_val - r)
verts, faces, _, _ = marching_cubes(phi, level=0, spacing=spacing)
verts[:, 0] += xx[0,0,0]
verts[:, 1] += yy[0,0,0]
verts[:, 2] += zz[0,0,0]
if len(verts) > 0and len(faces) > 0:
cap_verts, cap_faces = add_isocaps(x_coords, y_coords, z_coords, phi, level=0)
if len(cap_verts) > 0:
all_verts = np.vstack([verts, cap_verts])
all_faces = np.vstack([faces, cap_faces + len(verts)])
else:
all_verts, all_faces = verts, faces
vol, _ = stl_volume(all_verts, all_faces)
vol_frac = abs(vol)
else:
vol_frac = 0
if abs(vol_frac - vol_frac_target) < tol:
break
if vol_frac < vol_frac_target:
r += r_step
else:
r -= r_step
if r > 1or r < -1:
break
# 二分法精调
if vol_frac > vol_frac_target:
r1, r2 = r - r_step, r
else:
r1, r2 = r, r + r_step
for _ in range(max_iter):
r = (r1 + r2) / 2
if struct_type == 'net':
phi = F_val + r
else:
phi = -(F_val + r) * (F_val - r)
verts, faces, _, _ = marching_cubes(phi, level=0, spacing=spacing)
verts[:, 0] += xx[0,0,0]
verts[:, 1] += yy[0,0,0]
verts[:, 2] += zz[0,0,0]
if len(verts) > 0and len(faces) > 0:
cap_verts, cap_faces = add_isocaps(x_coords, y_coords, z_coords, phi, level=0)
if len(cap_verts) > 0:
all_verts = np.vstack([verts, cap_verts])
all_faces = np.vstack([faces, cap_faces + len(verts)])
else:
all_verts, all_faces = verts, faces
vol, _ = stl_volume(all_verts, all_faces)
vol_frac = abs(vol)
else:
vol_frac = 0
if abs(vol_frac - vol_frac_target) < tol:
break
if vol_frac > vol_frac_target:
r2 = r
else:
r1 = r
return r, vol_frac
# ==================== 主程序 ====================
defmain():
# ---------- 用户可调参数 ----------
n_final = 80# 最终网格分辨率（越高越精细）
n_calib = 60# 校准用分辨率
tpms = 'P'# TPMS类型: P, G, D, I, S, F, N
struct_type = "net"# "net" 或 "sheet"
vol_frac_target = 0.27# 目标体积分数 (0~1)
num_cell = 1# 每个方向上的晶胞个数
target_size = 1.0# 最终模型的边长
center = np.array([0, 0, 0]) # 模型中心坐标
# -------------------------------
print("=" * 50)
print("TPMS 晶格结构生成器")
print(f"类型: {tpms}, 结构: {struct_type}, 目标体积分数: {vol_frac_target}")
print("=" * 50)
# 1. 校准水平集常数 r
print("正在校准水平集常数 r ...")
x_range = (0, 1)
y_range = (0, 1)
z_range = (0, 1)
xx, yy, zz, xc, yc, zc, sp = generate_grid(x_range, y_range, z_range, n_calib)
r, vol_frac = calibrate_r(xx, yy, zz, xc, yc, zc, sp, tpms, struct_type, vol_frac_target)
print(f"校准完成: r = {r:.6f}, 实际体积分数 = {vol_frac:.6f}")
# 2. 生成最终网格
print("正在生成最终网格 ...")
x_range_f = (0, num_cell)
y_range_f = (0, num_cell)
z_range_f = (0, num_cell)
xx_f, yy_f, zz_f, xc_f, yc_f, zc_f, sp_f = generate_grid(x_range_f, y_range_f, z_range_f, n_final)
F_final = tpms_function(xx_f, yy_f, zz_f, tpms)
if struct_type == 'net':
phi_f = F_final + r
else:
phi_f = -(F_final + r) * (F_final - r)
verts, faces, _, _ = marching_cubes(phi_f, level=0, spacing=sp_f)
verts[:, 0] += x_range_f[0]
verts[:, 1] += y_range_f[0]
verts[:, 2] += z_range_f[0]
cap_verts, cap_faces = add_isocaps(xc_f, yc_f, zc_f, phi_f, level=0)
if len(cap_verts) > 0:
all_verts = np.vstack([verts, cap_verts])
all_faces = np.vstack([faces, cap_faces + len(verts)])
else:
all_verts, all_faces = verts, faces
# 平移和缩放
all_verts -= center
all_verts *= target_size
vol, _ = stl_volume(all_verts, all_faces)
vol_frac_final = abs(vol) / (target_size ** 3)
print(f"最终模型体积分数: {vol_frac_final:.6f}")
write_binary_stl('lattice.stl', all_verts, all_faces)
print("STL文件已保存为 lattice.stl")
if __name__ == "__main__":
main()
代码使用说明
上述代码完全独立，只需安装依赖库即可运行：
pip install numpy scipy scikit-image
所有核心算法均内嵌于脚本中，无需任何外部几何处理软件。用户只需修改主函数开头的参数（如TPMS类型、结构形式、目标体积分数、晶胞数量和模型尺寸），运行后即可在当前目录下生成lattice.stl文件。该文件可直接导入3D建模软件（如Blender、MeshLab）或切片软件，用于增材制造。
技术要点深度解读
TPMS的数学本质是隐函数F(x,y,z)=c，其中c为水平集常数。当c=0时，曲面通过原点；改变c可连续调整曲面的体积分数。对于杆状网络结构（net），曲面定义为F(x,y,z)=r；对于片状壳结构（sheet），曲面定义为(F+r)(F-r)=0，相当于两个等值面的夹层区域。通过二分法自动寻找合适的r，使得模型体积分数逼近目标值，这是整个算法的关键所在。
在网格生成环节，Marching Cubes算法将隐函数离散化后提取三角形面片，但原始输出是开放的（模型边界处有孔洞）。为了获得可用于实体打印的封闭模型，代码额外为六个边界平面添加了盖面：对每个边界平面上的隐函数值进行二维等值线提取，再使用Delaunay三角剖分生成填充面。这种方法保证了最终模型的水密性，可直接用于体积计算和3D打印。
体积分数的精确控制依赖于封闭网格的体积计算。代码通过散度定理计算有向体积，避免了繁琐的四面体剖分，仅需一次遍历所有三角形即可获得高精度结果。校准过程中，粗调阶段采用变步长搜索快速逼近目标区间，精调阶段使用二分法进一步收敛，兼顾了效率与准确性。
该实现将数学定义、离散化、网格生成、边界处理、体积校准和文件输出整合为一条流畅的流水线，为TPMS晶格结构的快速设计提供了一套开箱即用的解决方案。无论是用于科研中的力学性能分析，还是工程中的轻量化设计，都能从中获得高效支持。
10个月宝宝每天需要喝多少奶粉？