CN109741452B - 一种地质体3d打印自支撑结构自动生成方法 - Google Patents

Abstract

本发明公开了一种地质体3D打印自支撑结构自动生成方法。该方法针对没有自支撑性的地质体3D打印模型无法表现地质体原始产状(倾向和倾角)的问题，自动生成保证平衡性、结构强度和最小可见性的地质体3D打印模型的自支撑结构。为实现上述目标，构建了顾及自支撑结构平衡性、结构强度和最小可见性的最优化模型。针对最优化模型的求解，利用搜索算法得到带平衡性约束的极小化可见性场，并通过优化支撑结构的半径，增强支撑结构的强度。本发明所公开的地质体3D打印自支撑结构自动生成方法，对于提升地质3D打印模型的实用性具有很好的应用前景。

Description

一种地质体3D打印自支撑结构自动生成方法

技术领域

本发明涉及三维地质建模和3D打印领域，特别涉及一种地质体3D打印自支撑结构自动生成方法。

背景技术

基于FDM的3D打印机因为材料的廉价而被各大科研院所、生产厂商使用。这种打印机对于悬空结构、非平衡的地质体模型需要添加额外的支撑结构来保证打印的顺利进行。

在地质领域中，通过三维建模软件得到的地质体模型，在计算机虚拟空间下进行观察时，往往需要旋转模型或缩放模型，以实现它的全方位的观察。从某种程度上来说，对地质体模型的观察造成了一定的不便。此外，地质体模型由于不规则、空间结构复杂而难以通过传统工艺实现模型的实物展示。3D打印通过层叠打印材料或者通过激光烧结等方式，对于内部结构复杂的模型具有良好的复原效果。同时，3D打印得到的地质体模型，可以直观地表示它在地表或地下空间的形状，更有利于地质人员的观察。但是，地质体相对于3D打印中常见物体具有如下的特殊性：

1)地质体与机械零件等物体不同，主要用于向用户展示，供用户分析，在展示时要求还原姿态(产状)；

2)地质体与面向展示的文化创意品、展览品等人造物体不同，其本身没有自支撑性。

发明内容

本发明目的是，针对没有自支撑性的地质体3D打印模型无法表现地质体原始产状(倾向和倾角)，提供一种自动生成保证平衡性、结构强度和最小可见性的地质体3D打印模型的自支撑结构的自动生成方法，从而在保证自支撑结构平衡性、结构强度和最小可见性的前提下，实现地质体3D打印模型对地质体原始产状的展现。

为了实现上述技术目的，本发明的技术方案是：

一种地质体3D打印自支撑结构自动生成方法，包括以下步骤：

步骤一，建立自支撑结构平衡性的度量：

给定地质体三维模型的任意自支撑结构模型，以地质体模型的质心投影落在模型与支撑平面接触点所组成的凸包的内部为目的，先计算地质体三维模型的质心位置，然后以质心投影与凸包边界的最近距离作为此自支撑结构平衡性的度量；

步骤二，建立自支撑结构强度的度量：

给定地质体三维模型的任意自支撑结构模型，通过有限元方法计算模型在地质体3D打印模型和自身重力作用下的von Mise应力，作为此自支撑结构在重力下结构强度的度量；

步骤三，计算反映自支撑结构隐藏程度的可见性场：

基于光线投射算法计算目标点到观察点之间的可见性，通过计算地质体模型所在的三维栅格空间在观察范围内可见性的期望值，表达空间中任意点在地质体遮挡下的可见性，得到反映自支撑结构隐藏程度的可见性场；

步骤四，建立自支撑结构最优化数学模型：

基于步骤一至步骤三分别得到的自支撑结构平衡性、结构强度和可见性度量，综合构建自支撑结构的最优化数学模型；

步骤五，建立自支撑结构的最小可见性候选集：

基于步骤四所建立的自支撑机构最优化模型，采用动态规划方法，通过在三维可见性场中逐层递推的策略，计算得到用于放置地质体3D打印模型的支撑平面任意一点处，满足自支撑结构最优化模型中最优可见性的自支撑结构，对支撑平面上所有点按上述方法计算最小可见性支撑柱，综合建立自支撑结构的最小可见性候选集；

步骤六，计算自支撑结构位置：

基于步骤五所得到的自支撑结构最小可见性候选集，根据自支撑结构最优化模型中的平衡性约束，采用分支限界策略，在四叉树索引加速下，从最小可见性候选集中搜索获得最优支撑结构所在位置；

步骤七，计算自支撑结构半径：

