弹性特性

关键词： 弹性 飞机 设计 要求

弹性特性（精选八篇）

弹性特性 篇1

弹性支撑对梁动力特性的影响已经得到了较为深入的研究[1,2,3,4].关于拱结构动力特性的研究,主要集中于刚性支撑的情况[5,6,7].具有弹性支撑的拱结构,也有学者进行了研究,但大多是关于内力分析[8]和稳定性分析方面[9,10,11,12,13,14,15].有必要对弹性支撑拱结构的动力特性进行系统研究.当水平弹性支撑的刚度为0时,拱是一个圆弧梁,随着刚度逐渐增加而越来越接近普通的拱结构.因此水平弹性支撑拱结构的动力特性将介于圆弧梁和普通的拱结构之间,寻求水平弹性支撑拱转化为梁和普通拱结构的临界刚度是本文重点解决的问题.

首先用柔度法建立结构自由振动方程,其中考虑了拱脚处集中质量的附加惯性力.用解析法求得结构的柔度矩阵.计算分析了圆弧梁、两铰圆拱和水平弹性支撑圆拱的振动特性.利用圆弧梁与两铰圆拱自振频率和振型内力的区别,提出了临界柔度系数的概念,得到了两铰圆拱与圆弧梁的划分标准.最后分析了水平支撑刚度对拱结构振动内力特性的影响,得出了若干有益的结论.

1 基本方程

1.1 基本假定

(1)基于平面变形假设,变形前垂直于中线的平面,变形后仍垂直于中线;(2)忽略剪切变形的影响;(3)材料是均匀、连续、各向同性的;(4)振动是弹性的.

1.2 振动方程

利用集中质量法将两铰圆拱结构简化为图1所示的质点系.在平面内自由振动时,各质点的切向振动和法向振动均可设为简谐振动,即

式中,ui为切向位移,以绕曲率中心顺时针转动为正;wi为法向位移,以指向极点为正,x为拱脚处集中质量点m的水平位移,向内为正.见图1.

设水平弹性刚度为L为拱的跨度,EI为抗弯刚度,a为系数.

根据柔度法可建立结构自由振动方程

将振动方程(1)代入式(2),可得下列形式

式中,M为质量矩阵,是一个对角矩阵,λ=1/ω2,ω为结构自振圆频率;Y为振型向量,Y=(U1U2···UnW1W2···WnX)T;δ为柔度矩阵.

其中

1.3 柔度矩阵的计算

利用结构力学的方法,可求得结构的柔度矩阵.解方程(3)即可得结构的自振频率和主振型.计算时取图2所示圆弧梁为基本结构,各力以图示方向为正.

图中Hiu,Hjw分别表示切向单位力Fupi=1、法向单位力Fwpj=1作用下弹性支撑拱的水平推力.FuAi和FwAj为A点处的支座反力,FuBi和FwBj为B处的支座反力.

式中,Muip和Nuip,Mwjp和Nwjp,M1p和N1p分别为切向单位力Fupi=1、法向单位力Fwpj=1和单位水平推力H=1作用时圆弧梁的弯矩和轴力;δ11表示单位水平推力H=1作用时圆弧梁支座B处沿H方向的水平位移.

2 水平弹性支撑圆拱振动特性分析

2.1 水平弹性支撑对圆拱自振频率的影响

当水平弹性系数ku→0时,拱脚处的水平推力将趋于0,结构为圆弧简支梁.当ku→∞时,结构为两铰圆拱.圆弧梁、两铰圆拱、弹性支座拱结构的最小频率分别可用下式计算

式中,µl,µg,µk分别为圆弧梁、两铰圆拱和弹性支座拱的自振频率系数.

图3中从上到下共3组曲线,分别表示拱、弹性支座拱和圆弧梁的自振频率系数(考虑轴向变形)与径厚比R/h、矢跨比f/L的关系,其中

由图3可以看出,当径厚比R/h 20时,圆弧梁的最小频率系数只与结构的矢跨比有关,随着矢跨比的增大逐渐减小.对于拱结构而言,当矢跨比f/L>0.1时,最小频率系数随矢跨比的增大逐渐减小,与径厚比没有关系;当0

弹性支撑拱结构的自振频率系数比拱结构小,由此可见水平弹性约束会使拱的自振频率减小.设∆µ=µg-µk为水平弹性约束对拱结构频率系数的影响量,结果见图4.

由图4可见,水平弹性系数ku一定时,对矢跨比和径厚比不同的拱影响不同.图中的曲线可分为上升段和下降段.下降段曲线基本重合,说明当矢跨比大于0.1时,弹性支座对不同径厚比的拱影响是相同的,均随着矢跨比的增大逐渐减小.上升段曲线不重合,矢跨比越小,影响越小.矢跨比一定时,径厚比越大,影响越大.

总的来说,对于工程中一般的拱结构,当矢跨比在0.1左右时,影响最大,对支座的弹性最为敏感.

2.2 两铰圆拱和圆弧梁的自振内力特点

分析了两铰圆拱、圆弧简支梁的振动内力特性.结构参数为:跨度L=20 m,矢高f=4 m,矩形截面A=b×h=0.5 m×0.5 m,弹性模量E=3.15×104N/mm2,钢筋混凝土容重取为2 500 kg/m3.

取集中质量数n=7即可满足小数点后两位的计算精度,计算得14个自振频率和主振型.将振型向量范数归一标准化,即若,则为切向振动,反之为法向振动.

在圆弧简支梁的14个振型中,第1到第8共8个振型以法向振动为主,第9到第14共6个振型以切向振动为主.

两铰圆拱法向振动个数与切向振动个数相等.前7个振型以法向振动为主,后7个振型以切向振动为主.

圆弧梁、拱结构法向、切向振动的内力特点结果见图5,图6所示.图中横坐标表示拱(梁)上的各质点号n,其中“4”为拱顶处.为了便于比较不同振型的内力特点,图中弯矩、轴力振型图的纵坐标分别为:M/ωi2,FN/ωi2.M,FN为在各点惯性力幅值向量miωi2Yi作用下产生的弯矩和轴力,根据叠加原理求得;Yi为第i振型对应的振型向量.

由于圆弧简支梁第1振型的弯矩远大于其他振型的弯矩,故不便在图中标出.

由图5可以看出,法向振动(第2,8振型)产生的弯矩远大于轴力,切向振动(第9,10振型)产生的轴力远大于弯矩.其他法向振动和切向振动也具有类似的特点.由此可以得出圆弧梁的法向振动以弯曲振动为主,而切向振动以拉压振动为主.

目前普遍认为[6]:圆拱结构的切向振动以拉压振动为主,法向振动则以弯曲振动为主.由图6可以看出,圆弧拱结构切向振动(第8,9振型)产生的轴力较大,弯矩相对小得多,因此切向振动是以拉压振动为主的;只有反对称的法向振动才以弯曲振动为主(第1,3振型),而对称的法向振动则是弯曲振动和拉压振动并重(第2振型),或仍以拉压振动为主(第4,5振型).因此笼统地说圆拱结构的法向振动是以弯曲振动为主是不确切的.

2.3 拱与圆弧梁的划分标准

根据前面的研究,拱与圆弧梁之间动力特性区别主要表现在几个方面:(1)拱结构的频率大于圆弧梁的频率;(2)拱结构的法向振动个数与切向振动个数相等,而圆弧梁的法向振动个数大于切向振动个数;(3)拱结构反对称法向振动以弯曲振动为主,而对称的法向振动则是拉压振动和弯曲振动并重.圆弧梁的法向振动以弯曲振动为主.

水平弹性支座拱结构介于两铰圆拱和圆弧简支梁之间,因此有

将叫做柔度系数.

当ξ=0时,为两铰圆拱;当ξ=1时,为拱形梁;0<ξ<1时为弹性拱结构;因此可以按照ξ的大小来划分拱与拱形梁.

经过大量计算发现对于一般的水平弹性支撑拱结构(f/L>0.1),当水平弹性约束使拱结构自振频率减小不超过1%时,具有拱的振动特点,定义这时的弹性系数为临界弹性系数,用kgcr0.01表示,对应的柔度系数为临界柔度系数,用ξgcr0.01表示.

类似地定义ξLcr0.01为圆弧梁临界柔度系数,即当弹性支座拱的柔度系数ξξLcr0.01时,具有圆弧梁的振动特点,按照圆弧梁计算自振频率,误差不超过1%.对应的临界弹性系数为kLcr0.01.

经过计算,对于一般的拱结构(f/L 0.1),临界系数如表1所示.

由表1可以看出,ξgcr0.01在0.067 1∼0.078 0之间变化,ξLcr0.01在0.972∼0.982之间变化,因此可以取ξgcr0.01=0.070,ξLcr0.01=0.980.

当0ξ0.070时,可以按照拱结构计算;当0.980ξ1时,按照圆弧梁计算;当0.070<ξ<0.980时,按照水平弹性支撑拱结构计算.

2.4 水平弹性支撑对拱结构振动内力特性的影响

选取的结构参数与前面相同.不同柔度系数拱结构计算结果见图7.图中横坐标表示拱(梁)上的各质点号n,其中“4”为拱顶处.

弹性支座不仅会降低拱结构的自振频率,同时会改变其振动形态.由图7(a)可以看出,随着柔度系数ξ的增大,即弹性支撑刚度逐渐减小,结构的振型弯矩会逐渐增大,使结构受力不利.由图7(b)看出,当ξ=0.9时,弹性支撑拱第2弯矩振型图与梁完全重合.

由图7(c)可以看出,第4振型拱结构以拉压振动为主.水平弹性支撑拱在不同柔度系数ξ=0.1∼0.9时的振型曲线与梁完全重合,均以弯曲振动为主.这就说明即使弹性支撑产生微小位移也会显著改变拱结构的振动形态,使其具有梁的振动特点.这在高阶振型中影响更加显著.

3 结论

本文主要分析了水平弹性支撑对圆拱结构动力特性的影响,得到如下结论:

(1)n质点体系的结构有2n个自振频率和主振型.两铰圆拱法向振动个数和切向振动个数相等,切向振动以拉压振动为主,反对称法向振动以弯曲振动为主,对称的法向振动则是弯曲振动和拉压振动并重,或仍以拉压振动为主.圆弧梁法向振动的个数大于切向振动的个数,切向振动以拉压振动为主,法向振动以弯曲振动为主.

(2)水平弹性支撑会使拱结构的自振频率减小,且支撑的弹性系数越小,频率减小得越多.对于径厚比R/h 20的拱结构,当矢跨比f/L=0.1左右时,影响最为显著.

(3)利用圆弧梁与圆拱结构振动特性的区别,提出了柔度系数ξ的概念,通过计算,得到了水平弹性支撑拱转化为拱结构和圆弧梁的临界柔度系数ξgcr0.01=0.070,ξLcr0.01=0.980,以及对应的临界弹性系数kgcr0.01,kLcr0.01.这为水平弹性支撑拱结构的动力设计带来了方便,即当ξ≤ξgcr0.01时,可以不用考虑弹性支撑的影响.

