身管固有频率高效全局灵敏度分析

陈光宋，钱林方，吉 磊

（南京理工大学机械工程学院，南京 210094）

身管固有频率高效全局灵敏度分析

陈光宋，钱林方，吉 磊

（南京理工大学机械工程学院，南京 210094）

为获得影响身管固有频率的关键参数，将身管固有频率写成其参数的混沌多项式形式，提出基于灵敏度的混沌多项式自适应展开，采用LASSO（Least Absolute Shrinkage and Selection Operator）自动选择重要项及其展开系数，根据混沌多项式的正交性直接由混沌多项式展开系数解析获得身管参数的全局灵敏度因子。算例1表明基于灵敏度的多项式自适应展开能够根据变量的重要性选择多项式的展开阶数，算例2通过一个标准模型验证提出方法的有效性和高效性。此外，身管算例表明对于低阶固有频率摇架前支点位置及对应的身管圆柱段直径和身管炮口处直径为核心关键参数。

身管；灵敏度分析；混沌多项式；LASSO

身管为火炮的重要部件之一，其尺寸参数不仅要根据膛压曲线变化规律由强度计算确定［1］，同时还要考虑其他因素，其中身管固有频率即为一个重要的因素，获取影响身管固有频率的关键参数对身管的设计有重要意义。

灵敏度分析被应用于各个领域，如机械工程［2－3］，土木工程［4］，武器系统［5］，损伤诊断［6］等，其功能是量化输入参数的不确定性对输出的影响。通常，灵敏度分析可分为局部灵敏度分析和全局灵敏度分析，其中全局灵敏度分析研究的是参数在完整域内变化对系统输出的影响［11］，大致可归为三类：非参数技术［7－8］，方差分析法［9－10］和矩独立重要性测度［11－12］。针对工程中的不同问题，国内外学者进行了研究和改进，例如李璐祎等［13］为分析参数区域的重要性，基于方差分析法，提出了基本变量区域重要性测度，并给出了相应的算法；张磊刚等［14］分析了高阶核函数的性质，并将其应用到灵敏度分析中，能够提高灵敏度分析的效率和精度；Ge等［15］针对大量输入参数问题，将基于准优化轨迹的基本因素法和Sobol'灵敏度分析方法相结合，有效地解决多输入参数情况下的全局灵敏度分析；Blatman等［16］采用混沌多项式展开，可直接解析地计算出Sobol'灵敏度因子等。

本文采用混沌多项式展开构建身管固有频率与其参数的近似模型，根据多项式展开的稀疏性原理，提出基于灵敏度的混沌多项式自适应展开，并采用LASSO方法自动获取重要项及其展开系数；基于Sobol'灵敏度分析方法［10］，根据混沌多项式的正交特性计算身管参数的全局灵敏度因子，获得影响身管固有频率的关键参数。

1 身管动力学方程

根据身管的结构及受力情况，将身管简化为变截面Timoshenko梁。任意时刻t，梁上距局部坐标原点x处的横向振动方程为［17］：

式中：L为梁长度，y和φ为梁的横向位移和不考虑截面剪切位移的截面转角，ρ为单位体积质量，A（x）为横截面积，I（x）为横截面惯量，kg为与截面形状有关的剪切系数，G为剪切模量，E为弹性模量为横向外载荷。

本文采用时间域上谱单元法进行解算，其具有谱方法的高精度、指数收敛和有限单元法的灵活性。与有限元法相似，谱单元法单元形函数可通过拉格朗日插值得到，即：

式中：插值点ζj不再沿插值方向均匀分布，而是以下等式的根，也称为GLL（Gauss-Lobatto-Legendre）点分布：

式中L′n－1（ζ）为n－1阶Legendre正交多项式的一阶导数。

谱单元法数值积分采用GLL积分方案，即积分点与GLL插值方案的插值点一致：

式中f（ζ）为被积函数，ζi为积分点，wi为权系数。

将插值函数式（3）代入到式（1）和式（2）可得运动方程的弱形式：

式中Me、Ke、δe、Pe分别称为单元质量矩阵、刚度矩阵、节点位移列阵和载荷列阵。

获得变截面Timoshenko梁单元矩阵后可使用与有限元相似的方法组装成系统矩阵，假定身管被分割成N个单元，则系统刚度矩阵、质量矩阵、载荷列阵和位移列阵分别为：

