木贼属植物仿生薄壁结构的耐撞性优化

刘飞明; 雷建银; 乔力; 刘志芳

doi:10.11858/gywlxb.20220516

木贼属植物仿生薄壁结构的耐撞性优化

doi: 10.11858/gywlxb.20220516

太原理工大学机械与运载工程学院应用力学研究所, 山西太原　030024

基金项目: 国家自然科学基金（11902215, 11772216）

详细信息

作者简介:
刘飞明（1995-），男，硕士研究生，主要从事薄壁结构的耐撞性研究. E-mail：1637752941@qq.com

通讯作者:
雷建银（1989-），男，博士，副教授，主要从事仿生结构的力学行为研究.E-mail：leijianyin@tyut.edu.cn

中图分类号: O341
计量
- 文章访问数: 266
- HTML全文浏览量: 168
- PDF下载量: 32
出版历程
- 收稿日期: 2022-02-18
- 修回日期: 2022-05-20
- 刊出日期: 2022-10-11

Crashworthiness Optimization of Horsetail-Bionic Thin-Walled Structures

Institute of Applied Mechanics, College of Mechanical and Vehicle Engineering, Taiyuan University of Technology, Taiyuan 030024, Shanxi, China

摘要

摘要: 利用ABAQUS有限元软件构建木贼属植物仿生薄壁结构（horsetail-bionic thin-walled structure，HBTS）在侧向冲击下的数值模型，分析了结构的壁厚、内径和肋骨数对其耐撞性能和变形模态的影响。结果表明：肋骨数和整体壁厚的增加会明显提高HBTS的比吸能和最大峰值载荷，HBTS不同部分的壁厚变化显著影响其变形模态和耐撞性能。为了综合考虑HBTS各部分的壁厚、肋骨数、内径5个参数对耐撞性能的影响，集成modeFRONTIER优化软件与ABAQUS软件，通过参数化建模在设计空间上构建有限元模型，建立比吸能和最大峰值载荷的Kriging代理模型。使用基于Kriging代理模型的多目标优化方法获取Pareto前沿，以同时实现比吸能最大化和峰值载荷最小化。最后分析了Pareto前沿面上各HBTS的设计参数分布情况，并验证了优化结果，该方法可为薄壁结构的优化设计提供新思路。
- 仿生薄壁结构 /
- 多目标优化 /
- modeFRONTIER /
- Pareto前沿
Abstract: Numerical model of horsetail-bionic thin-walled structure (HBTS) under the lateral impact was constructed using ABAQUS. The effects of wall thickness, inner diameter and number of ribs on the crashworthiness performance and deformation modes were analyzed. The results indicate that the specific energy absorption and peak load of HBTS can be significantly enhanced with the increase of the number of ribs and the overall wall thickness. The changes in the wall thickness of each part significantly affects its deformation mode and crashworthiness performance. Based on the above results, the optimization software modeFRONTIER and the finite element analysis software ABAQUS were integrated to explore the influence of five design parameters, in terms of wall thickness, number of ribs, inner diameter and so on. Finite element models were uniformly distributed on the design space through parametric modeling method, hence the Kriging surrogate models for the specific energy absorption and peak load were established. Then, Pareto front was obtained using Kriging surrogate model-based multi-objective optimization method for the specific energy absorption maximization and peak load minimization simultaneously in one model. Finally, the distribution of each HBTS’s design parameters on Pareto front was analyzed and the optimization results were verified. The method is expected to provide new thoughts for the optimization design of the thin-walled structure.
- bionic thin-walled structure /
- multi-objective optimization /
- modeFRONTIER /
- Pareto front

HTML全文

复合材料板具有较高的强度质量比、良好的耐腐蚀性和优异的可设计性，被广泛应用于航空航天和工业制造等领域^[1]。在实际使用中复合材料板经常受到不同形式的冲击荷载，从而产生振动和屈曲问题，因此冲击载荷作用下复合材料板的动力稳定性问题备受关注。

近年来关于复合材料板的研究越来越多，尤其是冲击荷载作用下复合材料板的动力稳定性问题研究^[2]，对工程部件结构设计和使用具有重要的意义。Sun等^[3]研究了在加热环境中应力波对功能梯度圆柱壳轴向冲击屈曲的影响；毛柳伟等^[4]对弹性直杆在应力波作用下的动力分叉屈曲进行了分析与探讨，提出了求解应力波作用下直杆动力屈曲的数值方法；Lepik^[5]讨论了在应力波影响下轴向压缩的弹塑性梁的屈曲；Zhang等^[6]分析了不确定初始几何缺陷对薄板屈曲的影响；Abdelaziz等^[7]利用双曲线剪切变形理论，分析了在各种边界条件下复合材料板的弯曲变形和屈曲；Kouchakzadeh等^[8]采用线性和旋转弹簧的均匀分布来模拟边界条件，对矩形层压复合板的屈曲进行了分析；Czapski等^[9]通过数值和实验方法，研究了残余应力对压缩至破坏期间薄壁层压板屈曲性能的影响。

在实际工程中，复合材料板多应用于振动环境，其在动力响应下的动态特性和振动分析是必不可少的，因而对该类材料的振动屈曲研究至关重要。Kuo^[10]研究了两种非均匀分布纤维复合材料板的振动屈曲问题，Villarreal等^[11]对典型正交异性板的本征频率和振动屈曲进行了理论分析，Eftekhari等^[12]提出了通过组合应用有限元方法和微分正交方法求解矩形板的振动屈曲问题，Rehman等^[13]探讨了壳体结构的缺陷和损坏对结构振动屈曲的影响，Sayyad等^[14]将三角剪切变形理论应用于复合板的变形和振动屈曲研究。

关于复合材料板的振动屈曲问题已开展了较多的研究，但大多未考虑应力波效应对振动屈曲的影响，而动力屈曲一般与应力波相联系且具有局部发生的特点，研究含初始缺陷的复合材料板能更好地揭示实际工程中复合材料板在不同工况下发生动力屈曲的机理。基于此，本研究利用Kirchhoff薄板理论和Hamilton原理，考虑应力波效应，建立含初始几何缺陷的四边简支复合材料板的振动控制方程，得到板的屈曲临界荷载表达式，在此基础上通过数值计算讨论初始几何缺陷、振型函数初相位、铺层角度、屈曲模态阶数和铺层层数对复合材料板振动屈曲临界荷载的影响，为工程实际提供理论依据。

