EI、Scopus 收录
中文核心期刊

海洋柔性立管防弯器结构非线性刚度拓扑优化

范志瑞, 杨志勋, 许琦, 苏琦, 牛斌, 赵国忠

范志瑞, 杨志勋, 许琦, 苏琦, 牛斌, 赵国忠. 海洋柔性立管防弯器结构非线性刚度拓扑优化. 力学学报, 2022, 54(4): 929-938. DOI: 10.6052/0459-1879-21-589
引用本文: 范志瑞, 杨志勋, 许琦, 苏琦, 牛斌, 赵国忠. 海洋柔性立管防弯器结构非线性刚度拓扑优化. 力学学报, 2022, 54(4): 929-938. DOI: 10.6052/0459-1879-21-589
Fan Zhirui, Yang Zhixun, Xu Qi, Su Qi, Niu Bin, Zhao Guozhong. Nonlinear stiffness topology optimization for the bend stiffener of flexible riser. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(4): 929-938. DOI: 10.6052/0459-1879-21-589
Citation: Fan Zhirui, Yang Zhixun, Xu Qi, Su Qi, Niu Bin, Zhao Guozhong. Nonlinear stiffness topology optimization for the bend stiffener of flexible riser. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(4): 929-938. DOI: 10.6052/0459-1879-21-589
范志瑞, 杨志勋, 许琦, 苏琦, 牛斌, 赵国忠. 海洋柔性立管防弯器结构非线性刚度拓扑优化. 力学学报, 2022, 54(4): 929-938. CSTR: 32045.14.0459-1879-21-589
引用本文: 范志瑞, 杨志勋, 许琦, 苏琦, 牛斌, 赵国忠. 海洋柔性立管防弯器结构非线性刚度拓扑优化. 力学学报, 2022, 54(4): 929-938. CSTR: 32045.14.0459-1879-21-589
Fan Zhirui, Yang Zhixun, Xu Qi, Su Qi, Niu Bin, Zhao Guozhong. Nonlinear stiffness topology optimization for the bend stiffener of flexible riser. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(4): 929-938. CSTR: 32045.14.0459-1879-21-589
Citation: Fan Zhirui, Yang Zhixun, Xu Qi, Su Qi, Niu Bin, Zhao Guozhong. Nonlinear stiffness topology optimization for the bend stiffener of flexible riser. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(4): 929-938. CSTR: 32045.14.0459-1879-21-589

海洋柔性立管防弯器结构非线性刚度拓扑优化

基金项目: 国家自然科学基金(U1906233, 52001088, 11732004和51975087), 山东省重点研发计划(2019 JZZY010801), 黑龙江省自然科学基金(LH2021 E050)和工业装备结构分析国家重点实验室开发基金(GZ20105)资助项目
详细信息
    作者简介:

    杨志勋, 副教授, 主要研究方向: 海洋工程柔性管缆结构多尺度分析及实验验证. E-mail: yangzhixun@hrbeu.edu.cn

  • 中图分类号: TH122