式中为广义特征值问题，解上式即可得到固有频率ω和振型qΔ。

2 高效全局灵敏度分析

2.1 身管固有频率自适应混沌多项式展开

混沌多项式展开近年来在不确定问题分析中逐渐得到应用。其基本思想是将系统的响应展开成有限的p阶混沌多项式和的形式［18］。如图1所示，考虑火炮身管、炮口制退器、炮尾的综合模型，其中炮口制退器和炮尾简化成集中质量。

图1 简化身管模型Fig.1 Model of barrel

令x＝｛L1，L2，L3，L4，L5，φ1，φ2，φ3，φ4，Mw，Mk｝T为相互独立的随机变量，为方便表述，下文中用xj代表x中的第j个变量，f为身管的固有频率，则f展开为x的混沌多项式形式为：

式中βi为多项式展开系数，s为展开项数，φui（x）为多元多项式，n为变量个数，ψuij（xj）是以xj的概率密度函数为权函数的混沌多项式，uij为多项式ψuij（xj）的阶数，ui为各多项式阶数的集合，满足

式（11）中φui（x）包含了所有阶数小于等于p的项，随着输入变量n和展开阶数p的增加，s将呈指数增长。为实现模型的缩减，根据参数的重要性、主项和交叉项重要性的不同，假定STj为变量xj的全局灵敏度，定义：

记tmin为t中的最小值，定义下式：

则满足式（16）的混沌多项式展开集为：

由式（14）可知，若系统对某个参数xj的灵敏度STj较大，则tj较小，进而由式（16）可得(uij)tj较小，因此含有xj的项就越容易被保留下来。

2.2 关键项的选择及展开系数的计算

上节中通过基于灵敏度的混沌多项式自适应展开获得了满足稀疏性原理的混沌多项式展开集Bn，pt，由此可得：

为进一步降低试验设计的规模，采用LASSO［19］方法同时识别重要的项并计算相应系数，定义为：

式中：F为输出样本，argmin（·）为对（·）求最小值，λ为调整参数。

采用最小角回归（LAR）算法［20］解算式（19），并用相对留一误差检验模型的精度：

2.3 基于混沌多项式展开的全局灵敏度

混沌多项式具有正交特性，即

式中δij为克罗尼克符号。

由式（22）可很容易解析地求得f的期望μf和方差为：

式中：E［·］表示求［·］均值，D［·］表示求［·］方差。

为计算参数的全局灵敏度因子，基于Sobol'灵敏度分析方法，可得混沌多项式展开中各项的灵敏度因子为：

则考虑变量xk单独作用的灵敏度Sk和总灵敏度STk分别为：

可见，获得混沌多项式系数后，可直接解析地计算出全局灵敏度因子，与Monte Carlo相比可显著提高计算效率，当获得灵敏度因子后可用于更新自适应展开参数。计算流程如图2所示。

图2 计算流程图Fig.1 Flowchart of computation process

3 算 例

3.1 不同混沌多项式展开比较

取两变量(x1，x2)模型，n＝2，p＝4，比较混沌多项式完备展开，双曲展开［21］和自适应展开。其中双曲展开方案中t＝0.9，自适应展开方案中ST1＝0.72，ST2＝0.3。不同方案的结果如图3、图4和图5所示，图中坐标轴表示变量的阶数。由图中可见完备多项式展开保留了所有小于p的项，双曲多项式展开剔除了对称高阶交叉项，自适应多项式展开根据灵敏度分析结果剔除了不重要变量的高阶项。

图3 完备多项式展开Fig.3 Truncated by complete scheme

图4 双曲多项式展开Fig.4 Truncated by hyperbolic scheme

3.2 解析测试函数

Ishigami函数是灵敏度分析的基准函数［22］：

式中x1、x2、x3为相互独立的随机变量，在[－π，π]上满足均匀分布。计算结果如表1所示，各参数的灵敏度因子结果如表2所示。结果显示采用本文的方法可显著减少试验设计的次数，提高计算效率。灵敏度因子的相对误差远小于Crude MC方法，且不会出现灵敏度因子为负的现象。

表1 Ishigam i函数的前2阶统计矩Tab.1 Ishigam i function and the first two moments

表2 参数全局灵敏度分析Tab.2 Global sensitivity indices of parameters

图5 自适应多项式展开Fig.5 Truncated by adaptive scheme

3.3 身管结构参数灵敏度分析

