CN111125942A

CN111125942A - 用于三维单元结构建模和拓扑优化的b样条高清晰度单元水平集方法和计算机存储介质

Info

Publication number: CN111125942A
Application number: CN201911031400.8A
Authority: CN
Inventors: 王煜
Original assignee: Hong Kong University of Science and Technology HKUST
Current assignee: Hong Kong University of Science and Technology HKUST
Priority date: 2018-10-31
Filing date: 2019-10-28
Publication date: 2020-05-08
Anticipated expiration: 2039-10-28
Also published as: CN111125942B; US20200134918A1; US10984593B2

Abstract

本发明提供了一种用于三维单元结构建模和拓扑优化的B样条高清晰度单元水平集方法和计算机介质，该方法包括：将目标全局结构的全局设计域划分为多个子域单元；针对所述多个子域单元中的每一个，使用隐式B样条函数定义用于表示子域单元的结构材料域的单元水平集函数，使得整组单元水平集函数包括整个全局设计域上的全局水平集函数。

Description

用于三维单元结构建模和拓扑优化的B样条高清晰度单元水平集方法和计算机存储介质

技术领域

本发明涉及拓扑技术领域，尤其涉及一种用于三维单元结构建模和拓扑优化的B样条高清晰度单元水平集方法和计算机存储介质。

背景技术

拓扑优化是一种数学和计算方法，最初是为优化材料分布而开发的，目的是在给定的设计领域、目标函数和约束条件下获得更好的实体结构的结构性能。经过近三十年的发展，拓扑优化已获得相当大的普及，并且该技术已成为可行的并且在许多领域中不可或缺的工具，例如用于轻型结构、多孔材料(cellular material)和医疗植入物的设计和优化的工程实践。从早期的均质化方法到各向同性固体材料惩罚方法(SIMP)、演化结构优化方法和水平集方法，已经提出了各种建模和优化方法。

用于形状和拓扑优化的水平集方法基于材料域和界面边界的隐式表示。在水平集方法中，结构界面由高维水平集函数的零水平隐式表示，高维水平集函数的零水平通过水平集函数的动态变化自然地适应材料域中的拓扑变化。在最一般的情况下，结构形状和拓扑变化在数学上由哈密顿雅可比偏微分方程控制。通过伴随灵敏度分析，可以计算形状导数并将其用于构建优化算法，该优化算法用于演化水平集函数并获得优化解。与基于材料密度的拓扑优化方法(例如SIMP)相比，该拓扑优化方法表示在一组固定的元素或体素上具有构造的密度变量的结构域，该水平集方法的主要特征是材料域之间的几何边界用指定的平滑度精确表示。几何表示的性质(尽管是隐式的)提供了与计算机辅助设计(CAD)系统和用于优化结构的生产和制造的过程设计直接集成的可能性。

在经典的水平集方法中，水平集函数是通过在整个设计域中在三个维度上定义的单个函数来全局定义的。最常用的函数是有符号距离函数，该函数几乎在任何地方都可以微分，并且其梯度满足程函方程。计算符号距离函数的算法包括高效的快速行进方法，快速扫掠方法和更通用的PDE方法(通常在计算机视觉中使用)。可替代地，可以用参数来表示水平集函数。最常选择三种类型的基函数来进行参数化：径向基函数(RBF)或紧凑支持的径向基函数(CSRBF)、B样条或NURBS插值函数、和作为增加插值的稀疏性的近似值的离散小波变换。然而，在用于全局水平集函数的这些参数化方案中，必须在整个设计域上构造密集或适度稀疏的全局插值矩阵。因此，当计算形状和拓扑优化的动态变化时，这些参数化方案将导致较高的计算成本。

拓扑优化技术也已被部署用于设计多孔材料和结构，多孔材料和结构也被称为建筑材料或超材料，用于自然界中通常不被发现的高度定制或极端物理特性，例如在低质量密度下的超高强度。多孔材料的计算设计中最方便的方法是使用相同单位晶胞的周期性阵列来定义材料中的微观结构。因此，设计问题被转换为优化单个单位晶胞设计空间内的材料分布。在所谓的逆均质化方法中，均质化理论被应用于近似材料的有效特性。尽管已经对周期性微结构的优化进行了深入研究并针对各种物理特性进行了研究，但具有空间变化的多孔结构和相关物理特性的优化的多孔材料和结构的潜力更大。当将单元结构设计为具有空间变化的微结构时，可以最好地实现并满足对材料或结构特性的局部要求。

单元结构的拓扑优化中的一个关键问题是在每个微结构可以具有自己的配置的同时，用于对可变微结构建模的表示方案。一种方便的方法是使用显式参数化预定义一类微结构，通常是基于桁架/梁或类似壳的结构。在给定的参数范围内，可以针对设计目标调整每个微结构，只要该微结构可以保持其连通性和完整性。因此，所得的单元结构在结构性能或在多功能方面只能具有有限的增益。

在设计具有空间可变微结构集合的单元结构时，至关重要的是要确保相邻单元之间的几何连通性。当简单地使用拓扑优化来优化每个单独的单元时，相邻的单元不必形成不可分割的部分，从而导致物理性能下降。任何未连接的单元都不能承受外部负载，并且会在生产过程中造成制造困难。解决这个关键问题的现有工作通常集中于在相邻的微结构单元上追寻几何连通性条件，例如施加几何约束、应用预定的连接器或采用顺序过滤。这些分段方法(piece-meal approach)实质上用大量的几何约束对设计空间施加了严格的限制，从而损害了优化设计的性能。因此，它们实质上降低了设计优化的能力(如果仍将其视为拓扑优化)，并且它们无法轻易扩展到三维单元结构。

并发优化方法在材料和结构的分层优化中提供了更大的灵活性。它基于尺度分离(scale separation)，其中均质化近似值用于计算可以单独优化的微结构的宏观等效机械性能。遗憾的是，由于假设尺度分离，所得到的优化结构通常无法制造。采用多尺度FE2模型来设计具有空间变化的微结构的结构，但是在优化的结构中，不能保证相邻微结构之间的几何连通性。对于现有的分层拓扑优化方法，这是一个共同的挑战。最近，提出了一种在单个单元上具有形状变态的水平集方法，以获得拓扑相似的、因此可连接的微结构，其连接面至少具有C0光滑度。这个想法在所谓的VCUT水平集方法中得到了进一步的推广。

