The stochastic dynamics of nanoscale mechanical oscillators immersed in a viscous fluid

Ornstein-Uhlenbeck噪声驱动的过阻尼广义Langevin方程的共振钟苏川;彭皓;张路【摘要】本文研究了受外部周期信号激励的线性过阻尼广义Langevin方程的随机共振现象.本文将系统内噪声建模为指数型关联Ornstein-Uhlenbeck噪声,系统外噪声建模为双态噪声,并利用随机平均法和积分变换算法推导出系统响应的一阶稳态矩和稳态响应振幅的解析表达式.对解析结果的分析表明,该线性过阻尼广义Langevin方程具有丰富的共振行为,即系统的稳态响应振幅随噪声的特征参数、周期激励信号的频率及部分系统参数的变化而出现广义随机共振.【期刊名称】《四川大学学报（自然科学版）》【年(卷),期】2018(055)004【总页数】7页(P698-704)【关键词】广义Langevin方程;线性过阻尼振子;Ornstein-Uhelnbeck噪声;随机共振【作者】钟苏川;彭皓;张路【作者单位】四川大学空天科学与工程学院空天信息处理与应用中心,成都610064;西南交通大学数学学院,成都611756;四川大学数学学院,成都610064【正文语种】中文【中图分类】O291 引言上世纪80年代，Benzi等[1-3]在研究第四纪冰川问题时提出了随机共振的概念.他们用双稳态系统对天气系统进行建模，将短期气候的随机波动用高斯白噪声来描述.在一定强度的噪声调节下，气候的周期响应与外加微弱周期驱动会实现同步,这种现象被称为随机共振(Stochastic Resonance，SR).此后,关于随机共振的理论和实验研究引起了人们极大的关注.较早的研究普遍认为，系统要产生随机共振必须满足非线性、周期激励和噪声这三个条件.随着研究的深入，即使上述三个条件不完全满足的系统也被观测到随机共振现象.例如，随机共振可以出现在由乘性色噪声驱动的线性系统中[4]，甚至是白噪声激励的具有随机阻尼项的欠阻尼线性系统中[5，6].传统的随机共振是指系统响应的信噪比随噪声强度的非单调变化现象.广义随机共振(Generalized Stochastic Resonance，GSR)的概念则是由Gitterman等提出的，系指系统响应的某些函数(如矩、自相关函数、功率谱或信噪比等)随系统的某些特征参数(如激励振幅、频率或噪声的强度、相关率等)非单调变化的现象[7，8]. 关于随机共振的研究大多关注外噪声情形.此时，噪声来源于外部环境，与系统自身的动力学特性无关.然而，对于很多小尺度系统(如活体细胞中的生化反应，纳米尺度晶体表面上的催化反应等)而言，系统自身的涨落起着关键作用，此时内噪声对系统动力学行为的影响必须加以考虑[9，10].特别地，内噪声如何影响系统的动力学行为，系统能否在内噪声的作用下出现随机共振或广义随机共振现象，及共振现象产生的充分必要条件等问题都有待解决.对于这些问题的研究也将为后续将随机共振理论用于参数估计和弱信号检测等提供基础.在以往的研究中,系统内噪声通常被建模为高斯白噪声,此时粒子的均方位移为:〈x2(t)〉=2Dt，即正常扩散.然而,当粒子在稠密的黏性流体、具有内部自由度的流体以及湍流中运动时,粒子的均方位移是时间t的幂律非线性函数,即〈x2(t)〉=2Dδtδ,δ≠1，即反常扩散[11，12].反常扩散现象可以通过广义Langevin方程刻画.实际上，认为噪声在不同时刻不相关只是一种近似，实际上真正的白噪声是不存在的.具有非零相关时间的噪声就是色噪声.一种最简单的色噪声是奥恩斯坦-乌伦贝克(Ornstein-Uhlenbeck，OU)噪声，它具有指数型关联函数[9，13，14].目前，尚没有文献研究指数关联噪声驱动的广义Langevin方程的共振行为.通过将外部环境扰动建模为能够解析求解的对称双态噪声，计算稳态响应振幅的解析表达式，我们讨论了OU噪声驱动的线性过阻尼广义Langevin方程的共振行为.进一步，根据获得的稳态响应振幅的解析表达式，我们观察到系统的内、外噪声都会诱导出广义随机共振现象.本文还给出了系统的稳态响应振幅关于外噪声的噪声强度能够诱导产生广义随机共振现象的充分必要条件.由于文中推导得到的系统响应的解析表达式是未作任何近似的精确表达式,因而结果不必限制在弱噪声和弱激励信号范围内,适用于任意的噪声强度和周期信号的振幅、频率,能解释更为广泛的物理现象.2 模型一个处于强噪声环境的线性过阻尼广义Langevin方程为[15,16]:A0sin(Ωt)+F(t)(1)其中,x(t)为单位质量粒子的位移,A0、Ω为周期激励信号的振幅和频率,ω2为势阱的固有频率,γ>0为阻尼系数.随机力Z(t)为固有频率ω2受到的随机扰动,本文将其视为外噪声,建模为对称双态噪声[17-19]:Z(t)能取两个离散值{>-a,a}，a>0，取值的概率满足Ps(-a)=Ps(a)=1/2，Z(t)的均值和方差为:〈Z(t)〉=0,〈Z(t)Z(s)〉=a2exp(-v|t-s|)(2)其中,a2为噪声强度,v为噪声相关率.双态噪声是一种简单而常见的色噪声,具有相对简单的性质,并且在极限条件(噪声相关率v→+∞)时,双态噪声将退化为高斯白噪声.此外,将双态噪声以乘性方式引入线性系统可得到一个精确可解的色噪声模型，故在双态噪声作用下,系统的动力学性质的研究具有普适的理论和实际意义[18，20，21 ].(1) 式中的K(t)为系统阻尼核函数,F(t)为粒子受到的热噪声,满足涨落耗散定理[9，22]:〈F(t)F(s)〉=kBT·K(t-s)(3)其中kB为波尔兹曼常数,T是介质温度.在均匀介质中,粒子的运动轨迹是一个马尔科夫过程,对应的随机力F(t)应建模为高斯白噪声:〈F(t)〉=0,〈F(t)F(s)〉=kBT·δ(t-s).此时,K(t)=δ(t)，从而广义Langevin方程退化为经典Langevin方程;而在非均匀介质中,扩散过程具有记忆性,粒子的运动轨迹是一个非马尔科夫过程,此时随机力F(t)应建模为色噪声.本文选取的内噪声F(t)为零均值、指数型高斯色噪声,即Orstein-Uhlenbeck(OU)噪声[9，13，14]:〈F(t)〉=0,(4)其中D为噪声强度,τ为噪声关联时间.需要说明的是，本文将双态噪声Z(t)视为系统外噪声,OU噪声F(t)视为系统内噪声,假设Z(t)和F(t)不相关,即〈F(t)Z(s)〉=0(5)3 系统响应的一阶稳态矩及稳态响应振幅将(1)式两边取平均,得到一阶矩〈x(t)〉所满足的方程:〈Z(t)x(t)〉=A0sin(Ωt)(6)其中，〈x(t)〉和满足(7)由于(6)式中出现了新的耦合项〈Z(t)x(t)〉，需要对该项解耦.利用Shapiro-Loginov公式[23]可得〈Z(t)x(t)〉所满足的微分方程为(8)此时(8)式中又出现了新的耦合项需要继续寻找其所满足的方程.将(1)式两边同时乘以Z(t)再取平均，获得最后一个封闭方程：(9)综上，(6)～(9) 式构成关于的方程组.由Laplace变换有(10)其中,x1(0)，x3(0)为初值条件，Xi(s)=xi(t)e-stdt,i=1,2,3,4.本文选取内噪声F(t)为OU噪声,利用涨落耗散定理(3)式,不难由OU噪声的自相关函数((4)式)推导出此时阻尼核函数K(t)的Laplace域表达式为:(11)其中称为广义噪声强度.由(10)式看出,系统响应的一阶矩所满足的方程在s域为一个四元非齐次线性方程组,解方程组可得Xi(s),i=1,2,3,4的解析表达式:(12)将(12)式作逆Laplace变换,得系统响应一阶矩解析表达式为:xi(t)=A0Hi1(t-u)sin(Ωu)du+[Hi2(t)x1(0)+Hi3(t)x3(0)](13)其中，为Hik(t)的Laplace变换,而可以通过(10)式确定.特别地，(14)为保证(13)式所给出的系统响应一阶矩的稳定性(意味着Hik(s)=0,k=1,2,3不能有正实根),系统参数需满足如下不等式[8]:(15)本文的后续都是基于(15)式成立的情况进行讨论的.当t→∞时，初始条件对系统响应的影响将逐渐消失,则系统响应的一阶稳态矩为:(x(t))as =(x(t))|t→∞ =A1 sin(Ωt+ϖ1 )(16)其中A1和φ1为系统响应一阶稳态矩〈x(t)〉as的振幅和相移,A1又称为系统稳态响应振幅,它们可根据的表达式获得：(17)其中4 讨论4.1 系统稳态响应振幅A1关于OU噪声(内噪声)参数的共振行为根据推导获得的稳态响应振幅A1((17)式)，我们在图1(a)和(b)中分别给出了A1作为OU噪声F(t)(内部噪声)的噪声强度D和关联时间τ的函数所对应的3D曲面图和对应到不同噪声强度(D=1,D=2,D=3)下的投影曲线图A1(τ).从图1可以看出，对于任意给定的内噪声强度D，随着内噪声关联时间τ的变化.稳态响应振幅的曲线表现出了典型的非单调变化，这意味着广义随机共振现象产生了.具体地，随着OU噪声相关时间τ的增加，稳态响应振幅A1(τ)增大到极大值(波峰相)随后又达到极小值(波谷相)，这种现象被称为波峰-波谷共振行为.进一步，随着OU噪声强度D的增大，这种波峰-波谷共振行为表现的更加明显.同时，随着OU噪声强度D的减小(如D=1)，波峰-波谷共振行为逐渐退化为传统的单峰共振行为.图1 (a) 稳态响应振幅A1作为OU噪声F(t)的噪声强度D和关联时间τ的函数所对应的3D曲面图，参数为A0=1,Ω=0.5,ω=1,kBT=1,γ=0.05,a2=0.2,v=3；(b) 给定OU噪声强度D取值下稳态响应振幅A1随OU噪声相关时间τ的变化曲线，其他参数同图1(a).Fig.1(a) The output amplitude A1 as the function of D and τ of the internal OU noise F(t)， the other parameters are A0=1,Ω=0.5,ω=1,kBT=1,γ=0.05,a2=0 .2,v=3; (b) The output amplitude A1 as the function of τ with different valu es of noise intensity D.4.2 系统稳态响应振幅A1关于双态噪声(外噪声)参数的共振行为在分析了系统稳态响应振幅A1关于系统内噪声F(t)参数D,τ的共振行为后，我们继续研究A1关于系统外噪声Z(t)(双态噪声)强度a2和相关率v的共振行为.图2 (a)～(c)分别给出了对应的3D曲面图和对应于噪声强度a2和相关率v的投影曲线图.从图2(a)可以看出，稳态响应振幅A1关于外噪声强度a2和相关率v都能表现出明显的非单调变化曲线，即系统能够出现广义随机共振现象.进一步地，根据图2 (b)和(c)可以看出，随着外噪声相关率v的增加(v=0.05,v=0.1,v=0.2)，A1(a2)的共振峰逐渐平缓，直至消失.这说明此时只有较小的噪声相关率能够诱导显著的单峰共振现象.同时，随着外噪声强度a2的增加(a2=0.1,a2=0.4,a2=0.8)，A1(v)的共振行为更加明显，共振峰的峰值A1,max(v)增加，对应峰值点位置右移. 图2 (a) 稳态响应振幅A1作为双态噪声Z(t)的噪声强度a2和相关率v的函数所对应的3D曲面图，参数为A0=1,Ω=0.2,ω=1,kBT=1,γ=2.5,D=1.5,τ=0.6；(b) 给定Z(t)噪声相关率v取值下稳态响应振幅A1(a2)变化曲线，其他参数同图2(a)；(c) 给定Z(t)噪声强度a2取值下稳态响应振幅A1(v)变化曲线，其他参数同图2(a).Fig.2(a) The output amplitude A1 as the function of a2 and v of the external dic hotomous noise Z(t)，the other parameters are: A0=1,Ω=0.2,ω=1,kBT=1,γ=2.5,D=1.5,τ=0.6; (b) T he output amplitude A1 as the function of a2 with different values of noise correlation rate v; (c) The output amplitude A1 as the function of v with dif ferent values of noise intensity a2.4.3 系统稳态响应振幅A1关于双态噪声(外噪声)噪声强度a2共振行为的进一步分析根据上面的讨论，系统稳态响应振幅A1关于外噪声强度a2的共振行为依赖于系统参数的选取.同时,考虑到本文所讨论的模型共振行为能够被控制的主要因素是调整噪声强度a2，因而在本节我们将根据(17)式所给出的系统稳态响应振幅A1的表达式分析系统稳态响应振幅A1关于外噪声强度a2能够诱导产生共振行为的充分必要条件.根据(17)式，系统稳态响应振幅A1(a2)在处达到极大值，(18)其中根据共振的稳定性条件((15)式)可推导得此时稳态响应振幅A1关于外噪声强度a2能够诱导产生广义随机共振的充分必要条件为(19)图3给出的是根据共振的充分必要条件(19)式在参数空间τ-γ上做的相图.图中灰色区域对应的参数τ和γ取值能够诱导系统的稳态响应振幅A1产生共振行为,白色区域对应的参数τ和γ取值组合无法诱导A1产生随机共振行为.从图3可以看出，稳态响应振幅A1关于外噪声强度a2能够产生随机共振现象非常依赖系统参数取值的选取.图3 系统稳态响应振幅A1关于外噪声强度a2能够诱导产生广义随机共振现象对应于τ-γ平面的相图，参数ω2=1,kBT=1，Ω=2,D=1,v=0.1.Fig.3The phase diagrams for SR phenomenon of A1(a2) in the τ-γ plane with the other parameters ω2=1,kBT=1,Ω=2,D=1,v=0.1.图4给出的是其他参数同图3所取时稳态响应振幅A1(a2)针对不同阻尼系数γ取值下的曲线图.此时，根据图3给出的稳态响应振幅A1(a2)能够诱导共振行为的阻尼系数γ的取值范围为[0.26,0.67].对应图4的模拟结果，当阻尼系数γ=0.3和γ=0.4(γ的取值落入能够诱导共振行为的参数范围)时，对应的稳态相应振幅曲线A1(a2)表现出先增大后减小的非单调变换过程，也即典型的随机共振行为.此外，当γ=0.1和γ=0.8(γ的取值未落入能够诱导共振行为的参数范围)时，从模拟结果看，对应的稳态相应振幅曲线A1(a2)呈现出严格单调的变化趋势，即此时系统不能诱导随机共振行为产生.从模拟结果看来，我们给出的稳态相应振幅曲线A1(a2)能够有诱导产生随机共振行为的相图是可靠的，这也为后续产生和利用随机共振现象提供了坚实的理论基础.图4 不同的阻尼系数γ取值下，系统稳态响应振幅A1关于外噪声强度a2的变化曲线图，对应参数为：ω2=1,kBT=1,Ω=2,D=1,v=0.1,τ=0.1， (a) γ=0.1; (b) γ=0.3; (c) γ=0.4;(d)γ=0.8.Fig.4The output amplitude A1 as the function of noise intensity a2 with differen t values of damping coefficient γ， the other parameters are: ω2=1,kBT=1,Ω=2,D=1,v=0.1,τ=0.1， (a) γ=0.1; (b) γ=0.3; (c) γ=0.4;(d)γ=0.8.5 结论本文研究了周期激励信号作用下的线性过阻尼广义Langevin方程受Ornstein-Uhlenbeck噪声激励下的共振行为.通过随机平均法和积分变换算法本文推导了系统响应一阶稳态矩和稳态响应振幅的精确表达式(17) 式，讨论了稳态响应振幅关于噪声强度能够诱导共振行为的充分必要条件(19)式，并验证了其有效性，从而为后续将该随机共振现象用于工程实际应用提供基础.参考文献:【相关文献】[1] Benzi R，Sutera A，Vulpiani A.The mechanism of stochastic resonance [J].J Phys A： Math Theor， 1981，14: 453.[2] Benzi R，Parisi G，Sutera A，et al.Stochastic resonance in climatic change [J].Tellus， 1982， 34: 10.[3] Nicolis C.Stochastic aspects of climatic transitions response to a periodic forcing [J].T ellus， 1982， 34: 1.[4] 王传毅，任芮彬，邓科，周期调制噪声驱动的具质量涨落的欠阻尼谐振子的随机共振[J].四川大学学报:自然科学版， 2016， 53: 1183.[5] Berdichevsky V，Gitterman M.Stochastic resonance in linear systems subject to multiplicative and additive noise [J].Phys Rev E， 1999， 60: 1494.[6] Li J H，Han Y X.Phenomenon of stochastic resonance caused by multiplicative asymmetric dichot omous noise [J].Phys Rev E， 2006， 74: 051115.[7] Gitterman M.Harmonic oscillator with fluctuating damping parameter [J].Phys Rev E ， 2004， 69: 041101.[8] Gitterman M.Classical harmonic oscillator with multiplicative noise [J].Physica A， 20 05， 352: 309.[9] 包景东.经典和量子耗散系统的随机模拟方法 [M].北京: 科学出版社， 2009.[10] 龚玉兵，侯中怀，辛厚文.色噪声性质对一氧化氮还原反应体系随机共振的影响 [J].高等学校化学学报，2005， 26: 2331.[11]Klafter J， Sokolov I M.Anomalous diffusion spreads its wings [J].Phys World， 2005， 18: 29.[12] Metzler R，Klafter J.The random walk’s guide to anomalous diffusion: a fractional dynamics approach [J].Phys Rep， 2000， 339: 1.[13]Peter R.Thermally activated escape with potential fluctuations driven by an Ornstein-Uhlenbeck process [J].Phys Rev E， 1995， 1579: 1.[14] Cabrera J L，Gorronogoitia J，de la Rubia F J.Coherence enhancement in nonlinear systems subject to multiplicative Orn stein-Uhlenbeck noise [J].Phys Rev E， 2002， 66， 022101.[15]Kou S C.Stochastic modeling in nanoscale biophysics: subdiffusion within proteins [J]. Ann Appl Stat， 2008， 2: 501.[16]Kou S C， Xie X S.Generalized Langevin equation with fractional gaussian noise: subdiffusi on within a single protein molecule [J].Phys Rev Lett， 2004， 93: 180603.[17]Gitterman M.Harmonic oscillator with multiplicative noise: nonmonotonic dependence on the strength and the rate of dichotomous noise [J].Phys Rev E， 2003， 67: 057103. [18] Laio F，Ridolfi L，D’Odorico P.Noise-induced transitions in state-dependent dichotomous processes [J].Phys Rev E， 2008， 78: 031137.[19] Xu Y，Wu J，Zhang H Q，et al.Stochastic resonance phenomenon in an underdamped bistable system driven by we ak asymmetric dichotomous noise [J].Nonlinear Dynam， 2012， 70: 531.[20] Bena I，Van den Broeck C，Kawai R，et al.Nonlinear response with dichotomous noise [J].Phys Rev E， 2002， 66: 045603. [21] Bena I.Dichotomous markov noise: exact results for out-of-equilibrium systems [J].Int J Mod Phys B， 2006， 20: 2825.[22] 张季谦，辛厚文.噪声在非线性化学体系中作用的研究现状[J].化学进展， 2001, 13: 241.[23] Shapiro V E，Loginov V M.Formulae of differentiation and their use for solving stochastic equations [J]. Physica A， 1978， 91: 563.。

