一种基于多微结构原型和多尺度有限元的多孔结构设计方法

1.本发明属于工程结构优化领域，具体涉及一种基于多微结构原型和多尺度有限元的多孔结构设计方法。


背景技术：

2.多孔结构是一种多尺度结构，具备多尺度特征，即最大特征尺寸与最小特征尺寸相差甚大。随制造技术高速发展，多孔材料(cellular structure)凭借轻量化、比强高、隔热降噪、减震吸能等优点被广泛应用于各行各业。
3.以具有两个几何特征尺度的多孔结构为例，将其两个尺度分别称为宏观尺度和微观尺度。这两个尺度是相互影响的。在微观尺度上，微观结构构型显著影响其宏观等效物理性质。在宏观尺度上，结构的边界条件和外部荷载分布会影响结构的局部应力应变状态。而尺度关联的设计方法通常会消耗大量计算资源。为减少多尺度结构计算量，均匀化方法、代表体元法、多尺度有限元法等相继提出。
4.目前，已经有许多方法来解决这个问题。代表性的是基于计算均匀化方法的多尺度结构优化设计算法。然而，由于均质化方法建立在尺度分离的假设基础上，因此不能保证优化后的多尺度结构的可制造性和相邻微结构之间的连续性。为了保证相邻微结构之间的平滑连接，人们做了很多工作，包括采用去均质技术、映射和插值技术、全尺寸分析技术、数据驱动技术。基于可变切割(vcut)水平集方法，liu等人和xia等人开发了一种多变量切割(m-vcut)水平集方法，用于优化相邻微观结构之间至少具有c0连续性的细胞结构。在该方法中，微结构构型的优化设计转化为相应粗尺度单元节点切割高度的变化，在保持微结构性能的同时，大大减少了设计变量的数量。而m-vcut水平集方法是基于全尺度网格进行结构分析和灵敏度分析，会消耗大量的计算资源和计算时间。这严重限制了该方法在实际工程中的应用，特别是对三维问题的应用。


技术实现要素：