在当前的商业CAD系统中，普遍使用边界表示(B-rep)。通常，3D实体目标的2个流形边界由一组(修剪的)表面表示，由此描绘了其体积。尽管这种B-rep方案已成为标准，但已广泛认可3D目标的完整体积表示将在整个设计周期中作为建模、分析和优化的单个几何模型提供较大的潜在优势。最近，在三维参数空间中，提出了显式B样条三变量模型，作为对传统参数B脊线和曲面的体积扩展。作为显式模型，此类B样条三变量仅限于长方体拓扑，不能表示一般的体积形状。它可以是定义体积的构造块，在其域内仅具有一组修剪的B样条曲面。因此，显式体积模型需要广泛使用布尔运算来构建实体模型。在此处提出的HD-CLIBS表示中，隐式描述了体积模型。

发明内容

在本公开中，目标是开发一个通用的几何模型，用于具有空间变化的微结构以及优化的三维(3D)形状和拓扑的一般实体和/或单元结构的拓扑优化。这个模型自然地保证了在任何相邻的微结构单元之间的边界上具有定义的平滑度的几何连通性。使用该模型，拓扑优化过程同时优化结构的实体域以及包括单元结构的所有结构单元的物理属性。该设计模型被称为高清晰度B样条单元水平集(或简称HD-CLIBS)，因为采用了水平集方法的基本概念，而其应用被推广到一组子域单元上，这些子域单元分割了正在优化的结构的整个设计域，潜在地产生了高清晰度单元/实体结构。

首先，将优化中的全局结构的整个设计域划分为的一组3D子域单元。在每个子域单元上，分别定义水平集函数，使得整组水平集函数包括整个全局域上的全局水平集函数。因此，子域单元上的这些水平集函数被称为单元水平集函数。这与通常的水平集方法形成对比，通常的水平集方法采用单个全局水平集函数捕获整个结构的内部和边界。

其次，在每个子域单元上，使用隐式B样条函数来定义用于表示单元的结构材料域的单元水平集函数。得益于隐式B样条的内在特性，这种参数表示对于单元水平的拓扑优化具有很大的灵活性。更重要的是，可以容易地定义相邻单元的B样条函数的参数系数的简单兼容性约束，从而保证它们的几何连通性达到指定阶的平滑度，例如一阶可微C¹。

第三，除了单元结构的建模之外，HD-CLIBS表示进一步被结合到拓扑优化方案中。基于形状导数，设计了最陡下降方案，用于更新B样条系数和演化所有单元的水平集函数，根据设计目标和指定约束优化所有单元的物理性质以及全局结构。借助于隐式B样条，递归地对B样条系数的优化变量进行更新，而不涉及线性代数系统的求解。因此，所得到的优化过程实质上是有效的，并且也容易在并行计算平台上实现。这里描述的HD-CLIBS方法通过在不同体积分数下优化体积模量的三维单元结构的设计得到了证实。

附图说明

图1示出了水平集表示方法，其中第2部分的结构由第1部分的水平集函数隐式表示。

图2示出了单元水平集表示方法，其中第1部分的全局水平集函数是子域水平集函数的集合，例如第3部分，全局结构是子域结构的集合，例如第4部分

图3示出了使用快速B样条插值的一般模型初始化过程，其中第5部分是输入结构几何形状，第6部分是相应符号距离函数的提取样本，第7部分是B样条基础函数，第8部分是重构水平集函数。

图4给出了有关初始化过程的2D示例。第9部分是由二进制图像表示的骨组织，可以通过扫描获得。第10部分是基于第9部分的边界计算出的符号距离函数。第11部分是通过快速插值算法计算出的隐式B样条函数，第12部分是由第11部分的隐式B样条函数的零水平集表示的复制骨组织。

图5给出了有关初始化过程的两个3D示例。第13部分和第15部分是原始骨结构，第14部分和第16部分是由隐式B样条函数的零水平集表示的重建结构。

图6示出了通过将第17部分和第19部分的两个目标直接连接成2D中的第18部分的单元结构的几何集合。

图7比较了连接两个目标时的几何连续性。第18部分示出了两个目标之间的不匹配，第20部分拥有C⁰连续性，第21部分拥有C²连续性。

图8示出了具有C²连续性的3D单元连接。通过迫使第22部分满足等式(16)来产生第23部分，通过迫使第25部分满足等式(16)来产生第26部分。第24部分是通过连接第23部分和第26部分产生的新单元。

图9示出了具有满足周期性条件的系数的3D单元。通过迫使第14部分满足周期性条件来生成第27部分。通过以2×2×2的方式复制第27部分来生成第28部分。

图10示出了结插入后B样条基函数的变化及其在单元分裂中的用途。

图11给出了关于单元分裂的2D示例，其中将第12部分的单元分裂为第30部分和第31部分，其中第29部分和第32部分是分裂后的隐式B样条函数。

图12给出了关于单元分裂的2D示例，其中第12部分的单元分裂为第34部分、第35部分、第38部分和第39部分，其中第33部分、第36部分、第37部分和第40部分是分裂后的隐式B样条函数。

图13给出了有关单元分裂的3D示例，其中将第14部分的单元划分为第22部分、第25部分、第41-46部分的8个较小部分。

图14示出了通过使用分割和连接操作从第47部分的均匀单元离散化到第48部分的分层单元离散化的变化。

图15示出了通过同时调节系数来进行腐蚀和膨胀操作。通过同时将第14部分的系数减小一定量来生成第49部分，并且通过同时将第14部分的系数增大一定量来生成第50部分。

图16示出了通过将第51部分的渐变图案覆盖到第14部分来生成第52部分的渐变单元。