(4)水平弹性支撑不仅使拱结构的自振频率减小,而且会显著改变其振动形态.即使柔度系数很小,比如ξ=0.1,水平弹性支撑拱也将不再具有拱的振动特点,尤其高阶振型将完全按照梁的特点振动,即会使拱结构由拉压振动变为弯曲振动

摘要：用柔度法建立了水平弹性支撑拱结构的自由振动方程,考虑了拱脚处集中质量的附加惯性力.计算分析了水平弹性支撑对两铰圆拱固有特性的影响,水平弹性支撑会使拱结构的自振频率减小,当拱结构的矢跨比为0.1左右时影响最为显著,同时还会改变拱结构的振动形态,尤其在高阶振型中将完全按照梁的特点振动.分析了圆弧梁与两铰圆拱的振动内力特点,提出了柔度系数的概念,经过计算得到了水平弹性支撑拱转化为两铰圆拱和圆弧梁的临界柔度系数以及对应的临界刚度系数.

弹性特性 篇2

三维机翼的静气动弹性特性数值模拟研究

发展了一种计算流体动力学(CFD)和计算结构动力学(CSD)的耦合计算方法,对三维机翼的静气动弹性进行了数值模拟研究.采用三维欧拉方程为控制方程基于直角网格计算气动力,并耦合结构静平衡方程进行静气动弹性数值模拟,设计了CFD/ CSD耦合计算的数据交换的方式.以M6机翼作为算例,进行了机翼静气动弹性的数值模拟,计算结果表明所发展的`三维机翼静气动弹性数值模拟方法是合理可行的,并可作为机翼结构优化设计和考虑结构弹性变形影响的气动外形优化设计的基础.

作 者：王晓江 夏露 詹浩 刘付龙 WANG Xiao-jiang XIA Lu ZHAN Hao LIU Fu-long 作者单位：西北工业大学,翼型叶栅空气动力学国防科技重点研究室,陕西,西安,710072刊 名：航空计算技术 ISTIC英文刊名：AERONAUTICAL COMPUTING TECHNIQUE年，卷(期)：40(2)分类号：V211.41关键词：静气动弹性 CFD/CSD 欧拉方程 数据交换 细平板样条法

弹性特性 篇3

摘要：首先考虑各组件陀螺效应推导了旋转状态下转轴一轮盘及叶片的连续体能量方程并进行离散化和组集处理，获得了柔性转子-叶片耦合系统的运动微分方程。然后分析了叶片数和轮盘位置对系统固有频率的影响。通过对比分析发现，由于叶片和转轴的耦合出现了由叶片前两阶模态主导的转子-叶片耦合模态和叶片-叶片耦合模态，其中包含Nb-2（Nb为叶片数目）个重复的叶片一叶片耦合模态，2个重复的转子-叶片耦合模态；随着叶片数的增加，转子模态固有频率线性减小，叶片-叶片耦合模态的固有频率不变，转子-叶片耦合模态的固有频率线性增大；转子一叶片耦合模态和转子前两阶模态对应的固有频率均受轮盘位置影响，影响趋势关于转轴中点对称；其中，转子-叶片耦合模态对应的一阶固有频率在轮盘从轴端到中点过程中逐渐增大，在轴中点时达到最大，二阶固有频率在轮盘从轴端到中点过程中先增大后减小，在轴中点时达到最小；转子模态固有频率变化则与转子-叶片耦合模态完全相反，但不同的是其基本不受叶片数影响。

关键词：转子-叶片耦合系统；固有频率；数值分析；叶片数目；轮盘位置

引言

在大多数旋转机械中，转子系统属于柔性系统，而转子-叶片更是存在耦合关系。较为常见的分析方法是单独对叶片或简化转子系统的动力学特性进行分析，针对转子-叶片耦合系统的研究较少。随着科研工作不断深入，国内外众多学者开始对转子-叶片耦合系统的振动问题进行研究。Omprakash等研究了叶片偏角和叶片扭转角对叶盘-叶片系统固有频率的影响。Huang等介绍了一种新的方法来分析轴扭转和叶片弯曲的旋转轴盘叶片单元之问的动态耦合，并且研究了转轴转速对模态频率的影响。Guo等提出了一种转子、轮系耦合扭转振动的分析方法，将原来大规模耦合振动系统降阶为低阶等效的小规模耦合系统，进而求出其频率和振型。Yang等研究了纵向支撑的柔性和叶片偏角对固有特性的影响。Turhan等研究了轴扭转和叶片弯曲之问的耦合振动问题。Yang等研究了转子系统中的叶片弯曲、叶盘变形和转轴扭转之问的耦合振动问题。wang等研究了油膜力作用下的叶片-转轴耦合振动系统的非线性动力学行为。chiu等建立了考虑转轴扭转的轴-轮盘-叶片耦合系统，研究了带有失谐叶片的转子系统中的耦合振动问题。Chiu等研究了多叶盘转子系统中叶片弯曲、转轴扭转之问的耦合振动问题。sinou等采用谐波平衡法和随机有限元法相结合的方法，研究了旋转轴横向裂纹的影响。chouksey等采用的分析方法，探讨了内部转子材料的阻尼和流体膜力对柔性转子-叶片系统模态特性的影响。Chiu等分析性地研究了一种带分组叶片的多盘转子系统中转轴扭转、叶片弯曲、拉筋之问耦合振动对耦合系统振动的影响。Ma等建立了一个考虑轴的横向和扭转变形、转子陀螺效应、离心刚化、旋转软化和叶片科氏力新的转子-叶片动态模型。

本文以柔性转子-叶片耦合系统为研究對象，推导出旋转状态下的柔性转子-叶片耦合系统的振动微分方程，分析了转子-叶片耦合系统的固有特性，并通过数值模拟研究了叶片数、转盘位置对系统固有频率的影响规律。

1.动力学模型的建立

后掠机翼跨音速静气动弹性特性研究 篇4

通过对弹簧近似方法[3,4]的改进, 使弹簧近似方法在保持原有计算效率的前提下, 网格变形能力和网格质量得到大大的提高。最后, 应用本文发展的非结构动网格生成方法, 以Euler方程为控制方程, 计算弹性后掠机翼飞行时所受气动力载荷, 再耦合结构静平衡方程, 用结构影响系数法, 对中等展弦比、大展弦比后掠机翼的真实弹性变形和真实载荷进行了求解, 并对计算的结果进行了分析。

1 弹簧近似方法

根据Blom[5]的定义, 弹簧近似方法有两种描述方法, 一种是弹簧的平衡长度为零, 称为顶点弹簧 (vertex springs) ;另一种是弹簧的平衡长度等于边界未运动和变形前边的原长, 称为棱边弹簧 (segment springs) 。

传统的顶点弹簧近似方法一般用于网格优化, 弹簧倔强系数取为常数, 网格点移动后的新位置坐标相当于周围点坐标的平均。对这一经典方法进行了修改, 使其更适于实现网格变形运动。下面给出修改后的顶点弹簧近似方法。对于顶点弹簧, 节点i, j间的弹簧张力可以令为:

式 (1) 中Kij为连接节点i, j的弹簧的倔强系数;ri, rj分别为节点i, j的位置矢量。对任意的一个节点i, 它的合力方程为:

式 (2) 中Ni为计算区域内所有与节点i相连的非结构网格的节点总数。不同于网格优化的情况, 修改后的顶点弹簧方法令计算区域内所用网格节点的受力都为边界没有发生运动和变形时所受的合力。对计算区域内所有的m个非结构网格的节点列上述合力方程得到初始状态下的该线性系统的矩阵表达式:

式 (3) 中Ni为包括边界点的网格节点数。当边界运动或变形后, 通过改变边界的位置矢量和适当的整理式 (3) 就得到一个m阶的线性方程组。由于该线性方程组的刚度矩阵对称正定, 则求解的GaussSeidel迭代法关于任意初始向量均收敛。其相应的分量迭代格式为:

当边界变形或运动到n+1时刻, 利用式 (4) 经过数次迭代便可以达到需要的精度, 得到n+1时刻的网格点的位置。

2 弹簧倔强系数的选择

2.1 标准方法

标准的弹簧近似方法大都选取网格边的边长为参考量来确定弹簧的倔强系数, 大量数值试验表明, 对于二维或三维网格, 弹簧倔强系数按下式取能较好保持网格的疏密特征。

式 (5) 中lij为网格边的边长。于是当lij→0时, kspringij→∞。这样便很好的保证了任意的网格节点i, j不会碰撞。

2.2 方法改进

标准的弹簧近似方法仅仅考虑了弹簧在直线方向的伸缩作用, 在二维或三维时并不能保证避免网格边的互相交叉。而且, 在用标准的弹簧近似方法处理相对大一些的边界运动和变形时, 就会产生体积为负的“无效”网格单元。为了防止产生这种“无效”网格单元, 就必须在弹簧倔强系数中引入单元几何形状的变量。因此, 在实际计算中我们对公式 (5) 稍作改变, 成为公式 (6) 。

式 (6) 中Vmin为拥有网格边ij的所有单元中体积最小的单元。当Vmin→0时, kspringij→∞。这样就保证了与网格边ij相对的所有单元中体积最小的单元不会成为体积为负的“无效”网格单元。

为了防止挤压已有的四面体单元及四面体中的三角形面元, 必须考虑挤压弹簧倔强系数。挤压弹簧倔强系数计算的具体公式为:

式 (7) 中θ为网格边ij所对应的三角形内角以及网格边ij所对应的四面体单元所夹内角, 如图1所示。当θ→0或θ→π时, ksquashij→∞。这样就避免了挤压已有四面体及三角形面元。

此外, 一些文献把采用弹簧近似方法调整网格点位置的方法与椭圆型结构网格生成方法比较认为, 此种方法也具有椭圆型方法的性质, 即局部的扰动只在局部产生影响。计算也证明了这一点, 最早出现网格变形失败的单元就是在运动边界附近。为减小这一性质的影响就必须使弹簧的倔强系数沿内边界到外边界的径向方向合理分布, 因此我们引进弹簧倔强系数的径向分布系数Kd, 称之为边界修正。

经过上述的改进和修正, 最后得到总的弹簧系数为:

式 (8) 中Kd=1/d, d为计算区域内的点到物面的最小距离。改进后的方法与标准方法的对比结果如图2所示, 图2 (a) 是标准方法的结果, NACA0012翼型绕1/4弦长点偏转35°, 尾部已经出现失败。图2 (b) 是对同一问题使用改进后的方法, 此时翼型偏转到50°, 网格质量仍然较好。

3 气动力数值计算方法

3.1 主控方程

采用的控制方程是ALE (arbitrary lagrangianEulerian) 描述下的三维可压缩积分形式的Euler方程, 其在直角坐标系下积分形式为:

Ω是控制体, 是控制体边界, n为控制体边界外法向单位向量, 守恒项和对流项分别为:

式 (10) 中U={u, v, w}, ρ、u、v、w、e分别为流体的密度, x、y和z轴方向的流体速度和单位体积流体的总内能, nx、ny、nz是n的3个分量。

3.2 数值求解方法

采用中心有限体积法进行空间离散[6], 引入人工黏性后的离散表达式为:

式 (11) 中Vi是网格单元体积, Qi是网格单元内各物理量平均值, Ci为网格单元通量项, Di为人工黏性项。引入当地时间步长和隐式残值光顺的加速收敛措施后, 采用四步龙格-库塔方法对式 (3) 进行时间推进求解[7]。计算的物理空间生成非结构四面体网格, 采用弹簧近似的动网格方法耦合机翼弹性变形[8,9]。

4 机翼弹性变形的求解方法

采用结构影响系数法[10]求解机翼结构系统弹性变形。首先, 通过结构有限元分析方法得到机翼的结构柔度影响系数。

在小变形的假设下, 机翼表面结构上任一点 (x, y) 处的变形可写为:

式中W (x, y) 为 (x, y) 点处的法向变形, CZZ (x, y, ξ, η) 为结构柔度影响系数, F (ξ, η) 为作用于 (ξ, η) 处的载荷, 主要可分为气动力和惯性力。

式中L (ξ, η) 为作用于 (ξ, η) 处的气动力, N是过载系数, m (ξ, η) 是 (ξ, η) 的质量, g为重力加速度。

写成矩阵为:

机翼表面的坐标加上相应位置的变形, 就得到新的表面坐标。

静气动弹性变形的求解是一个迭代过程, 通过求解Euler方程得到气动力, 再耦合结构静平衡方程计算变形, 计算变形后再重新生成网格计算气动力, 然后再计算变形, 如此循环计算直到变形和气动力俊收敛, 则得到飞行器最后的弹性平衡状态。

5 算例及结果分析

5.1 中等展弦比后掠机翼静气动弹性变形计算

以一中等展弦比后掠机翼为例, 进行了刚性时气动力、弹性时气动力及静气动变形的计算, 该机翼半展长为7.41 m, 展弦比3.80, 梢根比0.561 8, 前缘后掠角29.91°, 刚轴位于机翼展向横截面27%弦长处。计算状态为:Ma=0.84, α=3.06°。

5.2 大展弦比后掠机翼静气动弹性变形计算

本文所计算的某无人机大展弦比后掠机翼的展弦比为8, 半展长为11.348 6 m, 前缘后掠角为27.66°, 刚轴位于机翼横截面37.27%弦长处。计算的状态为:Ma=0.80, α0=3.00°。

5.3 结果分析

从图3本文计算的机翼刚性状态下的展向升力系数与文献[11]上风洞试验结果和全速势方程计算结果比较可以看出, 本文方法计算结果好于全速势方程计算结果, 更接近于风洞试验结果, 表明本文发展的气动力计算方法是比较精确的。

弹性后掠机翼展向各剖面的扭转角分布见图6、图8, 图5、图7为后掠机翼刚轴展向弹性变形, 图4、图9为后掠刚性/弹性机翼沿展向升力系数分布, 图10为弹性变形前后的弦线平面, 虚线表示变形后的弦线平面, 可以清楚的看到, 静气动弹性变形使机翼产生向上的挠度而向上弯曲, 同时产生负的扭转角, 因而顺气流方向迎角减少了, 使得后掠弹性机翼的剖面升力系数小于刚性机翼的剖面升力系数, 从图4、图9可以清楚的看到, 而且, 静气动弹性变形产生的挠度和负的扭转角沿展向从翼根到翼梢逐渐增大, 这与图5~图8的计算结果也是比较相符的。

6 结论

本文采用了改进的弹簧近似方法, 以Euler方程为控制方程, 耦合结构静平衡方程, 用结构影响系数法, 对中等展弦比、大展弦比后掠机翼的气动力载荷及静变形进行了求解:并对后掠机翼的气动特性进行了研究, 得到以下结论: (1) 本文所发展的弹簧近似方法, 在不降低计算效率的同时, 可以大大提高网格的变形能力和网格质量, 很好地解决了非结构动网格技术中关键的网格变形问题。 (2) 后掠机翼静气动弹性变形使机翼产生向上的挠度而向上弯曲, 同时产生负的扭转角, 使弹性后掠机翼的总升力减小, 因次, 在后掠机翼设计时需要沿展向设计一个正的预转角, 这样后掠机翼在飞行中才能减小静气动弹性的影响, 提高后掠机翼的气动特性。 (3) 本文发展的静气动变形的求解方法, 可以应用到后掠弹性机翼气动/结构一体化设计中。

参考文献

[1] Flomenhoft H I.Aeroelasticity and dynamic loads-From 1903 to the supersonic era.AIAA Paper 2000—1597, 2000

[2] Manojk B, Kapania B K, Reichenbach E.A CFD/CSD interaction methodology for aircraft wings.AIAA Paper 98—4673, 1998

[3] Batina J T.Implicit flux-split Euler schemes for unsteady aerodynamic analysis involving unstructured dynamic meshes.N91—10918, 1990

[4] Batina J T.Unsteady Euler airfoil Solutions using unstructured dynamic meshes.AIAA Paper 89—0115, 1989

[5] Blom F J.Considerations on the spring analogy.Journal of Aircraft, 2000;32:647—668

[6] Delanaye M.Quadratic-reconstruction finite volume scheme for compressible flows on unstructured adaptive grids.AIAA Journal.1977;35 (4) :631—639

[7] Jameson N, Schmit W, Turkel E.Numerical solution of the Euler equations by finite volume methods using runge-kutta time stepping scheme.AIAA Paper 1981—1259, 1981

[8] 史爱明.非结构动网格下非定常气动力计算和跨音速嗡鸣研究.西安:西北工业大学, 2003Shi Aiming.Computing the unsteady aerodynamic loads and investigating transonic buzz by using the unstructured dynamic mesh.Xi’an:Northwestern Polytechnical University, 2003

[9] Gaitonde F.A three-dimensional moving mesh method for the calculation of unsteady transonic flows.The Aeronautical Journal, 1995;99 (984) :150—160

[10] 杨永年, 气动弹性教学素材.西安, 西北工业大学, 1986Yang Yongnian.Aeroelasticity teaching matetials.Xi’an:Northwestern Polytechnical University, 1986

弹性特性 篇5

膜片(碟形)弹簧广泛应用于汽车等机电产品中,一般均进行喷丸和强压等强化处理。由于强化处理对弹簧的寿命有很大影响,故许多研究者给予关注[1,2,3,4,5,6]。实际弹簧的载荷变形关系与目前教科书中给出的Almen-Laszlo(简写A-L)推导的计算公式(以下称A-L公式)计算结果有较大的差别,为了减小计算误差,一些研究者对弹簧的载荷变形关系进行了研究[1,6,7,8,9,10]。这些研究中,一种解决问题的方法是根据试验统计数据对A-L公式进行修正[1,9],另外一种方法是采用有限元法计算载荷变形关系。另外,20世纪90年代之前,还有研究者推导过其他形式的计算公式[1,10,11]。研究表明[1,5,6],这些方法都存在着各自的问题。特别是这些研究都没有考虑到喷丸和强压等强化处理对载荷变形关系的影响,以及如何在设计计算中考虑强化处理问题。本文研究了包含喷丸和强压等强化处理残余应力信息的膜片(碟形)弹簧载荷变形关系设计方法。

1 膜片弹簧的试验研究

1.1 试验准备

为了讨论喷丸和强压对膜片弹簧载荷变形关系的影响,首先进行膜片弹簧样件载荷变形相关数据的测量。检测在专用的试验设备上进行,试验检测夹具按照弹性样件设计要求设计制作。测试过程中,测试系统自动记录弹性元件的载荷与变形。选择常用结构型式的3种规格的膜片弹簧为试件(A类试件、B类试件、C类试件),图1为试件的结构示意图,图中符号的含义与文献[1]相同。3种规格膜片弹簧试件的区别是:B类试件的外径较大,C类试件的外径较小;D型孔形式不同(A类试件为长圆形的D型孔,B类试件和C类试件为T形的D型孔);膜片弹簧的锥角大小也有较大差别。

1.2 载荷变形特性试验结果

图2所示为试验样件的载荷变形特性关系曲线。图2中,未强化处理过是指试件未做喷丸或强压处理;强压处理过是指试件进行过喷丸或强压处理。计算膜片弹簧载荷变形关系的A-L公式为

$\begin{array}{l}F{}_1=\frac{\text{π}Et\lambda {}_1\ln (R/r)}{6(1-\mu {}^2)(L-l){}^2}[(h-\\ \lambda {}_1\frac{R-r}{L-l})(h-0.5\lambda {}_1\frac{R-r}{L-l})+t{}^2](1)\end{array}$

式中,R、r分别为碟形部分的外半径和内半径;L、l分别为外支撑半径和内支撑半径;F1、λ1分别为载荷和变形;t为板厚;E、μ分别为弹性模量和泊松比。

从图2可以看出,强化作用对试验样件的载荷变形特性有明显的影响,采用A-L公式计算的载荷峰值误差与谷值误差中的较大值超过了20%。根据经验,膜片弹簧载荷变形测试中可以出现更大的峰值载荷误差与谷值载荷误差。由于A、B、C三类试件是任意抽取的,故强化作用对载荷变形特性的影响是具有一定普遍性的。

1.3 残余应力检测

由于试件强化作用的结果是以残余应力的形式体现的,故笔者进行了残余应力检测。目前测量金属材料构件残余应力的方法有多种[12],技术上成熟且实用的是X射线衍射法,故笔者采用X-350型X射线应力测定仪检测残余应力。

膜片弹簧的结构是轴对称的,采用极坐标分析其应力、应变更加方便。所以,检测残余应力时,也按照极坐标法,即检测径向残余应力、切向残余应力,以及它们沿着深度方向的分布。图3为试件碟形部分的子午截面及检测位置示意图,在弹簧的上下两个表面(凸面、凹面)沿径向布置3个测量点,在3个表面测量点位置再沿板厚方向布置2个测量点,检测结果如图4、图5所示。

2 强化处理影响载荷变形关系的解决方法

2.1 影响机理分析

如上所述,强化处理后,膜片弹簧载荷变形关系中采用A-L公式计算的载荷与实际载荷有大的误差。A-L公式中的几何结构参数、材料特性等参数,若因工艺因素而改变,则将影响计算所得载荷的正确性。如果喷丸与强压使这些参数改变,则A-L公式计算结果将出现误差或很大的误差。

喷丸和强压强化结果的一个表现形式是膜片弹簧的几何结构参数发生变化,改变的结构参数主要有内径、外径及锥角。但是,几何结构参数对计算载荷变形关系的影响可以消除。因为用A-L公式设计计算时,已经充分考虑了几何结构参数。几何结构参数唯一可能带来的影响是几何结构的不规则性,这一点文献[6]已经做了充分的讨论。对于制作比较好的膜片(碟形)弹簧,几何结构参数产生的误差很小。在选择样件的时候,不选择几何结构不规则程度比较严重的样件,以排除几何结构参数对载荷变形关系的影响。

喷丸和强压强化结果的另一个表现形式是残余应力。根据图4、图5所示的结果,试验得到的切向残余应力与文献[2]的形式及数量级一致,径向残余应力与切向残余应力是同数量级的数值。然而,文献[1,2,3,4,12,13,14]仅讨论了喷丸和强压对疲劳寿命的影响,而几乎没有讨论过喷丸、强压对载荷变形特性的影响。实际上弹簧起作用时,残余应力自然也将成为总应力的一部分,这将使材料抵抗变形的能力发生变化。另外,笔者用有限元法试算了弹簧的应力,其中的径向应力不大于100MPa。因此,弹簧使用过程中的径向应力比径向残余应力小很多。总之,在喷丸与强压后,表层材料的特性发生了改变(主要是材料的应力应变关系发生了变化),这使试件抵抗变形能力改变,从而导致A-L公式计算结果产生了大的误差。

2.2 解决方法

有几种方法可以解决喷丸与强压引起的弹簧载荷变形关系设计计算误差的问题。

(1)建立准确的材料模型。

这种方法原则上是可以解决问题的,但需要进行严密的理论推导和大量的模型试验研究,目前,该方法实现起来还有诸多困难。

(2)采用修正系数方法。

包含残余应力的材料的数学模型难以推导,但可以采用修正系数的方法校正由喷丸与强压处理带来的计算误差。笔者根据图2中的试验数据,提出在A-L公式中加入修正系数f来修正计算的方法,即在式(1)的右侧乘以修正系数f。修正系数的具体表达形式为[7]

f=a+k(h-λ1) (2)

其中,a、k为求解修正系数的拟合系数,与喷丸、强压工艺参数有关。

在A-L公式中增加修正系数后,对图2中A类、B类、C类试件载荷变形关系的峰值和谷值进行计算,试件误差中有一个为5.8%,其余的在5%以内。该方法的推广使用,需要更多的试验数据支持。

(3)采用多种材料模型组合的有限元法计算弹簧载荷变形关系。

这种方法解决问题的基本思路是:首先将膜片或碟形弹簧实体分为多个区域,具体做法是沿着板厚方向分为多层,同时沿着径向也分为多层;而后,对不同区域的材料使用不同的材料模型在此基础上,将研究对象离散化为实体单元,并施加约束及位移加载;最后,进行计算。由于不同部位的残余应力情况不同,若要真实地反映实际情况,需要将弹簧实体划分为无限多的区域,这在实际中无法实现。但可以将弹簧划分为较多的区域来近似代替无限多区域,同样可以获得比较理想的结果。

根据以上思路,笔者对图2所示的A、B、C三类试件进行了试算。考虑膜片弹簧的结构对称性,取弹簧1/(2n)(n为分离指个数)的扇形区域为研究对象。考虑建模处理的简便性,将模型划分为较少的区域。将研究对象沿着板厚方向均匀分为3层,凸面层和凹面层沿着径向分别划分为m个小的区域,则研究对象共计被分为2m+1个小的区域。由于中间层受喷丸与强压的影响小,中间层小区域的材料使用试件材料本身的参数。

根据图4、图5所示残余应力的检测结果,给出不同小区域的材料模型参数。喷丸与强压处理后,材料在各个区域的残余应力的形式及大小有较大的差别。由图4、图5所示的检测结果可知,残余应力基本上都是压应力。考虑弹簧变形过程中,材料若为单向应力,则总应力为外载产生的应力与残余应力的合成,其应力与应变关系为

σ=σ0+E ε=E′ε (3)

式中,σ0为单向残余应力;E为材料的弹性模量;σ、ε分别为应力与相应的应变;E′为考虑了残余应力影响后的当量弹性模量,与残余应力有关。

式(3)中的残余应力为负值,则应变不是很大情况下的总应力均比较小,可以认为是当量弹性模量较小。在有限元法计算过程中,可以根据残余应力试验结果设置不同小区域的材料模型的弹性模量,即按照本小节的分析思路处理弹性模量。膜片弹簧某部位的当量弹性模量为

E′=E+aE+bE2+cE3 (4)

其中,a、b、c为拟合系数,根据实测的有残余应力、无残余应力状态下的膜片弹簧载荷变形关系,和残余应力的数值以及分布确定。

采用多种材料模型组合的有限元法对图2所示实例进行计算。根据膜片弹簧的对称性,对A、B、C三类试件分别取1/36、1/30、1/24的扇形区域进行分析计算。将膜片弹簧沿着板厚度方向划分为3个区域,即凸面层、中间层、凹面层。将中间层的弹性模量设为200GPa。再将凸面层和凹面层沿着径向分别划分为4个小区域,根据式(4)确定各个区域的弹性模量。则计算对象被划分为9个小区域,有限元建模时,共使用9种材料模型。采用三维实体单元对研究对象进行单元划分:A类试件被划分为13 785个单元、26 883个节点,B类试件被划分为7503个单元、18 409个节点,C类试件被划分为11 630个单元、27 260个节点。根据力学关系施加约束,图6为B类试件的分区与有限元模型。在有限元建模时,对模型进行了简化处理:去掉了小的圆弧和对膜片弹簧载荷变形关系影响微小的分离指部分[5]。

计算得到膜片弹簧载荷的峰值和谷值如表1所示。从表1可以看出,采用多种材料模型组合的有限元法计算峰值误差和谷值误差的绝对值中,有一项为5.8%,其余的均在5%以内,比A-L公式法计算精度高。C类试件的有限元法计算谷值误差略大,原因之一是在用试验数据拟合求解当量弹性模量时,残余应力使用的是平均值,各个试件的实际强化处理会有一定的偏差。

注:计算误差时以试验值为标准。

需要说明的是,喷丸与强压后的残余应力在弹簧的不同区域有较大的差异,若要更加接近实际情况,需要将膜片弹簧实体划分为更多的小区域,并且使用更多种类的材料模型。

3 结论

(1)喷丸与强压强化处理后,膜片(碟形)弹簧载荷变形的弹性特性发生了改变,弹簧在相同变形量下的载荷变小,而且载荷变小的数值是非线性的。

(2)喷丸与强压强化处理后,膜片(碟形)弹簧的材料产生了残余应力。残余应力将与正常使用时的应力共同起作用,导致弹簧在同样变形量的情况下的载荷发生变化,即残余应力改变了弹簧抵抗外界载荷变形的能力。

(3)目前可以采用两种方法解决由于喷丸与强压引起的膜片(碟形)弹簧载荷变形关系计算误差问题:一是在A-L公式中增加修正系数的方法;二是采用多种材料模型的有限元法计算载荷变形关系。

摘要：研究了强化处理后的膜片(碟形)弹簧载荷变形关系的计算方法。首先,进行样件载荷变形关系的测试,而后,采用X射线应力仪检测样件的残余应力数值及分布规律。在此基础上,分析喷丸与强压处理对弹簧载荷变形关系的影响机理,发现载荷变形关系与喷丸和强压有关。最后,基于试验数据拟合,得到了基于Almen-Laszlo公式的修正公式和采用多种材料模型有限元法求解膜片(碟形)弹簧载荷变形关系的方法。

双刚度弹性驱动器力学特性频域分析 篇6

对步行机器人的研究是机器人研究领域中一个活跃而又重要的分支。现在, 大多数步行机器人基本上处于实验室阶段, 特别是其在行走的速度、稳定性及地面适应能力方面仍不是十分理想[1]。与工业机械臂不同, 步行机器人在实际步行过程中, 摆动腿在落地时, 足部对地面瞬间会产生较大冲击, 造成了各关节的剧烈振动, 影响了机器人行走的稳定性[2]。特别是在快速运动过程当中, 极易对机体造成损害, 以致会造成机载设备的损坏。有时不得不以牺牲运动速度来换取机器人行走的稳定与安全。一方面, 步行过程中所产生的冲击能量通过机器人结构传递到各关节驱动器之中;另一方面, 关节驱动器作为动力源为各个关节提供动力。因此, 驱动器工作条件更为恶劣, 更容易损坏。如何减轻机器人在高速动态行走时的对地冲击, 避免关节和驱动器的损坏是机器人研究的一个重点。

步行机器人的步行技术实际上是一种仿生技术, 现在步行机器人的步态、姿态的参数以及结构等多是直接从对人或动物的研究中获得的[3]。在长期自然进化中, 人和其他动物适应各种复杂的环境和地貌, 其行走肌体具有良好的抗冲击性和对振动进行缓冲的能力。在这些生物肢体中, 软组织 (如肌、腱、软骨) 是加强刚度的组织。当负荷和变形增加时, 相应的刚度也会随着增加[4]。人类肌肉的最大刚度和肌肉放松状态下的刚度比约为5[5]。生物肌肉的这种变刚度特性使得由肌肉组成的关节运动具有很好的柔顺性和适应性。生物肌肉的刚度可以随着张力的不同而变化, 这一点是目前大多数驱动器所做不到的。

本文基于生物肌肉的变刚度这一特点, 在文献[6,7,8]的单刚度弹性驱动器的设计基础上, 针对机器人连续高速动态稳定步行的要求, 提出了一种双刚度弹性驱动器, 它可以更好地模拟生物体肌肉的变刚度特性。

1 模型的提出与假设

如图1、图2所示, 与传统的单刚度弹性驱动器相比, 改进后的双刚度弹性驱动器的结构特点是在动力源与负载之间装有带间隙的分段线性的非线性弹簧。这样, 驱动器输出力经过非线性弹簧的缓冲带动负载, 从而在驱动过程中根据负载大小可以被动调节刚度。

为降低公式的复杂程度, 以便集中精力研究驱动器的性质本身, 假定如下:

(1) 在机器人足部触地发生碰撞瞬间, 由于是承受冲击, 其带宽很大, 驱动的电机或液压等动力源由于受启动与制动过程以及功率密度等限制, 难以在瞬间完成相应的调整, 所以, 反馈控制系统根本来不及作用, 本文分析时不考虑反馈系统影响。

(2) 机器人在步行过程中, 步行足分为着地和腾空2个阶段。从动力学角度看, 由于机器人步行足着地时间短暂, 几乎是一瞬间, 所以在分析地面对机器人腿部的作用力时, 可用一个单自由度弹簧质量系统来描述, 即一个质量块m以初速度v0飞过来 (图3、图4) , 要用两驱动器分别对质量块进行制动, 则瞬态冲击功率P1由弹簧压缩速度v和质量块所受弹簧反作用力F1决定, 即

P1=F1v (1)

2 驱动器动力学分析

2.1单刚度弹性驱动器

如图3所示, 质量块m以初速度v0碰撞到驱动器, ks、x分别为该系统的弹簧刚度系数和弹簧压缩量。由动量原理[9]建立系统的运动微分方程为

质量块与驱动器接触瞬时, 设时间t=0, 有初始位移x (0) =0与初速度$\dot{x}(0)=v{}_0$, 解微分方程得

式中, ω1为该系统的固有频率, $\omega {}_1=\sqrt{k{}_{\text{s}}/m}$。

对式 (3) 求导, 得到弹簧压缩速度

对式 (3) 进行初始条件为0的拉氏变换得

另外, 质量块与驱动器碰撞过程中, 质量块所受弹簧反作用力为

F1= ksx (6)

所以, 由式 (1) 、式 (3) 、式 (4) 和式 (6) 得

对式 (7) 进行初始条件为0的拉氏变换得

2.2 双刚度弹性驱动器

如图4所示, k1和k2分别为该系统的弹性刚度系数, e为模型中质量块m与弹簧k2之间的距离。同样, 由动量原理建立系统的运动微分方程为

F (x) 为分段线性非线性弹性力 (图5) , 其表达式为

用等价线性化方法[10,11,12]求非线性方程的解, 首先建立一个与非线性振动方程相对应的等价线性化振动方程, 即

式中, me为等价质量;ke为等价弹簧刚度。

在一次近似的情况下, 方程的近似解为

式中, A为振幅;ωe为等效固有频率。

非线性弹性力在一次近似情况下可改写为

式中, φe为间隙e所对应的相位角, $\phi {}_{\text{e}}=\arcsin \frac{3}{A}$。

该系统的等价弹簧刚度为

将x的值代入式 (13) , 并进行分段积分, 可求得

因为e/A<1, 可将arcsin (e/A) 和$\frac{e}{A}\sqrt{1-(\frac{e}{A}){}^2}$展为幂级数, 于是有

该系统的等价质量为

me=m (16)

由式 (11) 、式 (15) 和式 (16) 得等价线性化振动方程, 即

$\begin{array}{l}m\ddot{x}+k{}_{2}x+k{}_{2}x\{1-\frac{4}{\text{π}}\frac{e}{A}[1-\\ \frac{1}{6}(\frac{e}{A}){}^2-\frac{1}{40}(\frac{e}{A}){}^4]\}=0(17)\end{array}$

解式 (17) 并对其进行初始条件为0的拉氏变换, 分别得

式中, ω2为双刚度弹性驱动器固有频率, $\omega {}_2=\sqrt{k{}_{\text{e}}/m}$。

质量块与驱动器碰撞过程中, 质量块受弹簧反作用力为

F1=kex (20)

同样, 由式 (20) 得到瞬态冲击功率P1并对其进行初始条件为0的拉氏变换, 分别得

3 模型评价

式 (5) 、式 (8) 、式 (19) 和式 (22) 中有两部分, 第一部分v0主要由外部环境 (如机器人奔跑速度、步行姿态等) 决定。因此, 为了关注两型驱动器 (单刚度弹性驱动器、双刚度弹性驱动器) 的抗冲击能力及弹簧伸缩量, 我们主要研究上述各式的第二部分, 即两型驱动器本身。

由式 (8) 、式 (22) 的第二部分可以看出, 影响瞬态冲击功率P1的参数有弹簧刚度k和系统固有频率ω。对于给定的机器人, 制造完成之后其整体质量随之确定下来。那么, 影响瞬态冲击功率的参数只有弹簧刚度k, 且瞬态冲击功率随着弹簧刚度的增大而增大。式 (8) 、式 (22) 表明, 降低瞬态冲击功率的因素是较小的驱动器弹簧刚度 (较软的弹簧) 和较低的系统固有频率 (较大的质量) 。同样, 由式 (5) 、式 (19) 的第二部分可看出, 影响驱动器弹簧压缩量的因素也只有弹簧刚度k。由控制理论可知, 在低频段弹簧压缩量随其自身固有频率ω的降低而增大, 而在高频段弹簧压缩量则趋于定值。

所以, 在选择弹簧刚度时, 首先依据机器人作业需要最大输出力以及驱动器最大伸缩量的要求, 确定驱动器弹簧刚度k的下限, 然后按照瞬态冲击功率最小原则, 确定驱动器弹簧刚度k的上限。

4 实例计算与结果分析

为了便于对两型驱动器的抗冲击特性和对弹簧压缩量进行分析, 首先假定单刚度弹性驱动器弹簧刚度ks与双刚度弹性驱动器在间隙e=0时的两并联弹簧合成刚度相等, 即k1+ k2= ks。选择文献[8]中提供的系统参数:m=128kg, ks=286N/mm, 选取k1= 200N/mm, k2=86N/mm, e/A=0.5, 描绘出在不同频率下两者冲击功率的对比曲线 (图6) 。图6显示外界激励频率γ处于低频段时, 两曲线基本重合;随着γ上升, 双刚度弹性驱动器首先到达共振区 (因为ke< ks) ;激励频率γ在高频段时, 双刚度弹性驱动器受到冲击功率明显低于单刚度弹性驱动器。图7为当质量块以初速度v0=10m/s分别向两型驱动器冲击时, 它们的弹簧压缩量对比曲线。由于双刚度弹性驱动器的等价刚度系数ke小于单刚度弹性驱动器弹簧系数ks, 因此在低频段前者的弹簧压缩量大于后者的弹簧压缩量;但随着外界激励频率γ的增大, 幅值趋于一致。

1.双刚度弹性驱动器 2.单刚度弹性驱动器

1.双刚度弹性驱动器 2.单刚度弹性驱动器

从图6、图7中看出, 对于双刚度弹性驱动器来说, 在结构上采用有间隙的分段线性的非线性弹簧后, 当外部冲击力幅值变化很小时, 其抵抗冲击的能力与单刚度弹性驱动器相当, 而伸缩量较大;当外部冲击力幅值变化很大时, 抗冲击能力显著增强, 驱动器的伸缩量又不会过大。所以, 与单刚度弹性驱动器相对比, 双刚度弹性驱动器提高了在高速动态下对外部环境的适应能力。

在m、k1、k2和v0不变的情况下, 双刚度弹性驱动器的非线性弹簧间隙e与瞬态冲击功率P1的关系曲线如图8所示。随着间隙e的增加, 非线性弹簧的等价刚度下降, 受到的冲击功率P1也显著下降。而当质量块与驱动器接触后某一时刻 (t= 0.5s) 时, 弹簧压缩量随着间隙e的增大而增大, 如图9所示。图8、图9说明, 随着间隙e的增大, 驱动器受到的冲击主要由系统弹簧k1承受, 冲击力就不会太大, 而相应的驱动器伸缩量则增大。

在m、k1和v0不变的情况下, 双刚度弹性驱动器的弹簧系数k2与冲击功率P1之间呈现线性特征, 瞬态冲击功率P1随着k2的增大而增大, 如图10所示;与此同时, 在质量块与驱动器接触后某一时刻 (t=0.5s) , 弹簧压缩量随着k2增大而减小, 如图11所示。图10、图11表明, 随着驱动器弹簧k1和k2之间的刚度系数差值增大, 受到外界冲击之后, 弹簧k1很快被压缩, 之后就由弹簧k1和k2共同承受冲击, 这样驱动器的伸缩量就不会太大。

从上述分析看出, 双刚度弹性驱动器的非线性弹簧组参数k1、k2和e影响着驱动器向外输出力和抵抗冲击的能力。那么, 对这三个参数的选取就是求设计变量向量z=[k1k2e]T为坐标所组成的设计空间, 使得在一定约束条件下目标函数, 即瞬态冲击功率P1最小。

具体的约束条件需根据实际要求设定, 如有时需限定一定外界冲击功率下的最大压缩量, 有时需限定两弹簧的刚度比值等。这样, 对k1、k2和e等参数的选取就转化为多变量优化问题, 最终通过优化设计得到驱动器弹簧各参数。

5 结论

本文从仿生学的角度出发, 提出双刚度弹性驱动器的结构。建立了单刚度和双刚度弹性驱动器的力学模型, 在频域内对两者的抗冲击功率性能和驱动器伸缩量进行了分析。通过分析, 本文给出了频域内单双刚度弹性驱动器冲击功率和压缩量的对比图, 同时通过定量分析给出了双刚度驱动器间隙e和刚度k1、k2等参数与压缩量、冲击功率之间的关系曲线, 这为实际驱动器的设计提供了一种依据。最后, 本文给出了双刚度弹性驱动器主要设计参数的概括性选取原则与方法。

参考文献

[1]谢涛, 徐建峰, 张永学, 等.仿人机器人的研究历史、现状及展望[J].机器人, 2002, 24 (4) :367-374.

[2]沈继红, 丁二华.一种类人型机器人的步行稳定性研究[J].哈尔滨工程大学学报, 2004, 25 (4) :536-539.

[3]马建旭, 马培荪, 杨保忠, 等.四足步行机器人中一种新型腿结构缓冲特性[J].上海交通大学学报, 1999, 33 (7) :847-850.

[4]余联庆, 吴昌林, 马世平.基于仿生研究的步行机缓冲型腿机构设计[J].华中科技大学学报 (自然科学版) , 2005, 33 (6) :105-107.

[5]隋立明, 王祖温, 包钢.气动肌肉与生物肌肉的力学特性对比研究[J].机床与液压, 2004 (6) :22-24.

[6]Robinson D W.Design and Analysisof Series Elasticityin Closed-loop Actuator Force Control[D].Massa-chusetts:MIT, 2000.

[7]Pratt J E, Kupp B T.Series Elastic Actuators for Leg-ged Robots[J].Proceedings of SPIE-the InternationalSociety for Optical Engineering, 2004, 5422:135-144.

[8]Robinson D W, Pratt J E, Paluska DJ.Series ElasticActuator Development for a Biomi metic Walking Robot[J].IEEE/ASME International Conference on Ad-vanced Intelligent Mechatronics, 1999 (6) :561-568.

[9]冯登泰.应用非线性振动力学[M].北京:中国铁道出版社, 1982.

[10]闻邦椿, 李以农, 韩清凯.非线性振动理论中的解析方法及工程应用[M].沈阳:东北大学出版社, 2001.

[11]马利娥, 葛文杰, 黄则兵.腿型仿生跳跃机器人运动机理的研究[J].机器人技术与应用, 2004 (6) :28-31.

弹性特性 篇7

自从Kneller 等[1]提出采用软磁/硬磁(H/S)交换耦合的纳米双相交换弹性磁体(Exchange-spring magnets)体系可得到120 MGOe(≈955 kJ/m3)(大约是普通永磁体磁能积的3倍)的磁能积[2]以来,由软磁/硬磁组成的双层、多层铁磁或反铁磁结构体系已成为磁性与磁性材料领域的前沿课题之一[3,4]。

近年来,由于这种双层、多层结构体系比纳米混合物更易于控制结构参数以及其磁特性的可裁剪性,在实验上双层、多层膜体系更易于实现理论所预言的超高磁能积,这种结构已成为实现“理想”各向异性交换弹性磁体的最佳方式[5]。已有多个研究团队对纳米双层或多层膜体系的反磁化机理、矫顽力性质进行了研究[6,7,8,9,10]。研究表明,在铁磁或反铁磁交换耦合软/硬磁双层、多层膜体系中,存在一个软磁层的临界厚度tc[11,12],当软磁层厚度小于tc时,软磁层严格耦合于硬磁层,软、硬磁层在同一形核场反转,形成一矩形的磁滞回线;软磁层厚度大于tc的交换弹性体系,反磁化过程通过两步完成,当外磁场达到软磁层的形核场时,软磁层首先发生偏转,由于界面处磁矩被来自于硬磁层的交换耦合所钉轧,在软磁层中沿着膜厚度方向形成螺旋式的磁结构,这个过程即所谓的交换弹性反磁化过程(Exchange-spring process),当外磁场进一步增大到硬磁层的反转场时,硬磁层发生不可逆的磁反转。人们普遍认为软磁层中发生的交换弹性反磁化过程是一个可逆过程,且在强交换耦合作用下,软磁层的形核场Hb与软磁层厚度的平方成反比[13,14]。

本课题组基于一维原子链模型研究了软磁/硬磁反铁磁交换耦合双层结构体系的交换弹性反磁化过程。结果表明,软磁层存在一个临界厚度tc,软磁层厚度小于tc时,其反磁化过程呈现典型的可逆交换弹性反磁化过程,而软磁层厚度大于tc时,软磁层中的交换弹性反磁化过程转变为不可逆过程。

1 理论分析及数值模拟

对于纳米尺度软/硬磁反铁磁交换耦合双层膜结构,薄膜厚度远小于薄膜大小,薄膜面内的退磁场可以忽略不计(可认为薄膜平面无穷大),绝大部分磁矩的运动将被限制在磁层面内。因此可以假设磁矩在面内转动,无垂直于膜面的分量,磁矩方向平行于膜面且沿薄膜厚度方向一维连续变化,θ为磁矩与易磁化轴间的夹角。软磁层、硬磁层厚度分别为ts和th,软硬磁层沿z轴增长,并且软硬磁层具有相同的易磁化方向,都沿+x方向,界面处为坐标轴z轴原点中心。外磁场方向与易磁化方向相同。根据微磁学理论,总自由能为:

$\begin{array}{l}E{}_{\text{t}\text{o}\text{t}}={\displaystyle {\int }_0^{t{}_s}\left[A{}_{\text{s}}(\frac{\text{d}\theta }{\text{d}z}){}^2+k{}_{\text{s}}\sin {}^2\theta -{\rm H}{\rm M}{}_{\text{s}}\cos \theta \right]}\text{d}z+\\ {\displaystyle {\int }_{-t{}_{\text{h}}}^0\left[A{}_{\text{h}}(\frac{\text{d}\theta }{\text{d}z}){}^2+k{}_{\text{h}}\sin {}^2\theta -{\rm H}{\rm M}{}_{\text{h}}\cos \theta \right]}\text{d}z+\\ A{}_{\text{s}\text{h}}\cos (\theta {}_{0\text{s}}-\theta {}_{0\text{h}})(1)\end{array}$

其中Ms(Mh),As (Ah) 和ks(kh)分别是软磁层(硬磁层)的饱和磁化强度、交换常数和磁各向异性常数,H为沿易磁化方向的外场。Ash是在软硬磁层相邻界面之间的交换耦合常数Ash>0。θ0s与θ0h分别为软磁和硬磁层在界面处磁化矢量与易磁化轴间的夹角。根据能量最小原理,平衡态时总自由能的一阶变分为零δE=0。因此平衡条件为:当0 ≤z ≤ts时,-2As(d2θ/dz2)+2kssinθcosθ+HMssinθ=0;而当-th≤z≤0时,-2Ah(d2θ/dz2)+2khsinθcosθ+HMhsinθ=0。由于软/硬磁层间强烈的反铁磁耦合,在软磁层边界z=ts与硬磁层的边界z=-th处, dθ/dz=0,而在软磁层表面,dθ/dz=-Ashsin(θ0s-θ0h)/2As, 硬磁层表面dθ/dz=-Ashsin(θ0s-θ0h)/2Ah。因此,磁层厚度与磁矩变化角度θ将存在如下积分关系。

$t{}_{\text{s}}={\displaystyle {\int }_{\theta {}_{0\text{s}}}^{\theta {}_{\text{s}}}\left(\frac{A{}_{\text{s}}}{{\rm K}{}_{\text{s}}\sin {}^2\theta -{\rm K}{}_{\text{s}}\sin {}^2\theta {}_{\text{s}}-{\rm H}{\rm M}{}_{\text{s}}\cos \theta +{\rm H}{\rm M}{}_{\text{s}}\cos \theta {}_{\text{s}}}\right)}{}^{1/2}\text{d}\theta (2)$

$-t{}_{\text{h}}={\displaystyle {\int }_{\theta {}_{0\text{h}}}^{\theta {}_{\text{h}}}\left(\frac{A{}_{\text{h}}}{{\rm K}{}_{\text{h}}\sin {}^2\theta -{\rm K}{}_{\text{h}}\sin {}^2\theta {}_{\text{h}}-{\rm H}{\rm M}{}_{\text{h}}\cos \theta +{\rm H}{\rm M}{}_{\text{h}}\cos \theta {}_{\text{h}}}\right)}{}^{1/2}\text{d}\theta (3)$

式(2)、式(3)能很好地描述磁矩变化与磁层厚度间的关系,却很难得到Hb和θ的解析表达式。考虑软磁层界面被硬磁层完全钉轧,且ks=0的特殊情况,则θ0s =π,θ0h=0,解上述方程可得θ=2arccos[cos(θs/2)sn[-(zπ/2ts)·(H/Hb)1/2],其中sn[-(zπ/2ts)·(H/Hb)1/2]为Jacobi椭圆正弦函数,其模数为cos(θs/2),θs为软磁层表面磁矩与易轴间的夹角。因此,当H>Hb时,磁矩才开始发生偏转,在越靠近软磁层表面处, 磁矩偏转越大。Hb即是形核场,Hb=Asπ2/(2Msts2), 与软磁层的厚度成平方反比关系,这与Guslienko等[15]的研究一致。

对于ks≠0情况,软/硬磁双层结构的反磁化特性很难进行解析描述,但可通过数值模拟进行诠释。为此,基于一维原子链模型对这种结构的反磁化过程进行了数值模拟。一维原子链模型是假设磁性层由一系列的原子层沿一维方向叠加而成,每一原子层内的磁化矢量是一致的,可用一个磁矩来表示。这一模型已被广泛地用于研究磁性多层膜体系,且取得了与实验符合很好的结果[10,12,14]。根据一维原子链模型假设,体系总能量可表示为:

$\begin{array}{l}E{}_{\text{t}\text{o}\text{t}}=-{\displaystyle \sum _{i=1}^{{\rm N}{}_{\text{s}}-1}\frac{A{}_{\text{s}}}{d{}^2}}\overrightarrow{m}{}_i\cdot \overrightarrow{m}{}_{i+1}-{\displaystyle \sum _{i=1}^{{\rm N}{}_{\text{s}}}k}{}_{\text{s}}(\overrightarrow{m}{}_i\cdot e{}_i){}^2-\\ {\displaystyle \sum _{i=1}^{{\rm N}{}_{\text{s}}}{\rm M}}{}_{\text{s}}\overrightarrow{{\rm H}}\cdot \overrightarrow{m}{}_i-\frac{A{}_{\text{s}\text{h}}}{d{}^2}\overrightarrow{m}{}_{{\rm N}{}_{\text{s}}}\cdot \overrightarrow{m}{}_{{\rm N}{}_{\text{s}}+1}-{\displaystyle \sum _{i={\rm N}{}_{\text{s}}+1}^{{\rm N}-1}\frac{A{}_{\text{h}}}{d{}^2}}\overrightarrow{m}{}_i\cdot \overrightarrow{m}{}_{i+1}-\\ {\displaystyle \sum _{i={\rm N}{}_{\text{s}}+1}^{{\rm N}}k}{}_{\text{h}}(\overrightarrow{m}{}_i\cdot e{}_i){}^2-{\displaystyle \sum _{i={\rm N}{}_{\text{s}}+1}^{{\rm N}}{\rm M}}{}_{\text{h}}\overrightarrow{{\rm H}}\cdot \overrightarrow{m}{}_i(4)\end{array}$

式中:mi为第i个原子磁化强度与饱和磁化强度的比值,ei为易磁化轴方向单位矢量,d是相邻两个原子层之间的距离,Ns和N分别是软磁层和体系的原子层数。根据能量最小原理,可求得任意外场所对应的平衡态磁矩分布,以及所需的退磁曲线。数值计算参数为:Ms=557 emu/cm3、As=1.46×106 erg/cm、ks=1×106 erg/cm3、Mh=1085 emu/cm3、Ah=1.46×106 erg/cm、kh=1×108 erg/cm3、d=2.5×10-8 cm(这些参数值都是来自于YFe2和DyFe2实验参数[16])。Ash=1.46×106 erg/cm。在计算过程中硬磁层的厚度固定为30个原子层;外场方向沿易磁化方向,最大偏离易磁化轴为1°。双层结构的初始磁化状态是软、硬磁磁矩严格反平行排列,而各磁层内磁矩为一致排列。

2 结果与讨论

利用一维原子链模型首先模拟了软磁层厚度(Ns)分别为25和 60个原子层的软磁/硬磁双层结构体系的反磁化过程。当软磁层厚度较薄(Ns=25)时,体系退磁过程是典型的可逆交换弹性反磁化过程。当外磁场反向增大到Hb1时,磁化强度均匀分布的初始磁化开始不稳定,随后各原子层的磁矩开始发生偏转,由于界面反铁磁交换耦合能和塞曼能间的相互竞争,导致软磁层界面处的磁矩处偏转小,而表面原子层磁矩角度偏转较大,软磁原子层内磁矩连续偏离易磁化轴方向,而形成螺旋旋转结构或称为布洛赫畴壁(如图1(a)所示)。随着外场增大,磁矩的偏转角度也逐渐增大,畴壁逐渐向硬磁层方向压缩。图2(a)为软磁层厚度为25个原子层体系相应的退磁曲线,可以看出在Hb1处,磁矩开始发生偏转,总的磁化强度呈现连续下降趋势,最后达到极值。当外场逐渐减小返回零时,系统的磁化矢量沿着原路径返回初始状态,是一个完全可逆的过程。图3(a)给出了软磁层界面第1个原子和表面第Ns个原子磁矩与易磁化方向的夹角θ1和θNs随外场的变化曲线,可以看出,外场超过形核场Hb1后,随外场的增大偏转角逐渐增大,表明螺旋式畴壁结构的逐渐形成和压缩过程。这些都充分证明,在软磁场较薄时,体系的这种交换弹性反磁化过程是一可逆过程。图1(b)展示了软磁层厚度为Ns=60体系内各原子层磁矩在Hb1处与易磁化轴间的夹角。在软磁层界面附近少数原子磁矩由于来自于硬磁层的反铁磁耦合的钉轧作用,磁矩仅偏转一个非常小的角度,而别的大多数的原子磁矩发生了大的偏转,表面处的磁矩与易磁化轴间的夹角接近180°。从图2(b)相应的退磁曲线也可以明显看出在形核场Hb1处磁化强度的突然减小,表明这是个一级磁相变过程。此外这一反磁化过程是不可逆的,减小外场,Hb1后形成的螺旋结构在另一较小的偏转场Hb2处才回到初始的均匀分布磁结构状态。这种不可逆的交换弹性反磁化过程也可以从图3(b)θ1和θNs随外场的变化曲线中获得。在Hb1处,θ1和θNs突然发生大的偏转,表明Hb1处突然形成一个大角度的螺旋结构,导致出现不可逆交换弹性反磁化过程。

从能量的角度来分析,根据一维原子链模型,体系的能量主要由塞曼能Eze、交换能Eex、界面交换耦合能EAex以及磁晶各向异性能Ean决定。交换弹性反磁化过程是可逆还是不可逆过程主要取决于在形核场Hb1前后,软磁层中交换能和磁晶各向异性能的变化ΔEex+ΔEan和塞曼能变化ΔEze间的竞争。在Hb1处,软磁层内磁矩开始发生偏转,若此时塞曼能减小能弥补形成大角度畴壁所需畴壁能(交换能和磁晶各向异性能)的增加,则畴壁即可形成,交换弹性反磁化为不可逆过程;反之,反磁化是可逆过程;但通常这种临界条件很难得到解析表达式,需进行数值计算。但可确定的是存在一个利于大角度畴形成的临界软磁层厚度tc。软磁层厚度大于tc,体系更易形成一大角度畴壁,形成磁滞现象,导致不可逆交换反磁化过程发生。若软磁层厚度接近或超过180°畴壁厚度时,外场在Hb1处,一个180°畴壁立即形成,导致不可逆的交换弹性过程。因为界面处软磁层磁矩更接近与初始态反方向的易磁化方向,软磁层的磁晶各向异性将对这些磁矩形成钉轧作用,从而形成磁滞现象,产生不可逆的交换弹性反磁化过程。

为了更深入地理解随着软磁层厚度增大,可逆交换弹性反磁化过程向不可逆转变。图4给出了当Ash=1.46×106 erg/cm时,在形核场Hb1处软磁层表面和界面磁矩与易磁化轴间夹角随软磁层厚度的变化曲线。由图4可以看出,软磁层存在一个临界软磁层厚度(数值计算tc=8.75 nm),当ts<tc时,θ1(Hb1) 和θNs(Hb1)都为零,表明交换弹性反磁化过程是可逆过程,当ts>tc时,θ1(Hb1) 和θNs(Hb1) 突然增大,表明不可逆交换弹性反磁化过程出现,这一可逆过程伴随着一个大角度畴壁的突然形成。值得注意的是,软磁层的临界厚度tc远小于180°布洛赫畴壁厚度(δs=π(As/ks)1/2 =37.9 nm)。当软磁层厚度满足tc<ts<δs时,在形核场Hb1处,软磁层中立即形成一个大角畴壁,且随着软磁层厚度的增大,θ1增大得越快。当软磁层的厚度接近δs时,θ1在Hb1处立即达到180°,软磁层磁矩形成一个完整的180°布洛赫畴壁。模拟也发现,软磁层的临界厚度随着ks的减小而增大。当ks→0时,δs→∞,因此不可能发生不可逆交换弹性反磁化过程。

导致磁结构不稳定的形核场Hb是软磁/硬磁反铁磁耦合双层膜体系的反磁化特性的另一个重要参量,只有在某些特殊情况下,才可获得Hb的解析表达式。当软磁层界面磁矩完全被硬磁层所钉轧,ks=0情况下,形核场Hb=π2As/(2Mst${}_{\text{s}}^2$)。一般情况下,一个简单的经验公式Hb=Hb0/tsn(Hb0和n为与磁参量有关的常数)常常用来描述交换耦合双层膜体系的形核场随软磁层厚度的变化关系[13,14],当软磁层非常厚时,n近似等于1.75。基于一维原子链模型,在Ash=1.46×106 erg/cm时,得到了Hb1随ts的变化关系,发现这个简单的经验公式并不能完全描述Hb随ts变化的规律,存在一个临界软磁层厚度tc。图5展示了形核场Hb1 和Hb2随ts的变化曲线,图6为lnHb-lnts的关系曲线,可以看出在软磁层厚度较小时,lnHb-lnts具有非常好的线性关系,经验公式可很好地描述Hb随ts的变化,发生可逆的交换弹性反磁化过程。随着ts的增大,Hb开始一分为二,表明发生不可逆交换反磁化过程,Hb1 和Hb2 不再满足经验公式,lnHb-lnts呈非线性关系,且Hb1和Hb2随ts的变化满足不同的规律。这也进一步证明了体系的反磁化过程由可逆的交换弹性反磁化转变为不可逆过程存在一临界软磁层厚度tc。

3 结论

基于微磁学理论对反铁磁耦合软磁/硬磁双层结构的反磁化过程进行了研究,当界面交换耦合强度非常强以至于界面软磁磁矩完全被硬磁层所钉轧,ks=0的情况下,形核场Hb=π2As/(2Mst${}_{\text{s}}^2$)。利用一维原子链模型模拟了这种结构体系的反磁化过程。研究发现,反铁磁耦合软磁/硬磁双层磁结构体系软磁层存在一个临界厚度tc,当软磁层厚度小于这个临界尺寸时,软磁/硬磁双层结构体系的反磁化过程呈现典型的可逆交换弹性反磁化过程,而软磁层厚度大于tc时,交换弹性反磁化过程为不可逆过程。研究了形核场随软磁层厚度的变化关系,仅当软磁层厚度较小, 反磁化过程为可逆的交换弹性反磁化过程时, 形核场Hb可由经验公式Hb=Hb0/tsn描述, lnHb随lnts呈线性减小。

摘要：基于微磁学理论研究了软磁/硬磁反铁磁交换耦合双层结构体系的反磁化特性,利用一维原子链模型模拟了其反磁化过程。研究表明:当考虑了软磁层的磁晶各向异性能后,随着软磁层厚度的增大,交换弹性反磁化过程从可逆过程转变为不可逆过程。存在一个临界的软磁层厚度tc,当软磁层厚度ts<tc时,交换弹性反磁化过程为可逆过程;而当ts>tc时,交换弹性反磁化过程为不可逆过程。形核场Hb随软磁层厚度的变化仅当体系的反磁化过程为可逆的交换弹性反磁化过程时才满足经验公式Hb=Hb0/tsn。

弹性特性 篇8

盾构法隧道施工具有自动化程度高、工作效率高、安装性强、受地面环境限制小、较少破坏人文自然景观、对城市交通影响小等优点,是地铁建设中区间施工常用方法之一[1]。岩石由于结构、岩性、所含矿物及环境不同,其抗风化能力也有差异,在花岗岩地区,当完整岩石中有交叉节理发育时,节理把岩石分割成棱角形块,风化作用在棱角部位加速,棱角部位岩石首先风化成土,被节理包围岩块变圆,形成球状风化体(孤石),其大小、规模、分布具有随机性和隐蔽性,但主要分布于全风化带和强风化带中。孤石给工程建设带来众多问题。在城市轨道建设中,孤石严重影响盾构法隧道施工,盾构过程中如遇到孤石,孤石会使盾构机刀盘磨损严重、刀座变形,严重时刀盘由于受力不均匀导致主轴承受损,较大孤石可能导致盾构转向,偏离隧道轴线。因此,为了确保工程顺利进行,减少经济损失,在盾构机通过之前必须查明孤石分布范围,提前清除孤石障碍[2]。

目前,孤石探测技术已在国际上引起重视并得到广泛研究。研究表明,地震和电磁层析成像技术采用射线“照射”目标体的观测方式及成像计算方法能获得较高的探测分辨率[3],已经广泛应用于矿体探测和工程勘探。两种方法综合运用于孤石探测及探测效果分析的参考资料较少,为给类似工程提供应用经验,本文从方法原理、工作方法、数据处理及解释做了详细论述。

1 孤石波速和衰减系数特征

孤石赋存在花岗岩全、强风化层中,由于所处环境和风化条件不同,相对于风化残积土孤石完全保留了花岗岩的物理化学性质,岩石参数良好[4]。表1[5]是某地孤石、强风化层、全风化层物性参数。

上述数据说明孤石与周围介质存在明显的波速差异和对电磁波的吸收差异,这为孤石探测提供了良好的应用地球物理基础。理想条件下微风化体和全强风化存在很大波速差异,差异比约为3.1;其衰减系数差异更大,差异比约为7.5。实际情况下二者差异没那么大,例如:对于风化程度相当、裂隙不发育、弱含水的孤石,其电性差异较小,波速差异占主导;对于风化程度相当、裂隙发育、富含水的孤石,其波速差异不明显,其电性差异占主导。因此,在对孤石发育特征不甚了解情况下,综合运用两种方法从孤石的波速和电性参数探测其发育范围更可靠,综合分析两种资料有助于孤石边界划定。以实例探测结合基坑开挖资料,针对两种探测技术对孤石的探测效果进行客观分析和评价,可为以后类似工程提供参考依据。

2 基本理论

2.1 跨孔层析成像

跨孔层析成像又称井—井CT(Computer Tomography),是在常规物探方法基础上发展起来的新技术,它通过在探测介质一侧发射某种信号(例如电磁波、弹性波、电压等),在另外一侧利用仪器记录透射信号,根据被测介质对信号的影响特征,利用计算机层析成像技术重构介质内部结构,再分析重构的某个物理量图像与介质的对应关系达到层析成像技术应用目的[6]。该技术具有很多显著优点,首先,数据采集时,信号发射源和接收器均放入钻孔中,接收器远离地面直达目的层,采集信号为目的层直接反应,提高了数据信噪比;其次,跨孔探测反应的资料为两钻孔之间区域,相比地表勘探单位面积数据量大,为反演技术应用提供了数据基础,探测精度高于地面物探方法;最后,跨孔探测是在孔中进行,占地少,对场地条件要求不高,更适合于城市工程物探。

跨孔层析成像工作时至少需要两个钻孔,一孔作为信号源发射孔,另一孔或多孔作为信号接收孔,然后按照某种观测系统进行信号采集,目前国内外发展起来的观测系统有以下三种[7]:①共发射点观测,在发射孔固定信号发射点,在一个或多个接收孔中同时布置多个接收器进行采集,工作示意图如图1所示;②共接收点观测,在接收孔固定接收点,在发射孔不同深度依次发射物理信号;③平行观测,观测时发射点、接收点分别在发射孔和接收孔平行等间距移动,顺序由深向浅或由浅向深。

跨孔弹性波CT工作时,在激发孔采用震源发射器发射弹性波,在接收孔利用传感器接收经过介质散射后的弹性波,接收信息传输到地震仪并记录,通过分析记录到弹性波运动学(走时、射线路径)和动力学(波形、振幅、相位、频率)特征,达到地质勘探目的(见图1)。而走时层析成像主要利用弹性波中的直达波,是对拾取的直达波时间(又称初至)进行反演,达到孔间弹性波场重建目的。按成像方式可分为基于射线方程的层析技术和基于波动方程的层析技术,按重建参数可分为波速层析成像和衰减系数层析成像。基于射线理论的弹性波CT由于旅行时具有信噪比高、受激发条件影响小、信号稳定、且无论是柱面波还是球面波其走时规律均相同等优点,相对于其他层析成像方法发展较早,也较为成熟,是目前地震层析成像的主要方法(1992,zhou Bing)[8]。

跨孔电磁波CT的工作方式与跨孔弹性波CT类似,它是在一钻孔放置发射天线,另一钻孔放置接收天线。发射天线发射固定频率的电磁信号,接收天线感应经地层传播来的电磁波场,衰减后的电磁波经专门仪器记录。通过分析接收到的电磁波场强和相位变化可达到地质勘探目的。由于硬件技术限制,目前相位反演尚处于研究阶段,相比之下场强应用较为广泛。电磁波CT工作时应避免采用金属套管保护,前人研究成果表明金属套管会对电磁波传播产生较大影响[9]。

2.2 正演理论

弹性波和电磁波在介质中传播具有相同的运动学性质[10],均遵从射线理论,因此传播路径正演追踪具有相同算法。数据分析区别在于反演算法,以弹性波为发射源的层析成像技术,可以对能量、旅行时和相位反演,以电磁波为发射源的层析成像技术,也可以对以上三种物理量进行反演,但是,基于当前硬件发展现状,弹性波勘探一般对旅行时进行层析成像,电磁波勘探一般对衰减能量进行层析成像。

弹性波传播过程中遵从如下公式:

式中:i为传播路径编号;j为剖分网格编号;Ti为第i条射线的传播时间;Li为第i条射线的传播路径长度;dl为第i条射线穿过第j个网格的弧长微元;Vj(x,y)为第j网格的速度,为方便计算,有时用慢度Sj(x,y)表示,且;由(1)式知弹性波的传播时间就是慢度沿传播路径的积分。

电磁波在传输工程中,电磁波振幅(能量)的衰减遵从如下公式[11]:

式中:Ai第i条射线的衰减后能量;Li为第i条射线的传播路径长度;Lij为第i条射线在j网格里面的弧长微元;βj为第j网格的衰减系数;A0为发射天线的初始辐射常数;fS(θS)、fR(θR)分别为发射天线和接收天线与垂线之间角度的函数。为讨论方便,公式推导时忽略发射和接收天线之间的角度影响,相当于把发射场认为球形源,简化后的公式为:

两边取对数之后再进行数学变换后,式(3)变为式(4):

简写为:

式中:Ei为接收天线的接收到衰减后的电磁场强数值。

对比式(1)和式(5),从数学分析两个公式一致,只是代表物理量不同而已,因此数据经变换之后,弹性波和电磁波层析成像所归结数学问题是完全相同的。

弹性波正演是以介质的波速分布为基础,按照波的传播理论追踪传播路径,计算从发射点到接收点的传播时间,射线追踪过程即是正演计算。正演计算的精度和速度决定着成像的分辨率和可靠程度。弹性波射线追踪理论是在高频近似条件下,波场的主能量沿射线轨迹传播。追踪算法分为直线追踪和弯曲射线追踪两类[12]。本试验采用射线理论中的直射线追踪算法,计算过程如下。

首先将孔间成像区域网格化成M×N个单元,每个网格用慢度(速度的倒数)Sj(j=1,2,3…M×N)表示(见图2),网格单元内图像函数看作常量,因此,图像函数沿第i条射线的积分值即投影函数为:

式中:i为通过图像区射线号;j为图像区单元号;Lij为第i条射线通过第j个单元的路径长度;Sj为图像区函数值,对于弹性波代表慢度,对于电磁波正演代表衰减系数[13]。

对于图2模型,第i条弹性波射线在成像区域中的旅行时间由式(1)和式(6)推导为:

由此可见,成像区域内介质慢度是旅行时在射线方向的投影函数。假如有N条射线依次通过图像重建区时,式(7)可列成如下线性方程组:

基于上述正演算法,地质模型给定后,由于每个网格单元慢度以及发射点和接收点坐标已知,每个网格切割每条射线的长度Lij则可求取。再由式(7)计算出每条射线的理论旅行时Ti,最后根据式(8)建立关于每条射线通过成像区域后旅行时矩阵方程。这就是已知慢度求理论旅行时的正演计算过程。正演计算的目的是确定弹性波旅行时路径,求取理论旅行时,为反演计算做准备。

2.3 反演理论

地球物理反演是依据观测信息和前向物理模型求解或推断地质体的赋存状态和物性参数的计算过程,最终归结为求解一个大型的﹑稀疏的﹑超定的线性方程组。常用算法有:标准层析重建方法(如:代数重建ART和联合迭代重建SIRT),最小二乘法,迭代矩阵法,共轭梯度法等。在工程中以SIRT应用最为广泛[14,15]。

SIRT(Simultaneous Iterative Reconstruction Techniques)算法首先由Givert提出。该算法把第K轮所有射线对某一网格单元的修正值保存下来,直到本轮迭代结束后求所有射线对该网格的修正平均,并作为下一次迭代的初值。SIRT算法对超定或欠定方程组均能求解,且不会出现解的奇异现象,同时具有平稳收敛、计算简便快捷、占用内存少等特点。详细计算步骤如下:

(1)初始化方程组中慢度值Sj,根据式(7)计算M条射线的旅行时估算值Ti(i=1,2…M),其中Lij由射线追踪算法求取。

(2)求每条射线实测旅行时ti与估算值Ti之差Δti=ti-Ti(i=1,2…M),按照加权算法计算第i条射线在第j个网格中的旅行时差Δti分摊量,并计算第j个网格中所有射线旅行时差Δti分摊量之和∑Δti。

(3)计算每个网格单元切割每条射线的长度之和∑L,求每个网格的慢度修正量。为直观起见,上述步骤的算术表达式如下:

式中:Mj表示穿过第j个网格的射线数;N表示射线数;M表示模型离散化后的网格数。为了加快计算速度,并突出短射线的作用,其平均化算法更为实用,即:

式中:和分别表示旅行时变化量和第i条射线的总长度;Ni表示第i条射线穿过成像区域的网格数。为了让求解结果更加平稳,可对每次迭代求取的慢度修正量ΔSj进行平滑。

(4)把步骤(3)中计算的慢度修正量ΔSj(k)与Sj(k-1)相加,作为模型的下一次迭代的初值。重复上述步骤,直至满足限定误差条件为止。

3 实例分析

试验地位于广州东北部,经勘察上覆地层为人工堆填、冲积、洪积、坡积及残积的第四系,主要包括全新统(Q4)和全—上更新统(Q3+4)。下伏基岩为中生代燕山期侵入岩,属燕山晚期第一阶段的萝岗岩体,以细、中、粗粒斑状黑云母二长花岗岩为主体岩石结构。在花岗岩全、强风化层中存在大小不等且随机分布的孤石,其中一处经钻探揭露,钻孔柱状描述为:0~2.6m为人工堆积粘性土;2.6~7.0m为坡积而成的粉质粘土;7.0~9.8m为花岗岩风化残积而成的砂质粘性土;9.8~12.2m为花岗岩微风化体(孤石);12.2~16.5m为花岗岩风化残积而成的砂质粘性土;16.5~33.6m为全风化,局部夹强风化花岗岩;稳定水位2.5m。本例以此钻孔揭露孤石为探测对象,探测方法综合应用上述介绍的跨孔弹性波和电磁波CT技术。

3.1 试验方案

试验时在孤石两侧各布置一个钻孔,钻孔深度超过现存孤石深度,两钻孔平面距离孤石各5m,通过设计使得探测孤石位于成像区域中部。跨孔弹性波CT和电磁波CT工作参数为:激发点距1m,接收点距0.5m,采用共发射点观测系统,单点24道接收。原始记录结果见图3和图4。

注:激发深度10m;接收深度8.5~20m

注:发射点深度10m;接收点深度3~18m

图3跨孔弹性波CT单炮记录表明激发震源能量较强,有效波信噪比高,首波清晰,初至时间易于拾取,记录中15~18道初至有超前现象,推测为孤石异常所致;图4跨孔电磁波CT定点发射记录表明,电磁波随着距离的增加衰减量增加,曲线平滑,外界无电磁干扰,但曲线的曲率具有变化,说明该段岩石电性具有差异性。

3.2 数据反演及对比分析

数据处理采用上文介绍的正反演方法,为了提高计算速度和精度,程序设计时对数据格式和正、反演算法进行了部分优化,如图5是两种方法的处理结果。跨孔弹性波通过对拾取的初至反演得到孔间区域内波速分布图,跨孔电磁波通过对接收到的场强反演得到孔间区域内电磁波衰减系数分布图,剖面图中纵坐标为钻孔高程,横坐标为孔间距,图中数值表示波速数值和电磁波衰减系数,白色竖条为揭露孤石钻孔位置。

图5(a)是跨孔弹性波CT反演结果,剖面图中显示探测区域内介质波速大小存在差异,大部分区域为低速介质,个别点为高速体,高、低速分布形态具有规律性,高速体推断为孤石发育区,异常体的波速数值大小介于1900~2600m/s,横向距离ZK02孔3~6.5m,发育深度9~11m;图5(b)是跨孔电磁波反演结果,剖面图中显示探测区域内介质对电磁波的吸收能力存在差异,大部分区域对电磁波强吸收,个别区域弱吸收,吸收系数大小分布形态具有规律性,弱吸收异常区推断为孤石发育区,异常区衰减系数数值大小1~2.5m/s,平面距ZK01孔2.5~6.5m,发育深度9~11.5m;对比两种方法所得结果,对异常体的发育范围推测基本一致,通过与既有钻孔资料对比,垂向深度相差约0.7m,结果基本一致;图5(c)是所探孤石开挖后照片,经实地测量孤石平面大小为3m×5m,本次跨孔探测剖面平行于3m方向,因此两种层析成像方法对所探异常平面大小和深度大小基本吻合。

岩石是矿物集合体,具有一定的孔隙度,其物性跟岩石骨架(矿物颗粒)、胶结物及孔隙流体相互作用等因素有关[16],物探资料反演的数值是岩石某个物性的综合体现,因此对反演剖面进行地质解释时,异常体的推断应以数值的相对大小为准则。跨孔弹性波CT探测结果反映了异常体与围岩之间的波速差异,跨孔电磁波CT反映了异常体与围岩之间的电性差异,两种跨孔物探方法综合应用,对异常体双性定位,增加对异常体解析的可靠性,且能起到相互补充和相互印证作用。本文综合运用了弹性波和电磁波两种技术对孤石探测,均取得成效,但勘探结果均与实际位置存在偏差,二者结果联立解释后减小了偏差,因此,本研究成果揭示跨孔弹性波CT和跨孔电磁波CT两种资料联合处理是可行的,并初步判断联合处理成果要优于单一处理结果。

4 结语

本文简要介绍了孤石成因及其给城市轨道工程带来的问题,说明孤石探测的必要性,实验数据表明孤石与周围介质之间存在波速差异和对电磁波的吸收差异,说明可以利用跨孔弹性波CT和跨孔电磁波CT层析成像技术进行探测。本文重点介绍了跨孔层析成像技术的工作方法和数据处理算法,说明二者数据数理具有一致性,并进行了详细公式推导,工程实践证明采用相同算法对跨孔弹性波和电磁波CT数据处理得到理想效果,并对试验成果作了详细分析,说明两种方法综合利用对探测孤石具有应用前景,也为两种资料的联合反演提供参考。

摘要：在工程建设过程中常遇到花岗岩球状风化体不良地质问题,城市轨道交通建设尤其明显,严重影响着盾构法隧道施工的效率和安全,因此,发展孤石探测技术意义重大,但目前仍是一项技术难题。本文以孤石探测技术为研究目的,提出从孤石的波速特征和对电磁波的衰减特性进行综合探测,文中详细论述了数据处理的正、反演理论,并重点介绍了波速层析和衰减层析两种算法的高度一致性,通过实例分析说明两种探测技术均对孤石探测存在成效,研究成果说明综合利用两种物性条件进行孤石探测可增加探测精度,成果相互补充,相互印证,最后以本研究为基础展望了基于波速和衰减系数的联合反演技术。