非线性弹性材料（共7篇）
非线性弹性材料 篇1
高空长航时无人机(全球鹰等)为了追求更长的巡航时间、更大的航程,要求升阻比大、结构轻质的机翼,大展弦比的复合材料机翼往往是其首选。这种细长的机翼具有较大的柔性,在气动载荷作用下会产生显著的弹性变形,并且有明显的几何非线性效应。此外,复合材料的各向异性也增加了气弹分析难度。但是,通过合理设计复合材料结构的刚度分布,就可以充分利用气动弹性效应,达到减轻结构重量、提高飞行性能的目的。为此,就需要建立简单、精确的机翼气弹模型和高效的气弹优化方法。
气动弹性优化涉及结构、气动、优化等多个领域。1990年以来,国内外学者在这些领域进行了大量研究工作。在结构建模方面,以Patil为代表的一些国外学者将大展弦比复合材料机翼模拟成低阶、精确的梁模型;这些模型从早期的线性模型[1,2]发展到考虑大变形的非线性模型[3,4]。在气动建模方面,基于升力面理论的偶极子网格模型被广泛应用[3,5];这种模型比早期基于片条理论的模型精度更高,而计算量远远小于CFD方法。复合材料机翼的剪裁优化主要通过各种优化方法对复合材料的铺层厚度、角度以及铺层数进行调整以达到减小机翼质量或提高气弹性能(发散速度、颤振速度)的目标[5,6]。优化方法主要有数学规划法和准则法。数学规划法寻优过程复杂,在处理多学科耦合的气弹优化问题时计算耗时,令人难以接受;而准则法的收敛速度快,更适用于气弹工程优化。
鉴于问题的复杂性,从工程设计的角度考虑,要求在保证合理精度的前提下尽量减小计算量。因此,选择的结构模型、气动模型及优化方法需要简单、精确,并且在复杂程度和计算耗费上相匹配。
1 气弹分析及优化方法
本文采用一种非线性盒形梁模型和偶极子网格法相结合的气弹模型,使用均匀导数优化准则法对大展弦比复合材料机翼进行气弹优化剪裁。
1.1 气弹分析方法
机翼被模拟成小应变中等变形的几何非线性盒形梁结构。应用Hamilton变分原理,在变形后坐标系下建立机翼运动平衡方程,机翼变形前后坐标系如图1。
式(1)中δU、 δT、 δW分别是机翼的应变能变分、动能变分和外力虚功。
考虑截面翘曲和剪切变形,位移和应变的几何非线性关系[7]:
其中λT是截面的翘曲函数。
由于结构是小应变,本构方程还是线性关系:
将以上控制方程组转换到变形前坐标系,采用有限元离散求解,最终得到机翼离散方程:
式(6)中M是机翼质量矩阵,其非线性影响很小,视作常值;K(q)是和结构位移向量q相关的刚度矩阵,它包含了几何非线性的影响;Faero是移置到节点上的气动载荷向量,它也与结构变形有关。
机翼梁单元的划分方法如图2,每个梁单元有5个节点,共19个自由度,其中两端节点各有8个自由度,3个内节点各一个自由度。
u—轴向位移, v—横向弯曲挠度, v″—横向弯曲转角, w—垂向弯曲挠度,w″—垂向弯曲转角,
首先求解机翼静气弹响应,忽略方程(6)的前一项,并将刚度项的非线性部分提出,作为等效非线性载荷:
将FNL进行一阶泰勒展开,保留其中的线性项:
式(8)中∂FNL/∂q是非线性几何关系确定的Jacobian矩阵。对于线性气动模型:
式(9)中A是气动系数矩阵。将式(8)、式(9)代入式(7),采用牛顿-拉弗森算法:
迭代得到收敛的静响应q。
在进行颤振分析时,认为机翼是在静平衡位置附近做小扰动。根据静响应结果得到机翼静变形状态下的刚度矩阵,采用模态法近似节点位移,将系统响应方程转换到广义模态空间;采用偶极子网格法求解非定常气动载荷。应用p-k法进行颤振分析,这种方法无需对马赫数进行迭代,比V-g法更适合于优化计算,其颤振平衡方程为:
方程中,MG、KG为广义质量和刚度矩阵;
求解颤振方程(12),得到不同速度下各阶模态自激振动频率和衰减率,绘出V-γ和V-ω图,V-γ曲线和横坐标的交点就是颤振发生点[8]。
1.2 优化方法
优化设计常用的方法有非线性数学规划法和优化准则法。非线性规划法通用性好,但方法复杂,计算耗时长;优化准则法基于工程经验获得的优化设计所必须满足的准则[9],方法简单,计算量小,适用于一些特定问题。常见的优化准则有满应力准则、位移准则、均匀导数准则等。本文使用均匀导数准则对机翼复合材料铺层进行优化剪裁,流程如图3所示。
以颤振速度Vf为约束条件,复合材料各铺层厚度tj为设计变量,优化机翼质量。由于优化目标是机翼质量而优化变量是铺层厚度,所以首先要将设计变量转换为铺层质量mj,它们组成了机翼的总质量。Vf对各设计变量的偏导数∂Vf/∂mj称为敏度。当Vf对各变量的敏度不相等时,应增大敏度大的变量,同时减小敏度小的变量,从而在Vf不变的情况下减轻结构总重。当Vf对各变量敏度相等时,优化计算收敛,我们就认为获得了最优的设计。优化递推公式[10]如下:
式(13)中,mjN是新的mj值,ND是设计变量的总数,ΔVf是希望经过优化后获得的Vf增量。本文使用数值差分来近似获得颤振速度对设计变量的敏度。对以下两组设计变量:
所以差分法计算敏度可得:
式(14)中,Vf (m″)和Vf (m″)为p-k法求得的颤振速度。
2 算例分析
2.1 气弹模型验证
盒形梁机翼一端悬臂固支,如图4所示。复合材料对称铺层,铺层角为 [0/±45/90]s。结构的几何尺寸和复合材料特性如表1所示。采用1.1节的方法将机翼离散为10个梁单元,共41个节点,118个自由度。
机翼气动面沿弦向2等分,展向40等分,共80个气动网格(图5)。
对以上机翼模型进行初始分析,选取包括一阶扭转在内的前五阶模态进行颤振分析,计算结果如图6所示。计算结果显示,扭转模态是机翼颤振的危险模态,因此结构的扭转刚度对颤振速度影响很大。优化的关键在于保证结构刚度的前提下寻求质量最轻的设计方案。
为验证模型正确性,采用Nastran软件建立机翼壳单元结构模型(包含240个单元,246个节点和1230个自由度),并进行颤振计算。两种模型的计算结果如表2所示,机翼临界颤振速度和颤振频率的误差均小于10%,发生颤振的都是一阶扭转模态,颤振形式是一阶扭转和三阶弯曲耦合颤振。由于本文模型的单元数目和复杂程度都要远远小于Nastran壳单元模型,所以在保证了计算准确性的前提下大大提高了计算效率。
2.2 气弹剪裁
复合材料铺层结构如图7所示。通过调整各方向铺层的厚度t0、t45和t90,使得机翼在满足颤振速度约束条件下,达到质量最轻。机翼飞行条件如表3所示。颤振速度下限200 m/s,使用均匀颤振速度导数准则进行优化,优化初始值如表4。
由于优化目标是机翼质量而优化变量是铺层厚度,所以首先要将优化变量转换为0o、45o和90o铺层的质量。机翼由这三个质量部件组成,其优化过程及结果如图8～图11所示。
在使用均匀颤振速度导数准则迭代了四轮以后,对迭代收敛条件进行判定,如表5所示。此时颤振速度处于约束边界上,对各部件质量敏度误差小于5%,基本相等。所以认为优化计算已经收敛。
最终的优化结果如表6所示,机翼质量从0.695 2 kg降低到了0.444 6 kg,下降了36%。45o铺层的厚度增加以保证扭转刚度,使结构满足颤振速度约束条件;0o和90o铺层由于对结构刚度贡献较小,所以减小其厚度以减轻结构总质量,其中90o铺层厚度减小最多,0o铺层次之。优化前后结构固有频率如表7所示,可以发现一阶扭转频率显著增加,而弯曲频率均有所降低。机翼依旧在一阶扭转分支发生颤振。
3 结论
本文将大展弦比复合材料机翼模拟成非线性盒形梁结构进行气弹分析,并使用优化准则法对机翼复合材料进行了气动弹性优化剪裁,得到以下结论:
(1)本文的气弹模型可以模拟柔性机翼的非线性气弹特性,具有精度好,计算量小的特点,适用于气弹分析及其优化的工程计算。
(2)对于本文研究的盒形梁复合材料机翼,随着速度的增加,一阶扭转模态总是首先失去稳定,导致颤振发生。优化计算中发现,颤振速度对复合材料±45o铺层的敏度最大,0o铺层次之,90o铺层最小;所以增大±45o铺层保证结构扭转刚度以满足约束条件,减小0o和90o铺层优化机翼质量。
(3)均匀导数准则法进行气弹优化具有方法简单,计算效率高的优点。
摘要：将复合材料机翼模拟成盒形梁。考虑变形的几何非线性,建立简单、精确的气弹模型进行优化裁剪研究。以颤振速度为约束条件,运用均匀导数准则对复合材料各铺层厚度进行优化,以获得质量最轻的设计。研究表明:机翼一阶扭转模态是颤振危险模态。颤振速度对复合材料45o铺层敏度最大。优化结果中增加了45o铺层厚度以增加结构扭转刚度,减小了0o和90o铺层厚度以优化质量。将低阶非线性梁模型与均匀导数准则相结合的方法具有精度合理、收敛快的优点,适用于复合材料机翼的气弹工程优化。
关键词：非线性梁,复合材料机翼,气弹优化,均匀导数准则
参考文献
[1] Weisshaar T A.Aeroelastic tailoring of forward swept compositewings.AIAA—80—0795,1980
[2] Taylor J M,Butler R.Optimum design and validation of flat compos-ite beams subject to frequency constraints.AIAA—96—1583—CP,1996
[3] Patil M J,Hodges D H.On the importance of aerodynamic and struc-tural geometrical nonlinearities in aeroelastic behavior of high aspectratio wings.Journal of Fluids and Sturctures,2004;19:905—915
[4] Smith M J,Patil M J,Hodges D H.CFD-based analysis of nonlinearaeroelastic behavior of high-aspect-ratio wings.AIAA—2001—1582,2001
[5]万志强,杨超.大展弦比复合材料机翼气动弹性优化.复合材料学报,2005;22:145—149
[6] Patil M J.Aeroelastic tailoring of composite box beams.AIAA—1997—0015,1997
[7] Bir.G,Chopra.I.UMARC theory manual.Maryland:Center forRotorcraft Education and Research University of Maryland CollagePark,1994
[8]陈桂彬,邹丛青,杨超.气动弹性设计基础.北京:北京航空航天大学出版社,2004;106—108
[9]管德.飞机气动弹性力学手册.北京:航空工业出版社,1994;151—155
非线性弹性材料 篇2
关键词：能量释放率,能量图形面积法,加载方式,线弹性材料,非线性弹性材料
引言
断裂力学教材[1,2,3,4]讲述平面裂纹构件的能量释放率时，一般只针对线弹性材料，加载方式或者是固定位移，或者是固定载荷.著名的Irwin--Kies公式也是在上述条件下推出的，于是能量释放率可以通过构件的柔度对裂纹长度的导数求得，与加载方式无关.但是，如果对既不是固定位移也不是固定载荷的任意加载方式，或者更进一步，材料又是非线性弹性的，能量释放率如何计算，以及是否与加载方式有关，有必要弄清楚.参考文献[5]用解析方法给出相关的全面的公式推导，以便于进一步的理论推导使用，但是要求读者有相当的数学基础.本文从不同的视角出发，采用比较形象的能量图形面积法加以推导，理解容易且使用方便，证明即使对非线性弹性材料，虽然不同加载方式下的能量释放率的计算公式在形式上有所不同，但在极限求导时它们的数值大小是一样的，差别只是高阶小量.参考文献[5]中的公式与本文的公式如果取极限写成解析公式，本质上应该是一致的.从2012年起，为了推进国际化教学，我们在本科生和研究生的断裂力学的教学中启用了美国著名的康奈尔大学(Cornell University)的最新原版教材[6]，本文吸纳了该英文教材的优点，是讲授双语课程的教学体会之一.
1 线弹性材料
能量释放率是势能对裂纹面积的变化率[6]，如式(1)所示
其中，W是外力的功，U是弹性应变能，A是裂纹面积，A=Ba,a是裂纹长度，B是厚度，B=1时，A=a.∆代表增量，很小的增量也可以用d代替∆.
1.1 固定位移加载(A→B，载荷必随da改变)
如图1所示，由于位移δ1固定，外力不做功，故∆W=0，∆W-∆U=(1).其中：“初”表示裂纹长度是a的状态，“末”表示裂纹长度是a+da的状态.∆U为末状态(裂纹扩展da后)下弹性应变能与初状态下弹性应变能的差,即,(2)–((1)+(2))=–(1)；符号(1)和(2)分别表示对应图形的面积，具有功和能的量纲，即力乘长度，于是
d P表示固定位移下，与da相对应的P的增量.
1.2 固定载荷加载(A→C，位移必随da改变)
如图2所示，加载载荷P不变，初末状态如A,C所示.
外力随位移的改变而做功
(1),(2),(3),(4)分别表示对应图形的面积.于是
1.3 任意加载(A→D，位移和载荷都随da改变)
如图3所示，既不固定位移也不固定载荷
则有：∆W-∆U=((3)+(4))–((3)–(1))=(1)+(4)，于是
dδ表示固定载荷或者任意加载下，与da相对应δ的增量，固定载荷时dδ=δ2-δ1，任意加载时dδ=δ3-δ1.
对于后两种情况都比第一种情况的∆W-∆U=(1)多出(4)，但总有面积(4)1/2(d P·dδ)，而面积(1)=1/2(d P·δ1),因为dδδ1，故有面积(4)面积(1)，故以上3种情况算出的G数值上相等，都是G=(1)/da.注意到在参考文献[4]的第56页有一条译者注：“原文为‘这一面积是与加载方式无关的’，疑有误.”，此书的中文译者把英文原文的意思改为：“这一面积的大小是取决于加载方式的.”其实原文并没有原则性的错误.根据上述分析，固定载荷或者任意加载只比固定位移算的有关面积略大一点，在极限(即求导)意义上此微量趋于零(见式(1)的第3个等式).
总结以上3种加载情况可知，计算能量释放率G的关键量∆W-∆U就等于裂纹长度为a和a+da的载荷-位移曲线与加载路径所围三角形的面积，在图1，图2和图3中，此三角形的两边是载荷-位移曲线，而第3边是加载路径.G为状态量，3种加载情况都能得到相同的G，而按固定位移方式求解最方便.在极限求导的意义下，G值与加载路径无关，都等于固定位移下的G值.
2 非线性弹性材料
同理可以证明，对非线性材料裂纹构件的载荷-位移关系式曲线，但仍可用式(1)求其能量释放率，上述结论仍然成立，唯一的区别是相应于线性材料的直边三角形现在变成了曲边三角形，故面积要用积分来算.
2.1 对恒位移加载(A→B，载荷必随da改变)
如图4所示，位移δ固定，外力不做功，故∆W=0，于是
从而
其中P即为P(a,δ)曲线.
2.2 对恒载荷加载(A→C，位移必随da改变)
如图5所示，加载P固定.于是
则∆W-∆U=(3)+(4)–[((2)+(3))-((1)+(2))]=(1)+(4)
即
从而
其中δ即为δ(a,p)函数.
2.3 任意加载(A→D，位移和载荷都随da改变)
如图6所示，任意加载途径下，位移、载荷都将随裂纹长度改变，此时
则∆W-∆U=(3)+(4)–[((2)+(3))–((1)+(2))]=(1)+(4)
即∆W-∆U是计算能量曲边三角形OA D=(1)+(4)的面积，可有两种算法如下
OA D=OA B+A B D(算法1)
OA D=OA C-A C D(算法2)
相应的能量释放率也有两种算法，即:
(1)算法1
其中P(δ)是加载路径.
(2)算法2
其中δ(p)是加载路径.
分别对(7a)和(7b)用积分中值定理有：
(1)算法1
其中P=[p(δ)-p(a+da,δ)]，P0是位于P1和P2之间的某一个值，P 0是位于P1和P3之间的某一个值，易知他们都是有限值.
(2)算法2
其中，δ=[δ(a+da,p)-δ(p)],δ0是位于δ1和δ2之间的某一个值，δ0是位于δ1和δ3之间的某一个值，易知它们都是有限值.
从图6可以看出，由于(δ3-δ1)δ1，式(8a)中的第2项要比第1项小得多；同理由于(p1-p3)p1，式(8b)中的第2项也要比第1项小得多.因此，若略去高阶小量，则式(8a)与式(5)的计算值相等.而式(8b)与式(6)的计算值相等.事实上，当da→0时，相对于三角形OA B而言，三角形A B D,A C D以及A B C都是高阶小量.因此，对于任意加载路径(图6中的A→D)，由于起始点A(δ1,P1)是相同的，无论用式(5)或式(6)或式(7)算出的G都是一样的，只要知道P(a)-δ(a)，也即P(a,δ)或δ(a,p)曲线关系即可计算出G.
2.4 退化到Irwin--Kies公式
Irwin的重要贡献之一就是推出Irwin--Kies公式，即G与柔度C对裂纹长度的导数有关
可以证明，对线弹性材料，分别利用式(5)和式(6)，以及，δ=PC(其中C(a)只与a有关)就可以推出式(9)，由于式(9)对任意加载路径都适用，这也表明G的计算与加载路径无关.
3 结论
(1)从能量释放率G的定义式(1)出发，用图解法，即能量图形面积法，推导出能量释放率G的计算公式.用于裂纹扩展的能量就是某一个三角形的面积，此三角形的两边分别是裂纹长度为a和a+da的载荷-位移曲线，第3边则是加载路径.线性弹性材料不同于非线性弹性材料之处仅在于，与a和a+da对应的那两边是直线，而非线性弹性材料对应a和a+da的两边是曲线.见图1∼图3和式(2)∼式(4).
(2)对非线性弹性材料，P(a)--δ(a)曲线下的面积是非线性弹性应变能，它是能量状态函数，在恒位移或恒载荷条件下，求出它关于裂纹长度a的导数就可以求出G.
(3)对非线性弹性材料在恒位移、恒载荷和任意加载的能量释放率公式分别是(5)∼(7)，经证明得出，在极限求导时(da→0)三个公式的计算值相等.
(4)当式(5)和式(6)退化到线弹性材料时，就得到式(9)，这就是著名的Irwin--Kies公式，只要对裂纹体的柔度求导即可得到G，它可以对任意加载路径适用[1,2,3,4].
参考文献
[1] Broek D.工程断裂力学基础.王克仁等译.北京:科学出版社,1980
[2] 沈成康.断裂力学.上海:同济大学出版社,1996
[3] 程靳,赵树山.断裂力学.北京:科学出版社,2006
[4] Lawn B.脆性固体力学(第2版).龚江宏译.北京:高等教育出版社,2010
[5] 张作启,刘彬.任意加载模式下含裂纹超弹性体的能量释放率.力学学报,2013,45(1):129-133
非线性弹性材料 篇3
岩石物理实验在当今的地球物理勘探和地震资料解释中越来越凸显其举足轻重的地位, 对于岩石物理这一基础学科的研究有助于从根本上了解岩石本身性质及其内部结构, 对地震波速度异常的解释具有十分重要的作用, 实验作为岩石物理基础内容, 对其研究也就有了现实指导意义。
按照胡克定律, 我们可知岩石体应变与体应力之间的关系:σ=Eε, 也即是体应变与体应力之间存在线性关系。非线性弹性[1], 即是指在实际岩石物理实验中, 宏观表现上体应变和体应力之间不是线性关系, 这种非线性弹性关系存在使得岩石物理实验存在诸多变数。对岩石的非线性弹性的研究将有助于提高岩石物理测定精度和操作规范性。对该类问题的研究主要从两个方面着手: (1) 是岩石内部的空隙, 岩石内部普遍存在诸多空隙, 这些空隙按照加压后的“反应”分为两类, 第一类是加压后就逐渐闭合的空隙, 称之为“软空隙”, 第二类是加压后 (本实验取30Mpa) 不易闭合的空隙, 称之为“硬空隙”。[2]这些空隙的存在是岩石非线性弹性的根本原因, 在本实验中, 测试对象致密砂岩的声波速度主要与空隙相关, 对其声波速度的测试能直观反映空隙的变化情况; (2) 滞后效应[3], 是指事物的发展需要经历由微到著, 由潜到显的变化过程, 致使在复杂事物的因果链条中, 事物的发生原因与事物发生的实际响应存在一定时间的间隔, 我们将这种延迟了较长时间的现象称为滞后效应。滞后效应在岩石物理实验中具体表现为压力已到达目标值, 岩石样品内部压力未达到目标值从而使得宏观表现滞后于压力的变化。由于滞后效应的存在, 导致了实践结果的多样性与不确定性。
1 实验及结果分析
在目前的岩石物理实验中, 测定岩石在不同压力情况下岩石速度实验极为常见, 在此实验中, 滞后效应的影响不容忽视, 它作为岩石的一种重要性质, 对滞后效应的研究对岩石物理实验具有十分现实的意义。本文所选用的实验对象为来自CQ地区的一批致密砂岩, 实验前均将实验样品至于80°干燥箱中干燥24小时, 确保足够干燥状态, 此外, 本实验仪器选用的是SCMS-E高温高压岩心多参数测试仪, 实验采用的是脉冲穿透法, 纵波发射频率为700KHZ, 按照5mpa为梯度, 实验获取不同压力 (0-60MPA) 下岩石的纵波速度。
实验结果图:
实验结果分析:
图1是选取CQ地区的9快致密砂岩样品测试成果图, 不同的岩石样品纵波速度随压力变化趋势图;图2则是一块岩石样品的压力-纵波速度曲线图, 在本图中, 分段拟合了30Mpa前和30Mpa后的变化趋势, 从图中我们可以看到: (1) 0-30MPA实验结果拟合成指数曲线, 即纵波速度与压力的增加不呈线性关系, 而是表现为非线性关系, 岩石物理性质上表现为非线性弹性; (2) 对30-60MPA的实验结果进行拟合, 发现可以用线性直线很好的拟合实验结果, 即表现为线性关系, 在岩石物理性质上为线性弹性; (3) 实验等待稳压的时间越长, 实验结果与拟合曲线更为相近。
岩石内部空隙的“软空隙”在0-30MPA阶段, 随着压力的增加, 逐渐由“开”变为“关”, 随着压力的增加, 关闭的速率逐渐放慢, 可关闭的“软空隙”越来越少, 岩石等效密度逐渐变大, 变大速率逐渐变小, 岩石的速度也随着逐渐变大, 变大速率逐渐变小;在30-60MPA阶段, 此时“软空隙”已全部关系, 由“硬空隙”和岩石基质骨架组成的部分变现出线性弹性性质, 即遵循胡克定律中应变与应力关系:σ=Eε, 此时岩石等效密度与应力也呈线性关系, 宏观变现为岩石速度与压力呈现线性增加关系。对比同一块岩石不同稳压等待时间可以看到, 随着等待时间的增长, 岩石速度与围压的关系更加符合拟合结果, 这表明等待稳压的时间越久, 岩石声波速度测定的精度越高。
2 意见和建议
岩石的非线性弹性性质是一种极为常见的岩石特征, 特别是对于岩石物理实验这一类高精度实验项目而言, 更是一个不可忽略的因素, 因此, 作者就此对岩石物理实验提出几点建议和改进意见: (1) 岩石物理实验需本着实事求是的精神, 认真记录每一次观测到的实验结果, 确保第一手资料的真实可靠; (2) 岩石物理测试过程中, 不同的测试人员需统一标准, 特别是在选取初至的时候, 统一选择某一个特征清晰明显的起跳点作为初至点, 一次实验过程和全部实验过程都保持这一标准, 可尽可能使得测试的结果准确, 避免人为的误差; (3) 岩石物理实验的加压和泄压过程都是一个动态的过程, 因为滞后效应的存在, 从而使得在该动态过程中, 在压力达到某一个测试目标压力值附近时, 在实际情况允许范围内尽可能选择长的等待稳定时间, 且保证每次等待时间大致相同, 尽量减少滞后效应对每单次测量结果的影响; (4) 建立滞后效应模板, 对某一沉积稳定, 岩性稳定的区域, 可以选择进行岩石样品的长时间滞后效应观测试验, 结合地震和实际录井资料, 绘制滞后效应模板图, 从而对以后测得的岩石速度值进行校正。
参考文献
[1]李廷, 席道瑛, 徐松林.动荷载作用下岩石非线性弹性响应研究[J].地学前言, 2006, 13 (3) :207-212.
[2]Franklin Ruiz and Arthur Cheng, OHM Rock Solid Images, A rock physics model for tight gas sand[J].THE LEADING EDGE.DECEMEER 2010, 1484.
非线性弹性材料 篇4
如何设计既能大幅度连续光滑变形, 又有足够刚度和强度的轻质可变形蒙皮[1,2]结构已成为智能变形飞行器的关键技术之一[3,4,5]。目前, 对于弹性蒙皮[1]的研究主要集中在以下两个方面:一种是波纹板式复合材料弹性蒙皮[5], 该蒙皮利用波纹扩张或收缩产生的变形累积效应实现沿波纹方向的大变形, 但是该蒙皮表面并不连续平滑, 同时垂直于波纹方向的纵向承载能力很弱;另一种是采用柔性较大的橡胶类材料制作的蒙皮, 这种蒙皮虽能满足机翼变形和气密性的要求, 但驱动方式和变形控制很复杂, 而且机翼的整体承载能力不高[6]。本文提出了基于柔性蜂窝结构的超弹性蒙皮结构设计理念, 以期得到变形柔度、承载刚度、轻质三方面性能俱佳的可变形蒙皮结构。
蜂窝芯柔性大变形问题的研究是柔性蜂窝芯超弹性蒙皮研制的基础。在单一规则蜂窝结构的研究中, 宏观上将其视为匀质连续材料, 选取其中的一个单元进行宏观结构上的等效力学和变形分析, 进而得到结构的等效弹性模量。Gibson等[7]运用Euler梁模型, 对蜂窝结构的力学行为进行分析研究, 得出了蜂窝结构弹性模量的经典计算公式 (Gibson公式) 。富明慧等[8]在考虑蜂窝壁板伸缩变形的情况下, 对Gibson公式进行了修正。Warren等[9]依据蜂窝结构中单元周期性重复排列的特点, 取其代表单元进行分析, 建立了宏观的弹性本构方程, 得出了W-K公式。王飞等[10]根据均匀化理论, 并结合有限元方法得出了不同相对密度下蜂窝结构的等效弹性参数。以上研究多以小变形性为前提假设, 但小变形性难以满足超弹性蒙皮的大变形要求[7,11]。柯映林等[12]将蜂窝壁板的变形归属于薄板大挠度变形问题, 得出了具有非线性响应的蜂窝芯材料等效面内弹性模量。Hu等[13]通过分析蜂窝壁板的扭转变形, 发现弹性模量随应变的变化已经不再是常数。祝涛等[14]考虑蜂窝壁板面内荷重对壁板弯曲的影响, 对Gibson公式进行了修正, 提出了一种蜂窝芯层的非线性等效拉伸弹性模量的拟合方法。目前, 对蜂窝结构非线性大变形问题的研究尚无成熟的理论体系形成, 实验研究则更为有限。
本文设计并加工了内六角形负泊松比蜂窝结构, 通过理论建模、有限元仿真和力学实验3种方式对这种蜂窝结构大变形下的面内弹性模量进行了分析, 建立了弹性模量的理论计算公式, 得出了等效弹性模量的非线性特性及相同方向和不同方向弹性模量的变化特性。
1 蜂窝结构变形力学建模
1.1 Y方向单向压缩变形
图1为超弹性柔性蜂窝芯蒙皮结构原理示意图和蜂窝单元图。其中, a为横壁板长, b为斜壁板长, t为斜壁板厚度, w为横壁板厚度, θ为斜壁板与横壁板的夹角。
假设横壁板为刚性 (w≥2t) , 蜂窝结构的变形由斜壁板的弹性变形引起。选取斜壁板OCB为研究对象, 将其视为一端固支 (O点) , 一端限制其转角 (B点) 的细长杆 (t/b≤0.015) , 其受力变形图见图2a。由于整个壁板的受力变形成反对称性, 故取其一半建立柔性悬臂梁模型进行分析 (图2b) 。图2b中γ为任意截面的转角, s*为斜壁板变形后的弧长, α为载荷P (与Y轴平行) 与X*轴的夹角。
悬臂梁的挠曲线方程为
式中, Es为材料本身的弹性模量;I为斜壁板的惯性矩。
引入, s=s*/b (Pcr=π2EsI/b2, 0≤s≤0.5) , 并代入边界条件 (C点的弯矩为零) 可得
令, γ1为壁板的最大转角 (即C点处的转角) , 则
对s进行积分和坐标转换, 可得其在X、Y方向上的投影为
1.2 Y方向单向拉伸变形
Y方向单向拉伸时, 斜壁板受力变形如图3所示。
悬臂梁的挠曲线方程为
采用上一节同样的方法进行推导可得量纲一杆长s在X、Y方向上的投影为
对于斜壁板X方向的受力变形情况, 可依照斜壁板Y方向受力分析的方法进行分析。
2 蜂窝芯等效弹性模量
利用上节的结论, 可求得相应大变形的格林应变, 进而求得等效弹性模量。
分别用上标1、2表示蜂窝材料单向压缩和拉伸的受力状态, 下标x、y表示受力方向分别为X、Y向, 则等效弹性模量公式为
3 理论计算、ANSYS仿真、模型实验结果分析
采用65Mn弹簧钢加工出单个单元实验模型, 在弹性变形范围内进行了X方向和Y方向的拉伸、压缩实验, 模型受力变形如图4所示。通过测量载荷力和位移计算出应力和应变, 进而求得等效弹性模量。选择满足超弹性大挠度和大应变能力的Soild187单元对蜂窝结构进行有限元仿真分析。根据实验模型的材料、尺寸参数和载荷, 设定有限元模型的材料、尺寸参数和载荷量级, 通过数值模拟计算出应变, 进而求得等效弹性模量。模型的参数为:Es=196.2 GPa、a=55 mm、b=43mm、t=0.4 mm、h=18.5 mm (蜂窝芯层的厚度) 、θ=70°。并将理论计算 (式 (7) 计算结果) 、ANSYS有限元仿真、模型实验、线性计算的结果[7,15]进行比较分析, 给出应力应变、等效弹性模量应变关系, 如图5~图8所示。线性计算公式为
由图5~图8可得如下结论:①图5~图8表明, 理论、仿真、实验三条曲线吻合程度良好, 证明了等效弹性模量计算公式的正确性和有效性;②从图5、图6可以看出, Y方向的等效弹性模量非线性特征明显。单向压缩、拉伸应变趋于零时等效弹性模量趋于线性计算的固定值;③从图7、图8可以看出, X方向等效弹性模量近似于线性变化, 应变趋于零时等效弹性模量趋于线性计算的固定值, 随应变的增大等效弹性模量和线性计算值的差值增大, 因此, 大变形条件下柔性蜂窝芯X方向等效弹性模量的计算应选用更为精确的非线性计算公式;④图5~图8表明, 同一方向单向拉伸和压缩的等效弹性模量随应变的增大变化趋势相反, 应变越大两者差值越大, 因此, 柔性蜂窝芯等效弹性模量的计算应根据受力方向和受力状态选择相应的计算公式进行计算。
4 结束语
非线性弹性材料 篇5
非线性波动的研究对于解决物理学、化学、生物学和地球物理学中遇到的复杂现象和问题有着极其重要的意义.孤立子理论的建立是非线性波动理论发展中的一项重大成就.在非线性物理许多领域,已经发现一大批非线性演化方程具有孤立子解.孤立子的典型特征是在其传播过程中伴随有能量的集聚,且孤立子间相互作用时表现出犹如粒子弹性碰撞一样的行为.这些特征在流体动力学、光纤通讯、等离子体及生物学等领域已开始获得了一些应用[1].近年来,对于固体结构中的应变孤立波的研究也取得了某些进展[2,3].本文将研究涉及固液耦合作用的埋置于弹性地基内的充液压力管道中非线性波的传播,寻找其孤立波解.
充液管中波的传播一直是一个十分活跃的研究领域,管壁的膨胀性和流体的压缩性可导致压力波沿管的轴向传播.依管的刚度和流体的可压缩性不同,所涉及问题大致可分为两类.如果管壁是刚性的,且流体的压缩性是不可忽略的,这类问题构成了水锤理论的基础,它在化工、核电等工业部门有着重要的作用.如果管壁刚性相对较弱,而流体几乎不可压缩时,这类问题相应于诸如血管中和输油管中脉冲传播的力学模型.本项研究属于后一类问题.充液管中压力波的传播研究至少有两百年的历史,较早的研究只能处理线性问题,如充液弹性薄壁管中波的传播速度和周期波的传播以及预应力管、黏弹性管和变截面管中波的传播问题.20世纪80年代以后,由于非线性科学发展的推动,开始了充液管中的非线性波、特别是孤立波的研究.过去二十多年中,根据所关心的问题不同,考虑不同效应的耦合,采用不同的力学模型,得到了一批有意义的成果.文献[4]是最早研究充液弹性薄管中非线性压力波的文献之一,基于曲面理论,借助张量工具导出了管壁运动方程.另外,特别值得一提的是H.Demiray及其合作者在这一领域发表了大量论文[5,6,7,8],在他们的模型中,大多数是假定管壁中有预加应力.
本文研究埋置于弹性地基内充液压力管道中的非线性波,对模拟地下输运管线或包裹在肌肉内的血管中的血液流动等有潜在的应用背景.本研究中,假设流体为不可压缩理想流体,管外地基反力采用Winkler模型,认为地基每单位面积的反力Pe与管壁的法向位移ω成正比,即Pe=kw.对于管壁的分析,使用了类似于文献[4,9,10]在充液管中压力波的研究中所采用的一个假定,即不计轴向应力的影响.正如文献[4]评述的那样,这个假定对于描述某种现象可能不够充分,然而对血管中脉冲波传播的绝大多数特性均能从这个模型得到合理的解释.另外,还假设管壁是橡胶类材料或软组织材料,在其变形过程中考虑了半径和壁厚的变化,从而导致壁厚运动方程是非线性的.初始时,假定在内压P0的作用下,管壁和流体均处于静止平衡状态;管壁中的动力位移场及流体流速和压力的变化是叠加在这个平衡状态上的微小扰动.在长波近似条件下,流体的压力和速度沿半径方向进行平均,仅沿轴向是变化的,即流体流动是一维的.基于该物理模型,建立其非线性数学模型,并利用约化摄动法求解该模型.
1 控制方程
1.1 管内流体的质量守恒方程
考虑一维流动,管中流体的质量守恒方程为
式中,ρf为流体密度,A为管内横截面积,V为流体轴向流动速度,x和t分别为轴向坐标和时间坐标.若流体是不可压缩的,则式(1)变为
1.2 管内流体的动量守恒方程
考虑沿x方向的流体动量守恒,对于理想流体有
式中,Pi是管内流体压力,Pi和截面积A之间的关系可以是相当复杂的,为了简单起见,引入一个简化假定[9],即认为面积A仅依赖于过剩压力Pi-Pe,即
其中,Pe为管壁外部压力.
1.3 管壁的法向运动方程
为清楚起见,将管壁的状态或构形分为3种情况(如图1).当Pi=0时,从而Pe=0,管处于无应力自然状态,称其为原始构形,此时管的半径和厚度分别记作R0和H0;当内压Pi=P0时,管壁的法向位移为w0,地基反力为Pe=kw0,此时的管半径和厚度分别记作r0和h0,流体速度V=0,这是一个静力平衡状态,称其为中间构形,或参考构形;从t=0时刻开始,系统进入扰动状态,我们称此后时刻t的构形为瞬时构形或现实构形.在瞬时构形中,内压Pi=Pi(x,t),Pi=Pi-P0为内压的扰动,流体速度为V(x,t),管的半径为R,厚度为h,法向挠度为w0+w,w=w(x,t)为从参考构形算起的挠度,地基反力为Pe(x,t)=k(w0+w).对于软组织材料,认为管壁不可压缩,则有R0H0=r0h0=Rh.根据这3种构形,容易得到
以原始构形为基准度量环向应变,于是有
式中,和εθ分别为中间构形和瞬时构形中的环向应变,为中间构形的环向伸长比.相应的环向应力为
Δσθ为从中间构形到现实构形环向应力增量,E为弹性模量.
在图2所示的现实构形中,考虑动力平衡,有
式中,P=Pi,ρω是管壁的密度.利用几何关系(6)和弹性关系(7),经适当整理后,式(8)可写成
式(9)是关于半径R为变量的运动方程.在中间构形上w=0,R=r0,式(9)退化
这正是参考构形的静力平衡方程(Laplace定律).
为了与方程(9)相协调,方程(2)也应变成关于R的方程,为此将A=πR[2]代入方程(2),可得
2 方程的综合及其求解
综合方程(3),(9)和(11),经适当整理可得关于埋置于弹性地基内充液压力管道中非线性波的动力学方程组
其中,第2式中用P代替了Pi.方程(12)是关于3个变量P,V和R的方程组.除了流体的对流项引起的非线性,上式第3个方程表明管壁方程也是非线性的.线性化表明,方程组(12)具有弱弥散性,应用约化摄动法求解.为此,采用以下G-M变换[11]
相应地,有
将式(14)作用于方程组(12),得到
假定方程中的各场量可以表达为如下渐进展开
式中ε为小参数,将式(16)代入到(15),得到关于ε的一阶、二阶方程组.O(ε)阶的方程为
O(ε[2])阶的方程为
对式(17)进行直接积分,取积分常数为零(由于ξ→∞,V1=R1=P1=0),可得
可以看出,一阶摄动得到了R1,P1和V1之间的关系,同时也给出了c0的表达式,要c0为实数,则要求β>0.将式(20)代入到二阶的方程(18),有
其中,.从式(21)中消去P2,R2,V2,得到
其中,.式中左边第1项为非定常项,第2项为非线性项,第3项为弥散项,α为非线性系数,γ为弥散系数.若令u=αU,则可得以下标准形式的KdV方程
u和U只是幅值相差α倍.KdV方程(22)有如下形式的孤立波解
可以看出,波速c与波的幅值有关,这是非线性波的特点.若要求c>0,则γ>0.
3 结论与讨论
(1)本文在长波近似情况下,研究了埋置于弹性地基内充液压力管中非线性波的传播,借助约化摄动法从流固耦合的动力学方程组中得到了KdV方程.结果表明,该流固耦合系统在一定条件下会有孤立波解存在.由于孤立波有许多重要特性,因此得到的结果会在生物医学工程或其它工业部门的相关问题研究中有一定参考价值.
(2)从解的表达式(24)看出,方程(22)中的参数α和γ对孤立波的形成与传播有重要作用.α与波的幅值成反比,波的宽度与γ的平方根成正比.显然,α和γ的值越小,波就会越陡峻(steeping)和峰化(peaking),这些现象应该是孤立波应用的重要指标.另外,波的幅值与波的传播速度成正比,这是非线性波的基本特征.如在G-M变换中取c0>0,则方程(21)中的a>0,从而有γ>0,α和γ同号就要求,进一步分析可得孤立波存在的条件是
在上式令k=0,并利用式(10),可以给出
式(26)为无弹性地基存在时的孤立波存在的条件.分析式(25)可以看出,弹性地基的存在,扩大了式(26)对ω0的要求范围.顺便指出,地基系数k过大,相当于大大增加了管壁的刚度,当其到达一定程度时,流体压缩性便不可略去,此种情况已不是本文研究所关心的问题.
(3)在已有相关研究中,描述管内流体运动的方程基本相同,大都是一维流动,至多再引入流体的黏性.对于描述管壁的模型,类型就比较多.一般认为管壁有软材料制成而不计其弯曲刚度,又考虑到周围基体的轴向束缚,认为轴向位移x=0.本文主要考虑基体法向束缚而引入了弹性地基,对管壁采用了类似于文献[4,9,10]的处理,略去了轴向膜应力作用,而考虑到软材料的特性,通过考虑半径和壁厚的变化引入了非线性.
摘要：研究了埋置于弹性地基内充液压力管道中非线性波的传播.假设管壁是线弹性的,地基反力采用Winkler线性地基模型,管中流体为不可压缩理想流体.假定系统初始处于内压为PO的静力平衡状态,动态的位移场及内压和流速的变化是叠加在静力平衡状态上的扰动.基于质量守恒和动量定理,建立了管壁和流体耦合作用的非线性运动方程组;进而用约化摄动法,在长波近似情况下得到了KdV方程,表征着系统有孤立波解.
关键词：弹性地基,充液圆管,约化摄动法,孤立波
参考文献
[1] Michel Remoissenet.Waves Called Solitons.2nd ed.New York:Springer-Varlag Berlin Heidelberg,1999
[2] Zhang Shanyuan,Zhuang Wei.The strain solitary waves in a nonlinear rod.Acta Mechanica Solida Sinica,1987, 1(3):62-72
[3] Zhang Shanyuan,Liu Zhifang.Three kinds of nonlinear dispersive wave in elastic rods with finite deformation.Appl Math Mech,2008,29(7):909-917
[4] Hashizume Y.Nonlinear pressure wave in a fluid-filled elastic tube.J Phys Soc Japan,1985,54(9):3305-3312
[5] Hilmi Demiray.Solitary waves in prestressed elastic tubes. Bulletin of Mathematical Biology,1996,58(5):939-955
[6] Hilmi Demiray.Propagation of weakly nonlinear waves in fluid-filled thin elastic tubes.Applied Mathematics and Computation,2002,133:29-41
[7] Hilmi Demiray.Nonlinear waves in a fluid-filled in homogeneous elastic tube with variable radius.International Journal of Non-linear Mechanics,2008,43:241-245
[8] Nalan Antar.The KdV-Burgers hierarchy in fluid-filled elastic tubes.International Journal of Engineering Science, 2002,40:1179-1198
[9] Fung YC.Biomechanics:Circulation.2nd ed.New York: Springer,1997
[10] James Lighthill.Waves in Fluids.London:Combridge University Press,1978
线性时变周期系统的弹性控制 篇6
线性时变周期(LTVP)系统是一类复杂而又非常有研究的价值系统,在物理以及工程技术中有着广泛应用[1]。近些年来,国内外许多研究学者对线性时变周期系统的研究作出过很多有意义的贡献。例如,文献[2]致力于离散的周期系统的模型匹配问题的研究,而张雪峰[3]研究了线性时变周期系统的稳定性和鲁棒控制,并且成功地给出了基于LMI的LTVP系统稳定性的充要条件。
尽管近年来,国内外学者对LTVP系统展开了多方面的研究,其中,不乏许多重要且独特的见解。但是,在设计LTVP系统的控制器时,关于控器的鲁棒性的问题的研究并不多见。
弹性控制[4]是指系统能够忍受控制装置的结构和参数、系统的结构和设计参数引起的波动,以及在非正常情况下,系统可以维持状态已知,以及可以接受一定恢复程度的运行。而这些非正常情况对系统的影响,在一定程度上,可以转化为控制器的鲁棒性。
因此,本文在研究弹性控制器的设计时,就考虑到了控器的鲁棒性,致力于LTVP系统的弹性控制的研究。并且,成功地给出了LTVP系统镇定的充分条件,以及相应的弹性控制器的设计方法。
二、LTVP系统简介
考虑一个LTVP系统
其中:x(t)∈Rn,u(t)∈Rm,A(t):R→Rr×n,C(t):R→Rr×n是以T为周期的分段函数,B:R→Rn×m为常数矩阵。Rn表示n维实数。
本文利用鲁棒控制的相关思想,可以把系统(1)简化成一种鲁棒扰动形式,而该鲁棒扰动形式是范数有界的。
对系统(1)而言,A(t)是范数有界且分段连续的,即存在一个常数ρ(0<ρ<∞),使得AT(t)A(t)≤ρ2I为不失一般性,设A(t)=A+MF(t)E,其中A,M,E是具有适当维数的矩阵,且FT(t)F(t)≤I。下面给出几个本文所需要涉及的引理。
引理1[5]LTVP系统(1)被称为是渐近稳定的,如果存在一个对称正定矩阵P,使得
引理2[5]对于任意给定的实矩阵S1,S2,S3(其中S1=S1T,S3>0)有S1+S2S3-1S2T<0,当且仅当
引理3[6]对于给定的具有适当维数的复数矩阵Ω,Γ,∑,且Ω是对称矩阵,则Ω+ΓF(t)∑+∑TFT(t)ΓT<0对任意的F(t)满足FT(t)F(t)≤I,当且仅当存在一个正数ε>0,使得Ω+εΓΓT+ε-1∑T∑<0。
下面给出Riccati稳定的概念,该稳定渐近稳定更强。
定义1系统(1)可以被称为是Riccati稳定,如果对于任意给定的正交矩阵F(t),FT(t)F(t)=I存在一个对称正定矩阵P使得
三、弹性控制器的设计
考虑如下形式的状态反馈控制律:
其中,K表示控制器的增益,△K表示增益的摄动。
考虑如下两种形式的摄动:
其中,FiY(t)Fi(t)≤I(i=1,2,I)(为适当维数的单位矩阵)。
定理1考虑具有加法摄动(i)的闭环LTVP系统(7),若存在对称正定和矩阵以及正数ε1,ε2使得
成立,其中“*”是由矩阵的对称性而省略的部分,Ω=XAT+AX+WTBT+BW,
则K=WX-1,u(t)=B(K+△K)x(t)为无记忆弹性控制律,使得闭环系统
渐进稳定。
证明:A(t)=A+MF(t)E,u(t)=(K+M1F1(t)N1)x(t)由引理1可知,为使系统(7)渐近稳定,有
左右同乘P-1,并令X=P-1,W=KX,可得
再连续两次使用引理3,可得Ω+ε1MMT+ε2BM1M1T+ε2BM1M1TBT+ε1-1XN1TN1X<0,
其中Ω=XAT+AX+WTBT+BW,再由引理2,可得出不等式(6),定理得证。
定理2考虑具有乘法摄动(ii)的闭环LTVP系统(7),若存在对称正定X和矩阵W以及正数ε3,ε4使得
成立,其中,“*”是由矩阵的对称性而省略的部分,
则K=WX-1,u(t)=(K+△K)x(t)为无记忆弹性控制律,使得闭环系统(7)渐进稳定。
证明与定理1类似,此略。
四、仿真算例
例1对于定理1,我们考虑具有如下参数的LTVP系统:
通过求解定理1中的线性矩阵不等式(6)可得:
例2对于定理2,我们考虑具有如下参数的LTVP系统:
M2=[0.5],N2=[0.1 0.1]T,其余参数矩阵与例1所给出的相同。求解定理2中所给出的线性矩阵不等式(8),可得:
五、结语
本文针对线性时变周期系统,给出了一个新的二次稳定方法,并且给出了基于LMI的弹性控制器的设计方法,同时,本文考虑了弹性控制器的两种摄动。这就较为全面地解决了传统的稳定控制器,未考虑控制器的鲁棒性,而相对脆弱的问题。最后,本文给出数值仿真的算例,验证了所提出的方法的有效性。当然,其他相关问题,还需要进一步的研究。
摘要：线性时变周期(LTVP)系统是一类非常有研究价值的系统,在研究该系统的镇定性问题时,一般需要设计状态反馈控制器,但是,传统的状态反馈控制器往往没有考虑到控制器的鲁棒性,从而导致控制器相对脆弱。而本文为了解决该问题,设计出弹性控制器,成功解决了该问题。而且给出基于LMI的弹性控制器的设计方法,利用仿真实例验证了所提出方法的有效性。
关键词：线性时变周期系统,弹性控制,线性矩阵不等式,鲁棒性
参考文献
[1]秦元勋,王慕秋,王联.运动稳定性理论与应用[M].北京:科学出版社,1981
[2]Colaneri P,Kucera V.The model matching problem for periodic discrete-time system[J].IEEE Trans on Automatic Control,1997,42(10):1472~1476
[3]张雪峰,杨明.基于LMI的线性时变周期系统的稳定性及鲁棒控制[J].控制与决策,2012,27(2):291~294
[4]马静,叶泳,贾秋生.弹性控制综述[J].信息与控制,2015,44(1):67~75
[5]张庆灵,张雪峰,翟丁.控制理论基础[M].北京:高等教育出版社,2008
非线性弹性材料 篇7
目前, 有关空间机器人系统的动力学分析及智能控制的研究已得到各国科学研究人员的广泛关注, 并已取得了一定的成果。但值得注意的是, 大多数研究都建立在空间机器人系统结构中各分体均为纯刚性体的基本假设上, 并将空间机器人视为一个多刚体系统[1,2,3]。然而, 在实际的空间应用中, 空间机器人的机械臂与装配在关节处驱动该机械臂运动的电机转子之间的连接不可能为绝对刚性。同时, 在空间机器人轻型化的要求下, 带有柔性的机械臂已经得到了越来越广泛的运用。因此, 空间机器人系统实际上为刚-柔性耦合系统, 而具有柔性关节和柔性臂的空间机器人模型为最接近实际的空间机器人模型。值得注意的是, 空间机器人结构中柔性因素的存在是一把“双刃剑”。一方面, 柔性臂能够减轻空间机器人的质量, 降低能量消耗, 使机械臂获得较大的操作空间和较高的工作效率;柔性关节能够吸收空间机器人在运动过程中发生意外碰撞时受到的冲击力, 降低空间机器人的损伤。另一方面, 柔性关节和柔性臂在运动过程中所引起的弹性变形和弹性振动会对系统的控制精度和稳定性造成不利的影响。而随着空间机器人技术朝着轻质、高速、高精度的方向发展, 空间机器人的大位移刚性运动与柔性关节和柔性臂的小位移弹性变形之间的耦合作用已不容忽视, 目前, 有关刚-柔性耦合的空间机器人系统的动力学分析和智能控制方法研究已成为了科学研究的重点, 但是大多数的研究对象为刚性关节-柔性臂空间机器人系统[4,5], 或柔性关节-刚性臂空间机器人系统[6,7]。虽然有少数研究同时考虑了柔性关节和柔性臂对系统的影响, 但是其研究对象主要为地面机器人[8,9,10], 而对柔性关节-柔性臂空间机器人系统的研究仍然比较少见。尤其对于载体自由漂浮的漂浮基空间机器人系统, 该系统呈现出的非线性和强耦合性使得空间机器人的动力学建模过程比固定基的地面机器人更加复杂[11], 进而又使得惯常用于地面机器人的一些控制方法无法直接应用和推广到空间机器人的控制中。因此, 有关漂浮基柔性关节-柔性臂空间机器人的研究难度较大, 同时也更具有挑战性。
基于以上讨论, 本文同时考虑了柔性关节和柔性臂对漂浮基空间机器人系统的影响, 首先利用动量、动量矩守恒关系和拉格朗日-假设模态法建立系统的动力学方程。接着基于奇异摄动法, 将系统分解为“刚性关节-刚性臂”慢变子系统、“柔性关节-刚性臂”快变子系统和“刚性关节-柔性臂”快变子系统, 并分别针对这三个子系统设计相应的控制律来实现系统期望运动轨迹的渐近跟踪和关节、臂弹性振动的抑制。最后通过仿真实验证明所提出的混合控制律的有效性。
1 漂浮基柔性关节-柔性臂空间机器人系统动力学建模
为不失一般性, 考虑如图1所示的双柔性关节、单柔性臂的漂浮基空间机器人系统。该系统由载体B0、刚性臂B1和柔性臂B2组成。Oi (i=1, 2) 为各分体与电机转子连接的关节转动铰。建立惯性坐标系 (OXY) 及各分体Bj (j=0, 1, 2) 的主轴坐标系 (Ojxjyj) 。假设各分体在 (OXY) 平面内作平面运动。
1.1 柔性关节的简化模型
根据Spong所提出的“转子-扭簧系统”简化模型[12]:在小变形的情况下, 柔性关节可简化为一个用来连接电机转子和机械臂的刚度系数为ki的无惯量线性扭转弹簧, 结构如图2所示。此时, 当关节Oi处的电机转子转过角度φi时, 与其相连接的机械臂Bi由于受到扭转弹簧弹性力的作用, 其实际的转动角度为qi=φi-σi, 其中σi为扭转弹簧引起的关节弹性变形偏差角。而关节Oi处电机转子与机械臂之间相互作用的弹性力大小可表示为ki (φi-qi) =kiσi。
1.2 柔性臂的简化模型
假设柔性臂B2为细长杆, 忽略其在运动过程中的剪切变形和转动惯量的影响, 于是可将B2视为Euler-Bemoli梁, 并利用假想模态法[13], 将B2的横向弹性变形v (x2, t) 表示为
式中, n为保留模态数;фi (x2) 为第i阶模态函数;δi (t) 为与фi (x2) 相对应的模态坐标。
本文对前二阶模态进行分析, 于是有:n=2, v (x2, t) =ф1 (x2) δ1 (t) +ф2 (x2) δ2 (t) 。
1.3 系统动力学模型
如图1所示, 令载体B0的质心Oc0与O0重合, 其相对于O的矢径为r0, 刚性臂B1的质心Oc1相对于O的矢径为r1, 柔性臂B2上任意一点相对于O的矢径为r2, 系统总质心C相对于O的矢径为rc。Oc0与O1之间的距离为l0, Oc1与O1之间的距离为d1, Bi的长度为li。B0的质量为m0, B1的质量为m1, B2的线密度为ρ, 电机的质量可忽略不计[12], 于是系统的总质量M=m0+m1+ρl2。B2的抗弯刚度为EI。根据系统几何位置关系和总质心定理, 各分体矢径rj及其一阶导数可分别表示为
式中, ei (i=0, 1, 2) 为系统各分体主轴坐标系xi轴的轴向基矢量;e3为柔性臂B2的主轴坐标系y2轴的轴向基矢量;Rj0、Rj1、Rj2、Rj3、Rj4为系统惯性参数的组合函数。
为不失一般性, 设定漂浮基柔性关节-柔性臂空间机器人系统不受外力作用, 系统相对惯性坐标系 (OXY) 满足动量、动量矩守恒关系。假设系统的初始动量、动量矩为零, 于是系统的动量、动量矩守恒关系可表示为
式中, w0、w1分别为B0、B1的转动角速度矢量;wφi为电机转子的自转角速度矢量;J0为B0的转动惯量;J1为B1的转动惯量;Jφi为关节Oi处电机转子的自转惯量。
柔性关节的存在使得在对漂浮基柔性关节-柔性臂空间机器人系统进行动力学分析时不能再将电机转子与机械臂简化为一整体进行分析, 而需要分别对电机转子的动力学和由载体、刚柔机械臂组成的空间机器人的动力学进行分析。于是, 系统的总动能T为电机转子的动能Tφ和空间机器人的动能Tq之和, 即
忽略宇宙中微弱的重力作用, 系统的总势能U为柔性关节简化扭转弹簧的弹性变形势能Uφ和柔性臂的弯曲应变势能Uq之和, 即
基于以上的讨论, 并结合拉格朗日方程, 可获得载体位置、姿态均不受控的漂浮基柔性关节-柔性臂空间机器人系统完全驱动形式的动力学方程:
其中, φ为由电机转子的转角φ1、φ2组成的向量, φ=[φ1φ2]T∈R2;q为由机械臂的转角q1、q2组成的向量, q=[q1q2]T∈R2;δ为由柔性臂的模态坐标δ1、δ2组成的向量, δ=δ1[δ2]T∈R2;θ=q[0qδ]T∈R5, q0为载体姿态角;δ为柔性关节引起的关节弹性变形偏差角组成的向量, σ=φ-q∈R2;Jφ为电机的对角正定惯量矩阵, Jφ=diag (Jφ1, Jφ2) ∈R2×2;D (θ) 为空间机器人的对称正定惯量矩阵, D (θ) ∈R4×4·;为包含科氏力和离心力的列向量, ;Kφ为柔性关节刚度系数矩阵, Kφ=diag (k1, k2) ∈R2×2;Kδ为柔性臂刚度系数矩阵, Kδ=diag (kδ1, kδ2) ∈R2×2;τ为由关节O1、O2处电机的输出力矩τ1、τ2组成的向量, τ=[τ1τ2]T∈R2。
式 (8) 可分解为电机转子和空间机器人两部分, 即
其中, D11∈R2×2、D12=DT21∈R2×2和D22∈R2×2均为D (θ) ∈R4×4的子矩阵;C1∈R2×1和C2∈R2×1均为C (θ, θ·, qa·) ∈R4的子矩阵。显然, 式 (9) 为电机转子的动力学方程, 式 (10) 为空间机器人的动力学方程。
2 系统动力学奇异摄动分解及控制律设计
柔性关节和柔性臂在运动过程中所引起的弹性变形会影响系统的控制精度, 所引起的弹性振动会影响系统的稳定性。为了克服这些问题, 本文基于奇异摄动法, 将漂浮基柔性关节-柔性臂空间机器人系统分解为表示系统刚性运动部分的慢变子系统和表示系统柔性运动部分的快变子系统, 并分别为子系统设计控制律。其中, 慢变子系统控制律τs用来实现系统期望运动轨迹的渐近跟踪;快变子系统控制律τf用来主动抑制柔性关节和柔性臂的双重弹性振动, 保证系统的稳定性。于是, 系统的总控制律可表示为τ=τs+τf。
2.1“刚性关节-刚性臂”慢变子系统
由于对称、正定惯性矩阵D (θ) 可逆, 于是由式 (9) 、式 (10) 可解得
定义奇异摄动因子ε2=1/min (kδ1, kδ2) , 变量zδ=δ/ε2、zσ=σ/ε2和矩阵。将它们代入式 (11) ~式 (13) , 可得到
为了获得“刚性关节-刚性臂”慢变子系统的动力学方程, 需要消除系统中的柔性变量。于是令ε=0, 则可将式 (14) ~式 (16) 重新写为
式中, C1σ和C2σ分别为消除矩阵C1和C2中有关柔性关节的变量 (即令) 后得到的新矩阵;上横线“-”表示消除有关柔性臂的变量后获得的新矩阵和新变量。
由式 (17) 可解得, 将其代入式 (19) 可得到
再将所求得的zσ和zδ代入式 (18) , 便可得到“刚性关节-刚性臂”慢变子系统的动力学方程:
定义qd=[qd1qd2]T为慢变子系统的期望输出向量, 则其与实际输出向量q之间的输出误差向量, ei=qi-qdi;对时间的一阶导数, 即速度误差向量。
考虑如下的非线性滑模面:
式中, α、β为正常数;p1、u1、p2、u2为正奇数, 且p1>u1, p2>u2。
由式 (22) 可看出, 该非线性滑模超曲面实际上为常规的线性滑模面和终端滑模面的组合和改进。以往研究表明:当系统状态远离平衡点时, 线性滑模的收敛速度优于终端滑模;而当系统状态在平衡点附近时, 终端滑模的收敛速度优于线性滑模。因此, 为了使系统从任意初始状态到达平衡点的过程中能够始终获得较高的收敛速度, 本文将线性滑模面和终端滑模面进行合理结合。同时注意到, 式 (22) 还对线性滑模面进行了改进:如果令p2>u2, 则此时e的指数大于1, 从而进一步加快了远离平衡点的系统状态的收敛速度。
将s对时间求导, 可得
为了消除滑模自身的抖振并提高趋近速度, 选取如下的双幂次趋近律:
式中, si为s的第i个元素。
于是结合式 (23) 和式 (24) 可得到慢变子系统如下的滑模控制律:
定理对式 (21) 所描述的慢变子系统, 滑模控制律式 (25) 可保证:当系统到达滑模面后, 对给定的任意初始状态e (0) , 系统都将保持稳定并在有限时间内到达平衡点。
证明在滑模面上, 令s=0, 则由式 (22) 可得到系统误差的收敛速度表达式:
选取如下形式的Lyapunov函数:
将V对时间求导, 并结合式 (26) 可得
由于α、β为正常数, p1、u1、p2、u2为正奇数, 故。于是由Lyapunov稳定性定理可知:系统渐近稳定。
假设系统误差的初始状态ei (0) >1, 则可将系统从初始状态收敛到达平衡点的过程分为两个阶段。
阶段1:从初始状态收敛到ei (t) =1的过程。此时e的收敛速度主要取决于式 (26) 中的改进线性滑模部分, 即
对式 (29) 两边取积分, 得完成该阶段所用的时间为
阶段2:从ei (t) =1收敛到平衡点的过程。此时e的收敛速度主要取决于式 (26) 中的终端滑模部分, 即
对式 (31) 两边取积分, 得到完成该阶段所用的时间为
注意到, 在求解t1和t2时都分别忽略了式 (26) 中的其中一项, 即采用了小于实际收敛速度的运动速度来计算到达时间, 因此, 系统从初始状态收敛到平衡点的总时间应该为
2.2 快变子系统
由于柔性关节和柔性臂都会引起系统的弹性振动, 而且振动级别不一定相同。因此, 我们考虑将快变子系统再次分解为两个子系统:描述柔性关节引起的系统弹性振动的“柔性关节-刚性臂”快变子系统和描述柔性臂引起的系统弹性振动的“刚性关节-柔性臂”快变子系统。对“柔性关节-刚性臂”快变子系统设计控制律τf1来抑制柔性关节引起的系统弹性振动;对“刚性关节-柔性臂”快变子系统设计控制律τf2来抑制柔性臂引起的系统弹性振动。于是, 快变子系统的总控制律可写为τf=τf1+τf2。
2.2.1“柔性关节-刚性臂”快变子系统
为了获得该子系统的动力学方程, 由上文的分析, 消除式 (14) ~式 (16) 中有关柔性臂的变量, 有
设计如下的基于转角速度差值的反馈控制律:
式中, K2为正定对角矩阵。
于是由式 (37) 可看出该控制律的控制原理如下:根据反馈回的电机转子和机械臂的转动角速度的差值来不断调节参数Kf, 从而保证系统的稳定性。
将式 (37) 代入式 (34) 可得到“柔性关节-刚性臂”快变子系统如下形式的动力学方程:
2.2.2“刚性关节-柔性臂”快变子系统
为了获得“刚性关节-柔性臂”快变子系统的动力学方程, 令系统动力学方程式 (14) ~式 (16) 中φ=q、, 消去系统中有关柔性关节的变量, 可得到
由式 (39) 可解得, 并将其代入式 (40) 、式 (41) , 整理后可得
引入快变时标tf=t/ε及边界层修正项。因为在快变系统中d于是结合式 (43) , 并令ε=0, 可得到“刚性关节-柔性臂”快变子系统如下形式的动力学方程:
由于式 (44) 为线性完全能控系统, 因此, 本文采用线性二次型最优控制器 (LQR) 来将系统状态ζ调节到零, 从而抑制柔性臂引起的系统弹性振动。若选取最优控制的性能泛函为 (Qf∈R4×4为对称正定常值矩阵, Rf∈R2×2为对称半正定常值矩阵) , 则可将LQR控制器设计为如下形式:
式中, P为Ricatti方程 (-PAf-AfTP+PBfRf-1BfTP-Qf=0) 的唯一解。
综上, 式 (21) 、式 (38) 和式 (44) 描述的即为漂浮基柔性关节-柔性臂空间机器人系统的奇异摄动模型。而式 (25) 、式 (37) 和式 (45) 分别为各子系统的控制律。
3 仿真算例
利用本文提出的慢变子系统的滑模控制律 (式 (25) ) 、快变子系统的速度差值反馈控制律 (式 (37) ) 及LQR控制器 (式 (45) ) 对图1所示的作平面运动的漂浮基柔性关节-柔性臂空间机器人系统进行数值仿真实验。系统的惯性参数为:l0=1.5m, l1=3m, l2=2.5m, d1=2m, m0=40kg, m1=2kg, J0=34.17kg·m2, J1=3kg·m2, Ja1=Ja2=0.53kg·m2, ρ=1kg/m, EI=300N·m2, Kφ=diag (200, 200) 。假设空间机器人的机械臂转角期望运动轨迹为qd1=0.5π (0.1t-0.5sin (0.2πt) /π) , qd2=0.5π (1-0.1t+0.5sin (0.2πt) /π) 。仿真初始值:q0 (0) =0, q (0) =0[.05 1.6]Trad, φ (0) =0[.05 1.6]Trad。仿真时间:t=10s。为了对比, 采用常规的线性滑模面与本文提出的非线性滑模面式 (22) 进行比较。仿真结果如图3~图8所示。其中, 图3为空间机器人的机械臂转角运动的轨迹跟踪对比图, 虚线表示采用基于线性滑模面的控制方法得到的空间机器人的机械臂转角qc的运动轨迹, 实线表示采用基于本文提出的非线性滑模面的控制方法得到的空间机器人的机械臂转角q的运动轨迹;点线表示空间机器人的机械臂期望转角qd的运动轨迹。图4所示为柔性关节引起的关节弹性变形偏差角。图5所示为关闭“柔性关节-刚性臂”快变系统控制律τf1后关节弹性变形偏差角。图6所示为柔性臂模态坐标变量δ的变化曲线。图7所示为柔性臂末端变形曲线。图8所示为关闭“刚性关节-柔性臂”快变子系统控制律τf2后柔性臂模态坐标变量δ的变化曲线。
从图3可看出, 本文所提出的混合控制方法能保证空间机器人的机械臂转角的运动精确且稳定地跟踪上期望运动轨迹, 保证了控制系统的精度和稳定性。而且与常规的基于线性滑模面的控制方法比较来看, 本文提出的基于非线性滑模面的控制方法的趋近速度得到了提高。从图4可看出, 柔性关节所引起的关节弹性变形偏差角σ虽然不为零, 但是被限制在一个非常小的范围内, 足以保证系统的控制精度。而从图5可看出, 当关闭τf1后关节弹性变形偏差迅速变大, 从而证明了τf1对抑制柔性关节引起的系统弹性振动的有效性。从图6和图7可看出, 柔性臂的振动得到了有效的抑制, 柔性臂末端的变形也很小。而从图8可看出, 当关闭τf2后柔性臂的二阶模态在3.5s的时候就变得很大, 进而仿真失效, 这说明此时柔性臂的振动无法得到抑制, 从而证明了τf2对抑制柔性臂的振动的有效性。
4 结束语
在考虑漂浮基空间机器人系统的机械臂和关节都存在柔性的情况下, 本文利用系统动量、动量矩守恒关系和拉格朗日-假设模态法建立了漂浮基柔性关节-柔性臂空间机器人的动力学模型, 并基于奇异摄动法提出了由控制系统运动的非线性滑模控制、抑制关节柔性振动的速度差值反馈控制和抑制臂柔性振动的LQR控制组成的混合控制律。仿真实验表明, 所提出的混合控制律能够补偿系统的关节转角弹性变形偏差, 保证空间机器人快速、精确、稳定地完成期望运动轨迹的渐近跟踪, 且能够有效地抑制柔性关节和柔性臂引起的系统弹性振动, 保证系统的稳定性, 体现了该混合控制律的良好控制品质。