图17示出了通过将第54部分的渐变图案覆盖到第53部分的原始均匀单元结构上来生成第55部分的渐变单元结构。

图18示出了对第56部分的单元的局部修改，以生成第57部分的新单元。

图19示出了混合操作。第58部分和第12部分是两个需要混合的单元。第59部分是通过对第58部分和第12部分的系数进行最小值运算生成的。第61部分是通过对第58部分的系数进行逆运算生成的。第62部分是通过对第58部分和第61部分的系数进行最小值运算生成的。第60部分是通过对第59部分和第62部分的系数进行最大值运算生成的。

图20示出了在第63部分的简化叶片和第64部分的松质组织上进行的混合操作，所生成的单元结构是第65部分。

图21示出了在第15部分的股骨和第28部分的松质组织上的混合操作。第66部分和第67部分是所产生的单元结构的两个半部分。

图22示出了基于HD-CLIBS模型的拓扑优化流程图。

图23示出了具有C0等式约束的HD-CLIBS模型的拓扑优化结果。第69部分是第68部分给出的边界和载荷条件的解，第71部分是第70部分给出的边界和载荷条件的解。

图24示出了在零件72中给出的桥梁荷载作用下零件73的解决方案。

图25给出了第74部分的边界和载荷条件，以及第75部分的初始设计，用于HD-CLIBS模型的拓扑优化。

图26示出了部分76的单比例设计，其中单元之间只有C0相等，而部分77是部分76的1/4。

图27示出了具有周期性条件的部分78的单元结构设计，而部分79是一个单个的周期性单元。

图28示出了具有部分80的单元格结构，在单元上分层平均，而部分81是部分80的1/4。

具体实施方式

A.高清晰度单元水平集模型

在经典的水平集方法中，一般的三维结构由单个全局水平集函数表示，如图1所示，其中该结构的几何模型由下式隐式获得：

其中x＝{x，y，z}∈R³表示设计域D内一个点的坐标；t表示水平集函数演化的伪时间；

表示结构的边界；Ω表示结构的材料域。通过使用全局水平集函数，结构边界的演化被隐式转换为水平集函数的零曲线(对于2D问题)或曲面(对于3D问题)的更新。在传统的水平集方法中，通过求解以下汉密尔顿-雅各比(Hamilton-Jacobi)偏微分方程来实现水平集函数的演化：

其中V＝dx/dt，Vⁿ是表示水平集动态演化的法线速度。符号距离函数最常用于函数φ，该函数几乎处处都可微分，并且其梯度满足程函方程(eikonal equation)。用于计算符号距离函数的有效算法包括快速行进方法和窄带方法，从而在整个设计域D上构造单个函数(图1)。

在本公开中，定义了单元水平集表示的概念，与使用单个全局水平集函数的常规水平集模式相反。首先，将全局设计域D和全局结构域Ω分别分割为M个子域单元D_s和微结构Ω_s(s＝1，…，M)，从而

其中，

是微结构的边界。如图2所示，在每个开放的有界子域单元D_s上，定义水平集函数φ_s来表示该单元内的微结构，即

根据水平集概念，子域水平集函数φ_s可以通过控制H-J方程在其子域单元D_s上进行动态演化：

其中

是子域水平集函数φ_s的向外法向速度场。因此，子域水平集函数的总集包括整个全局域D上的全局水平集函数。因此，子域单元上的这些水平集函数被称为单元水平集函数(cellular level set function)。利用这种设计域离散化，可以将子域单元的特征尺寸细化到超细规模，从而生成具有高清晰度和特别细节的设计。

B.高清晰度B样条单元水平集(HD-CLIBS)

在利用单个全局函数φ的经典水平集方法中，用离散方案求解Hamilton-Jacobi方程将遇到很多问题，例如重新初始化、速度扩展和数值稳定性的CFL条件。如果将水平集函数例如用径向基函数(RBF)参数化，则通常会导致在更新相关联的参数化参数的过程中求解线性系统，通常稀疏性有限，因此计算量大。为了解决使用基于RBF的参数化方法的大规模问题，提出了一种想法，即通过统一分区(POU)将全局设计域划分为一些重叠的局部子域。尽管POU技术降低了计算复杂性，但是在没有实际可行的指导原则的情况下，确定重叠区域的大小成为另一个问题。通常，单元水平集函数的概念可以结合任何可行的表示方案，例如经典的符号距离函数、基于反应扩散的模型和分段恒定水平集函数。但是，通过结合适当的参数化方案，发现可以充分利用单元水平集表示所提供的灵活性，来构建用于建模、设计和优化多孔材料和结构的通用且强大的平台。使用流行的基于参数建模的功能，可以在当前的CAD基础架构中自然地实现这种建模和优化方法。

在本公开中，我们选择将单元水平集函数φ_s定义为隐式B样条函数，它表示为B样条基函数的线性组合。对于一般的3D模型，x方向上次数p、y方向上次数q和z方向上次数r的张量积B样条函数是以下形式的三变量分段多项式函数：

其中

形成用于水平集函数的由(n+1)×(m+1)×(l+1)个B样条系数组成的数组，

是以下张量积形式的三变量B样条函数：

而

和

是分别对x、y、和z方向上的节点向量进行定义的B-样条基函数。B样条的节点向量位于横跨固定子域单元D_s的网格上，分别在三个方向中的每个方向上具有规则间距h。因此，用系数

作为参数来对单元水平集函数φ_s进行参数化，而其边界由子域单元上的水平集方程来描绘：

其中x∈Γ_s＝θΩ_s，s＝1，…，M(8)

隐式三变量B样条模型是通用的，因为系数

定义了用于子域单元上参数化建模和设计的变量。尽管这些系数通常不与网格的任何空间坐标相关联，但可以将它们视为所谓的格雷维尔横坐标(Grevilleabscissae)控制点处的φ_s的控制参数。