㊀第４９卷第４期煤炭科学技术Ｖｏｌ ４９㊀Ｎｏ ４㊀㊀２０２１年４月ＣｏａｌＳｃｉｅｎｃｅａｎｄＴｅｃｈｎｏｌｏｇｙ㊀Ａｐｒ．２０２１㊀移动扫码阅读邓志刚．动静载作用下煤岩多场耦合冲击危险性动态评价技术［Ｊ］．煤炭科学技术，２０２１，４９（４）：１２１－１３２ ｄｏｉ：１０ １３１９９／ｊ ｃｎｋｉ ｃｓｔ ２０２１ ０４ ０１５ＤＥＮＧＺｈｉｇａｎｇ．Ｍｕｌｔｉ－ｆｉｅｌｄｃｏｕｐｌｉｎｇｄｙｎａｍｉｃｅｖａｌｕａｔｉｏｎｍｅｔｈｏｄｏｆｒｏｃｋｂｕｒｓｔｈａｚａｒｄｃｏｎｓｉｄｅｒｉｎｇｄｙｎａｍｉｃａｎｄｓｔａｔｉｃｌｏａｄ［Ｊ］．ＣｏａｌＳｃｉｅｎｃｅａｎｄＴｅｃｈｎｏｌｏｇｙ，２０２１，４９（４）：１２１－１３２ ｄｏｉ：１０ １３１９９／ｊ ｃｎｋｉ ｃｓｔ ２０２１ ０４ ０１５动静载作用下煤岩多场耦合冲击危险性动态评价技术邓㊀志㊀刚１，２（１．煤炭科学技术研究院有限公司安全分院，北京㊀１０００１３；２．煤炭资源高效开采与洁净利用国家重点实验室，北京㊀１０００１３）摘㊀要：深部开采冲击地压灾害孕灾过程中既有静态基础量又有动态变化量，剧增的原岩应力与覆岩断裂㊁井下爆破等引起的动载扰动是诱发冲击地压灾害的源头，因此实现冲击危险性快速㊁高精度评价必须综合考虑动静载作用㊂笔者开展了典型煤岩霍普金森压杆试验及数值模拟，分析了动载对煤岩体破坏作用以及对应力场的影响，针对应力变化可以直接引起介质中震动波波速变化，且波速变化前的幅值与变化幅度均受应力场影响这一特性，掌握了震动场与应力场的耦合关系，建立了多场耦合冲击危险性动态评价技术：以原岩应力场表示煤岩孕灾过程的静态基础量，以采动应力场和震动场表示煤岩孕灾过程的动态变化量，以波速异常指数㊁波速梯度指数㊁应力异常指数㊁应力梯度指数为评价指标可实现煤岩冲击危险性动态评价㊂研究结果表明：动载作用下能量以震动波形式传递，造成应力场的重新分布，应力呈现分区传递特点，并且在能量达到某一阈值后引起煤岩损伤破坏，但无论动载直接作用在岩石上还是煤体上，岩石是能量传递路径，煤层是能量耗散㊁释放主体，破坏主要发生在煤体中㊂多场耦合冲击危险性评价技术在某工作面经现场应用，在工作面逐渐揭露断层过程中冲击危险性由强冲击危险性降低到中等冲击危险性，现场监测数据表明评价结果与现场实际相符㊂关键词：动静载荷；冲击危险性；震动场；多场耦合；动态评价中图分类号：ＴＤ３２４㊀㊀㊀文献标志码：Ａ㊀㊀㊀文章编号：０２５３－２３３６（２０２１）０４－０１２１－１２Ｍｕｌｔｉ－ｆｉｅｌｄｃｏｕｐｌｉｎｇｄｙｎａｍｉｃｅｖａｌｕａｔｉｏｎｍｅｔｈｏｄｏｆｒｏｃｋｂｕｒｓｔｈａｚａｒｄｃｏｎｓｉｄｅｒｉｎｇｄｙｎａｍｉｃａｎｄｓｔａｔｉｃｌｏａｄＤＥＮＧＺｈｉｇａｎｇ１，２（１．ＭｉｎｅＳａｆｅｔｙＴｅｃｈｎｏｌｏｇｙＢｒａｎｃｈ，ＣｈｉｎａＣｏａｌＲｅｓｅａｒｃｈＩｎｓｔｉｔｕｔｅ，Ｂｅｉｊｉｎｇ㊀１０００１３，Ｃｈｉｎａ；２．ＳｔａｔｅＫｅｙＬａｂｏｒａｔｏｒｙｏｆＣｏａｌＭｉｎｉｎｇａｎｄＣｌｅａｎＵｔｉｌｉｚａｔｉｏｎ，Ｂｅｉｊｉｎｇ㊀１０００１３，Ｃｈｉｎａ）收稿日期：２０２０－１２－０２；责任编辑：朱恩光基金项目：国家科技重大专项资助项目（２０１６ＺＸ０５０４５００３－００６－００２）；国家自然科学基金面上资助项目（５１６７４１４３）作者简介：邓志刚（１９８１ ），男，吉林长春人，研究员，博士，中国煤炭科工集团三级首席科学家㊂Ｔｅｌ：０１０－８４２６１５８１，Ｅ－ｍａｉｌ：ｄｅｎｇｚｈｉｇａｎｇ２００４＠１６３．ｃｏｍＡｂｓｔｒａｃｔ：Ｓｔａｔｉｃｂａｓｉｃｑｕａｎｔｉｔｙａｎｄｄｙｎａｍｉｃｖａｒｉａｔｉｏｎｑｕａｎｔｉｔｙｅｘｉｓｔｉｎｔｈｅｐｒｏｃｅｓｓｏｆｒｏｃｋｂｕｒｓｔｉｎｄｅｅｐｍｉｎｉｎｇ．Ｄｙｎａｍｉｃｌｏａｄｄｉｓｔｕｒｂａｎｃｅａｎｄｔｈｅｉｎｃｒｅａｓｉｎｇｏｆｉｎ－ｓｉｔｕｓｔｒｅｓｓｆｉｅｌｄａｒｅｔｈｅｓｏｕｒｃｅｏｆｒｏｃｋｂｕｒｓｔ．Ｔｈｅｒｅｆｏｒｅ，ｔｈｅｄｙｎａｍｉｃａｎｄｓｔａｔｉｃｌｏａｄｍｕｓｔｂｅｃｏｎｓｉｄｅｒｅｄｃｏｍｐｒｅｈｅｎｓｉｖｅｌｙｉｎｔｈｅｆａｓｔａｎｄｈｉｇｈ－ｐｒｅｃｉｓｉｏｎｅｖａｌｕａｔｉｏｎｏｆｒｏｃｋｂｕｒｓｔｈａｚａｒｄ．Ｈｏｐｋｉｎｓｏｎｐｒｅｓｓｕｒｅｂａｒｅｘｐｅｒｉｍｅｎｔｓａｎｄｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎｓｗｅｒｅｃａｒｒｉｅｄｏｕｔｔｏａｎａｌｙｚｅｔｈｅｉｎｆｌｕｅｎｃｅｓｏｆｄｙｎａｍｉｃｌｏａｄｏｎｔｈｅｄａｍａｇｅａｎｄｓｔｒｅｓｓｆｉｅｌｄｏｆｔｈｅｃｏａｌｒｏｃｋ．Ｉｎｖｉｅｗｏｆｔｈｅｆａｃｔｔｈａｔｔｈｅｃｈａｎｇｅｏｆｓｔｒｅｓｓｃｏｕｌｄｄｉ⁃ｒｅｃｔｌｙｃａｕｓｅｔｈｅｃｈａｎｇｅｏｆｖｉｂｒａｔｉｏｎｗａｖｅｖｅｌｏｃｉｔｙａｎｄｔｈｅａｍｐｌｉｔｕｄｅｂｅｆｏｒｅａｎｄａｆｔｅｒｔｈｅｃｈａｎｇｅｏｆｗａｖｅｖｅｌｏｃｉｔｙｗｅｒｅａｆｆｅｃｔｅｄｂｙｔｈｅｓｔｒｅｓｓｆｉｅｌｄ，ｔｈｅｃｏｕｐｌｉｎｇｒｅｌａｔｉｏｎｓｈｉｐｂｅｔｗｅｅｎｖｉｂｒａｔｉｏｎｆｉｅｌｄａｎｄｓｔｒｅｓｓｆｉｅｌｄｗａｓｍａｓｔｅｒｅｄａｎｄｔｈｅｍｕｌｔｉ－ｆｉｅｌｄｃｏｕｐｌｉｎｇｄｙｎａｍｉｃｅｖａｌｕａｔｉｏｎｍｅｔｈｏｄｏｆｒｏｃｋｂｕｒｓｔｈａｚａｒｄｗａｓｅｓｔａｂｌｉｓｈｅｄ．Ｉｎｔｈｅｐｒｏｃｅｓｓｏｆｃａｔａｓｔｒｏｐｈｅ，ｔｈｅｉｎ－ｓｉｔｕｓｔｒｅｓｓｆｉｅｌｄｒｅｐｒｅｓｅｎｔｓｔｈｅｓｔａｔｉｃｆｏｕｎｄａｔｉｏｎｑｕａｎｔｉｔｙ，ａｎｄｔｈｅｍｉｎｉｎｇｓｔｒｅｓｓｆｉｅｌｄａｎｄｔｈｅｖｉｂｒａｔｉｏｎｆｉｅｌｄｒｅｐｒｅｓｅｎｔｔｈｅｄｙｎａｍｉｃｃｈａｎｇｅｑｕａｎｔｉｔｙ．Ｔｈｅｗａｖｅｖｅｌｏｃｉｔｙａｎｏｍａｌｙｉｎｄｅｘ，ｗａｖｅｖｅｌｏｃｉｔｙｇｒａｄｉｅｎｔｉｎ⁃ｄｅｘ，ｓｔｒｅｓｓａｎｏｍａｌｙｉｎｄｅｘａｎｄｓｔｒｅｓｓｇｒａｄｉｅｎｔｉｎｄｅｘａｒｅｕｓｅｄａｓｅｖａｌｕａｔｉｏｎｉｎｄｅｘｅｓｔｏｒｅａｌｉｚｅｄｙｎａｍｉｃｅｖａｌｕａｔｉｏｎｏｆｒｏｃｋｂｕｒｓｔｈａｚａｒｄ．Ｔｈｅｒｅ⁃ｓｕｌｔｓｓｈｏｗｔｈａｔｔｈｅｅｎｅｒｇｙｉｓｔｒａｎｓｍｉｔｔｅｄｉｎｔｈｅｆｏｒｍｏｆｖｉｂｒａｔｉｏｎｗａｖｅｕｎｄｅｒｄｙｎａｍｉｃｌｏａｄ，ｒｅｓｕｌｔｉｎｇｉｎｔｈｅｒｅｄｉｓｔｒｉｂｕｔｉｏｎｏｆｓｔｒｅｓｓｆｉｅｌｄ．Ｔｈｅｓｔｒｅｓｓｐｒｅｓｅｎｔｓｔｈｅｃｈａｒａｃｔｅｒｉｓｔｉｃｓｏｆｚｏｎａｌｔｒａｎｓｍｉｓｓｉｏｎ，ａｎｄｃａｕｓｅｓｔｈｅｄａｍａｇｅｏｆｃｏａｌａｎｄｒｏｃｋｗｈｅｎｔｈｅｅｎｅｒｇｙｒｅａｃｈｅｓａｃｅｒｔａｉｎｔｈｒｅｓｈｏｌｄ．Ｈｏｗｅｖｅｒ，ｎｏｍａｔｔｅｒｔｈｅｄｙｎａｍｉｃｌｏａｄｄｉｒｅｃｔｌｙａｃｔｓｏｎｔｈｅｒｏｃｋｏｒｔｈｅｃｏａｌ，ｔｈｅｒｏｃｋｉｓｔｈｅｅｎｅｒｇｙｔｒａｎｓｆｅｒｐａｔｈ，ｔｈｅｃｏａｌｓｅａｍｉｓｔｈｅｍａｉｎｂｏｄｙｏｆｅｎｅｒｇｙｄｉｓｓｉｐａｔｉｏｎａｎｄｒｅｌｅａｓｅ，ａｎｄｔｈｅｆａｉｌｕｒｅｍａｉｎｌｙｏｃｃｕｒｓｉｎｔｈｅｃｏａｌ．Ｔｈｅｍｕｌｔｉ－ｆｉｅｌｄｃｏｕｐｌｉｎｇｄｙｎａｍｉｃｅｖａｌｕａｔｉｏｎｍｅｔｈｏｄｏｆｒｏｃｋｂｕｒｓｔ１２１２０２１年第４期煤炭科学技术第４９卷ｈａｚａｒｄｗａｓａｐｐｌｉｅｄｏｎａｃｅｒｔａｉｎｗｏｒｋｉｎｇｆａｃｅ．Ｔｈｅｒｏｃｋｂｕｒｓｔｈａｚａｒｄｗａｓｒｅｄｕｃｅｄｆｒｏｍｓｔｒｏｎｇｔｏｍｅｄｉｕｍｉｎｔｈｅｐｒｏｃｅｓｓｏｆｇｒａｄｕａｌｌｙｅｘｐｏｓｉｎｇｆａｕｌｔｓ．Ｔｈｅｆｉｅｌｄｍｏｎｉｔｏｒｉｎｇｄａｔａｓｈｏｗｅｄｔｈａｔｔｈｅｅｖａｌｕａｔｉｏｎｒｅｓｕｌｔｓｗｅｒｅｃｏｎｓｉｓｔｅｎｔｗｉｔｈｔｈｅａｃｔｕａｌｓｉｔｕａｔｉｏｎ．Ｋｅｙｗｏｒｄｓ：ｄｙｎａｍｉｃａｎｄｓｔａｔｉｃｌｏａｄｓ；ｒｏｃｋｂｕｒｓｔｈａｚａｒｄ；ｖｉｂｒａｔｉｏｎｆｉｅｌｄ；ｍｕｌｔｉ－ｆｉｅｌｄｃｏｕｐｌｉｎｇ；ｄｙｎａｍｉｃｅｖａｌｕａｔｉｏｎ０㊀引㊀㊀言我国多数矿井进入深部开采阶段，冲击地压灾害频度㊁强度显著增加［１］，冲击地压防治工作任重道远㊂２０１８年８月１日，国家煤矿安全监察局印发的‘防治煤矿冲击地压细则“开始实施，规定： 开采具有冲击倾向性的煤层必须进行冲击危险性评价 ， 开采冲击地压煤层必须进行采区㊁采掘工作面冲击危险性评价 ， 当评估煤层有冲击倾向性时，应当进行冲击危险性评价 ，并且以冲击危险性评价结果作为冲击地压监测㊁卸压等工作开展的依据㊂目前冲击危险性评价方法较多㊂一类是以冲击地压主要诱因为切入点的冲击危险性静态评价技术，如窦林名等［２］提出的综合指数法，综合考虑了岩体结构㊁力学特性㊁地质因素等条件㊂姜福兴等［３］采用模糊数学的方法，用垂直应力与煤体单轴抗压强度的比值㊁弹性能量指数２个指标评价煤体的冲击危险性，且根据应力叠加原理建立了冲击危险性评价模型，后又在此基础上提出了冲击地压分类评价技术手段㊂张科学等［６］综合考虑开采深度㊁冲击倾向性㊁煤层顶底板性质㊁地质构造㊁开采技术提出了基于层次分析法的煤层冲击危险性模糊综合评价模型㊂张宏伟等［７］应用地质动力区划方法对煤矿冲击危险进行评价㊂邓志刚［１０］基于三维地应力场反演技术开展了相关研究，综合考虑构造应力㊁采动影响等因素，实现了对采区宏观区域的冲击危险评价㊂欧阳振华等［１１］考虑瓦斯作用，将煤层气属性㊁抽采效果分析作为一类地质因素㊁开采技术条件，提出一种含瓦斯煤冲击危险性改进型综合指数评价方法㊂但是由于冲击地压致灾机理不清，灾害孕育㊁发展㊁发生的过程中影响因素繁杂，以及复杂多变的采掘及地质条件，致使静态评价方法主要是宏观上为煤层开采前的防冲工作提供一定参考，缺少对于采掘过程中因局部区域地质及开采条件变化㊁卸压措施等因素引起的冲击危险性动态变化的量化能力，因此，另一类基于现场监测数据的冲击危险性动态评价技术是当前研究工作的重点，如刘少虹等［１２］基于地音与电磁波ＣＴ探测数据提出的冲击危险性层次化评价方法；李宏艳等［１４］基于微震监测数据建立的考虑响应能量和无响应时间的冲击危险性动态评价技术㊂姜福兴等［１５］应用矿压观测法观测冲击地压工作面支架压力㊁立柱压缩量，判断工作面顶板来压规律，结合巷道的变形及其围岩应力分布进行观测，评价及预测冲击危险性㊂何学秋等［１７］采用电磁辐射法评价冲击危险性，主要参数为电磁辐射强度和脉冲数㊂曹民远等［１９］采用数值模拟和理论计算的方法分析了采掘工作面应力扰动叠加的影响，提出了近直立煤层动态权重评价法的计算体系㊂但是冲击地压的孕灾过程中既有静态基础量，又有动态变化量，因此目前仅依靠单一理论或方法快速㊁高精度的进行冲击危险性评价难度较大㊂我国煤矿进入深部开采后，剧增的原岩应力场成为冲击地压灾害发生的必要条件㊂覆岩断裂㊁井下爆破等带来的强动载扰动易成为诱发冲击灾变的充分条件，但目前冲击危险性评价的研究工作中少有兼顾动静载综合作用的理论或方法㊂为此，笔者以震动场㊁采动应力场表示孕灾过程中动态变化量，以原岩应力场表示孕灾过程中静态基础量㊂提出了波速异常指数㊁波速梯度指数㊁应力异常指数㊁应力梯度指数４个冲击危险性评价指标，并在此基础上建立了多场耦合冲击危险性动态评价技术以实现井下高精度冲击危险性动态评价㊂１㊀煤岩动载破坏试验分析１．１㊀典型煤岩动载破坏霍普金森压杆试验分离式霍普金森压杆（ＳＨＰＢ）试验系统（图１）由压杆系统㊁测量系统和数据采集处理系统３个部分组成㊂图１㊀ＳＨＰＢ试验装置Ｆｉｇ．１㊀ＳＨＰＢｅｘｐｅｒｉｍｅｎｔａｌｄｅｖｉｃｅ当动载试块受到不同气压后获得不同初速度撞击入射杆，在杆内产生入射脉冲εｉ，试件在该应力作用下产生高速变形，同时产生反射脉冲εｒ和透射脉冲εｔ㊂如图２所示㊂选取强冲击倾向性煤样试件４２２１邓志刚等：动静载作用下煤岩多场耦合冲击危险性动态评价技术２０２１年第４期个，中砂岩试件４个，尺寸均为ø５０ｍｍˑ１００ｍｍ㊂本次试验煤岩样取样点分别为某典型冲击地压矿井３－１煤回风大巷ＨＦ６导点处顶板和３１１１０２工作面煤层㊂煤岩物理力学参数见表１㊂分别采用气压０．２㊁０．４㊁０．６㊁０．８ＭＰａ发射子弹，撞击入射杆，记录其入射㊁反射和透射波曲线㊂图２㊀ＳＨＰＢ试验原理Ｆｉｇ．２㊀ＰｒｉｎｃｉｐｌｅｏｆＳＨＰＢｅｘｐｅｒｉｍｅｎｔａｌ㊀㊀煤样㊁岩样入射波㊁反射波和透射波曲线如图３㊁图４所示，仅出示驱动应力为０．２㊁０．４㊁０．８ＭＰａ时的结果㊂对比分析可知，随着撞击杆驱动应力增加，入射波波速幅值㊁入射波波速变化率均有所增加，反射波和透射波波峰和波谷增高，透射波持续时间缩短，这也和冲击地压发生的突然㊁猛烈性质一致㊂１．２㊀典型煤岩动载破坏数值模拟采取有限元方法对煤岩霍普金森压杆试验进行模拟，进一步分析动载作用下煤岩体损伤破坏机理㊂数值模型如图５所示㊂模拟试件分为煤样㊁岩样㊁煤－岩组合样，岩－煤组合样，其中煤－岩组合样是指震动波入射端在煤上，岩－煤组合样是指震动波入射端在岩石上㊂煤样㊁岩样尺寸为ø５０ｍｍˑ１００ｍｍ，煤岩组合样中煤㊁岩样尺寸均为ø５０ｍｍˑ５０ｍｍ㊂入射杆㊁透射杆材料参数按钢材设定［２０］，密度为７７９４ｋｇ／ｍ３，弹性模量为２１１ＧＰａ，泊松比为０．２８５㊂表１㊀煤岩物理力学参数Ｔａｂｌｅ１㊀Ｐｈｙｓｉｃａｌａｎｄｍｅｃｈａｎｉｃａｌｐａｒａｍｅｔｅｒｓｏｆｃｏａｌａｎｄｒｏｃｋ试样密度／（ｋｇ㊃ｍ－３）单轴抗压强度／ＭＰａ弹性模量／ＧＰａ泊松比抗拉强度／ＭＰａ内摩擦角／（ʎ）黏聚力／ＭＰａ煤样１３２５．４０３８．７６２３．４７４０．２８２２．４９３１８．５２１３．８９４岩样２１１１．９８４０．４３４７．３９５０．２２２２．８３９３５．６０１５．５２５图３㊀煤样不同气压下的波形Ｆｉｇ．３㊀Ｗａｖｅｆｏｒｍｓｏｆｃｏａｌｕｎｄｅｒｄｉｆｆｅｒｅｎｔａｉｒｐｒｅｓｓｕｒｅ图４㊀岩样不同气压下的波形Ｆｉｇ．４㊀Ｗａｖｅｆｏｒｍｓｏｆｒｏｃｋｕｎｄｅｒｄｉｆｆｅｒｅｎｔａｉｒｐｒｅｓｓｕｒｅ３２１２０２１年第４期煤炭科学技术第４９卷图５㊀霍普金森试验数值模型Ｆｉｇ．５㊀ＳＨＰＢｅｘｐｅｒｉｍｅｎｔｎｕｍｅｒｉｃａｌｍｏｄｅｌ煤岩物理力学参数见表２㊂加载在入射杆端部的震动波信号为ＳＨＰＢ试验中不同气压驱动子弹记录的入射杆应变波信号㊂不同震动波作用下煤岩体应力㊁损伤分布如图６ 图９所示，限于篇幅煤样㊁岩样仅出示驱动应力为０．２㊁０．４㊁０．８ＭＰａ时的结果，煤岩组合样仅出示驱动应力为０．２ＭＰａ和０．８ＭＰａ时的结果㊂分析可知，震动波作用引起煤岩应力重新分布，应力传递呈现分区传递特点，即存在应力传递优势面㊂在震动波波速峰值㊁波速变化率较低时，震动波对煤岩介质表２㊀数值模拟参数Ｔａｂｌｅ２㊀Ｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎｐａｒａｍｅｔｅｒｓ试样弹性模量／ＧＰａ泊松比密度／（ｋｇ㊃ｍ－３）屈服强度／ＭＰａ单轴抗压强度／ＭＰａ内摩擦角／（ʎ）黏聚力／ＭＰａ煤样３．４７４０．３２１３２０１７．２５２４．６０１８．５２１３．８９０岩样７．６８３０．２３２５１９４４．９７５０．２７４３．１０１０．６５６图６㊀煤样应力与损伤分布情况Ｆｉｇ．６㊀Ｓｔｒｅｓｓａｎｄｄａｍａｇｅｄｉｓｔｒｉｂｕｔｉｏｎｏｆｃｏａｌｓｐｅｃｉｍｅｎ没有破坏作用，即震动波对煤岩介质的破坏与损伤存在阈值㊂煤岩体发生破坏的位置同时是单元受拉损伤㊁受压损伤极值位置，因此震动波作用下煤岩体破坏模式为拉压复合破坏㊂无论震动波直接作用在岩石上还是煤上，煤岩组合试件的破坏主要发生在煤体上，说明岩石是能量传播的路径，煤体是能量耗４２１邓志刚等：动静载作用下煤岩多场耦合冲击危险性动态评价技术２０２１年第４期图７㊀岩样应力与损伤分布情况Ｆｉｇ．７㊀Ｓｔｒｅｓｓａｎｄｄａｍａｇｅｄｉｓｔｒｉｂｕｔｉｏｎｏｆｒｏｃｋｓｐｅｃｉｍｅｎ图８㊀煤－岩样应力与损伤分布情况Ｆｉｇ．８㊀Ｓｔｒｅｓｓａｎｄｄａｍａｇｅｄｉｓｔｒｉｂｕｔｉｏｎｏｆｃｏａｌ－ｒｏｃｋｓｐｅｃｉｍｅｎ５２１２０２１年第４期煤炭科学技术第４９卷图９㊀岩－煤样应力与损伤分布情况Ｆｉｇ．９㊀Ｓｔｒｅｓｓａｎｄｄａｍａｇｅｄｉｓｔｒｉｂｕｔｉｏｎｏｆｒｏｃｋ－ｃｏａｌｓｐｅｃｉｍｅｎ散㊁释放的主体，这也符合冲击地压主要发生在煤层中的事实㊂１．３㊀震动场与煤岩冲击危险性的关联依据采煤工作面和掘进工作面煤岩体破坏失稳主要形式，结合ＳＨＰＢ试验和数值模拟研究结果，煤岩体震动场与冲击危险性的关系总结如下：①震动波是能量传递的载体，震动波所具有的能量超过一定阈值时可引起煤岩破坏，易诱发冲击地压灾害㊂②震动波传递引起应力分布变化，应力传递沿优势面进行㊂随着震动波能量增加，优势面周围易出现煤岩损伤破坏，引起煤岩冲击灾变㊂③当震源位于岩层时，能量传递速度较快，在煤岩界面发生衰减，煤体在震动波作用下发生破坏；当震源位于煤层时，煤体对震动波传递速度相对较慢，能量多耗散在煤层中，主要诱发煤体破坏，对岩层造成的破坏较小㊂２㊀煤岩动㊁静载冲击危险性评价指标考虑动静载作用煤岩冲击危险性评价指标包括应力场相关指标和震动场相关指标，其中静载作用主要表现为应力场的变化，动载作用主要引起震动场的变化㊂２．１㊀应力场冲击危险性评价指标基于煤矿冲击地压应力控制理论［２１］，煤岩体冲击破坏是应力作用的结果，一是取决于应力绝对值大小，二是应力梯度变化㊂因此，建立应力异常指数和应力梯度指数㊂应力异常指数表征一定区域内不同位置应力差异的指标，计算公式为γσ＝σｒ－σｍｉｎσｍａｘ－σｍｉｎˑ１０（１）式中：γσ为应力异常指数；σｒ为监测区域某点应力，ＭＰａ；σｍａｘ㊁σｍｉｎ分别为监测区域内实时应力最大值和最小值，ＭＰａ㊂应力梯度指数是表征一定区域内不同位置应力变化速度差异的指标，计算公式为ｇσ＝ｇσｒ－ｇσｍｉｎｇσｍａｘ－ｇσｍｉｎˑ１０（２）式中：ｇσ为应力梯度异常指数；ｇσｒ为监测区域内某一点的应力场梯度；ｇσｍａｘ㊁ｇσｍｉｎ分别为监测区域内应力最大㊁最小梯度㊂２．２㊀震动场冲击危险性评价指标综上，震动场波速绝对值㊁变化速率对煤岩破坏有显著影响㊂因此，提出表征震动波波速的波速异常指数和表征震动波波速变化速率的波速梯度指数，作为２个基于震动场的冲击危险性动态评价指数㊂波速异常指数表征一定区域内不同位置震动波波速的差异，计算公式为γθ＝θｒ－θｍｉｎθｍａｘ－θｍｉｎˑ１０（３）式中：γθ为波速异常指数；θｒ为监测区域某点震动波６２１邓志刚等：动静载作用下煤岩多场耦合冲击危险性动态评价技术２０２１年第４期波速，ｍ／ｓ；θｍａｘ㊁θｍｉｎ分别为监测区域内震动波波速最大值和最小值，ｍ／ｓ㊂波速梯度指数ｇθ是通过震动场波速变化速率表征煤岩体发生冲击地压的危险程度，计算公式为ｇθ＝ｇθｒ－ｇθｍｉｎｇθｍａｘ－ｇθｍｉｎˑ１０（４）式中：ｇθ为波速梯度异常指数；ｇθｒ为监测区域内某一点的震动波波速梯度；ｇθｍａｘ㊁ｇθｍｉｎ为监测区域内震动波波速最大㊁最小梯度㊂３㊀煤岩多场耦合冲击危险性动态评价技术结合笔者以往研究［２２］和上述研究成果可知，一方面煤岩应力场改变可以直接引起介质中震动波波速变化，且波速变化前的幅值与变化幅度均与应力场大小相关；另一方面，震动场传递会造成煤岩应力场的重新分布㊂因此，考虑动㊁静载作用开展煤岩冲击危险性动态评价关键在于分析震动场－应力场的耦合作用㊂煤炭开采之前，煤岩体处于重力和构造应力组成的原岩应力场之中；开采过程中，煤岩体形成采动应力场；原岩应力场和采动应力场相互作用，煤岩体损伤变形，震动产生，以弹性波的形式向外传播形成震动场㊂冲击地压是原岩应力场㊁采动应力场和震动场综合作用的结果，煤岩体中多场耦合关系如图１０所示㊂图１０㊀煤岩体中多场耦合关系Ｆｉｇ．１０㊀Ｆｉｅｌｄｉｎｃｏａｌｒｏｃｋｍａｓｓａｎｄｉｔｓｃｏｕｐｌｉｎｇｒｅｌａｔｉｏｎｓｈｉｐ为了准确描述煤岩体中各种场的关系，从冲击危险性评价角度建立统一数学模型Ｒ（ｔｉ，ｓ；ｍｊ）＝０㊀㊀（ｉ，ｊ＝１，２，３， ）（５）式中：ｔｉ为场的变量，一般情况下有多个，既可以是标量也可以是矢量；ｓ为场的源或者汇，通常只有一个；ｍｊ为煤岩体的物理性质变量，如弹性模量㊁泊松比㊁剪切模量㊁波速等多个变量㊂基于该函数煤岩体中的３种场的冲击危险性评价具体表达式如下：１）原岩应力场为Ｙ（ｈ，ｃ，ｆ；ρ，μ）＝０（６）式中：ｈ为采深；ｃ为地应力；ｆ为体积力；ρ为煤岩体密度；μ为泊松比㊂２）震动场为Ｓ（ｘ，ｙ，ｚ，ｔ，Ｅ，ｆ；ρ，μ）＝０（７）式中：ｘ㊁ｙ㊁ｚ为震源的位置坐标；ｔ为发震时间；Ｅ为震源能量㊂３）采动应力场为Ｆ（ｕ，ｆ；ρ，μ）＝０（８）式中：ｕ为位移㊂３．１㊀原岩应力场与采动应力场（ＲＭ）耦合冲击危险性评价模型㊀㊀原岩应力场冲击危险性评价指标见表３㊂原岩应力场冲击危险性指数定义为Ｒ＝（Ｒ１＋Ｒ２＋Ｒ３＋Ｒ４）／４（９）其中，Ｒ１㊁Ｒ２㊁Ｒ３㊁Ｒ４为不同评价指标得分㊂原岩应力冲击危险性反映煤岩体自身发生冲击地压的固有属性，其数值大小反映了煤岩体采动后，发生自发型冲击地压的可能性和危险性㊂原岩应力场冲击危险性指数取值与冲击危险等级关系见表４㊂表３㊀原岩应力场冲击危险性评价指标Ｔａｂｌｅ３㊀Ｒｏｃｋｂｕｒｓｔｈａｚａｒｄｅｖａｌｕａｔｉｏｎｉｎｄｅｘｓｏｆｉｎ－ｓｉｔｕｓｔｒｅｓｓｆｉｅｌｄ变量影响因素阈值分值Ｒ１开采深度ｈｈɤ４００ｍ１４００ｍ＜ｈɤ６００ｍ２６００ｍ＜ｈɤ８００ｍ３ｈ＞８００ｍ４Ｒ２向落差大于３ｍ的断层推进的工作面或巷道，工作面或掘进工作面至断层的距离ＬｄＬｄȡ１００ｍ１５０ｍɤＬｄ＜１００ｍ２２０ｍɤＬｄ＜５０ｍ３Ｌｄ＜２０ｍ４Ｒ３向背斜或向斜推进的工作面或巷道，工作面或掘进工作面与之距离ＬｚＬｚȡ５０ｍ１２０ｍɤＬｚ＜５０ｍ２１０ｍɤＬｚ＜２０ｍ３Ｌｚ＜１０ｍ４Ｒ４同一水平煤层冲击地压发生次数ｎｎ＝０１ｎ＝１２２ɤｎ＜３３ｎȡ３４㊀㊀采动应力冲击危险指标包括：应力异常指数和应力梯度指数㊂二者取值与冲击危险等级之间的关系见表５㊁表６㊂７２１２０２１年第４期煤炭科学技术第４９卷表４㊀原岩应力场冲击危险性等级划分标准Ｔａｂｌｅ４㊀Ｒｏｃｋｂｕｒｓｔｈａｚａｒｄｃｌａｓｓｉｆｉｃａｔｉｏｎｃｒｉｔｅｒｉａｂａｓｅｄｏｎｉｎ－ｓｉｔｕｓｔｒｅｓｓｆｉｅｌｄ阈值冲击危险性评价指数冲击危险等级Ｒɤ１１无１＜Ｒ＜２２弱２ɤＲ＜３３中等Ｒȡ３４强表５㊀应力异常指数冲击危险性等级划分标准Ｔａｂｌｅ５㊀Ｒｏｃｋｂｕｒｓｔｈａｚａｒｄｃｌａｓｓｉｆｉｃａｔｉｏｎｃｒｉｔｅｒｉａｂａｓｅｄｏｎｓｔｒｅｓｓａｎｏｍａｌｙｉｎｄｅｘ阈值冲击危险性评价指数冲击危险等级γσɤ１１无１＜γσ＜３２弱３ɤγσ＜５３中等γσȡ５４强表６㊀应力梯度指数冲击危险性等级划分标准Ｔａｂｌｅ６㊀Ｒｏｃｋｂｕｒｓｔｈａｚａｒｄｃｌａｓｓｉｆｉｃａｔｉｏｎｃｒｉｔｅｒｉａｂａｓｅｄｏｎｓｔｒｅｓｓｇｒａｄｉｅｎｔｉｎｄｅｘ阈值冲击危险性评价指数冲击危险等级ｇθɤ１１无１＜ｇθ＜３２弱３ɤｇθ＜５３中等ｇθȡ５４强㊀㊀基于原岩应力场与采动应力场耦合的冲击危险性评价模型为ＤＲＭ＝ａ１Ｒ＋ｂ１γσ＋ｃ１ｇσ（１０）㊀㊀其中：ＤＲＭ是原岩应力场与采动应力场耦合的冲击危险性评价指数；ａ１，ｂ１，ｃ１分别为原岩应力场和采动应力场耦合冲击危险性评价权重系数，不同矿井取值不同㊂原岩应力场与采动应力场耦合的冲击危险性指数取值与冲击危险等级之间的关系见表７㊂表７㊀原岩应力场与采动应力场耦合冲击危险性等级划分标准Ｔａｂｌｅ７㊀Ｒｏｃｋｂｕｒｓｔｈａｚａｒｄｃｌａｓｓｉｆｉｃａｔｉｏｎｃｒｉｔｅｒｉａｂａｓｅｄｏｎｃｏｕｐｌｉｎｇｏｆｉｎ－ｓｉｔｕｓｔｒｅｓｓｆｉｅｌｄａｎｄｍｉｎｉｎｇｓｔｒｅｓｓｆｉｅｌｄ阈值冲击危险性评价指数冲击危险等级ＤＲＭɤ１１无１＜ＤＲＭ＜３２弱３ɤＤＲＭ＜５３中等ＤＲＭȡ５４强３．２㊀原岩应力场与震动场（ＲＳ）耦合冲击危险性评价模型㊀㊀震动场冲击危险性指标包括：波速异常指数和波速梯度指数㊂二者取值与冲击危险等级之间的关系见表８㊁表９㊂原岩应力场与震动场耦合的冲击危险性评价模型为ＤＲＳ＝ａ２Ｒ＋ｂ２γθ＋ｃ２ｇθ（１１）㊀㊀其中：ＤＲＳ为原岩应力场和震动场耦合的冲击危险性评价指数；ａ２，ｂ２，ｃ２为原岩应力场和震动场耦合冲击危险性评价权重系数，不同矿井取值不同㊂原岩应力场与震动场耦合的冲击危险性指数取值与冲击危险等级之间的关系见表１０㊂表８㊀波速异常指数冲击危险性等级划分标准Ｔａｂｌｅ８㊀Ｒｏｃｋｂｕｒｓｔｈａｚａｒｄｃｌａｓｓｉｆｉｃａｔｉｏｎｃｒｉｔｅｒｉａｂａｓｅｄｏｎｗａｖｅｖｅｌｏｃｉｔｙａｎｏｍａｌｙｉｎｄｅｘ阈值冲击危险性评价指数冲击危险等级γθɤ１１无１＜γθ＜３２弱３ɤγθ＜５３中等γθȡ５４强表９㊀波速梯度指数冲击危险性等级划分标准Ｔａｂｌｅ９㊀Ｒｏｃｋｂｕｒｓｔｈａｚａｒｄｃｌａｓｓｉｆｉｃａｔｉｏｎｃｒｉｔｅｒｉａｂａｓｅｄｏｎｗａｖｅｖｅｌｏｃｉｔｙｇｒａｄｉｅｎｔｉｎｄｅｘ阈值冲击危险性评价指数冲击危险等级ｇθɤ１１无１＜ｇθ＜３２弱３ɤｇθ＜５３中等ｇθȡ５４强表１０㊀原岩应力场与震动场耦合冲击危险性等级划分标准Ｔａｂｌｅ１０㊀Ｒｏｃｋｂｕｒｓｔｈａｚａｒｄｃｌａｓｓｉｆｉｃａｔｉｏｎｃｒｉｔｅｒｉａｂａｓｅｄｏｎｃｏｕｐｌｉｎｇｏｆｉｎ－ｓｉｔｕｓｔｒｅｓｓｆｉｅｌｄａｎｄｖｉｂｒａｔｉｏｎｆｉｅｌｄ阈值冲击危险性评价指数冲击危险等级ＤＲＳɤ１１无１＜ＤＲＳ＜３２弱３ɤＤＲＳ＜５３中等ＤＲＳȡ５４强３．３㊀采动应力场与震动场（ＭＳ）耦合冲击危险性评价模型㊀㊀采动应力场与震动场耦合冲击危险性评价模型为ＤＭＳ＝ａ３γσ＋ｂ３ｇσ＋ｃ３γθ＋ｄ３ｇθ（１２）㊀㊀其中：ＤＭＳ为采动应力场与震动场耦合冲击危险性评价指数；ａ３，ｂ３，ｃ３，ｄ３分别为应力异常指数，应力梯度指数，波速异常指数，波速梯度指数的权重系数，不同矿井取值不同㊂采动应力场与震动场耦合的冲击危险性指数取值与冲击危险等级之间的关系见表１１㊂３．４㊀多场耦合（ＲＭＳ）冲击危险性动态评价模型冲击地压发生的本质是煤岩体具有的冲击能量８２１邓志刚等：动静载作用下煤岩多场耦合冲击危险性动态评价技术２０２１年第４期超过围岩吸收能量的极限㊂应力场可以表现煤岩体表１１㊀采动应力场与震动场耦合冲击危险性等级划分标准Ｔａｂｌｅ１１㊀Ｒｏｃｋｂｕｒｓｔｈａｚａｒｄｃｌａｓｓｉｆｉｃａｔｉｏｎｃｒｉｔｅｒｉａｂａｓｅｄｏｎｃｏｕｐｌｉｎｇｏｆｍｉｎｉｎｇｓｔｒｅｓｓｆｉｅｌｄａｎｄｖｉｂｒａｔｉｏｎｆｉｅｌｄ阈值冲击危险性评价指数冲击危险等级ＤＭＳɤ１１无１＜ＤＭＳ＜３２弱３ɤＤＭＳ＜５３中等ＤＭＳȡ５４强未受扰动的地应力场和受采动影响而形成的采动应力场，是煤岩体承受应力的状态量㊂震动场主要表现煤岩体无法承受外部高应力差作用发生损伤破坏，在此过程中以震动形式释放出能量的时空域，可以表现煤岩体积蓄能量的过程㊂冲击地压的不仅发生在高应力区，也发生在煤岩体由低应力区向高应力区转化的过程中，采用煤岩体多场耦合的方法可以充分全面评价监测区域的冲击危险性㊂基于上述对ＲＭ耦合㊁ＲＳ耦合和ＭＳ耦合的冲击危险性评价模型，构建煤岩体多场耦合（ＲＭＳ）冲击危险性动态评价模型㊂冲击危险性指数算法如下Ｄ＝ＤＲＭ＋ＤＲＳ＋ＤＭＳ（１３）多场耦合冲击危险性评价指数Ｄ与冲击危险性等级的对应关系见表１２㊂表１２㊀多场耦合（ＲＭＳ）冲击危险性等级划分标准Ｔａｂｌｅ１２㊀Ｒｏｃｋｂｕｒｓｔｈａｚａｒｄｃｌａｓｓｉｆｉｃａｔｉｏｎｃｒｉｔｅｒｉａｂａｓｅｄｏｎｍｕｌｔｉ－ｆｉｅｌｄｃｏｕｐｌｉｎｇ阈值冲击危险性评价指数冲击危险等级Ｄɤ５１无５＜Ｄ＜１０２弱１０ɤＤ＜１５３中等Ｄȡ１５４强４㊀工程应用选取典型冲击地压矿井３１１２０２工作面为现场，开展相关应用㊂４．１㊀工作面概况３１１２０２工作面是该矿井１２盘区第２个回采工作面，是首个沿空回采工作面，位于１２盘区北部，为３１１２０１接续工作面，东部以１２盘区辅运大巷为界，西部至１２盘区西部边界，南部为实体煤，北部为正在回采的３１１２０１工作面，保护煤柱宽度６ｍ㊂该工作面采用走向长壁综合机械化一次采全高采煤法，采高５．２５ｍ，工作面倾斜长度２９９ｍ，走向长度３１４０ｍ，全部垮落法管理顶板，两回采巷道采用液压支架进行超前支护㊂工作面布置如图１１所示㊂图１１㊀３１１２０２工作面布置Ｆｉｇ．１１㊀ＬａｙｏｕｔｏｆＮｏ．３１１２０２ｍｉｎｉｎｇｆａｃｅ经鉴定，３－１煤及其顶底板均具有弱冲击倾向性，３－１煤层冲击危险等级为中等冲击危险㊂３１１２０２工作面所在地层构造形态总体为一向北西倾斜的单斜构造，倾向３００ʎ ３２０ʎ㊁倾角１ʎ ３ʎ，地层产状沿走向及倾向均有一定变化，沿走向发育有宽缓的波状起伏㊂３１１２０２工作面受ＤＦ１９㊁ＤＦ１８㊁Ｆ２２㊁Ｆ２４断层影响较大，其中ＤＦ１９断层影响最为显著，该断层走向长度约１２００ｍ，落差６．５ １０．０ｍ，预计影响３１１２０２工作面走向长度５６０ｍ，对生产过程中的冲击地压灾害影响最大㊂３１１２０２工作面主要断层情况见表１３，３１１２０２工作面煤层顶底板结构特征见表１４㊂表１３㊀３１１２０２工作面断层特征Ｔａｂｌｅ１３㊀ＦａｕｌｔｃｈａｒａｃｔｅｒｉｓｔｉｃｓｏｆＮｏ．３１１２０２ｍｉｎｉｎｇｆａｃｅ断层走向／（ʎ）倾向／（ʎ）倾角／（ʎ）性质落差／ｍＤＦ１８３１２４２７０正断层０ ５．０ＤＦ１９２９６２６４９正断层６．５ １０．０Ｆ２２２８５１５３０正断层１．１Ｆ２４３５７９０４６正断层０．３表１４㊀３１１２０２工作面煤层顶底板结构特征Ｔａｂｌｅ１４㊀ＳｔｒｕｃｔｕｒａｌｃｈａｒａｃｔｅｒｉｓｔｉｃｓｏｆｃｏａｌｓｅａｍｒｏｏｆａｎｄｆｌｏｏｒｉｎＮｏ．３１１２０２ｍｉｎｉｎｇｆａｃｅ顶底板岩性厚度／ｍ平均厚度／ｍ基本顶细粒砂岩９．２５ １９．７０１５．８４直接顶砂质泥岩２．２８ １２．８５８．５０直接底砂质泥岩４．６９ １２．９９７．６８基本底细粒砂岩５．２１ ２１．４５１４．８２４．２㊀多场耦合冲击危险性动态评价原岩应力场包括重力场和构造应力场，通过地应力测试及三维反演可得到㊂采动应力场通过应力在线监测系统监测得到㊂在３１１２０２回风巷生产帮安设应力在线监测系统，距离开切眼６０ｍ生产帮侧９２１２０２１年第４期煤炭科学技术第４９卷安设第１组应力测点，之后每隔４０ｍ安设一组，共布置１０组，主要监测工作面超前３００ｍ范围内回风巷一侧煤体采动应力分布情况；每组垂直于煤壁施工２个ø４４ｍｍ应力钻孔，孔深分别为１１ｍ和１６ｍ，钻孔间距１ｍ㊂当测点与工作面距离小于３０ｍ时开始回撤，随着工作面回采，测点依次前移，直至回采结束㊂测点布置方案如图１２所示㊂收集了３１１２０２工作面２０１９年５月至１１月的回风巷采动应力监测数据，并进行了分析和应用㊂图１２㊀应力在线监测测点布置Ｆｉｇ．１２㊀Ｌａｙｏｕｔｏｆｍｅａｓｕｒｉｎｇｐｏｉｎｔｓｆｏｒｏｎｌｉｎｅｓｔｒｅｓｓｍｏｎｉｔｏｒｉｎｇ工作面震动场数据由ＡＲＡＭＩＳＭ／Ｅ微震监测系统监测得到㊂３１１２０２工作面测站布置情况如图１３所示㊂井下布置４台微震拾震器（编号Ｓ９至Ｓ１２）和６个移动式监测探头（编号Ｔ１９至Ｔ２４），地面布置１台编号为Ａ２矿震测站组成联合监测网，对工作面进行全面监测㊂图１３㊀３１１２０２回采工作面微震监测系统测站布置Ｆｉｇ．１３㊀ＡｒｒａｎｇｅｍｅｎｔｏｆｔｈｅｓｔａｔｉｏｎｏｆｍｉｃｒｏｓｅｉｓｍｉｃｍｏｎｉｔｏｒｉｎｇｓｙｓｔｅｍｉｎＮｏ．３１１２０２ｍｉｎｉｎｇｆａｃｅ选取３１１２０２工作面回采至距离ＤＦ１９断层１０ｍ时，开始揭露ＤＦ１９断层时以及揭露ＤＦ１９断层２９５ｍ时，３个时间节点３１１２０２工作面超前１５０ｍ范围内的冲击危险性评价情况㊂回采至距离ＤＦ１９断层１０ｍ时，计算原岩应力场冲击危险性指数Ｒ，３－１煤层平均采深６２０ｍ，Ｒ１＝３；工作面距离断层１０ｍ，Ｒ２＝４；工作面前方无背斜或向斜，Ｒ３＝１；该区域未发生过冲击地压，Ｒ４＝１㊂根据式（９）计算得到Ｒ＝２．３㊂按照式（１）㊁式（２）计算得到γσ＝２．３，ｇσ＝３．３㊂３１１２０２工作面最大主应力与水平应力比约为１，取ａ１＝ｂ１＝ｃ１＝０．５，根据式（１０）计算得到ＤＲＭ＝４．０㊂同理计算出，揭露断层时ＤＲＭ＝５．０；揭露断层２９５ｍ时ＤＲＭ＝４．０㊂回采至距离ＤＦ１９断层１０ｍ时，Ｒ＝２．３；根据式（３）㊁（４）计算得到γθ＝３．４，ｇθ＝５．０；工作面最大主应力与水平应力比约为１，取ａ２＝ｂ２＝ｃ２＝０．５，根据式（１１）计算得到ＤＲＳ＝５．４㊂同理计算出，揭露断层时ＤＲＳ＝６．５；揭露断层２９５ｍ时ＤＲＳ＝４．５㊂回采至距离ＤＦ１９断层１０ｍ时，根据式（１）㊁式（２）计算得到γσ＝２．３，ｇσ＝３．３；根据式（３）㊁式（４）计算得到γθ＝３．４，ｇθ＝５．０㊂３１１２０２工作面最大主应力与水平应力比约为１，取ａ３＝ｂ３＝ｃ３＝ｄ３＝０．５，根据式（１２）计算得到ＤＭＳ＝７．０㊂同理计算出，揭露断层时ＤＭＳ＝９．２；揭露断层２９５ｍ时ＤＭＳ＝６．２㊂根据式（１３）计算得到，回采至距离ＤＦ１９断层１０ｍ时Ｄ＝１６．４，具有强冲击危险性；揭露断层时Ｄ＝２０．７，具有强冲击危险性；揭露断层２９５ｍ时Ｄ＝１４．７，具有中等冲击危险性㊂４．３㊀评价结果验证与对比依据３１１２０２工作面回采期间超前工作面３００ｍ范围内微震监测数据㊁钻孔应力监测数据平均值验证评价结果㊂在距离ＤＦ１９断层１０ｍ附近，当天微震释放总能量约为１９３００Ｊ，单次最大能量为７０００Ｊ，微震事件２６次；揭露断层时，当天微震释放总能量约为２２３００Ｊ，单次最大能量约为９０００Ｊ，微震事件１７次；揭露断层２９６ｍ附近，当天微震释放总能量约为７７００Ｊ，单次最大能量约为６０００Ｊ，微震事件６次㊂从微震事件能量㊁频次中可以看出冲击危险性降低㊂在距离断层１０ｍ附近㊁揭露断层附近以及揭露断层２９６ｍ附近选取３个煤层钻孔应力测点，３个测点应力监测数据如图１４所示㊂工作面推进过程中煤层应力数值增加，强冲击危险区域应力始终高于中等冲击危险区域㊂微震和煤层钻孔应力监测数据验证了冲击危险性动态评价结果的合理性㊂图１４㊀煤层钻孔应力监测数据平均值Ｆｉｇ．１４㊀Ａｖｅｒａｇｅｖａｌｕｅｓｏｆｓｔｒｅｓｓｍｏｎｉｔｏｒｉｎｇｄａｔａｓｉｎｃｏａｌｓｅａｍ０３１。