1. 复合材料板的控制方程

复合材料板在x = L_a处为固定边界条件，其余3边为简支边界条件，在 $z$ = 0的中性面上受x方向的面内阶跃荷载N作用，如图1所示。板在 $z$ 方向上含初始几何缺陷w₁，且x、y、 $z$ 方向的位移分别为u、v、w。复合材料板的长度、宽度和厚度分别为L_a、L_b、h，由n层单层板组成， $\theta$ 为单层板的铺层角度，即纤维材料铺设方向与x方向的夹角（见图1）。

图 1 复合材料板结构示意图

Figure 1. Schematic diagram of composite plate structure

下载: 全尺寸图片幻灯片

根据Kirchhoff薄板理论及经典弹性理论，复合材料板的位移与应变、弹性曲面的曲率和扭率的表达式为

${\begin{cases} {w = {w_1} + {w_0}} \\ {u = {u_0} - z\displaystyle\frac{{\partial \left( {w - {w_1}} \right)}}{{\partial x}}} \\ {v = {v_0} - z\displaystyle\frac{{\partial \left( {w - {w_1}} \right)}}{{\partial y}}} \end{cases}}$

(1)

${\begin{cases} {\varepsilon _x^0 = \displaystyle\frac{{\partial {u_0}}}{{\partial x}}} \\ {\varepsilon _y^0 = \displaystyle\frac{{\partial {v_0}}}{{\partial y}}} \\ {\gamma _{xy}^0 = \displaystyle\frac{{\partial {u_0}}}{{\partial y}} + \displaystyle\frac{{\partial {v_0}}}{{\partial x}}} \end{cases}}$

(2)

${\begin{cases} {{\kappa _x} = - \displaystyle\frac{{{\partial ^2}\left( {w - {w_1}} \right)}}{{\partial {x^2}}}} \\ {{\kappa _y} = - \displaystyle\frac{{{\partial ^2}\left( {w - {w_1}} \right)}}{{\partial {y^2}}}} \\ {{\kappa _{xy}} = - \displaystyle\frac{{2{\partial ^2}\left( {w - {w_1}} \right)}}{{\partial x\partial y}}} \end{cases}}$

(3)

${\begin{cases} {{\varepsilon _x} = \varepsilon _x^0 + z{\kappa _x}} \\ {{\varepsilon _y} = \varepsilon _y^0 + z{\kappa _y}} \\ {{\gamma _{xy}} = \gamma _{_{xy}}^0 + z{\kappa _{xy}}} \end{cases}}$

(4)

式中：u₀、v₀、w₀分别为复合材料板在x、y、 $z$ 方向上的中面位移， $\varepsilon _x^0$ 、 $\varepsilon _y^0$ 、 $\gamma _{xy}^0$ 分别为复合材料板中面应变分量， ${\kappa _x}$ 、 ${\kappa _y}$ 、 ${\kappa _{xy}}$ 为中面的曲率和扭率， ${\varepsilon _x}$ 、 ${\varepsilon _y}$ 、 ${\gamma _{xy}}$ 为复合材料板任意一点的应变。

板的内力（N_x、N_y、N_xy）与内力矩（M_x、M_y、M_xy）之间的关系为

$\left[ \begin{array}{l} {N_x} \\ {N_y} \\ {N_{xy}} \\ {M_x} \\ {M_y} \\ {M_{xy}} \end{array} \right] = \left[ {\begin{array}{*{20}{c}} {{A_{11}}}&{{A_{12}}}&{{A_{16}}}&{{B_{11}}}&{{B_{12}}}&{{B_{16}}} \\ {{A_{12}}}&{{A_{22}}}&{{A_{26}}}&{{B_{12}}}&{{B_{22}}}&{{B_{26}}} \\ {{A_{16}}}&{{A_{26}}}&{{A_{66}}}&{{B_{16}}}&{{B_{26}}}&{{B_{66}}} \\ {{B_{11}}}&{{B_{12}}}&{{B_{16}}}&{{D_{11}}}&{{D_{12}}}&{{D_{16}}} \\ {{B_{12}}}&{{B_{22}}}&{{B_{26}}}&{{D_{12}}}&{{D_{22}}}&{{D_{26}}} \\ {{B_{16}}}&{{B_{26}}}&{{B_{66}}}&{{D_{16}}}&{{D_{26}}}&{{D_{66}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\varepsilon _x^0} \\ {\varepsilon _y^0} \\ {\gamma _{xy}^0} \\ {{\kappa _x}} \\ {{\kappa _y}} \\ {{\kappa _{xy}}} \end{array}} \right]$

(5)

式中： ${A}_{ij}$ 、 ${B}_{ij}$ 、 ${D}_{ij} \left(i,j=1,2,6\right)$ 分别表示板的拉伸刚度、耦合刚度和弯曲刚度系数^[15]，表达式如下

${\begin{cases} {{A_{ij}} = \displaystyle\sum\limits_{k = 1}^n {{{\overline Q_{ij}^k }}\left( {{h_k} - {h_{k - 1}}} \right)} } \\ {{B_{ij}} = \displaystyle\frac{1}{2}\displaystyle\sum\limits_{k = 1}^n {{{\overline Q_{ij}^k }}( {{h^2_k} - h_{k - 1}^2} )} } \\ {{D_{ij}} = \displaystyle\frac{1}{3}\displaystyle\sum\limits_{k = 1}^n {{{\overline Q_{ij}^k }}( {{h^3_k} - h_{k - 1}^3} )} } \end{cases}}$

(6)

式中： ${\overline Q_{ij}^k }$ 为复合材料板第 $k$ 层的偏轴刚度系数

$\overline {\boldsymbol{Q}} = {{\boldsymbol{P}}^{ - 1}}{\boldsymbol{Q}}{({{\boldsymbol{P}}^{ - 1}})^{\rm{T}}}$

(7)

式中：P为坐标转换矩阵，Q为刚度矩阵。

${\boldsymbol{P}} = \left[ {\begin{array}{*{20}{c}} {{{\cos }^2}\theta }&{{{\sin }^2}\theta }&{2\sin \theta \cos \theta } \\ {{{\sin }^2}\theta }&{{{\cos }^2}\theta }&{ - 2\sin \theta \cos \theta } \\ { - \sin \theta \cos \theta }&{\sin \theta \cos \theta }&{{{\cos }^2}\theta - {{\sin }^2}\theta } \end{array}} \right]$

(8)

${\boldsymbol{Q}}{\rm{ = }}\left[ {\begin{array}{*{20}{c}} {{Q_{11}}}&{{Q_{12}}}&0 \\ {{Q_{12}}}&{{Q_{22}}}&0 \\ 0&0&{{Q_{66}}} \end{array}} \right]$

(9)

考虑材料为正交各向异性材料，设E₁、E₂、G₁₂、μ₁₂、μ₂₁分别为板材料x和y方向的拉压弹性模量、剪切弹性模量、主泊松比和副泊松比，则有

${Q_{11}}{\rm{ = }}\frac{{{E_1}}}{{1 - {\mu _{12}}{\mu _{21}}}},\;\;\;{Q_{22}}{\rm{ = }}\frac{{{E_2}}}{{1 - {\mu _{12}}{\mu _{21}}}},\;\;\;{Q_{12}}{\rm{ = }}\frac{{{\mu _{12}}{E_2}}}{{1 - {\mu _{12}}{\mu _{21}}}}{\rm{ = }}\frac{{{\mu _{21}}{E_1}}}{{1 - {\mu _{12}}{\mu _{21}}}},\;\;\;{Q_{66}}{\rm{ = }}{G_{12}}$

(10)

复合材料板在左端受面内冲击荷载N作用（见图1）时，应力波沿x方向在板内传播，其应力变化如图2所示。

图 2 应力波传播示意图

Figure 2. Schematic diagram of stress wave propagation

下载: 全尺寸图片幻灯片

当应力波传播至波阵面位置L_cr（临界长度）时，板发生振动屈曲，板的内力N_t和应力波波速c分别表示为

${N_{\rm{t}}} = {\begin{cases} N&{0 \leqslant x \leqslant {L_{{\rm{cr}}}}} \\ 0&{x > {L_{{\rm{cr}}}}} \end{cases}} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} c{\rm{ = }}\sqrt {\frac{{{A_{11}}}}{{\rho h}}}$

