一种结构拓扑优化模型的灵敏度分析方法及应用

1.本发明属于结构优化设计领域，更具体地，涉及一种结构拓扑优化模型的灵敏度分析方法及应用。


背景技术：

2.传统的工程结构设计，往往需要工程设计人员具有多年的实践经验，这种设计模式已无法满足当今工业产品高效、高精、柔性化等诸多方面的要求。基于数值分析技术的结构优化技术可通过计算机快速生成合理有效的结构设计方案，是工程结构智能设计的主要实现方式之一。尺寸优化、形状优化以及拓扑优化是实现工程结构优化设计的三种主要途径，并且由于拓扑优化具有更广阔的设计空间、内部结构拓扑形式可变等优点，使得拓扑优化成为最具潜力的结构优化设计工具。然而，基于单元变量(或节点变量)描述形式的结构拓扑优化模型本质上是一个0-1整数规划问题，极易出现棋盘格、网格依赖等一系列的数值不稳定问题。
3.过滤器技术可有效消除优化结果存在的棋盘格和网格依赖等数值不稳定性，是固体各向同性惩罚模型、渐进结构、离散变量等一系列变密度拓扑优化方法进行优化设计变量更新的基础。因此，基于过滤器技术的变密度拓扑优化方法可有效实现工程结构的优化设计，并逐渐应用于航空航天、汽车工程、建筑、医疗器械等领域。而结构拓扑优化模型的灵敏度分析作为结构拓扑优化的重要环节，它反映的是约束函数和目标函数对优化设计变量的导数信息，表征了优化设计变量对输出响应的影响程度和影响规律，因此研究结构拓扑优化模型的灵敏度分析方法具有重要意义；
4.基于过滤器技术的灵敏度分析方法使用单元附近区域的灵敏度加权平均替代对应优化设计变量的初始灵敏度，进而实现结构拓扑优化模型的正则化处理。无论是基于单元中心点距离的显式过滤器还是基于b样条跨单元特性的隐式过滤器，均需要对所有单元的权重矩阵进行计算和存储，以便能快速进行变密度拓扑优化方法的灵敏度分析。然而，当单元尺寸和过滤器半径较大时，以索引全部单元的方式进行过滤器矩阵权重值的计算和存储，将导致过滤器权重矩阵的计算时间过长、存储占用空间过大以及优化变量更新效率慢等一系列问题，进而严重影响了结构拓扑优化模型的灵敏度分析以及变密度拓扑优化方法的应用场景。因此，现急需一种非全局索引的结构拓扑优化模型的灵敏度分析方法来提升过滤器权重矩阵和优化变量更新的效率。


技术实现要素：

