﻿ 考虑结构自重的基于NURBS插值的3D拓扑描述函数法<sup>1)</sup>
 文章快速检索 高级检索
 0

### 引用本文 [复制中英文]

[复制中文]
Wang Xuan, Hu Ping, Zhu Xuefeng, Gai Yundong. TOPOLOGY DESCRIPTION FUNCTION APPROACH USING NURBS INTERPOLATION FOR 3D STRUCTURES WITH SELF-WEIGHT LOADS1)[J]. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(6): 1437-1445. DOI: 10.6052/0459-1879-16-145.
[复制英文]

### 基金项目

1) 国家自然科学基金（11272075，11302041）和中央高校基本科研业务费专项（DUT15RC（4）36）资助项目

### 通讯作者

2) 胡平, 教授, 主要研究方向:固体力学和车辆工程.E-mail:pinghu@dlut.edu.cn

### 文章历史

2016-05-27 收稿
2016-08-04 录用
2016-08-08网络版发表

*. 大连理工大学汽车工程学院, 大连 116024
†. 大连理工大学工业装备结构分析国家重点实验室, 大连 116024

1 基于NURBS插值的拓扑描述函数法 1.1 NURBS实体

 $V(\xi, \eta, \zeta ) = \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^m {\sum\limits_{k = 1}^l {P_{i, j, k} R_{i, j, k} } } } (\xi, \eta, \zeta )$ (1)

 $R_{i, j, k}^{p, q, r} (\xi, \eta, \zeta ) =\\ \quad \dfrac{N_{i, p} (\xi )M_{j, q} (\eta )L_{k, r} (\zeta )\omega _{i, j, k} }{\sum\limits_{\bar {i} = 1}^n {\sum\limits_{\bar {j} = 1}^m {\sum\limits_{\bar {k} = 1}^l {N_{\bar {i}, p} (\xi )M_{\bar {j}, q} (\eta )L_{\bar {k}, r} (\zeta )\omega _{\bar {i}, \bar {j}, \bar {k}} } } } }$ (2)

1.2 基于NURBS插值的拓扑描述函数法

 $\left.\begin{array}{ll} \phi (x) = 0, & \partial \varOmega \\ \phi (x) > 0, & \hbox{在} \ \varOmega \hbox{内} \\ \phi (x) < 0, & \hbox{在} \ D / \varOmega \hbox{内} \end{array} \!\!\right\}$ (3)

 $H(\phi (x)) = \left\{ \!\!\begin{array}{ll} 0, & \phi (x) < 0 \\ 1, & \phi (x) \geqslant 0 \end{array} \right.$ (4)

 $\phi ^h (x) = \sum\limits_{I = 1}^{n_{\rm cp}} R_I (\xi, \eta, \zeta )\phi _I$ (5)

 $\left.\begin{array}{l} E(x) = \Big [\sum\limits_{I = 1}^{n_{\rm cp}} R_I (\xi, \eta, \zeta )H(\phi _I ) \Big]^{g_1 }E_0 \\ \rho (x) = \Big [\sum\limits_{I = 1}^{n_{\rm cp}} R_I (\xi, \eta, \zeta )H(\phi _I ) \Big]^{g_2 }\rho _0 \end{array}\!\!\right\}$ (6)

 $\left. \begin{array}{l} \mathit{\boldsymbol{x}}(\xi, \eta, \zeta ) = \sum\limits_{I = 1}^{{n_{{\rm{cp}}}}} {{R_I}} (\xi, \eta, \zeta ){\mathit{\boldsymbol{\overline x}} _I} = \mathit{\boldsymbol{RX}}\\ \mathit{\boldsymbol{u}}(\xi, \eta, \zeta ) = \sum\limits_{I = 1}^{{n_{{\rm{cp}}}}} {{R_I}} (\xi, \eta, \zeta ){\mathit{\boldsymbol{\overline u}} _I} = \mathit{\boldsymbol{RU}} \end{array} \right\}$ (7)