已经注意到隐式B样条函数具有多种吸引人的特性，包括紧凑的局部支持和其梯度范数的严格限制。在图像处理领域，隐式B样条函数由于其固有属性而被用作分割和配准的图像模型。例如，已经开发了一种快速的B样条插值公式，可以根据体积数据进行有效的曲面重建。此外，局部支持B样条基函数，并且可以通过节点操作来调整支持域。因此，容易适应用户交互以控制重构表面的特性，例如局部平滑度。

在本节的以下部分中，将讨论四个主要的建模方面，即模型初始化和重建、单元连接性、单元划分和布尔运算，以示出HD-CLIBS模型作为具有体积表示功能的建模系统的几个优点。

首先简要描述与HD-CLIBS模型相关的基础的B样条特性，因为B样条和NURBS已成为CAD系统中的标准表示。在不失一般性的前提下，讨论一维单位单元内x方向上次数p的单变量B样条函数，

其中c_i是B样条系数，基函数B_i，p(x)可递归形成为：

其中x_i是节点序列X＝{x₀，x₁，x₂，...，x_N}的第i个值。如果节点值是唯一的，则它是一个简单的节点。当仅使用简单的节点时，所有基函数将共享相同的形状。当节点值重复S次时，它是具有重复度为S的重复节点。例如，对于S＝3，节点序列为x₀＝x₁＝x₂＝0。此外，重复度不能高于p+1并且N＝n+p+1。由于在方程(7)的3D情况下的张量积，所有这些性质都推广到三变量B样条基函数N_i，j，k(x)。

在HD-CLIBS模型中，定义了一个直线网格，横跨固定子域单元D_s，分别在三个方向中的每个方向上具有规则间距h。B样条的节点向量位于网格上。在每个方向上的第一个和最后一个节点处，将节点分别设置为以最高重复度完全重复，例如在x方向上为p+1。为了保持两个相邻单元在其公共面处的水平集函数之间的连通性，而进行了这种节点重复度规范，从而保证相邻单元之间的单元结构中所需的几何平滑度。

1)模型的初始化与重构

对于HD-CLIBS模型，其初始构建是通过确定三变量B样条系数来初始化其水平集函数的步骤。有两种通常的方法可以完成构建。如果结构的一组体积数据(例如CT扫描数据)可用，则可以计算该结构的HD-CLIBS模型。这本质上是模型重建的问题，已经在几何建模和计算机图形学领域中对此进行了广泛研究。另一方面，如果结构的实体或表面表示形式(例如CAD模型)存在，则可以像传统的水平集表示那样为模型建立符号距离函数。基于符号距离函数的离散化，可以为HD-CLIBS模型构建准备离散的体积模型。

对于HD-CLIBS模型的构建，可以基于两种快速技术，利用B样条的特定优势作为水平集表示的基础。对于立方次数即p＝3以下的三变量B样条，可以使用快速B样条插值公式，作为首选方案，可以将该公式进一步扩展到包括重复节点的情况以用于快速重建。对于高阶B样条，可以通过快速离散B样条卷积运算获得B样条系数。此外，由于三变量B样条函数是三维的张量积，因此该卷积可以作为简单一维卷积的序列执行。因此，特别是与任何L²投影或搭配方法相比，这种构造过程在计算上是有效的。

图3示出了重新初始化和重建过程的流程图。首先，利用给定的CAD模型，使用均匀的Euler网格计算符号距离场。在同一网格上，构建了由重复节点所支持的B样条基函数，并使用快速B样条内插法计算了它们的系数，从而得到了插入符号距离场的隐式B样条函数。

图5示出了一块松质骨样品，其边界用三角形网格建模。利用上述模型重建过程，使用快速插值算法对100×100×100离散样本进行操作，构建了三变量B样条函数。图5还给出了基于构建的HD-CLIBS表示的再生松质骨。

2)单元连接

对于任何单元模型，一个基本要求是指定和维持所表示的结构任何相邻单元之间在其相邻面处的几何连通性的能力。在HD-CLIBS模型中，具有隐式B样条函数的体积表示的性质允许在任何相邻单元之间具有指定平滑度的几何连通性。通过限制与连接两个相邻单元的边界面相关的单元水平集函数φ_s(x)中的三元系数

可以轻松实现这一点。在不失一般性的前提下，讨论了具有节点序列X＝{x₀，x₁，x₂，...，x_N}的1D单位单元内x方向上次数p的单元B样条函数。为了简化讨论，以单变量p＝3次数为例。两个相邻的单元分别用符号+和-简单表示。在一个单元的起始边界节点

处，它是重复节点。类似地，对于要在相同位置处连接的相邻单元，其结束边界结节点

也被重复。通过使用三次B样条函数，可以证明，当B样条系数之间分别满足以下关系时，两个相邻的1D单元在其连接位置处的几何连通性可达到二阶可微分：

用于C⁰连续性 (12)

用于C¹连续性 (13)

用于C²连续性 (14)

对于任何其他阶的B样条，也可以很容易地概括出相邻单元之间的几何连通性的这些条件。将通过以下示例说明这种功能。

如图6所示，构造两个直线2D目标，分别具有其单元水平集的双变量函数(如第17部分和第19部分所示)，并且域单元的大小相同。左单元第17部分是通过给定2D结构的腐蚀操作而构造的，而右单元第19部分是来自相同结构的膨胀操作的结果。在各个水平集函数的零集中，嵌入在第17部分和第19部分中的两个目标沿着选定的一侧进行组装，以形成新的单元结构，如第18部分所示。显然，组装后的目标在两个单元之间的连接处具有几何不匹配(即，台阶)，如第18部分所示。

如本示例所示，可以通过将在方程(12)(13)(14)中定义的几何连通性条件强加到所需的B样条系数，来消除几何不匹配的几何不连续性。例如，单元水平集函数的双变量系数可以更改为

和

表示的平均值，使得

这将确保连接面上的C⁰连续性。此外，如果选择新系数

和

也要满足相同的要求，