根据步骤六得到的最优支撑结构所在位置，搜索自支撑结构的半径，以满足自支撑结构最优化模型中的结构强度为约束，根据搜索得到的支撑结构半径，最终建立地质体3D打印自支撑结构三维模型。

所述的方法，所述的步骤一中，首先对地质体三维模型进行四面体网格剖分，然后采用以下方式计算地质体三维模型的质心位置：

其中，

为四面体的质量，

为对应的四面体质心，T_i为对当前的地质体三维模型进行四面体剖分后得到的四面体。

所述的方法，所述的步骤一中，计算质心投影与凸包边界的最近距离包括以下步骤：

计算质心G到所支撑平面上的投影G_proj，判断G_proj是否落在模型与支撑平面的接触点所组成的凸包的内部，即G_proj∈H，其中H为凸包，则有：

(P_i-G_proj)(P_(i+1)％k-G_proj)＜0，i＝0，l，…n-1

其中{P_i}为凸包H的顶点，即模型与支撑平面接触点，k为顶点的数量；寻找凸包H的边界上距离G_proj的最近点P*，计算G_proj到P*的距离d_G；如果G_proj落在凸包H的内部，则质心投影与凸包边界的最近距离为d_G，否则，令质心投影与凸包边界的最近距离为-d_G。

所述的方法，所述的步骤二中，利用有限元方法来计算在地质体3D打印模型及其支撑结构在自身重力作用下的von Mise应力包括以下步骤：

1)对地质体三维模型及其支撑结构模型进行四面体网格剖分；

2)通过求解刚度方程，得到自支撑结构在自身重力作用下的位移：

Kd＝f

其中K表示有限元刚度矩阵，d表示位移向量，f为表示模型自身重力的向量。

其中有限元刚度矩阵采用的是有限元通用的刚度矩阵。

3)基于支撑结构的位移，计算支撑结构的内应力σ：

σ＝Cε

其中，C为表达3D打印材料力学性质的弹性矩阵，

表示梯度算子。其中弹性矩阵是确定力学参数后，通过有限元法自动生成的。

4)对模型剖分得到的每个四面体单元，根据其所对应的内应力σ^e，计算单元的vonMise应力：

所述的方法，所述的步骤三包括以下步骤：

对于地质体所在三维空间M_V内任一目标点x和任一观察的位置即视点y，有x和y的可见性函数V(x，y)如下：

其中值为0时表示不可见，1表示可见；

以观察点集合作为上半球空间，则半球空间对目标点x的可见性V_x的期望通过下式表达：

E(V_x)＝∫_Ωp(y)V(x，y)dy， (1)

其中p(y)表示在半球上的分布；

将公式(1)离散化为每个观察点对当前目标点x的求和形式：

E(V_x)＝∑V(x，y)

通过计算空间中每点的可见性期望值，得到地质体所在空间M_V的可见性场，则M_V中支撑结构S的可见性V(S)表示为：

所述的方法，所述的步骤四中，最优化模型为：

其中，V(S)为根据步骤三计算得到的自支撑结构的可见性，D为质心投影点G_proj到凸包边e_i的最小距离，r为以质心投影点G_proj为圆心的圆半径，SM(S，R)表示支撑结构S在重力f下所得到的von Mise应力，该应力也是支撑结构半径R的函数，σ_t为导致自支撑结构断裂的最小von Mise应力值。

所述的方法，所述的步骤五中，通过动态规划算法求解自支撑结构最小可见性候选集包括以下步骤：

1)累计最小可见性值计算

给定步骤三得到的自支撑结构可见性场，对该空间按照z方向进行分层，从支撑平面出发计算每一层的最小可见性值，然后通过累计每一层最小可见性值，递推每一层的最小可见性值：

M(i，j，k)＝v(i，j，k)+min(M(i+1，j-1，k-1)，…，M(i+1，j，k)，…，M(i+1，j+1，k+1))

直到递推至最顶层(具有最大z值)为止，其中，M(i，j，k)表示累计最小可见性值，v(i，j，k)表示离散边界空间中某一点的可见性值，i表示高度上的分层的当前层数，j表示宽度上的分层的当前层数，k表示长度上的分层的当前层数；

2)回溯寻找最小可见性值的自支撑结构中心线

获取到累计最小可见性值后，通过回溯的方法，从最顶层开始回溯跟踪具有累计最小可见性的支撑柱的中心线：即在自上而下的最小可见性自支撑结构跟踪过程中，给定当前层跟踪点，反过来确定该点在下一层相邻九个点位置中最小的累计可见性值的点，最终回溯至支撑面位置，从而得到支撑结构的中心线。综合支撑平面上所有离散位置所产生的支撑结构的中心线，并记录其累计可见性值，以此作为支撑结构最小可见性候选集。