ornstein-uhlenbeck noise 原理 -回复Ornstein-Uhlenbeck noise (OU noise) is a stochastic process that is widely used in various fields, including physics, biology, and finance. This noise process is named after physicists Leonard Ornstein and George Eugene Uhlenbeck, who first introduced it in 1930 as a model for the velocity of a Brownian particle subject to friction.To understand the principles behind Ornstein-Uhlenbeck noise, let us first explore the basics of stochastic processes. In mathematics and st atistics, astochastic process is a collection of random variables indexed by a timeparameter. These random variables can represent the behavior of a system over time, where the value of the variable at a given time is uncertain.One common example of a stochastic process is Brownian motion. Brownian motion describes the random movement of particles suspended in a fluid, such as the movement of pollen grains in water. The particle's position changes randomly over time due to the continuous bombardment of water molecules.Ornstein-Uhlenbeck noise, on the other hand, models a system with some level of mean reversion or tendency to return to a central value. This characteristic isoften observed in real-world phenomena, such as the behavior of stock prices or the motion of a pendulum subject to damping.The Ornstein-Uhlenbeck process can be mathematically defined as a stochastic differential equation (SDE). The SDE representing Ornstein-Uhlenbeck noise is as follows:dX(t) = θ(μ- X(t))dt + σdW(t)In this equation, dX(t) represents the infinitesimal change in the value of the variable X at time t. The first term θ(μ- X(t))dt models the mean reversion behavior, where θrepresents the strength of mean reversion and μis themean or central value that the pr ocess tends to move towards. The second term σdW(t) captures the instantaneous random fluctuations, where dW(t) is an increment of standard Wiener process or Brownian motion, and σrepresents the intensity of these fluctuations.By solving the SDE, we can generate sample paths of the Ornstein-Uhlenbeck process. These sample paths exhibit a characteristic behavior where the process tends to revert to its mean value while still allowing for random fluctuations around it. The rate of mean reversion and the in tensity of fluctuations can be adjusted by changing the parameters θand σ.In practical applications,Ornstein-Uhlenbeck noise is often used as a tool for generating synthetic time series data that resemble certainreal-world phenomena. For example, in finance, the motion of stock prices is often modeled as a stochastic process with mean reversion to capture the observed tendency of prices to revert to their long-term average. Ornstein-Uhlenbeck noise provides a useful framework for simulating such proces ses. Similarly, in physics, this noise process has been used to model the motion of charged particles subject to random forces.Furthermore, Ornstein-Uhlenbeck noiseis used as a regularization technique in machine learning and optimization problems. By a dding this noise to the gradient updates during training, it introduces a form of randomness that can help prevent the model from getting stuck in local optima and improve its generalization performance.。

