Poisson方程设计Voronoi图力学结构

在结构力学、土木工程、建筑学等领域, 经常会涉及到这样一个问题: 如何获得某种几何结构或图案, 在满足一定强度等条件下尽可能的节省材料!

Topology Optimization

我最先想到的, 也是最自然的方法——拓扑优化（topology optimization）：根据给定的载荷、边界约束以及其他限制条件，在给定的区域内对材料分布进行优化的数学方法. 正好以前在读博士期间用 COMSOL Multiphysis 求解过很多问题, 其中包括一些拓扑优化的测试, 不得不佩服这个软体, 结果计算得真的很漂亮, 非常推荐大家学一下这个有限元软体.

A test about discontinuous Galerkin (dG) method for Poisson Equation?

cn.comsol.com

Convection dominated Convection-Diffusion Equation by upwind discontinuous Galerkin (dG) method?

cn.comsol.com

（自己之前在官网上上传的模型案例）

二维拓扑优化图片来源：http://cn.comsol.com/blogs/how-to-use-topology-optimization-results-as-model-geometries/ （侵删）

三维拓扑优化图片来源：http://cn.comsol.com/blogs/how-to-use-topology-optimization-results-as-model-geometries/ （侵删）

但不得不提的是, 拓扑优化在边界设置上很是繁琐, 而且计算也需要一定时间(PC机), 此外对于不同的载荷方向或者边界设置则需要再次重新计算, 使用起来并不是特别方便.

Voronoi Diagram

于是, 我调整了一下想法，如果我能提供个大概的稀疏孔洞的几何结构, 然后再验证在这种几何结构下是否满足我所需的结构强度不也可以么？因此该设计问题的关键就变成如何设计靠近最优解的稀疏孔洞模型. 事实上, 自然界给过很多启示: 肥皂泡的结构; 蜂巢的结构; 蜻蜓的翅膀; 细胞的结构等等, 貌似轻薄的自然物都具有这类结构.

自然界的某种启示图片来源网路（侵删）

肥皂泡的结构的成因为该状态下表面张力的内能最低,处于稳定的状态,细胞的结构和其类似,都可以看成空间的一种划分; 蜂巢的结构成因为该结构状态下密合度最高、所需材料最简、可使用空间最大,并且具有很高的结构强度; 蜻蜓的翅膀结构成因: 主脉里面是血管和神经，生长中越靠近脉的部分伸展越快，抻拉中自然就长成了这样.然而这些结构精致的美感，里面隐藏著巧妙的数学机制，即被称作Voronoi图的几何构造. 当前的太空梭、人造卫星、宇宙飞船在内部大量采用Voronoi图(蜂窝)结构，卫星的外壳也几乎全部是该结构; 在建筑中, 水立方则是模拟肥皂泡的Voronoi图结构等. 因此我们在做轻量化结构设计时也考虑使用Voronoi图结构的模型!

基于Voronoi图构造孔洞模板(Pattern)

Voronoi图的数学定义：一种基于距离的平面划分方法。在平面上有个若干不重合种子点(可以看成结构上需要去掉的孔洞位置)，把平面分为若干个区域，使得每个区域内的点到它所在区域的种子点的距离比到其它区域种子点的距离近。每个区域称为该种子点的Voronoi区域。Voronoi图是Delaunay三角剖分的对偶图。Voronoi图的每条边是由相邻种子点的垂直平分线构成，在边上的点到两个种子点的距离相等.

周期Voronoi图的生成方法：可以分为随机型和均匀型：

随机分布型: 为根据需要的孔洞密度, 在单位正方形中随机洒上若干个种子点(透气孔中心); 可以很方便地使用MATLAB中的random函数得到. 其最主要的特点是满足结构强度的前提下, 图样好看, 多变, 适合人们的审美; 但缺点是需要通过一定的有限元分析验证!

随机散点型