5.本发明主要目标是提出一个针对多孔材料的拓扑优化框架。使用水平集函数(level set method，lsm)来描述微结构边界及拓扑形状，并选取了一种多尺度计算手段即扩展多尺度有限元法(emsfem)来计算微结构性能，构建离线数据库，通过在线的拓扑优化分析设计得到宏观最优拓扑结构，最后重构细尺度模型。该框架被称为基于数据驱动的多孔材料结构拓扑优化方法，可以节省大量计算资源和时间。
6.本发明采用如下技术方案：
7.一种基于多微结构原型和多尺度有限元的多孔结构设计方法，包括以下步骤：
8.步骤s1.定义多孔结构微结构原型的水平集函数，并用多种切割平面切割这些水平集函数获得一系列具有不同相对密度的基本微结构；
9.步骤s2.将这些基本微结构进行相互组合生成一系列真实微结构，计算所述真实微结构的等效单元刚度矩阵和相对密度；
10.步骤s3.建立所述真实微结构等效刚度矩阵、微结构相对密度、设计变量之间的映射关系，形成离线数据库；
11.步骤s4.基于粗尺度网格进行优化，基于离线数据库中存储的数据，采用插值和数值求导技术进行结构分析和灵敏度分析，粗尺度优化后，得到设计变量的分布；
12.步骤s5.根据步骤s4中基本微结构设计变量分布和步骤s1的水平集函数，对原尺寸多孔结构进行几何表示和重构，最终得到全尺寸多孔结构。
13.进一步地，步骤s1中基本微结构用如下水平集函数表示：
[0014][0015]
其中，φm为微结构原型，hm是设计变量，即2d问题中的切割曲面的高度及3d问题中切割曲面的等值线。
[0016]
进一步地，所述步骤s2中用扩展多尺度有限元法计算真实微结构等效刚度矩阵，包括如下子步骤：
[0017]
1)数值基函数的构造，采用周期边界条件来构造数值基函数；
[0018]
2)建立了真实微结构的数值基函数后，通过下式计算真实微结构的等效单元刚度矩阵：
[0019]
kc＝n
tksn[0020]
其中，n为数值基函数，ks为子网格上的总体刚度矩阵，kc为真实微结构的等效单元刚度矩阵。
[0021]
进一步地，所述步骤s2中用利用数值积分法计算微结构的相对密度，具体为：微结构采用规则的固定背景网格进行网格划分，其中包括三种单元，即空单元、实体单元和混合单元，采用数值积分法计算微结构的相对密度，对于空单元和实体单元，使用2
dim
个积分点，对于混合单元，积分点数设置为6
dim
，其中dim表示微结构的维度。
[0022]
进一步地，所述步骤s3中建立所述真实微结构等效刚度矩阵、微结构相对密度、设计变量之间的映射关系为kc(h1,h2,
…
,hm)和ρm(h1,h2,
…
,hm)。
[0023]
进一步地，所述步骤s4中，进行粗尺度优化，以材料用量约束下的最小柔顺度优化为例，用如下优化列式进行粗尺度优化:
[0024][0025]
这里j为目标函数即柔度；v表示材料使用量；表示材料的使用上限；上标
[0026]
‘
min’和
‘
max’表示设计变量的下限和上限；为设计变量,表示第e个粗单元中第m个水平集函数的切割曲面的高度或切割曲面的等值线；fc和uc分别为结构在粗网格上的外力向量和位移向量.uc可以通过在粗尺度网格上求解平衡方程kcuc＝fc确定，p为水平集函数个数。
[0027]
进一步地，所述步骤s4中，灵敏度分析具体如下：
[0028]
由于外力向量与设计变量无关，因此目标函数相对于设计变量的导数为
[0029][0030]
由于设计变量与外力向量无关，下标
‘
m’代表相应的基本微结构；kc为总体刚度矩阵；可以由如下公式计算:
[0031][0032]cc
,ce为扩充矩阵，从已建立的离线数据库中通过数值求导或插值的手段获得，然后，将公式带回上式，得：
[0033][0034]
这里ue和ke分别为第
‘
e’个粗单元的节点位移向量和刚度矩阵；
[0035]
材料的使用量v可写为：
[0036][0037]
其中，vc和分别为粗单元的体积和相应相应微结构的相对密度，则对于第
‘
e’个粗单元，v关于设计变量的导数为：
[0038][0039]
可以从已建立的离线数据库中真实微结构的ρm和切割曲面之间的关系图通过数值求导及插值技术来获得。
[0040]
进一步地，所述步骤s5中，为了保证界面处相邻微结构的平滑连续性，首先将每个粗元的设计变量通过加权平均策略转换到节点上，然后，利用节点处的设计变量和双线性形状函数插值技术，可形成多个变高度切割曲面，将生成的变高度切割曲面切割对应微结构原型的水平集函数，得到由对应基本微结构组成的梯度多孔结构。
[0041]
与现有技术相比，本发明的有益效果为：
[0042]
1)数据驱动思想。所有离线数据都是在优化前计算的。离线数据库只需要构建一次即可进行多种模型或工况的优化问题。
[0043]
2)计算效率高。与传统的m-vcut水平集方法相比，该方法在设计变量较少的情况下，可获得略低的目标函数值。此外，对于二维悬臂梁问题，该方法的计算效率比传统的m-vcut水平集方法提高了近141倍。可以想象，在三维问题的计算效率将大大提高。
[0044]
3)该框架可进行多微结构原型和多材料的复合多孔结构设计，不止限于本发明所提到的微结构，基本微结构种类多样，可自由选取组合。
附图说明
[0045]
图1为本发明实例中的水平集法原理示意图，用一个三维水平集函数表示一个二维微结构ω。
[0046]
图2为本发明实例中使用多种微结构原型形成的真实微结构示例图；其中，(a)由
方形和
‘
x’形微结构原型生成；(b)由
‘
＝'和
‘
//'形状的微结构原型生成；(c)由
‘
||’，
‘
//’和
‘
\\'形状的微结构原型生成；(d)由3d
‘
cubic’和3d
‘
x’形微结构原型生成。
[0047]
图3为本发明实例中在emsfem中构造数值基函数所用到的位移边界条件图。
[0048]
图4为本发明实例中子网格中三种单元的示意图。
[0049]
图5为本发明实例中2d真实微结构的kc[1,1]和设计变量之间的关系图；其中，(a)仅用
‘
box’形基本微结构生成的真实微结构；(b)仅用
‘
x’形基本微结构生成的真实微结构；(c)采用
‘
box’和
‘
x’形基本微结构生成的真实微结构。
[0050]
图6为本发明实例中用
‘
＝’，
‘
||’，
‘
//'，
‘
\\’四种基本微结构生成的2d真实微结构的kc[1,1]和设计变量之间的关系图；其中，(a)～(f)分别表示hr＝0,0.2,0.4,0.6,0.8,1时的关系图。
[0051]
图7为本发明实例中3d真实微结构的kc[1,1]和设计变量之间的关系图；其中，(a)仅用3d
‘
cubic’形基本微结构生成的真实微结构；(b)仅用3d
‘
x’形基本微结构生成的真实微结构；(c)采用3d
‘
cubic’和3d
‘
x’形基本微结构生成的真实微结构。
[0052]
图8为本发明实例中2d真实微结构的相对密度ρm和设计变量之间的关系；其中，(a)仅用
‘
box’形基本微结构生成的真实微结构；(b)仅用
‘
x’形基本微结构生成的真实微结构；(c)采用
‘
box’和
‘
x’形基本微结构生成的真实微结构。
[0053]
图9为本发明实例中用
‘
＝’，
‘
||’，
‘
//'，
‘
\\’四种基本微结构生成的2d真实微结构的ρm和设计变量之间的关系图；其中，(a)～(f)分别表示hr＝0,0.2,0.4,0.6,0.8,1时的关系图。
[0054]
图10为本发明实例中3d真实微结构的ρm和设计变量之间的关系图；其中，(a)仅用3d
‘
cubic’形基本微结构生成的真实微结构；(b)仅用3d
‘
x’形基本微结构生成的真实微结构；(c)采用3d
‘
cubic’和3d
‘
x’形基本微结构生成的真实微结构。
[0055]
图11为本发明实例中重构过程示意图。
[0056]
图12为本发明实例中基于emsfem与mvcut水平集法的数据驱动结构优化的计算流程图。
[0057]
图13为本发明实例1中的悬臂梁模型图；其中，(a)展示了边界条件及粗网格划分(包含40
×
20个粗单元)；(b)为所使用的四种基本微结构。
[0058]
图14为本发明实例1中使用本发明提出的方法优化的结果图。其中，(a)为hh的分布图；(b)为基于hh的分布重构出的
‘
＝’形基本微结构的结构图；(c)为hv的分布图；(d)为基于hv的分布重构出的
‘
||’形基本微结构的结构图；(e)为h
l
的分布图；(f)为基于h
l
的分布重构出的
‘
//’形基本微结构的结构图；(g)为hr的分布图；(h)为基于hr的分布重构出的
‘
\\’形基本微结构的结构图。
[0059]
图15为本发明实例1中优化出的多孔结构图；其中，(a)为使用m-vcut法在全尺度网格上的优化结果；(b)为使用本发明提出的方法在粗尺度网格上优化得到的结果。
[0060]
图16为本发明实例1中目标函数和约束函数的变化历史图；其中，(a)为使用m-vcut法在全尺度网格上的优化结果；(b)为使用本发明提出的方法在粗尺度网格上优化得到的结果。
[0061]
图17为本发明实例2中的mbb梁及其边界条件示意图。
[0062]
图18为本发明实例2中用提出的方法对优化列式的优化结果。其中，(a)为hb的分布图；(b)为
‘
box’形基本微结构的全尺度多孔结构图；
[0063]
图19为本发明实例2中用提出的方法对优化列式的优化结果。其中，(a)为h
x
的分布图；(b)为
‘
x’形基本微结构的全尺度多孔结构图；
[0064]
图20为本发明实例2中用提出的方法对优化列式的优化结果。(a)为hb的分布图；(b)为h
x
的分布图；(c)为
‘
x’、
‘
box’形基本微结构的全尺度多孔结构图；
[0065]
图21为本发明实例2中用提出的方法对优化列式的优化结果。其中，(a)为hh的分布图；(b)为hv的分布图；(c)为h
l
的分布图；(d)为hr的分布图；(e)为
‘
＝’、
‘
||’、
‘
//’、
‘
\\’形基本微结构的全尺度多孔结构图。
[0066]
图22为本发明实例3中的l形模型及其边界条件示意图。
[0067]
图23为本发明实例3中的用提出的方法对优化列式的优化结果。其中，(a)为的分布图；(b)为3d
‘
cubic’形基本微结构的全尺度多孔结构图。
[0068]
图24为本发明实例3中的用提出的方法对优化列式的优化结果。其中，(a)为的分布图；(b)为3d
‘
x’形基本微结构的全尺度多孔结构图。
[0069]
图25为本发明实例3中的用提出的方法对优化列式的优化结果。其中，(a)为的分布图；(b)为的分布图；(c)为3d
‘
cubic’和3d
‘
x’形基本微结构的全尺度多孔结构图。
[0070]
图26为其他三维桁架状和壳状基本微结构示意图。
具体实施方式
[0071]
为了使本发明实施例所要解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
[0072]
需要理解的是，术语“长度”、“宽度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明实施例和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
[0073]
在本发明的描述中，除非另有说明，术语“连接”应做广义理解，例如可以是固定连接，也可以是可拆卸连接，或一体地连接。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
[0074]
实施例：
[0075]
本发明提供一种多孔结构的拓扑优化设计方法，包括如下步骤：
[0076]
步骤一：定义多孔结构微结构原型的水平集函数并进行多种切割、组合生成一系列真实微结构，利用周期边界条件构造了数值基函数，并利用emsfem建立了微结构等效刚度矩阵与设计变量之间的映射关系。此外，本发明还采用数值积分技术计算了微结构的相对密度。所有离线数据都是在优化前计算的。
[0077]
在该步骤中，首先，确定微结构原型。众所周知，任何结构的几何结构都可以用高
维水平集函数隐式表示。如图1所示，一个二维结构ω表示为一个三维水平集函数即:
[0078][0079]
其中x表示域d中任意点的坐标向量，而表示结构的边界。
[0080]
在本实施例中，一个真实的微结构是由几个基本的微结构组成的。如表1所示，选择六种二维和两种三维基本微结构来形成真实的微结构。每个基本微观构可以用对应的基本水平集函数表示。其中，φs，φ
x
，φh，φv，φ
l
，φr，分别对应
‘
box’，
‘
x’，
‘
＝’，
‘
||’，
‘
//’，
‘
\\’,3d
‘
cubic’和3d
‘
x’形的微结构原型。
[0081]
表1.基本微结构原型表
[0082][0083]
另外，相应的基本微结构可以分别用如下水平集函数表示：
[0084][0084][0085]
这里，hb，h
x
，hh，hv，h
l
，hr，是设计变量，即2d问题中的切割曲面的高度及
3d问题中切割曲面的等值线。通过调整切割面高度，可以得到一系列具有不同相对密度的基本微结构。这些基本微观结构可以单独构成真实的微观结构，也可以将它们组合起来得到真实的微观结构。如图2所示，通过将基本微结构进行某种组合生成了四种真实的微结构，其中，(a)由方形和
‘
x’形微结构原型生成；(b)由
‘
＝'和
‘
//'形状的微结构原型生成；(c)由
‘
||’，
‘
//’和
‘
\\'形状的微结构原型生成；(d)由3d
‘
cubic’和3d
‘
x’形微结构原型生成。
[0086]
利用emsfem计算微结构的等效单元刚度矩阵：
[0087]
如上所述，对多孔结构的直接进行全尺寸有限元分析计算极其昂贵。为了提高计算效率，扩大计算规模，在结构优化设计前，采用emsfem对各种真实微结构进行了等效刚度矩阵的计算。
[0088]
在本项工作中，由于我们的微结构都具有周期性的几何特征，因此采用周期边界条件来构造数值基函数。以粗单元dc节点1在x方向的数值基函数的构造过程为例，位移边界条件如图3所示，下标
‘
left’,
‘
right’，
‘
lower’，
‘
upper’分别代表子网格上位于左、右、下、上边界的节点。
[0089]
对于周期边界条件，在两个相对的边界上施加相对位移约束:
[0090][0091]
这里，u3＝0，v3＝0分别用来固定节点3的x和y方向的位移。在施加如上位移条件后，可以在子网格上通过求解弹性力学问题来获得数值基函数，即:这里ks代表子网格上的总体刚度矩阵，us是子网格内的位移向量，fs表示四个边界上指定的位移约束所引起的外力向量。则数值基函数可以表示为:
[0092][0093]
这里，和分别为计算所得位移矢向量的x方向和y方向的分量。同样的，其他节点的数值基函数也可通过施加相应的位移约束计算得到。最终的数值基函数矩阵n组成为:
[0094][0095]
这里表示考虑微观结构内部材料的泊松效应等因素而引起的附加耦合项。利用构造的数值基函数，子网格位移矢量与单个粗单元节点位移矢量的关系可表示为
[0096][0097]
这里us,vs分别表示子网格内x和y方向的位移向量，分别为粗单元dc在x和y方向的节点位移。建立了真实微结构的数值基函数后，通过下式计算真实微结构的等效单
元刚度矩阵：
[0098]
kc＝n
tksn[0099]
其中，n为数值基函数，ks为子网格上的总体刚度矩阵，kc为真实微结构的等效单元刚度矩阵。
[0100]
利用数值积分法计算微结构的相对密度：
[0101]
为便于计算，如图4所示，微结构采用规则的固定背景网格进行网格划分，其中包括三种单元，即空单元、实体单元和混合单元。采用数值积分法计算微结构的相对密度。对于空单元和实体单元，使用2
dim
个积分点，其中dim表示微结构的维度。对于混合单元，积分点数设置为6
dim
，以提高计算精度。
[0102]
建立映射关系数据库：
[0103]
如前所述，一个真实的微结构可以由单一的基本微结构或多个基本微结构组合而成。在本项工作中，为了生成不同的二维和三维真实微结构，共考虑了七种组合形式，即:
[0104]
(1)仅单个的
‘
box’形基本微结构；(2)仅单个的
‘
x’形基本微结构；(3)
‘
box’和
‘
x’形基本微结构的组合；(4)
‘
＝’,
‘
||’,
‘
//’,
‘
\\’四种基本微结构的组合；(5)仅单个的3d
‘
cube’形基本微结构；(6)仅单个的3d
‘
x’形基本微结构；(7)3d
‘
cube’形和3d
‘
x’形基本微结构的组合。
[0105]
在实际工程应用中，也可以构造其他形状的基本微结构。将基本微观组织相结合可以得到真实的微观结构。
[0106]
对于每种组合，计算了相应真实微结构的等效刚度矩阵和相对密度。通过大量的离线计算，建立了不同组合微结构的等效刚度矩阵、相对密度及设计变量之间的映射关系数据库。由于kc矩阵中元素太多，使得等效刚度矩阵与设计变量之间的关系曲线或曲面不容易可视化。在本节中，我们只给出矩阵的第一个元素kc与设计变量之间关系的曲线或曲面。如图5(a),(b),(c)所示，通过第1，2，3个组合形式，我们可以分别建立映射关系kc(hb)，kc(h
x
)和kc(hb,h
x
)，类似的，图6中的映射关系kc(hh,hv,h
l
,hr)对通过
‘
＝’，
‘
||’，
‘
//’，
‘
\\’四种基本微结构组合而成的真实微结构，即第4种组合形式得到。对于3d问题，映射关系和也可以分别通过第5，6,和第7种真实微结构的组合方式得到，如图7(a),(b),(c)所示。
[0107]
同时也用如上所述的七种组合方式建立了相对密度与设计变量的映射关系，分别如图8(a),图8(b),图8(c),图9,图10(a),图10(b),和图10(c)所示。
[0108]
步骤二：离线数据库准备好后，可以基于粗尺度网格进行优化，节省大量计算时间和资源。在粗尺度优化中，基于离线数据库中存储的数据，采用插值和数值求导技术进行结构分析和灵敏度分析。
[0109]
针对不同组合产生的真实微结构，提出了七种不同的优化列式:
[0110]
仅单个的
‘
box’形基本微结构:
[0111][0112]
仅单个的
‘
x’形基本微结构:
[0113][0114]
‘
box’和
‘
x’形基本微结构的组合:
[0115][0116]
‘
＝’,
‘
||’,
‘
//’,
‘
\\’四种基本微结构的组合:
[0117][0118]
仅单个的3d
‘
cube’形基本微结构:
[0119][0120]
仅单个的3d
‘
x’形基本微结构:
[0121][0122]
3d
‘
cube’形和3d
‘
x’形基本微结构的组合:
[0123][0124]
这里ji(i＝1～7)为目标函数；v表示材料使用量；表示材料的使用上限；上标
‘
min’和
‘
max’表示设计变量的下限和上限；fc和uc分别为结构在粗网格上的外力向量和位移向量.uc可以通过在粗尺度网格上求解平衡方程kcuc＝fc确定。
[0125]
灵敏度分析：
[0126]
由于外力向量与设计变量无关，因此目标函数相对于设计变量的导数为
[0127]
[0128]
这里表示第
‘
e’个粗单元的设计变量；由于设计变量与外力向量无关，上标
‘
m’代表表1中相应的基本微结构；kc是总体刚度矩阵；可以由如下公式计算:
[0129][0130]cc
,ce为扩充矩阵，ne为粗单元数目。对于所提出的优化列式，这里的可以从如图5，图6和图7已建立的离线数据库中通过数值求导或插值的手段获得。然后，将公式带回上式，得：
[0131][0132]
这里ue和ke分别为第
‘
e’个粗单元的节点位移向量和刚度矩阵。
[0133]
材料的使用量v可写为：
[0134][0135]
这里vc和分别为粗单元的体积和相应相应微结构的相对密度。则对于第
‘
e’个粗单元，v关于设计变量的导数为：
[0136][0137]
这里可以从如图8，图9和图10等不同优化列式下的离线数据库中通过数值求导及插值技术来获得。
[0138]
此外，采取最优准则法(oc)算法用于更新设计变量，为消除棋盘格图案和网格依赖关系的影响，所有设计变量都使用了密度过滤技术。
[0139]
步骤三：粗尺度优化后，可以得到设计变量的分布。然后，根据步骤二基本微结构设计变量分布和步骤一的水平集函数，对原尺寸多孔结构进行几何表示和重构，最终得到全尺寸多孔结构。
[0140]
在粗尺度优化后，可以得到设计变量的分布。最后，从相应微结构的原型水平集函数中减去设计变量值，即可表示出全尺寸多孔结构的几何形状。
[0141]
重构的过程如图11所示，其中有四个粗单元，每个粗单元考虑两个设计变量h
l
和hr。为了保证界面处相邻微结构的平滑连续性，首先将每个粗元的设计变量通过加权平均策略转换到节点上，如图11所示。然后，利用节点处的设计变量和双线性形状函数插值技术，可形成多个变高度切割曲面。将生成的变高度切割曲面减去对应微结构原型的水平集函数，得到由对应基本微结构组成的梯度多孔结构。如图11所示，通过将这些包含基本微观结构的多孔结构叠加在一起，可以形成最终的多孔结构。
[0142]
为了更清楚地理解本发明提出的方法，我们提供了一个计算流程图，如图12所示，其中包括三个主要部分，即离线数据库的建立、粗尺度优化和全尺寸多孔结构重建。在离线数据库中，利用emsfem和数值积分技术，得到了微结构等效刚度矩阵、微结构相对密度和设计变量之间的映射关系。在粗尺度优化中，采用基于粗尺度网格的emsfem完成结构分析，采
用基于离线数据库的数值求导和插值的手段完成灵敏度分析。因此，整个优化迭代只需要在粗尺度上进行，可以节省大量的计算时间和计算资源。在粗尺度优化后，生成了全尺寸多孔结构的几何形状。
[0143]
实施例1
[0144]
在本例中，悬臂梁模型及其边界条件如图13(a)所示，其中设计域的大小设为l＝40，w＝20，外力设为f＝1。此外，在粗尺度网格中有40
×
20的粗单元。本例中考虑了第四种优化组合，对应的四种基本微观结构如图13(b)所示。材料使用体积与设计域体积之比的上限设为0.5。
[0145]
通过使用所提出的基于emsfem的数据驱动方法，粗尺度上优化后的hh，hv，h
l
,和hr分布分别对应
‘
＝’,
‘
||’，
‘
//’，和
‘
\\’,如图14(a)，(c)，(e)和(g)所示。相应的，也分别显示了仅其中单一基本微观结构重建的结构，如图14(b)，(d)，(f)和(h)所示。如图15(b)所示，将仅含单一基本微结构的所有重构结构叠加在一起即可生成最终的多孔结构。
[0146]
为了便于比较，这里还提供了使用传统m-vcut水平集方法优化的多孔结构，如图15(a)所示。可以看出，本发明提出的基于emsfem的数据驱动方法优化的多孔结构构型与传统的m-vcut水平集方法基本相同，说明了本发明方法的有效性。此外，这两种方法的目标函数和约束函数的迭代历史也绘制在图16中，这里用移动渐近线法(mma)更新设计变量。表2中，n
design
表示设计变量个数，j4表示目标函数，t
total
表示总计算时间，n
step
为迭代步数，ne为网格中用于求解弹性问题和敏感度分析的单元数。
[0147]
表2优化结果数据表
[0148][0149]
通过与传统的m-vcut水平集方法的比较，我们可以看到，本实施例方法在设计变量较少的情况下，可以得到略低的目标函数值。在计算效率方面，我们可以看到，在迭代次数基本相同的情况下，本实施例方法比传统的m-vcut水平集方法消耗的计算时间要少得多。在本例中，与传统的m-vcut方法相比，所提方法的计算效率提高了约99.3％。
[0150]
实施例2
[0151]
如图17所示，这里研究了一个mbb模型，并考虑了四种不同的优化列式如图17所示，这里研究了一个mbb模型，并考虑了四种不同的优化列式和如图17所示的设计域的大小设置为l＝180和w＝30，外力设置为f＝1。模型由180
×
30的粗单元组成。材料使用体积与设计域体积之比的上限设置为0.3。另外，设计变量的下限设为0，上限设为1。
[0152]
对于只具有
‘
box’基本微观结构的优化列式采用本发明方法得到的hb分布和重构的多孔结构分别如图18(a)和图18(b)所示。从局部放大的图像可以看出，相邻微结构之间的平滑连接得到了保证。此外，对于只含’x’形基本微观结构的优化列式优化结果如图19所示。此外，对于优化列式与
‘
box’和
‘
x’形基本微观结构，设计变量hb和h
x
的分布分别显示在图20(a)和图20(b)中，重构的全尺寸多孔结构如图20(c)所示，从图中我们可以看到，基本微观结构的不同组合被安排在不同的位置，使结构的性能更好。同样，
我们可以在优化列式中发现同样的现象，如图21所示。优化列式的优化结果的目标函数值分别为j1＝200.5，j2＝162.3，j3＝158.8，j4＝145.9。可以看出，随着基本微结构数量的增加，目标函数值逐渐变小。当基本微结构数量增多时，由基本微结构生成的真实微结构种类数也将增加，即真实微结构构型的设计空间变大，使得目标函数值变小，计算结果更优。
[0153]
实施例3
[0154]
如图22所示，为了验证本发明方法对三维问题的有效性，我们设置了一个三维l形模型。模型的尺寸参数分别设为l＝30，l1＝20，w＝10。均匀外力设为q＝1。此外，粗尺度网格如图22所示，每个粗尺度单元的大小是1
×1×
1。材料使用体积与设计域体积之比的上限设置为0.3。设计变量的上下限分别设为0.1和0.8。
[0155]
对于优化列式只有三维
‘
cubic’形状的基本微观结构，设计变量的分布如图23(a)所示，对应重构的全尺度多孔结构如图23(b)所示。对于优化列式只带有3d的
‘
x’形基本微结构，优化结果如图24所示。类似地，在图25中给出了优化列式与3d
‘
cubic’和3d
‘
x’形状的基本微结构。
[0156]
优化列式和分别为j5＝3928.8，j6＝4096.2和j7＝3824.5。由此可见，在优化过程中，随着基本微结构数量的增加，目标函数值逐渐变小。这一结论已在算例2二维问题中给出，并在本例中对三维问题进行了验证。
[0157]
可以想象，三维问题的效率将有更大程度的提高。通过计算发现，不同的基本微结构组合布置在不同的位置，使得结构的性能更好。随着基本微结构数量的增加，设计空间变大，导致目标函数值变小。以上仅为本发明较佳的实施例，并非因此限制本发明的实施方式及保护范围，如图26所示的其他类型的基本微观结构，也可以选择在提出的方法框架下优化壳结构或晶格结构。
[0158]
对于本领域技术人员而言，应当能够意识到凡运用本发明说明书内容所作出的等同替换和显而易见的变化所得到的方案，比如在步骤一离线数据库的构建中，可以是扩展多尺度有限元法，也可以是代表体元法或计算均匀化法等等多尺度计算方法。这些方案均应当包含在本发明的保护范围内。