ＪｏｕｒｎａｌｏｆＭｅｃｈａｎｉｃａｌＳｔｒｅｎｇｔｈ２０２３，４５（３）：５１９⁃５２６ＤＯＩ：１０ １６５７９／ｊ．ｉｓｓｎ．１００１ ９６６９ ２０２３ ０３ ００２∗２０２１０９１３收到初稿，２０２１１１２０收到修改稿㊂国家自然科学基金项目（１２１７２３２１），河北省自然科学基金项目（Ａ２０２０２０３００７）资助㊂∗∗曹天笑，男，１９９６年生，河北石家庄人，汉族，燕山大学在读硕士研究生，主要研究方向为非线性磁弹性振动㊂∗∗∗胡宇达（通信作者），男，１９６８年生，黑龙江东宁人，汉族，燕山大学教授，博士研究生导师，主要研究方向为非线性振动和磁弹性力学㊂谐变磁力作用轴向运动铁磁板的磁弹性主共振∗ＭＡＧＮＥＴＯ⁃ＥＬＡＳＴＩＣＰＲＩＮＣＩＰＡＬＲＥＳＯＮＡＮＣＥＯＦＡＸＩＡＬＬＹＭＯＶＩＮＧＦＥＲＲＯＭＡＧＮＥＴＩＣＰＬＡＴＥＵＮＤＥＲＨＡＲＭＯＮＩＣＭＡＧＮＥＴＩＣＦＯＲＣＥ曹天笑∗∗１，２㊀胡宇达∗∗∗１，２（１．燕山大学建筑工程与力学学院，秦皇岛０６６００４）（２．燕山大学河北省重型装备与大型结构力学可靠性重点实验室，秦皇岛０６６００４）ＣＡＯＴｉａｎＸｉａｏ１，２㊀ＨＵＹｕＤａ１，２（１．ＳｃｈｏｏｌｏｆＣｉｖｉｌＥｎｇｉｎｅｅｒｉｎｇａｎｄＭｅｃｈａｎｉｃｓ，ＹａｎｓｈａｎＵｎｉｖｅｒｓｉｔｙ，Ｑｉｎｈｕａｎｇｄａｏ０６６００４，Ｃｈｉｎａ）（２．ＨｅｂｅｉＫｅｙＬａｂｏｒａｔｏｒｙｏｆＭｅｃｈａｎｉｃａｌＲｅｌｉａｂｉｌｉｔｙｆｏｒＨｅａｖｙＥｑｕｉｐｍｅｎｔｓａｎｄＬａｒｇｅＳｔｒｕｃｔｕｒｅｓ，ＹａｎｓｈａｎＵｎｉｖｅｒｓｉｔｙ，Ｑｉｎｈｕａｎｇｄａｏ０６６００４，Ｃｈｉｎａ）摘要㊀基于哈密顿变分原理，建立机械载荷和磁化产生的谐变磁力共同作用下铁磁板的非线性磁弹性振动方程㊂针对两长边简支约束边界条件，利用伽辽金积分法得到变量分离后的横向振动微分方程㊂应用多尺度法和李雅普诺夫稳定性理论求解电磁激发下的１阶主共振问题，得到稳态响应下的幅频方程和定常解的稳定性判据㊂通过算例，得到铁磁板的幅频特性㊁振幅⁃磁场强度幅值和振幅⁃速度的变化曲线图㊂曲线分析结果表明，稳定解部分的振幅随磁场强度幅值的增大而增大；非线性刚度随速度的增大而增大，硬弹簧特性增强，非线性特征更为显著㊂关键词㊀铁磁板㊀轴向运动㊀主共振㊀谐变磁力㊀多尺度法中图分类号㊀Ｏ３２２㊀Ｏ４４２㊀㊀㊀㊀Ａｂｓｔｒａｃｔ㊀ＢａｓｅｄｏｎＨａｍｉｌｔｏｎｐｒｉｎｃｉｐｌｅ，ｔｈｅｎｏｎｌｉｎｅａｒｍａｇｎｅｔｏ⁃ｅｌａｓｔｉｃｖｉｂｒａｔｉｏｎｅｑｕａｔｉｏｎｏｆｆｅｒｒｏｍａｇｎｅｔｉｃｐｌａｔｅｉｓｅｓｔａｂｌｉｓｈｅｄｕｎｄｅｒｔｈｅｃｏｍｂｉｎｅｄａｃｔｉｏｎｏｆｍｅｃｈａｎｉｃａｌｌｏａｄａｎｄｈａｒｍｏｎｉｃｍａｇｎｅｔｉｃｆｏｒｃｅｉｎｄｕｃｅｄｂｙｍａｇｎｅｔｉｚａｔｉｏｎ．Ｆｏｒｔｈｅｂｏｕｎｄａｒｙｃｏｎｄｉｔｉｏｎｓｗｉｔｈｔｗｏｌｏｎｇｓｉｍｐｌｙｓｕｐｐｏｒｔｅｄｅｄｇｅｓ，ｔｈｅｔｒａｎｓｖｅｒｓｅｖｉｂｒａｔｉｏｎｄｉｆｆｅｒｅｎｔｉａｌｅｑｕａｔｉｏｎｓａｆｔｅｒｖａｒｉａｂｌｅｓｓｅｐａｒａｔｉｏｎｉｓｏｂｔａｉｎｅｄｂｙＧａｌｅｒｋｉｎｍｅｔｈｏｄ．Ｔｈｅｍｕｌｔｉ⁃ｓｃａｌｅｍｅｔｈｏｄａｎｄＬｙａｐｕｎｏｖｓｔａｂｉｌｉｔｙｔｈｅｏｒｙａｒｅｕｓｅｄｔｏｓｏｌｖｅｔｈｅｆｉｒｓｔ⁃ｏｒｄｅｒｐｒｉｎｃｉｐａｌｒｅｓｏｎａｎｃｅｐｒｏｂｌｅｍｕｎｄｅｒｅｌｅｃｔｒｏｍａｇｎｅｔｉｃｅｘｃｉｔａｔｉｏｎ，ｔｈｅａｍｐｌｉｔｕｄｅ⁃ｆｒｅｑｕｅｎｃｙｅｑｕａｔｉｏｎａｎｄｔｈｅｓｔａｂｉｌｉｔｙｃｒｉｔｅｒｉｏｎｏｆｔｈｅｓｔｅａｄｙ⁃ｓｔａｔｅｓｏｌｕｔｉｏｎｓａｒｅｏｂｔａｉｎｅｄ．Ｔｈｒｏｕｇｈｎｕｍｅｒｉｃａｌｅｘａｍｐｌｅｓ，ｔｈｅａｍｐｌｉｔｕｄｅ⁃ｆｒｅｑｕｅｎｃｙｃｈａｒａｃｔｅｒｉｓｔｉｃｃｕｒｖｅｓ，ｔｈｅａｍｐｌｉｔｕｄｅ⁃ｍａｇｎｅｔｉｃｆｉｅｌｄｉｎｔｅｎｓｉｔｙｃｕｒｖｅｓａｎｄｔｈｅａｍｐｌｉｔｕｄｅ⁃ｖｅｌｏｃｉｔｙｖａｒｉａｔｉｏｎｃｕｒｖｅｓｏｆｔｈｅｆｅｒｒｏｍａｇｎｅｔｉｃｐｌａｔｅａｒｅｏｂｔａｉｎｅｄ．Ｔｈｅｒｅｓｕｌｔｓｓｈｏｗｔｈａｔｔｈｅａｍｐｌｉｔｕｄｅｏｆｓｔａｂｌｅｓｏｌｕｔｉｏｎｓｉｎｃｒｅａｓｅｓｗｉｔｈｔｈｅｉｎｃｒｅａｓｅｏｆｍａｇｎｅｔｉｃｆｉｅｌｄｉｎｔｅｎｓｉｔｙａｍｐｌｉｔｕｄｅ．Ｔｈｅｎｏｎｌｉｎｅａｒｓｔｉｆｆｎｅｓｓｉｎｃｒｅａｓｅｓｗｉｔｈｔｈｅｉｎｃｒｅａｓｅｏｆｔｈｅｖｅｌｏｃｉｔｙ，ｔｈｅｈａｒｄ⁃ｓｐｒｉｎｇｃｈａｒａｃｔｅｒｉｓｔｉｃｓａｒｅｅｎｈａｎｃｅｄａｎｄｔｈｅｎｏｎｌｉｎｅａｒｃｈａｒａｃｔｅｒｉｓｔｉｃｓａｒｅｍｏｒｅｓｉｇｎｉｆｉｃａｎｔ．Ｋｅｙｗｏｒｄｓ㊀Ｆｅｒｒｏｍａｇｎｅｔｉｃｐｌａｔｅ；Ａｘｉａｌｌｙｍｏｖｉｎｇ；Ｐｒｉｎｃｉｐａｌｒｅｓｏｎａｎｃｅ；Ｈａｒｍｏｎｉｃｍａｇｎｅｔｉｃｆｏｒｃｅ；Ｍｕｌｔｉ⁃ｓｃａｌｅｍｅｔｈｏｄＣｏｒｒｅｓｐｏｎｄｉｎｇａｕｔｈｏｒ：ＨＵＹｕＤａ，Ｅ⁃ｍａｉｌ：ｈｕｙｕｄａ０３＠１６３．ｃｏｍ，Ｔｅｌ：＋８６⁃３３５⁃８０５７１０１，Ｆａｘ：＋８６⁃３３５⁃８０５７１０１ＴｈｅｐｒｏｊｅｃｔｓｕｐｐｏｒｔｅｄｂｙｔｈｅＮａｔｉｏｎａｌＮａｔｕｒａｌＳｃｉｅｎｃｅＦｏｕｎｄａｔｉｏｎｏｆＣｈｉｎａ（Ｎｏ．１２１７２３２１），ａｎｄｔｈｅＮａｔｕｒａｌＳｃｉｅｎｃｅＦｏｕｎｄａｔｉｏｎｏｆＨｅｂｅｉＰｒｏｖｉｎｃｅ（Ｎｏ．Ａ２０２０２０３００７）．Ｍａｎｕｓｃｒｉｐｔｒｅｃｅｉｖｅｄ２０２１０９１３，ｉｎｒｅｖｉｓｅｄｆｏｒｍ２０２１１１２０．０㊀引言㊀㊀多场环境下轴向运动结构具有重要工程应用背景，如磁悬浮运输车体悬力结构，高精轴向运动板带的多因素轧制调质工艺和车辆ＣＶＴ传动钢带等㊂当铁磁结构在多重耦合场中工作时，会不可避免地受到激励，发生复杂的振动行为㊂利用铁磁材料对磁场极其敏感的特性，通过外加磁场可实现对铁磁结构振动特㊀５２０㊀机㊀㊀械㊀㊀强㊀㊀度２０２３年㊀性的电磁控制㊂因此，对处于磁场中轴向运动结构磁弹性振动特性的研究具有重要的理论和应用价值㊂针对轴向运动结构，ＬＩＮＣＣ［１］研究了二维轴向运动平板的稳定性问题㊂ＺＨＯＵＹＦ等［２］分析了物理参数和速度对厚度抛物线变化下轴向运动黏弹性矩形板振动特性的影响㊂ＭＯＨＡＭＡＤＩＡ等［３］研究了受黏性结构阻尼影响的轴向运动简支薄圆柱壳的线性自由振动问题㊂ＬＩＨＹ等［４］分析了流固耦合影响下浸水轴向运动板的动力学特性和稳定性㊂ＦＡＲＳＨＢＡＦＺＲ等［５］对具有中间非线性支撑的轴向运动简支黏弹性梁进行了非线性振动分析㊂ＨＵＹＤ等［６］应用改进多尺度法对周期外载作用下轴向运动导电条形板的强非线性振动及混沌运动问题进行了研究㊂黄建亮等［７⁃８］采用多元Ｌ⁃Ｐ方法研究了轴向运动体系横向内共振和联合共振问题㊂殷振坤等［９］应用增量谐波平衡法研究了轴向运动薄板横向非线性振动特性及其稳定性㊂李中华等［１０］研究了轴向运动黏弹性夹层板的模态耦合横向振动问题㊂对于磁弹耦合场中的振动问题，李晶等［１１］分析了不同情况下磁场强度对矩形导电薄板内共振特性的影响㊂ＨＵＹＤ等［１２］对磁场中的轴向运动载流梁的１ʒ３主内共振问题进行了研究㊂高原文等［１３⁃１４］研究了铁磁梁式板在横向磁场和脉冲磁场作用下的磁弹性动力响应特征和动力失稳现象㊂徐浩然等［１５］研究了处于平行共轴三线圈和球形载流线圈产生磁场中的导电圆板的固有振动问题㊂ＨＡＳＡＮＹＡＮＤ等［１６］分析了磁场和导电率对轴向磁场中有限导电板振动特性的影响㊂ＷＡＮＧＸ等［１７］研究了具有磁弹性相互作用和磁阻尼的铁磁梁式板在横向磁场中的动态稳定性问题㊂ＨＵＹＤ等［１８］研究了磁场中导电薄板的磁弹性组合共振和次谐波共振问题㊂将轴向运动铁磁板置于横向谐变磁场中，不仅需要考虑轴向速度对薄板振动时非线性特征的影响，同时，由于铁磁材料会被磁化产生谐变磁力作用于铁磁板，因此会产生由磁场环境激发的振动行为㊂本文研究谐变磁力作用下轴向运动铁磁板主共振问题，分析几何和物理参量对振幅和非线性振动特性的影响㊂１㊀横向磁场中轴向运动铁磁板的振动方程㊀㊀考虑横向磁场环境（Ｈｎ１为上表面磁场强度，Ｈｎ２为下表面磁场强度）中，沿着ｘ方向以速度Ｖ０ｘ做轴向匀速运动，并受到边界面内拉力Ｆ０ｘ和均布横向机械载荷Ｐｚ作用的各向同性铁磁条形薄板㊂如图１所示，板长为ｌ；板厚为ｈ；薄板的弹性模量为Ｅ；泊松比为μ；密度为ρ㊂１ １㊀电磁力㊀㊀各向同性线性软铁磁介质所受到的磁体力和边界图１㊀磁场中的轴向运动铁磁板Ｆｉｇ．１㊀Ａｘｉａｌｌｙｍｏｖｉｎｇｆｅｒｒｏｍａｇｎｅｔｉｃｐｌａｔｅｉｎｍａｇｎｅｔｉｃｆｉｅｌｄ磁力［１９］分别为磁体力㊀ｆｅｍ＝▽Ｂ()Ｍ＝μ０μｒχｍ▽Ｈ()２（１）边界磁力㊀Ｆｅｍ＝－μ０χｍ（μｒ＋１）２Ｈｔ()２ｎ（２）式中，Ｂ为磁感应强度矢量，Ｂ＝μ０μｒＨ；Ｈ为磁场强度矢量；Ｍ为磁化矢量，Ｍ＝χｍＨ；μ０为真空磁导率；μｒ为相对磁导率；χｍ为材料的磁化率，χｍ＝μｒ－１；ｎ为铁磁介质表面的单位法向矢量；▽＝∂∂ｘｉ＋∂∂ｚｋ，ｉ和ｋ分别为ｘ方向和ｚ方向上的单位矢量㊂将磁体力和边界磁力向中面简化后得到等效横向磁力为Ｆｚ＝ʏｈ２－ｈ２ｆｅｍｚ（ｘ，ｚ）ｄｚ＋Ｆｅｍｚ（ｘ，ｈ２）－Ｆｅｍｚ（ｘ，－ｈ２）＝μ０μｒχｍ２Ｈｎ（ｘ，ｈ２）éëêêùûúú２－Ｈｎ（ｘ，－ｈ２）éëêêùûúú２{}－μ０χｍ２Ｈｔ（ｘ，ｈ２）éëêêùûúú２－Ｈｔ（ｘ，－ｈ２）éëêêùûúú２{}（３）式中，Ｈｎ为铁磁板表面法线方向磁场强度；Ｈｔ为铁磁板表面切线方向磁场强度㊂设薄板内部磁场沿ｚ轴线性分布［２０］，则Ｂ０ｚ＝μ０μｒＨｎ１＋Ｈｎ２２＋Ｈｎ２－Ｈｎ１ｈｚæèçöø÷（４）㊀㊀轴向运动铁磁板在横向磁场中所受洛伦兹电磁力为ｆｘ＝ＪｙＢ０ｚ＝－σ０Ｂ２０ｚＶ０ｘ－ｚｄｄｔ∂ｗ∂ｘæèçöø÷éëêêùûúú（５）式中，Ｊｙ为铁磁板在磁场中运动所产生的电流密度；σ０为电导率；Ｖ０ｘ为轴向运动速度；ｗ为中面横向位移㊂将式（４）代入式（５）并对ｚ沿板厚方向积分得电磁力矩为ｍｘ＝ʏｈ２－ｈ２ｆｘｚｄｚ＝σ０μ２０μ２ｒｈ３ｄｄｔ∂ｗ∂ｘæèçöø÷㊃㊀㊀（Ｈｎ１＋Ｈｎ２）２４８＋éëêê（Ｈｎ２－Ｈｎ１）２８０ùûúú＋㊀㊀σ０μ２０μ２ｒｈ２Ｖ０ｘ（Ｈ２ｎ２－Ｈ２ｎ１）１２（６）㊀第４５卷第３期曹天笑等：谐变磁力作用轴向运动铁磁板的磁弹性主共振５２１㊀㊀１ ２㊀动能和势能㊀㊀对于沿ｙ方向上无限长的铁磁板，做仅随坐标ｘ变化的横向振动时，系统的动能表示为Ｔ＝∬ρ２Ｖ２ｄｘｄｚ＝ʏＶ２０ｘ＋ｄｗｄｔæèçöø÷２éëêêùûúúｄｘ＋ρｈ３２４ʏｄｄｔ∂ｗ∂ｘæèçöø÷éëêêùûúú２ｄｘ（７）式中，Ｖ为板内任一点的绝对速度矢量；ｔ为时间变量；ｄ／ｄｔ＝Ｖ０ｘ∂／∂ｘ＋∂／∂ｔ㊂轴向运动铁磁板发生变形时，势能Ｕ包括弯曲应变势能㊁中面应变势能和边界面内拉力Ｆ０ｘ引起的应变势能㊂根据弹性薄板理论，势能表达式为Ｕ＝ＤＭ２ʏ∂２ｗ∂ｘ２æèçöø÷２ｄｘ＋１２ʏＮｘεｘｄｘ＋ʏＦ０ｘεｘｄｘ（８）式中，ＤＭ为弯曲刚度，ＤＭ＝Ｅｈ３／［１２（１－μ２）］；ＤＮ为拉伸强度，ＤＮ＝Ｅｈ／（１－μ２）；Ｎｘ为中面应力，Ｎｘ＝ＤＮ（εｘ＋μεｙ）；εｘ㊁εｙ均为中面应变分量，εｘ＝（∂ｗ／∂ｘ）２／２，εｙ＝（∂ｗ／∂ｙ）２／２㊂１ ３㊀磁弹性横向振动方程㊀㊀根据哈密顿变分原理，有ʏｔ１ｔ２（δＴ－δＵ＋δＷＰ＋δＷ）ｄｔ＝０（９）式中，ｔ１㊁ｔ２分别为两固定时刻；δＷＰ为机械载荷Ｐｚ所做虚功；δＷ为电磁力虚功㊂将动能㊁势能变分后和电磁力虚功代入式（９），整理得到忽略面内位移的轴向运动铁磁板的非线性磁弹性横向振动方程为－ＤＭ∂４ｗ∂ｘ４＋３２ＤＮ∂ｗ∂ｘæèçöø÷２∂２ｗ∂ｘ２＋Ｆ０ｘ∂２ｗ∂ｘ２＋Ｆｚ＋㊀㊀∂ｍｘ∂ｘ＋Ｐｚ＝ρｈＶ２０ｘ∂２ｗ∂ｘ２＋２Ｖ０ｘ∂２ｗ∂ｘ∂ｔ＋∂２ｗ∂ｔ２æèçöø÷－㊀㊀㊀ρｈ３１２Ｖ２０ｘ∂４ｗ∂ｘ４＋æèç２Ｖ０ｘ∂４ｗ∂ｘ３∂ｔ＋∂４ｗ∂ｘ２∂ｔ２öø÷（１０）２㊀轴向运动铁磁板的主共振分析２ １㊀振动方程的伽辽金离散㊀㊀对于两长边简支的轴向运动铁磁板，其满足边界条件为ｗ＝０，∂２ｗ∂ｘ２＝０，ｘ＝０ｗ＝０，∂２ｗ∂ｘ２＝０，ｘ＝ｌìîíïïïï（１１）㊀㊀在考虑３阶模态的情况下，设满足边界条件式（１１）的位移函数为ｗ＝ð３ｎ＝１ｐｎ（ｔ）Ｗｎ（ｘ）＝ð３ｎ＝１ｐｎ（ｔ）ｓｉｎｎπｘｌ（１２）㊀㊀设铁磁板上㊁下表面的磁场强度分别为Ｈｎ１＝Ｈ０ｃｏｓ（Ω１ｔ），Ｈｎ２＝Ｈ０ｓｉｎ（Ω１ｔ）（１３）式中，Ｈ０为磁场强度幅值；Ω１为磁场强度频率㊂将式（１３）代入式（３）中，且仅考虑横向磁场，得到横向谐变磁力为Ｆｚ＝μ０μｒχｍ２Ｈ２０ｃｏｓ（２Ω１ｔ）（１４）同时，机械载荷按谐变规律给定：Ｐｚ＝Ｐ０ｃｏｓ（Ω２ｔ）（１５）式中，Ｐ０为载荷幅值；Ω２为载荷频率㊂令磁场强度频率和机械载荷频率关系为２Ω１＝Ω２＝Ω，并将式（１４）和式（１５）代入式（１０）中，进行伽辽金积分，推导得无量纲化后的横向振动微分方程为ｑ㊆１＋ω２１ｑ１＝η２１Ｈ２０１２０ｓｉｎ（Ω－τ）＋Ｈ２０３０éëêêùûúúｑ２＋φ２１ｑ㊃２＋㊀㊀㊀ξ１Ｈ２０１２０ｓｉｎ（Ω－τ）éëêê＋Ｈ２０３０ùûúúｑ㊃１＋α１Ｓ１１ｑ３１＋Ｓ４１ｑ１ｑ２２＋(㊀㊀㊀Ｓ５１ｑ１ｑ２３＋Ｓ７１ｑ２１ｑ３＋Ｓ８１ｑ２２ｑ３)＋ｆ１ｃｏｓ（Ω－τ）ｑ㊆２＋ω２２ｑ２＝η１２Ｈ２０１２０ｓｉｎ（Ω－τ）＋Ｈ２０３０éëêêùûúúｑ１＋φ１２ｑ㊃１＋㊀㊀㊀φ３２ｑ㊃３＋η３２Ｈ２０１２０ｓｉｎ（Ω－τ）＋éëêêＨ２０３０ùûúúｑ３＋㊀㊀㊀ξ２Ｈ２０１２０ｓｉｎ（Ω－τ）＋Ｈ２０３０éëêêùûúúｑ㊃２＋α２Ｓ２２ｑ３２＋(㊀㊀㊀Ｓ６２ｑ２１ｑ２＋Ｓ９２ｑ２ｑ２３＋Ｓ０２ｑ１ｑ２ｑ３)ｑ㊆３＋ω２３ｑ３＝η２３Ｈ２０１２０ｓｉｎ（Ω－τ）＋Ｈ２０３０éëêêùûúúｑ２＋φ２３ｑ㊃２＋㊀㊀㊀ξ３Ｈ２０１２０ｓｉｎ（Ω－τ）＋éëêêＨ２０３０ùûúúｑ㊃３＋α３Ｓ１３ｑ３１＋Ｓ３３ｑ３３＋(㊀㊀㊀Ｓ４３ｑ１ｑ２２＋Ｓ７３ｑ２１ｑ３＋Ｓ８３ｑ２２ｑ３)＋ｆ３ｃｏｓ（Ω－τ）ìîíïïïïïïïïïïïïïïïïïïïïïïïïïïïï（１６）式中ｑｎ＝ｐｎｈ，Ω＝ΩωＮ，τ＝ωＮｔ，ηｎｉ＝σ０μ２０μ２ｒｈ３Ｖ０ｘＤｎｉＡｉｉω２Ｎ，ξｉ＝σ０μ２０μ２ｒｈ３ＣｉｉＡｉｉωＮ，φｎｉ＝ρｈ３Ｖ０ｘＤｎｉ－１２ρｈＶ０ｘＢｎｉ６ＡｉｉωＮ，αｉ＝３ＤＮｈ２２Ａｉｉω２Ｎ，ｆｉ＝Ｇｉω２ＮＡｉｉｈμ０μｒχｍ２Ｈ２０＋Ｐ０æèçöø÷，ωＮ＝３ｇ１ｇ２ｇ３，ω１＝ｇ１ωＮ，ω２＝ｇ２ωＮ，ω３＝ｇ３ωＮ，ｇ２１＝１２ＤＭＥ１１－１２Ｆ０ｘＣ１１＋１２ρｈＶ２０ｘＣ１１－ρｈ３Ｖ２０ｘＥ１１１２Ａ１１，ｇ２２＝１２ＤＭＥ２２－１２Ｆ０ｘＣ２２＋１２ρｈＶ２０ｘＣ２２－ρｈ３Ｖ２０ｘＥ２２１２Ａ２２，㊀５２２㊀机㊀㊀械㊀㊀强㊀㊀度２０２３年㊀ｇ２３＝１２ＤＭＥ３３－１２Ｆ０ｘＣ３３＋１２ρｈＶ２０ｘＣ３３－ρｈ３Ｖ２０ｘＥ３３１２Ａ３３㊂式中，Ａｎｉ㊁Ｂｎｉ㊁Ｃｎｉ㊁Ｄｎｉ㊁Ｅｎｉ㊁Ｇｉ㊁Ｓｂｉ（ｎ＝１，２，３；ｉ＝１，２，３；ｂ＝０，１， ，９）均为积分式㊂２ ２㊀多尺度法求解㊀㊀研究系统主共振问题，经验证，系统为弱非线性，引入小参数ε，式（１６）可写为ｑ㊆１＋ω２１ｑ１＝εη－２１Ｈ２０１２０ｓｉｎ（Ω－τ）＋Ｈ２０３０éëêêùûúúｑ２＋εφ－２１ｑ㊃２＋㊀㊀㊀εξ－１Ｈ２０１２０ｓｉｎ（Ω－τ）éëêê＋Ｈ２０３０ùûúúｑ㊃１＋εα－１Ｓ１１ｑ３１＋Ｓ４１ｑ１ｑ２２＋(㊀㊀㊀Ｓ５１ｑ１ｑ２３＋Ｓ７１ｑ２１ｑ３＋Ｓ８１ｑ２２ｑ３)＋εｆ－１ｃｏｓ（Ω－τ）ｑ㊆２＋ω２２ｑ２＝εη－１２Ｈ２０１２０ｓｉｎ（Ω－τ）＋Ｈ２０３０éëêêùûúúｑ１＋εφ－１２ｑ㊃１＋㊀㊀εφ－３２ｑ㊃３＋εη－３２Ｈ２０１２０ｓｉｎ（Ω－τ）＋éëêêＨ２０３０ùûúúｑ３＋㊀㊀εξ－２Ｈ２０１２０ｓｉｎ（Ω－τ）＋Ｈ２０３０éëêêùûúúｑ㊃２＋εα－２Ｓ２２ｑ３２＋(㊀㊀Ｓ６２ｑ２１ｑ２＋Ｓ９２ｑ２ｑ２３＋Ｓ０２ｑ１ｑ２ｑ３)ｑ㊆３＋ω２３ｑ３＝εη－２３Ｈ２０１２０ｓｉｎ（Ω－τ）＋Ｈ２０３０éëêêùûúúｑ２＋εφ－２３ｑ㊃２＋㊀㊀εξ－３Ｈ２０１２０ｓｉｎ（Ω－τ）＋éëêêＨ２０３０ùûúúｑ㊃３＋εα－３Ｓ１３ｑ３１＋Ｓ３３ｑ３３＋(㊀㊀Ｓ４３ｑ１ｑ２２＋Ｓ７３ｑ２１ｑ３＋Ｓ８３ｑ２２ｑ３)＋εｆ－３ｃｏｓ（Ω－τ）ìîíïïïïïïïïïïïïïïïïïïïïïïïïïïï（１７）式中，η－ｎｉ＝ηｎｉε，ξ－ｉ＝ξｉε，φ－ｎｉ＝φｎｉε，α－ｉ＝αｉε，ｆ－ｉ＝ｆｉε（ｎ＝１，２，３；ｉ＝１，２，３）㊂为了定量表示主共振发生时激励力频率Ω－与第１阶固有频率ω１之间的接近程度，引入调谐参数σ，并令Ω－＝ω１＋εσ（１８）㊀㊀应用多尺度法进行１阶近似求解，近似解设为ｑ１（τ，ε）＝ｑ１１（Ｔ０，Ｔ１）＋εｑ１２（Ｔ０，Ｔ１）ｑ２（τ，ε）＝ｑ２１（Ｔ０，Ｔ１）＋εｑ２２（Ｔ０，Ｔ１）ｑ３（τ，ε）＝ｑ３１（Ｔ０，Ｔ１）＋εｑ３２（Ｔ０，Ｔ１）ìîíïïïï（１９）其中，Ｔｎ＝εｎτ㊂将式（１９）代入式（１７）中，展开后令ε的同次幂项系数相等，将ε０同次幂项系数方程组的解写为复数形式：ｑ１１＝Ａ１（Ｔ１）ｅｉω１Ｔ０＋Ａ－１（Ｔ１）ｅ－ｉω１Ｔ０ｑ２１＝Ａ２（Ｔ１）ｅｉω２Ｔ０＋Ａ－２（Ｔ１）ｅ－ｉω２Ｔ０ｑ３１＝Ａ３（Ｔ１）ｅｉω３Ｔ０＋Ａ－３（Ｔ１）ｅ－ｉω３Ｔ０ìîíïïïï（２０）㊀㊀将式（１８）和式（２０）代入ε１同次幂项系数方程组中，可得消除久期项的条件为－２ω１ｉＤ１Ａ１＋ξ－１ω１ｉＡ１Ｈ２０３０＋３α－１Ｓ１１Ａ２１Ａ－１＋㊀㊀２α－１Ｓ４１Ａ１Ａ２Ａ－２＋２α－１Ｓ５１Ａ１Ａ３Ａ－３＋ｆ－１２ｅｉσＴ１＝０－２ω２ｉＤ１Ａ２＋ξ－２ω２ｉＡ２Ｈ２０３０＋３α－２Ｓ２２Ａ２２Ａ－２＋㊀㊀２α－２Ｓ６２Ａ２Ａ１Ａ－１＋２α－２Ｓ９２Ａ２Ａ３Ａ－３＝０－２ω３ｉＤ１Ａ３＋ξ－３ω３ｉＡ３Ｈ２０３０＋３α－３Ｓ３３Ａ２３Ａ－３＋㊀㊀２α－３Ｓ７３Ａ３Ａ１Ａ－１＋２α－３Ｓ８３Ａ３Ａ２Ａ－２＝０（２１）ìîíïïïïïïïïïïïïïïïï式中，Ｄ１＝∂∂Ｔ１㊂将复函数Ａ写为Ａｋ（Ｔ１）＝１２ａｋｅｉθｋ（ｋ＝１，２，３），代入式（２１）中，并令γ１＝σＴ１－θ１，分离虚部和实部得ａ㊃１＝Ｈ２０６０ξ－１ａ１＋ｆ－１２ω１ｓｉｎγ１ａ１γ㊃１＝ａ１σ＋３α－１８ω１Ｓ１１ａ３１＋α－１４ω１Ｓ４１ａ１ａ２２＋㊀㊀㊀α－１４ω１Ｓ５１ａ１ａ２３＋ｆ－１２ω１ｃｏｓγ１ａ㊃２＝Ｈ２０６０ξ－２ａ２ａ２θ㊃２＝－３α－２８ω２Ｓ２２ａ３２－α－２４ω２Ｓ６２ａ２ａ２１－α－２４ω２Ｓ９２ａ２ａ２３ａ㊃３＝Ｈ２０６０ξ－３ａ３ａ３θ㊃３＝－３α－３８ω３Ｓ３３ａ３３－α－３４ω３Ｓ７３ａ３ａ２１－α－３４ω３Ｓ８３ａ３ａ２２ìîíïïïïïïïïïïïïïïïïïïïïïï（２２）㊀㊀由于ξ－２＜０㊁ξ－３＜０，可见ａ２㊁ａ３都将衰减，共振不被激发㊂对于系统的稳态响应，令ａ㊃１＝０，γ㊃１＝０，得到１阶主共振幅频响应式为Ｈ４０３６００ξ－２１ａ２１＋ａ１σ＋３α－１８ω１Ｓ１１ａ３１æèçöø÷２＝ｆ－２１４ω２１（２３）２ ３㊀稳定性分析㊀㊀主共振发生时，分析系统稳态运动中解的稳定性条件，设ａ１＝ａ１０＋ａ１１，γ１＝γ１０＋γ１１（２４）式中，ａ１０㊁γ１０均为系统的稳态解；ａ１１㊁γ１１均为小的扰动量㊂㊀第４５卷第３期曹天笑等：谐变磁力作用轴向运动铁磁板的磁弹性主共振５２３㊀㊀将式（２４）代入式（２２），根据Ｌｙａｐｕｎｏｖ稳定性理论得到判断稳态解稳定性的本征方程为λ２＋ｃ１１λ＋ｃ１２＝０（２５）式中，ｃ１１＝－Ｈ２０３０ξ－１，ｃ１２＝σ２＋Ｈ４０３６００ξ－２１＋３α－１２ω１σＳ１１ａ２１０＋２７α－２１６４ω２１Ｓ２１１ａ４１０㊂根据Ｒｏｕｔｈ⁃Ｈｕｒｗｉｔｚ判据，在ξ－１＜０情况下可得系统稳态解稳定的充要条件为ｃ１２＞０（２６）３　算例分析㊀㊀针对处于交变磁场中马氏体钢材料制成的轴向运动铁磁板进行算例分析㊂主要物理参数为：密度ρ＝７８００ｋｇ／ｍ３；泊松比μ＝０ ３；材料电导率σ０＝２ ３ˑ１０６（Ω㊃ｍ）－１；相对磁导率μｒ＝１０００；弹性模量Ｅ＝２００ＧＰａ；板长ｌ＝０ ５ｍ；面内拉力Ｆ０ｘ＝１０ｋＮ／ｍ㊂图２（ａ）为不同板厚下轴向运动铁磁板的幅频特性曲线图㊂图２（ａ）中曲线表明，共振区域（εσʈ０）幅值明显增大㊂图２（ａ）中点线代表非稳定解，实线代表稳定解（下同）㊂当调谐值一定时，随着板厚的增加，稳定解部分的共振幅值随之减小；脊骨线附近曲线随着板厚的增加逐渐内缩㊂同时，调谐值由负数增大到一定值时，共振幅值出现多值现象㊂点画线所包夹区域（ｃ１２＜０）为非稳定解区域，即板厚变化时，幅频特性曲线中非稳定解部分所在区域㊂图２（ｂ）为拾取不同厚度下振幅出现多值解的临界调谐值所绘的调谐值⁃厚度多值解分岔点曲线㊂该曲线表明了在不同板厚下开始出现多值解时的调谐值㊂图２（ｂ）中分岔点曲线下方区域对应单值解，上方区域对应多值解㊂图２㊀不同板厚下幅频特性曲线及多值解分岔点曲线Ｆｉｇ．２㊀Ａｍｐｌｉｔｕｄｅ⁃ｆｒｅｑｕｅｎｃｙｃｈａｒａｃｔｅｒｉｓｔｉｃｃｕｒｖｅｓａｔｄｉｆｆｅｒｅｎｔｐｌａｔｅｔｈｉｃｋｎｅｓｓｅｓａｎｄｂｉｆｕｒｃａｔｉｏｎｐｏｉｎｔｓｏｆｍｕｌｔｉ⁃ｖａｌｕｅｄｓｏｌｕｔｉｏｎｓ图３㊀不同物理参量下幅频特性曲线Ｆｉｇ．３㊀Ａｍｐｌｉｔｕｄｅ⁃ｆｒｅｑｕｅｎｃｙｃｈａｒａｃｔｅｒｉｓｔｉｃｃｕｒｖｅｓａｔｄｉｆｆｅｒｅｎｔｐｈｙｓｉｃａｌｐａｒａｍｅｔｅｒｓ㊀㊀图３所示为板厚ｈ＝０ ００９ｍ时不同物理参量下的幅频特性曲线㊂激发主共振的激励力由磁化力Ｆｚ和机械载荷Ｐｚ两部分组成，并且磁化力幅值与磁场强度幅值的平方Ｈ２０成正比㊂因此，在图３（ａ）和图３（ｂ）中，当调谐值一定时，随着磁场强度幅值Ｈ０和机械载荷幅值Ｐ０的增大，共振幅值会随之增大，并且两图中的多值解分岔点向右上方移动㊂在图３（ｃ）中，速度的增大使铁磁板的非线性刚度增大和硬弹簧特性增强，进而幅频特性曲线向右弯曲程度加深，使得上支曲线出现交点，在交点之前振幅随着速度的增大而增大，在交点之后振幅随着速度的增大而减小；在下支曲线的稳定解区域中，振幅随着速度的增大而增大㊂图４（ａ）为不同调谐值下振幅－磁场强度幅值曲线图㊂该曲线图由类椭圆闭合曲线和上支凹型曲线组㊀５２４㊀机㊀㊀械㊀㊀强㊀㊀度２０２３年㊀成，曲线关于磁场强度幅值（Ｈ０＝０）对称㊂椭圆闭合曲线下半部分随着调谐值的增加而减小，上支凹型曲线随着调谐值的增大而增大㊂磁场强度幅值为零（Ｈ０＝０）时所对应振幅不为零，是由于此时共振由机械载荷Ｐｚ激发的㊂共振幅值随着磁场强度幅值的增大，将由多值解转化为单值解㊂多值现象消失的原因在于随着磁场强度幅值的增大，幅频特性曲线图中的脊骨线附近曲线发生外拓，使得调谐值所对应的振幅由多值解转向单值解㊂由于非稳定解只出现于类椭圆闭合曲线部分，在图４（ｂ）中，对类椭圆闭合曲线部分单独分析，点画线为拾取不同调谐值下振幅⁃磁场强度幅值多值解分岔点所绘，其上部区域（ｃ１２＜０）为非稳定解区域，即调谐值变化时，类椭圆闭合曲线中非稳定解部分所在区域㊂图４（ｃ）为拾取不同调谐值下振幅变为单值解的临界磁场强度幅值所绘的磁场强度幅值⁃调谐值多值解分岔点曲线图㊂该曲线表明在不同调谐值下多值现象消失时的磁场强度幅值㊂图４（ｃ）中分岔点曲线下方区域对应单值解，上方区域对应多值解㊂㊀㊀图５为调谐值εσ＝０ ０５时不同几何和物理参量下的振幅⁃磁场强度幅值曲线图㊂由图５（ａ）和图５（ｂ）可知，稳定解部分的共振幅值随着板厚的减小和机械载荷幅值的增大而增大㊂在图５（ｃ）中，随着速度的增大，椭圆闭合曲线稳定解部分共振幅值增大，上支凹型曲线共振幅值减小㊂图６（ａ）为不同调谐值下振幅⁃速度曲线图㊂该曲线图由呈子弹状的下支曲线和上支下凹型曲线两部分组成㊂该曲线图表明，随着调谐值的减小，上支曲线的振幅减小而下支曲线稳定解部分的振幅增大，并且整个下支曲线向内收缩㊂同时，共振幅值随着速度的增大，将由多值解转化为单值解㊂多值现象消失的原因是随着速度的增大，幅频特性曲线图中的脊骨线附近曲线向右弯曲程度加深，使得调谐值所对应的振幅由多值解转向单值解㊂由于非稳定解只出现于下支曲线部分，因此，对图６（ｂ）中下支曲线部分单独分析，点画线为拾取不同调谐值下振幅⁃速度多值解分岔点所绘，其上部区域（ｃ１２＜０）为非稳定解区域，即调谐值变化时，下支曲线中非稳定解部分所在区域㊂图４㊀不同调谐值下振幅⁃磁场强度幅值曲线及多值解分岔点曲线Ｆｉｇ．４㊀Ａｍｐｌｉｔｕｄｅ⁃ｍａｇｎｅｔｉｃｉｎｔｅｎｓｉｔｙａｍｐｌｉｔｕｄｅｃｕｒｖｅｓａｔｄｉｆｆｅｒｅｎｔｔｕｎｉｎｇｖａｌｕｅｓａｎｄｂｉｆｕｒｃａｔｉｏｎｐｏｉｎｔｓｃｕｒｖｅｓｏｆｍｕｌｔｉ⁃ｖａｌｕｅｄｓｏｌｕｔｉｏｎｓ图５㊀不同几何和物理参量下的振幅⁃磁场强度幅值曲线Ｆｉｇ．５㊀Ａｍｐｌｉｔｕｄｅ⁃ｍａｇｎｅｔｉｃｉｎｔｅｎｓｉｔｙａｍｐｌｉｔｕｄｅｃｕｒｖｅｓａｔｄｉｆｆｅｒｅｎｔｇｅｏｍｅｔｒｉｃａｎｄｐｈｙｓｉｃａｌｐａｒａｍｅｔｅｒｓ㊀第４５卷第３期曹天笑等：谐变磁力作用轴向运动铁磁板的磁弹性主共振５２５㊀㊀㊀㊀图６（ｃ）所示为拾取不同调谐值下振幅转换为单值解时的临界速度所绘的速度⁃调谐值多值解分岔点曲线㊂图６（ｃ）中分岔点曲线下方区域对应多值解，上方区域对应单值解㊂图７为调谐值εσ＝０ ０７时不同几何和物理参量下振幅⁃速度曲线图㊂在图７（ａ）中，稳定解部分的共振幅值随着板厚的增大而减小㊂在图７（ｂ）和图７（ｃ）中，稳定解部分的共振幅值随着磁场强度幅值和机械载荷幅值的增大而增大㊂图６㊀不同调谐值下振幅⁃速度曲线及多值解分岔点曲线Ｆｉｇ．６㊀Ａｍｐｌｉｔｕｄｅ⁃ｖｅｌｏｃｉｔｙｃｕｒｖｅｓａｔｄｉｆｆｅｒｅｎｔｔｕｎｉｎｇｖａｌｕｅｓａｎｄｂｉｆｕｒｃａｔｉｏｎｐｏｉｎｔｓｃｕｒｖｅｓｏｆｍｕｌｔｉ⁃ｖａｌｕｅｄｓｏｌｕｔｉｏｎｓ图７㊀不同几何和物理参量下振幅⁃速度曲线Ｆｉｇ．７㊀Ａｍｐｌｉｔｕｄｅ⁃ｖｅｌｏｃｉｔｙｃｕｒｖｅｓａｔｄｉｆｆｅｒｅｎｔｇｅｏｍｅｔｒｉｃａｎｄｐｈｙｓｉｃａｌｐａｒａｍｅｔｅｒｓ㊀㊀选取参数ｈ＝０ ０１ｍ，Ｖ０ｘ＝２０ｍ／ｓ，Ｐ０＝１ｋＮ／ｍ２，Ｈ０＝１００Ａ／ｍ，绘制幅频特性曲线［图８（ａ）］，由曲线可得调谐值εσ＝－０ ０２和εσ＝０ ０６时的振幅稳定解Ｓ１㊁Ｓ２和Ｓ３㊂为了对解析结果进行数据仿真验证，选取相同参数和相应的调谐值，直接对式（１６）数值求解，得到调谐值εσ＝－０ ０２和εσ＝０ ０６时的稳定状态下的响应图［图８（ｂ）㊁图８（ｃ）］㊂经过对比可见，数值结果与解析结果基本一致㊂图８㊀幅频特性曲线及稳定解Ｓ１㊁Ｓ２和Ｓ３的响应图Ｆｉｇ．８㊀Ａｍｐｌｉｔｕｄｅ⁃ｆｒｅｑｕｅｎｃｙｃｈａｒａｃｔｅｒｉｓｔｉｃｃｕｒｖｅｓａｎｄｒｅｓｐｏｎｓｅｄｉａｇｒａｍｓｏｆｓｔａｂｌｅｓｏｌｕｔｉｏｎｓＳ１，Ｓ２ａｎｄＳ３㊀５２６㊀机㊀㊀械㊀㊀强㊀㊀度２０２３年㊀４㊀结论㊀㊀本文针对轴向运动铁磁板，考虑谐变磁化力和运动效应，得到非线性主共振的近似解析解，确定了共振幅值随几何和物理参量的变化规律㊂结果表明：１）当激励力频率接近条形板第一阶固有频率时，系统发生主共振，共振区域内的振幅明显增加，并且受板厚㊁磁场强度幅值和速度等参量的显著影响㊂２）铁磁板磁化产生的谐变磁力与机械载荷共同组成激励力，使磁场强度不单单影响阻尼项，同时出现于激励力项中，使系统呈现更加复杂的非线性振动特征㊂参考文献（Ｒｅｆｅｒｅｎｃｅｓ）［１］㊀ＬＩＮＣＣ．Ｓｔａｂｉｌｉｔｙａｎｄｖｉｂｒａｔｉｏｎｃｈａｒａｃｔｅｒｉｓｔｉｃｓｏｆａｘｉａｌｌｙｍｏｖｉｎｇｐｌａｔｅｓ［Ｊ］．ＩｎｔｅｒｎａｔｉｏｎａｌＪｏｕｒｎａｌｏｆＳｏｌｉｄｓａｎｄＳｔｒｕｃｔｕｒｅｓ，１９９７，３４（２４）：３１７９⁃３１９０．［２］㊀ＺＨＯＵＹＦ，ＷＡＮＧＺＭ．Ｖｉｂｒａｔｉｏｎｓｏｆａｘｉａｌｌｙｍｏｖｉｎｇｖｉｓｃｏｅｌａｓｔｉｃｐｌａｔｅｗｉｔｈｐａｒａｂｏｌｉｃａｌｌｙｖａｒｙｉｎｇｔｈｉｃｋｎｅｓｓ［Ｊ］．ＪｏｕｒｎａｌｏｆＳｏｕｎｄａｎｄＶｉｂｒａｔｉｏｎ，２００８，３１６（１／２／３／４／５）：１９８⁃２１０．［３］㊀ＭＯＨＡＭＡＤＩＡ，ＳＨＡＨＧＨＯＬＩＭ，ＡＳＨＥＮＡＩＧＦ．Ｆｒｅｅｖｉｂｒａｔｉｏｎａｎｄｓｔａｂｉｌｉｔｙｏｆａｎａｘｉａｌｌｙｍｏｖｉｎｇｔｈｉｎｃｉｒｃｕｌａｒｃｙｌｉｎｄｒｉｃａｌｓｈｅｌｌｕｓｉｎｇｍｕｌｔｉｐｌｅｓｃａｌｅｓｍｅｔｈｏｄ［Ｊ］．Ｍｅｃｃａｎｉｃａ，２０１９，５４（１４）：２２２７⁃２２４６．［４］㊀ＬＩＨＹ，ＬＡＮＧＴＹ，ＬＩＵＹＪ，ｅｔａｌ．Ｎｏｎｌｉｎｅａｒｖｉｂｒａｔｉｏｎｓａｎｄｓｔａｂｉｌｉｔｙｏｆａｎａｘｉａｌｌｙｍｏｖｉｎｇｐｌａｔｅｉｍｍｅｒｓｅｄｉｎｆｌｕｉｄ［Ｊ］．ＡｃｔａＭｅｃｈａｎｉｃａＳｏｌｉｄａＳｉｎｉｃａ，２０１９，３２（６）：７３７⁃７５３．［５］㊀ＦＡＲＳＨＢＡＦＺＲ，ＲＥＺＡＥＥＭ，ＬＯＴＦＡＮＳ．ＮｏｎｌｉｎｅａｒｖｉｂｒａｔｉｏｎａｎｄｓｔａｂｉｌｉｔｙａｎａｌｙｓｉｓｏｆｖｉｓｃｏｅｌａｓｔｉｃＲａｙｌｅｉｇｈｂｅａｍｓａｘｉａｌｌｙｍｏｖｉｎｇｏｎａｆｌｅｘｉｂｌｅｉｎｔｅｒｍｅｄｉａｔｅｓｕｐｐｏｒｔ［Ｊ］．ＩｒａｎｉａｎＪｏｕｒｎａｌｏｆＳｃｉｅｎｃｅａｎｄＴｅｃｈｎｏｌｏｇｙ，ＴｒａｎｓａｃｔｉｏｎｓｏｆＭｅｃｈａｎｉｃａｌＥｎｇｉｎｅｅｒｉｎｇ，２０２０，４４（４）：８６５⁃８７９．［６］㊀ＨＵＹＤ，ＨＵＰ，ＺＨＡＮＧＪＺ．Ｓｔｒｏｎｇｌｙｎｏｎｌｉｎｅａｒｓｕｂｈａｒｍｏｎｉｃｒｅｓｏｎａｎｃｅａｎｄｃｈａｏｔｉｃｍｏｔｉｏｎｏｆａｘｉａｌｌｙｍｏｖｉｎｇｔｈｉｎｐｌａｔｅｉｎｍａｇｎｅｔｉｃｆｉｅｌｄ［Ｊ］．ＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌａｎｄＮｏｎｌｉｎｅａｒＤｙｎａｍｉｃｓ，２０１５，１０（２）：０２１０１０（１⁃１２）．［７］㊀黄建亮，陈树辉．轴向运动体系非线性振动分析的多元Ｌ⁃Ｐ方法［Ｊ］．中山大学学报（自然科学版），２００４，４３（４）：１１５⁃１１７．ＨＵＡＮＧＪｉａｎＬｉａｎｇ，ＣＨＥＮＳｈｕＨｕｉ．ＭｕｌｔｉｖａｒｉａｔｅＬ⁃Ｐｍｅｔｈｏｄｆｏｒｎｏｎｌｉｎｅａｒｖｉｂｒａｔｉｏｎａｎａｌｙｓｉｓｏｆａｘｉａｌｌｙｍｏｖｉｎｇｓｙｓｔｅｍｓ［Ｊ］．ＡｃｔａＳｃｉｅｎｔｉａｒｕｍＮａｔｕｒａｌｉｕｍＵｎｉｖｅｒｓｉｔａｔｉｓＳｕｎｙａｔｓｅｎｉ，２００４，４３（４）：１１５⁃１１７（ＩｎＣｈｉｎｅｓｅ）．［８］㊀黄建亮，陈树辉．轴向运动体系横向非线性振动的联合共振［Ｊ］．振动工程学报，２００５，１９（１）：２４⁃２８．ＨＵＡＮＧＪｉａｎＬｉａｎｇ，ＣＨＥＮＳｈｕＨｕｉ．Ｊｏｉｎｔｒｅｓｏｎａｎｃｅｏｆｌａｔｅｒａｌｎｏｎｌｉｎｅａｒｖｉｂｒａｔｉｏｎｏｆａｘｉａｌｌｙｍｏｖｉｎｇｓｙｓｔｅｍｓ［Ｊ］．ＪｏｕｒｎａｌｏｆＶｉｂｒａｔｉｏｎＥｎｇｉｎｅｅｒｉｎｇ，２００５，１９（１）：２４⁃２８（ＩｎＣｈｉｎｅｓｅ）．［９］㊀殷振坤，陈树辉．轴向运动薄板非线性振动及其稳定性研究［Ｊ］．动力学与控制学报，２００７，５（４）：３１４⁃３１９．ＹＩＮＺｈｅｎＫｕｎ，ＣＨＥＮＳｈｕＨｕｉ．Ｎｏｎｌｉｎｅａｒｖｉｂｒａｔｉｏｎａｎｄｓｔａｂｉｌｉｔｙｏｆａｘｉａｌｍｏｖｉｎｇｔｈｉｎｐｌａｔｅ［Ｊ］．ＪｏｕｒｎａｌｏｆＤｙｎａｍｉｃｓａｎｄＣｏｎｔｒｏｌ，２００７，５（４）：３１４⁃３１９（ＩｎＣｈｉｎｅｓｅ）．［１０］㊀李中华，李映辉．轴向运动黏弹性夹层板的多模态耦合横向振动［Ｊ］．复合材料学报，２０１２，２９（３）：２１９⁃２２５．ＬＩＺｈｏｎｇＨｕａ，ＬＩＹｉｎｇＨｕｉ．Ｍｕｌｔｉｍｏｄａｌｃｏｕｐｌｅｄｔｒａｎｓｖｅｒｓｅｖｉｂｒａｔｉｏｎｏｆａｎａｘｉａｌｌｙｍｏｖｉｎｇｖｉｓｃｏｅｌａｓｔｉｃｓａｎｄｗｉｃｈｐｌａｔｅ［Ｊ］．ＡｃｔａＭａｔｅｒｉａｅＣｏｍｐｏｓｉｔａｅＳｉｎｉｃａ，２０１２，２９（３）：２１９⁃２２５（ＩｎＣｈｉｎｅｓｅ）．［１１］㊀李㊀晶，胡宇达．横向磁场中矩形导电薄板的内共振特性分析［Ｊ］．机械强度，２０１７，３９（６）：１２５５⁃１２６３．ＬＩＪｉｎｇ，ＨＵＹｕＤａ．Ａｎａｌｙｓｉｓｏｆｉｎｔｅｒｎａｌｒｅｓｏｎａｎｃｅｃｈａｒａｃｔｅｒｉｓｔｉｃｓｏｆｒｅｃｔａｎｇｕｌａｒｃｕｒｒｅｎｔ⁃ｃｏｎｄｕｃｔｉｎｇｔｈｉｎｐｌａｔｅｉｎｔｒａｎｓｖｅｒｓｅｍａｇｎｅｔｉｃｆｉｅｌｄ［Ｊ］．ＪｏｕｒｎａｌｏｆＭｅｃｈａｎｉｃａｌＳｔｒｅｎｇｔｈ，２０１７，３９（６）：１２５５⁃１２６３（ＩｎＣｈｉｎｅｓｅ）．［１２］㊀ＨＵＹＤ，ＷＡＮＧＪ．Ｐｒｉｎｃｉｐａｌ⁃ｉｎｔｅｒｎａｌｒｅｓｏｎａｎｃｅｏｆａｎａｘｉａｌｌｙｍｏｖｉｎｇｃｕｒｒｅｎｔ⁃ｃａｒｒｙｉｎｇｂｅａｍｉｎｍａｇｎｅｔｉｃｆｉｅｌｄ［Ｊ］．ＮｏｎｌｉｎｅａｒＤｙｎａｍｉｃｓ，２０１７，９０（１）：６８３⁃６９５．［１３］㊀高原文，周又和，郑晓静．横向磁场激励下铁磁梁式板的混沌运动分析［Ｊ］．力学学报，２００２，４６（１）：１０１⁃１０８．ＧＡＯＹｕａｎＷｅｎ，ＺＨＯＵＹｏｕＨｅ，ＺＨＥＮＧＸｉａｏＪｉｎｇ．Ｃｈａｏｔｉｃｍｏｔｉｏｎａｎａｌｙｓｉｓｏｆｆｅｒｒｏｍａｇｎｅｔｉｃｂｅａｍｐｌａｔｅｕｎｄｅｒｔｒａｎｓｖｅｒｓｅｍａｇｎｅｔｉｃｆｉｅｌｄｅｘｃｉｔａｔｉｏｎ［Ｊ］．ＣｈｉｎｅｓｅＪｏｕｒｎａｌｏｆＴｈｅｏｒｅｔｉｃａｌａｎｄＡｐｐｌｉｅｄＭｅｃｈａｎｉｃｓ，２００２，４６（１）：１０１⁃１０８（ＩｎＣｈｉｎｅｓｅ）．［１４］㊀高原文，缑新科，周又和．脉冲磁场激励下铁磁梁式板动力响应特征研究［Ｊ］．振动工程学报，２００５，１９（３）：３１４⁃３１７．ＧＡＯＹｕａｎＷｅｎ，ＧＯＵＸｉｎＫｅ，ＺＨＯＵＹｏｕＨｅ．Ｄｙｎａｍｉｃｒｅｓｐｏｎｓｅｃｈａｒａｃｔｅｒｉｓｔｉｃｓｏｆｆｅｒｒｏｍａｇｎｅｔｉｃｂｅａｍｐｌａｔｅｕｎｄｅｒｐｕｌｓｅｄｍａｇｎｅｔｉｃｆｉｅｌｄｅｘｃｉｔａｔｉｏｎ［Ｊ］．ＪｏｕｒｎａｌｏｆＶｉｂｒａｔｉｏｎＥｎｇｉｎｅｅｒｉｎｇ，２００５，１９（３）：３１４⁃３１７（ＩｎＣｈｉｎｅｓｅ）．［１５］㊀徐浩然，胡宇达，李文平．载流线圈中导电圆板的磁弹性固有振动［Ｊ］．机械强度，２０１９，４１（６）：１２７１⁃１２７７．ＸＵＨａｏＲａｎ，ＨＵＹｕＤａ，ＬＩＷｅｎＰｉｎｇ．Ｍａｇｎｅｔｏｅｌａｓｔｉｃｎａｔｕｒａｌｖｉｂｒａｔｉｏｎｏｆｃｏｎｄｕｃｔｉｖｅｃｉｒｃｕｌａｒｐｌａｔｅｉｎｃｕｒｒｅｎｔ⁃ｃａｒｒｙｉｎｇｃｏｉｌｓ［Ｊ］．ＪｏｕｒｎａｌｏｆＭｅｃｈａｎｉｃａｌＳｔｒｅｎｇｔｈ，２０１９，４１（６）：１２７１⁃１２７７（ＩｎＣｈｉｎｅｓｅ）．［１６］㊀ＨＡＳＡＮＹＡＮＤ，ＬＩＢＲＥＳＣＵＬ，ＱＩＮＺ，ｅｔａｌ．Ｎｏｎｌｉｎｅａｒｖｉｂｒａｔｉｏｎｏｆｆｉｎｉｔｅｌｙ⁃ｅｌｅｃｔｒｏｃｏｎｄｕｃｔｉｖｅｐｌａｔｅｓｔｒｉｐｓｉｎａｎａｘｉａｌｍａｇｎｅｔｉｃｆｉｅｌｄ［Ｊ］．Ｃｏｍｐｕｔｅｒｓ＆Ｓｔｒｕｃｔｕｒｅｓ，２００５，８３（１５／１６）：１２０５⁃１２１６．［１７］㊀ＷＡＮＧＸ，ＬＥＥＪＳ．Ｄｙｎａｍｉｃｓｔａｂｉｌｉｔｙｏｆｆｅｒｒｏｍａｇｎｅｔｉｃｂｅａｍ⁃ｐｌａｔｅｓｗｉｔｈｍａｇｎｅｔｏ⁃ｅｌａｓｔｉｃｉｎｔｅｒａｃｔｉｏｎａｎｄｍａｇｎｅｔｉｃｄａｍｐｉｎｇｉｎｔｒａｎｓｖｅｒｓｅｍａｇｎｅｔｉｃｆｉｅｌｄｓ［Ｊ］．ＪｏｕｒｎａｌｏｆＥｎｇｉｎｅｅｒｉｎｇＭｅｃｈａｎｉｃｓ，２００６，１３２（４）：４２２⁃４２８．［１８］㊀ＨＵＹＤ，ＬＩＪ．Ｔｈｅｍａｇｎｅｔｏ⁃ｅｌａｓｔｉｃｓｕｂｈａｒｍｏｎｉｃｒｅｓｏｎａｎｃｅｏｆｃｕｒｒｅｎｔ⁃ｃｏｎｄｕｃｔｉｎｇｔｈｉｎｐｌａｔｅｉｎｍａｇｎｅｔｉｃｆｉｌｅｄ［Ｊ］．ＪｏｕｒｎａｌｏｆＳｏｕｎｄａｎｄＶｉｂｒａｔｉｏｎ，２００９，３１９（３／４／５）：１１０７⁃１１２０．［１９］㊀周又和，郑晓静．电磁固体力学［Ｍ］．北京：科学出版社，１９９９：８２⁃８９．ＺＨＯＵＹｏｕＨｅ，ＺＨＥＮＧＸｉａｏＪｉｎｇ．Ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃｓｏｌｉｄｍｅｃｈａｎｉｃｓ［Ｍ］．Ｂｅｉｊｉｎｇ：ＳｃｉｅｎｃｅＰｒｅｓｓ，１９９９：８２⁃８９（ＩｎＣｈｉｎｅｓｅ）．［２０］㊀胡宇达．轴向运动导电薄板磁弹性耦合动力学理论模型［Ｊ］．固体力学学报，２０１３，３４（４）：４１７⁃４２５．ＨＵＹｕＤａ．Ｍａｇｎｅｔｏｅｌａｓｔｉｃｃｏｕｐｌｅｄｄｙｎａｍｉｃｓｔｈｅｏｒｅｔｉｃａｌｍｏｄｅｌｏｆａｘｉａｌｌｙｍｏｖｉｎｇｃｕｒｒｅｎｔ⁃ｃｏｎｄｕｃｔｉｎｇｔｈｉｎｐｌａｔｅｓ［Ｊ］．ＣｈｉｎｅｓｅＪｏｕｒｎａｌｏｆＳｏｌｉｄＭｅｃｈａｎｉｃｓ，２０１３，３４（４）：４１７⁃４２５（ＩｎＣｈｉｎｅｓｅ）．。