所述的方法，所述的步骤六中，计算自支撑结构位置包括以下步骤：

1)在支撑结构最小可见性候选集中选择具有最小可见性支撑柱的支撑点作为全局可见性值最小点g，并以此为支撑结构的第一个支撑点g₀；然后，计算自支撑结构的后续支撑点位置的可行域；根据全局可见性值最小点g和地质体模型M的质心在平面上的投影点G_proj，建立最小可见性即全局最小点与质心投影点G_proj组成的直线gG的参数方程：

gG＝g+t(G_proj-g)/||G_proj-g||，

2)建立圆心为G半径为r的约束圆，则以g为起点到该圆的两条切线l₁和l₂的参数方程为：

利用切线l₁和l₂以及支撑平面矩形区域的边界，圈定包含约束圆的封闭区域，则该区域以外的区域即为可行域U₁；在可行域U₁中，从自支撑结构最小可见性候选集中选择具有最小可见性支撑柱的支撑点作为第二个支撑点g₁，求得g₁与给定圆的切线l₃和l₄的切线方程为：

根据得到的切线l₃和l₄以及支撑平面矩形区域的边界圈定封闭区域，以该封闭区域外与U₁相交的部分作为第三个支撑点的可行域U₂；在可行域U₂中，从自支撑结构最小可见性候选集中选择具有最小可见性支撑柱的支撑点作为第三个支撑点g₂；最终，以g₀，g₁和g₂分别作为最优支撑结构3条支撑柱的所在位置；

在上述从可行域U₁和U₂中搜索最小可见性支撑柱过程中，采用四叉树索引加分支限界的方法实现快速搜索，步骤如下：

1)四叉树索引构建

首先，在支撑平面上，划定包含支撑结构最小可见性候选集中所有支撑结构中心线支撑点的矩形域。然后，将矩形域作为当前四叉树的根结点，将当前根结点进行四分，得到四个子结点，并依次递归每个子结点，最终完成整个区域的四叉树索引构建；在构建四叉树过程中，在每个结点中存储自支撑结构最小可见性候选集在当前区域中的最小可见性值；

2)分支限界算法加速问题求解

利用分支限界算法求解最小可见性值问题，其中在进行支撑点g₁和g₂搜索时，预先设定一个较小值，按至顶向下的顺序访问落入可行域U₁和U₂中的四叉树结点，在访问过程中，如果当前结点存储的最小可见性值没有比预设的值更小，则直接终止对其子结点所对应区域的搜索，从而实现搜索的剪枝。

所述的方法，所述的步骤七中，计算自支撑结构半径包括以下步骤：

1)给定步骤六得到的最优支撑结构的所在位置，采用初始支撑结构半径R₀，建立基本支撑结构；对基本支撑结构进行四面体剖分，记录四面体剖分后的四面体体元总数total；

2)执行有限元分析，得到每个四面体单元的应力分布；

3)计算每个四面体单元的von Mise应力，并判断超出给定阈值σ_t四面体单元的个数n，计算n与四面体总数的比率ratio＝n/total；

4)如果ratio＞0或者迭代次数小于设定值，执行R＝R+α，搜索满足结构强度要求的支撑结构半径R，其中α表示每次迭代对R的递增步长；如果迭代次数超过指定值或n小于预设的最大个数，则结束支撑半径R的搜素，否则转步骤2)。

本发明的技术效果在于，在通过加设支撑结构来展示地质体模型的原始产状的基础上，优化了支撑结构的可见性和平衡性，并通过优化支撑结构的半径提升了支撑结构的结构强度。通过研究地质体在3D打印中的平衡问题，实现地质体模型的自支撑，对于提升地质体3D打印模型的实用性，具有很好的现实意义，同时对于3D打印技术在地质中的应用具有重要的意义。

下面结合附图对本发明作进一步说明。

附图说明

图1为通过动态规划算法求解式(10)得到的支撑结构的中心线；其中(a)为I矿区计算得到的支撑柱中心线；(b)为II-2矿区计算得到的支撑柱中心线；(c)为III矿区计算得到的支撑柱中心线；(d)为IV矿区计算得到的支撑柱中心线；

图2为通过动态规划算法求解式(12a)得到的支撑结构的中心线，其中(a)为参数β取1，(b)为参数β取2；

图3为带约束下的可行域U₁和U₂对自支撑结构位置求解示意图；

图4为IV矿区四叉树索引构建示意图；