则它们满足方程(12)(13)(14)的所有连续性条件，从而产生C²连续性。几何约束对最终连接的单元目标中的几何连续性的这些影响如图7所示，可以发现，对于C⁰连续性，图7第18部分中的几何不匹配在图7第20部分中得到了消除，对于C²连续性，图7第18部分中的几何不匹配在图7第21部分中得到了消除。应当注意，对于B样条系数，方程(12)(13)(14)的几何连续性条件可以以其他选择的形式来满足，并且这里使用简单的形式进行说明。

在下一个示例中，说明3D情况。为了便于说明，从图5的第14部分中选择了两个完全不连接的单元，一个在左上角，另一个在右下角，如图8的第22部分和第25部分所示。如果这两个单元直接并排连接，在这两个单元的连接面上可能会存在几何不匹配。同样，将如上面的示例中的用于C²连续性要求的约束条件方程(16)应用到所选的用于其连接的每个单元的各个面(图7的第21部分)上。由于HD-CLIBS模型B样条系数的变化，这两个单元的几何形状发生了变化，如图8的第23部分和第26部分所示。最后，这两个单元在选定的面连接，产生在相邻面之间具有C²连续性的组合单元结构(图8的第24部分)。

3)单元周期性

接下来，讨论如何对由周期性分布在结构中的单个单元组成的单元结构进行建模，正如在有关超材料的文献中所研究的许多情况一样。在HD-CLIBS模型中，当将方程(12)(13)(14)的几何连续性条件沿周期重复方向应用到单元的相对面时，很容易表示周期单元。值得一提的是，不仅提供了周期性，而且借助HD-CLIBS模型也很容易适应周期性连接中的更高阶连续性(例如C²)。

如图9所示，根据第14部分的松质骨CT扫描的体积数据对一个单元进行建模。分别沿三个方向中的每一个方向，将方程(16)的C²连通性条件应用于与相对面相关联的第一个节点和最后一个节点的B样条系数，简单地进行如下更改：

(c₀)^*＝(c₁)^*＝(c₂)^*＝(c_n)^*＝(c_n-1)^*＝(c_n-2)^*＝(c₀+c_n)/2 (17)

然后，将原始单元更改为第27部分的单元。将这个新单元重复成2×2×2的阵列的单元后，将构造一个新的周期单元结构作为第28部分，该结构在连接单元的边界处具有C²连续性。图9中结构的俯视图描绘了周期性结构公共边界上的平滑连接。该操作简单有效。

4)单元划分

在HD-CLIBS模型中，任何子域单元都可以分为两个或多个单元。使用三变量B样条函数，可以通过在节点序列中插入足够多的特定节点来实现此目的。当结重复p+1次时，效果是此单元域在相应的方向上被划分为两个子域，并且相应的单变量B样条中的所有基函数都分为两组。在细分的新子域中支持每组中的基函数，并且没有可以跨越插入的节点的基函数，如图10所示。这意味着，一组中的B-样条函数的系数只能控制该子域中B样条的形状，它们对另一个子域中的B样条的形状没有影响，反之亦然。在这方面，原始的B样条函数基本上被分为两个独立的B样条函数。使用这种技术，可以将一个单元划分为多个较小的单元，这可以为在几何图形表示上进行局部细化提供有效的工具。

基本单元划分过程如图10所示，为了简单起见，这里再次使用一维三次B样条曲线。在原始单元中，存在一系列系数c_i。然后，在中间插入一个新的节点，重复4次。因此，创建了一系列新的基函数，并生成了它们的新系数。所得的两个划分的新单元分别用符号+和-表示。如上图所示，给出针对左单元新创建的单元的B样条系数，

对于右单元

以类似的方式，并借助HD-CLIBS模型中三变量B样条基函数的张量积，可以轻松地扩展单元划分并将其应用于2D和3D情况，2D和3D情况分别对应c_i，j(在2D模式下)和c_i，j，k(在3D模式下)，具有相似的分配。

图11所示的2D示例用于说明。第12部分的原始单元在水平方向上均等地分成第30部分和第31部分两个半单元。根据单元划分算法，在新的细分单元中分别生成了两个新的隐式B样条函数，作为表示细分的2D目标的新HD-CLIBS模型，如第29部分和第32部分所示。随后可将该划分操作应用在另一方向(垂直方向)，这将导致将两个半单元被划分为四单元，如图12所示。

也可以使用该划分过程对3D单元进行分割。如图13所示，将3D小梁骨单元平均分为2×2×2立方体的八个单元(oct-cell)。在实践中，划分过程可以为内部结构检查提供一种方法，并且可以在不影响相邻区域的情况下方便进行局部几何修改。

除了物理上划分和连接单元外，划分和连接功能还为用户提供了有效而灵活的工具，以在逻辑上更改单元结构内部的单元布局和层次结构。如图14所示，将2D中的原始单元结构均匀地分成8×8个单元。通常，离散化在单元结构设计期间是固定的。现在，利用单元划分和连接操作，可以通过操纵节点和系数轻松调整单元布局，并且对于相同的单元结构，可以同时存在具有不同单元大小的分层离散化。

5)单元操作

HD-CLIBS模型的隐式三变量B样条函数实质上是水平集函数的参数化。通过调整单元的B样条函数的系数

可以操纵和修改结构。另外，由于基函数是局部支持的，因此系数的影响将限制在由节点区间和次数p确定的域中。在这里，讨论一些有用的操作，旨在证明HD-CLIBS模型的有效性。

首先，基于初始化的B样条函数，可以通过同时在系数上添加或减去正常数来轻松更改全局特征厚度。如图15所示，从系数中减去正常数将侵蚀特征并导致单元密度较低。相反，对系数增加正常数将扩大特征并增加单元密度。可替代地，可以以渐变(graded)的方式调整系数，以便实现单元内部的渐变的过渡。在图16中，内部渐变过渡单元是通过将第51部分的渐变图案添加到第14部分中的原始系数而获得的。

渐变图案也可以应用于由多个单元组成的单元结构。如图17所示，当将第54部分的渐变图案分配给由4乘以4单元构成的第53部分的单元结构时，生成第55部分的渐变单元结构，该渐变单元结构在每个单元之间完全连接。