第９卷㊀第２期２０２４年３月气体物理ＰＨＹＳＩＣＳＯＦＧＡＳＥＳＶｏｌ.９㊀Ｎｏ.２Ｍａｒ.２０２４㊀㊀ＤＯＩ:１０.１９５２７/ｊ.ｃｎｋｉ.２０９６￣１６４２.１０９８Ｍａｃｈ数和壁面温度对ＨｙＴＲＶ边界层转捩的影响章录兴ꎬ㊀王光学ꎬ㊀杜㊀磊ꎬ㊀余发源ꎬ㊀张怀宝(中山大学航空航天学院ꎬ广东深圳５１８１０７)ＥｆｆｅｃｔｓｏｆＭａｃｈＮｕｍｂｅｒａｎｄＷａｌｌＴｅｍｐｅｒａｔｕｒｅｏｎＨｙＴＲＶＢｏｕｎｄａｒｙＬａｙｅｒＴｒａｎｓｉｔｉｏｎＺＨＡＮＧＬｕｘｉｎｇꎬ㊀ＷＡＮＧＧｕａｎｇｘｕｅꎬ㊀ＤＵＬｅｉꎬ㊀ＹＵＦａｙｕａｎꎬ㊀ＺＨＡＮＧＨｕａｉｂａｏ(ＳｃｈｏｏｌｏｆＡｅｒｏｎａｕｔｉｃｓａｎｄＡｓｔｒｏｎａｕｔｉｃｓꎬＳｕｎＹａｔ￣ｓｅｎＵｎｉｖｅｒｓｉｔｙꎬＳｈｅｎｚｈｅｎ５１８１０７ꎬＣｈｉｎａ)摘㊀要:典型的高超声速飞行器流场存在着复杂的转捩现象ꎬ其对飞行器的性能有着显著的影响ꎮ针对ＨｙＴＲＶ这款接近真实高超声速飞行器的升力体模型ꎬ采用数值模拟方法ꎬ研究Ｍａｃｈ数和壁面温度对ＨｙＴＲＶ转捩的影响规律ꎮ采用课题组自研软件开展数值计算ꎬＭａｃｈ数的范围为３~８ꎬ壁面温度的范围为１５０~９００Ｋꎮ首先对γ￣Ｒｅ~θｔ转捩模型和ＳＳＴ湍流模型进行了高超声速修正:将压力梯度系数修正㊁高速横流修正引入到γ￣Ｒｅ~θｔ转捩模型ꎬ并对ＳＳＴ湍流模型闭合系数β∗和β进行可压缩修正ꎻ然后开展了网格无关性验证ꎬ通过与实验结果对比ꎬ确认了修正后的数值方法和软件平台ꎻ最终开展Ｍａｃｈ数和壁面温度对ＨｙＴＲＶ边界层转捩规律的影响研究ꎮ计算结果表明ꎬ转捩区域主要集中在上表面两侧㊁下表面中心线两侧ꎻ增大来流Ｍａｃｈ数ꎬ上下表面转捩起始位置均大幅后移ꎬ湍流区大幅缩小ꎬ但仍会存在ꎬ同时上表面层流区摩阻系数不断增大ꎬ下表面湍流区摩阻系数不断减小ꎻ升高壁面温度ꎬ上下表面转捩起始位置先前移ꎬ然后快速后移ꎬ最终湍流区先后几乎消失ꎮ关键词:转捩ꎻＨｙＴＲＶꎻ摩阻ꎻＭａｃｈ数ꎻ壁面温度㊀㊀㊀收稿日期:２０２３￣１２￣１３ꎻ修回日期:２０２４￣０１￣０２基金项目:国家重大项目(ＧＪＸＭ９２５７９)ꎻ广东省自然科学基金－面上项目(２０２３Ａ１５１５０１００３６)ꎻ中山大学中央高校基本科研业务费专项资金(２２ｑｎｔｄ０７０５)第一作者简介:章录兴(１９９８ )㊀男ꎬ硕士ꎬ主要研究方向为高超声速空气动力学ꎮＥ￣ｍａｉｌ:184****8082＠１６３.ｃｏｍ通信作者简介:张怀宝(１９８５ )㊀男ꎬ副教授ꎬ主要研究方向为空气动力学ꎮＥ￣ｍａｉｌ:ｚｈａｎｇｈｂ２８＠ｍａｉｌ.ｓｙｓｕ.ｅｄｕ.ｃｎ中图分类号:Ｖ２１１ꎻＶ４１１㊀㊀文献标志码:ＡＡｂｓｔｒａｃｔ:Ｔｈｅｒｅｉｓａｃｏｍｐｌｅｘｔｒａｎｓｉｔｉｏｎｐｈｅｎｏｍｅｎｏｎｉｎｔｈｅｆｌｏｗｆｉｅｌｄｏｆａｔｙｐｉｃａｌｈｙｐｅｒｓｏｎｉｃｖｅｈｉｃｌｅꎬｗｈｉｃｈｈａｓａｓｉｇｎｉｆｉ￣ｃａｎｔｉｍｐａｃｔｏｎｔｈｅｐｅｒｆｏｒｍａｎｃｅｏｆｔｈｅｖｅｈｉｃｌｅ.ＴｈｅｅｆｆｅｃｔｓｏｆＭａｃｈｎｕｍｂｅｒａｎｄｗａｌｌｔｅｍｐｅｒａｔｕｒｅｏｎｔｈｅｔｒａｎｓｉｔｉｏｎｏｆＨｙＴＲＶｗｅｒｅｓｔｕｄｉｅｄｂｙｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎｍｅｔｈｏｄｓ.Ｔｈｅｓｅｌｆ￣ｄｅｖｅｌｏｐｅｄｓｏｆｔｗａｒｅｏｆｔｈｅｒｅｓｅａｒｃｈｇｒｏｕｐｗａｓｕｓｅｄｔｏｃａｒｒｙｏｕｔｎｕ￣ｍｅｒｉｃａｌｃａｌｃｕｌａｔｉｏｎｓ.ＴｈｅｒａｎｇｅｏｆＭａｃｈｎｕｍｂｅｒｗａｓ３~８ꎬａｎｄｔｈｅｒａｎｇｅｏｆｗａｌｌｔｅｍｐｅｒａｔｕｒｅｗａｓ１５０~９００Ｋ.Ｆｉｒｓｔｌｙꎬｔｈｅｈｙｐｅｒｓｏｎｉｃｃｏｒｒｅｃｔｉｏｎｓｏｆｔｈｅγ￣Ｒｅ~θｔｔｒａｎｓｉｔｉｏｎｍｏｄｅｌａｎｄｔｈｅＳＳＴｔｕｒｂｕｌｅｎｃｅｍｏｄｅｌｗｅｒｅｃａｒｒｉｅｄｏｕｔ.Ｔｈｅｐｒｅｓｓｕｒｅｇｒａｄｉｅｎｔｃｏｅｆｆｉｃｉｅｎｔｃｏｒｒｅｃｔｉｏｎａｎｄｔｈｅｈｉｇｈ￣ｓｐｅｅｄｃｒｏｓｓ￣ｆｌｏｗｃｏｒｒｅｃｔｉｏｎｗｅｒｅｉｎｔｒｏｄｕｃｅｄｉｎｔｏｔｈｅγ￣Ｒｅ~θｔｔｒａｎｓｉｔｉｏｎｍｏｄｅｌꎬａｎｄｔｈｅｃｏｍ￣ｐｒｅｓｓｉｂｉｌｉｔｙｃｏｒｒｅｃｔｉｏｎｓｏｆｔｈｅｃｌｏｓｕｒｅｃｏｅｆｆｉｃｉｅｎｔｓβ∗ａｎｄβｏｆｔｈｅＳＳＴｔｕｒｂｕｌｅｎｃｅｍｏｄｅｌｗｅｒｅｃａｒｒｉｅｄｏｕｔ.Ｔｈｅｎꎬｔｈｅｇｒｉｄｉｎ￣ｄｅｐｅｎｄｅｎｃｅｖｅｒｉｆｉｃａｔｉｏｎｗａｓｃａｒｒｉｅｄｏｕｔꎬａｎｄｔｈｅｍｏｄｉｆｉｅｄｎｕｍｅｒｉｃａｌｍｅｔｈｏｄａｎｄｓｏｆｔｗａｒｅｐｌａｔｆｏｒｍｗｅｒｅｃｏｎｆｉｒｍｅｄｂｙｃｏｍ￣ｐａｒｉｎｇｗｉｔｈｅｘｐｅｒｉｍｅｎｔａｌｒｅｓｕｌｔｓ.ＦｉｎａｌｌｙꎬｔｈｅｅｆｆｅｃｔｓｏｆＭａｃｈｎｕｍｂｅｒａｎｄｗａｌｌｔｅｍｐｅｒａｔｕｒｅｏｎｔｈｅｔｒａｎｓｉｔｉｏｎｌａｗｏｆｔｈｅＨｙＴＲＶｂｏｕｎｄａｒｙｌａｙｅｒｗｅｒｅｓｔｕｄｉｅｄ.Ｔｈｅｒｅｓｕｌｔｓｓｈｏｗｔｈａｔｔｈｅｔｒａｎｓｉｔｉｏｎａｒｅａｉｓｍａｉｎｌｙｃｏｎｃｅｎｔｒａｔｅｄｏｎｂｏｔｈｓｉｄｅｓｏｆｔｈｅｕｐｐｅｒｓｕｒｆａｃｅａｎｄｔｈｅｃｅｎｔｅｒｌｉｎｅｏｆｔｈｅｌｏｗｅｒｓｕｒｆａｃｅ.ＷｉｔｈｔｈｅｉｎｃｒｅａｓｅｏｆｔｈｅｉｎｃｏｍｉｎｇＭａｃｈｎｕｍｂｅｒꎬｔｈｅｓｔａｒｔｉｎｇｐｏｓｉｔｉｏｎｏｆｔｒａｎｓｉｔｉｏｎｏｎｔｈｅｕｐｐｅｒａｎｄｌｏｗｅｒｓｕｒｆａｃｅｓｉｓｇｒｅａｔｌｙｂａｃｋｗａｒｄꎬａｎｄｔｈｅｔｕｒｂｕｌｅｎｔｚｏｎｅｉｓｇｒｅａｔｌｙｒｅｄｕｃｅｄꎬｂｕｔｉｔｓｔｉｌｌｅｘ￣ｉｓｔｓ.Ａｔｔｈｅｓａｍｅｔｉｍｅꎬｔｈｅｆｒｉｃｔｉｏｎｃｏｅｆｆｉｃｉｅｎｔｏｆｔｈｅｌａｍｉｎａｒｆｌｏｗｚｏｎｅｏｎｔｈｅｕｐｐｅｒｓｕｒｆａｃｅｉｎｃｒｅａｓｅｓｃｏｎｔｉｎｕｏｕｓｌｙꎬａｎｄｔｈｅｆｒｉｃｔｉｏｎｃｏｅｆｆｉｃｉｅｎｔｏｆｔｈｅｔｕｒｂｕｌｅｎｔｚｏｎｅｏｎｔｈｅｌｏｗｅｒｓｕｒｆａｃｅｄｅｃｒｅａｓｅｓ.Ａｓｔｈｅｗａｌｌｔｅｍｐｅｒａｔｕｒｅｉｎｃｒｅａｓｅｓꎬｔｈｅｓｔａｒｔｉｎｇｐｏｓｉ￣ｔｉｏｎｏｆｔｒａｎｓｉｔｉｏｎｏｎｔｈｅｕｐｐｅｒａｎｄｌｏｗｅｒｓｕｒｆａｃｅｓｓｈｉｆｔｓｆｏｒｗａｒｄꎬｔｈｅｎｒａｐｉｄｌｙｓｈｉｆｔｓｂａｃｋｗａｒｄꎬａｎｄｆｉｎａｌｌｙｔｈｅｔｕｒｂｕｌｅｎｔｚｏｎｅａｌｍｏｓｔｄｉｓａｐｐｅａｒｓ.气体物理２０２４年㊀第９卷Ｋｅｙｗｏｒｄｓ:ｔｒａｎｓｉｔｉｏｎꎻＨｙＴＲＶꎻｆｒｉｃｔｉｏｎꎻＭａｃｈｎｕｍｂｅｒꎻｗａｌｌｔｅｍｐｅｒａｔｕｒｅ引㊀言高超声速飞行器具有突防能力强㊁打击范围广㊁响应迅速等显著优势ꎬ正逐渐成为各国空天竞争的热点[１]ꎮ高超声速飞行器边界层转捩是该类飞行器气动设计中的重要问题[２]ꎮ在边界层转捩过程中ꎬ流态由层流转变为湍流ꎬ飞行器的表面摩阻急剧增大到层流时的３~５倍ꎬ严重影响飞行器的气动性能与热防护系统ꎬ转捩还会导致飞行器壁面烧蚀㊁颤振加剧㊁飞行姿态控制难度大等一系列问题ꎬ对飞行器的飞行安全构成严重的威胁[３￣５]ꎬ开展高超声速飞行器边界层转捩研究具有十分重要的意义ꎮ影响边界层转捩的因素很多ꎬ例如ꎬＭａｃｈ数㊁Ｒｅｙｎｏｌｄｓ数㊁湍流强度㊁表面传导热等ꎮ在高超声速流动条件下ꎬ强激波㊁强逆压梯度㊁熵层等高超声速现象及其相互作用ꎬ会使得转捩流动的预测和研究难度进一步增大[６]ꎮ目前高超声速飞行器转捩数值模拟方法主要有直接数值模拟(ＤＮＳ)㊁大涡模拟(ＬＥＳ)和基于Ｒｅｙｎｏｌｄｓ平均Ｎａｖｉｅｒ￣Ｓｔｏｋｅｓ(ＲＡＮＳ)的转捩模型方法ꎬ由于前两种计算量巨大ꎬ难以推广到工程应用ꎬ基于Ｒｅｙｎｏｌｄｓ平均Ｎａｖｉｅｒ￣Ｓｔｏｋｅｓ的转捩模型在工程实践中应用最为广泛ꎬ其中γ￣Ｒｅ~θｔ转捩模型基于局部变量ꎬ与现代ＣＦＤ方法良好兼容ꎬ目前已经有多项研究尝试从一般性的流动问题拓展到高超声速流动转捩模拟[６￣９]ꎮ目前高超声速流动转捩的研究对象主要是结构相对简单的构型ꎮＭｃＤａｎｉｅｌ等[１０]研究了扩口直锥在高超声速流动条件下的转捩现象ꎮＰａｐｐ等[１１]研究了圆锥在高超声速流动条件下的转捩特性ꎮ美国和澳大利亚组织联合实施的ＨＩＦｉＲＥ计划[１２]ꎬ研究了圆锥形状的ＨＩＦｉＲＥ１和椭圆锥形的ＨＩＦｉＲＥ５的转捩问题ꎮ杨云军等[１３]采用数值模拟方法ꎬ分析了椭圆锥的转捩影响机制ꎬ并研究了Ｒｅｙｎｏｌｄｓ数对转捩特性的影响规律ꎮ另外ꎬ袁先旭等[１４]于２０１５年成功实施了圆锥体ＭＦ￣１航天模型飞行试验ꎮ以上对高超声速流动的转捩研究ꎬ都取得了比较理想的结果ꎬ然而所采用的模型都是圆锥㊁椭圆锥等简单几何外形ꎬ这与真实高超声速飞行器有较大差异ꎬ较难反映真实的转捩特性ꎮ为了有效促进对真实高超声速飞行器的转捩问题研究ꎬ中国空气动力研究与发展中心提出并设计了一款接近真实飞行器的升力体模型ꎬ即高超声速转捩研究飞行器(ｈｙｐｅｒｓｏｎｉｃｔｒａｎｓｉｔｉｏｎｒｅｓｅａｒｃｈｖｅｈｉｃｌｅꎬＨｙＴＲＶ)[１５]ꎬ模型详细的参数见参考文献[１６]ꎮＨｙＴＲＶ外形如图１所示ꎬ其整体外形较为复杂ꎬ不同区域发生转捩的情况也不尽相同ꎮ对ＨｙＴＲＶ的转捩问题研究能够显著提高对真实高超声速飞行器转捩特性的认识水平ꎮＬｉｕ等[１７]采用理论分析㊁数值模拟和风洞实验３种方法对ＨｙＴＲＶ的转捩特性进行了研究ꎻ陈坚强等[１５]分析了ＨｙＴＲＶ的边界层失稳特征ꎻＣｈｅｎ等[１８]对ＨｙＴＲＶ进行了多维线性稳定性分析ꎻＱｉ等[１９]在来流Ｍａｃｈ数６㊁攻角０ʎ的条件下对ＨｙＴＲＶ进行了直接数值模拟ꎻ万兵兵等[２０]结合风洞实验与飞行试验ꎬ利用ｅＮ方法预测了ＨｙＴＲＶ升力体横流区的转捩阵面形状ꎮ目前ꎬ相关研究主要集中在ＨｙＴＲＶ的稳定性特征及转捩预测两个方面ꎬ而对若干关键参数ꎬ特别是Ｍａｃｈ数和壁面温度对转捩的影响研究还比较少ꎮ(ａ)Ｆｒｏｎｔｖｉｅｗ(ｂ)Ｓｉｄｅｖｉｅｗ㊀㊀㊀图１㊀ＨｙＴＲＶ外形Ｆｉｇ.１㊀ＳｈａｐｅｏｆＨｙＴＲＶ基于此ꎬ本文采用数值模拟方法ꎬ应用课题组自研软件开展Ｍａｃｈ数和壁面温度对ＨｙＴＲＶ转捩流动的影响规律研究ꎮ１㊀数值方法１.１㊀控制方程和数值方法控制方程为三维可压缩ＲＡＮＳ方程ꎬ采用结构网格技术和有限体积方法ꎬ变量插值方法采用２阶ＭＵＳＣＬ格式ꎬ通量计算采用低耗散的通量向量差分Ｒｏｅ格式ꎬ黏性项离散采用中心格式ꎬ时间推进方法采用ＬＵ￣ＳＧＳ格式ꎮ壁面采用等温㊁无滑移壁面条件ꎬ入口采用Ｒｉｅｍａｎｎ远场边界条件ꎬ出口采用零梯度外推边界条件ꎮ１.２㊀γ￣Ｒｅ~θｔ转捩模型γ￣Ｒｅ~θｔ转捩模型是Ｍｅｎｔｅｒ等[２１ꎬ２２]于２００４年提０１第２期章录兴ꎬ等:Ｍａｃｈ数和壁面温度对ＨｙＴＲＶ边界层转捩的影响出的一种基于拟合公式的间歇因子转捩模型ꎬ在２００９年公布了完整的拟合公式及相关参数[２３]ꎮ许多学者也开发了相应的程序ꎬ并进行了大量的算例验证[２４￣２８]ꎬ证明了该模型具有较好的转捩预测能力ꎬ预测精度较高ꎻ通过合适的标定ꎬγ￣Ｒｅ~θｔ转捩模型可以适用于多种情况下的转捩模拟ꎮ该模型构建了关于间歇因子γ的输运方程和关于转捩动量厚度Ｒｅｙｎｏｌｄｓ数Ｒｅ~θｔ的输运方程ꎮ具体来说ꎬγ表示该位置是湍流流动的概率ꎬ取值范围为０<γ<１ꎮ关于γ的控制方程为Ə(ργ)Əｔ＋Ə(ρｕｊγ)Əｘｊ＝Ｐγ－Ｅγ＋ƏƏｘｊμ＋μｔσｆæèçöø÷ƏγƏｘｊéëêêùûúú其中ꎬＰγ为生成项ꎬＥγ为破坏项ꎮ关于Ｒｅ~θｔ的输运方程为Ə(ρＲｅ~θｔ)Əｔ＋Ə(ρｕｊＲｅ~θｔ)Əｘｊ＝Ｐθｔ＋ƏƏｘｊσθｔ(μ＋μｔ)ƏＲｅ~θｔƏｘｊéëêêùûúú其中ꎬＰθｔ为源项ꎬ其作用是使边界层外部的Ｒｅ~θｔ等于Ｒｅθｔꎬ定义式为Ｐθｔ＝ｃθｔρｔ(Ｒｅθｔ－Ｒｅ~θｔ)(１.０－Ｆθｔ)Ｒｅθｔ采用以下经验公式Ｒｅθｔ＝１１７３.５１－５８９ ４２８Ｔｕ＋０.２１９６Ｔｕ２æèçöø÷Ｆ(λθ)ꎬＴｕɤ０.３Ｒｅθｔ＝３３１.５０(Ｔｕ－０.５６５８)－０.６７１Ｆ(λθ)ꎬＴｕ>０.３ìîíïïïïＦ(λθ)＝１＋(１２.９８６λθ＋１２３.６６λ２θ＋４０５.６８９λ３θ)ｅ－(Ｔｕ１.５)１.５ꎬ㊀λθɤ０Ｆ(λθ)＝１＋０.２７５(１－ｅ－３５.０λθ)ｅ－(Ｔｕ０.５)ꎬλθ>０ìîíïïïï在实际计算中ꎬ通过γ￣Ｒｅ~θｔ转捩模型获得间歇因子ꎬ再通过间歇因子来控制ＳＳＴｋ￣ω湍流模型中湍动能的生成ꎮγ￣Ｒｅ~θｔ转捩模型与ＳＳＴｋ￣ω湍流模型耦合为Ə(ρｋ)Əｔ＋Ə(ρｕｊｋ)Əｘｊ＝γｅｆｆτｉｊƏｕｉƏｘｊ－ｍｉｎ(ｍａｘ(γｅｆｆꎬ０.１)ꎬ１.０)ρβ∗ｋω＋ƏƏｘｊμ＋μｔσｋæèçöø÷ƏｋƏｘｊéëêêùûúúƏ(ρω)Əｔ＋Ə(ρｕｊω)Əｘｊ＝γｖｔτｉｊƏｕｉƏｘｊ－βρω２＋ƏƏｘｊ(μ＋σωμｔ)ƏωƏｘｊéëêêùûúú＋２ρ(１－Ｆ１)σω２１ωƏｋƏｘｊƏωƏｘｊ模型中具体参数定义见文献[２３]ꎮ１.３㊀高超声速修正原始ＳＳＴ湍流模型及γ￣Ｒｅ~θｔ转捩模型都是基于不可压缩流动发展的ꎬ为了更好地预测高超声速流动转捩ꎬ本节引入了３种重要的高超声速修正方法ꎮ１.３.１㊀压力梯度修正压力梯度对边界层转捩的影响较大ꎬ在高Ｍａｃｈ数情况下ꎬ边界层厚度较大ꎬ进而影响压力梯度的大小ꎬ因此在模拟高超声速流动时应该考虑Ｍａｃｈ数对压力梯度的影响ꎮ本文采用张毅峰等[２９]提出的压力梯度修正方法ꎬ具体修正形式如下λᶄθ＝λθ１＋γᶄ－１２Ｍａｅæèçöø÷其中ꎬＭａｅ为边界层外缘Ｍａｃｈ数ꎬγᶄ为比热比ꎮ１.３.２㊀高速横流修正在原始γ￣Ｒｅ~θｔ转捩模型中ꎬ没有考虑横流不稳定性对转捩的影响ꎬ对于横流模态主导的转捩ꎬ原始转捩模型计算的结果并不理想ꎮＬａｎｇｔｒｙ等[３０]在２０１５年对γ￣Ｒｅ~θｔ转捩模型进行了低速横流修正ꎬ向星皓等[９]在Ｌａｎｇｔｒｙ低速横流修正的基础上ꎬ对高超声速椭圆锥转捩ＤＮＳ数据进行了拓展ꎬ提出了高速横流转捩判据ꎬ本文直接采用向星皓提出的高速横流转捩方法ꎮＬａｎｇｔｒｙ将横流强度引入转捩发生动量厚度Ｒｅｙｎｏｌｄｓ数输运方程中Ə(ρＲｅ~θｔ)Əｔ＋Ə(ρｕｊＲｅ~θｔ)Əｘｊ＝Ｐθｔ＋ＤＳＣＦ＋ƏƏｘｊσθｔ(μ＋μｔ)ƏＲｅ~θｔƏｘｊéëêêùûúú式中ꎬＤＳＣＦ为横流源项ꎬＬａｎｇｔｒｙ低速横流修正为ＤＳＣＦ＝ｃθｔρｔｃｃｒｏｓｓｆｌｏｗｍｉｎ(ＲｅＳＣＦ－Ｒｅ~θｔꎬ０.０)Ｆθｔ２其中ꎬＲｅＳＣＦ为低速横流判据ＲｅＳＣＦ＝θｔρＵｌｏｃａｌ０.８２æèçöø÷μ＝－３５.０８８ｌｎｈθｔæèçöø÷＋３１９.５１＋ｆ(＋ΔＨｃｒｏｓｓｆｌｏｗ)－ｆ(－ΔＨｃｒｏｓｓｆｌｏｗ)其中ꎬｈ为壁面粗糙度高度ꎬθｔ为动量厚度ꎬ１１气体物理２０２４年㊀第９卷ΔＨｃｒｏｓｓｆｌｏｗ是横流强度抬升项ꎮ向星皓提出的高速横流转捩判据ꎬ其中高速横流源项ＤＳＣＦ￣Ｈ为ＤＳＣＦ￣Ｈ＝ｃＣＦρｍｉｎ(ＲｅＳＣＦ￣Ｈ－Ｒｅ~θｔꎬ０)ＦθｔＲｅＳＣＦ￣Ｈ＝ＣＣＦ￣１ｌｎｈｌμ＋ＣＣＦ￣２＋(Ｈｃｒｏｓｓｆｌｏｗ)其中ꎬＣＣＦ￣１＝－９.６１８ꎬＣＣＦ￣２＝１２８.３３ꎻｌμ为粗糙度参考高度ꎬｌμ＝１μｍꎻｆ(Ｈｃｒｏｓｓｆｌｏｗ)为抬升函数ｆ(Ｈｃｒｏｓｓｆｌｏｗ)＝６００００.１０６６－ΔＨｃｒｏｓｓｆｌｏｗ＋５００００(０.１０６６－ΔＨｃｒｏｓｓｆｌｏｗ)２其中ꎬΔＨｃｒｏｓｓｆｌｏｗ与Ｌａｎｇｔｒｙ低速横流修正中保持一致ꎮ１.３.３㊀ＳＳＴ可压缩修正高超声速流动具有强可压缩性ꎬ所以在进行高超声速计算时ꎬ应该对湍流模型进行可压缩修正ꎮＳａｒｋａｒ[３１]提出了膨胀耗散修正ꎬ对ＳＳＴ湍流模型中的闭合系数β∗ꎬβ进行了可压缩修正ꎬＷｉｌｃｏｘ[３２]在Ｓａｒｋａｒ修正的基础上考虑了可压缩生成项产生时的延迟效应ꎬ使得可压缩修正在湍流Ｍａｃｈ数较小的近壁面关闭ꎬ在湍流Ｍａｃｈ数较大的自由剪切层打开ꎬ本文采用Ｗｉｌｃｏｘ提出的可压缩性修正β∗＝β∗０[１＋ξ∗Ｆ(Ｍａｔ)]β＝β０－β∗０ξ∗Ｆ(Ｍａｔ)其中ꎬβ０ꎬβ∗均为原始模型中的系数ꎬξ∗＝１.５ꎮＦ(Ｍａｔ)＝[Ｍａｔ－Ｍａｔ０]Ｈ(Ｍａｔ－Ｍａｔ０)Ｍａｔ０＝１/４ꎬＨ(ｘ)＝０ꎬｘɤ０１ꎬｘ>０{其中ꎬＭａｔ＝２ｋ/ａ为湍流Ｍａｃｈ数ꎬａ为当地声速ꎮ２㊀网格无关性验证及数值方法确认２.１㊀网格无关性验证计算采用３套网格ꎬ考虑到ＨｙＴＲＶ的几何对称性ꎬ生成３套半模网格ꎬ第１层网格高度为１ˑ１０－６ｍꎬ确保ｙ＋<１ꎬ流向ˑ法向ˑ周向的网格数分别为:网格１是３０１ˑ２０１ˑ２０１ꎬ网格２是３０１ˑ３０１ˑ２０１ꎬ网格３是４０１ˑ３８１ˑ２８１ꎮ全模下表面如图２所示ꎬ选取ｙ/Ｌ＝０中心线和ｘ/Ｌ＝０.５处ꎬ对比３套网格的表面摩阻系数ꎬ计算结果如图３所示ꎮ采用网格１时ꎬ表面摩阻系数分布与另外两个结果存在明显差异ꎻ而采用网格２和网格３时ꎬ表面摩阻系数曲线基本重合ꎬ表明在流向㊁法向和周向均满足网格无关性ꎬ后续数值计算采用网格２ꎮ图２㊀截取位置示意图Ｆｉｇ.２㊀Ｓｃｈｅｍａｔｉｃｄｉａｇｒａｍｏｆｔｈｅｉｎｔｅｒｃｅｐｔｉｏｎｌｏｃａｔｉｏｎ(ａ)Ｓｕｒｆａｃｅｆｒｉｃｔｉｏｎａｔｙ/Ｌ＝０(ｂ)Ｓｕｒｆａｃｅｆｒｉｃｔｉｏｎａｔｘ/Ｌ＝０.５图３㊀采用３套网格计算得到的摩阻对比Ｆｉｇ.３㊀Ｃｏｍｐａｒｉｓｏｎｏｆｔｈｅｆｒｉｃｔｉｏｎｄｒａｇｃａｌｃｕｌａｔｅｄｕｓｉｎｇｔｈｒｅｅｓｅｔｓｏｆｇｒｉｄｓ２.２㊀数值方法和自研软件的确认采用修正后的转捩模型对ＨｙＴＲＶ开展计算ꎬ计算工况为Ｍａ＝６ꎬ来流温度Ｔɕ＝９７Ｋꎬ单位２１第２期章录兴ꎬ等:Ｍａｃｈ数和壁面温度对ＨｙＴＲＶ边界层转捩的影响Ｒｅｙｎｏｌｄｓ数为Ｒｅ＝１.１ˑ１０７/ｍꎬ攻角α＝０ʎꎬ来流湍流度ＦＳＴＩ＝０.８％ꎬ壁面温度Ｔ＝３００Ｋꎮ为方便对比分析ꎬ计算结果与参考结果均采用上下对称形式布置ꎬ例如ꎬ图４是模型下表面计算结果与实验结果对比:对于下表面两侧转捩的起始位置ꎬ高超声速修正前的转捩位置在ｘ＝０.６８ｍ附近ꎬ高超声速修正后的计算结果与实验结果吻合良好ꎬ均在ｘ＝０.６０ｍ附近ꎬ并且湍流边界层区域形状基本一致ꎬ说明修正后的转捩模型能够较好地预测ＨｙＴＲＶ转捩的位置ꎮ(ａ)Ｃａｌｃｕｌａｔｉｏｎｏｆｔｈｅｆｒｉｃｔｉｏｎｄｉｓｔｒｉｂｕｔｉｏｎ(ｂｅｆｏｒｅｈｙｐｅｒｓｏｎｉｃｃｏｒｒｅｃｔｉｏｎ)(ｂ)Ｃａｌｃｕｌａｔｉｏｎｏｆｔｈｅｆｒｉｃｔｉｏｎｄｉｓｔｒｉｂｕｔｉｏｎ(ａｆｔｅｒｈｙｐｅｒｓｏｎｉｃｃｏｒｒｅｃｔｉｏｎ)(ｃ)Ｅｘｐｅｒｉｍｅｎｔａｌｒｅｓｕｌｔｓｏｆｔｈｅｈｅａｔｆｌｕｘｄｉｓｔｒｉｂｕｔｉｏｎ[１７]图４㊀下表面计算结果和实验结果对比Ｆｉｇ.４㊀Ｃｏｍｐａｒｉｓｏｎｏｆｔｈｅｃａｌｃｕｌａｔｅｄａｎｄｅｘｐｅｒｉｍｅｎｔａｌｒｅｓｕｌｔｓｏｎｔｈｅｌｏｗｅｒｓｕｒｆａｃｅ３㊀ＨｙＴＲＶ转捩的基本流动特性计算工况采用Ｍａ＝６ꎬ攻角α＝０ʎꎬ来流湍流度ＦＳＴＩ＝０.６％ꎬ分析ＨｙＴＲＶ转捩的基本流动特性ꎮ从图５可以看出ꎬ模型两侧和顶端均出现高压区ꎬ高压区之间为低压区ꎬ横截面上存在周向压力梯度ꎬ流动从高压区向低压区汇集ꎬ从而在下表面中心线附近和上表面两侧腰部区域均形成流向涡结构(见图６)ꎬ沿流动方向ꎬ高压区域逐渐扩大ꎬ流向涡结构的影响范围也越大ꎮ在流向涡结构的边缘位置ꎬ壁面附近的低速流体被抬升到外壁面区域ꎬ外壁面区域的高速流体又被带入到近壁面区域ꎬ进而导致流向涡结构边缘处壁面的摩阻显著增加ꎬ最终诱发转捩ꎬ这些流动特征与文献[１５]的结果一致ꎮ图７显示了上下表面摩阻的分布情况ꎬ其中上表面两侧区域在ｘ/Ｌ＝０.８０附近ꎬ摩阻显著增加ꎬ出现明显的转捩现象ꎬ转捩区域分布在两侧边缘位置ꎻ而下表面两侧区域在ｘ/Ｌ＝０.７５附近ꎬ也出现明显的转捩ꎬ转捩区域相对集中在中心线两侧ꎮ图５㊀不同截面位置处的压力云图Ｆｉｇ.５㊀Ｐｒｅｓｓｕｒｅｃｏｎｔｏｕｒｓａｔｄｉｆｆｅｒｅｎｔｃｒｏｓｓ￣ｓｅｃｔｉｏｎｌｏｃａｔｉｏｎｓ图６㊀不同截面位置处的流向速度云图Ｆｉｇ.６㊀Ｓｔｒｅａｍｗｉｓｅｖｅｌｏｃｉｔｙｃｏｎｔｏｕｒｓａｔｄｉｆｆｅｒｅｎｔｃｒｏｓｓ￣ｓｅｃｔｉｏｎｌｏｃａｔｉｏｎｓ３１气体物理２０２４年㊀第９卷(ａ)Ｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀(ｂ)Ｌｏｗｅｒｓｕｒｆａｃｅ图７㊀上下表面摩阻分布云图Ｆｉｇ.７㊀Ｆｒｉｃｔｉｏｎｃｏｅｆｆｉｃｉｅｎｔｃｏｎｔｏｕｒｓｏｎｔｈｅｕｐｐｅｒａｎｄｌｏｗｅｒｓｕｒｆａｃｅｓ４㊀不同Ｍａｃｈ数对ＨｙＴＲＶ转捩的影响保持来流湍流度ＦＳＴＩ＝０.６％不变ꎬＭａｃｈ数变化范围为３~８ꎮ图８是不同Ｍａｃｈ数条件下ＨｙＴＲＶ上下表面的摩阻分布云图ꎬ从图中可知ꎬ随着Ｍａｃｈ数的增加ꎬ上下表面的湍流区域均逐渐减少ꎬ其中上表面两侧转捩起始位置由ｘ/Ｌ＝０.５６附近后移至ｘ/Ｌ＝０.９２附近ꎬ下表面两侧转捩起始位置由ｘ/Ｌ＝０.４８附近后移至ｘ/Ｌ＝０.９９附近ꎬ上下表面两侧转捩起始位置均大幅后移ꎬ说明Ｍａｃｈ数对ＨｙＴＲＶ转捩的影响很大ꎮｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀ｌｏｗｅｒｓｕｒｆａｃｅ(ａ)Ｍａ＝３ｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀ｌｏｗｅｒｓｕｒｆａｃｅ(ｂ)Ｍａ＝４ｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀ｌｏｗｅｒｓｕｒｆａｃｅ(ｃ)Ｍａ＝５４１第２期章录兴ꎬ等:Ｍａｃｈ数和壁面温度对ＨｙＴＲＶ边界层转捩的影响ｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀ｌｏｗｅｒｓｕｒｆａｃｅ(ｄ)Ｍａ＝６ｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀ｌｏｗｅｒｓｕｒｆａｃｅ(ｅ)Ｍａ＝７ｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀ｌｏｗｅｒｓｕｒｆａｃｅ(ｆ)Ｍａ＝８图８㊀不同Ｍａｃｈ数条件下摩阻系数分布云图Ｆｉｇ.８㊀ＦｒｉｃｔｉｏｎｃｏｅｆｆｉｃｉｅｎｔｃｏｎｔｏｕｒｓａｔｄｉｆｆｅｒｅｎｔＭａｃｈｎｕｍｂｅｒｓ上表面选取图７中ｚ/Ｌ＝０.１２的位置ꎬ下表面选取ｚ/Ｌ＝０.１０的位置进行分析ꎮ从图９中可以分析出ꎬ随着Ｍａｃｈ数的增加ꎬ上表面转捩起始位置不断后移ꎬ当Ｍａｃｈ数增加到７时ꎬ由于湍流区的缩小ꎬ此处位置不再发生转捩ꎬ此外ꎬＭａｃｈ数越高层流区摩阻系数越大ꎻ下表面转捩起始位置也不断后移ꎬ当Ｍａｃｈ数增加到８时ꎬ此处位置不再发生转捩ꎬ此外ꎬＭａｃｈ数越高ꎬ湍流区的摩阻系数越小ꎬ这些结论与关于来流Ｍａｃｈ数对转捩位置影响的普遍研究结论一致ꎮ(ａ)Ｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀(ｂ)Ｌｏｗｅｒｓｕｒｆａｃｅ图９㊀不同位置摩阻系数随Ｍａｃｈ数的变化Ｆｉｇ.９㊀ＶａｒｉａｔｉｏｎｏｆｆｒｉｃｔｉｏｎｃｏｅｆｆｉｃｉｅｎｔｗｉｔｈＭａｃｈｎｕｍｂｅｒａｔｄｉｆｆｅｒｅｎｔｌｏｃａｔｉｏｎｓ５１气体物理２０２４年㊀第９卷５㊀不同壁面温度对ＨｙＴＲＶ转捩的影响保持来流湍流度ＦＳＴＩ＝０.６％及Ｍａ＝６不变ꎬ壁面温度的变化范围为１５０~９００Ｋꎮ图１０是不同壁面温度条件下ＨｙＴＲＶ上下表面的摩阻分布云图ꎬ可以看出随着壁面温度的增加ꎬ上表面两侧湍流区域先是缓慢扩大ꎬ在壁面温度为５００Ｋ时湍流区域快速缩小ꎬ增加到９００Ｋ时ꎬ已无明显湍流区域ꎻ下表面两侧湍流区域先是无明显变化ꎬ同样当壁面温度升高到５００Ｋ时ꎬ湍流区域快速缩小ꎬ当壁面温度升高到７００Ｋ时ꎬ两侧已经无明显的湍流区域ꎬ相比上表面两侧湍流区域ꎬ下表面湍流区域消失得更早ꎮ由此可以得出壁面温度对转捩的产生有较大的影响ꎬ壁面温度增加到一定程度将导致ＨｙＴＲＶ没有明显的转捩现象ꎮｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀ｌｏｗｅｒｓｕｒｆａｃｅ(ａ)Ｔ＝１５０Ｋｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀ｌｏｗｅｒｓｕｒｆａｃｅ(ｂ)Ｔ＝２００Ｋｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀ｌｏｗｅｒｓｕｒｆａｃｅ(ｃ)Ｔ＝３００Ｋｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀ｌｏｗｅｒｓｕｒｆａｃｅ(ｄ)Ｔ＝５００Ｋ６１第２期章录兴ꎬ等:Ｍａｃｈ数和壁面温度对ＨｙＴＲＶ边界层转捩的影响ｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀ｌｏｗｅｒｓｕｒｆａｃｅ(ｅ)Ｔ＝７００Ｋｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀ｌｏｗｅｒｓｕｒｆａｃｅ(ｆ)Ｔ＝９００Ｋ图１０㊀不同壁面温度条件下摩阻系数分布云图Ｆｉｇ.１０㊀Ｆｒｉｃｔｉｏｎｃｏｅｆｆｉｃｉｅｎｔｃｏｎｔｏｕｒｓａｔｄｉｆｆｅｒｅｎｔｗａｌｌｔｅｍｐｅｒａｔｕｒｅｃｏｎｄｉｔｉｏｎｓ上表面选取ｚ/Ｌ＝０.１２５的位置ꎬ下表面选取ｚ/Ｌ＝０.１００的位置进行分析ꎮ从图１１中可以分析出ꎬ随着壁面温度的增加ꎬ上表面转捩起始位置先前移ꎬ当壁面温度增加到５００Ｋ时ꎬ转捩起始位置后移ꎬ转捩区长度逐渐增加ꎬ层流区域的摩阻系数逐渐增加ꎬ当壁面温度增加到７００Ｋ时ꎬ该位置已不再出现转捩ꎻ下表面转捩起始位置先小幅后移ꎬ当壁面温度增加到３００Ｋ时ꎬ转捩起始位置开始后移ꎬ当壁面温度增加到７００Ｋ时ꎬ由于湍流区域的减小ꎬ该位置不再发生转捩ꎮ(ａ)Ｕｐｐｅｒｓｕｒｆａｃｅ㊀㊀㊀㊀㊀(ｂ)Ｌｏｗｅｒｓｕｒｆａｃｅ图１１㊀不同位置摩阻系数随壁面温度的变化Ｆｉｇ.１１㊀Ｖａｒｉａｔｉｏｎｏｆｆｒｉｃｔｉｏｎｃｏｅｆｆｉｃｉｅｎｔｗｉｔｈｗａｌｌｔｅｍｐｅｒａｔｕｒｅａｔｄｉｆｆｅｒｅｎｔｌｏｃａｔｉｏｎｓ为进一步分析壁面温度的影响ꎬ本文分别在上下表面湍流区选取一点(０.９ꎬ０.０２９ꎬ０.１４)ꎬ(０.９７ꎬ－０.３４ꎬ０.１２)ꎬ分析边界层湍动能剖面ꎬ结果如图１２所示ꎮ从图中可以看到ꎬ随着壁面温度升高ꎬ边界层厚度先略微变厚ꎬ再变薄ꎬ当壁面温度升高到７００Ｋ时ꎬ边界层厚度迅速降低ꎮ这些结果与转捩位置先前移再后移的结论相符合ꎬ因为边界层厚度会影响不稳定波的时间和空间尺度ꎬ边界层厚度低时ꎬ不稳定波增长速度变慢ꎬ延迟转捩发生ꎮ需要指出的是ꎬ仅采用当前使用的方法ꎬ无法从更深层７１气体物理２０２４年㊀第９卷次揭示转捩反转的流动机理ꎬ而须另外借助稳定性分析方法ꎬ例如ꎬ使用ｅＮ方法开展基于模态的稳定性研究ꎮ文献[３３]采用该手段研究了大掠角平板钝三角翼随壁温比变化出现转捩反转的内在机理:壁温比升高促进横流模态和第１模态扰动增长ꎬ抑制第２模态发展ꎬ在第１㊁２模态联合作用影响下ꎬ出现转捩反转现象ꎮ我们将在后续开展进一步研究ꎮ(ａ)Ｕｐｐｅｒｓｕｒｆａｃｅ(ｂ)Ｌｏｗｅｒｓｕｒｆａｃｅ图１２㊀不同位置湍动能剖面随壁面温度的变化Ｆｉｇ.１２㊀Ｖａｒｉａｔｉｏｎｏｆｔｕｒｂｕｌｅｎｔｋｉｎｅｔｉｃｅｎｅｒｇｙｗｉｔｈｗａｌｌｔｅｍｐｅｒａｔｕｒｅａｔｄｉｆｆｅｒｅｎｔｌｏｃａｔｉｏｎｓ６㊀结论针对ＨｙＴＲＶ转捩问题ꎬ在Ｍａｃｈ数Ｍａ＝３~８ꎬ壁面温度Ｔ＝１５０~９００Ｋ的条件下ꎬ基于课题组自研软件ꎬ对γ￣Ｒｅ~θｔ转捩模型和ＳＳＴ湍流模型进行了高超声速修正ꎬ研究了Ｍａｃｈ数和壁面温度对ＨｙＴＲＶ转捩的影响ꎬ得出以下结论:１)经过高超声速修正后的γ￣Ｒｅ~θｔ转捩模型和ＳＳＴ湍流模型能够较为准确地预测ＨｙＴＲＶ转捩位置ꎬ并且湍流边界层区域形状与实验结果基本一致ꎻＨｙＴＲＶ存在多个不同的转捩区域ꎬ上表面两侧转捩区域分布在两侧边缘位置ꎬ下表面两侧转捩区域分布在中心线两侧ꎮ２)Ｍａｃｈ数的增加会导致上下表面转捩起始位置均大幅后移ꎬ湍流区大幅缩小ꎬ但当Ｍａｃｈ数增加到８时ꎬ湍流区仍然存在ꎬ并没有消失ꎻ上表面层流区摩阻不断增加ꎬ下表面湍流区摩阻不断减小ꎮ３)壁面温度的增加会导致上下表面转捩起始位置先前移ꎬ再后移ꎬ这与边界层厚度变化规律一致ꎬ当壁面温度增加到７００Ｋ时ꎬ下表面湍流区已经基本消失ꎬ当壁面温度增加到９００Ｋ时ꎬ上表面湍流区也基本消失ꎻ上表面在层流区域的摩阻系数逐渐增大ꎬ在湍流区的摩阻系数逐渐减小ꎮ致谢㊀感谢中国空气动力研究与发展中心和空天飞行空气动力科学与技术全国重点实验室提供的ＨｙＴＲＶ模型数据和实验数据ꎮ参考文献(Ｒｅｆｅｒｅｎｃｅｓ)[１]㊀ＯｂｅｒｉｎｇＩＩＩＨꎬＨｅｉｎｒｉｃｈｓＲＬ.Ｍｉｓｓｉｌｅｄｅｆｅｎｓｅｆｏｒｇｒｅａｔｐｏｗｅｒｃｏｎｆｌｉｃｔ:ｏｕｔｍａｎｅｕｖｅｒｉｎｇｔｈｅＣｈｉｎａｔｈｒｅａｔ[Ｊ].Ｓｔｒａ￣ｔｅｇｉｃＳｔｕｄｉｅｓＱｕａｒｔｅｒｌｙꎬ２０１９ꎬ３(４):３７￣５６. [２]ＣｈｅｎｇＣꎬＷｕＪＨꎬＺｈａｎｇＹＬꎬｅｔａｌ.Ａｅｒｏｄｙｎａｍｉｃｓａｎｄｄｙｎａｍｉｃｓｔａｂｉｌｉｔｙｏｆｍｉｃｒｏ￣ａｉｒ￣ｖｅｈｉｃｌｅｗｉｔｈｆｏｕｒｆｌａｐｐｉｎｇｗｉｎｇｓｉｎｈｏｖｅｒｉｎｇｆｌｉｇｈｔ[Ｊ].ＡｄｖａｎｃｅｓｉｎＡｅｒｏｄｙｎａｍｉｃｓꎬ２０２０ꎬ２(３):５.[３]ＢｅｒｔｉｎＪＪꎬＣｕｍｍｉｎｇｓＲＭ.Ｆｉｆｔｙｙｅａｒｓｏｆｈｙｐｅｒｓｏｎｉｃｓ:ｗｈｅｒｅｗｅᶄｖｅｂｅｅｎꎬｗｈｅｒｅｗｅᶄｒｅｇｏｉｎｇ[Ｊ].ＰｒｏｇｒｅｓｓｉｎＡｅｒｏｓｐａｃｅＳｃｉｅｎｃｅｓꎬ２００３ꎬ３９(６/７):５１１￣５３６. [４]陈坚强ꎬ涂国华ꎬ张毅锋ꎬ等.高超声速边界层转捩研究现状与发展趋势[Ｊ].空气动力学学报ꎬ２０１７ꎬ３５(３):３１１￣３３７.ＣｈｅｎＪＱꎬＴｕＧＨꎬＺｈａｎｇＹＦꎬｅｔａｌ.Ｈｙｐｅｒｓｎｏｎｉｃｂｏｕｎｄａｒｙｌａｙｅｒｔｒａｎｓｉｔｉｏｎ:ｗｈａｔｗｅｋｎｏｗꎬｗｈｅｒｅｓｈａｌｌｗｅｇｏ[Ｊ].ＡｃｔａＡｅｒｏｄｙｎａｍｉｃａＳｉｎｉｃａꎬ２０１７ꎬ３５(３):３１１￣３３７(ｉｎＣｈｉｎｅｓｅ).[５]段毅ꎬ姚世勇ꎬ李思怡ꎬ等.高超声速边界层转捩的若干问题及工程应用研究进展综述[Ｊ].空气动力学学报ꎬ２０２０ꎬ３８(２):３９１￣４０３.ＤｕａｎＹꎬＹａｏＳＹꎬＬｉＳＹꎬｅｔａｌ.Ｒｅｖｉｅｗｏｆｐｒｏｇｒｅｓｓｉｎｓｏｍｅｉｓｓｕｅｓａｎｄｅｎｇｉｎｅｅｒｉｎｇａｐｐｌｉｃａｔｉｏｎｏｆｈｙｐｅｒｓｏｎｉｃｂｏｕｎｄａｒｙｌａｙｅｒｔｒａｎｓｉｔｉｏｎ[Ｊ].ＡｃｔａＡｅｒｏｄｙｎａｍｉｃａＳｉｎｉｃａꎬ２０２０ꎬ３８(２):３９１￣４０３(ｉｎＣｈｉｎｅｓｅ).[６]ＺｈａｎｇＹＦꎬＺｈａｎｇＹＲꎬＣｈｅｎＪＱꎬｅｔａｌ.Ｎｕｍｅｒｉｃａｌｓｉ￣８１第２期章录兴ꎬ等:Ｍａｃｈ数和壁面温度对ＨｙＴＲＶ边界层转捩的影响ｍｕｌａｔｉｏｎｓｏｆｈｙｐｅｒｓｏｎｉｃｂｏｕｎｄａｒｙｌａｙｅｒｔｒａｎｓｉｔｉｏｎｂａｓｅｄｏｎｔｈｅｆｌｏｗｓｏｌｖｅｒｃｈａｎｔ２.０[Ｒ].ＡＩＡＡ２０１７￣２４０９ꎬ２０１７. [７]ＫｒａｕｓｅＭꎬＢｅｈｒＭꎬＢａｌｌｍａｎｎＪ.Ｍｏｄｅｌｉｎｇｏｆｔｒａｎｓｉｔｉｏｎｅｆｆｅｃｔｓｉｎｈｙｐｅｒｓｏｎｉｃｉｎｔａｋｅｆｌｏｗｓｕｓｉｎｇａｃｏｒｒｅｌａｔｉｏｎ￣ｂａｓｅｄｉｎｔｅｒｍｉｔｔｅｎｃｙｍｏｄｅｌ[Ｒ].ＡＩＡＡ２００８￣２５９８ꎬ２００８. [８]ＹｉＭＲꎬＺｈａｏＨＹꎬＬｅＪＬ.Ｈｙｐｅｒｓｏｎｉｃｎａｔｕｒａｌａｎｄｆｏｒｃｅｄｔｒａｎｓｉｔｉｏｎｓｉｍｕｌａｔｉｏｎｂｙｃｏｒｒｅｌａｔｉｏｎ￣ｂａｓｅｄｉｎｔｅｒｍｉｔ￣ｔｅｎｃｙ[Ｒ].ＡＩＡＡ２０１７￣２３３７ꎬ２０１７.[９]向星皓ꎬ张毅锋ꎬ袁先旭ꎬ等.Ｃ￣γ￣Ｒｅθ高超声速三维边界层转捩预测模型[Ｊ].航空学报ꎬ２０２１ꎬ４２(９):６２５７１１.ＸｉａｎｇＸＨꎬＺｈａｎｇＹＦꎬＹｕａｎＸＸꎬｅｔａｌ.Ｃ￣γ￣Ｒｅθｍｏｄｅｌｆｏｒｈｙｐｅｒｓｏｎｉｃｔｈｒｅｅ￣ｄｉｍｅｎｓｉｏｎａｌｂｏｕｎｄａｒｙｌａｙｅｒｔｒａｎｓｉｔｉｏｎｐｒｅｄｉｃｔｉｏｎ[Ｊ].ＡｃｔａＡｅｒｏｎａｕｔｉｃａｅｔＡｓｔｒｏｎａｕｔｉｃａＳｉｎｉｃａꎬ２０２１ꎬ４２(９):６２５７１１(ｉｎＣｈｉｎｅｓｅ).[１０]ＭｃＤａｎｉｅｌＲＤꎬＮａｎｃｅＲＰꎬＨａｓｓａｎＨＡ.Ｔｒａｎｓｉｔｉｏｎｏｎｓｅｔｐｒｅｄｉｃｔｉｏｎｆｏｒｈｉｇｈ￣ｓｐｅｅｄｆｌｏｗ[Ｊ].ＪｏｕｒｎａｌｏｆＳｐａｃｅｃｒａｆｔａｎｄＲｏｃｋｅｔｓꎬ２０００ꎬ３７(３):３０４￣３０９.[１１]ＰａｐｐＪＬꎬＤａｓｈＳＭ.Ｒａｐｉｄｅｎｇｉｎｅｅｒｉｎｇａｐｐｒｏａｃｈｔｏｍｏｄｅｌｉｎｇｈｙｐｅｒｓｏｎｉｃｌａｍｉｎａｒ￣ｔｏ￣ｔｕｒｂｕｌｅｎｔｔｒａｎｓｉｔｉｏｎａｌｆｌｏｗｓ[Ｊ].ＪｏｕｒｎａｌｏｆＳｐａｃｅｃｒａｆｔａｎｄＲｏｃｋｅｔｓꎬ２００５ꎬ４２(３):４６７￣４７５.[１２]ＪｕｌｉａｎｏＴＪꎬＳｃｈｎｅｉｄｅｒＳＰ.ＩｎｓｔａｂｉｌｉｔｙａｎｄｔｒａｎｓｉｔｉｏｎｏｎｔｈｅＨＩＦｉＲＥ￣５ｉｎａＭａｃｈ￣６ｑｕｉｅｔｔｕｎｎｅｌ[Ｒ].ＡＩＡＡ２０１０￣５００４ꎬ２０１０.[１３]杨云军ꎬ马汉东ꎬ周伟江.高超声速流动转捩的数值研究[Ｊ].宇航学报ꎬ２００６ꎬ２７(１):８５￣８８.ＹａｎｇＹＪꎬＭａＨＤꎬＺｈｏｕＷＪ.Ｎｕｍｅｒｉｃａｌｒｅｓｅａｒｃｈｏｎｓｕｐｅｒｓｏｎｉｃｆｌｏｗｔｒａｎｓｉｔｉｏｎ[Ｊ].ＪｏｕｒｎａｌｏｆＡｓｔｒｏｎａｕｔｉｃｓꎬ２００６ꎬ２７(１):８５￣８８(ｉｎＣｈｉｎｅｓｅ).[１４]袁先旭ꎬ何琨ꎬ陈坚强ꎬ等.ＭＦ￣１模型飞行试验转捩结果初步分析[Ｊ].空气动力学学报ꎬ２０１８ꎬ３６(２):２８６￣２９３.ＹｕａｎＸＸꎬＨｅＫꎬＣｈｅｎＪＱꎬｅｔａｌ.ＰｒｅｌｉｍｉｎａｒｙｔｒａｎｓｉｔｉｏｎｒｅｓｅａｒｃｈａｎａｌｙｓｉｓｏｆＭＦ￣１[Ｊ].ＡｃｔａＡｅｒｏｄｙｎａｍｉｃａＳｉｎｉｃａꎬ２０１８ꎬ３６(２):２８６￣２９３(ｉｎＣｈｉｎｅｓｅ).[１５]陈坚强ꎬ涂国华ꎬ万兵兵ꎬ等.ＨｙＴＲＶ流场特征与边界层稳定性特征分析[Ｊ].航空学报ꎬ２０２１ꎬ４２(６):１２４３１７.ＣｈｅｎＪＱꎬＴｕＧＨꎬＷａｎＢＢꎬｅｔａｌ.Ｃｈａｒａｃｔｅｒｉｓｔｉｃｓｏｆｆｌｏｗｆｉｅｌｄａｎｄｂｏｕｎｄａｒｙ￣ｌａｙｅｒｓｔａｂｉｌｉｔｙｏｆＨｙＴＲＶ[Ｊ].ＡｃｔａＡｅｒｏｎａｕｔｉｃａｅｔＡｓｔｒｏｎａｕｔｉｃａＳｉｎｉｃａꎬ２０２１ꎬ４２(６):１２４３１７(ｉｎＣｈｉｎｅｓｅ).[１６]陈坚强ꎬ刘深深ꎬ刘智勇ꎬ等.用于高超声速边界层转捩研究的标模气动布局及设计方法.中国:１０９９６９３７４Ｂ[Ｐ].２０２１￣０５￣１８.ＣｈｅｎＪＱꎬＬｉｕＳＳꎬＬｉｕＺＹꎬｅｔａｌ.Ｓｔａｎｄａｒｄｍｏｄｅｌａｅｒｏ￣ｄｙｎａｍｉｃｌａｙｏｕｔａｎｄｄｅｓｉｇｎｍｅｔｈｏｄｆｏｒｈｙｐｅｒｓｏｎｉｃｂｏｕｎｄａｒｙｌａｙｅｒｔｒａｎｓｉｔｉｏｎｒｅｓｅａｒｃｈ.ＣＮꎬ１０９９６９３７４Ｂ[Ｐ].２０２１￣０５￣１８(ｉｎＣｈｉｎｅｓｅ).[１７]ＬｉｕＳＳꎬＹｕａｎＸＸꎬＬｉｕＺＹꎬｅｔａｌ.Ｄｅｓｉｇｎａｎｄｔｒａｎｓｉｔｉｏｎｃｈａｒａｃｔｅｒｉｓｔｉｃｓｏｆａｓｔａｎｄａｒｄｍｏｄｅｌｆｏｒｈｙｐｅｒｓｏｎｉｃｂｏｕｎｄａｒｙｌａｙｅｒｔｒａｎｓｉｔｉｏｎｒｅｓｅａｒｃｈ[Ｊ].ＡｃｔａＭｅｃｈａｎｉｃａＳｉｎｉｃａꎬ２０２１ꎬ３７(１１):１６３７￣１６４７.[１８]ＣｈｅｎＸꎬＤｏｎｇＳＷꎬＴｕＧＨꎬｅｔａｌ.Ｂｏｕｎｄａｒｙｌａｙｅｒｔｒａｎ￣ｓｉｔｉｏｎａｎｄｌｉｎｅａｒｍｏｄａｌｉｎｓｔａｂｉｌｉｔｉｅｓｏｆｈｙｐｅｒｓｏｎｉｃｆｌｏｗｏｖｅｒａｌｉｆｔｉｎｇｂｏｄｙ[Ｊ].ＪｏｕｒｎａｌｏｆＦｌｕｉｄＭｅｃｈａｎｉｃｓꎬ２０２２ꎬ９３８(４０８):Ａ８.[１９]ＱｉＨꎬＬｉＸＬꎬＹｕＣＰꎬｅｔａｌ.Ｄｉｒｅｃｔｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎｏｆｈｙｐｅｒｓｏｎｉｃｂｏｕｎｄａｒｙｌａｙｅｒｔｒａｎｓｉｔｉｏｎｏｖｅｒａｌｉｆｔｉｎｇ￣ｂｏｄｙｍｏｄｅｌＨｙＴＲＶ[Ｊ].ＡｄｖａｎｃｅｓｉｎＡｅｒｏｄｙｎａｍｉｃｓꎬ２０２１ꎬ３(１):３１.[２０]万兵兵ꎬ陈曦ꎬ陈坚强ꎬ等.三维边界层转捩预测ＨｙＴＥＮ软件在高超声速典型标模中的应用[Ｊ].空天技术ꎬ２０２３(１):１５０￣１５８.ＷａｎＢＢꎬＣｈｅｎＸꎬＣｈｅｎＪＱꎬｅｔａｌ.ＡｐｐｌｉｃａｔｉｏｎｓｏｆＨｙＴＥＮｓｏｆｔｗａｒｅｆｏｒｐｒｅｄｉｃｔｉｎｇｔｈｒｅｅ￣ｄｉｍｅｎｓｉｏｎａｌｂｏｕｎｄａｒｙ￣ｌａｙｅｒｔｒａｎｓｉｔｉｏｎｉｎｔｙｐｉｃａｌｈｙｐｅｒｓｏｎｉｃｍｏｄｅｌｓ[Ｊ].ＡｅｒｏｓｐａｃｅＴｅｃｈｎｏｌｏｇｙꎬ２０２３(１):１５０￣１５８(ｉｎＣｈｉｎｅｓｅ). [２１]ＭｅｎｔｅｒＦＲꎬＬａｎｇｔｒｙＲＢꎬＬｉｋｋｉＳＲꎬｅｔａｌ.Ａｃｏｒｒｅｌａｔｉｏｎ￣ｂａｓｅｄｔｒａｎｓｉｔｉｏｎｍｏｄｅｌｕｓｉｎｇｌｏｃａｌｖａｒｉａｂｌｅｓ ＰａｒｔⅠ:ｍｏｄｅｌｆｏｒｍｕｌａｔｉｏｎ[Ｊ].ＪｏｕｒｎａｌｏｆＴｕｒｂｏｍａｃｈｉｎｅｒｙꎬ２００６ꎬ１２８(３):４１３￣４２２.[２２]ＬａｎｇｔｒｙＲＢꎬＭｅｎｔｅｒＦＲꎬＬｉｋｋｉＳＲꎬｅｔａｌ.Ａｃｏｒｒｅｌａｔｉｏｎ￣ｂａｓｅｄｔｒａｎｓｉｔｉｏｎｍｏｄｅｌｕｓｉｎｇｌｏｃａｌｖａｒｉａｂｌｅｓ ＰａｒｔⅡ:ｔｅｓｔｃａｓｅｓａｎｄｉｎｄｕｓｔｒｉａｌａｐｐｌｉｃａｔｉｏｎｓ[Ｊ].ＪｏｕｒｎａｌｏｆＴｕｒ￣ｂｏｍａｃｈｉｎｅｒｙꎬ２００６ꎬ１２８(３):４２３￣４３４.[２３]ＬａｎｇｔｒｙＲＢꎬＭｅｎｔｅｒＦＲ.Ｃｏｒｒｅｌａｔｉｏｎ￣ｂａｓｅｄｔｒａｎｓｉｔｉｏｎｍｏｄｅｌｉｎｇｆｏｒｕｎｓｔｒｕｃｔｕｒｅｄｐａｒａｌｌｅｌｉｚｅｄｃｏｍｐｕｔａｔｉｏｎａｌｆｌｕｉｄｄｙｎａｍｉｃｓｃｏｄｅｓ[Ｊ].ＡＩＡＡＪｏｕｒｎａｌꎬ２００９ꎬ４７(１２):２８９４￣２９０６.[２４]孟德虹ꎬ张玉伦ꎬ王光学ꎬ等.γ￣Ｒｅθ转捩模型在二维低速问题中的应用[Ｊ].航空学报ꎬ２０１１ꎬ３２(５):７９２￣８０１.ＭｅｎｇＤＨꎬＺｈａｎｇＹＬꎬＷａｎｇＧＸꎬｅｔａｌ.Ａｐｐｌｉｃａｔｉｏｎｏｆγ￣Ｒｅθｔｒａｎｓｉｔｉｏｎｍｏｄｅｌｔｏｔｗｏ￣ｄｉｍｅｎｓｉｏｎａｌｌｏｗｓｐｅｅｄｆｌｏｗｓ[Ｊ].ＡｃｔａＡｅｒｏｎａｕｔｉｃａｅｔＡｓｔｒｏｎａｕｔｉｃａＳｉｎｉｃａꎬ２０１１ꎬ３２(５):７９２￣８０１(ｉｎＣｈｉｎｅｓｅ).[２５]牟斌ꎬ江雄ꎬ肖中云ꎬ等.γ￣Ｒｅθ转捩模型的标定与应用[Ｊ].空气动力学学报ꎬ２０１３ꎬ３１(１):１０３￣１０９.ＭｏｕＢꎬＪｉａｎｇＸꎬＸｉａｏＺＹꎬｅｔａｌ.Ｉｍｐｌｅｍｅｎｔａｔｉｏｎａｎｄｃａｌｉｂｅｒａｔｉｏｎｏｆγ￣Ｒｅθｔｒａｎｓｉｔｉｏｎｍｏｄｅｌ[Ｊ].ＡｃｔａＡｅｒｏｄｙ￣ｎａｍｉｃａＳｉｎｉｃａꎬ２０１３ꎬ３１(１):１０３￣１０９(ｉｎＣｈｉｎｅｓｅ). [２６]郭隽ꎬ刘丽平ꎬ徐晶磊ꎬ等.γ￣Ｒｅ~θｔ转捩模型在跨声速涡轮叶栅中的应用[Ｊ].推进技术ꎬ２０１８ꎬ３９(９):１９９４￣２００１.９１气体物理２０２４年㊀第９卷ＧｕｏＪꎬＬｉｕＬＰꎬＸｕＪＬꎬｅｔａｌ.Ａｐｐｌｉｃａｔｉｏｎｏｆγ￣Ｒｅ~θｔｔｒａｎ￣ｓｉｔｉｏｎｍｏｄｅｌｉｎｔｒａｎｓｏｎｉｃｔｕｒｂｉｎｅｃａｓｃａｄｅｓ[Ｊ].ＪｏｕｒｎａｌｏｆＰｒｏｐｕｌｓｉｏｎＴｅｃｈｎｏｌｏｇｙꎬ２０１８ꎬ３９(９):１９９４￣２００１(ｉｎＣｈｉｎｅｓｅ).[２７]郑赟ꎬ李虹杨ꎬ刘大响.γ￣Ｒｅθ转捩模型在高超声速下的应用及分析[Ｊ].推进技术ꎬ２０１４ꎬ３５(３):２９６￣３０４.ＺｈｅｎｇＹꎬＬｉＨＹꎬＬｉｕＤＸ.Ａｐｐｌｉｃａｔｉｏｎａｎｄａｎａｌｙｓｉｓｏｆγ￣Ｒｅθｔｒａｎｓｉｔｉｏｎｍｏｄｅｌｉｎｈｙｐｅｒｓｏｎｉｃｆｌｏｗ[Ｊ].ＪｏｕｒｎａｌｏｆＰｒｏｐｕｌｓｉｏｎＴｅｃｈｎｏｌｏｇｙꎬ２０１４ꎬ３５(３):２９６￣３０４(ｉｎＣｈｉ￣ｎｅｓｅ).[２８]孔维萱ꎬ阎超ꎬ赵瑞.γ￣Ｒｅθ模式应用于高速边界层转捩的研究[Ｊ].空气动力学学报ꎬ２０１３ꎬ３１(１):１２０￣１２６.ＫｏｎｇＷＸꎬＹａｎＣꎬＺｈａｏＲ.γ￣Ｒｅθｍｏｄｅｌｒｅｓｅａｒｃｈｆｏｒｈｉｇｈ￣ｓｐｅｅｄｂｏｕｎｄａｒｙｌａｙｅｒｔｒａｎｓｉｔｉｏｎ[Ｊ].ＡｃｔａＡｅｒｏｄｙ￣ｎａｍｉｃａＳｉｎｉｃａꎬ２０１３ꎬ３１(１):１２０￣１２６(ｉｎＣｈｉｎｅｓｅ). [２９]张毅锋ꎬ何琨ꎬ张益荣ꎬ等.Ｍｅｎｔｅｒ转捩模型在高超声速流动模拟中的改进及验证[Ｊ].宇航学报ꎬ２０１６ꎬ３７(４):３９７￣４０２.ＺｈａｎｇＹＦꎬＨｅＫꎬＺｈａｎｇＹＲꎬｅｔａｌ.ＩｍｐｒｏｖｅｍｅｎｔａｎｄｖａｌｉｄａｔｉｏｎｏｆＭｅｎｔｅｒᶄｓｔｒａｎｓｉｔｉｏｎｍｏｄｅｌｆｏｒｈｙｐｅｒｓｏｎｉｃｆｌｏｗｓｉｍｕｌａｔｉｏｎ[Ｊ].ＪｏｕｒｎａｌｏｆＡｓｔｒｏｎａｕｔｉｃｓꎬ２０１６ꎬ３７(４):３９７￣４０２(ｉｎＣｈｉｎｅｓｅ).[３０]ＬａｎｇｔｒｙＲＢꎬＳｅｎｇｕｐｔａＫꎬＹｅｈＤＴꎬｅｔａｌ.Ｅｘｔｅｎｄｉｎｇｔｈｅγ￣Ｒｅθｔｃｏｒｒｅｌａｔｉｏｎｂａｓｅｄｔｒａｎｓｉｔｉｏｎｍｏｄｅｌｆｏｒｃｒｏｓｓｆｌｏｗｅｆｆｅｃｔｓ(Ｉｎｖｉｔｅｄ)[Ｒ].ＡＩＡＡ２０１５￣２４７４ꎬ２０１５. [３１]ＳａｒｋａｒＳ.Ｔｈｅｐｒｅｓｓｕｒｅ￣ｄｉｌａｔａｔｉｏｎｃｏｒｒｅｌａｔｉｏｎｉｎｃｏｍｐｒｅｓｓｉ￣ｂｌｅｆｌｏｗｓ[Ｊ].ＰｈｙｓｉｃｓｏｆＦｌｕｉｄｓＡꎬ１９９２ꎬ４(１２):２６７４￣２６８２.[３２]ＷｉｌｃｏｘＤＣ.Ｄｉｌａｔａｔｉｏｎ￣ｄｉｓｓｉｐａｔｉｏｎｃｏｒｒｅｃｔｉｏｎｓｆｏｒａｄｖａｎｃｅｄｔｕｒｂｕｌｅｎｃｅｍｏｄｅｌｓ[Ｊ].ＡＩＡＡＪｏｕｒｎａｌꎬ１９９２ꎬ３０(１１):２６３９￣２６４６.[３３]马祎蕾ꎬ余平ꎬ姚世勇.壁温对钝三角翼边界层稳定性及转捩影响[Ｊ].空气动力学学报ꎬ２０２０ꎬ３８(６):１０１７￣１０２６.ＭａＹＬꎬＹｕＰꎬＹａｏＳＹ.Ｅｆｆｅｃｔｏｆｗａｌｌｔｅｍｐｅｒａｔｕｒｅｏｎｓｔａｂｉｌｉｔｙａｎｄｔｒａｎｓｉｔｉｏｎｏｆｈｙｐｅｒｓｏｎｉｃｂｏｕｎｄａｒｙｌａｙｅｒｏｎａｂｌｕｎｔｄｅｌｔａｗｉｎｇ[Ｊ].ＡｃｔａＡｅｒｏｄｙｎａｍｉｃａＳｉｎｉｃａꎬ２０２０ꎬ３８(６):１０１７￣１０２６(ｉｎＣｈｉｎｅｓｅ).０２。