(11)

板发生振动屈曲时的变形能可以表示为

$\begin{split} U = &\displaystyle\frac{1}{2}\displaystyle\int_0^{{L_{\rm{cr}}}} {\displaystyle\int_0^{{L_{\rm{b}}}} {\left( {{N_x}\varepsilon _x^0 + {N_y}\varepsilon _y^0 + {N_{xy}}\gamma _{xy}^0 + {M_x}{\kappa _x} + {M_y}{\kappa _y} + {M_{xy}}{\kappa _{xy}}} \right)} } {\rm{d}}x{\rm{d}}y {\rm{ = }}\\ &\displaystyle\frac{1}{2}\displaystyle\int_0^{{L_{\rm{cr}}}} {\displaystyle\int_0^{{L_{\rm{b}}}} {\left[ {{N_x}\frac{{\partial {u_0}}}{{\partial x}} + {N_y}\frac{{\partial {v_0}}}{{\partial y}} + {N_{xy}}\left( {\frac{{\partial {u_0}}}{{\partial x}} + \frac{{\partial {v_0}}}{{\partial y}}} \right) + {M_x}\left( { - \frac{{{\partial ^2}{w_0}}}{{\partial {x^2}}}} \right) + {M_y}\left( { - \frac{{{\partial ^2}{w_0}}}{{\partial {y^2}}}} \right) + {M_{xy}}\left( { - \frac{{2{\partial ^2}{w_0}}}{{\partial x\partial y}}} \right)} \right]} } {\rm{d}}x{\rm{d}}y \end{split}$

(12)

发生振动屈曲时的动能（考虑转动惯量）可以表示为

$\begin{split} T =& \displaystyle\frac{1}{2}\displaystyle\int_{ - {h}/{2}}^{{h}/{2}} {\int_0^{{L_{\rm{cr}}}} {\int_0^{{L_{\rm{b}}}} {{\rho _{\left( k \right)}}\left[ {{{\left( {\displaystyle\frac{{\partial u}}{{\partial t}}} \right)}^2} + {{\left( {\displaystyle\frac{{\partial v}}{{\partial t}}} \right)}^2} + {{\left( {\displaystyle\frac{{\partial {w_0}}}{{\partial t}}} \right)}^2}} \right]} } } {\rm{d}}x{\rm{d}}y{\rm{d}}z {\rm{ = }}\\ &\displaystyle\frac{1}{2}\int_0^{{L_{\rm{cr}}}} {\int_0^{{L_{\rm{b}}}} {\left[ {{I_0}{{\left( {\displaystyle\frac{{\partial {u_0}}}{{\partial t}}} \right)}^2} + {I_0}{{\left( {\displaystyle\frac{{\partial {v_0}}}{{\partial t}}} \right)}^2} + {I_0}{{\left( {\displaystyle\frac{{\partial {w_0}}}{{\partial t}}} \right)}^2} - 2{I_1}\displaystyle\frac{{\partial {u_0}}}{{\partial t}}\displaystyle\frac{{{\partial ^2}{w_0}}}{{\partial x\partial t}} - 2{I_1}\displaystyle\frac{{\partial {v_0}}}{{\partial t}}\displaystyle\frac{{{\partial ^2}{w_0}}}{{\partial y\partial t}} + {I_2}{{\left( {\displaystyle\frac{{{\partial ^2}{w_0}}}{{\partial x\partial t}}} \right)}^2} + {I_2}{{\left( {\displaystyle\frac{{{\partial ^2}{w_0}}}{{\partial y\partial t}}} \right)}^2}} \right]} } {\rm{d}}x{\rm{d}}y{\rm{ }} \end{split}$

(13)

式中： $\left( {{I_0},\; {I_1},\; {I_2}} \right) = \displaystyle\sum\limits_{k = 1}^{{N_k}} {\int_{{h_{k - 1}}}^{{h_k}} {{\rho _{\left( k \right)}}} } \left( {1, \;z, \;{z^2}} \right){\rm{d}}z$ ， ${{\rho _{\left( k \right)}}}$ 为第k层材料的密度。

板发生振动屈曲时的外力功可以表示为

$W = \frac{1}{2}\int_0^{{L_{\rm{cr}}}} {\int_0^{{L_{\rm{b}}}} {{N_{\rm{t}}}{{\left( {\frac{{\partial w}}{{\partial x}}} \right)}^2}} } {\rm{d}}x{\rm{d}}y$

(14)

考虑Hamilton变分原理，即

${\rm{\delta}} \int_{{{{t}}_1}}^{{{{t}}_2}} {\left( {T - U + W} \right)} {\rm{ d}}t = 0$

(15)

将式(12)～式(14)代入式(15)并进行变分计算，可得

$\frac{{\partial {N_x}}}{{\partial x}} + \frac{{\partial {N_{xy}}}}{{\partial x}} - {I_0}\frac{{{\partial ^2}{u_0}}}{{\partial {t^2}}} + {I_1}\frac{{{\partial ^3}{w_0}}}{{\partial x\partial {t^2}}}{\rm{ = }}0$

(16)

$\frac{{\partial {N_y}}}{{\partial y}} + \frac{{\partial {N_{xy}}}}{{\partial y}} - {I_0}\frac{{{\partial ^2}{v_0}}}{{\partial {t^2}}} + {I_1}\frac{{{\partial ^3}{w_0}}}{{\partial y\partial {t^2}}}{\rm{ = }}0$

(17)

$\frac{{{\partial ^2}{M_x}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{M_y}}}{{\partial {y^2}}} - 2\frac{{{\partial ^2}{M_{xy}}}}{{\partial x\partial y}} - {N_{\rm t}}\left( {\frac{{{\partial ^2}{w_0}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{w_1}}}{{\partial {x^2}}}} \right) = {I_0}\frac{{{\partial ^2}{w_0}}}{{\partial {t^2}}}{\kern 1pt} + {I_1}\left( {\frac{{{\partial ^3}{u_0}}}{{\partial x\partial {t^2}}} + \frac{{{\partial ^3}{v_0}}}{{\partial x\partial {t^2}}}} \right) - {I_2}\left( {\frac{{{\partial ^4}{w_0}}}{{\partial {x^2}\partial {t^2}}} + \frac{{{\partial ^4}{w_0}}}{{\partial {y^2}\partial {t^2}}}} \right)$

(18)

对于正交各向异性正规对称正交铺设的复合材料板，其刚度矩阵满足^[16]

${D}_{16}={D}_{26}=0,\;\;{A}_{16}={A}_{26}=0,\;\;{B}_{ij}=0$

(19)

根据 Kirchhoff 薄板理论及经典弹性理论，薄板中面在变形过程中没有伸长变形，将板的本构关系代入式(16)～式(18)中，略去含 ${u_0}$ 和 ${v_0}$ 的项，得到复合材料板在面向阶跃荷载激励下的控制方程

${D_{11}}\frac{{{\partial ^4}{w_0}}}{{\partial {x^4}}} + \left( {2{D_{12}} + 4{D_{66}}} \right)\frac{{{\partial ^4}{w_0}}}{{\partial {x^2}\partial {y^2}}} + {D_{22}}\frac{{{\partial ^4}{w_0}}}{{\partial {y^4}}} + {N_{\rm t}}\left( {\frac{{{\partial ^2}{w_0}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{w_1}}}{{\partial {x^2}}}} \right) - {I_0}\frac{{{\partial ^2}{w_0}}}{{\partial {t^2}}} + {I_2}\left( {\frac{{{\partial ^4}{w_0}}}{{\partial {x^2}\partial {t^2}}} + \frac{{{\partial ^4}{w_0}}}{{\partial {y^2}\partial {t^2}}}} \right) = 0$

(20)