身管模型的基本参数为：E＝2.1×1011N／m2，G＝1／2.6E，kg＝9／10，ρ＝7 850 kg／m3。L1、L2、L3、L4、L5、φ1、φ2、φ3、φ4、Mw、Mp均服从高斯分布，其均值为现有身管的名义值，不确定参数如表3所示。混沌多项式的最大阶数pmax＝7。前5阶固有频率如表4所示，图6～图10为输入参数对各阶固有频率的灵敏度因子。

身管计算表明对于前五阶固有频率摇架前支点位置L3都是重要关键参数，除四阶频率外φ2对各阶频率都有重要影响，身管炮口处直径φ4对一阶频率和三阶频率有重要影响，φ3对二阶频率和三阶频率较大影响，L1仅对四阶频率有较大影响，φ1对四阶频率有较大影响，对其余各阶频率影响较弱，L5对三阶频率有影响，L2、L4、Mw和Mk对前五阶固有频率均无影响。

表3 输入随机参数Tab.3 Parameters of barrel

图6 一阶固有频率各参数灵敏度Fig.6 Sensitivity indices of parameters for first frequency

图7 二阶固有频率各参数灵敏度Fig.7 Sensitivity indices of parameters for second frequency

图8 三阶固有频率各参数灵敏度Fig.8 Sensitivity indices of parameters for third frequency

表4 前5阶固有频率Tab.4 First five frequencies

图9 四阶固有频率各参数灵敏度Fig.9 Sensitivity indices of parameters for forth frequency

图10 五阶固有频率各参数灵敏度Fig.10 Sensitivity indices of parameters for fifth frequency

4 结 论

本文将火炮身管固有频率参数写成其结构参数的混沌多项式形式，提出基于灵敏度的混沌多项自适应展开，并采用LASSO萃取出关键项保留在最终的模型中，通过基于混沌多项式的全局灵敏度分析获得身管结构的关键参数，算例表明身管不同阶固有频率对应的关键参数不同，对于低阶固有频率摇架前支点位置L3及对应的身管圆柱段直径φ2和身管炮口处直径φ4为核心关键参数，算例结果可为火炮身管的优化设计提供依据。

［1］谈乐斌，张相炎，管红根，等.火炮概论［M］.北京：北京理工大学出版社，2005.

［2］黄毅，刘辉，项昌乐，等.车辆传动系统非线性平移扭转耦合振动响应灵敏度研究［J］.振动与冲击，2014，33（23）：92－99.

HUANG Yi，LIU Hui，XIANG Chang-le，et al.Response sensitivity of nonlinear translation-torsional vibration coupled model of a vehicle transmission system［J］.Journal of Vibration and Shock，2014，33（23）：92－99.

［3］余衍然，李成，姚林泉，等.基于傅里叶幅值检验扩展法的轨道车辆垂向模型全局灵敏度分析［J］.振动与冲击，2014，33（6）：77－81.

YU Yan-ran，LI Cheng，YAO Lin-quan，et al.Global sensitivity analysis on verticalmodel of railway vehicle based on extended Fourier amplitude sensitivity test［J］.Journal of Vibration and Shock，2014，33（6）：77－81.

［4］Hasofer A M.Modern sensitivity analysis of the CESARERisk computer fire model［J］.Fire Safety Journal，2009，44（3）：330－338.

［5］李皓川，孙志礼，王海.基于降维可视化技术的结构可靠性灵敏度分析［J］.兵工学报，2014，35（11）：1876－1882.

LI Hao-chuan，SUN Zhi-li，WANG Hai.Analysis of structural reliability sensitivity based on dimensionality reduction and visualization technique［J］.Acta Armamentar，2014，35（11）：1876－1882.

［6］傅奕臻，魏子天，吕中荣，等.基于时域响应灵敏度分析的板结构损伤识别［J］.振动与冲击，2015，34（4）：117－120.

FU Yi-zhen，WEI Zi-tian，LÜZhong-rong，et al.Damage identication of a plate based on response sensitivity analysis in time domain［J］.Journal of Vibration and Shock，2015，34（3）：117－120.

［7］Helton J C，Johson JD，Sallaberry C J，et al.Survey of sampling-based methods for uncertainty and sensitivity analysis［J］.Reliability Engineering＆System Safety，2006，91（10／11）：1175－1209.