除了在整个单元上进行操作外，还可以通过改变三变量系数来实现局部调整，因为每个基函数均受局部支持。如图18所示，如果想要扩大中心的特征，则仅应增加第56部分的有限区域中的系数。第57部分的结果表明，只有该区域中的特征被扩大，位于外部的特征保持不变。

此外，可以通过处理三变量系数将不同的单元融合在一起。如图19所示，给出了两个初始单元，包括模型系数为C₁的2D股骨样品和系数为C₂的松质骨样品。首先，通过对两组系数进行最小值操作(min operation)生成内部部分，

C_inside＝min(C₁，C₂) (20)

然后，通过另一个最小值操作生成边界皮肤部分，

C_bound＝min(C₁，-(C₁-r)) (21)

其中r是控制皮肤厚度的参数。最后，用最大值操作(max operation)将这两个部分结合起来，将得到

C_final＝max(C_inside，C_bound) (22)

使用C_final系数，保留模型1形状的混合单元由填充模型2的细节填充，并被具有由r控制的厚度的皮肤覆盖。

此混合操作可以轻松扩展到3D情况。在图20中，简化的叶片填充有由2×1×2周期性单元组成的骨组织。在图21中，使用真实的股骨样品进行说明，其中实体区域填充有由2×2×2个周期单元组成的骨组织。此外，如果需要更多细节，则可以使用更多的周期性单元构建骨骼样品。还可以根据实际需要局部调整内部微结构的相对密度。通过这种混合操作，很容易建立由复杂的下层微结构组成的单元结构。

C.使用HD-CLIBS优化单元结构的拓扑

本公开的主要动机是开发一种基于参数水平集的建模和优化方法，以克服经典水平集方法中的主要困难。上面描述的HD-CLIBS模型在分割全局设计域的一组单元子域上运行，而每个单元内的结构的材料域由其自身的水平集函数捕获。与在经典水平集方法中使用单个全局水平集函数相反，现在在每个子域单元上分别并行地执行水平集函数的动态演化。此外，在使用三变量B样条函数进行参数化的过程中，HD-CLIBS模型的内在属性可以在结构拓扑优化过程中得到广泛利用。因此，本节中描述的优化方法具有以下特征。

(1)HD-CLIBS参数模型使得能够将拓扑优化问题的解计算为将变分问题在B样条横跨的空间上的限制，换句话说，在B样条空间上的投影。因此，可以根据B样条系数的迭代独立地直接获得优化解，而无需求解H-J方程或具有插值矩阵的任何线性系统。

(2)快速B样条构造方案很容易部署用于计算B样条系数的灵敏度。这项技术将像在常规水平集方法中通常处理的那样，将灵敏度积分转换为离散的卷积运算。此外，由于B样条基函数在三个维度上是可分离的，因此作为张量积，可以将3D情况下的卷积作为简单的三个1D卷积序列执行。

(3)单元建模(即将设计空间分割为单元子域)和具有隐式三变量B样条基函数的水平集表示的结合为并行处理(尤其是在高性能上具有CPU或GPU多核的计算平台)中的高效计算开辟了前景。应该注意的是，HD-CLIBS模型完美而轻松地保持了微结构单元之间的几何连续性。这意味着HD-CLIBS模型在单元、晶格和功能梯度材料和结构的设计和优化方面具有巨大潜力。在这种情况下，可以以较小的规模选择或细化子域单元的大小，从而以可能是当前最新技术的数百倍高清晰度产生具有精细结构特征和细节的设计。这样的设计是对参数系数方面可能存在数百万或数十亿个设计变量(按参数系数计)的优化问题的解决方案。因此，在B样条曲线实现中的高清晰度单元水平集将结构拓扑优化扩展到超级计算机环境。

在本公开中，提出了均值一致性问题的拓扑优化，目的是说明该方法的实现方式和潜在的好处。基于水平集描述，可以给出结构均值一致性问题的的数学公式：

最小化_φ：J(u，φ)＝∫_Dε_ij(u)C_ijklε_kl(u)H(φ)dΩ

约束

其中J(u，φ)是目标函数，u是位移场，ε_ij是由二阶应变张量表示的应变场，C_ijkl是实体材料的四阶弹性张量，H是Heaviside函数，G是具有表示固体材料的最大允许体积V_max的体积约束方程。静态平衡方程以能量双线性形式a(u，v，φ)和载荷线性形式l(v，φ)的形式以其弱变分形式给出，其中v表示表示运动学上允许的位移场空间U中的虚拟位移场。

众所周知，在经典的基于水平集的拓扑优化方法中，拉格朗日函数很容易构造，其导数简单表示为

G表示由下式定义的梯度密度函数，

g(x)＝--ε_ij(u(x，t))C_ijklε_kl(u(x，t))+λ (25)

其中λ是拉格朗日乘数，用于处理方程(23)中的体积约束。

使用HD-CLIBS模型，拓扑优化将在子域单元D_s的水平进行。根据汉密尔顿-雅各比方程(方程(5))，将法向速度场

表示为：

其中s＝1，…，M (26)

从方程(8)，很容易获得

其中x∈Γ_s＝θΩ_s (27)

因此，将方程(27)代入(26)得到

其中s＝1，…，M (28)

将方程(28)中的法向速度代入式(24)，拉格朗日导数变为

与形状敏感性分析的一般情况一样，此导数定义为结构边界上的边界积分。通常，通过使用亥维赛函数(Heaviside function)H(φ)及其导数狄拉克德尔塔函数δ(φ)，将导数重新形式化为固定设计域D上的域积分。然后，将其转换为

对于最小化过程，此拉格朗日导数应保持为负。因此，可以使用最速下降法，得到所需导数系数

其中s＝1，…，M (31)

其中W(x)＝-g(x)δ(φ(x))。

在这里，使用了以下事实：全局结构域D被分割为M个子域单元D_s(s＝1，…，M)，使得

和

严格来说，在结构的边界

上定义等式(29)中的形状导数，这意味着只有每个包含边界的一部分的子域单元D_s的速度场