图5为给定圆半径r＝35得到的最小的可见性值的平面视点分布示意图，其中三角形为组成的地面凸包，圆表示圆心在地质体质心投影点，半径为r的圆，点表示为自支撑结构候选集在支撑平面上的(支撑)点；左上角为I矿区平面视点，左下角为III矿区平面视点，右上角为II-2矿区平面视点，右下角为IV矿区平面视点；

图6为实施例中某矿床岩体三维模型图；

图7为实施例中模型的观察半球上采样点及可见性场示意图；

图8为实施例中地质体的潜在支撑结构的集合及其可见性值示意图；

图9为实施例中地质体的在支撑底面上的潜在支撑点示意图；

图10为实施例中地质体表达支撑点及其所在支撑柱可见性的四叉树结构图；

图11为实施例中地质体四叉树索引构建示意图；

图12为实施例中地质体支撑柱半径的调整示意图；至上向下分别为支撑柱半径R＝7，14，21时的模型应力分布；

图13为实施例中地质体模型最终支撑结构及打印结果图。

具体实施方式

本发明包括如下步骤：

①建立自支撑结构平衡性度量

通过添加外支撑结构来保证地质体模型的站立，解决地质体模型的平衡性问题。为了保证打印后地质体模型的平衡性，地质体模型M的质心G的投影G_proj应当落在模型与地面点接触部分所组成的凸包H的内部。因此，约束条件如下：

G_proj∈H， (1)

在实现中可分为如下两步：

1)求解质心

在求解三角网格模型的质心前，对当前的三角网格模型进行四面体化，可以得到模型的每个四面体

然后可以根据四面体化后的三角网格模型获取到每个四面体的质量

和对应的质心

通过得到的每个四面体的质心和它们的质量，得到当前地质体模型G的质心：

通过垂直投影G到所在的平面即可得到G_proj。

2)计算G_proj与凸包的关系(平衡约束)

假设凸包H有k个顶点{P_i}(i＝0，1，2，...，n-1)，它们按照顺时针方向排列于所在平面上。那么当前的平衡约束可以描述如下：

(P_i-G_proj)(P_(i+1)％k-G_proj)＜0，i＝0，1，…n-1 (3)

其中{P_i}为凸包H的顶点，即模型与支撑平面接触点，k为顶点的数量。通过寻找凸包H的边界上距离G_proj的最近点P*，计算G_proj到P*的距离d_G。如果G_proj落在凸包H的内部，则质心投影与凸包边界的最近距离为d_G，否则，令质心投影与凸包边界的最近距离为-d_G。

②建立自支撑结构强度度量

有限元法用来计算结构的刚度问题，它的刚度方程如下：

Kd＝f， (4)

其中K表示刚度矩阵，d表示位置偏移量，在上面的代数方程组中，d即为我们所要求解的变量。等式右边的f为我们所施加的外力大小。f可以指定在任意位置上或三角网格的面上，在后面的小节中将会根据具体应用来设定不同的应力大小。

③建立自支撑结构可见性度量

首先计算两点之间的可见性，如公式(7)。然后通过两点之间的可见性计算结果，计算任意一点在空间中的可见性结果，随后对将所有点的可见性结果进行积分(公式(8))，并将公式(8)中的积分结果转换为相应的期望表达(公式(9))，得到了最终的可见性场。具体过程如下：

给定当前的观察空间S_p和地质体块体空间M_V，对于给定一目标点x(x∈M_V)和人所观察的位置y(y∈S_p)，那么定义x和y的可见性函数V(x，y)如下：

其中值为0时表示不可见，1表示可见。

将观察点集合看成是上半球空间，在给定目标点集合后，半球空间对目标点集合的可见性可以看成是每个观察点对当前目标点的期望的积分，所以有如下公式：

E(V_x)＝∫_Ωp(y)V(x，y)dy， (8)

其中p(y)表示视线在半球上的分布，E(V_x)表示目标点x处可见性V_x的期望。

将公式(8)离散化为每个观察点对当前目标点x的求和形式：

E(V_x)＝∑V(x，y)

通过计算空间中每点的可见性期望值，得到地质体所在空间M_V的可见性场，则M_V中某条支撑结构S的可见性V(S)表示为：

通过计算地质体所在空间中每点的可见性期望值，得到空间的可见性场，在后续模型中用于引导支撑结构的定位与结构展布。

④自支撑结构自动生成最优化模型

给定的地质体模型，为保证公式(1)的平衡性约束，以及应满足公式(4)所设计的刚度平衡约束，构建如下优化问题：

其中，SM(S，R)表示支撑结构S在给定外力f下所得到的von Mise应力，这里SM(S，f)表示四面体网格中的每一个四面体单元所得到的von Mise应力