1.3 Heaviside函数的光滑性处理

 $H_\varepsilon (\phi ) = \left\{ \begin{array}{ll} \delta, & {\phi \leqslant - \varepsilon } \\ {[1 + \sin (\pi \phi / 2\varepsilon )] / 2}, & { - \varepsilon < \phi < \varepsilon } \\ 1, & {\phi \geqslant \varepsilon } \end{array} \right.$ (10)

 ${\boldsymbol F}_e = \bar {\boldsymbol F}_e \int_{V_e } {\rho (x)g} {\rm{d}} V_e$ (13)

 ${\boldsymbol F}_e^i = \left\{ \begin{array}{cl} 0, & {i = 2 j - 1} \\ - \dfrac{1}{n_{\rm cp}}, & {i = 2 j} \end{array} \right.$ (14)

 ${\boldsymbol F}_e^i = \left\{\begin{array}{cl} 0 & {i = 3 j - 1, \ \ 3 j - 2} \\ - \dfrac{1}{n_{\rm cp}} & {i = 3j} \end{array} \right.$ (15)

 $V_e = \int_{V_e } \sum\limits_{I = 1}^{n_{\rm cp}} R_I H_\varepsilon (\phi _I ) {\rm{d}} V_e$ (16)
2.2 灵敏度分析

 $\left. \begin{array}{l} \frac{{\partial E(x)}}{{\partial {\phi _I}}} = {g_1}{[\sum\limits_{I = 1}^{{n_{cp}}} {{R_I}{H_\varepsilon }({\phi _I})}]^{{g_1} - 1}}{R_I}{H_{\varepsilon '}}({\phi _I}){E_0}\\ \frac{{\partial \rho (x)}}{{\partial {\phi _I}}} = {g_2}{[\sum\limits_{I = 1}^{{n_{cp}}} {{R_I}{H_\varepsilon }({\phi _I})}]^{{g_2} -1}}{R_I}{H_{\varepsilon '}}({\phi _I}){\rho _0} \end{array} \right\}$ (17)

 $\left. \begin{array}{l} \frac{{\partial {\mathit{\boldsymbol{K}}_e}}}{{\partial {\phi _I}}} = \int_{{V_e}} {\frac{{\partial E(x)}}{{\partial {\varphi _I}}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} {\mathit{\boldsymbol{C}}_0}\mathit{\boldsymbol{B}}{\rm{d}}{V_e}\\ \frac{{\partial {{\bf{F}}_e}}}{{\partial {\phi _I}}} = {{\bar F}_e}\int_{{V_e}} {\frac{{\partial \rho (x)}}{{\partial {\phi _I}}}g} {\rm{d}}\\ {V_e}\frac{{\partial {V_e}}}{{\partial {\phi _I}}} = \int_{{V_e}} {{R_I}{H_{\varepsilon '}}({\phi _I})} {\rm{d}}V \end{array} \right\}$ (18)

 $\dfrac{\partial c}{\partial \phi _I } = 2\dfrac{\partial {\boldsymbol F}^{\rm T}}{\partial \phi _I }{\boldsymbol U}-{\boldsymbol U}^{\rm T}\dfrac{\partial {\boldsymbol K}}{\partial \phi_I }{\boldsymbol U}$ (19)

 $\dfrac{\partial c}{\partial \phi _I } = 2\sum\limits_{Ie = 1}^{N_{Ie} } {\dfrac{\partial {\boldsymbol F}_{Ie}^{\rm T} }{\partial \phi _I }{\boldsymbol U}_{Ie} } - \sum\limits_{Ie = 1}^{N_{Ie} } {U_{Ie}^{\rm T} \dfrac{\partial {\boldsymbol K}_{Ie} }{\partial \phi _I } {\boldsymbol U}_{Ie} }$ (20)

 $\dfrac{\partial V}{\partial \phi _I } = \sum\limits_{Ie = 1}^{N_{Ie} } {\dfrac{\partial V_{Ie} }{\partial \phi _I }}$ (21)

2.3 优化算法

(1)构建设计域：定义设计域几何参数，使用NURBS构建优化模型的几何区域.

(2)定义模型参数：定义实体材料杨氏模量和质量密度及材料容许的体积分数等模型参数.

(3)等几何分析：运用等几何分析方法求解位移场.

(4)灵敏度分析：计算目标函数值，及通过式(20)和式(21)分别计算目标函数和材料体积对设计变量的敏度值.

(5)优化模型求解：用MMA算法求解优化模型，进行设计变量更新.

(6)收敛性判断：本文以最大的迭代次数来判断优化过程是否停止；若收敛，则迭代停止，输出优化结果；若不收敛，则返回步骤(3)，继续迭代直到收敛为止.

3 数值算例与讨论

3.1 固定外载荷的最小柔顺性问题：L形梁

 图 1 L形梁设计域和边界条件 Figure 1 Design domain and boundary condition for L-shaped beam with curvature

 图 2 L形梁的最优拓扑构型 Figure 2 Optimal topology of L-shaped beam

 图 3 柔顺性和拓扑构型的迭代历程 Figure 3 Iteration histories of the compliance and topology
3.2 2D自重拱设计

 图 4 算例2的模型和优化结果 Figure 4 Model and optimal topology of example 2
 图 5 算例2柔顺性和体积分数的进化历史 Figure 5 Evolution histories of the compliance and volume fraction of example 2

 图 6 约束体积分数与实际体积分数关系图 Figure 6 Relation between the constrained volume fraction and the real volume fraction
3.3 3D自重拱设计

 图 7 算例3的设计域和边界条件 Figure 7 Design domain and boundary condition of example 3

 图 8 3D自重拱的最优拓扑 Figure 8 Optimal topology of 3D arch

 图 9 3D拱形桥的最优拓扑 Figure 9 Optimal topology of 3D arched bridge

 图 10 不同边界条件下柔顺性和体积分数的迭代历程 Figure 10 Iteration histories of the compliance and volume fraction for different boundary conditions
4 结论

(1)高阶NURBS基函数被同时用于三维NURBS实体片中的几何场、位移场及设计变量场插值，实现了几何模型、分析模型和优化模型的有效统一，确保了位移场及设计变量场的高阶连续性；

(2)无需对刚度和质量密度的插值格式进行修正，本文方法能够有效克服考虑自重载荷作用的优化问题中出现的材料的寄生效应及目标函数的非单调性问题；

(3)研究表明，本文方法优化迭代过程稳健，收敛性快.