数学: 科学的王后和仆人Mathematics: Queen and Servant of Science北京理工大学叶其孝本文的题目是已故的美国科学院院士、著名数学家、数学史学家和科普作家Eric Temple Bell(贝尔, 1883, 02, 07 ~ 1960, 12, 21)于1951年写的一本书的书名Mathematics: Queen and Servant of Science (数学: 科学的王后和仆人). 该书主要是为大学生和非数学领域的人士写的, 介绍纯粹和应用数学的各个方面, 更着重在说明数学科学的极端重要性.The Mathematical Association of America, 1996, 463 pages实际上这是他1931年写的The Queen of the Sciences (科学的王后)和1937年写的The Handmaiden of the Sciences (科学的女仆)这两本通俗数学论著的合一修订扩大版.Eric Temple Bell Alexander Graham Bell (1847 ~ 1922) 按常识的理解, 女王是优美、高雅、无懈可击、至尊至贵的, 在科学中只有纯粹数学才具有这样的特点, 简洁明了的数学定理一经证明就是永恒的真理, 极其优美而且无懈可击;另一方面, 科学和工程的各个分支都在不同程度上大量应用数学, 这时数学科学就是仆人, 这些仆人是否强有力, 用起来是否得心应手是雇佣这些仆人的主人最为关心的事. 事实上, servant这个字本身就有“供人们利用之物, 有用的服务工具”的意思. 毫无疑问, 我们的目的不是为数学争一个好的名分, 而是想说明数学是怎样通过数学建模来解决各种实际问题的; 数学(数学建模)的极端重要性, 以及探讨正确认识和理解数学科学的作用对于发展我国科学技术、经济以及教育, 从而争取在21世纪把我国真正建设成为屹立于世界民族之林的强国，乃至个人事业发展的至关重要性. 当然, 我们也希望说明王后和仆人集于一身并不矛盾. 历史上, 很多特别受人尊敬的科学家, 不仅仅是由于他们的科学成就, 更因为他们的科学成就能够服务于人类.数学是科学的王后, 算术是数学的王后. 她常常放下架子为天文学和其他科学效劳, 但是在所有情况下, 第一位的是她(数学)应尽的责任. (高斯)Mathematics is the Queen of the Sciences, and Arithmetic the Queen of Mathematics. She often condescends to render service to astronomy and other natural sciences, but under all circumstance the first place is her due.— Carl Friedrich Gauss (卡尔·弗里德里希·高斯, 1777, 4, 30 ~ 1855, 2, 23)From: Bell, Eric T., Mathematics: Queen and Servant of Science, MAA, 1951, p.1;Men of Mathematics, Simon and Schuster, New York, 1937, p. xv.***************************************************自古以来，数学的发展始终与科学技术的发展紧密相连，反之亦然. 首先, 我们来看一下导致我们现在这个飞速发展的信息社会的19、20世纪几乎所有重大科学理论的发展和完善过程中数学(数学建模)所起到的不可勿缺的作用.数学研究的成果往往是重大科学发明的催生素(仅就19、20世纪而言, 流体力学、电磁理论、相对论、量子力学、计算机、信息论、控制论、现代经济学、万维网和互联网搜索引擎、生物学、CT、甚至社会政治学领域等). 但是20世纪上半世纪, 数学虽然也直接为工程技术提供一些工具, 但基本方式是间接的: 先促进其他科学的发展, 再由这些科学提供工程原理和设计的基础. 数学是幕后的无名英雄.现在, 数学无处不在, 数学和工程技术之间,在更广阔的范围内和更深刻的程度上, 直接地相互作用着, 极大地推动了科学和工程科学的发展, 也极大地推动了技术的发展. 数学不仅是幕后的无名英雄, 很多方面开始走向“前台”. 但是对数学的极端重要性迄今尚未有共识, 取得共识对加强一个国家的竞争力来说是至关重要的.硬能力―一位美国朋友谈及对未来中国人的看法: 20年后, 中国年轻人会丢了中国人现在的硬能力, 他们崇拜各种明星, 不愿献身科学, 不再以学术研究为荣, 聪明拔尖的学生都去学金融、法律等赚钱的专业; 而美国人因为认识到其硬能力(例如数学)不行, 进行教育改革, 20年后, 不但保持了其软实力即非专业能力的优势, 而且在硬能力上赶上中国人.‖“正在丢失的硬实力”, 鲁鸣, 《青年文摘》2011年第5期动向：美国很多州新办STEM高中, 一些大学开始开设STEM课程等.STEM = Science + Technology + Engineering + Mathematics2012年2月7日公布的美国总统科技顾问委员会给总统的报告，参与超越：培养额外的100万具有科学、技术、工程和数学学位的大学生(Engage to Excel: Producing One Million Additional College Graduates with Degrees in Science, Technology, Engineering, and Mathematics)The Mathematical Sciences in 2025, the National Academies Press, 2013人们使用的数学科学思想、概念和方法的范围在不断扩大的同时，数学科学的用途也在不断扩展. 21世纪的大部分科学与工程将建立在数学科学的基础上.This major expansion in the uses of the mathematical sciences has been paralleled by a broadening in the range of mathematical science ideas and techniques being used. Much of twenty-first century science and engineering is going to be built on a mathematical science foundation, and that foundation must continue to evolve and expand.数学科学是日常生活的几乎每个方面的组成部分.互联网搜索、医疗成像、电脑动画、数值天气预报和其他计算机模拟、所有类型的数字通信、商业和军事中的优化问题以及金融风险的分析——普通公民都从支撑这些应用功能的数学科学的各种进展中获益，这样的例子不胜枚举.The mathematical sciences are part of almost every aspect of everyday life. Internet search, medical imaging, computer animation, numerical weather predictions and othercomputer simulations, digital communications of all types, optimization in business and the military, analyses of financial risks —average citizens all benefit from the mathematical science advances that underpin these capabilities, and the list goes on and on.调查发现:数学科学研究工作正日益成为生物学、医学、社会科学、商业、先进设计、气候、金融、先进材料等许多研究领域不可或缺的重要组成部分. 这种研究工作涉及最广泛意义下数学、统计学和计算综合，以及这些领域与潜在应用领域的相互作用. 所有这些活动对于经济增长、国家竞争力和国家安全都是至关重要的，而且这种事实应该对作为整体的数学科学的资助性质和资助规模产生影响. 数学科学的教育也应该反映数学科学领域的新的状况.Finding: Mathematical sciences work is becoming an increasingly integral and essential component of a growing array of areas of investigation in biology, medicine, social sciences, business, advanced design, climate, finance, advanced materials, and many more. This work involves the integration of mathematics, statistics, and computation in the broadest sense and the interplay of these areas withareas of potential application. All of these activities are crucial to economic growth, national competitiveness, and national security, and this fact should inform both the nature and scale of funding for the mathematical sciences as a whole. Education in the mathematical sciences should also reflect this new stature of the field.****************************************************************为了以下讲述的方便, 我们先来了解一下什么是数学建模.数学模型(Mathematical Model)是用数学符号对一类实际问题或实际发生的现象的(近似的)描述.数学建模(Mathematical Modeling)则是获得该模型并对之求解、验证并得到结论的全过程.数学建模不仅是了解基本规律, 而且从应用的观点来看更重要的是预测和控制所建模的系统的行为的强有力的工具.数学建模是数学用来解决各种实际问题的桥梁.↑→→→→→→→→↓↑↓↑↓↓↑↓←←←←←通不过↓↓通过)定义：数学建模就是上述框图多次执行的过程数学建模的难点观察、分析实际问题, 作出合理的假设, 明确变量和参数, 形成明确的数学问题. 不仅仅是翻译的问题; 涉及的数学问题可能是复杂、困难的, 求解也许涉及深刻的数学方法. 如何作出正确的判断, 寻找合适、简洁的(解析或近似) 解法; 如何验证模型.简言之:合理假设、模型建立、模型求解、解释验证.记住这16个字, 将会终生受用.数学建模的重要作用：源头创新当然数学建模也有局限性, 不能单独包打天下, 因为实际问题是非常复杂的, 需要多学科协同解决.在图灵(A. M. Turing)的文章: The Chemical Basis of Morphogenesis (形态生成的化学基础), Philosophical Transactions of the Royal Society of London (伦敦皇家学会哲学公报), Series B (Biological Sciences),v.237(1952), 37-72.1. 一个胚胎的模型. 成形素本节将描述一个正在生长的胚胎的数学模型. 该模型是一种简化和理想化, 因此是对原问题的篡改. 希望本文论述中保留的一些特征, 就现今的知识状况而言, 是那些最重要的特征.1. A model of the embryo. MorphogensIn this section a mathematical model of the growing embryo will be described. This model will be asimplification and an idealization, and consequently a falsification. It is to be hoped that the features retained for discussion are those of greatest importance in the present state of knowledge.想单靠数学建模本身来解决重大的生物学问题是不可能的，另一方面，想仅仅依靠实验来获得对生物学的合理、完整的理解也是极不可能的. There is no way mathematical modeling can solve major biological problems on its own. On the other hand, it ishighly unlikely that even a reasonably complete understanding could come solely from experiment.—— J. D. Murray, Why Are There No 3-Headed Monsters? Mathematical Modeling in Biology, Notices of the AMS,v. 59 (2012), no. 6, p.793.自古以来公平、公正的竞赛都是培养、选拔人才的重要手段, 科学和数学也不例外.中学生IMO (国际数学奥林匹克(International Mathematical Olympiad), 1959 ～)北美的大学生Putnbam数学竞赛(1938 ～)全国大学生数学竞赛(2010 ～)Mathematical Contest in Modeling (MCM, 1985 ～)美国大学生数学建模竞赛Interdisciplinary Contest in Modeling (ICM, 1999～)美国大学生跨学科建模竞赛China Undergraduate Mathematical Contest in Modeling (CUMCM, 1992～) 中国大学生数学建模竞赛中国大学生参加美国大学生数学建模竞赛情况中国大学生数学建模竞赛情况在以下讲述中涉及物理方面的具体的数学模型 (问题)的叙述和初步讨论可参考《物理学与偏微分方程》, 李大潜、秦铁虎编著, (上册, 1997; 下册, 2000), 高等教育出版社.Seven equations that rule your world (主宰你生活的七个方程式), by Ian Stewart, NewScientist, 13 February 2012.Fourier transformation 2ˆ()()ix f f x e dx πξξ∞--∞=⎰Wave equation 22222u u c t x ∂∂=∂∂ Ma xwell‘s equation110, , 0, H E E E H H c t c t∂∂∇⋅=∇⨯=-∇⋅=∇⨯=∂∂Schrödinger‘s equation ˆψH ψi t∂=∂Ian Stewart, In Pursuit of the Unknown:17 Equations That Changed the World (追求对未知的认识：改变世界的17个方程), Basic Books, March 13, 2012.目录(Contents)Why Equations? /viii1. The squaw on the hippopotamus ——Pythagoras‘sTheorem/12. Shortening the proceedings —— Logarithms/213. Ghosts of departed quantities —— Calculus/354. The system of the world ——Newton‘s Law ofGravity/535. Portent of the ideal world —— The Square Root ofMinus One/736. Much ado about knotting ——Euler‘s Formula forPolyhedra/837. Patterns of chance —— Normal Distribution/1078. Good vibrations —— Wave Equation/1319. Ripples and blips —— Fourier Transform/14910. The ascent of humanity —— Navier-StokesEquation/16511. Wave in the ether ——Maxwell‘s Equations/17912. Law and disorder —— Second Law ofThermodynamics /19513. One thing is absolute —— Relativity/21714. Quantum weirdness —— Schrödinger Equation/24515. Codes, communications, and computers ——Information Theory/26516. The imbalance of nature —— Chaos Theory/28317. The Midas formula —— Black-Scholes Equation/195Where Next?/317Notes/321Illustration Credits/330Index/331相对论Albert Einstein(1879, 3, 14 ～1955, 4, 18)20世纪最伟大的科学成就莫过于Einstein(爱因斯坦)的狭义和广义相对论了, 但是如果没有Minkowski (闵可夫斯基)几何、Riemann(黎曼)于1854年发明的Riemann几何, 以及Cayley(凯莱), Sylvester(西勒维斯特)和Noether(诺特)等数学家发展的不变量理论, Einstein的广义相对论和引力理论就不可能有如此完善的数学表述. Einstein自己也不止一次地说过.早在1905年, 年仅26岁的爱因斯坦就已提出了狭义相对论. 狭义相对论推倒了牛顿力学的质量守恒、能量守恒、质量能量互不相关、时空永恒不变的基本命题. 这是一场真正的科学革命.为了导出狭义相对论，爱因斯坦作出了两个假设：运动的相对性(所有匀速运动都是相对的)和光速为常数(光的运动例外, 它是绝对的). (1)狭义相对性原理,即在所有惯性系中, 物理学定律具有相同的数学表达形式;(2)光速不变原理,真空中光沿各个方向传播的速率都相等，与光源和观察者的运动状态无关.时空不是绝对独立的.由此可以导出一些推论: 相对论坐标变换式和速度变换式, 同时的相对性, 钟慢尺缩效应和质能关系式等.他的好友物理学家P.Ehrenfest指出实际上还蕴涵着第三个假设, 即这两个假设是不矛盾的. 物体运动的相对性和光速的绝对性, 两者之间的相互制约和作用乃是相对论里一切我们不熟悉的时空特征的根源.(部分参阅李新洲:《寻找自然之律--- 20世纪物理学革命》, 上海科技教育出版社, 2001.)1907 年德国数学家H. Minkowski (1864 ～1909) 提出了―Minkowski 空间‖,即把时间和空间融合在一起的四维空间1,3R. Minkowski 几何为Einstein 狭义相对论提供了合适的数学模型.“没有任何客观合理的方法能够把四维连续统分离成三维空间连续统和一维时间连续统. 因此从逻辑上讲, 在四维时空连续统(space- time continuum)中表述自然定律会更令人满意. 相对论在方法上的巨大进步正是建立在这个基础之上的, 这种进步归功于闵可夫斯基(Minkowski).”—Albert Einstein, The Meaning of Relativity, 1922, Princeton University Press. 中译本, 阿尔伯特·爱因斯坦著, 相对论的意义, (普林斯顿科学文库(Princeton Science Library) 1), 郝建纲、刘道军译, 上海科技教育出版社, 2001, p. 27.有了Minkowski 时空模型后, Einstein 又进一步研究引力场理论以建立广义相对论. 1912 年夏他已经概括出新的引力理论的基本物理原理, 但是为了实现广义相对论的目标, 还必须寻求理论的数学结构, Einstein 为此花了 3 年的时间, 最后, 在数学家M. Grossmann 的介绍下学习掌握了发展相对论引力学说所必需的数学工具—以Riemann几何和Ricci, Levi - Civita的绝对微分学, 也就是Einstein 后来所称的张量分析.“根据前面的讨论, 很显然, 如果要表达广义相对论, 就需要对不变量理论以及张量理论加以推广. 这就产生了一个问题, 即要求方程的形式必须对于任意的点变换都是协变的. 在相对论产生以前很久, 数学家们就已经建立了推广的张量演算理论. 黎曼(Riemann)首先把高斯(Gauss)的思路推广到了任意维连续统, 他很有预见性地看到了……进行这种推广的物理意义. 随后, 这个理论以张量微积分的形式得到了发展, 对此里奇(Ricci)和莱维·齐维塔(Tulio Levi-Civita, 1873～1941)做出了重要贡献. ”—阿尔伯特·爱因斯坦著, 相对论的意义, 郝建纲、刘道军译, 上海科技教育出版社, 2001, p. 57.从数学建模的角度看, 广义相对论讨论的中心问题是引力理论, 其基础是以下两个假设: 1. (等效原理)惯性力场与引力场的动力学效应是局部不可分辨的，(或说引力和非惯性系中的惯性力等效)；2. (广义相对性原理) 一切参考系都是平权的，换言之，客观的真实的物理规律应该在任意坐标变换下形式不变——广义协变性(即一切物理定律在所有参考系[无论是惯性的或非惯性的]中都具有相同的形式)。

