TOPOLOGICAL OPTIMIZATION DESIGN METHOD OF LAYER-WISE GRADED LATTICE STRUCTURES WITH HIGH LOAD-BEARING
-
摘要: 随着增材制造技术的迅速发展, 点阵结构由于其高比强度、高比刚度等优异力学性能受到广泛关注, 但其单胞分布设计大多基于均布式假设, 导致其承载能力相对较差. 基于拓扑优化技术提出了一种梯度分层的点阵结构设计方法. 首先, 基于水平集函数建立点阵单胞几何构型的显式描述模型, 引入形状插值技术实现点阵单胞的梯度构型生成; 其次, 构建基于Kriging的梯度点阵单胞宏观等效力学属性预测模型, 建立宏观有限单元密度与微观点阵单胞等效力学属性的内在联系; 然后, 以点阵结构刚度最大为优化目标, 结构材料用量和力学控制方程为约束条件, 构建点阵结构的梯度分层拓扑优化模型, 并采用OC算法进行数值求解. 算例结果表明, 所提方法可实现点阵结构的最优梯度分层设计, 充分提高了点阵结构的承载性能, 同时可保证不同梯度点阵单胞之间的几何连续性. 最后, 开展梯度分层点阵结构与传统均匀点阵结构和线性梯度点阵结构的准静态压缩仿真分析, 仿真结果表明, 与传统均匀点阵结构和线性梯度点阵结构相比, 梯度分层点阵结构的承载能力明显提高. 研究结果可为高承载点阵结构设计提供理论参考.Abstract: With the rapid development of additive manufacturing technology, lattice structures have attracted extensive attention due to their excellent mechanical properties, such as high specific strength and high specific stiffness. However, the designs of lattice structures are mostly based on the assumption of uniform distribution, resulting in a relatively poor load-bearing capacity. This paper proposes a layer-wise graded lattice structure design method based on a topology optimization technology. Firstly, an explicit description model of lattice geometric configuration is established by using the level set function, and a shape interpolation technology is employed to generate the graded configurations of lattice cells. Secondly, a prediction model of macro effective mechanical property for these graded lattice cells is constructed based on the Kriging metamodel, achieving the essential relationship between the effective density of macro element and the effective mechanical property of micro lattice cell. Then, with the maximum stiffness of lattice structures as the optimization objective, the allowable material usage amount and structural system equilibrium equation as the constraint conditions, a layer-wise graded topology optimization model of lattice structures is established, which is solved numerically by using the OC algorithm. The numerical results indicate that the proposed method can realize the optimal layer-wise graded design of lattice structures, which not only fully improve the load bearing performance of lattice structures, but also ensure the geometric connectivity between different graded lattice cells. Finally, the quasi-static compression simulation analyses of the layer-wise graded lattice structures, the traditional uniform lattice structures and the linear graded lattice structures are carried out and discussed. The simulation results show that, compared with the traditional uniform lattice structures and the linear graded lattice structures, the loading capacity of the layer-wise graded lattice structures is significantly improved. The proposed method provides a theoretical reference for the design of high loading lattice structures.
-
Keywords:
- lattice structure /
- graded layer /
- topology optimization /
- shape interpolation /
- Kriging model
-
引 言
点阵结构是由杆件组成的点阵单胞, 按一定规则重复排列构成的空间桁架结构. 点阵结构因质量轻、比强度高、比刚度高和能量吸收性能优异等力学优势[1-2], 在航空航天、海军舰船和生物医学等诸多领域有着广泛的应用前景[3-5]. 但点阵微结构多变[6]与高孔隙率[7]的特点常导致其不易制备. 近年来, 随着增材制造技术的迅速发展, 极大地方便了具有复杂微观构型的点阵结构制备, 使得具有高力学性能的点阵结构设计与优化成为了研究热点.
点阵结构的力学性能受其微结构构型的影响, 通过优化其微结构的尺寸和形状参数可实现点阵结构力学性能的改进, 这是一种相对简单且有效的方式[8]. Bai等[9]为解决体心立方(body centered cubic, BCC)晶格在节点处应力集中的问题, 通过调整点阵单胞的节点处半径, 实现了单胞杆件的变截面设计, 有限元仿真和实验测试证明新构型点阵的力学性能优势. El-Sayed等[10]采用响应面和方差统计分析方法设计实验, 研究设计参数对八面体点阵结构力学性能的影响, 通过对晶格支柱直径、长度和角度进行尺寸优化, 得到与人类骨骼相似的点阵结构. 王书恒等[11]设计出一种力学性能可调的组合桁架点阵结构, 考虑材料的各向异性, 基于代表体元法, 推导出点阵结构弹性性能的表达式以及弹性各向同性条件. 徐世鹏等[12]提出了一种考虑时变刚度特性的复合材料微结构拓扑优化设计方法, 在不同降解时间步利用均匀化方法计算中间结构的力学性能, 使复合材料为微结构满足时变刚度特性.
上述研究中点阵单胞大多是均匀分布. 但在实际应用中, 点阵结构通常会承受不同的载荷作用. 因此与均匀点阵结构相比, 具有不同梯度分布的点阵结构设计可具有更优的力学性能[13-15]. 如图1的点阵分层结构, 从纵向看, 点阵结构呈梯度分布, 从横向看, 每一层点阵结构密度相同, 以分层的策略设计点阵结构. Dumas等[16]提出了一种适用于骨替代应用的孔隙度分层点阵结构设计方法, 控制结构特性, 如孔径和支柱厚度, 以提供所需的孔隙率分布和结构的力学性能. Liu等[17]研究了生物多孔支架在物理和力学特性的梯度变化要求, 基于三重周期极小曲面提出了Gyroid和Diamond晶格关于密度、异质结构和单元尺寸梯度三种模式下的梯度分层点阵结构. Sanjairaj等[18]提出了一种对三重周期极小曲面片状支架的优化方法, 设计其单元尺寸和壳厚度, 生成梯度分层点阵结构, 以同时满足杨氏模量、孔隙率和孔径的要求. Bai等[19]提出了一种尺寸梯度分层点阵结构设计策略, 利用不同尺寸形成的晶胞单元组装成梯度结构, 通过压缩实验证明尺寸梯度结构连接稳定, 且具有模量高、强度高和能量吸收性能好等优点. Chamini等[20]设计出具有更优坍塌模式、抗冲击以及能量吸收性能的体心立方型双向梯度点阵结构.
然而, 在上述研究中点阵结构的梯度分布形式大多基于设计者的经验确定, 难以实现点阵结构梯度分布与外载荷条件的有效匹配, 以充分发挥材料潜力, 提升结构承载性能. 此外, 在实际工程中, 点阵结构会面临各种复杂载荷工况, 导致其最优梯度分布设计难以通过经验设计获得[21].
拓扑优化通过寻求设计域内材料的最优分布, 力求在满足约束的条件下, 获得性能最优的结构拓扑形式, 是一种先进的结构创新设计方法[22]. 基于此, 本文引入拓扑优化技术, 开展点阵结构的梯度分层设计方法研究. 廖中源等[23]基于均匀化方法计算拟合得到了单胞等效力学属性函数, 一定程度减轻了计算成本, 采用参数化建模实现点阵结构拓扑优化设计; 蔡金虎等[24]采用变密度法对结构进行一次拓扑优化, 得到结构大致的相对密度与应力分布, 并建立了密度、应力与胞元支柱截面尺寸的映射关系, 实现结构的点阵化设计; 赵芳垒等[25]在此基础上提出了一种局部相对密度映射的设计方法, 使得变密度结构具有较好的连接性, 更加有效地保留了拓扑优化信息. 相较于变密度点阵结构优化, 分层设计的策略会使设计自由度降低, 但计算效率大幅提高.
然而, 基于拓扑优化的梯度分层点阵结构优化设计的研究仍然存在一些挑战. 一方面, 计算大量梯度点阵单胞的宏观等效属性需要重复调用均匀化方法(homogenization method, HM)[26], 导致计算成本高昂; 另一方面, 梯度点阵单胞几何构型之间存在阶跃连接, 易产生应力集中, 是导致结构提前失效的主因[27-28]. 因此, 本文构建基于Kriging模型的梯度点阵单胞宏观等效力学属性预测模型, 实现梯度点阵单胞等效属性的高效预测; 建立基于水平集函数的点阵单胞几何构型显式描述模型, 并引入形状插值技术保证梯度点阵单胞几何构型的光滑连接; 构建基于拓扑优化的点阵结构梯度分层优化模型, 并采用OC算法进行求解, 实现点阵结构的最优梯度分层设计. 最后通过数值算例证明所提方法的正确性和有效性, 并与相同约束下的均匀点阵和线性梯度点阵进行对比以验证梯度分层点阵结构的承载优势.
1. 梯度分层点阵结构拓扑优化模型
拓扑优化是在一定的约束条件和设计域内, 寻找材料最优分布的一种方法. 本文以点阵单胞的单元密度为设计变量, 为满足点阵结构的刚度需求, 以柔度最小化作为优化目标函数, 同时满足结构体积约束, 建立的梯度点阵结构拓扑优化模型为
$$ \left.\begin{split} & {\text{Find:}}\qquad{{\boldsymbol{\rho}} }{\text{ = }}\left( {{\rho _1},{\rho _2},\cdots,{\rho _q}} \right),q = 1,2,\cdots,y \\ & {\text{Minimize: }}{{\boldsymbol{C}}}\left( {{\boldsymbol{\rho}} } \right) = {{{\boldsymbol{F}}}^{\text{T}}}{{\boldsymbol{U}}}{\text{ = }}{{{\boldsymbol{U}}}^{\text{T}}}{{\boldsymbol{KU}}} \\ & {\text{Subject to: }}{{{\boldsymbol{G}}}^{\text{S}}}\left( {{\boldsymbol{\rho}} } \right) = {{\boldsymbol{v}}}\left( {{\boldsymbol{\rho }}} \right) - f{v_0} \leqslant 0, \\ & \qquad\qquad\quad {{\boldsymbol{F}} = {\boldsymbol{KU}}}, \\ & \qquad\qquad\quad{\text{0}} \leqslant {\rho _{\min }} \leqslant {\rho _q} \leqslant 1 \end{split}\right\} $$ (1) 式中,
$ {\rho _q} $ 为设计变量(单元相对密度);$ {\rho _{\min }} $ 为最小相对密度; y为设计变量的个数; C为总柔度; F为总体载荷矩阵; U为在载荷F作用下产生的总体位移矩阵; K为点阵结构的总体刚度矩阵;${{\boldsymbol{v}}}\left( {{\boldsymbol{\rho}} } \right)$ , v0分别为优化后实体体积以及设计域总体积; f为允许的材料体积分数.2. 点阵单胞几何构型的数学描述
2.1 二维点阵结构
二维点阵结构由水平集函数[29-30]表示, 用于获得若干密度不同的点阵单元, 二维点阵结构由图2所示的构件组成, 其水平集函数为
$$ \left.\begin{split} &{\boldsymbol{\phi} }_{{\rm{2D}}}\left({{\boldsymbol{x}}}\right) > 0, \;{\rm{if}}\;{{\boldsymbol{x}}}\in {\boldsymbol{\varOmega}} \\ &{\boldsymbol{\phi} }_{{\rm{{\rm{2D}}}}}\left({{\boldsymbol{x}}}\right) = 0, {\rm{if}}\;{{\boldsymbol{x}}}\in \partial {\boldsymbol{\varOmega}} \\ &{\boldsymbol{\phi} }_{{\rm{2D}}}\left({{\boldsymbol{x}}}\right) < 0, {\rm{if}}\;{{\boldsymbol{x}}}\in {{\boldsymbol{D}}}\backslash {\boldsymbol{\varOmega}}\\ &{\boldsymbol{\phi} }_{{\rm{2D}}}\left({{\boldsymbol{x}}}\right) = \mathrm{max}\left({\boldsymbol{\phi} }_{{\rm{2D}},{r}}\left(x,y\right),{\boldsymbol{\phi} }_{{\rm{2D}},c1}\left(x,y\right),{\boldsymbol{\phi} }_{{\rm{2D}},c2}\left(x,y\right)\right)\\ &{\boldsymbol{\phi} }_{{\rm{2D}},{r}}\left(x,y\right) = \mathrm{min}\left({\boldsymbol{\phi} }_{{\rm{2D}},r1}\left(x,y\right),{\boldsymbol{\phi} }_{{\rm{2D}},r2}\left(x,y\right)\right)\\ &{\boldsymbol{\phi} }_{{\rm{2D}},{r}1}\left(x,y\right) = 1-\frac{{\left[-\mathrm{sin}{\theta}_{{\rm{2D}}}\cdot \left(x-{x}_{0}\right) + \mathrm{cos}{\theta}_{{\rm{2D}}}\cdot \left(y-{y}_{0}\right)\right]}^{2}}{{\left({t}_{{\rm{2D}}}/2\right)}^{2}}\\ &{\boldsymbol{\phi} }_{{\rm{2D}},{r}2}\left(x,y\right) = 1-\frac{{\left[\mathrm{cos}{\theta}_{{\rm{2D}}}\cdot \left(x-{x}_{0}\right) + \mathrm{sin}{\theta}_{{\rm{2D}}}\cdot \left(y-{y}_{0}\right)\right]}^{2}}{{\left({L}_{{\rm{2D}}}/2\right)}^{2}}\\ &\mathrm{sin}{\theta}_{{\rm{2D}}} = \frac{{y}_{2}-{y}_{1}}{{L}_{{\rm{2D}}}}, \mathrm{cos}{\theta}_{{\rm{2D}}} = \frac{{{x}}_{2}-{x}_{1}}{{L}_{{\rm{2D}}}}\\ &{x}_{0} = \frac{{{x}}_{1} + {x}_{2}}{2},{y}_{0} = \frac{{{y}}_{1} + {y}_{2}}{2},{L}_{{\rm{2D}}} = \sqrt{{\left({{x}}_{2}-{x}_{1}\right)}^{2} + {\left({y}_{2}-{y}_{1}\right)}^{2}}\\ &{\boldsymbol{\phi} }_{{\rm{2D}},c1}\left(x,y\right) = 1-\frac{{\left(x-{x}_{1}\right)}^{2} + {\left(y-{y}_{1}\right)}^{2}}{{\left({t}_{{\rm{2D}}}/2\right)}^{2}}\\ &{\boldsymbol{\phi} }_{{\rm{2D}},c2}\left(x,y\right) = 1-\frac{{\left(x-{x}_{2}\right)}^{2} + {\left(y-{y}_{2}\right)}^{2}}{{\left({t}_{{\rm{2D}}}/2\right)}^{2}}\end{split}\right\} $$ 式中,
${{\boldsymbol{D}}}$ 代表二维空间的结构设计域,${{\boldsymbol{\varOmega}} }$ 代表存在二维棒状实体的区域,$\partial {{\boldsymbol{\varOmega}} }$ 则为实体区域的边界, 棒状部分由一个矩形和两个圆组成, 水平集函数分别为${{\boldsymbol{\phi}} _{{\rm{2D}},r}}\left( {x,y} \right),$ ${{\boldsymbol{\phi}} _{{\rm{2D}},{{c}}1}}\left( {x,y} \right)$ 和${{\boldsymbol{\phi}} _{{\rm{2D}},{{c2}}}}\left( {x,y} \right),$ ${{\boldsymbol{\phi}} _{{\rm{2D}}}}\left( {x,y} \right)$ 取各组成部分${\boldsymbol{\phi}}$ 的最大值; 矩形的两端以$ \left( {{x_1},{y_1}} \right) $ 和$ \left( {{x_2},{y_2}} \right) $ 为中心, 这两个点的坐标也是两个圆的圆心坐标, 即矩形的厚度t2D与圆的直径相等.$ {\theta _{{\rm{2D}}}} $ 为矩形的倾斜角, 角度从水平面逆时针开始测量.$ \left( {{x_0},{y_0}} \right) $ 和L2D分别为矩形的中点坐标和长度.通过基于水平集函数的显式描述模型描述点阵单胞, 便于宏观与微观两个尺度的转化, 图3是BCC晶格的三维水平集函数以及体积分数为0.5的水平面. 对水平集函数进行不同水平面的插值, 得到100个BCC单胞的样本数据, 用于后续Kriging模型构建. 基于水平集函数的形状插值技术使所有梯度点阵均具有相似的拓扑构型, 保证了点阵单胞间的连接性[31].
2.2 三维点阵结构
三维点阵结构由图4所示的构件组成, 其水平集函数由式(3)表示为
$$ \left.\begin{split} &{\boldsymbol{\phi} }_{{\rm{3D}}}\left({{\boldsymbol{x}}}\right) > 0\text{, }{\rm{if}}\;{{\boldsymbol{x}}}\in {\boldsymbol{\varOmega}}\\ &{\boldsymbol{\phi} }_{{\rm{{\rm{3D}}}}}\left({{\boldsymbol{x}}}\right) = 0\text{, }{\rm{if}}\;{{\boldsymbol{x}}}\in \partial {\boldsymbol{\varOmega}}\\ &{\boldsymbol{\phi} }_{{\rm{3D}}}\left({{\boldsymbol{x}}}\right) < 0\text{, }{\rm{if}}\;{{\boldsymbol{x}}}\in {{\boldsymbol{D}}}\backslash {\boldsymbol{\varOmega}}\\ &{\boldsymbol{\phi} }_{{\rm{3D}}}\left({{\boldsymbol{x}}}\right) = \mathrm{max}\left({\boldsymbol{\phi} }_{{\rm{3D}},r}\left(x,y,z\right),{\boldsymbol{\phi} }_{{\rm{3D}},c1}\left(x,y,z\right),{\boldsymbol{\phi} }_{{\rm{3D}},c2}\left(x,y,z\right)\right)\\ &{\boldsymbol{\phi} }_{{\rm{3D}},r}\left(x,y,z\right) = \mathrm{min}\left({\boldsymbol{\phi} }_{{\rm{3D}},r1}\left(x,y,z\right),{\boldsymbol{\phi} }_{{\rm{3D}},r2}\left(x,y,z\right)\right)\\ &{\boldsymbol{\phi} }_{{\rm{3D}},r1}\left(x,y,z\right) = {\left({L}_{{\rm{3D}}}/2\right)}^{2}-{\left(\mathrm{cos}{\theta}_{{\rm{3D}}}\cdot {L}_{{\rm{dis}}}\right)}^{2}\\ &{\boldsymbol{\phi} }_{{\rm{3D}},r2}\left(x,y,z\right) = {\left({t}_{{\rm{3D}}}/2\right)}^{2}-{\left(\mathrm{sin}{\theta}_{{\rm{3D}}}\cdot {L}_{\mathrm{{\rm{dis}}}}\right)}^{2}\\ &\mathrm{cos}{\theta}_{{\rm{3D}}} = \sqrt{{\left(\frac{{d}_{{x}_{2}}{d}_{x} + {d}_{{y}_{2}}{d}_{y} + {d}_{{z}_{2}}{d}_{z}}{{L}_{{\rm{dis}}}\sqrt{{d}_{{x}_{2}}{}^{2} + {d}_{{y}_{2}}{}^{2} + {d}_{{z}_{2}}{}^{2}}}\right)}^{2}}\\ &\mathrm{sin}{\theta}_{{\rm{3D}}} = \sqrt{1-{\mathrm{cos}}^{2}{\theta}_{{\rm{3D}}}},{L}_{{\rm{dis}}} = \sqrt{{d}_{x}{}^{2} + {d}_{y}{}^{2} + {d}_{z}{}^{2}}\\ &{d}_{{x}_{2}} = {x}_{2}-{x}_{0},{d}_{{y}_{2}} = {y}_{2}-{y}_{0},{d}_{{z}_{2}} = {z}_{2}-{z}_{0}\\ &{d}_{x} = x-{x}_{0},{d}_{y} = y-{y}_{0},{d}_{{z}} = z-{z}_{0},\\ &{x}_{0} = \frac{{{x}}_{1} + {x}_{2}}{2},{y}_{0} = \frac{{{y}}_{1} + {y}_{2}}{2},{z}_{0} = \frac{{z}_{1} + {z}_{2}}{2}\\ &{L}_{{\rm{3D}}} = \sqrt{{\left({{x}}_{2}-{x}_{1}\right)}^{2} + {\left({y}_{2}-{y}_{1}\right)}^{2} + {\left({z}_{2}-{z}_{1}\right)}^{2}}\\ &{\boldsymbol{\phi} }_{{\rm{3D}},c1}\left(x,y,z\right) = {\left({t}_{{\rm{3D}}}/2\right)}^{2}-\left[{\left(x-{x}_{1}\right)}^{2} + {\left(y-{y}_{1}\right)}^{2} + {\left(z-{z}_{1}\right)}^{2}\right]\\ &{\boldsymbol{\phi} }_{{\rm{3D}},c2}\left(x,y,z\right) = {\left({t}_{{\rm{3D}}}/2\right)}^{2}-\left[{\left(x-{x}_{2}\right)}^{2} + {\left(y-{y}_{2}\right)}^{2} + {\left(z-{z}_{2}\right)}^{2}\right]\end{split}\right\} $$ (3) 式中,
${{\boldsymbol{D}}}$ 代表三维空间的结构设计域,${{\boldsymbol{\varOmega }}}$ 代表存在三维棒状实体的区域,$\partial {{\boldsymbol{\varOmega}} }$ 为三维实体区域的边界, 棒状部分由一个圆柱体和两个球体组成, 水平集函数分别为${\boldsymbol{\phi} _{{\rm{3D}},r}}\left( {x,y,z} \right)$ ,${\boldsymbol{\phi} _{{\rm{3D}},{c} 1}}\left( {x,y,z} \right)$ 和${\boldsymbol{\phi} _{{\rm{3D}},{c} 2}}\left( {x,y,z} \right)$ ,${\boldsymbol{\phi} _{{\rm{3D}}}}\left( {x,y,z} \right)$ 取各组成部分$ \boldsymbol{\phi} $ 的最大值; 圆柱体的末端以$ \left( {{x_1},{y_1},{z_1}} \right) $ 和$ \left( {{x_2},{y_2},{z_2}} \right) $ 为中心, 这两个点同样是两个球体的球心, 即圆柱体和两个球体的直径具有相同的值t3D.$ \left( {{x_0},{y_0},{z_0}} \right) $ 和L3D分别为圆柱体的中心点坐标和长度. 图5是BCC晶格的四维水平集函数以及体积分数为0.5的等值面.本文通过插值点阵单胞的水平集函数构建样本点阵单胞
$$ {{\boldsymbol{\phi}} ^{{e}}}\left( {{\boldsymbol{x}}} \right) = {{\boldsymbol{\phi}} ^{{{\rm{pro}}} }}\left( {{\boldsymbol{x}}} \right) - {{\varphi }} $$ (4) 式中,
$\varphi \in \left[\mathrm{min}\left({{\boldsymbol{\phi}} }^{\text{pro}}\left({{\boldsymbol{x}}}\right)\right), \text{max}\left({{\boldsymbol{\phi}} }^{\text{pro}}\left({{\boldsymbol{x}}}\right)\right)\right]$ ,$ \varphi $ 为形状插值系数, 可由二分法计算获得.${{\boldsymbol{\phi}} ^{{\text{pro}}}}$ 为点阵单胞的水平集函数, 而${{\boldsymbol{\phi}} ^e}$ 为第e层点阵单元的水平集函数. 随着${\boldsymbol{\phi}}$ 的变化, 点阵单元的结构域也是变化的, 这意味着点阵单元的体积分数与${\boldsymbol{\phi}}$ 相关. 考虑到水平集函数和形状插值适用于各种点阵, 本文选取等效密度为设计变量. 考虑到增材制造的实际精度限制, 本文规定点阵单元密度的范围从0.05到1. 样本点阵单元的等效密度是一个等差数列, 并且不同密度的点阵单元模型参数都可以通过式(4)获得.2.3 基于Kriging模型的梯度点阵单胞宏观等效属性预测
对于大规模变密度拓扑优化问题, 常将实体划分为均匀点阵单胞, 当采用传统的HM求解点阵单胞的等效力学属性时, 需对不同密度下的点阵单元进行分析, 因而导致其计算效率较低. 本文通过构建Kriging模型[32-34]预测不同密度下点阵单元的有效弹性模量, 可大幅提高计算效率.
本文作者在以往的工作中[35-36]验证过当梯度出现时, 基于周期性边界条件得到的等效张量是准确的. 根据HM, 点阵结构的等效弹性张量计算公式为
$$ \begin{split} & {\left( {{{\boldsymbol{D}}}_{ijkl}^{\text{H}}} \right)_e} = \frac{1}{{\left| {{{{\boldsymbol{\varOmega}} }_e}} \right|}} \times \\ & \qquad \int_{{{{\boldsymbol{\varOmega}} }_e}} {\left( {{{\boldsymbol{\varepsilon}} }_{pq}^{0\left( {ij} \right)} - {{\boldsymbol{\varepsilon}} }_{pq}^ * \left( {{{\boldsymbol{u}}}_e^{ij}} \right)} \right)} {{{\boldsymbol{D}}}_{{{pqrs}}}}\left( {{{\boldsymbol{\varepsilon}} }_{rs}^{0\left( {kl} \right)} - {{\boldsymbol{\varepsilon}} }_{rs}^ * \left( {{{\boldsymbol{u}}}_e^{kl}} \right)} \right){{ {\boldsymbol{H}}}}\left( {{{{\boldsymbol{\phi}} }_e}} \right){\rm{d}}{{{\boldsymbol{\varOmega}} }_e} \end{split} $$ (5) 式中, i, j, k和l的值为正整数1, 2, ···, d, d为空间维度; Dpqrs为固体材料的弹性张量;
$ {{{\boldsymbol{\varOmega}} }_e} $ 为设计域,$ \left| {{{{\boldsymbol{\varOmega}} }_e}} \right| $ 为第e层点阵结构的面积或体积. 单位脉冲函数${{ {\boldsymbol{H}}}}\left( {{{{\boldsymbol{\phi}} }_e}} \right)$ 表示设计域的不同部分. 在给定的宏观应变${{\boldsymbol{\varepsilon}} }_{{{pq}}}^{0\left( {ij} \right)}$ 下, 位移矩阵对应的应变矩阵${{\boldsymbol{\varepsilon}} }_{{{pq}}}^ * \left( {{{\boldsymbol{u}}}_e^{\left( {ij} \right)}} \right)$ 可以在微观尺度上求解以下公式来计算$$ \left.\begin{split} & \int_{{{{{{\boldsymbol{\varOmega}}}} }_e}} {\left( {{{\boldsymbol{\varepsilon}} }_{{{pq}}}^{0\left( {ij} \right)} - {{\boldsymbol{\varepsilon}} }_{{{pq}}}^ * \left( {{{\boldsymbol{u}}}_e^{ij}} \right)} \right)} {{{\boldsymbol{D}}}_{{{pqrs}}}}{{\boldsymbol{\varepsilon}} }_{{{rs}}}^ * \left( {{{\boldsymbol{v}}}_e^{kl}} \right){{ {\boldsymbol{H}}}}\left( {{{{\boldsymbol{\phi}} }_e}} \right){\rm{d}}{{{\boldsymbol{\varOmega}} }_e} = 0 \\ & \forall {{\boldsymbol{v}}}_e^{kl} \in \overline {{\boldsymbol{U}}} \left( {{{{\boldsymbol{\varOmega}} }_e}} \right) \end{split}\right\}$$ (6) 式中,
${{\boldsymbol{v}}}_e^{kl}$ 为虚拟位移矩阵,$\overline {{\boldsymbol{U}}} \left( {{{{\boldsymbol{\varOmega}} }_e}} \right)$ 为在周期性边界条件下符合运动学规律的第e层点阵结构的位移场. 计算出若干个点阵结构的等效弹性张量, 将其存储为训练数据, 基于这些训练数据, Kriging模型可以表示为$$ {{{\boldsymbol{D}}}^{\text{H}}}\left( \rho \right) = {{\boldsymbol{f}}}{\left( \rho \right)^{\text{T}}}{{\boldsymbol{\beta}} + {\boldsymbol{Z}}}\left( \rho \right),\;\;\rho \in \left[ {0.05,1} \right] $$ (7) 式中,
${{\boldsymbol{f}}}{\left( \rho \right)^{\text{T}}}{{\boldsymbol{\beta}} }$ 与${{\boldsymbol{Z}}}\left( \rho \right)$ 分别是全局响应估计值和局部偏差;${{\boldsymbol{f}}}\left( \rho \right)$ 为一个ne × 1的回归函数向量, 其中ne表示样本点阵结构的数量;${{\boldsymbol{\beta}} }$ 是一阶多项式回归模型的回归系数向量;$ \rho $ 表示预测点的样本点阵结构等效密度;${{{\boldsymbol{D}}}^{\text{H}}}\left( \rho \right)$ 为样本点阵结构的等效弹性张量. 最后, 得到基于Kriging模型计算点阵结构的等效弹性张量计算公式, 即$$ {{\hat {\boldsymbol{D}}}^{\text{H}}} \sim N\left( {{{{\hat {\boldsymbol{D}}}}^{\text{H}}}\left( \rho \right),{{{\hat {\boldsymbol{\sigma}} }}_{{{{\hat {\boldsymbol{D}}}}^{\text{H}}}}}\left( \rho \right)} \right) $$ (8) $$ {{\boldsymbol{D}}}_{ijkl}^{\text{H}}\left( \rho \right) = {{\hat {\boldsymbol{D}}}^{\text{H}}}\left( \rho \right) $$ (9) Kriging模型可以通过点阵结构的等效密度获得相应密度下的等效弹性张量, 避免在每次迭代中调用HM对点阵结构有效性能进行优化和运算, 使设计方法可以在较低计算负担下实现梯度分层点阵拓扑优化设计.
2.4 梯度分层点阵结构灵敏度分析
点阵结构的刚度矩阵可以表示为
$$ {{\boldsymbol{K}}} = \sum\limits_{q = 1}^{Ne} {\int_{{{{\boldsymbol{\varOmega}} }_q}} {{{{\boldsymbol{B}}}^{\text{T}}}} } {{\boldsymbol{D}}}_{ijkl}^{{\text{ele}}}\left( {\rho _q^{{\text{ele}}}} \right){{\boldsymbol{B}}}{\rm{d}}{{\boldsymbol{\varOmega}} }_q^{{\text{ele}}} $$ (10) 式中, B为应变−位移矩阵,
${{\boldsymbol{\varOmega}} }_q^{{\text{ele}}}$ 代表第q个单元的定义域, 其含有Ne个单元,${{\boldsymbol{D}}}_{ijkl}^{{\text{ele}}}\left( {\rho _q^{{\text{ele}}}} \right)$ 表示第q个单元在该单元密度下的弹性张量, 可以用如下插值法求得$$ {{\boldsymbol{D}}}_{ijkl}^{{\text{ele}}}\left( {\rho _q^{{\text{ele}}}} \right) = \rho _q^{{\text{ele}}}{\tilde {\boldsymbol{D}}}_m^0{\left( {\rho _q^{{\text{ele}}}} \right)_{ijkl}} $$ (11) 式中,
${\tilde {\boldsymbol{D}}}_m^0$ 是用于识别第m个单元的等效力学属性的变量.基于HM,
${{\boldsymbol{D}}}_{ijkl}^{{\text{ele}}}\left( {\rho _q^{{\text{ele}}}} \right)$ 可以表示为$$ {{\boldsymbol{D}}}_{ijkl}^{{\text{ele}}}\left( {\rho _q^{{\text{ele}}}} \right) = {{\boldsymbol{D}}}_{ijkl}^{\text{H}}\left( {\rho _q^{{\text{ele}}}} \right) $$ (12) 式中,
${{\boldsymbol{D}}}_{ijkl}^{\text{H}}\left( {\rho _q^{{\text{ele}}}} \right)$ 可以由式(7)和式(8)预测所得. 根据式(1)和式(10),${\tilde {\boldsymbol{D}}}_m^0{\left( {\rho _q^{{\text{ele}}}} \right)_{ijkl}}$ 可以由下式得到$$ {\tilde {\boldsymbol{D}}}_m^0{\left( {\rho _q^{{\text{ele}}}} \right)_{ijkl}} = {{{{{\hat {\boldsymbol{D}}}}^{\text{H}}}\left( {\rho _q^{{\text{ele}}}} \right)} \mathord{\left/ {\vphantom {{{{{\hat {\boldsymbol{D}}}}^{\text{H}}}\left( {\rho _q^{{\text{ele}}}} \right)} {\rho _q^{{\text{ele}}}}}} \right. } {\rho _q^{{\text{ele}}}}} $$ (13) 值得注意的是,
${\tilde {\boldsymbol{D}}}_m^0{\left( {\rho _q^{{\text{ele}}}} \right)_{ijkl}}$ 的作用在于更新$ \rho _q^{{\text{ele}}} $ , 并按单元由式(11)计算得到.研究关于柔度的点阵结构密度优化问题, 目标函数对设计变量的灵敏度[33]可以表示为
$$ \frac{{\partial {{\boldsymbol{C}}}\left( {{\boldsymbol{\rho }}} \right)}}{{\partial {{{\boldsymbol{\rho}} }_q}}} = mean\left( {{{\boldsymbol{\zeta}} },h} \right) $$ (14) $$ {{\boldsymbol{\zeta}} } = - {{{\boldsymbol{U}}}^{\text{T}}}\frac{{\partial {{\boldsymbol{K}}}\left( {{\boldsymbol{\rho }}} \right)}}{{{{{\boldsymbol{\rho}} }_q}}}{{\boldsymbol{U}}} $$ (15) 式中,
${{\boldsymbol{\zeta}} }$ 表示点阵结构的柔度对单元等效密度的灵敏度,$mean\left( {{{\boldsymbol{\zeta}} },h} \right)$ 表示同一元素层内灵敏度的平均值, 其中h表示元素层的方向. 基于构建的Kriging模型, 全局刚度矩阵的导数可以联立式(1)、式(10) ~ 式(12)得到$$ \frac{{\partial {{\boldsymbol{K}}}\left( {{\boldsymbol{\rho }}} \right)}}{{{{{\boldsymbol{\rho}} }_q}}} = \sum\limits_{q = 1}^y {\int_{{{{\boldsymbol{\varOmega}} }_q}} {{{{\boldsymbol{B}}}^{\text{T}}}} } \left( {\frac{{{{{\hat {\boldsymbol{D}}}}^{\text{H}}}\left( {{{{\boldsymbol{\rho}} }_q}} \right)}}{{{{{\boldsymbol{\rho}} }_q}}}} \right){{\boldsymbol{B}}}{\rm{d}}{{{\boldsymbol{\varOmega}} }_q} $$ (16) 因此, 式(14)可以改写为
$$ {{\boldsymbol{\zeta}} } = \frac{{\partial {{\boldsymbol{C}}}\left( {{\boldsymbol{\rho}} } \right)}}{{\partial {{{\boldsymbol{\rho}} }_q}}} = - {{{\boldsymbol{U}}}^{\text{T}}}\sum\limits_{q = 1}^y {\int_{{{{\boldsymbol{\varOmega}} }_q}} {{{{\boldsymbol{B}}}^{\text{T}}}} } \left( {\frac{{{{{\hat {\boldsymbol{D}}}}^{\text{H}}}\left( {{{{\boldsymbol{\rho}} }_q}} \right)}}{{{{{\boldsymbol{\rho}} }_q}}}} \right){{\boldsymbol{B}}}{\rm{d}}{{{\boldsymbol{\varOmega}} }_q}{{\boldsymbol{U}}} $$ (17) 最后可以计算得到体积约束对单元等效密度
${{{\boldsymbol{\rho }}}_q}$ 的灵敏度$$ \frac{{\partial {{{\boldsymbol{G}}}^{\text{S}}}\left( {{\boldsymbol{\rho}} } \right)}}{{\partial {{{\boldsymbol{\rho}} }_q}}} = {{V}} $$ (18) 3. 算法优化流程
梯度分层点阵结构的拓扑优化设计流程如图6所示. 由流程图可得知, 实现过程主要分为两大步骤: 步骤一是构建Kriging模型预测梯度点阵单胞宏观等效力学属性; 步骤二是梯度分层点阵结构等效密度优化. 最终迭代的结果需要对梯度点阵结构连接问题进行优化: 首先将点阵密度转化为节点密度, 利用加权平均式优化节点密度, 保证连通性[23], 然后根据式(4)得到对应等效密度下的点阵单元, 最后输出点阵结构尺寸.
4. 数值算例及结果分析
本文分别对二维和三维结构开展了优化设计, 以点阵结构刚度最大为优化目标, 用Matlab编程实现所提出的梯度点阵结构分层设计方法. 本文所采用板状二维模型、MBB三维模型、柱状三维模型是为了说明所提方法适用范围广泛, 进而分别从二维结构、三维结构以及对不同载荷方向来验证方法的有效性. 已有的拓扑优化梯度分层点阵设计主要有均匀、线性梯度等形式, 故本文选取均匀与线性梯度两种典型点阵结构, 在同一体积分数条件下, 与拓扑优化所得的梯度结构进行对比.
在本节所考虑的数值算例中, 实体材料的杨氏模量E0设置为20.1 MPa, 材料的泊松比
$ \mu $ = 0.3, 灵敏度过滤半径为1.5, 惩罚因子为1, 设计变量初始值设置为目标体积分数. 当连续两次迭代之间各单元等效密度的最大变化小于0.001或者达到最大迭代步数时, 将停止优化过程.4.1 二维梯度点阵结构
二维梯度点阵结构的设计域、边界条件以及载荷如图7所示. 设计域尺寸为: L = 30 cm, H = 15 cm, 均布载荷F的大小为10 kN, 方向竖直向下, 施加在设计域上表面, 设计域下表面为固定约束. 将二维设计域离散成450个四节点单元, 规定体积分数为0.3, 进行梯度分层点阵结构的拓扑优化迭代, 目标函数、体积分数随迭代步数的变化曲线如图8所示. 可以看到, 优化在第20步达到收敛, 且在迭代9次以后, 已经较为接近收敛的目标函数值, 收敛过程稳定且快速. 将优化得到的梯度结构与相同体积约束下的线性梯度结构和均匀结构进行比较, 其结果如表1所示. 经过有限元分析, 可得拓扑优化梯度点阵结构柔度相比均匀结构和线性梯度结构分别降低了51.27%和31.85%.
取设计域L的值为15 cm, 设计域 L与H的比值将变为1∶1, 保持其余参数和边界条件不变. 将设计域离散成225个四节点单元, 初始均匀点阵结构的柔度为293.04 N·cm. 经过20次迭代后, 结构柔度收敛至172.34 N·cm, 结构柔度降低了约41.19%. 根据表2所示的三种点阵构型比较可知, 拓扑优化点阵结构呈现出明显的梯度分布, 材料分布在主传递力路径上, 承载性能得到极大提升.
表 1 设计域30 cm × 15 cm结构图及柔度Table 1. The structure and compliance of design domain 30 cm × 15 cmUniform Linear graded Topology graded Relative density structure compliance 124.24 N·cm 88.84 N·cm 60.55 N·cm 表 2 设计域15 cm × 15 cm结构图及柔度Table 2. The structure and compliance of design domain 15 cm × 15 cmUniform Linear graded Topology graded Relative density structure compliance 293.04 N·cm 229.45 N·cm 172.34 N·cm 4.2 MBB三维梯度点阵结构
利用所提方法对如图9所示的MBB三维结构进行优化. 设计域尺寸为: L = 40 cm, H = 15 cm, W = 10 cm, 线载荷F的大小为100 kN, 方向沿Z轴方向加载在上表面中线, 设计域下表面的四个端点为固定约束. 将设计域离散成6000个六面体单元, 体积约束为0.3, 迭代收敛条件与上小节相同.
目标函数和体积约束的迭代过程如图10所示. 可以看到, 目标函数和体积约束曲线经过15次迭代后收敛, 在迭代5次以后, 已经较为接近收敛的目标函数值, 收敛过程快速稳定. 将优化结果与相同体积约束下的线性梯度结构和均匀结构作比较, 结果如表3所示. 在相同体积分数条件下, 拓扑优化梯度相较于均匀分布和线性梯度的结构柔度分别降低51.04%和63.81%. 经过优化设计的拓扑优化点阵结构, 呈现明显的梯度分层特征且具有良好的连续性, 同时证明了所提方法对于三维结构同样有效.
表 3 设计域40 cm × 15 cm × 10 cm结构图及柔度Table 3. The structure and compliance of design domain 40 cm × 15 cm × 10 cmUniform Linear graded Topology graded Relative density structure compliance 10171.05 N·cm 13759.11 N·cm 4979.70 N·cm 4.3 三维柱状梯度点阵结构
采用BCC单胞对如图11所示的三维柱状结构进行拓扑优化设计, 设计域尺寸为: L = 15 cm, H = 20 cm, W = 10 cm, 在上表面施加大小为100 kN的斜向面载荷, 平行于XOY平面且与X轴夹角呈45°, 下表面为固定约束, 体积约束设置为0.2. 经过迭代后得到点阵密度, 基于水平集函数的形状插值函数使梯度点阵拓扑结构相似, 利用节点密度优化点阵结构连接, 保证了不同单元密度点阵之间的连接, 得到的立体构型以及微结构如图12所示.
利用点阵单胞几何构型显式描述模型实现点阵结构的生成, 分别得到相同体积分数的均匀点阵、线性梯度点阵和拓扑优化梯度点阵的立体构型图, 结果如表4所示.
表 4 设计域15 cm × 20 cm × 10 cm结构图与柔度Table 4. The structure and compliance of design domain 15 cm × 20 cm × 10 cmUniform Linear graded Topology graded Relative density structures compliance 8165.04 N·cm 7708.78 N·cm 6945.11 N·cm 经过有限元分析, 可以得知所提方法的拓扑优化梯度点阵结构在均布斜载荷下的柱状设计域中得到的柔度值最小, 相较于均匀、线性梯度点阵结构分别降低14.94%和9.91%. 该算例表明, 所提方法可以实现基于不同载荷方向对三维结构的优化设计, 同时满足结构材料用量和力学控制方程. 但显然, 设计域和载荷方向同样影响梯度结构性能优化的效果.
4.4 计算效率对比
对于梯度分层点阵结构的柔度最小化问题, 采用Kriging模型与HM的方法计算迭代步时间, 结果如图13所示. 计算迭代时间所使用的计算机CPU为intel i5-8300 H 2.30 GHz, 内存为8 GB Micron 2667 MHz. 可以看出, 基于Kriging模型的方法每一步迭代所用的时间均少于HM, 在设计域30 cm × 15 cm的算例下, 基于HM的方法迭代总时间为59.83 s, 而基于Kriging模型的方法迭代总时间为24.02 s, 计算效率提高了约59.88%, 证明了采用Kriging大幅提高了所提方法的计算效率, 减轻了设计域单元数的增加带来计算量急剧上升的负担.
值得注意的是, 设计域尺寸15 cm × 20 cm × 10 cm基于HM的迭代步时间略高于设计域40 cm × 15 cm × 10 cm的迭代步时间, 而采用Kriging模型的结果却相反. 在迭代过程中, 影响HM迭代时间的因素主要为设计变量的个数, 而构建基于Kriging的梯度点阵单胞宏观等效力学属性预测模型, 大幅缩短了反复调用HM的计算时间, 故影响迭代时间的因素主要为单胞数量, 所以导致上述结果的出现.
5. 仿真算例
为了进一步验证本文提出的拓扑优化梯度点阵结构设计方法的有效性, 分别利用二维和三维算例对三种不同梯度的点阵结构在COMSOL中进行仿真分析. 在下述各例中, 假设模型为线弹性材料, 材料属性初始屈服应力为235 MPa, 各向同性切线模量为6100 MPa.
对表2的三种点阵结构进行压缩试验仿真, 模型的尺寸按比列放大为750 mm × 750 mm, 其中每个点阵单胞的尺寸为50 mm × 50 mm, 模型下表面为固定约束, 改变上表面应力, 对点阵结构进行应变求解, 分别得到三种方案的应力云图如图14所示. 对上表面应变求平均值可得其应力应变曲线如图15所示. 结果显示, 拓扑优化梯度点阵相较于其他两种梯度点阵整体受力均匀, 更符合力学性能要求, 由应力应变曲线可得, 拓扑优化梯度在同一应变下, 承载能力明显优于其余两种梯度, 与有限元分析结果一致, 验证了本文所提出优化方法对二维点阵结构的有效性.
取设计域尺寸为3 cm × 4 cm × 2 cm的三维点阵结构进行压缩试验仿真, 模型的尺寸等比例放大为90 mm × 120 mm × 60 mm, 其中每个点阵单胞的尺寸为30 mm × 30 mm × 30 mm, 模型下表面为固定约束, 对上表面施加斜载荷, 对点阵结构进行应力求解, 分别得到三种梯度点阵的应力云图如图16所示. 对上表面沿受力方向的应变计算可得应力应变曲线如图17所示. 通过有限元分析计算得出, 均匀点阵柔度为22258.44N·cm, 线性梯度点阵为20906.23 N·cm, 而拓扑优化点阵柔度为18792.69 N·cm, 柔度分别降低了15.57%和10.11%. 由应力应变曲线可得, 三种梯度在小应变下, 应力应变曲线相近, 随着应力不断增加, 拓扑优化梯度显示出更好的承载能力, 与有限元分析结果一致, 验证了本文所提出的拓扑优化理论方法对三维点阵结构的有效性.
6. 结论
(1)本文基于拓扑优化技术提出了点阵结构的梯度分层设计方法, 以点阵结构的刚度最大为优化目标, 结构材料用量和力学控制方程为约束条件, 构建梯度分层点阵结构的拓扑优化模型, 采用OC算法求解, 获得了具有最优梯度分层的点阵结构设计. 数值算例和仿真结果表明, 相比传统均匀梯度与线性梯度点阵结构设计, 梯度分层点阵结构整体刚度得到提升, 显现出显著的力学性能优势, 为点阵结构设计提供了理论参考.
(2)建立了基于水平集函数的点阵单胞几何构型显式描述模型, 采用基于水平集函数的形状插值技术实现了点阵单胞的梯度构型生成, 同时所有梯度点阵均具有相似的拓扑特征, 保证了彼此之间构型的连接性.
(3)构建了基于Kriging的梯度点阵单胞宏观等效力学属性高效预测模型, 避免了均匀化方法计算所有梯度点阵单胞宏观等效力学属性的高昂计算成本, 有效提高了所提方法的计算效率.
(4)本文以体心立方点阵单胞验证了所提方法的有效性, 但所提方法同样适用于其他点阵单胞构型; 同时本文只研究了点阵结构的高刚度承载性能, 因此, 可在本文研究基础上进一步考虑点阵结构减振隔振、冲击性能、抗弯曲等力学性能, 从而实现点阵结构的超轻质、高性能和多功能设计.
-
表 1 设计域30 cm × 15 cm结构图及柔度
Table 1 The structure and compliance of design domain 30 cm × 15 cm
Uniform Linear graded Topology graded Relative density structure compliance 124.24 N·cm 88.84 N·cm 60.55 N·cm 表 2 设计域15 cm × 15 cm结构图及柔度
Table 2 The structure and compliance of design domain 15 cm × 15 cm
Uniform Linear graded Topology graded Relative density structure compliance 293.04 N·cm 229.45 N·cm 172.34 N·cm 表 3 设计域40 cm × 15 cm × 10 cm结构图及柔度
Table 3 The structure and compliance of design domain 40 cm × 15 cm × 10 cm
Uniform Linear graded Topology graded Relative density structure compliance 10171.05 N·cm 13759.11 N·cm 4979.70 N·cm 表 4 设计域15 cm × 20 cm × 10 cm结构图与柔度
Table 4 The structure and compliance of design domain 15 cm × 20 cm × 10 cm
Uniform Linear graded Topology graded Relative density structures compliance 8165.04 N·cm 7708.78 N·cm 6945.11 N·cm -
[1] Jia Z, Liu F, Jiang X, et al. Engineering lattice metamaterials for extreme property, programmability, and multifunctionality. Journal of Applied Physics, 2020, 127(15): 150901 doi: 10.1063/5.0004724
[2] Jihong Z, Han Z, Chuang W, et al. A review of topology optimization for additive manufacturing: Status and challenges. Chinese Journal of Aeronautics, 2020, 34(1): 91-110
[3] 易长炎, 柏龙, 陈晓红等. 金属三维点阵结构拓扑构型研究及应用现状综述. 功能材料, 2017, 48(10): 10055-10065 (Yi Changyan, Bai long, Chen Xiaohong, et al. Review on the metal three-dimensional lattice topology configurations research and application status. Journal of Functional Materials, 2017, 48(10): 10055-10065 (in Chinese) [4] 雷红帅, 赵则昂, 郭晓岗等. 航天器轻量化多功能结构设计与制造技术研究进展. 宇航材料工艺, 2021, 51(4): 10-22 (Lei Hongshuai, Zhao Zegang, Guo Xiaogang, et al. Research progress on the design and manufacture technology of lightweight multifunctional spacecraft structures. Aerospace Materials & Technology, 2021, 51(4): 10-22 (in Chinese) doi: 10.12044/j.issn.1007-2330.2021.04.002 [5] 陶斯嘉, 王小锋, 曾婧等. 点阵材料及其3D打印. 中国有色金属学报, 2022, 32(2): 416-444 (Tao Sijia, Wang Xiaofeng, Zeng Jing, et al. Lattice materials and its fabrication by 3D printing: A review. The Chinese Journal of Nonferrous Metals, 2022, 32(2): 416-444 (in Chinese) [6] Shuheng W, Yongbin M, Zichen D. Two-node method for the effective elastic modulus of periodic cellular truss materials and experiment verification via stereolithography. European Journal of Mechanics - A/Solids, 2020, 87: 104201
[7] Jung A, Diebels S. Microstructural characterisation and experimental determination of a multiaxial yield surface for open-cell aluminium foams. Materials & Design, 2017, 131: 252-264
[8] Nazir A, Abate KM, Kumar A, et al. A state-of-the-art review on types, design, optimization, and additive manufacturing of cellular structures. The International Journal of Advanced Manufacturing Technology, 2019, 104(9): 3489-3510
[9] Bai L, Yi C, Chen X, et al. Effective design of the graded strut of bcc lattice structure for improving mechanical properties. Materials, 2019, 12(13): 2192 doi: 10.3390/ma12132192
[10] El-Sayed MA, Essa K, Ghazy M, et al. Design optimization of additively manufactured titanium lattice structures for biomedical implants. The International Journal of Advanced Manufacturing Technology, 2020, 110(9): 2257-2268
[11] 王书恒, 戴时, 吴鑫伟等. 考虑材料各向异性的熔丝制造PLA点阵结构弹性各向同性设计. 力学学报, 2022, 54(5): 1291-1302 (Wang Shuheng, Dai Shi, Wu Xinwei, et al. Design of elastically isotropic PLA lattice strucrure in fused filament fabrication considering material anisotropy. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(5): 1291-1302 (in Chinese) doi: 10.6052/0459-1879-22-031 [12] 徐世鹏, 丁晓红, 段朋云等. 考虑时变刚度特性的复合材料微结构拓扑优化设计方法. 力学学报, 2022, 54(1): 134-146 (Xu Shipeng, Ding Xiaohong, Duan Pengyun, et al. Topology optimization of composite material microstructure considering time-changeable stiffness. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(1): 134-146 (in Chinese) doi: 10.6052/0459-1879-21-395 [13] Seharing A, Azman AH, Abdullah S. A review on integration of lightweight gradient lattice structures in additive manufacturing parts. Advances in Mechanical Engineering, 2020, 12(6): 1-21
[14] Yu S, Sun J, Bai J. Investigation of functionally graded TPMS structures fabricated by additive manufacturing. Materials & Design, 2019, 182: 108021
[15] Peng Z, Dexing Q, Rui X, et al. Mechanical design and energy absorption performances of rational gradient lattice metamaterials. Composite Structures, 2021, 277: 114606 doi: 10.1016/j.compstruct.2021.114606
[16] Dumas M, Terriault P, Brailovski V. Modelling and characterization of a porosity graded lattice structure for additively manufactured biomaterials. Materials & Design, 2017, 121: 383-392
[17] Liu F, Mao Z, Zhang P, et al. Functionally graded porous scaffolds in multiple patterns: New design method, physical and mechanical properties. Materials & Design, 2018, 160: 849-860
[18] Sanjairaj V, Zhang L, Zhang S, et al. Triply periodic minimal surfaces sheet scaffolds for tissue engineering applications: An optimization approach toward biomimetic scaffold design. ACS Applied Bio Materials, 2018, 1(2): 259-269 doi: 10.1021/acsabm.8b00052
[19] Bai L, Gong C, Chen X, et al. Mechanical properties and energy absorption capabilities of functionally graded lattice structures: Experiments and simulations. International Journal of Mechanical Sciences, 2020, 182: 105735 doi: 10.1016/j.ijmecsci.2020.105735
[20] Chamini R, Shanqing X, Yvonne D, et al. Crushing behavior of functionally graded lattice. JOM, 2021, 73(12): 4130-4140 doi: 10.1007/s11837-021-04946-x
[21] Li H, Luo Z, Gao L, et al. Topology optimization for functionally graded cellular composites with metamaterials by level sets. Computer Methods in Applied Mechanics and Engineering, 2018, 328: 340-364 doi: 10.1016/j.cma.2017.09.008
[22] Sigmund O, Maute K. Topology optimization approaches. Structural and Multidisciplinary Optimization, 2013, 48(6): 1031-1055 doi: 10.1007/s00158-013-0978-6
[23] 廖中源, 王英俊, 王书亭. 基于拓扑优化的变密度点阵结构体优化设计方法. 机械工程学报, 2019, 55(8): 65-72 (Liao Zhongyuan, Wang Yingjun, Wang Shuting, et al. Graded-density lattice structures optimization design based on topology optimization. Journal of Mechanical Engineering, 2019, 55(8): 65-72 (in Chinese) doi: 10.3901/JME.2019.08.065 [24] 蔡金虎, 王春洁. 基于映射的梯度点阵结构设计方法. 振动与冲击, 2020, 39(20): 74-81 (Cai Jinhu, Wang Chunjie. A graded lattice structures design method based on mapping progress. Journal of Vibration and Shock, 2020, 39(20): 74-81 (in Chinese) doi: 10.13465/j.cnki.jvs.2020.20.010 [25] 赵芳垒, 敬石开, 刘晨燕. 基于局部相对密度映射的变密度多孔结构设计方法. 机械工程学报, 2018, 54(19): 121-128 (Zhao Fanglei, Jing Shikai, Liu Chenyan. Variable density cellular structure design method base on local relative density mapping. Journal of Mechanical Engineering, 2018, 54(19): 121-128 (in Chinese) doi: 10.3901/JME.2018.19.121 [26] 侯淑娟, 梁慧妍, 汪全中等. 基于迭代法的非线性弹性均质化研究. 力学学报, 2018, 50(4): 837-846 (Hou Shujuan, Liang Huiyan, Wang Quanzhong, et al. Study on nonlinear elastic homogenization with iterative method. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(4): 837-846 (in Chinese) doi: 10.6052/0459-1879-18-039 [27] Zadpoor AA. Mechanical performance of additively manufactured meta-biomaterials. Acta Biomaterialia, 2018, 85: 41-59
[28] Lei Y, Massimiliano F, Raya M, et al. An investigation into the effect of gradients on the manufacturing fidelity of triply periodic minimal surface structures with graded density fabricated by selective laser melting. Journal of Materials Processing Tech, 2019, 275: 116367
[29] Chu S, Gao L, Xiao M, et al. Design of sandwich panels with truss cores using explicit topology optimization. Composite Structures, 2018, 210: 892-905
[30] 付君健, 舒正涛, 田启华等. 功能梯度多孔结构拓扑优化的混合水平集方法. 机械工程学报, 2022, 48: 1-12 (Fu Junjian, Shu Zhengtao, Tian Qihua, et al. A hybrid level set method for topology optimization of functionally graded cellular structures. Journal of Mechanical Engineering, 2022, 48: 1-12 (in Chinese) [31] 郭旭, 赵康. 基于拓扑描述函数的连续体结构拓扑优化方法. 力学学报, 2004, 36(5): 520-526 (Guo Xu, Zhao Kang. A new topology description function based approach for structural topology optimization. Chinese Journal of Theoretical and Applied Mechanics, 2004, 36(5): 520-526 (in Chinese) doi: 10.3321/j.issn:0459-1879.2004.05.002 [32] 赵丹阳, 刘韬, 李红霞等. 可降解聚合物血管支架结构优化设计. 力学学报, 2017, 49(6): 1409-1417 (Zhao Danyang, Liu Tao, Li Hongxia, et al. Optimization design of degraable polymer vascular stent structure. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(6): 1409-1417 (in Chinese) doi: 10.6052/0459-1879-17-214 [33] Zhang Y, Li H, Xiao M, et al. Concurrent topology optimization for cellular structures with nonuniform microstructures based on the kriging metamodel. Structural and Multidisciplinary Optimization, 2019, 59(4): 1273-1299 doi: 10.1007/s00158-018-2130-0
[34] Zhang Y, Zhang L, Ding Z, et al. A multiscale topological design method of geometrically asymmetric porous sandwich structures for minimizing dynamic compliance. Materials & Design, 2022, 214: 110404
[35] Mi X, Xiliang L, Yan Z, et al. Design of graded lattice sandwich structures by multiscale topology optimization. Computer Methods in Applied Mechanics and Engineering, 2021, 384: 113949 doi: 10.1016/j.cma.2021.113949
[36] Xiliang L, Liang G, Mi X, et al. Kriging-assisted design of functionally graded cellular structures with smoothly-varying lattice unit cells. Computer Methods in Applied Mechanics and Engineering, 2022, 390: 114466 doi: 10.1016/j.cma.2021.114466