的集合。在公式(11)中的优化问题基础上，通过约束质心投影G_proj到凸包边界的距离，以保证质心投影G_proj与凸包的边界尽可能远。

给定地质体模型的质心投影点G_proj，以质心投影点为圆心，给定半径r，以及投影平面上的凸包H中的边e_i，i＝1，2，3…n，首先计算G_proj到凸包边e_i的最小距离D：

D＝min||G_proj-C_i||，

其中，C_i表示G_proj在每个边e_i垂直投影点。其次，通过计算得到质心投影到凸包边的最短距离D，为保证当前半径的圆在凸包H的内部，需要满足下面的约束：

r≤D， (11)

因此，式(11)的优化问题转变为：

④自支撑结构最小可见性候选集的求解

首先，将式(12)中的优化退化为对于不考虑约束下的极小化可见性场的问题(即式(12a))，利用动态规划算法进行求解。通过动态规划算法求解最小可见性值集合主要分为以下几个步骤：

1)可见性场计算

此步骤实现支撑结构的可见性计算通过求和支撑结构每个结点的可见性值，最终实现支撑结构的生成。在计算可见性场的过程中，主要应用了光线投射算法计算目标点到观察点之间的可见性结果，得到离散化的地质体边界空间的可见性场，这里约束支撑结构不可能在地质体模型的内部，这样避免了支撑结构从地面向地质体模型的内部延伸或穿过地质体模型的内部。

2)累计最小可见性值计算

给定地质体边界空间，对该边界空间按照z方向进行分层，假设边界空间最小z值为地面点，从地面出发计算每一层的最小可见性值，然后累计每一层最小可见性值，计算所有层的最小可见性值。下面是动态规划算法求解式(10)最小可见性值的递推公式：

M(i，j，k)＝v(i，j，k)+min(M(i+1，j-1，k-1)，…，M(i+1，j，k)，…，M(i+1，j+1，k+1)) (13)

这里，M(i，j，k)表示累计最小可见性值，v(i，j，k)表示离散边界空间中某一点的可见性值，由上面的可见性场得到，后面的公式表示最小化当前点上一层的相邻九个点的最小可见性值，其中i表示高度上的分层的当前层数，j表示宽度上的分层的当前层数，k表示长度上的分层的当前层数。这里的公式表示从最高层数的下一层开始计算每一层的最小可见性值。

3)回溯寻找最小可见性值集合

在上一步中，通过计算每一层中当前点的上一层中相邻九个点的最小可见性值，从而获取了每个点的累计最小可见性值M(i，j，k)。获取到累计最小可见性值后，通过回溯的方法得到每一层中的累计最小可见性值，最终确定支撑结构的中心线。具体方法是将自上而下的累计可见性值计算过程反过来从第一层开始确定第一层中的累计最小可见性值，利用之前计算方法，得到该点的上一层相邻九个点位置的最小的累计最小可见性值，最终回溯到边界空间的顶部位置。图1和图2分别给出了通过动态规划算法求解式(10)和式(12a)得到的支撑柱中心线集合。

⑥自支撑结构位置的求解

首先，计算自支撑结构位置的可行域，给定全局可见性值最小点g和给定的地质体模型M的质心在平面上的投影点G_proj，首先需要求解全局最小点(最小可见性值对应的视点，后面简称为全局最小点)与质心投影点G_proj组成的直线gG的参数方程。这里定义直线方程为：

gG＝g+t(G_proj-g)/||G_proj-g||， (15)

通过直线gG，在给定约束圆的半径r，求解以全局最小点为起点到该圆的两条切线。首先计算两条切线t₁和t₂与gG的夹角：

利用旋转矩阵可以得到gG沿顺时针和逆时针旋转θ之后，得到的直线l₁和l₂。这里求解结果点的x值的旋转矩阵为：

求解结果点的y值的旋转矩阵为：

因此，通过旋转得到的两个结果点分别记为res₁和res₂，所以直线l₁和l₂的参数方程为：

利用切线l₁和l₂，可以圈定矩形区域中包含给定圆的范围，根据r≤D约束，可以知道在两条切线所包含的区域外即为可行域U₁。同理，在可行域U₁中，选择第二点g₁，利用上面的方法可求得g₁与给定圆的切线l₃和l₄，它们的切线方程为：

根据得到的切线l₃和l₄，可以确定第三点的可行域U₂。

采用四叉树索引加分支限界的方法快速搜索，步骤如下：

1)四叉树索引构建