均匀分布型则: 根据需要的孔洞密度, 在单位正方形中均匀洒上若干个种子点(透气孔中心); 通常可以借助先在单位正方形上打上均匀的Delaunay三角网格, 以这些网格顶点作为种子点即可. 其主要特点是结构应力分布更均匀, 一般不存在应力集中的情况, 但图形比较单一.

均匀散点型

按道理来说, 在得到透气孔位置分布后, 我们即对这些种子点进行Voronoi图划分. 如MATLAB中包含内置的voronoi.m 函数, 即可直接输入种子点的坐标得到Voronoi图的结果, 包括各多边形的顶点等. 但是, 此时得到的Voronoi图 (蓝色多边形) 并不具有周期性对称结构: 即如果把此图看成一块」瓷砖」, 将它们铺满平面的话得到的图在边界处是不连续的!

为了得到周期的、可镶嵌整个平面的Voronoi图, 还需要做如下步骤的处理: 将单位正方形内的随机或均匀种子点上下左右对角等8个方位平移, 得到周期性的大正方形(均匀分布型的种子点需要增加边界点), 然后将大正方形内的所有种子点输入voronoi.m 函数, 得到如下结果:

我们只需要将上面的大正方形中心的单位正方形考虑既可.

不难发现, 虽然我们得到了所需的Voronoi图结构, 但在实际的使用中, 如此窄的线段无法满足结构的强度和贴合度的需求, 我们需要得到一定厚度的边界, 如果如下图直接将Voronoi图的每根线段加粗虽然能够得到一定厚度边界的孔洞结构图, 但其边缘不光滑, 容易给人一种很「机械」、「毛躁」的感觉！有没有，有没有！

如何加宽Voronoi图的每根线而且还能保证光滑？

Poisson方程

在数理方程里应该都接触过这个微分方程, 我们考虑带有齐次边界调节的微分方程

$egin{equation} left{ egin{array}{c} -left(frac{partial^2 u}{partial x^2}+frac{partial^2 u}{partial y^2}+frac{partial^2 u}{partial z^2} ight)=1,~~~ ext{in}~~Omega \ u=0,~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ext{on}~~partialOmega \ end{array} ight. end{equation}$

该方程可描述大量的物理模型，热传导, 浓度扩散，电磁学，引力场等等. 应该还记得该函数的解具有很好的光滑性质，因此我们利用此Poisson方程解的某一条等高线(contour)作为孔洞的光滑边缘，并可方便利用等高线的取值来控制边的宽细！

因此, 使用上述方法使用在单位正方形附近Voronoi多边形构成的每一个多边形求解Poisson方程(利用有限元方法), 得到如下的结果: (分别随机型和均匀型)

Voronoi图各个区域Poisson方程的解(单位化后的解，即最大值为1)

从上图中很容易看出, 调整不同的等高线值, 0~1之间, 既可得到不同宽度的Voronoi图结构.

光滑孔洞Pattern

根据需要, 选取一定粗细的边缘, 然后对其进行黑白填充, 并截取单位正方形部分, 即得到如下结果:

总结过程

STEP1: 利用现有的三角网格剖分软体得到[0,1]x[0,1]区间上的网格划分(可以控制大小), 或者先将[0,1]x[0,1]的边界划分若干等分, 再利用区间的均匀分布随机生成内点;
STEP2: 先将单位区间的边界上下左右延拓为」井」字, 再将三角网格的划分内点(或者随机生成的内点), 在上下左右斜向等8个方位( x坐标+-1或不变, y坐标+-1或不变)平移, 得到[-1,2]x[-1,2]区间上的离散点分布, 这些点关于边界具有周期性;
STEP3: 根据[-1,2]x[-1,2]区间上的离散点得到其相应的voronoi图!
STEP4: 根据[-1,2]x[-1,2]区间上的voronoi图, 得到包围每个离散点的多边形, 筛选出单位区间[-1,2]x[-1,2]内以及附近的那些voronoi多边形;
STEP5: 对每一个筛出的多边形求解Poisson方程: