2. 中国舰船研究设计中心,湖北 武汉 430064
2. China Ship Development and Design Center, Wuhan 430064, China
结构动力学试验是获取结构振动性能并对其进行减振降噪设计的重要手段。对于多舱段这类大型复杂结构,结构整体试验对试验条件、试验技术要求很高[1]。由于多舱段结构的振动源主要集中在某些舱段,因此可以将其截断出来单独进行舱段试验以降低试验成本,并通过修正相邻舱段的影响,得到截断舱段在多舱段中的实际振动特性。
王路才[2]通过设置不同的声反射边界条件来研究水下结构舱段以外部分对舱段自身辐射声反射的影响,探讨了以舱段模型代替整艇模型进行噪声估算的可行性。丁宏[3]通过对含环肋圆柱壳的单舱段一端分别施加自由、周向径向简支约束、轴向约束以及三向约束等不同的边界条件,从模态和频响曲线两方面与全船模型的舱段进行比较。庞福振[4]提出一种利用局部舱段代替整船模型进行船舶结构噪声数值预报的模型简化方法,可先由行波法要求确定船舶舱段的最小截断尺度,再结合船舶舱段结构划分的实际,最终确定船舶舱段的截断尺度。王璟旻[5]提出舰船抗冲击分析的截断模型的2种动力学等效边界条件,其中集中刚度替代法要优于基于MPC法的梁-舱-梁有限元方法。
前面的研究多是对结构截断边界的近似替代或者忽略截断边界的影响,其精度仍有待提高。本文充分考虑相邻舱段的刚度和质量对截断舱段的影响,通过动态缩聚法将相邻舱段的阻抗矩阵缩聚到其与截断舱段的连接处,并将缩聚阻抗矩阵作为修正项加到截断舱段的动力学方程中,实现对截断舱段的精确修正。此外,还研究了利用自由界面模态综合法,将截断舱段与相邻舱段的模态结果综合起来,得到整体模型的动力学参数。
1 截断舱段的修正从多舱段艇体结构中截取出的独立舱段缺少相邻舱段的约束,因此其与多舱段的动力学特性有很大的区别。即使考虑舱段在截断处的相互作用力,采用一定的边界约束条件来修正单舱段,也不能完全替代多舱段结构的连续性条件[3]。而充分考虑边缘舱段的刚度K、质量M和阻尼C对中间舱段的影响,并利用动态缩减法[6]将边缘舱段的阻抗矩阵缩聚到其与中间舱段的连接处,再将缩聚矩阵与中间舱段进行耦合,可实现对中间舱段动力学特性的精确修正。
1.1 多级子结构法当边缘舱段的自由度数量较多时,若直接进行缩聚,会造成庞大的计算量而降低求解的速度。子结构法是解决这类问题的一种有效方法,其将一个大型的复杂结构划分为若干子结构,先确定各子结构的特性,然后再将子结构装配成整体结构,最后确定整体结构的特性。子结构与子结构之间的连接节点(或外部节点)就是子结构的主节点,对应的自由度就是主自由度。由于舱段含有多种周期性的构件,因此可以使用多级子结构法进一步减小部件的自由度来进行求解。
1.1.1 子结构划分对于含有几种典型构件的多舱段结构,可以有如下的划分:
1)一级子结构划分。如图1所示,根据舱壁位置,将多舱段结构离散为多个单舱段(中间舱段含有舱壁),得到一级子结构。
![]() |
图 1 一级子结构的划分示意图 Fig. 1 The schematic diagram of division of first-order substructures |
2)二级子结构划分。根据环肋位置,将各个舱段离散为多个圆柱壳子段和环肋,得到二级子结构。这里采用有限元模型来表示二级子结构。边缘舱段的二级子结构的划分如图2所示,包含圆柱壳子段和环肋。
![]() |
图 2 边缘舱段的二级子结构划分示意图 Fig. 2 The schematic diagram of division of second-order substructures in edge segment |
对于包含基座的中间舱段,为了后面的分析计算,这里仅将舱壁和基座与加筋圆柱壳体分离,其二级子结构的划分如图3所示,包含加筋圆柱壳、基座和舱壁。
![]() |
图 3 中间舱段的二级子结构划分示意图 Fig. 3 The schematic diagram of division of second-order substructures in central segment |
前面将多舱段结构逐级划分为二级子结构,下面通过2次缩聚与耦合,将多舱段结构的整体自由度缩聚到中间舱段上。
1)第1次缩聚。将二级子结构(圆柱壳子段、环肋、基座和舱壁),缩聚到与其他子结构(一级或二级)的连接处,得到各二级子结构的缩聚阻抗矩阵。如图4所示,红点是保留的节点。值得注意的是,中间舱段的加筋圆柱壳部分不进行缩聚(即保留全部节点自由度),而基座面板上还保留了载荷作用点。
![]() |
图 4 二级子结构的缩聚示意图 Fig. 4 The schematic diagram of condensation of second-order substructures |
2)第1次矩阵拼装。将一级子结构下属的二级子结构部件的缩聚阻抗矩阵按照对应位置拼装在一起,得到一级子结构的初步缩聚阻抗矩阵
![]() |
图 5 二级子结构的组合示意图 Fig. 5 The schematic diagram of combination of second order substructures |
$ {{\boldsymbol{Z}}_I} = \left[ {\begin{array}{*{20}{c}} {Z_{11}^{{h_1}}}&{Z_{12}^{{h_1}}}&{}&{}&{}&{} \\ {Z_{21}^{{h_1}}}&{Z_{22}^{{h_1}} + Z_{11}^{{h_2}} + {Z^{{r_1}}}}&{Z_{12}^{{h_1}}}&{}&{}&{} \\ {}&{Z_{21}^{{h_2}}}&{Z_{22}^{{h_2}} + Z_{11}^{{h_3}} + {Z^{{r_2}}}}&{}&{}&{} \\ {}&{}&{}& \ddots &{}&{} \\ {}&{}&{}&{}&{Z_{22}^{{h_{n - 1}}} + Z_{11}^{{h_n}} + {Z^{{r_{n - 1}}}}}&{Z_{12}^{{h_n}}} \\ {}&{}&{}&{}&{Z_{21}^{{h_n}}}&{Z_{11}^{{h_n}} + {Z^b}} \end{array}} \right]。$ | (1) |
式中:
3)第2次缩聚。将组装的一级子结构的初步缩聚阻抗矩阵,进一步缩聚到与其他一级子结构(舱壁或舱段)的连接处。
4)第2次矩阵拼装。将各一级子结构的二次缩聚矩阵拼装在一起。最终得到只含有中间舱段的加筋圆柱壳部分的全部节点自由度和基座上载荷作用点自由度的阻抗矩阵。
1.2 算例验证如图6所示的多舱段有限元模型,壳体轴向划分120份,周向划分96份,总的自由度数为16351×6。该多舱段的模型是在中间主舱段的左右各延伸一个与其等长度的舱段,并按照相同的肋距设置环肋。左右边缘舱段不含内部基座,其外壳厚度和环肋参数均与主舱段相同。具体尺寸:单舱长L=2 m,直径D=1.5 m,肋距lr=0.4 m,肋骨尺寸0.008 m×0.02 m,壳体厚度t1=0.006 m,舱壁厚度t2=0.03 m。基座位于主舱段底部中心位置,其长宽高为0.6 m×0.388 m×0.3 m,基座面板厚度t3=0.015 m、腹板厚度t4=0.010 m。所有构件均由钢制成,材料参数为:杨氏模量E=210 GPa,泊松比v=0.3,密度ρ=7 850 kg/m3,损耗因子η=0.01。激励力为4个机脚处的竖直向下的单位力,分析频段为1~400 Hz。柱坐标系的原点设置在中间主舱段的壳体中心处,取壳体表面上的测点为:P1(D/2,180,0.8),P2(D/2,180,0)。在Ansys采用Shell181单元分别建立各二级子结构(圆柱壳子段、环肋、舱壁等)的有限元模型,并提取质量矩阵和刚度矩阵进行多次缩聚和组合,最后进行谐响应计算。
![]() |
图 6 多舱段的有限元模型 Fig. 6 The finite element model of the multiple cabin segments |
为了验证多级子结构方法的正确性,在Ansys中直接建立整体三舱段的有限元模型,并完成谐响应计算(下称“直接法”)。2种方法得到的测点位移响应结果如图7所示。由图可知,子结构法与完全法的频响曲线完全吻合,验证了子结构法正确性。
![]() |
图 7 舱段的测点速度响应 Fig. 7 The velocity responses of the cabin segments |
对于本节计算的舱段模型,通过划分为两级子结构,计算机处理的最大矩阵维度为5184×5184,完成单个频点的两级缩聚和组合耗时30 s左右。如果直接对边缘舱段整体进行缩聚,需要处理的最大矩阵维度为25920×25920,完成单个频点的缩聚耗时180 s左右。因此,多级子结构法极大地加快了计算分析的效率。
2 子舱段模态综合前面是在物理坐标下,通过采用多级子结构法将相邻舱段的阻抗矩阵缩聚到与截断舱段的连接处,从而实现对截断舱段的动力学特性修正。三舱段结构的模态比单舱段结构丰富,其边缘舱段之间存在复杂的耦合作用,三舱段模型存在2个船体梁模态,这些模态在单舱段模型中不会出现,因此截断舱段不能完全替代多舱段结构整体的振动特性。而模态综合法[7]相比于前面的动态缩聚求解方法,不仅可以大幅减少系统的自由度,有效提高计算效率,还可以获得多舱段结构整体的振动特性。相比于固定界面模态综合法,自由界面模态综合法更接近于截断舱段试验时两边自由的状态,便于将试验模型与理论模型结合起来。
2.1 自由界面模态综合法将系统划分为若干个子结构,子结构的无阻尼动力学方程为:
$ {\boldsymbol{M}}\ddot x + {\boldsymbol{K}}x = f ,$ | (2) |
将x分为内部自由度
$ {x}_{j}={\boldsymbol B}_{j}x\text{,}f={\boldsymbol B}_{j}^{\rm T}{f}_{j} 。$ | (3) |
式中:
在简谐力激励下,子结构的响应为:
$ x\left( \omega \right) = {\boldsymbol{G}}\left( \omega \right)f = {\boldsymbol{G}}\left( \omega \right){\boldsymbol{B}}_j^{\rm{T}}{f_j} ,$ | (4) |
式中,G(ω)为结构的动柔度矩阵。
取子结构的前k阶模态作为保留主模态
$ \begin{split} & {\boldsymbol{G}}\left( \omega \right) = {{\boldsymbol{G}}_k}\left( \omega \right) + {{\boldsymbol{G}}_h}\left( \omega \right) ,\\ & {{\boldsymbol{G}}_k}\left( \omega \right) = {\varPhi _k}{\left( {{\varLambda _k} - {\omega ^2}{I_k}} \right)^{ - 1}}\varPhi _k^{\rm{T}} ,\\ & {{\boldsymbol{G}}_h}\left( \omega \right) = {\varPhi _h}{\left( {{\varLambda _h} - {\omega ^2}{{\boldsymbol{I}}_h}} \right)^{ - 1}}\varPhi _h^{\rm{T}} 。\end{split} $ | (5) |
式中:
定义剩余附着模态为:
$ {\varPsi _s} = {{\boldsymbol{G}}_h}\left( \omega \right){\boldsymbol{B}}_j^{\rm{T}} = {\boldsymbol{G}}\left( \omega \right){\boldsymbol{B}}_j^{\rm{T}} - {\varPhi _k}{\left( {{\varLambda _k} - {\omega ^2}{{\boldsymbol{I}}_k}} \right)^{ - 1}}\varPhi _k^{\rm{T}}{\boldsymbol{B}}_j^{\rm{T}} ,$ | (6) |
由于高阶剩余模态Фh是未知的,将柔度矩阵在ω0处展开,可得结构的剩余动柔度为:
$ {{\boldsymbol{G}}_h}\left( {{\omega _0}} \right) = {\boldsymbol{G}}\left( {{\omega _0}} \right) - {\varPhi _k}{\left( {{\varLambda _k} - \omega _0^2{{\boldsymbol{I}}_k}} \right)^{ - 1}}\varPhi _k^{\rm{T}} 。$ | (7) |
式(7)将ω0处的柔度矩阵作为所有频率下的柔度矩阵,因此也称准剩余动柔度。于是子结构的剩余附着模态表达式为:
$ {\varPsi _s} = {\boldsymbol{G}}\left( {{\omega _0}} \right){\boldsymbol{B}}_j^{\rm{T}} - {\varPhi _k}{\left( {{\varLambda _k} - \omega _0^2{{\boldsymbol{I}}_k}} \right)^{ - 1}}\varPhi _k^{\rm{T}}{\boldsymbol{B}}_j^{\rm{T}}。$ | (8) |
式(8)考虑了
将保留主模态
$ \varPhi = \left[ {\begin{array}{*{20}{c}} {{\varPhi _k}}&{{\varPsi _s}} \end{array}} \right],$ | (9) |
于是,可以将子结构的物理坐标x用模态坐标p表示为:
$ x = \left[ {\begin{array}{*{20}{c}} {{\varPhi _k}}&{{\varPsi _s}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{p_k}} \\ {{f_j}} \end{array}} \right\} = \varPhi \cdot p 。$ | (10) |
式中:
将式(10)代入式(2)可得,子结构在模态坐标下的动力学方程为:
$ \bar M\ddot p + \bar Kp = g 。$ | (11) |
式中:
如图8所示,可以将三舱段结构截断离散为3个子舱段结构,即中间含内部基座的舱段(β)、两边的舱段(α)和(γ)。
![]() |
图 8 三舱段的子结构示意图 Fig. 8 The schematic diagram of the substructures of the three cabin segments |
将3个子结构联立,得到整个系统的自由振动方程为:
$\begin{split} &\left[ {\begin{array}{*{20}{c}} {{}^\alpha \bar M}&{}&{} \\ {}&{{}^\gamma \bar M}&{} \\ {}&{}&{{}^\beta \bar M} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{}^\alpha \ddot p} \\ {{}^\gamma \ddot p} \\ {{}^\beta \ddot p} \end{array}} \right\} + \left[ {\begin{array}{*{20}{c}} {{}^\alpha \bar K}&{}&{} \\ {}&{{}^\gamma \bar K}&{} \\ {}&{}&{{}^\beta \bar K} \end{array}} \right]\times\\ &\qquad \left\{ {\begin{array}{*{20}{c}} {{}^\alpha p} \\ {{}^\gamma p} \\ {{}^\beta p} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{}^\alpha g} \\ {{}^\gamma g} \\ {{}^\beta g} \end{array}} \right\} 。\end{split}$ | (12) |
根据界面力平衡条件
$\begin{split} \left\{ {\begin{array}{*{20}{c}} {{}^\alpha p} \\ {{}^\gamma p} \\ {{}^\beta p} \end{array}} \right\} = &{\left[ {\begin{array}{*{20}{c}} I&0&0&0&0&0&0 \\ 0&0&I&0&0&0&0 \\ 0&0&0&0&I&0&0 \\ 0&{ - I}&0&0&0&I&0 \\ 0&0&0&{ - I}&0&0&I \end{array}} \right]^{\text{T}}} \times\\ &\left[ {\begin{array}{*{20}{c}} I \\ {L_1^{ - 1}{L_2}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{}^\alpha {p_k}} \\ {{}^\gamma {p_k}} \\ {{}^\beta {p_k}} \end{array}} \right\} = T \cdot q 。\end{split}$ | (13) |
式中:
用式(13)对式(12)作坐标转换,并且由于整体不受外力,可以证明式(12)的右端项经过坐标转换后为零向量,所以有
$ \tilde {\boldsymbol{M}}\ddot q + \tilde {\boldsymbol{K}}q = 0 。$ | (14) |
式中:
$ \tilde {\boldsymbol{M}} = {{\boldsymbol{T}}^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} {{}^\alpha \bar M}&{}&{} \\ {}&{{}^\gamma \bar M}&{} \\ {}&{}&{{}^\beta \bar M} \end{array}} \right]{\boldsymbol{T}},\tilde {\boldsymbol{K}} = {{\boldsymbol{T}}^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} {{}^\alpha \bar K}&{}&{} \\ {}&{{}^\gamma \bar K}&{} \\ {}&{}&{{}^\beta \bar K} \end{array}} \right]{\rm{T}}。$ | (15) |
求解式(14)得到整个系统的固有频率和在广义坐标q下的振型,根据式(13)和式(10),可将结果返回到物理坐标上。例如,对于子结构α,有
$ {}^\alpha x = {}^\alpha \Phi {}^\alpha {\boldsymbol{BT}}q 。$ | (16) |
式中,
前面的研究是基于有限元建模的理论舱段,而在工程中试验舱段的质量矩阵和刚度矩阵是未知的,因此需要建立试验舱段与理论舱段的模态综合方法。考虑子结构的模态质量和模态刚度,即
$ \bar{\boldsymbol M}=\left[\begin{array}{ll} \bar{M}_{k} & \\ & \bar{M}_{s} \end{array}\right], \quad \bar{\boldsymbol K}=\left[\begin{array}{ll} K_{k} & \\ & \bar{K}_{s} \end{array}\right] 。$ | (17) |
式中:
$ \bar{\boldsymbol M}_{z}=\varPsi_{z}^{\rm T} {\boldsymbol{M}} \varPsi_{z},$ | (18) |
$ {\boldsymbol{R}}_{t}=\varPsi_{t}^{\rm T} {\boldsymbol{K}} \varPsi_{t} 。$ | (19) |
求解式(18)和式(19)需要知道子结构的质量矩阵M和刚度矩阵K。由于目前还没有可靠的方法通过试验手段获得结构的刚度矩阵和质量矩阵,因此需要对试验舱段进行特殊的处理[8]。
将式(6)代入式(19),可得
$ {\boldsymbol{R}}_{f}={\boldsymbol{B}}_{j} \varPhi_{h}\left(\varLambda_{h}-\omega^{2} {\boldsymbol{I}}_{A}\right)^{-1} \varLambda_{A}\left(\varLambda_{A}-\omega^{2} {\boldsymbol{I}}_{A}\right) \varPhi_{h}^{\rm T} {\boldsymbol{B}}_{j}^{\rm T},$ | (20) |
若采用ω0处的准剩余动柔度,由于ω0一般远小于截断频率,则
$ {\boldsymbol{R}}_{s} \approx {\boldsymbol{B}}_{j} \varPhi_{A}\left(\varLambda_{A}-\omega_{0}^{2} {\boldsymbol{I}}_{A}\right) \varPhi_{A}^{\rm T} {\boldsymbol{B}}_{j}^{\rm T}={\boldsymbol{B}}_{j} \varPsi_{s} 。$ | (21) |
将式(6)代入式(12),可推导出
$ \bar{\boldsymbol M}_{s}=0 。$ | (22) |
上述处理避免了刚度阵和质量阵的识别,适用于试验舱段参与模态综合。
2.3 数值验证为了验证子结构模态综合法计算整体模态的正确性,取1.2节的算例模型,在Ansys中直接对三舱段的整体模型进行模态计算(下称“直接法”),与模态综合法的计算结果进行对比。
首先通过理论舱段之间的模态综合,研究不同自由度数的界面力对模态综合结果的影响。表1给出了模态综合法与直接法计算的固有频率对比。从表中的结果可知,当界面力取完备自由度(六自由度)时,子结构模态综合法计算的固有频率与整体模型吻合良好。若减少界面力的自由度,忽略转动自由度对应的力矩,当仅保留x,y和z三个方向的力时,二者的计算结果仍能保持良好的吻合;而当界面力自由度只剩x和z两个方向时,模态综合法计算的固有频率误差增大,且开始出现部分模态的缺失。经检验,至少需要考虑x,y和z三个方向的界面力,才能保证模态综合结果的正确性。
![]() |
表 1 模态综合法与直接法计算的多舱段固有频率对比 Tab.1 Comparison of natural frequencies calculated by CMS method and Direct method |
进一步假设中间主舱段为试验舱段,参与模态综合时仅提供模态参数,而质量矩阵和刚度矩阵未知。左、右舱段仍按有限元模型参与模态综合,其质量矩阵和刚度矩阵已知。界面力取x,y和z三个方向的力。为进一步考察模态综合的准确性,基于综合的模态参数对结构作谐响应计算。激励力为垂直于主舱段的基座面板中心向上的单位简谐力,分析频段为1~300 Hz。设置A(D/2,0,0)和B(D/2,90,0) 两个测点,柱坐标原点在主舱段的几何中心位置。计算结果均与Ansys整体模型计算结果(下称“直接法”)进行对比。2种方法的计算结果吻合良好(见表1)。图9给出了模态综合法和直接法计算的振动响应对比。从图中可知,2种方法计算的测点响应结果吻合良好。结果表明,子结构模态综合法在实现试验舱段与理论舱段综合方面具有很强的适用性和可靠性。
![]() |
图 9 模态综合法和直接法计算的舱段振动响应对比 Fig. 9 Comparison of the velocity responses of the cabin segments calculated by CMS method and Direct method |
对于模态频率的截断,各子舱段的截断频率要保持一致或者接近,并不低于分析频率,即可保证足够精度的结果。对于本节算例,分析频率为1~300 Hz,主舱段截取到304.5 Hz,左、右舱段均截取到302.5 Hz。经过验证,若分析频率超过截断频率,超过部分的综合结果将不正确。
3 结 语本文通过研究子结构法在陆上截断舱段振动特性修正中的应用,得到以下结论:
1)通过动态缩聚方法得到边缘舱段在连接处的阻抗,可实现对截断舱段的精确修正。进一步利用多级子结构法,将舱段按实际需要划分为若干级子结构,可大大提高动态缩聚的速度。
2)采用自由界面模态综合法,将各子舱段的模态结果综合起来,可得到整体模型的动力学参数。至少需要考虑x,y和z三个自由度方向的力,才能保证模态综合结果的准确性。
3)通过构造剩余刚度和剩余质量,基于数值模拟计算,实现了试验舱段和理论舱段之间的模态综合。通过与整体舱段的结果进行对比,验证了方法的适用性和准确性。
进一步地,可以将多级子结构法和模态综合法结合起来使用。首先,利用多级子结构法对理论舱段进行自由度缩减,然后再与试验舱段进行模态综合,以实现对陆上截断舱段振动特性的快速修正。
[1] |
朴思扬, 祁峰, 张亚辉. 基于模态综合法的大型结构模型修正技术[J]. 应用数学和力学, 2018, 39(9): 989-998. |
[2] |
王路才, 周其斗, 纪刚, 等. 以舱段模型代替整艇模型进行噪声估算的可行性探讨[J]. 中国舰船研究, 2010, 5(6): 26-32. DOI:10.3969/j.issn.1673-3185.2010.06.006 |
[3] |
丁宏, 陈美霞, 魏建辉, 等. 舱段截断时边界条件的选取方法[J]. 舰船科学技术, 2010, 5(6): 26-32. DING H, CHEN M X, WEI J H, et al. Selection method of boundary conditions in the cabin segment truncation[J]. Ship Science and Technology, 2010, 5(6): 26-32. |
[4] |
庞福振. 船舶结构噪声截断模型数值预报方法研究[D]. 哈尔滨: 哈尔滨工程大学, 2012.
|
[5] |
王璟旻. 船舶结构抗水下爆炸分析简化模型研究[D]. 哈尔滨: 哈尔滨工程大学, 2017.
|
[6] |
GORDIS J H. An analysis of the improved reduced system (IRS) model reduction procedure[C]. Proceedings of the 10th International Modal Analysis Conference, 1992: 471−479.
|
[7] |
温争鸣. 自由界面模态综合法的研究与系统实现[D]. 武汉: 华中科技大学, 2014.
|
[8] |
董兴建, 孟光. 剩余质量阵和剩余刚度阵的近似计算方法[J]. 噪声与振动控制, 2009, 29(6): 11-14. DONG X J, MENG G. An approximated method for calculating residual mass matrix and residual stiffness matrix[J]. Noise and Vibration Control, 2009, 29(6): 11-14. |
[9] |
MACNEAL R H. A hybrid method of component mode synthesis[J]. Computers & Structures, 1971, 1(4): 581-601. |
[10] |
KIM J, BOO S, LEE P, et al. A dynamic condensation method with free interface substructuring[J]. Mechanical Systems and Signal Processing, 2019, 218-234. |
[11] |
WANG T, HE H, YAN W, et al. A model-updating approach based on the component mode synthesis method and perturbation analysis[J]. Journal of Sound and Vibration, 2018, 433: 349-365. DOI:10.1016/j.jsv.2018.07.026 |