5.针对现有技术的以上缺陷或改进需求，本发明提供了一种结构拓扑优化模型的灵敏度分析方法及应用，以解决现有灵敏度分析方法中结构拓扑优化过滤器权重矩阵计算时间长以及过滤器权重矩阵存储效率低的技术问题。
6.为了实现上述目的，第一方面，本发明提供了一种结构拓扑优化模型的灵敏度分析方法；其中，结构拓扑优化模型以工程结构分析网格中所有分析网格单元的应变能之和
为目标函数、所有分析网格单元的材料体积之和为体积约束函数构建得到，其优化设计变量为分析网格的b样条控制点坐标位置的材料密度；
7.上述灵敏度分析方法包括：
8.将工程结构各分析网格单元的应变能按照其在分析网格中的位置进行行排列整合，得到分析网格的应变能矩阵c；计算分析网格的应变能矩阵c关于单元中心点密度值的偏导，得到目标函数关于单元中心点密度值的偏导矩阵dc1；
9.将各分析网格单元的材料体积对其中心点密度的导数值按照其在分析网格中的位置进行排列整合，得到体积约束函数关于单元中心点密度值的偏导矩阵dv1；
10.采用基于张量积分解法转变后的目标函数关于优化设计变量的偏导矩阵表达式计算得到目标函数关于每个优化设计变量的偏导值；
11.采用基于张量积分解法转变后的体积约束函数关于优化设计变量的偏导矩阵表达式计算得到体积约束函数关于每个优化设计变量的偏导值；
12.其中，hm为分析网格在第m个参数方向上的子权重矩阵；m＝1,2,
…
,m；m为参数方向的数量，取值为2或3；hm的行数为第m个参数方向上b样条函数的个数，列数为第m个参数方向上的分析网格单元的个数；hm第r行第c列的元素hm(r,c)为第m个参数方向上的第r个b样条函数和第c个分析网格单元所对应的b样条函数值。
13.进一步优选地，上述基于张量积分解法转变后的目标函数关于优化设计变量的偏导矩阵表达式以及上述基于张量积分解法转变后的体积约束函数关于优化设计变量的偏导矩阵表达式的获取方法包括：
14.基于链式求导法，得到目标函数关于优化设计变量的偏导矩阵表达式h
×
dc1，以及体积约束函数关于优化设计变量的偏导矩阵表达式h
×
dv1；其中，h为单元中心点密度值关于优化设计变量的偏导矩阵，即分析网格的单元中心点b样条基函数值矩阵；
15.采用张量积分解法将偏导表达式h
×
dc1中的偏导矩阵h分解为分析网格不同参数方向上的子权重矩阵，从而将偏导表达式h
×
dc1转化为分析网格不同参数方向上的子权重矩阵与偏导矩阵dc1相乘的形式，即
16.采用张量积分解法将偏导表达式h
×
dv1中的偏导矩阵h同样分解为分析网格不同参数方向上的子权重矩阵，从而将偏导表达式h
×
dv1转化为分析网格不同参数方向上的子权重矩阵与偏导矩阵dv1相乘的形式，即
17.进一步优选地，上述子权重矩阵的获取方法包括：根据工程结构分析网格各参数方向上的b样条节点向量，得到每个参数方向上分析网格单元中心点的b样条函数值及b样条函数在该参数方向上的索引，从而得到每个参数方向上的子权重矩阵。
18.进一步优选地，分析网格的b样条控制点坐标位置基于工程结构的cad信息，结合b样条节点向量插入方法计算得到。
19.第二方面，本发明提供了一种工程结构的拓扑优化方法，包括以下步骤：
20.s11、根据工程结构的拓扑优化模型的实体材料体积约束，对工程结构的各分析网格单元的材料杨氏弹性模量进行初始化；
21.s12、基于各分析网格单元的材料杨氏弹性模量，得到工程结构的刚度矩阵；基于所得工程结构的刚度矩阵以及边界条件，求出每个分析网格单元的位移向量，进而求出各分析网格单元的应变能；
22.s13、以各分析网格单元的应变能之和为目标函数、各分析网格单元的材料体积之和为体积约束函数，构建结构拓扑优化模型；结构拓扑优化模型的优化设计变量为分析网格的b样条控制点坐标位置的材料密度；
23.s14、采用本发明第一方面所提供的灵敏度分析方法对上述结构拓扑优化模型进行灵敏度分析，得到目标函数和体积约束函数对每个优化设计变量的偏导值；
24.s15、将上述偏导值输入到优化求解算子得到各优化设计变量的更新值，进而对各分析网格单元的材料杨氏弹性模量进行更新；
25.s16、重复步骤s12-s15，直至相邻迭代步的目标函数相对变化值小于第一预设阈值或者分析网格单元的材料杨氏弹性模量所对应的相对材料密度的最大变化值小于第二预设阈值。
26.第三方面，本发明提供了一种结构拓扑优化模型的灵敏度分析方法，包括：其中，结构拓扑优化模型以工程结构分析网格中所有分析网格单元的应变能之和为目标函数、所有分析网格单元的材料体积之和为体积约束函数构建得到，其优化设计变量为分析网格单元中心点的材料密度；
27.上述灵敏度分析方法包括：
28.将工程结构各分析网格单元的应变能按照其在分析网格中的位置进行行排列整合，得到分析网格的应变能矩阵c；计算分析网格的应变能矩阵c关于优化设计变量的偏导，得到目标函数关于优化设计变量的偏导矩阵dc2；
29.将各分析网格单元的材料体积对其中心点密度的导数值按照其在分析网格中的位置进行排列整合，得到体积约束函数关于优化设计变量的偏导矩阵dv2，从而得到体积约束函数关于每个优化设计变量的偏导值；其中，
30.采用基于张量积分解法转变后的灵敏度修正公式对上述偏导矩阵dc2进行修正，得到修正后的目标函数关于每个优化设计变量的偏导值；
31.其中，fn为分析网格在第n个方向上的子权重矩阵；n＝1,2,
…
,n；n为分析网格几何维度方向的数量，取值为2或3；fn的行数和列数均为分析网格第n个方向上单元个数；fn第p行第q列的元素fn(p,q)为第n个方向上第p个单元和第q个单元之间的距离权重值。
32.进一步优选地，上述灵敏度修正公式的获取方法包括：
33.采用张量积分解法将灵敏度修正公式f
×
dc2中的过滤器权重矩阵f分解为分析网格不同方向上的子权重矩阵，从而将灵敏度修正公式f
×
dc2转化为分析网格不同方向上的子权重矩阵与偏导矩阵dc2相乘的形式，即
34.进一步优选地，过滤器权重矩阵f为单元中心点之间的权重矩阵，为所有分析网格在各方向上的子权重矩阵的kronecker乘积，即：
35.进一步优选地，上述子权重矩阵的获取方法包括：
36.根据工程结构的分析网格各方向上的网格单元的中心点之间的距离及最小过滤半径，计算每个方向上的子权重矩阵；
37.其中，r
min
为最小过滤半径，为第n个方向上第p个单元和第q个单元中心点在该方向上的距离。
38.第四方面，本发明提供了一种工程结构的拓扑优化方法，包括以下步骤：
39.s21、根据工程结构的拓扑优化模型的实体材料体积约束，对工程结构的各分析网格单元的材料杨氏弹性模量进行初始化；
40.s22、基于各分析网格单元的材料杨氏弹性模量，得到工程结构的刚度矩阵；基于所得工程结构的刚度矩阵以及边界条件，求出每个分析网格单元的位移向量，进而求出各分析网格单元的应变能；
41.s23、以各分析网格单元的应变能之和为目标函数、各分析网格单元的材料体积之和为体积约束函数，构建结构拓扑优化模型；结构拓扑优化模型的优化设计变量为分析网格单元中心点的材料密度；
42.s24、采用本发明第三方面所提供的灵敏度分析方法对上述结构拓扑优化模型进行灵敏度分析，得到目标函数和体积约束函数对每个优化设计变量的偏导值；
43.s25、将上述偏导值输入到优化求解算子得到各优化设计变量的更新值，进而对各分析网格单元的材料杨氏弹性模量进行更新；
44.s26、重复步骤s22-s25，直至网格单元的材料杨氏弹性模量所对应的相对材料密度的最大变化值小于第三预设阈值；对分析网格进行均匀细化；
45.s27、重复步骤s26，直至网格细化次数大于第四预设阈值。
46.第五方面，本发明提供了一种工程结构的拓扑优化方法系统，包括：存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时执行本发明第二方面或第四方面所提供的拓扑优化方法。
47.第六方面，本发明还提供了一种机器可读存储介质，所述机器可读存储介质存储有机器可执行指令，所述机器可执行指令在被处理器调用和执行时，所述机器可执行指令促使所述处理器实现本发明第一方面所提供的灵敏度分析方法、第二方面所提供的拓扑优化方法、第三方面所提供的灵敏度分析方法、以及第四方面所提供的拓扑优化方法中的一种或多种方法。
48.总体而言，通过本发明所构思的以上技术方案，能够取得以下有益效果：
49.1、本发明第一方面所提供的结构拓扑优化模型的灵敏度分析方法，基于张量积分解法转变后的目标函数关于优化设计变量的偏导矩阵表达式和采用基于张量积分解法转变后的体积约束函数关于优化设计变量的偏导矩阵表达式对结构拓扑优化模型进行灵敏度分析；本发明基于b样条基函数的张量分解特性，预先采用张量积分解法将偏导表达式h
×
dc1和偏导表达式h
×
dv1中偏导矩阵h分解为分析网格不同参数方向上的子权重矩阵，进而将偏导表达式转换成分析网格不同参数方向上的子权重矩阵与对应的偏导矩阵dc1和偏导矩阵dv1相乘的形式进行计算，故本发明仅需要存储每个方向上的b样条基函数的值，而不需要存储所有单元的b样条基函数的值，即可实现b样条函数值的等效表达，本发明无需对作用在每个单元的函数值及其索引进行存储，仅需要存储每个方向上的所有单元的函数值及其索引，存储效率较高，且计算速度也较快。
50.2、在本发明第二方面所提供的工程结构的拓扑优化方法的迭代更新过程中，采用本发明第一方面所提供的方法进行灵敏度分析，每一步灵敏度分析的计算效率和存储效率存在显著提高；在工程结构分析域的离散等几何分析网格的各方向上b样条单元个数均相等且为n时，工程结构的等几何分析权重矩阵的空间复杂度可以从o(n3)减少到o(n)，得到了较为显著的优化。
51.3、本发明第三方面所提供的结构拓扑优化模型的灵敏度分析方法，在获取目标函数和体积约束函数关于优化设计变量的偏导矩阵dc2和偏导矩阵dv2后，对目标函数关于优化设计变量的偏导矩阵dc2进行一步进行修正，本发明基于张量积的分解特性，预先采用张量积分解法将灵敏度修正公式f
×
dc2中的显式过滤器权重矩阵f分解为分析网格不同方向上的子权重矩阵，从而将灵敏度修正公式f
×
dc2转化为分析网格不同方向上的子权重矩阵与偏导矩阵dc2相乘的形式进行修正计算，故本发明仅需要存储分析网格不同方向上的子权重矩阵，而不需要存储显式过滤器权重矩阵f中的所有权重值，存储效率较高，且计算速度也较快。
52.4、在本发明第四方面所提供的工程结构的拓扑优化方法的迭代更新过程中，采用本发明第三方面所提供的方法进行灵敏度分析，每一步灵敏度分析的计算效率和存储效率存在显著提高。
附图说明
53.图1为本发明实施例1所提供的2
×
2大小的b样条等几何分析网格示意图；
54.图2为本发明实施例1提供的b样条等几何分析网格的ien列阵示意图；其中，(a)为x方向上b样条等几何分析网格的ien列阵示意图，(b)为y方向上b样条等几何分析网格的ien列阵示意图；
55.图3为本发明实施例2所提供的工程结构的拓扑优化方法流程图；
56.图4为本发明实施例3所提供的过滤器的结构示意图，其中，(a)为过滤半径为2的基于单元中心点欧式距离的传统过滤器，(b)为过滤半径为2的基于中心点轴向距离张量积结构的显式过滤器；
57.图5为本发明实施例4所提供的工程结构的拓扑优化方法流程图；
58.图6为本发明实施例7提供的基于本发明所提供的基于张量积分解的隐式过滤器
的灵敏度分析方法和基于传统基于隐式过滤器的灵敏度分析方法进行等几何结构拓扑优化时计算权重矩阵所需时间的比较结果示意图；
59.图7为本发明实施例7提供的基于本发明所提供的基于张量积分解的隐式过滤器的灵敏度分析方法和基于传统基于隐式过滤器的灵敏度分析方法进行等几何结构拓扑优化时权重矩阵所占内存的比较结果示意图；
60.图8为本发明实施例8提供的基于本发明所提供的基于张量积分解的显式过滤器的灵敏度分析方法和基于传统基于显式过滤器的灵敏度分析方法进行多分辨率结构拓扑优化时计算权重矩阵所需时间的比较结果示意图；
61.图9为本发明实施例8提供的基于本发明所提供的基于张量积分解的显式过滤器的灵敏度分析方法和基于传统基于显式过滤器的灵敏度分析方法进行多分辨率结构拓扑优化时权重矩阵所占内存的比较结果示意图。
具体实施方式
62.为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
63.实施例1、
64.本实施例提供了一种结构拓扑优化模型的灵敏度分析方法，该方法是基于隐式过滤器技术的灵敏度分析方法；其中，结构拓扑优化模型以工程结构分析网格中所有分析网格单元的应变能之和为目标函数、所有分析网格单元的材料体积之和为体积约束函数构建得到，其优化设计变量为分析网格的b样条控制点坐标位置的材料密度；其中，分析网格的b样条控制点坐标位置基于工程结构的cad信息，结合b样条节点向量插入方法计算得到。
65.具体地，本实施例所提供的灵敏度分析方法包括：
66.将工程结构各分析网格单元的应变能按照其在分析网格中的位置进行行排列整合，得到分析网格的应变能矩阵c；计算分析网格的应变能矩阵c关于单元中心点密度值的偏导，得到目标函数关于单元中心点密度值的偏导矩阵dc1；
67.将各分析网格单元的材料体积对其中心点密度的导数值按照其在分析网格中的位置进行排列整合，得到体积约束函数关于单元中心点密度值的偏导矩阵dv1；
68.采用基于张量积分解法转变后的目标函数关于优化设计变量的偏导矩阵表达式计算得到目标函数关于每个优化设计变量的偏导值；
69.采用基于张量积分解法转变后的体积约束函数关于优化设计变量的偏导矩阵表达式计算得到体积约束函数关于每个优化设计变量的偏导值；
70.其中，hm为分析网格在第m个参数方向上的子权重矩阵；m＝1,2,
…
,m；m为参数方向的数量，取值为2或3；hm的行数为第m个参数方向上b样条函数的个数，列数为第m个参数
方向上的分析网格单元的个数；hm第r行第c列的元素hm(r,c)为第m个参数方向上的第r个b样条函数和第c个分析网格单元所对应的b样条函数值。
71.具体地，上述子权重矩阵的获取方法包括：根据工程结构分析网格各参数方向上的b样条节点向量，得到每个参数方向上分析网格单元中心点的b样条函数值及b样条函数在该参数方向上的索引，从而得到每个参数方向上的子权重矩阵。
72.为了进一步说明本实施例所提供的灵敏度分析方法，下面以2
×
2大小的b样条等几何分析网格为例进行详述：
73.如图1所示为2
×
2大小的b样条等几何分析网格示意图，其中，b样条等几何分析网格基于工程结构的cad模型数据信息生成；通过对工程结构的b样条等几何分析网格的基函数的局部索引和b样条等几何分析网格中各等几何分析单元的单元编号进行组合排序，得到如图2所示的各个参数方向(这里为几何参数，即x方向和y方向)的b样条等几何网格的ien列阵；其中，(a)为x方向上b样条等几何分析网格的ien列阵示意图；(b)为y方向上b样条等几何分析网格的ien列阵示意图；x和y方向等几何分析单元的个数都为2个，每一个b样条单元在每一参数方向上均对应3个b样条基函数，即对应3个b样条基函数的局部索引。最后，基于每个方向的ien列阵，通过稀疏矩阵组装法对工程结构的等几何单元数据进行组装；在本例中，最终组装得到的子权重矩阵的行索引对应相关b样条函数的位置，列索引对应单元的位置。最终得到的结果如下：
[0074][0075]
其中，h1为分析网格在x方向上的子权重矩阵，h2为分析网格在y方向上的子权重矩阵。
[0076]
在获取了子权重矩阵后，进一步获取目标函数关于单元中心点密度值的偏导矩阵进而采用基于张量积分解法转变后的目标函数关于优化设计变量的偏导矩阵表达式，计算得到目标函数关于优化设计变量的偏导矩阵为：
[0077][0078]
需要说明的是，上述基于张量积分解法转变后的目标函数关于优化设计变量的偏导矩阵表达式以及上述基于张量积分解法转变后的体积约束函数关于优化设计变量的偏导矩阵表达式的获取方法包括：
[0079]
基于链式求导法，得到目标函数关于优化设计变量的偏导矩阵表达式h
×
dc1，以及体积约束函数关于优化设计变量的偏导矩阵表达式h
×
dv1；其中，h为单元中心点密度值关于优化设计变量的偏导矩阵，即分析网格的单元中心点b样条基函数值矩阵；
[0080]
采用张量积分解法将偏导表达式h
×
dc1中的偏导矩阵h分解为分析网格不同参数方向上的子权重矩阵，从而将偏导表达式h
×
dc1转化为分析网格不同参数方向上的子权重
矩阵与偏导矩阵dc1相乘的形式，即
[0081]
采用张量积分解法将偏导表达式h
×
dv1中的偏导矩阵h同样分解为分析网格不同参数方向上的子权重矩阵，从而将偏导表达式h
×
dv1转化为分析网格不同参数方向上的子权重矩阵与偏导矩阵dv1相乘的形式，即
[0082]
本实施例基于b样条基函数的张量分解特性，即b样条基函数值矩阵h中存在冗余，已知其各个参数方向上的b样条基函数的值，即可还原出矩阵h，因此，仅需要存储每个方向上的b样条基函数的值，而不需要存储空间中所有单元的b样条基函数的值。可以实现b样条函数值的等效表达，本发明无需对作用在每个单元的函数值及其索引进行存储，仅需要存储每个方向上的所有单元的函数值及其索引，存储效率较高，且计算速度也较快。
[0083]
实施例2、
[0084]
本发明提供了一种工程结构的拓扑优化方法，该方法为基于可分解隐式过滤器的结构拓扑优化方法，如图3所示，包括以下步骤：
[0085]
s11、根据工程结构的拓扑优化模型的实体材料体积约束，对工程结构的各分析网格单元的材料杨氏弹性模量进行初始化；
[0086]
s12、基于各分析网格单元的材料杨氏弹性模量，得到工程结构的刚度矩阵；基于所得工程结构的刚度矩阵以及边界条件，求出每个分析网格单元的位移向量，进而求出各分析网格单元的应变能；
[0087]
s13、以各分析网格单元的应变能之和为目标函数、各分析网格单元的材料体积之和为体积约束函数，构建结构拓扑优化模型；结构拓扑优化模型的优化设计变量为分析网格的b样条控制点坐标位置的材料密度；
[0088]
s14、采用本发明实施例1所提供的灵敏度分析方法对上述结构拓扑优化模型进行灵敏度分析，得到目标函数和体积约束函数对每个优化设计变量的偏导值；
[0089]
s15、将上述偏导值输入到优化求解算子得到各优化设计变量的更新值，进而对各分析网格单元的材料杨氏弹性模量进行更新；其中，所采用的优化求解算子可以为优化准则法、移动渐进线等；
[0090]
s16、重复步骤s12-s15，直至相邻迭代步的目标函数相对变化值小于第一预设阈值(本实施例中取值为10-6
)或者分析网格单元的材料杨氏弹性模量所对应的相对材料密度的最大变化值小于第二预设阈值(本实施例中取值为10-2
)。
[0091]
相关技术方案同实施例1、这里不做赘述。
[0092]
实施例3、
[0093]
本实施例提供了一种结构拓扑优化模型的灵敏度分析方法，该方法是基于显式过滤器技术的灵敏度分析方法；包括：其中，结构拓扑优化模型以工程结构分析网格中所有分析网格单元的应变能之和为目标函数、所有分析网格单元的材料体积之和为体积约束函数构建得到，其优化设计变量为分析网格单元中心点的材料密度；
[0094]
上述灵敏度分析方法包括：
[0095]
将工程结构各分析网格单元的应变能按照其在分析网格中的位置进行行排列整
合，得到分析网格的应变能矩阵c；计算分析网格的应变能矩阵c关于优化设计变量的偏导，得到目标函数关于优化设计变量的偏导矩阵dc2；
[0096]
将各分析网格单元的材料体积对其中心点密度的导数值按照其在分析网格中的位置进行排列整合，得到体积约束函数关于优化设计变量的偏导矩阵dv2，从而得到体积约束函数关于每个优化设计变量的偏导值；其中，体积约束函数关于每个优化设计变量的偏导值均为1；
[0097]
采用基于张量积分解法转变后的灵敏度修正公式对上述偏导矩阵dc2进行修正，得到修正后的目标函数关于每个优化设计变量的偏导值；并将目标函数关于每个优化设计变量的偏导值更新为修正后的目标函数关于每个优化设计变量的偏导值；
[0098]
其中，fn为分析网格在第n个方向上的子权重矩阵；n＝1,2,
…
,n；n为分析网格几何维度方向的数量，取值为2或3；fn的行数和列数均为分析网格第n个方向上单元个数；fn第p行第q列的元素fn(p,q)为第n个方向上第p个单元和第q个单元之间的距离权重值。
[0099]
具体地，上述子权重矩阵的获取方法包括：
[0100]
根据工程结构的分析网格各方向上的网格单元的中心点之间的距离及最小过滤半径，计算每个方向上的子权重矩阵；
[0101]
其中，r
min
为最小过滤半径，为第n个方向上第p个单元和第q个单元中心点在该方向上的距离。
[0102]
以n等于2为例，记第p个单元的中心点位置坐标为(i
p
,j
p
)，第q个单元的中心点位置坐标为(iq,jq)，则为x方向上第p个单元和第q个单元中心点在该方向上的距离，即|i
p-iq|，为y方向上第p个单元和第q个单元中心点在该方向上的距离，即|j
p-jq|。n等于3时与n等于2时情况类似，这里不做赘述。
[0103]
需要说明的是，上述灵敏度修正公式的获取方法包括：
[0104]
采用张量积分解法将灵敏度修正公式f
×
dc2中的过滤器权重矩阵f分解为分析网格不同方向上的子权重矩阵，从而将灵敏度修正公式f
×
dc2转化为分析网格不同方向上的子权重矩阵与偏导矩阵dc2相乘的形式，即
[0105]
进一步地，过滤器权重矩阵f为单元中心点之间的权重矩阵，为所有分析网格在各方向上的子权重矩阵的kronecker乘积，即：
[0106]
以二维过滤器为例，此时n为2，如图4所示为过滤器的结构示意图，其中，(a)为过
滤半径为2的基于单元中心点欧式距离的传统过滤器；(b)为过滤半径为2的基于中心点轴向距离张量积结构的显式过滤器，该过滤器的构造方式如下：
[0107][0108][0109][0110]
其中，表示第p个设计变量(即第p个单元)对第q个设计变量(即第q个单元)的权重影响值，可以由两个方向上的子权重值构成，因此可以构成基于张量积结构的显式过滤器。
[0111]
本实施例基于张量积的分解特性，显式过滤器权重矩阵可以分解为沿着每个方向上的子权重矩阵。相比于传统的显式过滤矩阵，本发明的存储效率和计算速度都有着显著的提升。
[0112]
实施例4、
[0113]
一种工程结构的拓扑优化方法，该方法为基于可分解显式过滤器的结构拓扑优化方法，也是一种多分辨率结构拓扑优化方法，如图5所示，包括以下步骤：
[0114]
s21、根据工程结构的拓扑优化模型的实体材料体积约束，对工程结构的各分析网格单元的材料杨氏弹性模量进行初始化；
[0115]
s22、基于各分析网格单元的材料杨氏弹性模量，得到工程结构的刚度矩阵；基于所得工程结构的刚度矩阵以及边界条件，求出每个分析网格单元的位移向量，进而求出各分析网格单元的应变能；
[0116]
s23、以各分析网格单元的应变能之和为目标函数、各分析网格单元的材料体积之和为体积约束函数，构建结构拓扑优化模型；结构拓扑优化模型的优化设计变量为分析网格单元中心点的材料密度；
[0117]
s24、采用本发明实施例3所提供的灵敏度分析方法对上述结构拓扑优化模型进行灵敏度分析，得到目标函数和体积约束函数对每个优化设计变量的偏导值；
[0118]
s25、将上述偏导值输入到优化求解算子得到各优化设计变量的更新值，进而对各分析网格单元的材料杨氏弹性模量进行更新；
[0119]
s26、重复步骤s22-s25，直至网格单元的材料杨氏弹性模量所对应的相对材料密度的最大变化值小于第三预设阈值(本实施例中取值为10-1
)；对分析网格进行均匀细化；
[0120]
需要说明的是，在对分析网格进行均匀细化后，下一次迭代过程中进行灵敏度分析时所采用的子权重矩阵也会发生变化，子权重矩阵会根据细化网格的中心点距离以及各维度方向的单元数目进行更新；优选地，对分析网格进行均匀细化后，可以结合细化程度将下一次灵敏度分析中的最小过滤半径r
min
更新为其当前值的k倍；优选地，k∈[1.5,2]。
[0121]
s27、重复步骤s26，直至网格细化次数大于第四预设阈值(本实施例中取值为3)。
[0122]
相关技术方案同实施例3、这里不做赘述。
[0123]
实施例5、
[0124]
一种工程结构的拓扑优化方法系统，包括：存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时执行本发明实施例2或实施例4所提供的拓扑优化方法。
[0125]
相关技术方案同实施例2或实施例4，这里不做赘述。
[0126]
实施例6、
[0127]
一种机器可读存储介质，所述机器可读存储介质存储有机器可执行指令，所述机器可执行指令在被处理器调用和执行时，所述机器可执行指令促使所述处理器实现本发明实施例1所提供的灵敏度分析方法、实施例2所提供的拓扑优化方法、实施例3所提供的灵敏度分析方法、以及实施例4所提供的拓扑优化方法中的一种或多种方法。
[0128]
相关技术方案同实施例1-实施例4，这里不做赘述。
[0129]
实施例7、
[0130]
为了进一步说明本发明实施例1所提供的基于可分解隐式过滤器的灵敏度分析方法的有效性，进行数值实验验证如下：
[0131]
在硬件条件为xeon(r)gold 5120cpu@2.20ghz 2.19ghz intel core and 64gb ram且软件环境由windows 10操作系统和matlab r2021a构成的计算机中，对本发明实施例1所提供的基于张量积分解的隐式过滤器的灵敏度分析方法的有效性进行验证。分别采用本发明所提供的灵敏度分析方法和传统基于隐式过滤器的灵敏度分析方法在等几何分析中来进行比较，其中用于表达设计域材料密度分布的b样条函数阶次均为10。实验结果如图6和图7所示：其中，图6为基于本发明所提供的基于张量积分解的隐式过滤器的灵敏度分析方法和基于传统基于隐式过滤器的灵敏度分析方法进行等几何结构拓扑优化时计算权重矩阵所需时间的比较结果示意图(横坐标为等几何分析网格的分辨率，纵坐标为进行灵敏度分析时计算权重矩阵所需要的时间)；图7为基于本发明所提供的基于张量积分解的隐式过滤器的灵敏度分析方法和基于传统基于隐式过滤器的灵敏度分析方法进行等几何结构拓扑优化时权重矩阵所占内存的比较结果示意图(横坐标为等几何分析网格的分辨率，纵坐标为进行灵敏度分析时权重矩阵所占的内存)；从图6可以看出，当将本发明所提供的可分解隐式过滤器替换传统隐式过滤器时，传统方法计算权重矩阵所需的时间与本发明实施例1计算权重矩阵所需的时间的比值最大可达到5300。于此同时，从图7中可以看出，传统隐式过滤器与本发明所提供的可分解隐式过滤器的权重矩阵所占内存之比的最大值为4976。由此可见本发明所提供的基于可分解隐式过滤器的灵敏度分析方法的内存效率和预处理效率都得到了显著提升。
[0132]
实施例8、
[0133]
为了进一步说明本发明实施例3所提供的基于可分解显式过滤器的灵敏度分析方法的有效性。在与实施例7相同的硬件和软件环境下，进行如下数值实验：
[0134]
分别采用本发明所提供的基于张量积分解的显式过滤器和传统显式过滤器来进行比较，其中每种过滤器中的初始过滤半径值都是8。实验结果如图8和图9所示：其中，图8为基于本发明所提供的基于张量积分解的显式过滤器的灵敏度分析方法和基于传统基于显式过滤器的灵敏度分析方法进行多分辨率结构拓扑优化时计算权重矩阵所需时间的比
较结果示意图(横坐标为网格的分辨率，纵坐标为进行灵敏度分析时计算权重矩阵所需要的时间)；图9为基于本发明所提供的基于张量积分解的显式过滤器的灵敏度分析方法和基于传统基于显式过滤器的灵敏度分析方法进行多分辨率结构拓扑优化时的权重矩阵所占内存的比较结果示意图(横坐标为网格的分辨率，纵坐标为进行灵敏度分析时权重矩阵所占的内存)；从图8可以看出，当将本发明所提供的基于张量积分解的显式过滤器替换传统显式过滤器时，传统方法计算权重矩阵所需的时间与本发明实施例3计算权重矩阵所需的时间之比的最大值为2217.1。于此同时，从图9中可以看出，传统显式过滤器与本发明所提供的可分解显式过滤器的权重矩阵所占内存之比的最大值为727.5。由此可见本发明所提供的基于张量积分解的显式过滤器的预处理效率和内存效率都有着非常显著的提升。
[0135]
综上所述，本发明实施例1基于b样条函数张量积分解的特性，将传统的多维度隐式过滤器权重矩阵分解为几个单一维度的子权重矩阵，从而大大减少权重矩阵的存储空间和计算时间。本发明实施例3基于工程结构所划分网格的每个方向上的网格数量以及最小过滤半径，得到每个参数单元方向上的子权重矩阵。各个方向上的子权重矩阵通过kronecker矩阵乘积法可以得到整体的矩阵，从而实现权重矩阵的等效表达。基于上述方法，本发明无需对每一个单元的敏度过滤权重矩阵进行计算，仅需要计算各个维度方向上的子权重矩阵即可，存储效率高，且计算效率较快。
[0136]
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。