2. 中国航发沈阳发动机研究所, 沈阳 110015
2. AECC Shenyang Engine Research Institute, Shenyang 110015, China
螺栓连接是一种典型的机械连接固定形式,由于具有可以承受较大载荷并且可重复装配拆卸的优点,被广泛应用于各种机械结构中,如机床、航天器、钢架桥等工程结构的主体受力部分均包含螺栓结合部。已有的研究表明:结合部对于机械结构整体的动态性能具有较大影响[1-2],因此在对结构进行动力学建模时不能简单地把结合部等效于刚性连接,而需要对其进行更准确的建模。然而目前结合部的作用机理尚不明确,因此很难通过解析的方法直接得到准确的结合部动力学特性参数,在实际工程应用中需要对结合部的动力学模型进行简化,再通过有限元仿真、实验和参数辨识方法对结合部的等效动力学参数进行辨识来完成结合部的动力学建模。
目前已有许多结合部动力学建模、参数辨识方法等方面的研究[3]。螺栓结合部常用的建模方法有子结构方法[4-5]、弹簧阻尼法[6]、虚拟材料法[7-8]等,而参数辨识方法通常采用有限元模型与实验模态数据[9]或频响函数[10]相结合,辨识出结合部的动力学参数,最终对比仿真和实验结果加以验证。
另一方面,螺栓连接本身具有不连续、接触、干摩擦等非线性特征,导致螺栓结合部表现出一定的非线性动力学特性。现有的研究大多考虑在特定的螺栓连接状态下如何更准确地进行动力学建模和参数辨识的问题,采用的测试方法也以力锤测试等宽频激励为主,较少考虑螺栓连接非线性对于动力学特性影响。而实际工况中的载荷通常以固定频率载荷为主,如果系统具有较为明显的非线性则会在固定频率载荷的激发下引起等效动力学参数的变化,导致利用线性方法测量和建立的动力学模型与结构实际受载过程中表现出的动力学特性之间存在较大误差。因此,为了建立更加准确的结构动力学模型,有必要研究非线性特性对螺栓结合部动力学特性的影响。
目前已经有一些螺栓结合部非线性动力学特性的研究。陈学前等[11]通过正弦激振实验对螺栓法兰结构的单自由度非线性振动模型的动力学参数进行了辨识。蔡力钢等[12]等结合Tikhonov正则化与迭代算法提高了力状态映射法的辨识精度并通过算例对非线性结合部等效动力学参数进行了辨识。孙伟等[13]等利用力状态映射法辨识了螺栓联接梁结构结合部刚度和阻尼特性参数。以上研究均采用了力状态映射法辨识结合部的非线性动力学参数,但是均是基于特定的螺栓连接状态,没有考虑到螺栓连接状态对非线性特性的影响。螺栓预紧力作为实际螺栓连接结构装配中一个重要的可控因素,对于结合部承载能力、连接刚度阻尼等有着重要影响,另外预紧力的变化也可能引起结合部非线性程度的改变。
为研究螺栓预紧力对结合部动力学的特性,特别是非线性动力学特性的影响情况,本文选择了一种双螺栓结合部试件作为研究对象,通过理论推导对结构在一阶模态附近进行简化并建立了等效动力学模型;然后通过模态实验获得了不同螺栓预紧力下结构的一阶模态参数,得到了线性等效动力学参数随预紧力的变化情况;最后考虑到非线性因素,探究了结合部在激振状态下的动力学特性,并利用力状态映射法对结合部在不同预紧力下非线性等效动力学参数进行了辨识,为进一步研究螺栓结合部的非线性动力学特性以及动力学建模提供了依据。
1 建立模型为便于研究结合部本身,突出结合部对于整体动力学特性的影响,本文选择一种典型的双螺栓结合部试件作为研究对象,结构如图 1所示。试件通过2个M12螺栓连接。
为了快速对试件的动力学特性有一个预测,这里通过ANSYS Workbench进行简单有限元分析,将结合面处设置为bonded,得到试件的前6阶模态如表 1所示。
为了保证测量精度并且便于建模,这里主要以第一阶模态作为研究对象,其对应的振型为扭转,主要表征螺栓结合部的切向动力学特性。通过模态分析可以看到,试件一阶自由模态主要是扭转,若在一阶模态频率附近进行正弦激振,则试件振动主要由扭转振动构成,根据集中参数法,此时系统可以简化为单自由度振动系统,如图 2所示。
由于结合面的刚度和阻尼为主要影响成分,并且改变预紧力不改变系统原有的质量分布、材料特性和振型,只影响结合面的接触特性和等效动力学参数,因此可以认为单自由度系统振动方程中质量为不变量。并且由振型可知,其扭转振动的模态质量为试件对Z轴的旋转惯量、阻尼和刚度为待辨识量。自由条件下,对模型进行受力分析,则有:
$ {I_{z1}}{{\ddot \theta }_1} = {C_\theta }\left( {{{\dot \theta }_2} - {{\dot \theta }_1}} \right) + {K_\theta }\left( {{\theta _2} - {\theta _1}} \right), $ | (1) |
$ {I_{z2}}{{\ddot \theta }_2} = {C_\theta }\left( {{{\dot \theta }_1} - {{\dot \theta }_2}} \right) + {K_\theta }\left( {{\theta _1} - {\theta _2}} \right) + {M_z}\left( t \right). $ | (2) |
考虑到模型对称性,有
$ {I_{z1}} = {I_{z2}} = {I_z}, $ | (3) |
并令
$ \theta = {\theta _2} - {\theta _1}. $ | (4) |
由式(1)和(2)得到该模型的振动微分方程
$ {I_z}\ddot \theta + 2{C_\theta }\dot \theta + 2{K_\theta }\theta = {M_z}\left( t \right). $ | (5) |
其中:Iz为单个试件对Z轴的旋转惯量,Cθ和Kθ为试件扭转状态下阻尼和刚度系数,Mz(t)为系统的输入力矩。当试件扭转转角比较小时,力矩和角加速度信号可以通过测量力信号和加速度信号与转动半径计算获得:
$ {M_z}\left( t \right) = F\left( t \right) \times {l_F}, $ | (6) |
$ \ddot \theta \left( t \right) = {{\ddot \theta }_2} - {{\ddot \theta }_1} \approx \frac{{{a_2}\left( t \right)}}{{{l_{a2}}}} - \frac{{{a_1}\left( t \right)}}{{{l_{a1}}}}. $ | (7) |
令
$ {K_{\theta eq}} = 2{K_\theta }, $ | (8) |
$ {C_{\theta eq}} = 2{C_\theta }. $ | (9) |
式(5)变为
$ {I_z}\ddot \theta + {C_{\theta eq}}\dot \theta + {K_{\theta eq}}\theta = {M_z}\left( t \right). $ | (10) |
根据线性振动理论,系统的固有频率f和阻尼比ζ满足
$ f = \frac{{{\omega _d}}}{{2{\rm{ \mathsf{ π} }}}} = \frac{{{\omega _n}}}{{2{\rm{ \mathsf{ π} }}}}\sqrt {1 - {\zeta ^2}} . $ | (11) |
其中:
$ {\omega _n} = \sqrt {\frac{{{K_{\theta eq}}}}{{{I_z}}}} , $ | (12) |
$ \zeta = \frac{{{C_{\theta eq}}}}{{2{\omega _n}{I_z}}}. $ | (13) |
ζ和f通过模态分析得到,由上两式即可计算得到系统一阶模态的等效阻尼和刚度系数Cθeq和Kθeq。
如果考虑非线性,可以引入位移项和速度项的多项式函数f(
$ {I_z}\ddot \theta + {C_{\theta eq}}\dot \theta + {K_{\theta eq}}\theta + f\left( {\dot \theta ,\theta } \right) = {M_z}\left( t \right). $ | (14) |
搭建图 3中的测试系统:通过弹性尼龙绳对试件进行悬吊,模拟自由状态,利用力锤敲击激励并测量试件在自由条件下响应,数据采集设备为LMS SCADAS,通过加速度传感器采集响应信号。试件在测量范围(0~1 kHz)内除刚体模态外主要有两阶模态,其中刚体模态远低于第一阶非刚体模态,因此可以忽略悬挂条件的影响认为试件处于自由状态,一阶模态为扭转,二阶模态为前后弯曲,模态振型如图 4所示,与有限元仿真结果一致,改变预紧力后模态振型无明显变化。
当2个螺栓预紧力相同时,得到一阶模态附近频响函数曲线如图 5所示。可以看到,随着预紧力降低,结构的模态频率和共振峰值也随着降低,并且在预紧力较低时模态频率下降越快。另外从频响函数图中可以看出,共振峰右侧具有明显的频响函数畸变,相关系数也出现了显著降低,因此系统可能存在非线性。
利用力矩扳手对2个螺栓施加不同预紧力,在TEST.LAB模态分析软件中获取不同预紧力工况下的系统固有频率和阻尼比,得到模态参数的变化规律,根据式(11)—(13)可以计算出等效线性刚度和阻尼值,实验结果如表 2所示。
|
||||
(5, 5) | 527.2 | 0.5 | 855 889 | 33.125 |
(5, 10) | 538.2 | 0.44 | 892 006 | 29.759 |
(5, 15) | 543.5 | 0.08 | 909 674 | 5.464 |
(5, 20) | 547.7 | 0.06 | 923 625 | 4.129 |
(10, 5) | 534.2 | 0.1 | 878 669 | 6.713 |
(10, 10) | 544.9 | 0.1 | 914 182 | 6.847 |
(10, 15) | 551.0 | 0.07 | 934 782 | 4.847 |
(10, 20) | 553.4 | 0.05 | 942 991 | 3.477 |
(15, 5) | 534.8 | 0.19 | 880 844 | 12.770 |
(15, 10) | 546.9 | 0.07 | 921 174 | 4.811 |
(15, 15) | 551.1 | 0.08 | 935 088 | 5.540 |
(15, 20) | 556.0 | 0.05 | 951 801 | 3.493 |
(20, 5) | 537.0 | 0.09 | 888 070 | 6.074 |
(20, 10) | 547.8 | 0.09 | 924 027 | 6.195 |
(20, 15) | 552.9 | 0.06 | 941 186 | 4.168 |
(20, 20) | 558.3 | 0.05 | 959 892 | 3.508 |
从实验结果可以看出,系统一阶模态的等效刚度随着螺栓预紧力减小而减小,等效阻尼随预紧力减小而增大,当预紧力低于10 N·m时,系统阻尼比显著增大;当预紧力均大于10 N·m时,系统阻尼比很小(小于0.1%)。这可能是由于预紧力较低时结合部的微观接触面积较低导致接触刚度下降,同时微观接触滑移增加导致阻尼增加。另外从2幅图的等势线可以看出,同一结合面的两个螺栓对等效刚度和阻尼的贡献趋势相同,因此同一模态参数下可能对应多组预紧力值。
2.2 激振实验
为了探究结合部在不同激励情况下的动力学特性,这里搭建了激振实验平台,如图 7所示。将试件进行悬吊,固定激振器于工作台上,通过加速度传感器测量其激振响应。
首先对试件采用随机猝发激励和定力幅步进正弦激励,步距为1 Hz,得到的频响函数在一阶模态频率附近与锤击获得的频响函数进行对比,图 8为预紧力为10 N·m时的对比图。
可以看到,10 N步进扫频激励得到的共振峰相比锤击和随机猝发激励的共振峰位置出现了明显的左偏,并且实验中发现正弦激振力幅越大,共振峰频率越低,这说明系统具有刚度渐软特性[14],即系统固有频率和等效刚度会随振动幅值增加而降低。具体来看,因为随机激励和锤击激励的能量较小,系统振幅较小,不易激起刚度渐软非线性,因此随机激励与锤击激励得到的共振峰位置吻合较好,也说明激振器的安装对试件动力学特性的影响较小;而正弦激励的能量集中,系统振幅较大,容易激发系统的刚度渐软非线性,因此造成系统固有频率降低,共振峰左移。
在一阶模态频率附近对试件进行定频激励,通过控制功率放大器和信号发生器的幅值,使输入的正弦激振力幅逐步增加,记录力幅稳定时传感器测得的加速度幅值并绘成曲线,图 9为20 N·m预紧在542 Hz处激振得到的曲线,可以看见加速度幅值随激振力幅值变化呈现出“S”形,这是由于系统具有刚度渐软非线性,频响函数在共振峰及共振峰左侧存在分叉[14],当激振力较小时,系统的非线性未被激发,加速度幅值随激振力幅近似线性变化(k1);当激振力达到一定程度后,系统非线性刚度部分逐渐增强,于是加速度/激振力比值逐渐增加,最终达到另一个激振力和加速度幅值均稳定增加的状态(k2),表现在频响函数上即从较低的分叉点跃变到较高的分叉点上。通过实验发现,激振频率较低时,Ff值较高,随着激振频率增加,Ff值逐渐减小,激振频率增加到一定程度时,图中的跃变现象消失。另外,Ff值与定力幅步进激励得到的共振峰有关,如图中542 Hz激振,Ff值大约为15 N,则以15 N步进激励得到的共振峰位于542 Hz左右,根据Ff值随激振频率的变化规律,也可以得知激振力越高,则共振峰频率越低,即验证了刚度软化特性的存在。
3 基于力状态映射法的参数辨识
力状态映射法是一种用于辨识非线性结构等效动力学参数的常用方法,最初由MARSI等[15]提出,并经CRAWLEY等[16-17]加以发展。
根据力状态映射法,把力表示为速度、位移、加速度的函数,将式(14)改写为
$ {M_z}\left( t \right) - {I_z}\ddot \theta = {C_{\theta eq}}\dot \theta + {K_{\theta eq}}\theta + f\left( {\dot \theta ,\theta } \right). $ | (15) |
获得加速度、速度和位移信号后,在时域中根据多元线性回归便可求得方程的待辨识参数。由于实际中仅测得了加速度的时域信号,需要通过积分获得速度和位移的信息,考虑到直接时域积分会产生趋势项而导致计算误差,而去趋势项会使得信号发生变形而导致计算误差,这里采用频域积分的方法。实际测量中发现,稳定振动的加速度响应信号主要包括基频分量和高次谐波,因此加速度响应可以展开为三角函数叠加的形式,根据三角函数的积分在频域中的对应关系,计算频域积分结果,再通过Fourier逆变换得到时域信号。
$ \ddot x = \sum\limits_{i = 1}^n {{a_i}\cos \left( {{\omega _i}t + {\varphi _i}} \right)} , $ | (16) |
$ \dot x = \int {\ddot x{\rm{d}}t} = \sum\limits_{i = 1}^n {\frac{{{a_i}}}{{{\omega _i}}}\cos \left( {{\omega _i}t + {\varphi _i} - \frac{{\rm{ \mathsf{ π} }}}{2}} \right)} , $ | (17) |
$ x = \int {\dot x{\rm{d}}t} = \sum\limits_{i = 1}^n {\frac{{{a_i}}}{{\omega _i^2}}\cos \left( {{\omega _i}t + {\varphi _i} - {\rm{ \mathsf{ π} }}} \right)} . $ | (18) |
通过激振器定频正弦激励,可以得到试件在共振峰附近的力和加速度时域信号,然后利用力状态映射法来对系统在不同预紧力下的动力学参数进行辨识。
为了使结构的非线性得到充分激发,这里选用共振峰附近10 N定力幅激振数据进行辨识,考虑到测量设备对于高频信号的采集精度,这里将待辨识的非线性项设置为3阶以内多项式非线性项并引入Coulomb摩擦项,待辨识方程:
$ \begin{array}{*{20}{c}} {f\left( {\dot \theta ,\theta } \right) = {a_0} + {a_{01}}\dot \theta + {a_{10}}\theta + {a_{20}}{\theta ^2} + }\\ {{a_{11}}\theta \dot \theta + {a_{02}}{{\left( {\dot \theta } \right)}^2} + {a_{30}}{\theta ^3} + {a_{21}}{\theta ^2}\dot \theta + }\\ {{a_{12}}\theta {{\left( {\dot \theta } \right)}^2} + {a_{03}}{{\left( {\dot \theta } \right)}^3} + {b_1}{\rm{sign}}\left( {\dot \theta } \right).} \end{array} $ | (19) |
辨识结果如表 3所示,其中贡献因子根据所在项对应的信号分量在总信号中的能量占比计算得到。
参数 | 5 N·m | 贡献因子 | 10 N·m | 贡献因子 | 15 N·m | 贡献因子 | 20 N·m | 贡献因子 |
a0 | 0.012 | 0.070 | 0.023 | 0.095 | 0.037 | 0.112 | 0.059 | 0.160 |
a01 | 9.548 | 3.705 | 5.905 | 2.214 | 5.159 | 1.894 | 4.695 | 1.694 |
a10 | 780 372 | 95.445 | 865 020 | 97.017 | 908 750 | 97.609 | 937 372 | 97.707 |
a02 | -2.655 | 0.086 | -2.798 | 0.116 | -2.515 | 0.137 | -3.326 | 0.196 |
a11 | 637.409 | 0.004 | -408.134 | 0.003 | 1 495.226 | 0.014 | 199.938 | 0.002 |
a03 | 76.706 | 0.216 | 40.275 | 0.195 | -0.156 | 0.001 | 13.640 | 0.138 |
a12 | -926 782 | 0.369 | -388 915 | 0.253 | -114 917 | 0.128 | -34 588 | 0.045 |
b1 | -0.018 | 0.106 | -0.026 | 0.107 | -0.035 | 0.106 | -0.022 | 0.059 |
注:未列出项贡献因子小于0.001 |
为了验证辨识结果的准确性,通过MATLAB Simulink建立仿真模型,将实际输入和仿真参数带入模型,利用Runge-Kutta法(ode4)计算得到仿真输出,并与实际输出进行对比,对比结果如图 10所示,其中MSE值为仿真输出和实测输出的均方误差。可以看到仿真预测结果与实验结果吻合很好,证明了辨识结果的准确性。
从辨识结果中可以看到,系统主要存在的非线性项为a12θ(
等效参数 | 测试条件 | 5 N·m | 10 N·m | 15 N·m | 20 N·m |
刚度/(105 N·rad-1) | 模态测试 | 8.56 | 9.14 | 9.35 | 9.60 |
正弦激励 | 7.80 | 8.65 | 9.09 | 9.37 | |
误差/% | 8.82 | 5.38 | 2.82 | 2.35 | |
阻尼比/% | 模态测试 | 0.50 | 0.10 | 0.08 | 0.05 |
正弦激励 | 0.152 | 0.089 | 0.076 | 0.068 | |
误差/% | 69.7 | 10.9 | 5.10 | 36.1 |
4 结论
本文对一种双螺栓结合部的一阶扭转模态进行了简化动力学建模,并通过模态测试和激振测试,利用锤击模态参数识别方法和力状态映射法得到了螺栓结合部在不同预紧力下的等效动力学特性。可以得到如下结论:
1) 螺栓预紧力对于结合部动态特性有重要影响,锤击模态刚度随预紧力降低而降低,模态阻尼比随预紧力降低而增加,并且预紧力较低时变化更加显著。
2) 对于同一低阶模态,结合面有多个螺栓作用时,每个螺栓预紧力对于模态参数变化的影响趋势相同,因此同一模态频率或阻尼比可能对应多组预紧力值。
3) 螺栓结合部存在一定程度的刚度渐软非线性,结合部的等效动力学参数会随着激励信号变化而变化,并且预紧力越低,结合部非线性程度越强,等效动力学参数受激励信号的影响越大,因此在预紧力较低时对结合部的准确建模不能忽略非线性因素的影响。
[1] |
IBRAHIM R A, PETTIT C L. Uncertainties and dynamic problems of bolted joints and other fasteners[J]. Journal of Sound and Vibration, 2005, 279(3-5): 857-936. DOI:10.1016/j.jsv.2003.11.064 |
[2] |
叶佩青, 王仁彻, 赵彤, 等. 机床整机动态特性研究进展[J]. 清华大学学报(自然科学版), 2012, 52(12): 1758-1763. YE P Q, WANG R C, ZHAO T, et al. Recent research advances of whole machine tool dynamics[J]. Journal of Tsinghua University (Science and Technology), 2012, 52(12): 1758-1763. (in Chinese) |
[3] |
蔡力钢, 王峰, 李玲, 等. 栓接结合部动态特性研究进展[J]. 机械工程学报, 2013, 49(9): 158-168. CAI L G, WANG F, LI L, et al. Review on dynamic properties of bolted joints[J]. Journal of Mechanical Engineering, 2013, 49(9): 158-168. (in Chinese) |
[4] |
ČELIČ D, BOLTEŽAR M. Identification of the dynamic properties of joints using frequency-response functions[J]. Journal of Sound and Vibration, 2008, 317(1-2): 158-174. |
[5] |
董冠华, 殷勤, 殷国富, 等. 机床结合部动力学建模与辨识方法的研究[J]. 机械工程学报, 2012, 52(5): 162-168. DONG G H, YIN Q, YIN G F, et al. Research on dynamics modeling and identification of machine tool joints[J]. Journal of Mechanical Engineering, 2012, 52(5): 162-168. (in Chinese) |
[6] |
张辉, 于长亮, 王仁彻, 等. 机床支撑地脚结合部参数辨识方法[J]. 清华大学学报(自然科学版), 2014, 54(6): 815-821. ZHANG H, YU C L, WANG R C, et al. Parameters identification method for machine tool support joints[J]. Journal of Tsinghua University (Science and Technology), 2014, 54(6): 815-821. (in Chinese) |
[7] |
张学良, 范世荣, 温淑花, 等. 基于等效横观各向同性虚拟材料的固定结合部建模方法[J]. 机械工程学报, 2017, 53(15): 141-147. ZHANG X L, FAN S R, WEN S H, et al. Modeling method of fixed joint interfaces based on equivalent transversely isotropic virtual material[J]. Journal of Mechanical Engineering, 2017, 53(15): 141-147. (in Chinese) |
[8] |
孙清超, 黄庆涛, 孙志勇, 等. 基于梯度虚拟材料的栓接结合部连接参数表征[J]. 机械工程学报, 2018, 54(11): 102-109. SUN Q C, HUANG Q T, SUN Z Y, et al. Interface parameter identification of bolted connections based on gradient virtual material[J]. Journal of Mechanical Engineering, 2018, 54(11): 102-109. (in Chinese) |
[9] |
李玲, 蔡力钢, 蔡安江, 等. 不同预紧力下栓接结合部切向等效特性[J]. 振动、测试与诊断, 2013, 33(1): 82-87. LI L, CAI L G, CAI A J, et al. Tangential equivalent properties of the bolted joints in different preload[J]. Journal of Vibration, Measurement & Diagnosis, 2013, 33(1): 82-87. DOI:10.3969/j.issn.1004-6801.2013.01.016 (in Chinese) |
[10] |
WANG M, WANG D, ZHENG G T. Joint dynamic properties identification with partially measured frequency response function[J]. Mechanical Systems and Signal Processing, 2012, 27: 499-512. DOI:10.1016/j.ymssp.2011.09.024 |
[11] |
陈学前, 杜强, 冯加权. 螺栓连接非线性振动特性研究[J]. 振动与冲击, 2009, 28(7): 196-198. CHEN X Q, DU Q, FENG J Q. Nonlinear vibrational characteristic of bolt-joints[J]. Journal of Vibration and Shock, 2009, 28(7): 196-198. DOI:10.3969/j.issn.1000-3835.2009.07.043 (in Chinese) |
[12] |
蔡力钢, 李玲, 郭铁能, 等. 基于力状态映射法辨识非线性结合部动态参数[J]. 机械工程学报, 2011, 47(7): 65-72. CAI L G, LI L, GUO T N, et al. Identification of nonlinear joint parameters with force-state mapping method[J]. Journal of Mechanical Engineering, 2011, 47(7): 65-72. (in Chinese) |
[13] |
孙伟, 李星占, 韩清凯. 螺栓联接梁结构结合部非线性特性参数辨识[J]. 振动工程学报, 2013, 26(2): 185-191. SUN W, LI X Z, HAN Q K. Nonlinear joint parameter identification for bolted beam structure[J]. Journal of Vibration Engineering, 2013, 26(2): 185-191. DOI:10.3969/j.issn.1004-4523.2013.02.005 (in Chinese) |
[14] |
WORDEN K, TOMLINSON G R. Nonlinearity in structural dynamics:Detection, identification, and modeling[M]. Bristol: Institute of Physics Publishing, 2001.
|
[15] |
MASRI S F, CAUGHEY T K. A nonparametric identification technique for nonlinear dynamic problems[J]. Journal of Applied Mechanics, 1979, 46(2): 433-447. DOI:10.1115/1.3424568 |
[16] |
CRAWLEY E F, O'DONNELL K J. Force-state mapping identification of nonlinear joints[J]. AIAA Journal, 1987, 25(7): 1003-1010. DOI:10.2514/3.9733 |
[17] |
CRAWLEY E F, AUBERT A C. Identification of nonlinear structural elements by force-state mapping[J]. AIAA Journal, 1986, 24(1): 155-162. DOI:10.2514/3.9236 |