NONLINEAR STIFFNESS TOPOLOGY OPTIMIZATION FOR THE BEND STIFFENER OF FLEXIBLE RISER

  • 摘要: 防弯器是海洋柔性立管过弯保护的关键附件, 对提高立管的结构安全性具有重要意义. 目前, 防弯器结构设计主要采用尺寸优化方法. 然而, 与拓扑优化方法相比, 该方法能提供的设计空间有限, 其在提高防弯器的力学性能, 以及探索防弯器创新构型方面具有很大不足. 本文在Dirichlet边界条件下, 以最大化弯曲刚度为目标, 对同时考虑材料和几何非线性的防弯器结构拓扑优化方法进行研究. 通过引入Helmholtz-PDE过滤和Heaviside惩罚, 以克服优化中出现的棋盘格现象和灰度单元等数值不稳定性问题. 与此同时, 研究引入了对称算子, 以提高往复性载荷作用下防弯器结构的承载能力和可制造性, 并且采用伴随法对优化问题的灵敏度进行了推导. 此外, 为了提高结构分析和优化的效率, 研究还基于PETSc库建立了并行程序框架. 数值算例中, 在材料体分比相同的情况下, 对防弯器结构分别进行了2D和3D非线性拓扑优化, 并对两种优化结果的承载能力进行了对比. 数值算例结果表明, 相比于防弯器2D拓扑优化结果, 在大部分波浪载荷方向下, 3D拓扑优化所给出的防弯器设计方案具有更为优越的结构性能. 本文所建立的3D非线性拓扑优化技术, 为深水恶劣海况下的高性能防弯器结构设计提供了新的理论模型和实现技术.
    Abstract: The bend stiffener, as the key over-bending protection accessory, is of significant importance to improve the safety of the flexible riser used in deep water. At present, the size optimization method is mainly used in the structural design of the bend stiffener. However, compared with the topology optimization, the design freedom provided by this method is limited, and it lacks capacities in sufficiently improving the mechanical performance and finding the novel configurations of the bend stiffener. In the present study, under Dirichlet boundary conditions, a topology optimization method considering the material and geometric nonlinearity is developed to maximize the structural stiffness of the bend stiffener. The Helmholtz-PDE filter and Heaviside projection are introduced to eliminate the numerical issues caused by the checkerboard pattern and the gray element phenomenon, respectively. The symmetry operator is employed to enhance the load bearing capability under the reciprocating ocean wave load and improve the manufacturability of the bend stiffener. Making use of the adjoint method, a sensitivity analysis is performed to enable a gradient-based algorithm for solving the optimization problems. Simultaneously, a parallel computational framework based on PETSc library is also utilized to improve the efficiency of the structural analysis and optimization. In the numerical examples, with the constant material volume fraction, 2D and 3D optimizations for the bend stiffener are performed to improve the stiffness of the bend stiffener, respectively. Based on that, the load carrying capacity of the two optimization results under different load directions are compared. The numerical examples show that, compared to the 2D optimized result, the 3D optimization can significantly improve the stiffness of the bend stiffener in most loading directions. The present 3D nonlinear topology optimization method provides the new theorical model and implementation technology for the high-performance bend stiffener with the severe water environment in the deep ocean.
  • 海洋柔性立管是连接海上浮体及海底井口的关键装备, 被誉为海洋浮式生产系统的“血管”[1-2]. 如图1所示, 在服役过程中, 由于风浪流及浮体运动的作用, 立管与浮体间的连接部位会因弯曲变形过大而造成结构的失效破坏. 因此, 防弯器在防止立管的过弯失效, 以及提高立管结构安全性等方面具有重要的意义[3].

    图  1  柔性立管与防弯器结构示意图
    Figure  1.  Schematic diagram about the flexible riser and bend stiffener

    国内外学者对防弯器结构的分析方法做了大量研究, 并提出了不同的结构分析模型. 如Deruntz[4]所建立的细长梁分析模型, Boet和Out[5]提出的变截面梁理论模型, Tong等[6]提出的斜拉力悬臂梁模型, 以及Drobyshevski [7]建立的异形梁模型等. 上述研究中多采用梁模型对防弯器结构开展性能分析, 可有效降低防弯器结构建模和分析的成本. 然而, 梁结构的分析精度依赖于实际防弯器结构的长细比, 特别对于深水应用的防弯器, 当长细比较小时, 采用梁结构可能会造成较大的分析误差. 为了能够更加准确地获得防弯器性能, 本文研究中将采用能够反映防弯器截面/整体真实构型的平面/实体有限元模型.

    结构优化是提高防弯器结构性能的有效设计方法. Tanaka等[8]以柔性立管的最大曲率, 以及防弯器尺寸、最大应变和端部最大弯矩为约束, 采用遗传算法对防弯器进行了轻量化设计. Yang和Wang[9]采用最优拉丁超立方方法建立了考虑材料、几何和载荷不确定性的柔性立管疲劳分析代理模型, 通过优化提高了防弯器结构的疲劳可靠性. Tang等[10]以质量和疲劳强度设置为目标, 以柔性立管的最大曲率为约束, 对防弯器的结构尺寸参数进行了多目标优化. 然而, 上述基于尺寸优化的结构设计所提供的设计自由度有限, 难以充分发挥防弯器结构的性能. 拓扑优化方法[11-14]能够在给定的目标函数和约束条件下, 获得设计域内最优的材料分布, 发现结构的创新构型, 从而相比传统的尺寸优化在更大程度上提高结构的性能. 因此, 本文研究中采用拓扑优化方法对防弯器进行结构优化.

    实际的海洋工程中, 防弯器由合成树脂材料构造而成, 同时在位运行中, 由于恶劣的海洋环境载荷往往发生几何大变形. 因此在防弯器结构分析和优化设计中需要考虑几何和材料的非线性影响[10]. 与线弹性结构的拓扑优化不同, 在结构的非线性分析中, 外载荷与位移响应之间表现为非线性映射关系. 这使得结构刚度的描述既可以采用切线刚度阵, 也可以采用割线刚度阵, 且两者所描述的结构刚度性能往往具有较大差异, 如图2所示. Wallin等[15]的研究表明以结构割线刚度最大化和切线刚度最大化为目标时, 优化结果具有很大的不同. Kemmler等[16]对非线性刚度优化中的目标函数进行了总结, 其包括割线刚度(end-compliance)[15,17]、切线刚度(end-stiffness) [15]、应变能[18]以及余能[19]等. 实际工程中应根据实际需要选取恰当的目标函数.

    图  2  切线刚度${{\boldsymbol{K}}_{\tan }}$与割线刚度${{\boldsymbol{K}}_{\cot }}$的定义, 其中${\boldsymbol{F}}$${\boldsymbol{u}}$分别为结构的载荷和位移
    Figure  2.  The definition of tangent stiffness ${{\boldsymbol{K}}_{\tan }}$ and secant stiffness ${{\boldsymbol{K}}_{\cot }}$, where ${\boldsymbol{F}}$ and ${\boldsymbol{u}}$ are the force and displacement of a structure, respectively

    综上所述, 本文在考虑材料和几何非线性情况下, 以结构刚度最大化为目标, 对柔性立管的防弯器结构进行拓扑优化. 并在此基础上, 对不同设计假设下所得优化结果进行对比, 探寻不同载荷方向对防弯器承载能力的影响.

    图1中柔性立管及防弯器的工作过程进行抽象, 可得图3所示简化的力学模型. 其中, 模型的上端固定, 在立管下端部施加波浪载荷. 本文研究中涉及2D和3D两种防弯器设计方案, 为了避免两种方案中复杂的载荷转换关系, 采用了Dirichlet边界条件来模拟波浪载荷的作用, 即在$y - z$平面内, 立管端部施加位移${{\boldsymbol{a}}_{\text{p}}}$, 且${{\boldsymbol{a}}_{\text{p}}}$$y$轴之间的夹角记为$\theta $.

    图  3  海洋柔性立管、防弯器力学模型及边界条件
    Figure  3.  The mechanical model and boundary conditions of ocean flexible riser and the bend stiffener

    关于Dirichlet边界条件下的线弹性结构拓扑优化问题, Niu等[20]采用了以结构支反力${{\boldsymbol{R}}_{\text{p}}}$与节点位移${{\boldsymbol{a}}_{\text{p}}}$之间的点积进行定义的目标函数$c$, 即

    $$ c = {{\boldsymbol{R}}_{\text{p}}} \cdot {\text{ }}{{\boldsymbol{a}}_{\text{p}}} $$ (1)

    然而, 如前所述, 当结构呈现非线性特征时, 结构刚度的描述比较复杂. 以下将对上式的力学意义进行进一步分析.

    图4给出了不同结构刚度下支反力随位移变化示意图. 由图可知, 当位移加载结束时, 曲线2的支反力大于曲线1, 即${R_{{\text{p2}}}} > {R_{{\text{p1}}}}$. 同时, 曲线2所对应的结构割线刚度也大于曲线1的对应值, 即${\alpha _2} > {\alpha _1}$. 因此, 就割线刚度角度而言, 曲线2所对应的结构刚度大于曲线1. 与此同时, 由于位移${{\boldsymbol{a}}_{\text{p}}}$已经给定, 式(1)所示目标函数的值与具体的加载过程无关, 而仅与加载结束时的支反力水平相关. 此时, 式(1)与结构割线刚度所反映的结构性能一致, 即支反力越大, 结构的割线刚度越大.

    图  4  式(1)所定义目标函数的力学含义
    Figure  4.  Nature of the objective defined by Eq. (1)

    此外, 在图4中, 曲线2始终位于曲线1的上方, 因此曲线2与坐标轴围成的面积大于曲线1, 即曲线2所对应的结构应变能大于曲线1. 并且, 在加载结束点处, 曲线2的切线斜率大于曲线1. 但是, 上述情况并不表示式(1)所示目标函数可以反映结构的应变能和切线刚度水平. 这是因为应变能水平与整体的加载过程相关, 而切线刚度则与加载结束点附近的邻域相关. 如图5图6所示, 曲线2的割线刚度大于曲线1, 但是前者的应变能和切线刚度均小于后者.

    图  5  应变能与式(1)所定义目标函数的不同
    Figure  5.  The Difference between the strain energy and the objective defined by Eq. (1)
    图  6  切线刚度与式(1)所定义目标函数的不同
    Figure  6.  The difference between the tangent stiffness and the objective defined by Eq. (1)

    综上分析, 本研究中采用Dirichlet边界条件下割线刚度作为优化目标, 立管防弯器结构的拓扑优化列式可描述为

    $$ \left. \begin{array}{l} \begin{array}{*{20}{c}} \begin{gathered} {\text{find}} \hfill \\ \min \hfill \\ {\text{s}}{\text{.t}}{\text{.}} \hfill \\ {\text{ }} \hfill \\ {\text{ }} \hfill \\ \end{gathered} &\begin{gathered} {\boldsymbol{z}} = \left\{ {{z_1},{z_1}, \cdots ,{z_{\text{n}}}} \right\} \hfill \\ c{\text{(}}{\boldsymbol{a}}{\text{, }}{\boldsymbol{z}}{\text{) = }} - {{\boldsymbol{R}}_{\text{p}}} \cdot {\text{ }}{{\boldsymbol{a}}_{\text{p}}} \hfill \\ {\boldsymbol{R}}({\boldsymbol{a}},{\boldsymbol{z}}) = {{\boldsymbol{0}}} \hfill \\ g = \xi - {\xi _0} \leqslant 0 \hfill \\ 0 < {{{\varepsilon }}_0} < {z_k} \leqslant 1,\;\;k \in \left\{ {1,2, \cdots ,{{n}}} \right\} \hfill \\ \end{gathered} \end{array} \end{array} \right\}$$ (2)

    式中, ${\boldsymbol{z}}$为所有设计域内单元密度设计变量集合, ${\boldsymbol{a}}$为节点位移场, $ {\boldsymbol{R}}({\boldsymbol{a}},{\boldsymbol{z}}) $为不平衡载荷向量, ${{n}}$为结构内单元的总数, $g$为材料用量体分比约束, ${\xi _0}$为给定的材料用量体分比上限, $\xi = V/{V_0}$为优化后的材料体分比. 其中, ${V_0}$为结构的总体积, $V$为优化后材料所占的体积.

    弹性体的虚功平衡方程可表示为

    $$ \delta \psi \left( {{\boldsymbol{u}},\delta {\boldsymbol{u}}} \right) = \int_{{\varOmega _0}}^{} {\delta {\boldsymbol{E}}} :{\boldsymbol{S}}{\rm{d}}v - \int_{\partial {\varOmega _0}}^{} {\delta {\boldsymbol{u}} \cdot {\boldsymbol{t}}{\rm{d}}s} $$ (3)

    式中, ${\boldsymbol{u}}$为弹性体的位移场, ${\boldsymbol{E}}$为Green-Lagrange应变张量, ${\boldsymbol{S}}$为二阶Piola-Kirchhoff应力张量, ${\boldsymbol{t}}$为作用在弹性体外边界的牵引力.

    对式(3)进行泰勒展开, 并取一阶截断, 可得

    $$ \delta \psi ({\boldsymbol{u}} + {\rm{d}}{\boldsymbol{u}},\delta {\boldsymbol{u}}) = \delta \psi ({\boldsymbol{u}},\delta {\boldsymbol{u}}) + {\rm{d}}\delta \psi ({\boldsymbol{u}},\delta {\boldsymbol{u}}) $$ (4)

    与此同时, 引入Galerkin方法对上式进行有限元离散化, 即

    $$ \delta {\boldsymbol{u}} = {\boldsymbol{N}}\delta {\boldsymbol{a}} $$ (5)

    式中, ${\boldsymbol{N}}$为形函数矩阵, ${\boldsymbol{a}}$为单元位移场${\boldsymbol{u}}$的节点值.

    对式(4)右端第一项进行推导可得非线性迭代的平衡条件

    $$ {\boldsymbol{R}}({\boldsymbol{a}}) = {{\boldsymbol{F}}_{{{\rm{int}}} }} - {{\boldsymbol{F}}_{{\text{ext}}}} $$ (6)

    式中, ${{\boldsymbol{F}}_{{{\rm{int}}} }}$${{\boldsymbol{F}}_{{\text{ext}}}}$分别为内外力载荷向量, 其表达式如下

    $$ {{\boldsymbol{F}}_{{{\rm{int}}} }} = \int_{{\varOmega _0}}^{} {{{\boldsymbol{B}}^{\text{T}}}} {\boldsymbol{S}}{\rm{d}}v $$ (7)
    $$ {{\boldsymbol{F}}_{{\text{ext}}}} = \int_{\partial {\varOmega _0}}^{} {{{\boldsymbol{N}}^{\text{T}}}} {\boldsymbol{t}}{\rm{d}}s $$ (8)

    对式(4)右端第二项进行推导可得切线刚度阵的表达式为

    $$ {\boldsymbol{K}} = \int_\varOmega ^{} {({{\boldsymbol{B}}^{\text{T}}}{\boldsymbol{DB}} + {{\boldsymbol{H}}^{\text{T}}}{\boldsymbol{RH}}){\rm{d}}v} $$ (9)

    式中, ${\boldsymbol{D}}$为本构矩阵, ${\boldsymbol{R}}$为结构应力相关矩阵, ${\boldsymbol{B}}$${\boldsymbol{H}}$均为单元的应变−位移矩阵. 关于非线性分析的具体理论详见文献[21].

    在本文研究中采用Dirichlet边界条件来模拟海洋波浪载荷的作用. 如图7所示, 结构的边界$\partial \varOmega $上施加位移条件的边界包含两种, 即固定边界$\partial {\varOmega _{\text{z}}}$和非零位移边界$\partial {\varOmega _{\text{p}}}$, 并且$\partial {\varOmega _{\text{u}}} = \partial {\varOmega _{\text{z}}} \cup \partial {\varOmega _{\text{p}}}$, $\partial {\varOmega _{\text{z}}} \cap \partial {\varOmega _{\text{p}}} = \emptyset $. 除Dirichlet边界外的其它自由边界记作$\partial {\varOmega _{\text{f}}}$, 即$\partial {\varOmega _{\text{f}}} = \partial \varOmega - \partial {\varOmega _{\text{u}}}$. 为了后续推导方便, 在此约定下标“${\text{u}}$”, “${\text{f}}$”, “${\text{z}}$”和“${\text{p}}$”分别表示与边界$\partial {\varOmega _{\text{u}}}$, $\partial {\varOmega _{\text{f}}}$, $\partial {\varOmega _{\text{z}}}$$\partial {\varOmega _{\text{p}}}$相关的物理量.

    图  7  Dirichlet边界条件结构分析示意图
    Figure  7.  Illustration of Dirichlet boundary conditions

    Dirichlet边界条件作用下的Newton-Raphson迭代过程中, 边界上无任何外载荷作用, 即

    $$ {{\boldsymbol{F}}_{{\text{ext}}}} = {{\bf{0}}} $$ (10)

    严格意义上而言, Newton-Raphson迭代过程中内外力的平衡条件${\boldsymbol{R}} = {{\bf{0}}}$是指所有位于边界$\partial {\varOmega _{\text{f}}}$上的节点所对应的不平衡载荷向量${{\boldsymbol{R}}_{\text{f}}}$取值为0, 即

    $$ {{\boldsymbol{R}}_{\text{f}}} = {{\boldsymbol{F}}_{{{\rm{int}}} ,{\text{f}}}} - {{\boldsymbol{F}}_{{\text{ext,f}}}} = {{\bf{0}}\;\; }{\text{on }}\;\;\partial {\varOmega _{\text{f}}} $$ (11)

    而对于施加位移约束的边界$\partial {\varOmega _{\text{u}}}$, 内外力的差值则不为0. 该差值即为结构的支反力

    $$ {{\boldsymbol{R}}_{\text{u}}} = {{\boldsymbol{F}}_{{{\rm{int}}} ,{\text{u}}}} - {{\boldsymbol{F}}_{{\text{ext,u}}}}\;\;{\text{ on }}\;\;\partial {\varOmega _{\text{u}}} $$ (12)

    结合式(10), Dirichlet边界下的支反力可进一步表示为

    $$ {{\boldsymbol{R}}_{\text{u}}} = {{\boldsymbol{F}}_{{{\rm{int}}} ,{\text{u}}}}\;\;{\text{ on }}\;\;\partial {\varOmega _{\text{u}}} $$ (13)

    此时, 上述支反力又可进一步分解为位于边界$\partial {\varOmega _{\text{z}}}$上的${{\boldsymbol{R}}_{\text{z}}}$和边界$\partial {\varOmega _{\text{p}}}$上的${{\boldsymbol{R}}_{\text{p}}}$. 相应地, 不平衡载荷向量${\boldsymbol{R}}$ 可分解为以下3部分

    $$ {\boldsymbol{R}} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{R}}_{\text{f}}^{\text{T}}}&{{\boldsymbol{R}}_{\text{p}}^{\text{T}}}&{{\boldsymbol{R}}_{\text{z}}^{\text{T}}} \end{array}} \right]^{\text{T}}} $$ (14)

    $ {\boldsymbol{R}}({\boldsymbol{a}} + {\rm{d}}{\boldsymbol{a}}) $进行泰勒展开, 并取一阶截断可得

    $$ {\boldsymbol{R}}({\boldsymbol{a}} + {\rm{d}}{\boldsymbol{a}}) \approx {\boldsymbol{R}}({\boldsymbol{a}}) + \frac{{\partial {\boldsymbol{R}}({\boldsymbol{a}})}}{{\partial {\boldsymbol{a}}}}{\rm{d}}{\boldsymbol{a}} $$ (15)

    式中, $ \partial {\boldsymbol{R}}({\boldsymbol{a}})/\partial {\boldsymbol{a}} $为结构的切线刚度阵K. 且根据式(14), 其可进一步写为如下分块矩阵格式

    $$ \begin{split} {\boldsymbol{K}} =& \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}}}&{\dfrac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {{\boldsymbol{a}}_{\text{p}}}}}}&{\dfrac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {{\boldsymbol{a}}_{\text{z}}}}}} \\ {\dfrac{{\partial {{\boldsymbol{R}}_{\text{p}}}}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}}}&{\dfrac{{\partial {{\boldsymbol{R}}_{\text{p}}}}}{{\partial {{\boldsymbol{a}}_{\text{p}}}}}}&{\dfrac{{\partial {{\boldsymbol{R}}_{\text{p}}}}}{{\partial {{\boldsymbol{a}}_{\text{z}}}}}} \\ {\dfrac{{\partial {{\boldsymbol{R}}_{\text{z}}}}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}}}&{\dfrac{{\partial {{\boldsymbol{R}}_{\text{z}}}}}{{\partial {{\boldsymbol{a}}_{\text{p}}}}}}&{\dfrac{{\partial {{\boldsymbol{R}}_{\text{z}}}}}{{\partial {{\boldsymbol{a}}_{\text{z}}}}}} \end{array}} \right] = \\ &\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{K}}_{{\text{ff}}}}}&{{{\boldsymbol{K}}_{{\text{fp}}}}}&{{{\boldsymbol{K}}_{{\text{fz}}}}} \\ {{{\boldsymbol{K}}_{{\text{pf}}}}}&{{{\boldsymbol{K}}_{{\text{pp}}}}}&{{{\boldsymbol{K}}_{{\text{pz}}}}} \\ {{{\boldsymbol{K}}_{{\text{zf}}}}}&{{{\boldsymbol{K}}_{{\text{zp}}}}}&{{{\boldsymbol{K}}_{{\text{zz}}}}} \end{array}} \right] \end{split}$$ (16)

    研究采用Neo-Hookean超弹性材料来模拟防弯器和立管的材料非线性行为, 其应变能密度函数表达式为

    $$ {\omega _0} = \frac{{{K}}}{2}\left( {\frac{{{J^2} - 1}}{2} - \ln J} \right) + \frac{{{G}}}{2}\left( {{J^{ - \frac{2}{3}}}{\text{tr}}{\boldsymbol{C}} - 3} \right) $$ (17)

    式中, $J = \det ({\boldsymbol{F}})$, ${\boldsymbol{F}}$为应变梯度张量, ${\boldsymbol{C}} = {{\boldsymbol{F}}^{\text{T}}}{\boldsymbol{F}}$为右Cauchy应变张量, ${{K}}$为体模量, ${{G}}$为剪切模量.

    根据二阶Piola-Kirchhoff应力张量${\boldsymbol{S}}$和刚度张量${\boldsymbol{D}}$的定义可得两者的表达式分别为

    $$ {S_{ij}} = \frac{{{K}}}{2}\left( {J - 1} \right){C_{ij}} + {{G}}{J^{ - \frac{2}{3}}}\left( {{\delta _{ij}} - \frac{{{C_{pp}}}}{3}C_{ij}^{ - 1}} \right) $$ (18)
    $$ \qquad\;\;\begin{split} {D_{ijkl}} =& {a_1}C_{ij}^{ - 1}C_{kl}^{ - 1} - {a_2}\left( {{\delta _{ij}}C_{kl}^{ - 1} + C_{ij}^{ - 1}{\delta _{kl}}} \right) + \\ & {a_3}\left( {C_{ik}^{ - 1}C_{jl}^{ - 1} + C_{il}^{ - 1}C_{jk}^{ - 1}} \right) \end{split} $$ (19)

    式中

    $$ {a_1} = {{K}}{J^2} + \frac{{2{\text{G}}}}{9}{J^{ - \frac{2}{3}}}{C_{pp}} $$ (20)
    $$ {a_2} = \frac{{2{{G}}}}{3}{J^{ - \frac{2}{3}}} $$ (21)
    $$ {a_3} = \frac{{{G}}}{3}{J^{ - \frac{2}{3}}}{C_{pp}} - \frac{{{K}}}{2}\left( {{J^2} - 1} \right) $$ (22)

    在实际运行中, 防弯器结构往往受到往复性波浪载荷的作用. 为了保证防弯器具有较好的承载能力, 通常要求防弯器为对称结构. 为了保证优化后的结构的对称性, 研究中引入了对称算子${\boldsymbol{M}}$, 此时初始设计变量${\boldsymbol{z}}$与变换后的设计变量${\boldsymbol{\bar z}}$之间满足如下关系

    $$ {\boldsymbol{\bar z}} = {\boldsymbol{M}} \cdot {\boldsymbol{z}} $$ (23)

    以下引入3自由度系统来简要阐述对阵矩阵的构造方法. 初始的设计变量所构成的向量${\boldsymbol{z}}$

    $$ {\boldsymbol{z}} = \left\{ {\begin{array}{*{20}{c}} {{z_1}}&{{z_2}}&{{z_3}} \end{array}} \right\} $$ (24)

    同时要求设计变量具有如下对称关系

    $$ {z_1} = {z_3} $$ (25)

    则对称矩阵可以定义为如下两种方式

    $$ \left[ {\begin{array}{*{20}{c}} {{{\bar z}_1}} \\ {{{\bar z}_2}} \\ {{{\bar z}_3}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {0.5}&0&{0.5} \\ 0&{1.0}&0 \\ {0.5}&0&{0.5} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{z_1}} \\ {{z_2}} \\ {{z_3}} \end{array}} \right] $$ (26)
    $$ \left[ {\begin{array}{*{20}{c}} {{{\bar z}_1}} \\ {{{\bar z}_2}} \\ {{{\bar z}_3}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {1.0} \\ 0 \\ {1.0} \end{array}}&{\begin{array}{*{20}{c}} 0 \\ {1.0} \\ 0 \end{array}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{z_1}} \\ {{z_2}} \end{array}} \right] $$ (27)

    在式(26)中, 对称算子的构造涉及所有设计变量, 并且${\boldsymbol{M}} = {{\boldsymbol{M}}^{\text{T}}}$. 在式(27)中, 由于${z_3}$变量取值总与${z_1}$相等, 因此在构造对称算子${\boldsymbol{M}}$时, 可略去$ {z}_{3} $. 相应地, 此时对称算子${\boldsymbol{M}}$为非对称矩阵. 需注意的是, 在实际优化中, ${\boldsymbol{M}}$通常为稀疏矩阵. 并且, 对称矩阵更容易实施转置等运算. 因此, 研究中采用了式(26)所定义的对称算子.

    在拓扑优化过程中, 常会遇到棋盘格、网格依赖以及灰度单元等数值不稳定问题. 对于棋盘格和网格依赖性问题, 拓扑优化中常采用的方法包括灵敏度过滤[22-23]和密度过滤[24]等. 密度过滤可有效克服优化过程中的棋盘格现象. 在本文研究中, 为了提高结构分析和优化的效率, 基于MPI (massage passing interface)进程通信协议及PETSc库[25-26]建立了并行分析和优化框架. 然而, 在并行计算过程中, 单元信息及设计变量会在不同的进程中存储. 当采用传统的密度过滤时, 需要频繁的进程通信来获取临近的单元信息, 从而造成大量的进程通信成本.

    为了克服传统密度过滤方法在并行计算中的不足, 本研究采用了基于Helmholtz-PDE方程的过滤方法[27], 以下简称PDE过滤. PDE过滤的控制方程可写为如下格式

    $$ {{\boldsymbol{K}}_{\text{f}}}{\boldsymbol{\rho }} = {{\boldsymbol{T}}_{\text{f}}}{\boldsymbol{\bar z}} $$ (28)

    式中, ${\boldsymbol{\bar z}}$为过滤前所有单元密度构成的向量, ${\boldsymbol{\rho }}$为过滤后所有节点密度构成的向量. ${{\boldsymbol{K}}_{\text{f}}}$${{\boldsymbol{T}}_{\text{f}}}$的定义如下

    $$ {{\boldsymbol{K}}_{\text{f}}} = \int_\varOmega ^{} {\nabla {\boldsymbol{N}}_{\text{f}}^{\text{T}}{{\boldsymbol{K}}_{\text{d}}}\nabla {{\boldsymbol{N}}_{\text{f}}}{\rm{d}}\varOmega + \int_\varOmega ^{} {{\boldsymbol{N}}_{\text{f}}^{\text{T}}{{\boldsymbol{N}}_{\text{f}}}} } {\rm{d}}\varOmega $$ (29)
    $$ {{\boldsymbol{T}}_{\text{f}}} = \int_\varOmega ^{} {{\boldsymbol{N}}_{\text{f}}^{\text{T}}{\rm{d}}\varOmega } $$ (30)
    $$ {{\boldsymbol{K}}_{\text{d}}} = \sum\limits_{i = 1}^{{d}} {r_i^2{{\boldsymbol{v}}_i}{\boldsymbol{v}}_i^{\text{T}}} $$ (31)

    式中, ${{\boldsymbol{N}}_{\text{f}}}$为形函数, ${{\boldsymbol{K}}_{\text{d}}}$为过滤半径相关的正定矩阵, ${{d}}$为过滤的维度, ${r_i}$为维度$i$的过滤半径; ${{\boldsymbol{v}}_i}$为过滤的方向矢量. 当各个方向的过滤半径相同, 且${{\boldsymbol{v}}_i}$为单位向量时, 各向异性的过滤则退化为各向同性过滤. 过滤后单元内各高斯积分点处的单元密度可通过下式求得

    $$ \rho = {{\boldsymbol{N}}_{\text{f}}}{\boldsymbol{\rho }} $$ (32)

    研究中, 采用SIMP方法对单元的材料属性进行插值. 在传统线弹性拓扑优化中, SIMP方法的插值格式通常表述为对材料弹性模量或者单元刚度矩阵的插值. 在非线性优化问题中则略有区别, 通常将其表述为对应变能密度函数的插值, 即

    $$ \bar \omega = \bar \rho \omega $$ (33)
    $$ \bar \rho = {\rho ^p} $$ (34)

    式中, $\omega $$\bar \omega $分别为惩罚前后的应变能密度函数, $p$为惩罚指数, $\bar \rho $为惩罚后的单元密度. 由于单元内PDE过滤后的单元密度是变化的, 此处的$\rho $$\bar \rho $为高斯积分点处的单元密度. 需注意的是, 对于各向同性的Neo-Hookean超弹性材料, 上述插值格式与对弹性模量的插值格式等效.

    为了改善优化过程中的灰度单元现象. 在此, 引入了Heaviside惩罚策略[28-30], 并且Heaviside函数表达式为

    $$ \tilde \rho = H\left( {\bar \rho } \right) = \frac{{\tanh \left( {\beta \gamma } \right) + \tanh \left[ {\beta \left( {\bar \rho - \gamma } \right)} \right]}}{{\tanh \left( {\beta \gamma } \right) + \tanh \left[ {\beta \left( {1 - \gamma } \right)} \right]}} $$ (35)

    式中, $\tilde \rho $为Heaviside惩罚后的单元密度, $\beta $$\gamma $为控制参数. 当$\beta \to \infty $时, 上述Heaviside函数将趋于理想阶跃函数$h = \left( {\bar \rho - \gamma } \right)$, 如图8所示. 严格意义上而言, 为了使得过滤前后的体积守恒, $\gamma $的取值需要通过一维搜索进行确定. 一般地, 当r = 0.5时, 可使得过滤前后的体积近似保持, 并且可保证优化迭代过程的平稳性.

    图  8  Heaviside函数的性质
    Figure  8.  The characteristic of Heaviside function

    综合上述PDE过滤、SIMP插值以及Heaviside惩罚, 最终单元应变能密度的插值格式为

    $$ \tilde \rho = H\left\{ {\bar \rho \left\{ {\rho \left[ {{\boldsymbol{\rho }}\left( {\bar z} \right)} \right]} \right\}} \right\} $$ (36)

    采用伴随法对优化灵敏度进行推导, 则引入伴随向量${\boldsymbol{\lambda }}$后的目标函数可写为

    $$ c\left( {{\boldsymbol{a}},{\boldsymbol{z}}} \right) = c\left( {{\boldsymbol{a}},{\boldsymbol{z}}} \right) - {\boldsymbol{\lambda }}_{\text{f}}^{\text{T}}{{\boldsymbol{R}}_{\text{f}}}\left( {{\boldsymbol{a}},{\boldsymbol{z}}} \right) $$ (37)

    上式两端同时对单元密度设计变量${\boldsymbol{\bar z}}$求偏导数, 可得

    $$ \begin{gathered} \frac{{\partial c\left( {{\boldsymbol{a}},{\boldsymbol{z}}} \right)}}{{\partial {\boldsymbol{\bar z}}}} = \frac{{\partial c}}{{\partial {\boldsymbol{\bar z}}}} + \frac{{\partial c}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}}\frac{{\partial {{\boldsymbol{a}}_{\text{f}}}}}{{\partial {\boldsymbol{\bar z}}}} - {\boldsymbol{\lambda }}_{\text{f}}^{\text{T}}\left( {\frac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {\boldsymbol{\bar z}}}} + \frac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}}\frac{{\partial {{\boldsymbol{a}}_{\text{f}}}}}{{\partial {\boldsymbol{\bar z}}}}} \right) \\ \end{gathered} $$ (38)

    在此需注意的是, 上式右端中$\partial c/\partial {\boldsymbol{\bar z}}$以及$\partial c/\partial {{\boldsymbol{a}}_{\text{f}}}$分别表示目标函数对${\boldsymbol{\bar z}}$${{\boldsymbol{a}}_{\text{f}}}$的显式导数, 即在求导过程中首先假设${{\boldsymbol{a}}_{\text{f}}}$${\boldsymbol{\bar z}}$是独立的, 而${\boldsymbol{\bar z}}$${{\boldsymbol{a}}_{\text{f}}}$的影响则体现在$\partial {{\boldsymbol{a}}_{\text{f}}}/\partial {\boldsymbol{\bar z}}$项. 同理, 关于不平衡载荷向量${{\boldsymbol{R}}_{\text{f}}}$的求导过程同样存在上述情况.

    对式(38)进一步整理可得

    $$ \begin{gathered} \frac{{\partial c\left( {{\boldsymbol{a}},{\boldsymbol{z}}} \right)}}{{\partial {\boldsymbol{\bar z}}}} = \frac{{\partial c}}{{\partial {\boldsymbol{\bar z}}}} - {\boldsymbol{\lambda }}_{\text{f}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {\boldsymbol{\bar z}}}} + \left( {\frac{{\partial c}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}} - {\boldsymbol{\lambda }}_{\text{f}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}}} \right)\frac{{\partial {{\boldsymbol{a}}_{\text{f}}}}}{{\partial {\boldsymbol{\bar z}}}} \\ \end{gathered} $$ (39)

    同时, 结合式(1), 可得

    $$ \frac{{\partial c}}{{\partial {\boldsymbol{\bar z}}}} = - {\boldsymbol{a}}_{\text{p}}^{\text{T}}\frac{{\partial {\boldsymbol{R}}_{\text{p}}^{}}}{{\partial {\boldsymbol{\bar z}}}} - {\boldsymbol{R}}_{\text{p}}^{\text{T}}\frac{{\partial {\boldsymbol{a}}_{\text{p}}^{}}}{{\partial {\boldsymbol{\bar z}}}} = - {\boldsymbol{a}}_{\text{p}}^{\text{T}}\frac{{\partial {\boldsymbol{R}}_{\text{p}}^{}}}{{\partial {\boldsymbol{\bar z}}}} $$ (40)
    $$ \frac{{\partial c}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}} = - {\boldsymbol{a}}_{\text{p}}^{\text{T}}\frac{{\partial {\boldsymbol{R}}_{\text{p}}^{}}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}} - {\boldsymbol{R}}_{\text{p}}^{\text{T}}\frac{{\partial {\boldsymbol{a}}_{\text{p}}^{}}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}} = - {\boldsymbol{a}}_{\text{p}}^{\text{T}}\frac{{\partial {\boldsymbol{R}}_{\text{p}}^{}}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}} $$ (41)

    将式(40)和式(41)代入式(39)可得

    $$ \begin{gathered} \frac{{\partial c\left( {{\boldsymbol{a}},{\boldsymbol{z}}} \right)}}{{\partial {\boldsymbol{\bar z}}}} = - {\boldsymbol{a}}_{\text{p}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{p}}}}}{{\partial {\boldsymbol{\bar z}}}} - {\boldsymbol{\lambda }}_{\text{f}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {\boldsymbol{\bar z}}}} - \left( {{\boldsymbol{a}}_{\text{p}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{p}}}}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}} + {\boldsymbol{\lambda }}_{\text{f}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}}} \right)\frac{{\partial {{\boldsymbol{a}}_{\text{f}}}}}{{\partial {\boldsymbol{\bar z}}}} \\ \end{gathered} $$ (42)

    为了消除对隐导数$\partial {{\boldsymbol{a}}_{\text{f}}}/\partial {\boldsymbol{\bar z}}$的直接求解, 可令

    $$ {\boldsymbol{a}}_{\text{p}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{p}}}}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}} + {\boldsymbol{\lambda }}_{\text{f}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {{\boldsymbol{a}}_{\text{f}}}}} = {{\boldsymbol{0}}} $$ (43)

    同时, 联合式(16)可知, 伴随向量${\boldsymbol{\lambda }}$可通过下式求得

    $$ {\boldsymbol{K}}_{{\text{ff}}}^{\text{T}}{{\boldsymbol{\lambda }}_{\text{f}}} = - {\boldsymbol{K}}_{{\text{pf}}}^{\text{T}}{{\boldsymbol{a}}_{\text{p}}} $$ (44)

    此时, 目标函数对单元密度设计变量的导数可写为

    $$ \frac{{\partial c\left( {{\boldsymbol{a}},{\boldsymbol{z}}} \right)}}{{\partial {\boldsymbol{\bar z}}}} = - {\boldsymbol{a}}_{\text{p}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{p}}}}}{{\partial {\boldsymbol{\bar z}}}} - {\boldsymbol{\lambda }}_{\text{f}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {\boldsymbol{\bar z}}}} $$ (45)

    考虑PDE过滤时, 上式可进一步整理为

    $$ \frac{{\partial c\left( {{\boldsymbol{a}},{\boldsymbol{z}}} \right)}}{{\partial {\boldsymbol{\bar z}}}} = \left( { - {\boldsymbol{a}}_{\text{p}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{p}}}}}{{\partial {\boldsymbol{\rho }}}} - {\boldsymbol{\lambda }}_{\text{f}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {\boldsymbol{\rho }}}}} \right){\boldsymbol{K}}_{\text{f}}^{ - 1}{{\boldsymbol{T}}_{\text{f}}} $$ (46)

    为了灵敏度分析的求解, 对上式中的物理量进行如下变换

    $$ {\boldsymbol{a}}_{\text{p}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{p}}}}}{{\partial {\boldsymbol{\rho }}}} = \left[ {\begin{array}{*{20}{c}} {0}&{{\boldsymbol{a}}_{\text{p}}^{\text{T}}}&{0} \end{array}} \right]{\boldsymbol{\nu }}_{\text{a}}^{\text{T}} $$ (47)
    $$ {\boldsymbol{\lambda }}_{\text{f}}^{\text{T}}\frac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {\boldsymbol{\rho }}}} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{\lambda }}_{\text{f}}^{\text{T}}}&{0}&{0} \end{array}} \right]{\boldsymbol{\nu }}_{\text{a}}^{\text{T}} $$ (48)

    式中

    $$ {{\boldsymbol{\nu }}_{\text{a}}} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {{\boldsymbol{R}}_{\text{f}}}}}{{\partial {\boldsymbol{\rho }}}}}&{\dfrac{{\partial {{\boldsymbol{R}}_{\text{p}}}}}{{\partial {\boldsymbol{\rho }}}}}&{\dfrac{{\partial {{\boldsymbol{R}}_{\text{z}}}}}{{\partial {\boldsymbol{\rho }}}}} \end{array}} \right] = \dfrac{{\partial {\boldsymbol{R}}}}{{\partial {\boldsymbol{\rho }}}} $$ (49)

    此时, 式(39)可进一步整理为

    $$ \frac{{\partial c\left( {{\boldsymbol{a}},{\boldsymbol{z}}} \right)}}{{\partial {\boldsymbol{\bar z}}}} = - {{\boldsymbol{\kappa }}^{\text{T}}}\frac{{\partial {{\boldsymbol{F}}_{{{\rm{int}}} }}}}{{\partial {\boldsymbol{\rho }}}}{\boldsymbol{K}}_{\text{f}}^{ - 1}{{\boldsymbol{T}}_{\text{f}}} $$ (50)

    式中

    $$ {{\boldsymbol{\kappa }}^{\text{T}}} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{\lambda }}_{\text{f}}^{\text{T}}}&{{\boldsymbol{a}}_{\text{p}}^{\text{T}}}&{0} \end{array}} \right] $$ (51)

    根据式(7)中内力${{\boldsymbol{F}}_{{{\rm{int}}} }}$的表达式, 可得

    $$ {{\boldsymbol{\kappa }}^{\text{T}}}\frac{{\partial {{\boldsymbol{F}}_{{{\rm{int}}} }}}}{{\partial {\boldsymbol{\rho }}}} = \int_\varOmega ^{} {{\varsigma ^{\text{T}}}\frac{{\partial {\boldsymbol{S}}}}{{\partial {\boldsymbol{\rho }}}}{\rm{d}}\varOmega } $$ (52)

    式中, $ {\boldsymbol{\varsigma }} = {\boldsymbol{B\kappa }} $.

    进一步地, 考虑优化中的SIMP插值以及Heaviside惩罚, 上式可进一步写为如下离散格式

    $$ {{\boldsymbol{\kappa }}^{\text{T}}}\frac{{\partial {{\boldsymbol{F}}_{{{\rm{int}}} }}}}{{\partial {\boldsymbol{\rho }}}} = \bigcup\nolimits_{i = 1}^{\text{n}} {\int_{{\varOmega _i}}^{} {\left[ {\frac{p}{{{{\tilde \rho }_i}}}{H^{'}}\left( {{\rho _i}} \right){\boldsymbol{\varsigma }}_i^{\text{T}}{{\boldsymbol{S}}_i}} \right]{\boldsymbol{N}}_{\text{f}}^{\text{T}}{\rm{d}}{v_i}} } $$ (53)

    联合式(44)、式(50)和式(53)即可求得目标函数对设计变量${\boldsymbol{\bar z}}$的导数. 进一步地, 根据式(23)可得目标函数对${\boldsymbol{z}}$的导数为

    $$ \frac{{\partial c\left( {{\boldsymbol{a}},{\boldsymbol{z}}} \right)}}{{\partial {\boldsymbol{\bar z}}}} = {\boldsymbol{M}} \cdot \frac{{\partial c\left( {{\boldsymbol{a}},{\boldsymbol{z}}} \right)}}{{\partial {\boldsymbol{\bar z}}}} $$ (54)

    采用上述所建立非线性分析和优化理论对防弯器结构进行优化. 算例中, 所有参数均采用无量纲化处理. 防弯器的弹性模量${{{E}}_{\text{1}}} = 210$, 柔性立管的弹性模量为${{{E}}_2} = 2100$, 两种材料的泊松比均为$ {{\mu }} = 0.3 $. 在Newton-Raphson迭代中, 采用均匀的载荷增量步. 为了保证迭代的稳定性, 所有算例均经过预先试算来确定恰当载荷增量步数目. 一般而言, 载荷增量步数目的设置要保证每个增量步只需要3至4次迭代即可收敛. 迭代收敛的条件为不平衡载荷向量的2-范数小于10−6, 即$\left\| {\boldsymbol{R}} \right\| \leqslant {10^{ - 6}}$.

    为了保证优化过程的稳定性, 采用了连续化惩罚策略. 在优化前期, 采用较小的惩罚指数, 以降低优化过程陷入局部最优的概率. 后期则采用较大的惩罚指数, 以使得优化结果逐渐趋于0/1化. 具体的惩罚策略如下: (1) 初始的SIMP惩罚指数$p$ 和Heaviside惩罚指数$\beta $均设置为1.0; (2) 在后续优化迭代中, SIMP惩罚指数每5步迭代增加0.05; (3) 当SIMP惩罚指数达到3.0后, 其保持不变; Heaviside惩罚指数将会每5步增加0.5, 直至其达到12.0.

    本算例假设防弯器在位运行中会受到沿各个方向的波浪载荷的作用(各方向概率相等), 因此要求防弯器为旋转对称结构. 此时, 防弯器的设计问题可退化为如图9所示的2D平面问题. 防弯器的尺寸为${{{R}}_1} = 1800$, ${{{R}}_2} = 300$, ${{{L}}_1} = 2400$, ${{{L}}_2} = 900$. 其左端固定, 右端截面上施加方向竖直向上、大小为${{{a}}_{\text{p}}} = 200$的指定位移. 为了保证防弯器与柔性立管的配合, 在防弯器与立管的接触区域设置厚度为30的不可设计域.

    图  9  算例1中悬臂梁结构的设计域及边界条件
    Figure  9.  Geometry of the admissible domain and boundary conditions about the cantilever beam considered in example 1

    有限元计算中, 该结构被离散为20 400个4节点平面单元, MPI进程数为10. 在结构优化中, 约束材料用量体分比上限${\xi _0} = 0.5$, PDE过滤半径为30. 此外, 采用式(23)和式(26)构建了对称算子${\boldsymbol{M}}$, 以保证优化结果的对称性.

    优化后防弯器的拓扑优化构型如图10所示. 对图10所示构型进行旋转操作, 即可得3D的防弯器构型, 如图11所示. 由图9可知, 施加指定位移场${{\boldsymbol{a}}_{\text{p}}}$后, 防弯器的变形主要以弯曲为主, 同时伴随面内的剪切. 因此, 在优化后的结果中, 材料主要分布在防弯器的外侧, 用于提高截面的惯性矩, 进而增强结构的弯曲刚度. 同时, 在防弯器的内部与立管相连接的区域出现了交叉的斜撑结构, 其主要用于抵抗防弯器内部出现的剪切变形.

    图  10  算例1中防弯器2D拓扑优化结果
    Figure  10.  The 2D optimized result of bend stiffener in example 1
    图  11  算例1中旋转操作后所得防弯器优化结果剖面图
    Figure  11.  Sectional view about the optimized results of bend stiffener in example 1 obtained by using rotation operation

    考虑到由于波浪、海流等环境载荷在特定海域、特定浮式平台处的方向性, 使得作用在防弯器上的载荷往往也呈现特定的方向性. 此时, 最优的立管防弯器不再具有轴对称特征, 必须按照3D结构开展非线性拓扑优化设计. 在本算例的结构分析和优化中, 仍然采用如图9所示几何尺寸. 由于波浪载荷的往复性, 此时要求优化后的结果具有$ x $$ z $平面的对称性, 如图12所示. 与算例1相同, 结构左端固定, 右端立管截面施加沿$ y $轴竖直向上的指定位移场${{\boldsymbol{a}}_{\text{p}}} $, 大小为200. 由于载荷、边界条件及结构变形也具有$x - y$平面的对称性, 仅选取防弯器和立管中$z \geqslant 0$的部分进行有限元分析和优化.

    图  12  算例2中所采用的结构有限元模型
    Figure  12.  The FEM model used in example 2

    有限元计算中, 结构整体被离散为494 080个8节点实体单元, 此例求解时设置MPI进程数为50. 需注意的是, 对图10所示2D优化结果进行旋转操作形成3D结构时, 由于各平面单元在旋转路径上扫略所产生的体积不同, 会使得图11所示3D结构的材料用量体分比与图10有所不同. 本算例中, 优化所采用的材料用量体分比约束的上限与图11所示3D构型的实测体分比保持一致. 此外, PDE过滤半径设为30. 同时, 与算例1类似, 为了保证防弯器与立管的配合, 图12$\sqrt {{y^2} + {z^2}} \leqslant 330$的区域(含柔性立管区域)设定为不可设计域.

    考虑载荷方向的防弯器优化结果如图13图14所示. 其中, 图14中所示构型为图13所得防弯器构型经$x - y$平面对称操作后得到. 由图14可知, 优化后的防弯器呈现出明显的“H”型构型. 与算例1不同, 该构型主要用于抵抗沿$y$轴方向的载荷. 其中, 位于上下两侧的翼板主要用于提高结构的弯曲刚度; 而中间的腹板则用于抵抗弯曲变形过程所产生的剪切变形.

    图  13  算例2中防弯器3D拓扑优化结果
    Figure  13.  The 3D optimized result of bend stiffener in example 2
    图  14  算例2中对称操作后所得优化结果
    Figure  14.  The optimized results obtained by symmetry operation in example 2

    以下将对图11图14所示优化结果在不同波浪载荷方向下的承载能力进行对比. 如第3节所述, 当${{\boldsymbol{a}}_{\text{p}}}$的大小一定时, 位移施加端的支反力越大, 则对应结构的刚度也越大. 因此, 选用支反力指标对两种设计方案进行对比. 在对比分析中, 立管右端$ y $$ z $平面内所施加位移${{\boldsymbol{a}}_{\text{p}}}$的幅值固定为200. 提取立管右端面所有节点的支反力, 并且以所提取支反力的合力大小为指标, 对结构的承载能力进行量化. 不同$\theta $取值情况下, 两种设计方法的支反力对比如图15所示.

    图  15  不同载荷方向($\theta $)下, 算例1和算例2优化结果的刚度对比
    Figure  15.  Stiffness comparison between the optimized results of example 1 and 2 with different loading directions ($\theta $)

    图15可知, 由于图11所示优化结果为旋转对称结构(各方向载荷概率相等), 因此其对各个方向的波浪载荷具有相同的承载能力. 而对于图14所示优化结果, 在不同波浪载荷方向下的承载能力则呈现出明显的不同. 具体而言, 如图15所示在相当大的角度范围内(除$\theta $位于$\left[ {80,100} \right]$$\left[ {260,280} \right]$), 3D拓扑优化结果的承载能力则显著优于2D的优化设计; 而当载荷方向位于刚度薄弱区域, 即$\theta \in \left[ {80,100} \right] \cup $$ \left[ {260,280} \right]$范围时, 其承载能力弱于图11所示优化结果. 此时, 对防弯器的安装提出了更严格的要求, 即在安装过程中, 需要对海域波浪载荷方向进行预先勘测, 避免波浪载荷方向应避开防弯器的刚度薄弱区域.

    本文建立了Dirichlet边界条件下考虑材料和几何非线性的防弯器结构拓扑优化数学模型及框架. 为了保证优化结果的对称性, 研究引入了对称算子, 并给出了对称算子的具体构建方法. 采用伴随法进行了灵敏度分析, 并且引入了并行计算方法, 提高了结构分析和优化的效率. 数值算例部分, 给出了基于2D和3D拓扑优化的防弯器创新设计方案. 通过对此两种方案的承载能力进行了对比, 主要得到以下结论:

    (1) 当不考虑特定的波浪载荷方向时, 可采用旋转对称假设的防弯器结构创新设计, 其能够保证在任意载荷方向下, 防弯器均具有较好的承载能力.

    (2) 当海域内的波浪载荷具有特定的方向时, 可采用平面对称假设的防弯器设计方案. 需注意的是, 安装过程中载荷方向应避开防弯器的刚度薄弱区域.

  • 图  1   柔性立管与防弯器结构示意图

    Figure  1.   Schematic diagram about the flexible riser and bend stiffener

    图  2   切线刚度${{\boldsymbol{K}}_{\tan }}$与割线刚度${{\boldsymbol{K}}_{\cot }}$的定义, 其中${\boldsymbol{F}}$${\boldsymbol{u}}$分别为结构的载荷和位移

    Figure  2.   The definition of tangent stiffness ${{\boldsymbol{K}}_{\tan }}$ and secant stiffness ${{\boldsymbol{K}}_{\cot }}$, where ${\boldsymbol{F}}$ and ${\boldsymbol{u}}$ are the force and displacement of a structure, respectively

    图  3   海洋柔性立管、防弯器力学模型及边界条件

    Figure  3.   The mechanical model and boundary conditions of ocean flexible riser and the bend stiffener

    图  4   式(1)所定义目标函数的力学含义

    Figure  4.   Nature of the objective defined by Eq. (1)

    图  5   应变能与式(1)所定义目标函数的不同

    Figure  5.   The Difference between the strain energy and the objective defined by Eq. (1)

    图  6   切线刚度与式(1)所定义目标函数的不同

    Figure  6.   The difference between the tangent stiffness and the objective defined by Eq. (1)

    图  7   Dirichlet边界条件结构分析示意图

    Figure  7.   Illustration of Dirichlet boundary conditions

    图  8   Heaviside函数的性质

    Figure  8.   The characteristic of Heaviside function

    图  9   算例1中悬臂梁结构的设计域及边界条件

    Figure  9.   Geometry of the admissible domain and boundary conditions about the cantilever beam considered in example 1

    图  10   算例1中防弯器2D拓扑优化结果

    Figure  10.   The 2D optimized result of bend stiffener in example 1

    图  11   算例1中旋转操作后所得防弯器优化结果剖面图

    Figure  11.   Sectional view about the optimized results of bend stiffener in example 1 obtained by using rotation operation

    图  12   算例2中所采用的结构有限元模型

    Figure  12.   The FEM model used in example 2

    图  13   算例2中防弯器3D拓扑优化结果

    Figure  13.   The 3D optimized result of bend stiffener in example 2

    图  14   算例2中对称操作后所得优化结果

    Figure  14.   The optimized results obtained by symmetry operation in example 2

    图  15   不同载荷方向($\theta $)下, 算例1和算例2优化结果的刚度对比

    Figure  15.   Stiffness comparison between the optimized results of example 1 and 2 with different loading directions ($\theta $)

  • [1] 阎军, 英玺蓬, 步宇峰等. 深水柔性立管结构技术进展综述. 海洋工程装备技术, 2019, 6(6): 745-749 (Yan Jun, Ying Xipeng, Bu Yufeng, et al. Summary of development of flexible riser structural technology in deep water. Ocean Engineering Equipment and Technology, 2019, 6(6): 745-749 (in Chinese)
    [2]

    Patel MH, Seyed FB. Review of flexible riser modelling and analysis techniques. Engineering Structures, 1995, 17(4): 293-304 doi: 10.1016/0141-0296(95)00027-5

    [3] 李兰, 张荫纳. 静态柔性软管弯曲加强器的设计. 海洋工程装备与技术, 2014, 1(3): 223-226 (Li Lan, Zhang Mengna. Design of static bend stiffener for flexible pipes. Ocean Engineering Equipment and Technology, 2014, 1(3): 223-226 (in Chinese) doi: 10.3969/j.issn.2095-7297.2014.03.007
    [4]

    Deruntz JA. End effect bending stresses in cables. Journal of Applied Mechanics, 1969, 36(4): 750-756 doi: 10.1115/1.3564766

    [5]

    Boet WJC, Out JMM. Analysis of a flexible-riser top connection with bend restrictor//Proceedings of the Annual Offshore Technology Conference, Houston, 1990: 131-142

    [6]

    Tong DJ, Low YM, Sheehan JM. Nonlinear bend stiffener analysis using a simple formulation and finite element method. China Ocean Engineering, 2011, 25(4): 577-590 doi: 10.1007/s13344-011-0047-0

    [7]

    Drobyshevski Y. Investigation into non-linear bending of elastic bars with application to design of bend stiffeners. Marine Structures, 2013, 31: 102-130 doi: 10.1016/j.marstruc.2012.10.008

    [8]

    Tanaka RL, da Silveira LMY, Novaes JPZ, et al. Bending stiffener design through structural optimization//Proceedings of the International Conference on Offshore Mechanics and Arctic Engineering - OMAE, American Society of Mechanical Engineers Digital Collection, Hawaii, 2010, 3: 411-418

    [9]

    Yang H, Wang A. Fatigue reliability based design optimization of bending stiffener. Journal of Ship Research, 2012, 56(2): 120-128 doi: 10.5957/jsr.2012.56.2.120

    [10]

    Tang M, Yan J, Chen J, et al. Nonlinear analysis and multi-objective optimization for bend stiffeners of flexible riser. Journal of Marine Science and Technology, 2015, 20: 591-603

    [11]

    Cheng GD, Olhoff N. An investigation concerning optimal design of solid elastic plates. International Journal of Solids and Structures, Pergamon, 1981, 17(3): 305-323 doi: 10.1016/0020-7683(81)90065-2

    [12]

    Bendsøe MP, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering, 1988, 71(2): 197-224 doi: 10.1016/0045-7825(88)90086-2

    [13]

    Ye HL, Dai ZJ, Wang WW, et al. ICM method for topology optimization of multimaterial continuum structure with displacement constraint. Acta Mechanica Sinica, 2019, 35(3): 552-562 doi: 10.1007/s10409-018-0827-3

    [14] 隋允康, 叶红玲, 彭细荣等. 连续体结构拓扑优化应力约束凝聚化的ICM方法. 力学学报, 2007, 39(4): 554-563 (Sui Yunkang, Ye Hongling, Peng Xirong, et al. The ICM method for contimuum structural topology optimization with condenstation of stress constraints. Chinese Journal of Theoretical and Applied Mechanics, 2007, 39(4): 554-563 (in Chinese) doi: 10.3321/j.issn:0459-1879.2007.04.019
    [15]

    Wallin M, Ivarsson N, Tortorelli D. Stiffness optimization of non-linear elastic structures. Computer Methods in Applied Mechanics and Engineering, 2018, 330: 292-307 doi: 10.1016/j.cma.2017.11.004

    [16]

    Kemmler R, Lipka A, Ramm E. Large deformations and stability in topology optimization. Structural and Multidisciplinary Optimization, 2005, 30(6): 459-476 doi: 10.1007/s00158-005-0534-0

    [17]

    Bruns TE, Tortorelli DA. Topology optimization of non-linear elastic structures and compliant mechanisms. Computer Methods in Applied Mechanics and Engineering, 2001, 190(26-27): 3443-3459 doi: 10.1016/S0045-7825(00)00278-4

    [18]

    Wallin M, Ristinmaa M, Askfelt H. Optimal topologies derived from a phase-field method. Structural and Multidisciplinary Optimization, 2012, 45(2): 171-183 doi: 10.1007/s00158-011-0688-x

    [19]

    Buhl T, Pedersen CBW, Sigmund O. Stiffness design of geometrically nonlinear structures using topology optimization. Structural and Multidisciplinary Optimization, 2000, 19(2): 93-104 doi: 10.1007/s001580050089

    [20]

    Niu F, Xu S, Cheng G. A general formulation of structural topology optimization for maximizing structural stiffness. Structural and Multidisciplinary Optimization, 2011, 43(4): 561-572 doi: 10.1007/s00158-010-0585-8

    [21]

    Krenk S. Non-linear Modeling and Analysis of Solids and Structures. Cambridge: Cambridge University Press, 2009

    [22]

    Sigmund O. On the design of compliant mechanisms using topology optimization. Mechanics of Structures and Machines, 1997, 25(4): 493-524 doi: 10.1080/08905459708945415

    [23]

    Sigmund O, Maute K. Sensitivity filtering from a continuum mechanics perspective. Structural and Multidisciplinary Optimization, 2012, 46(4): 471-475 doi: 10.1007/s00158-012-0814-4

    [24]

    Bruns TE. A reevaluation of the SIMP method with filtering and an alternative formulation for solid-void topology optimization. Structural and Multidisciplinary Optimization, 2005, 30(6): 428-436 doi: 10.1007/s00158-005-0537-x

    [25]

    Abhyankar S, Brown J, Constantinescu EM, et al. PETSc/TS: A modern scalable ODE/DAE solver library. arXiv:1806.01437v1

    [26]

    Aage N, Andreassen E, Lazarov BS. Topology optimization using PETSc: An easy-to-use, fully parallel, open source topology optimization framework. Structural and Multidisciplinary Optimization, 2015, 51(3): 565-572 doi: 10.1007/s00158-014-1157-0

    [27]

    Lazarov BS, Sigmund O. Filters in topology optimization based on Helmholtz-type differential equations. International Journal for Numerical Methods in Engineering, 2011, 86(6): 765-781 doi: 10.1002/nme.3072

    [28]

    Duan Z, Yan J, Zhao G. Integrated optimization of the material and structure of composites based on the heaviside penalization of discrete material model. Structural and Multidisciplinary Optimization, 2015, 51(3): 721-732 doi: 10.1007/s00158-014-1168-x

    [29]

    Xu S, Cai Y, Cheng G. Volume preserving nonlinear density filter based on heaviside functions. Structural and Multidisciplinary Optimization, 2010, 41(4): 495-505 doi: 10.1007/s00158-009-0452-7

    [30]

    Wang FW, Sigmund O, Jensen JS. Design of materials with prescribed nonlinear properties. Journal of the Mechanics and Physics of Solids, 2014, 69(1): 156-174

  • 期刊类型引用(5)

    1. 岳宝增,马伯乐,唐勇,刘峰. 液体大幅晃动情形的航天器刚液耦合动响应仿真分析. 宇航学报. 2022(02): 173-182 . 百度学术
    2. 孙船斌,沈民民,童宝宏,黄虎,张坤. 变工况下车辆燃油箱多孔挡板抑浪机理研究. 中国机械工程. 2022(11): 1377-1385 . 百度学术
    3. 王震,祝恒佳,陈晓宇,张云清. 基于交叉型双气室空气互联悬架的全地形车侧倾特性研究. 力学学报. 2020(04): 996-1006 . 本站查看
    4. 李海波,刘世兴,宋海燕,梁立孚. 非保守非线性刚-弹-液-控耦合分析动力学及其应用研究. 力学学报. 2020(04): 932-944 . 本站查看
    5. 司震,钱霙婧,杨晓东,张伟. 基于参激共振的受摄小行星悬停轨道设计方法. 力学学报. 2020(06): 1774-1788 . 本站查看

    其他类型引用(8)

图(15)
计量
  • 文章访问数:  1149
  • HTML全文浏览量:  483
  • PDF下载量:  134
  • 被引次数: 13
出版历程
  • 收稿日期:  2021-11-10
  • 录用日期:  2021-12-15
  • 网络出版日期:  2021-12-16
  • 刊出日期:  2022-04-17

目录

/

返回文章
返回