2021年4月农业机械学报第52卷第4期doi：10.6041/j.issn.1000-1298.2021.04.041考虑非线性摩擦的绳驱动连续体机器人动力学研究齐飞1余世刚1高书苑1陈柏2吴洪涛2（1.常州大学机械与轨道交通学院，常州213016；2.南京航空航天大学机电学院，南京210016）摘要：针对连续体机器人因自身的柔顺性及自由度的冗余性所带来的运动建模复杂和控制精度较低等问题，提出了一种考虑非线性摩擦的绳驱动连续体机器人动力学建模方法。

基于凯恩方法构建了包含惯性力、弹性力、重力及驱动力作用的动力学模型，并根据改进的Capstan方程对绳轮传动系统的力传递特性进行分析，提出了一种基于非线性摩擦补偿的前馈控制策略，最后通过对单节连续体机器人进行仿真和运动实验来验证所建模型的正确性。

结果表明，所建立的运动模型和控制策略能够较为准确地描述连续体机器人的变形特性。

”关键词：连续体机器人；动力学模型；传动摩擦；凯恩方法中图分类号：TP242文献标识码：A文章编号：1000-1298（2021）04-037549OSID：|Dynamic of Cable-driven Continuum Robot withNonlinear Friction ModelQI Fei1SHE Shigang1GAO Shuyuan1CHEN Bai2WU Hongtao2(1.School of Mechanical Engineering and Rail Transit,Changzhou University,Changzhou213016,China2.College of Mechanical and Electrical Engineering,Nanjing University of Aeronautics and Astronautics,Nanjing210016,China)Abstract:According to the complexity of motion modeling and low control accuracy of the continuum robot due to its flexibility and redundancy of freedom，a dynamic modeling method of the cable-driven continuum robot with nonlinear friction was proposed.Firstly，the kinematics model of the continuum robot was derived to analyze the relationship between the end pose and joint parameters,and the velocity kinematics model was presented for the convenience of mechanical analysis.Then，a dynamic model with considering of the inertial force，elastic force，gravity force and driving force，was established based on the Kane's method to analyze the mapping relationship between its motion and the force applied to the robot.Moreover,the force transmission characteristics of the cable—pulley system were analyzed by the improved Capstan equation，and a feedforward control strategy based on the established motion model and the nonlinear friction was presented to compensate the loss of the transmission force.Finally，the simulations and motion experiments were carried out to verify the established model and the proposed control method of the continuum robot.The results showed that the established dynamic model with nonlinear friction can well reflect the deformation characteristics of the continuum robot，which can be used as a reference for the shape reconstruction of the robot in future research.Key words:continuum robot;dynamic model;transmission friction；Kane's method0引言连续体机器人具有独特的柔顺性、安全性和灵活性,能在一些曲折、狭小的工作环境中作业，因此在工业管道、航空发动机及医疗手术等领域得到了广泛研究。