才有意义。因此，式(31)定义了H(φ)≠0的那些子域单元的B样条基函数系数的变化。对于其他单元，它们的B样条函数不会发生任何变化，并且它们的单元水平集函数对于改变或进化不是必需的。

因此，梯度下降法会导致B样条系数的变化，

其中Δt是迭代伪时间步长。在迭代过程中，B样条系数在每个子域单元上独立更新。更新后，可以通过将所有子域水平集函数φ_s拼凑在一起来获得全局域水平集函数φ。使用HD-CLIBS模型进行拓扑优化的计算复杂度非常低。因此，它可以直接应用于3D单元结构的大规模问题，克服了基于全局水平集公式的常规水平集方法中的许多困难。

D.数值实现

将进一步讨论在HD-CLIBS方法的计算实现中的许多数值特征。这些数值特征被充分利用以进行有效的实施。

1)用于梯度计算的离散可分卷积

在讨论中，假设子域单元是由分割全局设计域的直线单元组成的。对于每个子域单元，将B-样条基函数的次数在所有三个方向上均设置为p＝q＝r＝3，这意味着将部署三次B样条。而且，每个单元的边界面上的节点完全重复四次。在这方面，由边界节点所支持的基函数与由内部节点完全支持的基函数不同。在等间距为h的单元上的离散直线网格上，基于HD-CLIBS的优化的实际计算依赖于梯度密度函数g(x)或W(x)和B样条基函数的离散化。

对于完全由简单节点所支持的基函数，它们共享相同的形状轮廓，并且可以视为由具有单位空间的简单节点所支持的原始基函数β³(x))的放大和平移。因此，可以将这些函数表示为

其中对于三次样条，x_i是第i个节点的位置，h_x是节点区间或放大率，i∈[3，4，…，n-4，n-3]表示三次样条。因此，方程(31)的与这些基函数相关的B样条系数的导数在通常的3D中变为，

此外，特征函数W(x)可以沿着相同的序列离散化，通常会添加正则化过程，从而产生其离散版本W_h(x)。这种离散化方案最终获得以下形式的梯度密度和B样条曲线基函数的卷积，

其中

是β³的离散形式，

是基函数

的支持域。

由于B样条核是可分离的，因此可以将导数实质上计算为具有一维B样条核的梯度密度三个卷积的单个序列。这种卷积提供了一种计算系数导数的有效实现，从而通过方程(32)设置的单元集演变。

这种离散的卷积运算提供了额外的好处。使用B样条核，更新B样条系数的每个迭代步骤都可以看作是滤波操作。已知B样条核起低通滤波器的作用，并且在水平集演变和结构优化中引起固有的平滑。可以通过选择的B样条基函数的次数和节点间距来明确控制此平滑过程的效果。这些特性将成为为实际应用构建强大的单元结构设计和优化系统的有用工具。

2)重新规范化

三变量B样条表示法的另一个特征来自于水平集函数的边界属性。证明了B样条系数的l_∞范数提供了边界，

如所建议的，该边界属性可用于防止陡峭的梯度并避免重新初始化，这在经典的水平集实现中通常是必需的，并且已知是有问题的。在每个迭代步骤中，将以下重新规范化方案应用于获得的系数，

实际上，此操作将水平集函数的范围限制为[-1，1]。它还使水平集函数的梯度范数以某个值为界，从而防止了数值逼近中的不准确性。但是，此规范化操作不会阻止在优化过程中创建新的零级结构，从而保持了水平集表示的拓扑灵活性。

3)其他计算细节

在本文中，拉格朗日乘数λ是根据增强型拉格朗日乘数方案计算和更新的：

其中η是前N_R次迭代的参数集，而γ_i由下式更新，

γ_i+1＝min(γ_i+Δγ，γ_max) (38)

在优化过程中的进一步迭代中，其中Δγ和γ_max是参数γ的增量和上限。

在这里遵循现在的经典框架，针对亥维赛函数(尤其是狄拉克德尔塔函数)，将正则化近似定义为，

其中α是一个小的正数，而Δ描述了数值近似的宽度。

众所周知，当使用标准FEM时，结构状态方程的求解，特别是对于3D大规模或单元结构，在每次优化迭代中都非常耗时。在这项研究中，采用有限元方法(FCM)对结构进行线性弹性分析。在单元水平集方案中，全局设计域通常以直线元素被分割为一组子域单元。因此，由于我们的HD-CLIBS模型的下划线元素层次结构与FCM技术的要求兼容，因此实施有限元方法是非常自然的。在典型的FCM方法中，使用人造材料分布模型，其中空隙元素的有效弹性张量按比例缩小到很小的值，例如，在整个单元的封装中为10^-6。实际上，该策略已广泛应用于基于密度的拓扑优化方案中，例如SIMP方法。有限单元法允许采用更精确的方案，通过递归细分和基于单纯形的正交方案，使用多级方案评估修整后的单元的积分。尽管本文的重点不是结构分析的有限元建模，但在实现中结合了有限元方法的基本过程，以用于下面的数值示例。或者，可以考虑使用多节点扩展多尺度有限元方法(EMsFEM)来提高计算效率。图22示出了所提出的用于拓扑优化的HD-CLIBS方法的流程图。

E.数值示例

给出了三个3D示例，以证明该方法的构思和有效性。在第一示例中，将使用L＝4，H＝2和W＝2的直线设计域来优化固体悬臂梁，如图23所示。材料的杨氏模量为Es＝1泊松比为0.3。最大体积比设置为0.3。两种不同的加载条件被视为不同的优化案例(图23)。将设计域划分为4×2×2个单位单元，每个单元具有等间距的16×16×16网格。B样条的节点位于节点区间为0.0625的网格上。使用三次B样条。如上所述，在优化过程中，在连接单元之间应用了C⁰连续性条件。拓扑优化导致所谓的有机设计(图23的第69部分和第71部分)，这在以前的研究中是众所周知的。经过详细检查后，很容易确认结构的几何边界在任何相邻单元之间的公共边界上已完全连接。