首先，在支撑平面上，划定包含支撑结构最小可见性候选集中所有支撑结构中心线支撑点的矩形域。然后，将矩形域作为当前四叉树的根结点，利用树的递归特性，可以将当前根结点进行四分，得到四个子结点，并依次递归每个子结点，最终完成整个区域的四叉树索引构建。在构建四叉树过程中，在每个根结点中存储了当前区域最小的可见性值点，这样在进行第二点第三点的查找时，直接通过判断当前根结点存储的局部最小可见性值点即可判断是否需要进一步查找，从而节约了递归查找的时间。图4给出了IV矿区四叉树索引构建结果。

2)分支限界算法加速问题求解

利用分支限界算法求解最小可见性值问题时，可以在确定第一点后，在进行第二点和第三点搜索时，可以预先假定一个较小值，如果当前值没有比假定的值更小，那么后面四叉树结点的查找直接终止，从而实现剪枝的效果。

利用上述四叉树索引加分支限界算法可以从整个搜索域中求解最小可见性值和G。图5给出了不同地质体模型通过给定半径r后，利用搜索算法得到的最小可见性值的平面视点的分布。从图5可以发现，通过添加一定半径圆之后，质心到凸包边界的距离在一定程度上增加了，地质体模型的平衡性都得到了相应的改善。

图5给定圆半径r＝35得到的最小的可见性值的平面视点分布。图中黄色线段为组成的地面凸包，绿色圆表示给定半径的圆，白色为平面所有视点。左上角为I矿区平面视点，左下角为III矿区平面视点，右上角为II-2矿区平面视点，右下角为IV矿区平面视点。

⑦自支撑结构半径的求解

在求解式中的vonMise应力前，对地质体三维模型进行四面体化。四面体化后，应用开源有限元方法库(OOFEM)进行应力场的获取。

在进行支撑结构半径优化过程中，主要有下面几步：

1)根据3个支撑点g₀，g₁，g₂及其在支撑结构最小可见性候选集中的中心线，采用初始支撑结构半径R₀，建立初始支撑结构；对初始支撑结构。对初始支撑结构进行四面体剖分，记录四面体剖分后的四面体体元总数total；

2)执行有限元分析，得到每个四面体单元的应力分布；

3)计算每个四面体单元的von Mise应力，并判断超出给定阈值σ_t四面体单元的个数n，计算n与四面体总数的比率ratio＝n/total；

4)如果ratio＞0或者迭代次数小于设定值，执行R＝R+α。然后，重复上述2)-3)步，直到迭代次数超过指定值或n值小于预设值为止。

以下4个步骤描述地质体3D打印自支撑结构自动生成方法的实施，以某矿床岩体(如图6)为例。

步骤一：以地质体包围盒底部中心为球心，设置地质体观察半球及观察半球上的采样点，如图7所示，利用公式(9)，通过光线投射算法，计算得到地质体自支撑结构空间的可见性场。

步骤二：利用动态规划算法，极小化支撑结构可见性

计算得到地质体的潜在支撑结构的集合及其可见性值(如图8)。根据动态规划算法得到的累计最小可见性值，将每根支撑结构中心线上的顶点的可见性值累加到底面点上，得到支撑底面上的潜在支撑结构支撑点集合(如图9)。

步骤三，给定地质体的在支撑底面上的潜在支撑点，构建表达支撑点及其所在支撑柱可见性的四叉树结构(如图10)。然后在支撑底面二维空间上，寻找最优的3个潜在支撑点及其支撑柱g₁，g₂，g₃：

G＝min(g₁+g₂+g₃)， (16)

以此作为最优支撑结构(如图11)。

步骤四，基于有限元模拟结果，调整最优支撑点所对应支撑柱的半径，以保证自支撑结构的强度(如图12)，得到最终的实现3D打印的地质体自支撑结构(如图13)。

Claims (9)

1.一种地质体3D打印自支撑结构自动生成方法，其特征在于,包括以下步骤：

步骤一，建立自支撑结构平衡性的度量：

给定地质体三维模型的任意自支撑结构模型，以地质体模型的质心投影落在模型与支撑平面接触点所组成的凸包的内部为目的，先计算地质体三维模型的质心位置，然后以质心投影与凸包边界的最近距离作为此自支撑结构平衡性的度量；

步骤二，建立自支撑结构强度的度量：

给定地质体三维模型的任意自支撑结构模型，通过有限元方法计算模型在地质体3D打印模型和自身重力作用下的von Mise应力，作为此自支撑结构在重力下结构强度的度量；

步骤三，计算反映自支撑结构隐藏程度的可见性场：

基于光线投射算法计算目标点到观察点之间的可见性，通过计算地质体模型所在的三维栅格空间在观察范围内可见性的期望值，表达空间中任意点在地质体遮挡下的可见性，得到反映自支撑结构隐藏程度的可见性场；