［8］Saltelli A，Marivoet J.Non-parametric statistics in sensitivity analysis formodeloutput：a comparison of selected techniques ［J］.Reliability Engineering＆System Safety，1990，28 （2）：229－253.

［9］Iman R L，Hora S C.A robust measure of uncertainty importance for use in fault tree system analysis［J］.Risk Analysis，1990，10（3）：401－406.

［10］Sobol'I M.Global sensitivity indices for nonlinear mathematicalmodels and their Monte Carlo estimates［J］.Mathematics and Computers in Simulation，2001，55（1／2／3）：271－280.

［11］Borgonovo E.A new uncertainty importance measure［J］.Reliability Engineering＆System Safety，2007，92（6）：771－784.

［12］Liu H B，Chen W，Sudjianto A.Relative entropy based method for probabilistic sensitivity analysis in engineering design［J］.Journal of Mechanical Design，2006，128（2）：326－336.

［13］李璐祎，吕震宙.基本变量区域重要性测度及其稀疏网格解［J］.力学学报，2013，45（4）：569－579.

LILu-yi，LÜZhen-zhou.Regional importancemeasure of the basic variable and its sparse grid solution［J］.Chinese Journal of Theoretical and Applied Mechanics，2013，45（4）：569－579.

［14］张磊刚，吕震宙，吕召燕，等.高阶核函数的性质及其在灵敏度分析中的应用［J］.机械工程学报，2014，50（6）：27－32.

ZHANG Lei-gang，LÜZhen-zhou，LÜZhao-yan，etal.High order properties of kernel functions and their application in sensitivity analysis［J］.Journal of Mechanical Eengineering，2014，50（6）：27－32.

［15］Ge Q，Ciuffo B，Menendez M.Combining screening and metamodel-based methods：an efficient sequential approach for the sensitivity analysis of model outputs［J］.Reliability Engineering＆System Safety，2015，134：334－344.

［16］Blatman G，Sudret B.Efficient computation of global sensitivity indices using sparse polynomial chaos expansions ［J］.Reliability Engineering＆System Safety，2010，95（11）：1216－1229.

［17］Timoshenko S.Vibration Problem in Engineering［M］.New York：D.Van Nostrang Company，INC.，1937.

［18］Wiener N.The homogeneous chaos［J］.American Journal of Mathematics，1938，60：897－936.

［19］Tibshirani R.Regression shrinkage and selection via the Lasso［J］.Journal of the Royal Statistical Society：Series B （Methodological），1996，58（1）：267－288.

［20］Efron B，Hastie T，Johnstone I，et al.Least angle regression ［J］.The Annals of Statistics，2004，32（2）：407－499.

［21］Blatman G，Sudret B.Adaptive sparse polynomial chaos expansion based on least angle regression［J］.Journal of Computational Physics，2011，230（6）：2345－2367.

［22］Saltelli A，Chan K，Scott E M.Sensitivity Analysis［M］.New Jersey：Jhon Wildy and Sons，2000.

An effective global sensitivity analysismethod for natural frequencies of a barrel

CHEN Guang-song，QIAN Lin-fang，JILei

（School of Mechanical Engineering，Nanjing University of Science ang Technology，Nanjing 210094，China）

In order to obtain a barrel's key parameters affecting its natural frequencies，the polynomial chaos expansion was used to describe the relation of its natural frequencies and its structural parameters.The adaptive expansion strategy for the polynomial chaos based on the sensitivity was proposed to obtain polynomial chaos expansionswith different orders.The important terms and corresponding expansion coefficients were obtained automatically by using the least absolute shrinkage and selection operator（LASSO）.Based on the orthogonality of the polynomial chaos，the global sensitivity factors for structural parameters of the barrel were obtained directly with coefficients of the polynomial chaos expansion.Example 1 indicated that the proposed method can select expansion orders according to the importance of variables.A benchmark example was presented to demonstrate the effectiveness and higher efficiency of the proposed method.In addition，the numerical examples of barrels indicated that the position of the front fulcrum of the cradle，the corresponding cylinder diameters of barrels and the diameters at barrelsmuzzles are the key parameters affecting lower order natural frequencies of barrels.

barrel；sensitivity analysis；polynomial chaos；LASSO

TJ301

A

10.13465／j.cnki.jvs.2015.21.006

国家自然科学基金（11472137，51205207）

2015－02－02 修改稿收到日期：2015－04－14

陈光宋男，博士生，1987年生

钱林方男，教授，1961年生