接下来，如图24的第72部分中所示，在L＝16，H＝4和W＝4的设计域中考虑桥梁设计的一个相对复杂的情况。材料特性和最大体积比保持不变。桥梁在用L_s＝2的四个点固定在地面上。此外，在优化过程中，固定厚度为0.125的桥梁甲板，其中施加了一组均匀载荷P＝1。将设计域划分为16×4×4单位B样条单元，每个单元具有等间距的8×8×8网格。三次B样条的节点位于节点区间为0.125的网格上。由于模型的对称性，仅选择了设计域的1/4进行优化。在连接单元之间强制执行C⁰连续性约束。如图24的第73部分所示，通过拓扑优化获得了有机形状的桥。

在最后一个示例中，研究了具有周期性单元格和分层单元格的单元格结构的拓扑优化。如图25所示，设计域的尺寸为L＝8，W＝8，H＝4，最大体积比设置为0.3。整体结构在四个底角处固定在地面上，而集中载荷F＝4施加在顶面的中心。设计域被划分为8×8×4个单位单元，每个单元上具有等间距的16×16×16网格，用于三次B样条的节点分布，节点区间为0.0625。再次，应用C⁰连续性条件。由于结构对称，仅选择了设计域的1/4进行优化。

拓扑优化的初始结构设置如下。首先使用本文描述的模型初始化技术构建单位单元，使其具有0.53的相对体积分数(图25第75部分)。该晶胞在所有三个方向上都必须是周期性的，并且放置在HD-CLIBS模型的8×8×4子域单元中。换句话说，用于优化的初始结构是具有周期性重复的单元的周期性单元结构。然后，进行以下三种不同的优化。

一阶实体结构。作为常规的结构拓扑优化问题，将所有8×8×4个B样条单元的B样条系数视为设计变量，并且除连接单元之间的几何连续性必需的C⁰连续性条件外，不向其单元水平集函数施加任何其他限制。换句话说，B样条单元可以自由更改，只要它们保持C⁰连接即可。在优化过程的收敛过程中，发现如图26的第76部分所示，获得了一种具有常规有机结构优化已知特征的有机形状的单尺度实体结构。一些初始的单元变得空洞，而其他单元则在内部变成固体或被优化结构的边界表面切割。总材料体积比根据需要更改为0.30，其中材料主要在对角线方向上分布，以支撑来自顶部的载荷并将其传递到四个底部支撑。图26的第77部分示出了拓扑结构和有机形状的优化结构的四分之一。

周期性单元结构。如果要求B样条单元在所有三个方向上保持周期性，则优化过程将不允许消除任何单元，但是它将实质上优化单位单元以创建具有重复单元分布的优化单元结构。图27的第78部分示出了优化的单元结构，而图27的第79部分示出了优化的周期性单元。显然，优化单元中的材料会重新分布，主要沿垂直壁和对角壁在从负载的拐角到中心的方向分布。由于结构的四分之一模型用于优化，因此组装的整个单元结构由四个四分之一部分组成，对角线方向相应对齐，如图27的第78部分所示。

分层的单元结构。层状结构具有在一个或两个方向上周期性或重复的微结构。通过消除B样条单元在分层构造方向上的周期性连接性约束，可以将所提出的HD-CLIBS方法用于设计分层结构。因此，通过这样做，在该示例的图28的第80部分中获得了优化的分层单元结构。图28的第81部分示出了，由于将层中的单元约束为周期性，因此每个层中的单元均显示相同的几何构型。

Claims

1.一种用于三维单元结构建模和拓扑优化的B样条高清晰度单元水平集方法，包括：

将目标全局结构的全局设计域划分为多个子域单元；以及

针对所述多个子域单元中的每一个，使用隐式B样条函数定义用于表示子域单元的结构材料域的单元水平集函数，使得整组单元水平集函数包括整个全局设计域上的全局水平集函数。

2.根据权利要求1所述的方法，还包括：

定义相邻子域单元的隐式B样条函数的参数系数的兼容性约束，以保证它们的几何连通性达到指定阶的平滑度。

3.根据权利要求1所述的方法，还包括：

通过确定每个子域单元中的隐式B样条函数的系数来初始化所述单元水平集函数。

4.根据权利要求3所述的方法，其中，基于目标全局结构的体积数据和目标全局结构的实体或表面表示中的一种使用快速B样条差值法或快速离散B样条卷积运算获得隐式B样条函数的系数。

5.根据权利要求2所述的方法，还包括：将相邻子域单元的隐式B样条函数的参数系数的兼容性约束沿周期重复方向应用到子域单元的相对面，以表示周期单元。

6.根据权利要求1所述的方法，其中，通过在隐式B样条函数的节点序列中插入预定数量的特定节点来实现单元划分。

7.根据权利要求6所述的方法，其中，通过划分得到的子域单元中的隐式B样条函数的系数只能控制该子域单元中的B样条的形状，而对其他子域单元中的B样条的形状没有影响。

8.根据权利要求1所述的方法，其中，通过调整各个隐式B样条函数的参数系数来操纵和修改目标全局结构。

9.根据权利要求8所述的方法，其中，基于初始化的隐式B样条函数，通过同时在各个隐式B样条函数的系数上添加或减去正常数来更改全局特征厚度。

10.根据权利要求8所述的方法，其中，通过对隐式B样条函数的参数系数进行调整来将不同的子域单元融合在一起。

11.根据权利要求1所述的方法，其中，在各个子域单元上并行地执行隐式B样条函数的动态演化。

12.根据权利要求2所述的方法，其中，通过基于梯度的优化方法来更新隐式B样条函数的参数系数。

13.根据权利要求1所述的方法，其中，在子域单元的水平进行目标全局结构的拓扑优化。

14.一种计算机存储介质，其上存储有计算机程序，所述计算机程序由处理器执行时，使得处理器执行以下步骤：

将目标全局结构的全局设计域划分为多个子域单元；

15.根据权利要求14所述的计算机存储介质，其中，所述处理器还被配置为定义相邻子域单元的B样条函数的参数系数的兼容性约束，以保证它们的几何连通性达到指定阶的平滑度。