步骤四，建立自支撑结构最优化数学模型：

基于步骤一至步骤三分别得到的自支撑结构平衡性、结构强度和可见性度量，综合构建自支撑结构的最优化数学模型；

步骤五，建立自支撑结构的最小可见性候选集：

基于步骤四所建立的自支撑机构最优化模型，采用动态规划方法，通过在三维可见性场中逐层递推的策略，计算得到用于放置地质体3D打印模型的支撑平面任意一点处，满足自支撑结构最优化模型中最优可见性的自支撑结构，对支撑平面上所有点按上述方法计算最小可见性支撑柱，综合建立自支撑结构的最小可见性候选集；

步骤六，计算自支撑结构位置：

基于步骤五所得到的自支撑结构最小可见性候选集，根据自支撑结构最优化模型中的平衡性约束，采用分支限界策略，在四叉树索引加速下，从最小可见性候选集中搜索获得最优支撑结构所在位置；

步骤七，计算自支撑结构半径：

根据步骤六得到的最优支撑结构所在位置，搜索自支撑结构的半径，以满足自支撑结构最优化模型中的结构强度为约束，根据搜索得到的支撑结构半径，最终建立地质体3D打印自支撑结构三维模型。

2.根据权利要求1所述的方法，其特征在于，所述的步骤一中，首先对地质体三维模型进行四面体网格剖分，然后采用以下方式计算地质体三维模型的质心位置：

其中，n为四面体单元的个数，

为四面体的质量，

为对应的四面体质心，T_i为对当前的地质体三维模型进行四面体剖分后得到的四面体。

3.根据权利要求2所述的方法，其特征在于，所述的步骤一中，计算质心投影与凸包边界的最近距离包括以下步骤：

计算质心G到所支撑平面上的投影G_proj，判断G_proj是否落在模型与支撑平面的接触点所组成的凸包的内部，即G_proj∈H，其中H为凸包，则有：

(P_i-G_proj)(P_(i+1)％k-G_proj)＜0，i＝0,1,…n-1

其中{P_i}为凸包H的顶点，即模型与支撑平面接触点，k为顶点的数量；寻找凸包H的边界上距离G_proj的最近点P*，计算G_proj到P*的距离d_G；如果G_proj落在凸包H的内部，则质心投影与凸包边界的最近距离为d_G，否则，令质心投影与凸包边界的最近距离为-d_G。

4.根据权利要求1所述的方法，其特征在于，所述的步骤二中，利用有限元方法来计算在地质体3D打印模型及其支撑结构在自身重力作用下的von Mise应力包括以下步骤：

1)对地质体三维模型及其支撑结构模型进行四面体网格剖分；

2)通过求解刚度方程，得到自支撑结构在自身重力作用下的位移：

Kd＝f

其中K表示有限元刚度矩阵，d表示位移向量，f为表示模型自身重力的向量；

3)基于支撑结构的位移，计算支撑结构的内应力σ：

σ＝Cε

其中，C为表达3D打印材料力学性质的弹性矩阵，

表示梯度算子；

4)对模型剖分得到的每个四面体单元，根据其所对应的内应力σ^e,计算单元的von Mise应力：

5.根据权利要求1所述的方法，其特征在于，所述的步骤三包括以下步骤：

对于地质体所在三维空间M_V内任一目标点x和任一观察的位置即视点y，有x和y的可见性函数V(x,y)如下：

其中值为0时表示不可见，1表示可见；

以观察点集合作为上半球空间，则半球空间对目标点x的可见性V_x的期望通过下式表达：

其中p(y)表示在半球上的分布；

将公式(1)离散化为每个观察点对当前目标点x的求和形式：

E(V_x)＝∑V(x,y)

通过计算空间中每点的可见性期望值，得到地质体所在空间M_V的可见性场，则M_V中支撑结构S的可见性V(S)表示为：

6.根据权利要求1所述的方法，其特征在于，所述的步骤四中，最优化模型为：

其中，V(S)为根据步骤三计算得到的自支撑结构的可见性，D为质心投影点G_proj到凸包边e_i的最小距离，r为以质心投影点G_proj为圆心的圆半径，SM(S,R)表示支撑结构S在重力f下所得到的von Mise应力，该应力也是支撑结构半径R的函数，σ_t为导致自支撑结构断裂的最小von Mise应力值。

7.根据权利要求6所述的方法，其特征在于，所述的步骤五中，通过动态规划算法求解自支撑结构最小可见性候选集包括以下步骤：

1)累计最小可见性值计算