第 36 卷第 5 期2023 年10 月振 动 工 程 学 报Journal of Vibration EngineeringVol. 36 No. 5Oct. 2023多频未知时变扰动下的结构微振动鲁棒自适应控制方昱斌1，朱晓锦2，杨龙飞1，许志超1，田梦楚1，张小兵3（1.南京理工大学智能制造学院，江苏南京 210094； 2.上海大学机电工程及自动化学院，上海 200072；3.南京理工大学能源与动力工程学院，江苏南京 210094）摘要: 本文针对多频窄带未知和时变扰动，基于内模原理和Y⁃K参数化方法，提出一种反馈鲁棒自适应振动的主动控制算法。

该算法通过设计PID中央鲁棒控制器，有效解决了次级通道模型未知情况下的鲁棒控制器参数设计问题。

同时提出一种变步长最小均方（Variable Step Size Least Mean Square，VSSLMS）方法，可以在保证稳态误差的基础上大幅提升收敛速度，并通过系统辨识实验验证了所提VSSLMS方法相较于其他VSSLMS算法在收敛性能上的优越性。

通过结构微振动主动控制实时实验，对比验证了单独采用滤波x最小均方（Least Mean Square，LMS）自适应控制算法、基于LMS算法的鲁棒自适应控制算法和基于VSSLMS算法的鲁棒自适应控制算法的抑振效果。

实验结果表明，本文基于VSSLMS算法的鲁棒自适应控制算法在面向双频正弦窄带扰动以及其频谱、幅值突变情况时，都具有较好的收敛性和鲁棒性。

关键词:振动主动控制；Y⁃K参数化；变步长（VSS）； LMS算法；鲁棒自适应中图分类号: TB535；O327 文献标志码: A 文章编号: 1004-4523（2023）05-1309-09DOI：10.16385/ki.issn.1004-4523.2023.05.015引言在高分遥感探测领域中，卫星结构微振动会引发空间相机产生视线抖动和像移，继而降低成像质量和分辨率，其影响不能忽视［1⁃2］。