设三边简支和应力波传播到L_cr（ ${w_0}\left( {{L_{{\rm{cr}}}},y,t} \right) = {w'_0}\left( {{L_{{\rm{cr}}}},y,t} \right) = 0$ ）时的振型函数^[17]为

${w_0}\left( {x,y,t} \right) = {R_{ij}}\left[ {\sin \frac{{i{\text{π}} x}}{{{L_{{\rm{cr}}}}}} + \frac{i}{{i + 1}}\sin \frac{{\left( {i + 1} \right){\text{π}}x }}{{{L_{{\rm{cr}}}}}}} \right]\sin \frac{{j{\text{π}} {{y}}}}{{{L_{\rm{b}}}}}\sin \left( {\omega {\kern 1pt} t + \varphi } \right){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt}$

(21)

由屈曲模态确定的缺陷分布形式是板结构最有可能发生的屈曲形式，能够很好地确定结构的缺陷敏感性^[18]。对于复合材料板在制造过程中出现的初始几何缺陷，引入屈曲模态的 $\varepsilon$ 倍变形作为初始几何缺陷^[19]，可以表示为

${w_1}\left( {x,y} \right) = \varepsilon {R_{ij}}\left[ {\sin \frac{{i{\text{π}}x }}{{{L_{{\rm{cr}}}}}} + \frac{i}{{i + 1}}\sin \frac{{\left( {i + 1} \right){\text{π}} x}}{{{L_{{\rm{cr}}}}}}} \right]\sin \frac{{j{\text{π}} {{y}}}}{{{L_{\rm{b}}}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt}$

(22)

式中：i、j为屈曲模态阶数， $i,j = 1,2,3, \cdots$ ；R_ij为板的第(i, j)阶模态幅值； $\varepsilon$ 表示初始几何缺陷系数。

根据式(21)和式(22)，利用棣莫弗公式对控制方程式(20)中的各项进行求导并化简，得到

$\begin{split} &{D_{11}}{\left( {\displaystyle\frac{{\text{π}} }{{{L_{\rm{cr}}}}}} \right)^4}\left[ {{i^4} + {{\left( {i + 1} \right)}^4}} \right]\sin \left( {\omega t + \varphi } \right) + \left( {2{D_{12}} + 4{D_{66}}} \right){\left( {\displaystyle\frac{{\text{π}} }{{{L_{\rm{cr}}}}}} \right)^2}\left[ {{i^2} + {{\left( {i + 1} \right)}^2}} \right]{\left( {\displaystyle\frac{{j{\text{π}} }}{{{L_{\rm{b}}}}}} \right)^2}\sin \left( {\omega t + \varphi } \right) + \\ &2{D_{22}}{\left( {\displaystyle\frac{{j{\text{π}} }}{{{L_{\rm{b}}}}}} \right)^4}\sin \left( {\omega t + \varphi } \right) - {N_{\rm{t}}}{\left( {\displaystyle\frac{{\text{π}} }{{{L_{\rm{cr}}}}}} \right)^2}\left[ {{i^2} + {{\left( {i + 1} \right)}^2}} \right]\left[ {\sin \left( {\omega t + \varphi } \right) + \varepsilon } \right] + 2{I_0}{\omega ^2}\sin \left( {\omega t + \varphi } \right) + \\ &{I_2}{\omega ^2}\left\{ {{{\left( {\displaystyle\frac{{\text{π}} }{{{L_{\rm{cr}}}}}} \right)}^2}\left[ {{i^2} + {{\left( {i + 1} \right)}^2}} \right] + 2{{\left( {\displaystyle\frac{{j{\text{π}} }}{{{L_{\rm{b}}}}}} \right)}^2}} \right\}\sin \left( {\omega t + \varphi } \right) = 0 \end{split}$

(23)

根据式(23)可以得到N的表达式

$\begin{split} N{\rm{ = }}&\left\{ {{D_{11}}{{\left( {\displaystyle\frac{{\text{π}} }{{{L_{\rm{cr}}}}}} \right)}^4}\left[ {{i^4} + {{\left( {i + 1} \right)}^4}} \right] + \left( {2{D_{12}} + 4{D_{66}}} \right){{\left( {\displaystyle\frac{{\text{π}} }{{{L_{\rm{cr}}}}}} \right)}^2}\left[ {{i^2} + {{\left( {i + 1} \right)}^2}} \right]{{\left( {\displaystyle\frac{{j{\text{π}} }}{{{L_{\rm{b}}}}}} \right)}^2}} \right. + 2{D_{22}}{\left( {\displaystyle\frac{{j{\text{π}} }}{{{L_{\rm{b}}}}}} \right)^4} + 2{I_0}{\omega ^2} + \\ &{\left.{I_2}{\omega ^2}\left\{ {{{\left( {\displaystyle\frac{{\text{π}} }{{{L_{\rm{cr}}}}}} \right)}^2}\left[ {{i^2} + {{\left( {i + 1} \right)}^2}} \right] + 2{{\left( {\displaystyle\frac{{j{\text{π}} }}{{{L_{\rm{b}}}}}} \right)}^2}} \right\}\right\}} \displaystyle\frac{{\sin \left( {\omega t + \varphi } \right)}}{{{{\left( {\displaystyle\frac{{\text{π}} }{{{L_{\rm{cr}}}}}} \right)}^2}\left[ {{i^2} + {{\left( {i + 1} \right)}^2}} \right]\left[ {\sin \left( {\omega t + \varphi } \right) + \varepsilon } \right]}} \end{split}$

(24)

板发生屈曲时，临界条件为 $\omega = 0$ ^[20]，代入式(24)可得振动屈曲临界荷载为

${N_{\rm{cr}}}{\rm{ = }}{D_{11}}{\left( {\displaystyle\frac{{\text{π}} }{{{L_{\rm{cr}}}}}} \right)^2}\displaystyle\frac{{\left[ {{i^4} + {{\left( {i + 1} \right)}^4}} \right]\sin \varphi }}{{\left[ {{i^2} + {{\left( {i + 1} \right)}^2}} \right]\left( {\sin \varphi + \varepsilon } \right)}} + \left( {2{D_{12}} + 4{D_{66}}} \right){\left( {\displaystyle\frac{{j{\text{π}} }}{{{L_{\rm{b}}}}}} \right)^2}\displaystyle\frac{{\sin \varphi}}{{\left( {\sin \varphi + \varepsilon } \right)}} + 2{D_{22}}{\left( {\displaystyle\frac{{{j^2}{\text{π}} }}{{L_{\rm{b}}^2}}} \right)^2}L_{{\rm{cr}}}^2\displaystyle\frac{{\sin \varphi}}{{\left[ {{i^2} + {{\left( {i + 1} \right)}^2}} \right]\left( {\sin \varphi + \varepsilon } \right)}}$

(25)

2. 算例分析

利用MATLAB数值分析应力波未反射时初始几何缺陷、初相位、铺层角度、屈曲模态阶数、铺设厚度以及铺层层数对复合材料板振动屈曲临界荷载的影响，使用的材料参数见表1^[21]。

表 1 复合材料板参数^[21]

Table 1. Material parameters of composite plate^[21]

E₁/GPa	E₂/GPa	G₁₂/GPa	μ₁₂	L_a/m	L_b/m
140.0	8.6	5.0	0.35	0.60	0.50

下载: 导出CSV

| 显示表格

图3显示了复合材料板的模态阶数i为1、2、3，j为1时的屈曲模态。当板的x方向模态增大时，x方向的屈曲模态第一峰值增大且波数增加，而y方向的屈曲模态呈正对称分布。模态阶数的增加使板振动的屈曲模态变得更复杂。

图 3 x方向模态取值增大时板的屈曲模态

Figure 3. Buckling mode of composite plate with increasing mode value in x direction

下载: 全尺寸图片幻灯片

设置7组算例，分别以初始几何缺陷、初相位、铺层角度、x和y两个方向屈曲模态阶数、铺层层数及铺设厚度为变量进行算例分析，研究以上因素对板振动屈曲临界荷载的影响，算例参数见表2。

表 2 算例分析参数表

Table 2. Example analysis parameter table

Group	Initial defect coefficient	Order of mode		Laying angle/(°)	Initial phase	Number of layers laid	Thickness of the plate/m
Group	Initial defect coefficient	x direction	y direction	Laying angle/(°)	Initial phase	Number of layers laid	Thickness of the plate/m
A	Variable	i = 1	j = 1	[0, 0, 0, 0, 0]	${\text{π}}$ /2	5	0.01
B	0.1	Variable	j = 1	[0, 0, 0, 0, 0]	${\text{π}}$ /2	5	0.01
C	0.1	i =1	Variable	[0, 0, 0, 0, 0]	${\text{π}}$ /2	5	0.01
D	0.1	i =1	j = 1	Variable	${\text{π}}$ /2	5	0.01
E	0.1	i =1	j = 1	[0, 0, 0, 0, 0]	Variable	5	0.01
F	0.1	i =1	j = 1	[0, 0, 0, 0, 0]	${\text{π}}$ /2	Variable	0.01
G	0.1	i =1	j = 1	[0, 0, 0, 0, 0]	${\text{π}}$ /2	5	Variable

下载: 导出CSV

| 显示表格

将A组数据代入式(25)中，可以得到不同初始缺陷系数对复合材料板振动屈曲的影响，如图4所示。由L_cr-N_cr曲线可知：在应力波传播过程中，N_cr呈指数下降，分为敏感区和非敏感区。应力波在L_cr < 0.4 m区域传播时，N_cr的变化较陡峭，该区域为敏感区；应力波在L_cr > 0.4 m区域传播时，N_cr的变化趋于平缓，该区域为非敏感区，因此敏感分界点为0.4 m。当选取的初始缺陷系数增大时，临界荷载N_cr也随之增大。在敏感区，初始缺陷系数对临界荷载N_cr的影响较大，且随应力波传播呈减小趋势。此外，初始缺陷系数对非敏感区的影响较小。图4表明，初始几何缺陷系数越大，板越容易发生屈曲。

图 4 不同初始缺陷系数条件下N_cr与L_cr的关系曲线

Figure 4. Relationship between N_cr and L_cr under different initial defect coefficients

下载: 全尺寸图片幻灯片

将B组数据代入式(25)中，可以得到不同x方向模态阶数对复合材料板振动屈曲的影响，如图5所示。由L_cr-N_cr曲线可知：当选取的x方向模态阶数增大时，临界荷载N_cr随之明显增大。在敏感区，x方向模态阶数对临界荷载N_cr的影响很大，并随应力波的传播不断减小，到达非敏感区之后影响较小并趋于稳定。图5表明，x方向模态阶数的增加会显著增大板的屈曲临界荷载。

图 5 x方向模态阶数不同时N_cr与L_cr的关系曲线

Figure 5. Relationship between N_cr and L_cr with different order of modes in x direction

下载: 全尺寸图片幻灯片

将C组数据代入式(25)中，可以得到不同的y方向模态阶数对复合材料板振动屈曲的影响，如图6所示。由L_cr-N_cr曲线可知：当选取的y方向模态阶数增大时，临界荷载N_cr也随之变大。在应力波传播过程中，在敏感区y方向模态阶数对临界载荷基本没有影响，而在非敏感区有极小的影响。图6表明，y方向模态阶数的变化对板屈曲临界荷载基本没有影响。

图 6 y方向模态阶数不同时N_cr与L_cr的关系曲线

Figure 6. Relationship between N_cr and L_cr with different order of modes in y direction

下载: 全尺寸图片幻灯片

将D组数据代入式(25)中，得到不同铺层角度对复合材料板振动屈曲的影响，如图7所示。由L_cr-N_cr曲线可知：在敏感区，不同的铺设角度对临界荷载N_cr的影响较大，且随应力波的传播不断减小，到达非敏感区后趋于稳定。图7表明，铺层角度小的单层板的层数越多，板的临界屈曲荷载越大，说明复合材料板的铺设角度直接影响板的屈曲临界荷载。

图 7 不同铺层角度条件下N_cr与L_cr的关系曲线

Figure 7. Relationship between N_cr and L_cr under different laying angles

下载: 全尺寸图片幻灯片

将E组数据代入式(25)中，得到不同初相位对复合材料板振动屈曲的影响，如图8所示。由L_cr-N_cr曲线可知：振型函数的初相位越大，对应的临界荷载越大。在敏感区，振型函数的初相位对临界荷载N_cr的影响较小，且随应力波的传播不断减小；到达非敏感区之后，影响趋于平缓。图8表明，振型函数的初相位越大，板的屈曲临界荷载越大。

图 8 不同初相位条件下N_cr与L_cr的关系曲线

Figure 8. Relationship between N_cr and L_cr under the condition of initial phase of different mode functions

下载: 全尺寸图片幻灯片

将F组数据代入式(25)中，得到不同铺层层数对复合材料板振动屈曲的影响，如图9所示。由L_cr-N_cr曲线可知：当按照不同层数铺设时，敏感区的临界荷载N_cr的变化较大，且随应力波的传播不断减小；到达非敏感区之后变化较小并趋于平缓。图9表明，对于厚度固定且对称铺设的板，当铺设层数达到7时，其屈曲荷载随层数增加趋于稳定。

图 9 不同铺层层数下N_cr与L_cr的关系曲线

Figure 9. Relationship between N_cr and L_cr under different laying modes

下载: 全尺寸图片幻灯片

将G组数据代入式(25)中，得到不同铺设厚度对复合材料板振动屈曲的影响，如图10所示。由L_cr-N_cr曲线可知：板的铺设厚度越大，对应的临界荷载越大。在敏感区，不同的板厚对临界荷载N_cr的影响很大，且随应力波的传播不断减小，到达非敏感区后趋于稳定。图10表明，复合材料板的铺设厚度将直接决定板的屈曲临界荷载。

图 10 不同铺设厚度下N_cr与L_cr的关系曲线

Figure 10. Relationship between N_cr and L_cr under different thicknesses

下载: 全尺寸图片幻灯片

3. 结　论

基于Kirchhoff薄板理论和Hamilton变分原理，建立了具有初始几何缺陷的四边简支复合材料板的振动控制方程。采用伽辽金法，选取符合边界条件的振型函数求解控制方程，得到屈曲临界载荷表达式。数值计算结果表明：应力波在未发生反射前的传播过程中，复合材料板的振动屈曲临界载荷随着临界长度的增大、铺设厚度的减小、初始几何缺陷系数的增大、振型函数初相位的减小而减小；复合材料板的各层铺层角度与荷载作用方向的夹角越小，屈曲临界载荷越大，当对称铺设层数达7层时，临界荷载趋于稳定。研究结果可为工程中复合材料板的结构设计与应用提供一定的参考。

图 1 3种木贼属植物：(a)沼泽问荆，(b)野生问荆，(c)斑纹木贼^[11]

Figure 1. Three horsetails: (a) Marsh horsetail, (b) field horsetail and (c) variegated horsetail^[11]

下载: 全尺寸图片幻灯片

图 2 模型的载荷条件(a)和横截面(b)

Figure 2. Load condition (a) and cross section (b) of the model

下载: 全尺寸图片幻灯片

图 3 网格收敛结果：(a)载荷-位移曲线, (b)吸能量和计算时间

Figure 3. Mesh convergence results: (a) load-displacement curves, (b) energy absorption and computation time

下载: 全尺寸图片幻灯片

图 4 实验^[14]与有限元模拟的比较：(a) 变形模态，(b) 载荷-位移曲线

Figure 4. Comparison between experiment^[14] and finite element simulation: (a) deformation mode, (b) load-displacement curves

下载: 全尺寸图片幻灯片

图 5 n、d和t对 $S_{ \text{EA}}$ 的影响

Figure 5. Effects of n, d and t on $S_{ \text{EA}}$

下载: 全尺寸图片幻灯片

图 7 n、d和t对 $C_{\text{FE}}$ 的影响

Figure 7. Effects of n, d and t on $C_{\text{FE}}$

下载: 全尺寸图片幻灯片

图 6 n、d和t对 $F_{\text{max}}$ 的影响

Figure 6. Effects of n, d and t on $F_{\text{max}}$

下载: 全尺寸图片幻灯片

图 8 不同壁厚下HBTS的变形模态

Figure 8. Deformation modes of HBTS with different wall thicknesses

下载: 全尺寸图片幻灯片

图 9 多目标优化流程图

Figure 9. Multi-objective optimization flow chart

下载: 全尺寸图片幻灯片

图 10 S_EA和F_max的Pareto前沿

Figure 10. Pareto front of S_EA and F_max

下载: 全尺寸图片幻灯片

图 11 Pareto前沿面上各HBTS的参数分布：(a) 壁厚分布，(b) 内径和肋骨数分布

Figure 11. Parameter distribution of each structure on the Pareto front：(a) distribution of wall thickness, (b) distributions of inner diameter and number of ribs

下载: 全尺寸图片幻灯片

图 12 不同 ${{F}}_{\text{max}}$ 限制下最优HBTS的载荷-位移曲线

Figure 12. Load-displacement curves of optimal HBTSs with different ${{F}}_{\text{max}}$ limits

下载: 全尺寸图片幻灯片

表 1 铝合金AA6061的材料参数

Table 1. Parameters of aluminum alloy AA6061

Density/(g∙cm⁻³)	Poisson’s ratio	Elastic modulus/GPa	Yield stress/MPa	Ultimate strength/MPa
2.70	0.3	70.0	137.03	223.53

下载: 导出CSV

表 2 HBTS的耐撞性比较

Table 2. Comparison of crashworthiness of HBTS

Structure	${ {t} }{_{\text{R} }}$ /mm	${ {t} }{_{\text{I} }}$ /mm	${ {t} }{_{\text{r} }}$ /mm	M/g	${ {E} }{_{\text{A} }}$ /J	${ {S} }{_{ \text{EA} } }$ /(J∙g⁻¹)	${ {F} }{_{\text{max} }}$ /kN	${ {C} }{_{\text{FE} }}$ /%
a	0.800	0.100	0.100	110.2	70.16	0.637	1.841	82.8
b	0.158	0.800	0.158	110.2	28.02	0.254	1.125	54.1
c	0.271	0.271	0.800	110.2	49.11	0.446	1.401	76.2
d	1.500	0.800	0.800	308.1	361.61	1.174	9.095	86.4
e	0.858	1.500	0.858	308.1	390.85	1.268	9.691	87.7
f	0.971	0.971	1.500	308.1	315.82	1.025	8.110	84.7

下载: 导出CSV

表 3 壁厚与内径水平

Table 3. Wall thickness and inner diameter level

Level	${ {t} }{_{\mathrm{R} }}$ /mm	${ {t} }{_{\text{I} }}$ /mm	${ {t} }{_{\text{r} }}$ /mm	d/mm
1	0.10	0.10	0.10	10.0
2	0.45	0.45	0.45	17.5
3	0.80	0.80	0.80	25.0
4	1.15	1.15	1.15	32.5
5	1.50	1.50	1.50	40.0

下载: 导出CSV

表 4 Kriging代理模型精度

Table 4. Accuracy assessment of the Kriging surrogate models

${ {S} }{_{ \text{EA} } }$			${ {F} }{_{\text{max} } }$
${ {R} }{^{\text{2} }}$	δ_ARE/%	δ_RAAE/%	${ {R} }{^{\text{2} }}$	δ_ARE/%	δ_RAAE/%
0.989	1.52	1.33	0.996	1.68	0.67

下载: 导出CSV

表 5 不同 $F_{\text{max}}$ 限制下的最优HBTS

Table 5. Optimal HBTSs with different $F_{\text{max}}$ limits

Type	Design parameters					${ {S} }{_{ \text{EA} } }$ /(J∙g⁻¹)		$\delta_{S_{\rm EA}}$ /%	${ {F} }{_{\text{max} }}$ /(J∙g⁻¹)		$\delta_{F_{\max}}$ /%
Type	${ {t} }{_{\text{R} }}$ /mm	${ {t} }{_{\text{I} } }$ /mm	${ {t} }{_{\text{r} } }$ /mm	n	d/mm	Kriging	FE	$\delta_{S_{\rm EA}}$ /%	Kriging	FE	$\delta_{F_{\max}}$ /%
${ {F} }{_{\text{max} }}$ <5 kN	0.759	0.413	0.775	16	13.9	1.054	1.026	2.73	4.994	4.953	0.83
${ {F} }{_{\text{max} } }$ <10 kN	0.954	1.221	1.064	16	31.9	1.347	1.348	−0.07	9.952	10.111	−1.57

下载: 导出CSV

参考文献(20)

[1]	武和全, 毛鸿锋, 侯海彬. 复合材料仿竹薄壁管耐撞性和可靠性研究 [J]. 南京理工大学学报, 2017, 41(2): 186–190, 197. WU H Q, MAO H F, HOU H B. Study on crashworthiness and reliability of composite bamboo-like thin-walled tube [J]. Journal of Nanjing University of Science and Technology, 2017, 41(2): 186–190, 197.
[2]	ZHAO Z, HUANG W, LI B, et al. Synergistic effects of chiral morphology and reconfiguration in cattail leaves [J]. Journal of Bionic Engineering, 2015, 12(4): 634–642. doi: 10.1016/S1672-6529(14)60153-0
[3]	LIU Q, MA J, HE Z, et al. Energy absorption of bio-inspired multi-cell CFRP and aluminum square tubes [J]. Composites Part B: Engineering, 2017, 121: 134–144. doi: 10.1016/j.compositesb.2017.03.034
[4]	范晓文, 杨欣, 许述财, 等. 仿骨单位薄壁结构轴向和斜向耐撞性研究 [J]. 载人航天, 2020, 26(2): 142–151. doi: 10.3969/j.issn.1674-5825.2020.02.002 FAN X W, YANG X, XU S C, et al. Study on crashworthiness of thin-walled structure based on osteon under axial and oblique loads [J]. Manned Spaceflight, 2020, 26(2): 142–151. doi: 10.3969/j.issn.1674-5825.2020.02.002
[5]	TSANG H H, RAZA S. Impact energy absorption of bio-inspired tubular sections with structural hierarchy [J]. Composite Structures, 2018, 195: 199–210. doi: 10.1016/j.compstruct.2018.04.057
[6]	黄晗, 许述财, 杜雯菁, 等. 基于虾螯结构的仿生薄壁管吸能特性分析及优化 [J]. 北京理工大学学报, 2020, 40(3): 267–274. HUANG H, XU S C, DU W J, et al. Energy absorption analysis and optimization of a bionic thin-walled tube based on shrimp chela [J]. Transactions of Beijing Institute of Technology, 2020, 40(3): 267–274.
[7]	于鹏山, 刘志芳, 李世强. 新型仿竹薄壁圆管的设计与吸能特性分析 [J]. 高压物理学报, 2021, 35(5): 054205. doi: 10.11858/gywlxb.20210710 YU P S, LIU Z F, LI S Q. Design and energy absorption characteristic analysis of a new bio-bamboo thin-walled circular tube [J]. Chinese Journal of High Pressure Physics, 2021, 35(5): 054205. doi: 10.11858/gywlxb.20210710
[8]	HUANG J L, DENG X L, LIU W Y. Bionic design of bend-twist coupled thin-walled beam based on the structure of rice stem [J]. Mechanics of Advanced Materials and Structures, 2021: 1–14.
[9]	霍鹏, 许述财, 范晓文, 等. 鹿角骨单位仿生薄壁管斜向冲击耐撞性研究 [J]. 爆炸与冲击, 2020, 40(11): 124–135. HUO P, XU S C, FAN X W, et al. Oblique impact resistance of a bionic thin-walled tube based on antles osteon [J]. Explosion and Shock Waves, 2020, 40(11): 124–135.
[10]	白芳华, 张林伟, 白中浩, 等. 基于甲虫鞘翅的客车八边形仿生多胞薄壁管耐撞性研究 [J]. 振动与冲击, 2019, 38(21): 24–30. BAI F H, ZHANG L W, BAI Z H, et al. Crashworthiness of coach’s octagonal bionic mult-cell thin-walled tubes based on beetle elytra [J]. Journal of Vibration and Shock, 2019, 38(21): 24–30.
[11]	YIN H F, XIAO Y Y, WEN G L, et al. Crushing analysis and multi-objective optimization design for bionic thin-walled structure [J]. Materials & Design, 2015, 87: 825–834.
[12]	XIAO Y Y, YIN H F, FANG H B, et al. Crashworthiness design of horsetail-bionic thin-walled structures under axial dynamic loading [J]. International Journal of Mechanics and Materials in Design, 2016, 12(4): 563–576. doi: 10.1007/s10999-016-9341-6
[13]	YIN H F, XIAO Y Y, WEN G L, et al. Multi-objective robust optimization of foam-filled bionic thin-walled structures [J]. Thin-Walled Structures, 2016, 109: 332–343. doi: 10.1016/j.tws.2016.10.011
[14]	SUN G Y, TIAN X Y, FANG J G, et al. Dynamical bending analysis and optimization design for functionally graded thickness (FGT) tube [J]. International Journal of Impact Engineering, 2015, 78: 128–137. doi: 10.1016/j.ijimpeng.2014.12.007
[15]	尹华伟, 王陈凌, 段金曦, 等. 新型薄壁管耐撞性分析及优化设计 [J]. 高压物理学报, 2021, 35(3): 034202. doi: 10.11858/gywlxb.20200624 YIN H W, WANG C L, DUAN J X, et al. Crashworthiness analysis and optimization design of new thin-walled tube [J]. Chinese Journal of High Pressure Physics, 2021, 35(3): 034202. doi: 10.11858/gywlxb.20200624
[16]	郝文乾, 卢进帅, 黄睿, 等. 轴向冲击载荷下薄壁折纹管的屈曲模态与吸能 [J]. 爆炸与冲击, 2015, 35(3): 380–385. doi: 10.11883/1001-1455-(2015)03-0380-06 HAO W Q, LU J S, HUANG R, et al. Buckling and energy absorption properties of thin-walled corrugated tubes under axial impacting [J]. Explosion and Shock Waves, 2015, 35(3): 380–385. doi: 10.11883/1001-1455-(2015)03-0380-06
[17]	MEGUID S A, ATTIA M S, MONFORT A. On the crush behaviour of ultralight foam-filled structures [J]. Materials & Design, 2004, 25(3): 183–189.
[18]	穆雪峰, 姚卫星, 余雄庆, 等. 多学科设计优化中常用代理模型的研究 [J]. 计算力学学报, 2005, 22(5): 608–612. MU X F, YAO W X, YU X Q, et al. A survey of surrogate models used in MDO [J]. Chinese Journal of Computational Mechanics, 2005, 22(5): 608–612.
[19]	FANG J G, GAO Y K, SUN G Y, et al. Dynamic crashing behavior of new extrudable multi-cell tubes with a functionally graded thickness [J]. International Journal of Mechanical Sciences, 2015, 103: 63–73. doi: 10.1016/j.ijmecsci.2015.08.029
[20]	FU J, LIU Q, LIU-FU K, et al. Design of bionic-bamboo thin-walled structures for energy absorption [J]. Thin-Walled Structures, 2019, 135: 400–413. doi: 10.1016/j.tws.2018.10.003