给定步骤三得到的自支撑结构可见性场，对该空间按照z方向进行分层，从支撑平面出发计算每一层的最小可见性值，然后通过累计每一层最小可见性值，递推每一层的最小可见性值：

M(i,j,k)＝v(i,j,k)+min(M(i+1,j-1,k-1),…,M(i+1,j,k),…,M(i+1,j+1,k+1))直到递推至最顶层为止，其中，M(i,j,k)表示累计最小可见性值，v(i,j,k)表示离散边界空间中某一点的可见性值，i表示高度上的分层的当前层数，j表示宽度上的分层的当前层数，k表示长度上的分层的当前层数；

2)回溯寻找最小可见性值的自支撑结构中心线

获取到累计最小可见性值后，通过回溯的方法，从最顶层开始回溯跟踪具有累计最小可见性的支撑柱的中心线：即在自上而下的最小可见性自支撑结构跟踪过程中，给定当前层跟踪点，反过来确定该点在下一层相邻九个点位置中最小的累计可见性值的点，最终回溯至支撑面位置，从而得到支撑结构的中心线；综合支撑平面上所有离散位置所产生的支撑结构的中心线，并记录其累计可见性值，以此作为支撑结构最小可见性候选集。

8.根据权利要求1所述的方法，其特征在于，所述的步骤六中，计算自支撑结构位置包括以下步骤：

1)在支撑结构最小可见性候选集中选择具有最小可见性支撑柱的支撑点作为全局可见性值最小点g，并以此为支撑结构的第一个支撑点g₀；然后，计算自支撑结构的后续支撑点位置的可行域；根据全局可见性值最小点g和地质体模型M的质心在平面上的投影点G_proj，建立最小可见性即全局最小点与质心投影点G_proj组成的直线gG的参数方程：

gG(t)＝g+t(G_proj-g)/||G_proj-g||，

2)建立圆心为G半径为r的约束圆，则以g为起点到该圆的两条切线l₁和l₂的参数方程为：

利用切线l₁和l₂以及支撑平面矩形区域的边界，圈定包含约束圆的封闭区域，则该区域以外的区域即为可行域U₁；在可行域U₁中，从自支撑结构最小可见性候选集中选择具有最小可见性支撑柱的支撑点作为第二个支撑点g₁，求得g₁与给定圆的切线l₃和l₄的切线方程为：

根据得到的切线l₃和l₄以及支撑平面矩形区域的边界圈定封闭区域，以该封闭区域外与U₁相交的部分作为第三个支撑点的可行域U₂；在可行域U₂中，从自支撑结构最小可见性候选集中选择具有最小可见性支撑柱的支撑点作为第三个支撑点g₂；最终，以g₀,g₁和g₂分别作为最优支撑结构3条支撑柱的所在位置；

在上述从可行域U₁和U₂中搜索最小可见性支撑柱过程中，采用四叉树索引加分支限界的方法实现快速搜索，步骤如下：

1)四叉树索引构建

首先，在支撑平面上，划定包含支撑结构最小可见性候选集中所有支撑结构中心线支撑点的矩形域；然后，将矩形域作为当前四叉树的根结点，将当前根结点进行四分，得到四个子结点，并依次递归每个子结点，最终完成整个区域的四叉树索引构建；在构建四叉树过程中，在每个结点中存储自支撑结构最小可见性候选集在当前区域中的最小可见性值；

2)分支限界算法加速问题求解

利用分支限界算法求解最小可见性值问题，其中在进行支撑点g₁和g₂搜索时，预先设定一个较小值，按至顶向下的顺序访问落入可行域U₁和U₂中的四叉树结点，在访问过程中，如果当前结点存储的最小可见性值没有比预设的值更小，则直接终止对其子结点所对应区域的搜索，从而实现搜索的剪枝。

9.根据权利要求1所述的方法，其特征在于，所述的步骤七中，计算自支撑结构半径包括以下步骤：

1)给定步骤六得到的最优支撑结构的所在位置，采用初始支撑结构半径R₀，建立基本支撑结构；对基本支撑结构进行四面体剖分，记录四面体剖分后的四面体体元总数total；

2)执行有限元分析，得到每个四面体单元的应力分布；

3)计算每个四面体单元的von Mise应力，并判断超出给定阈值σ_t四面体单元的个数n，计算n与四面体总数的比率ratio＝n/total；

4)如果ratio＞0或者迭代次数小于设定值，执行R＝R+α，搜索满足结构强度要求的支撑结构半径R，其中α表示每次迭代对R的递增步长；如果迭代次数超过指定值或n小于预设的最大个数，则结束支撑半径R的搜素，否则转步骤2)。

Images