• 首页关于本刊投稿须知期刊订阅编委会期刊合作查询检索English
动三轴试验滞回曲线的椭圆拟合分析
动三轴试验滞回曲线的椭圆拟合分析
毛远凤1*) 叶庆东1) 沈宇鹏2)
1)中国地震局第一监测中心,天津 300180;
2)北京交通大学土木工程学院,北京 100044
 [收稿日期]: 2017-03-06
摘要

采用正弦波形循环荷载动三轴试验获取土的动剪切模量和阻尼比的数据处理过程中通常存在两个困难:一是因为试验中存在各种噪声使滞回曲线椭圆形态不明显;二是椭圆拟合会因方法不当出现不收敛或者误差大等情况。为了在一定程度上克服第一个困难,本文将滤波技术引入到对应力应变时间序列的处理中来;为了更好地拟合滞回曲线椭圆,本文结合了主成份分析技术与椭圆的几何拟合方法,通过计算椭圆焦点位置、长半轴长度等来确定拟合椭圆,结果表明,该方法容易线性化且相对稳定,值得在动三轴数据处理中进行推广。



引言

土的动剪切模量比和阻尼比是土动力学特征的重要参数(Hardin等,1968),在工程场地土层地震反应分析和地震安全性评价工作中不可缺少(郭婷婷等,2016)。通过描述土在动力荷载作用下的动应力应变关系(动本构关系)模型,将动剪切模量和阻尼比与动剪应变幅值的函数关系具体化。Seed等(1970)首先给出了砂土和黏性土的动剪切模量和阻尼比与动剪应变幅值的关系曲线。Hardin等(1972)给出了动剪切模量与动剪应变的关系(即Hardin-Drnevich双曲线模型)。Martin等(1982)在Hardin-Drnevich模型的基础上进行了改进,提出了三参数的Davidenkov模型和具有幂次形式的阻尼比拟合公式。陈国兴等(2005)基于Davidenkov模型,采用上限剪应变幅值作为分界点,推导出了修正的Davidenkov模型阻尼比计算公式。在所有模型中,Hardin-Drnevich模型需要求解的参数最少,只需要2个参数,因而较为简便,目前应用较为广泛。下文中讨论动剪切模量和阻尼比与动剪应变幅值关系就基于此模型。

土的动剪切模量和阻尼特性测定是地震等动荷载作用下土工建筑物设计和计算的基本依据(赵红芬等,2009)。通常,为获得土体的动剪切模量和阻尼比,最直接的方法是进行土动力学试验。在土的动力学特性与岩土地震工程试验研究中最常用的方法主要是动三轴和自振柱试验两种方法。其中,动三轴试验主要用于在室内测定粘性土的动强度、各类土的动模量和阻尼比、砂土抗液化强度等相关参数(吴世明等,2001张涛,2004申权等,2013郭可骍,2015)。基于前人对土动力特性的分析,同时考虑到试验测试及数据采集系统的特征,试验荷载一般采用正弦波形循环荷载(南京水利科学研究院,2003)。在动三轴试验中,数据处理是保证试验成果有效产出的关键,而应力应变滞回曲线椭圆拟合则是判断试验成果准确性的决定因素。

1 循环应力动三轴试验

在循环应力动三轴试验中,土试样直径选取39.1mm,高取为80mm,对试样施加相等的固结压应力σ1c=σ3c=σ0,此时固结应力比Kc,为等压固结(Hardin等,1972张茹等,2006)。动力特性试验荷载激振频率定为1Hz,逐级施加动荷,每级振次为10次,在排水条件下进行振动试验。在动三轴试验中,测定动模量和阻尼比的方法是将轴向动应力作用于试样,应力幅值由小逐级增大进行试验,当试样变形超过允许值时,试验终止。将试验过程中采集到的各级荷载下的动应力应变滞回曲线进行整理,用以确定不同动应变时的动模量和阻尼比,如图 1所示。动弹性模量Ed为椭圆长轴斜率,阻尼比λd定义为循环消耗的能量与施加的最大应变能之比的1/4π倍(李松林,1990),根据应变能定义,一个循环内损耗的能量为滞回圈椭圆的面积πab,其中ab为椭圆的长短轴;施加的最大应变能为Pdεd/2,即图中三角形AOB的面积,Pdεd为椭圆与长轴交点中坐标值较大的一点坐标。


图 1 确定动弹性模量和阻尼比示意图 Fig. 1 Schematic diagram of calculating dynamic elastic modulus and damping ratio
2 拟合滞回曲线椭圆的方法分析

试验采用正弦波循环荷载,即动应力可以表示为:

$ {{P}_{\rm{a}}}={{P}_{\rm{max}}}\sin \left(\omega t+{{\varphi }_{1}} \right)+{{P}_{0}} $ (1)

式中,Pa表示动应力,Pmax为动应力最大变化幅度,P0为静态基准应力,ω为荷载振动的角频率,t为时间,${{\varphi }_{1}}$为动应力的初相。由于存在阻尼,动应变与动应力不再满足简单的线性关系,而是存在一个相位差。

对应的动应变可以表示为:

$ {{\varepsilon }_{\rm{a}}}={{\varepsilon }_{\rm{max}}}\sin \left(\omega t+{{\varphi }_{2}} \right)+{{\varepsilon }_{0}} $ (2)

式中,εa为动应变,εmax为动应变最大变化幅度,ε0为静态基准应变。公式(1)和(2)中,角频率ω与时间t是相同的,但因为相位差的存在使得${{\varphi }_{1}}\ne {{\varphi }_{2}}$,因此,动应力-应变关系不再是直线,而是在直角坐标系中的长短轴均与坐标轴斜交的椭圆。用于计算阻尼比λdPdεd可直接从应力-应变图上读取(图 4(b)),但椭圆的长短轴必须通过拟合得到,因此椭圆拟合成为获得阻尼比的关键。


图 4 椭圆拟合迭代前(a)后(b)变化示意图 Fig. 4 Schematic diagrams of ellipse fitting iteration, before(a) and after(b)

理论上,试验过程中应力及应变响应均为时间的函数,是一条光滑的单频谐波曲线,但由于误差的存在,实际得到的曲线除了设定频率的主频外,通常含有其它频率成份的噪声,如图 2(a)(c)所示。此时,可以通过设置滤波器来滤掉主频以外的成份。为了保证不发生相位的移动,应该采用无相移的滤波器;同时考虑振幅失真最小,我们选用通带内具有最大平坦特性的巴特沃斯滤波器。图 2(b)2(d)是采用4阶无相移巴特沃斯2Hz低通滤波器滤波后的结果,可见经滤波后曲线光滑了很多,振幅与滤波前变化不大。原则上,可以在图 2(b)(d)中直接读取初相${{\varphi }_{1}}$${{\varphi }_{2}}$,然后通过三角函数和差化积相关公式求得待拟合椭圆的相关参数,但是这样做误差较大,因此多数情况下仍采用椭圆拟合的方法来得到椭圆的相关参数。


图 2 应力、应变初始((a)、(b))及滤波后((c)、(d))的时间曲线 Fig. 2 Initial time series of stress (a) and strain (b), filtered time series of stress (c) and strain (d)

椭圆拟合的方法有很多种,如蒙特卡罗法、非线性最小二乘法(许正文等,2008闫蓓等,2008)、参数拟合法(陈基伟,2007a2007b)和几何法(彭青玉,2003)等。蒙特卡罗法也称统计模拟法,需要借助计算机生成大量样本并对其特征进行统计分析,通常较为费时;非线性最小二乘法直观明了,但是需要求解5个参数,同时还要求二元二次函数满足椭圆的约束条件,事实上是一个二次规划问题,容易因陷入局部极小而不能很好地拟合观测数据;参数拟合法处理已知初相的问题较为方便,但对本文涉及的问题无法线性化,求解存在困难。彭青玉(2003)根据椭圆到两定点距离之和为常数的几何定义出发,提出了椭圆拟合的几何方法。与传统方法相比,该方法直接基于几何定义,因此相对稳定,并且很容易线性化,可操作性强,下文中将对该方法进行详细介绍。

椭圆拟合的几何方法中,需要给定半长、短轴以及长轴与横坐标的夹角。我们的做法是首先将滤波后的应变-应力曲线画在一张图上,然后手动删去明显异常的点,如图 3


图 3 手动删除异常点前(a)后(b)应力应变曲线 Fig. 3 Stress-strain curves before (a) and after (b) deleted abnormity point manually

由于椭圆拟合问题本质上是非线性问题,初始值越接近于拟合值越容易收敛,计算也越快。我们首先通过二维主成份分析确定长轴的方向。主成份分析通过正交变换将原坐标系进行旋转得到一个新的坐标系,使变量在新坐标系下方差最大,通过确定旋转角度θ就可以得到椭圆长短轴的方向。关于主成份分析理论及计算过程可参考文献(杨淑莹等,2015)。尽管主成份分析的数学推导过程较为复杂,但是利用Matlab仅需要一条命令就可以实现。椭圆中心初始值取所有点的质心(εcPc),然后可以确定椭圆长轴所在的直线方程。由于长短轴是正交的,且都过椭圆的中心,因此椭圆的短轴方向也随之确定了。椭圆半长轴长度取所有点到短轴距离的最大值a,半短轴长度取所有点到长轴的距离最大值b。椭圆函数是二次函数,可以从实际椭圆的外部或内部逼近拟合,经尝试后发现从内部逼近较容易收敛。为此,取初始半长轴为所有点到短轴距离最大值的0.8倍。初始焦距长度为$c=\sqrt{{{a}^{2}}-{{b}^{2}}} $,两个焦点(ε1P1)和(ε2P2)的位置分别为:

$ \begin{matrix} {{\varepsilon }_{1}}={{\varepsilon }_{\rm{c}}}+c\cos \theta \\ {{P}_{1}}={{P}_{\rm{c}}}+c\sin \theta \\ {{\varepsilon }_{2}}={{\varepsilon }_{\rm{c}}}-c\cos \theta \\ {{P}_{2}}={{P}_{\rm{c}}}-c\sin \theta \\ \end{matrix} $ (3)

理想情况下椭圆上的点到两焦点距离之和为常数2a,实际中由于误差的存在,椭圆拟合方程表示为:

$ \sqrt{{{\left(\varepsilon -{{\varepsilon }_{1}} \right)}^{2}}+{{\left(P-{{P}_{1}} \right)}^{2}}}+\sqrt{{{\left(\varepsilon -{{\varepsilon }_{2}} \right)}^{2}}+{{\left(P-{{P}_{2}} \right)}^{2}}}=2a+\upsilon $ (4)

式中υ为残差。设:

$ \begin{matrix} {{l}_{10}}=\sqrt{{{\left(\varepsilon -{{\varepsilon }_{1}} \right)}^{2}}+{{\left(P-{{P}_{1}} \right)}^{2}}} \\ {{l}_{20}}=\sqrt{{{\left(\varepsilon -{{\varepsilon }_{2}} \right)}^{2}}+{{\left(P-{{P}_{2}} \right)}^{2}}} \\ \end{matrix} $ (5)

l10ε1P1的偏导数为:

$ \begin{matrix} \frac{\partial {{l}_{10}}}{\partial {{\varepsilon }_{1}}}=-\frac{\varepsilon -{{\varepsilon }_{1}}}{{{l}_{10}}} \\ \frac{\partial {{l}_{10}}}{\partial {{P}_{1}}}=-\frac{P-{{P}_{1}}}{{{l}_{10}}} \\ \end{matrix} $ (6)

l20ε2P2的偏导数具有与公式(6)相同的形式。将(4)式在初值处进行泰勒展开可得:

$ \frac{\varepsilon -{{\varepsilon }_{1}}}{{{l}_{10}}}\Delta {{\varepsilon }_{1}}+\frac{P-{{P}_{1}}}{{{l}_{10}}}\Delta {{P}_{1}}+\frac{\varepsilon -{{\varepsilon }_{2}}}{{{l}_{20}}}\Delta {{\varepsilon }_{2}}+\frac{P-{{P}_{2}}}{{{l}_{20}}}\Delta {{P}_{2}}+2\Delta a={{l}_{10}}+{{l}_{20}}-2a-\upsilon $ (7)

其中,Δε1、ΔP1,Δε2、ΔP2分别为两个焦点位置的改正量,Δa为半长轴的改正量。对于(7)式,每一个已知点对应于一个方程,若有N个试验值就可得到一个Gm=d形式的线性方程组,其中GN×5的系数矩阵,m为焦点位置与长轴改正量这5个系数的列矢量,d为由初始值构成的N元素列矢量,采用最小二乘可以求得5个改正量。实际中,将上一次的初始值加上在该初始值下得到的改正量,作为下一次计算的初始值进行新一轮的计算。通过迭代,不断更新焦点位置和长轴的长度,直至相邻两次的值之差或残差小于一定的阈值,一般5—10步就可以收敛。实际中,为了保证系数矩阵的非奇异,可以加一个微小的阻尼,即采用阻尼最小二乘法。图 4(b)中的拟合椭圆是以图 4(a)的初始椭圆迭代6次的结果。

3 试验结果

采用残余变形(孔压)较小的若干级循环,在每级循环中选取有代表性的滞回曲线,用本文的方法拟合椭圆,并计算动剪切模量和阻尼比(图 4(b))。然后,采用Hardin-Drnevich双曲线模型(Hardin等,1972)对这一系列的值进行拟合得到动剪切模量G与动剪应变γd、阻尼比λd的关系,最终将某一粘土样的试验结果绘制成G/Gmax-γdλd-γd的关系曲线,如图 5所示(Gmax为最大动剪切模量),不同动剪应变γd对应的G/Gmax值、λd值如表 1所示。


图 5 不同动剪应变下的G/Gmax及动阻尼比λd曲线 Fig. 5 Curves of G/Gmax and λd at different dynamic shear strain levels
表 1 试样动三轴试验结果,不同动剪应变下的G/Gmaxλd Table 1 The dynamical tri-axial test results of G/Gmax and λd at different shear strain levels

从表中可以看出,不同剪应变下,G/Gmax及动阻尼比λd值在粘性土合理值范围内(孙静等,2004),满足工程项目要求。

图 5可以看出,本文中椭圆拟合的几何方法能够很好地拟合曲线,是能够应用于工程实践的实用方法。

4 结论

通过上述试验和分析,得到了不同动剪应变下该粘土样的动剪切模量和阻尼比随剪应变变化的关系曲线。可以看出,本文进行椭圆拟合的几何方法是一种相对稳定、可操作性强的方法。在实际工程中,采用人机交互可有效地控制数据质量。将该方法应用于同一批其他试验土样的滞回曲线拟合,均得到了较好的结果,进一步证明了该方法可以确保土动三轴试验数据结果的可靠性以及各土动力学特性参数取值的准确性。

参考文献
陈国兴, 庄海洋, 2005. 基于Davidenkov骨架曲线的土体动力本构关系及其参数研究[J]. 岩土工程学报, 27(8): 860-864.
陈基伟, 2007a. 椭圆直接拟合算法研究[J]. 工程勘察, (6): 49-51.
陈基伟, 2007b. 工程测量中一类参数曲线的拟合[J]. 大地测量与地球动力学, 27(1): 100-103.
郭可骍, 2015. 非饱和黄土的动三轴试验研究. 西安: 长安大学.
郭婷婷, 秦梅梅, 2016. 土动三轴试验参数选取的理论分析与计算[J]. 工程抗震与加固改造, 38(2): 144-149.
李松林, 1990. 动三轴试验的原理与方法. 北京: 地质出版社.
南京水利科学研究院土工研究所, 2003. 土工试验技术手册. 北京: 人民交通出版社.
彭青玉, 2003. 木星土星边缘的椭圆拟合[J]. 云南天文台台刊, (4): 43-48.
申权, 李明俊, 蒋文明, 等, 2013. 动三轴试验测试土阻尼的影响因素与不足[J]. 江西科学, 31(1): 84-89.
孙静, 袁晓铭, 孙锐, 2004. 土动剪切模量和阻尼比的推荐值和规范值的合理性比较[J]. 地震工程与工程振动, 24(2): 125-133.
吴世明, 周健, 杨挺, 2001. 土动力学理论与计算. 北京: 中国建筑工业出版社.
许正文, 姚连璧, 2008. 基于稳健估计的直接最小二乘椭圆拟合[J]. 大地测量与地球动力学, 28(1): 77-80.
闫蓓, 王斌, 李媛, 2008. 基于最小二乘法的椭圆拟合改进算法[J]. 北京航空航天大学学报, 34(3): 295-298.
杨淑莹, 张桦, 2015. 模式识别与智能计算-MATLAB技术实现. 3版. 北京: 电子工业出版社.
张茹, 何昌荣, 费文平, 等, 2006. 固结应力比对土样动强度和动孔压发展规律的影响[J]. 岩土工程学报, 28(1): 101-105.
张涛, 2004. 中型动三轴仪的研制及城市垃圾土动力特性试验研究. 大连: 大连理工大学.
赵红芬, 何昌荣, 王莉娜, 2009. 动模量阻尼的动三轴试验研究[J]. 路基工程, (4): 158-160.
Hardin B. O., Black W. L., 1968. Vibration modulus of normally consolidated clay[J]. American Society of Civil Engineers, 94(2): 353-370.
Hardin B. O., Drnevich V. P., 1972. Shear modulus and damping in soils[J]. Journal of the Soil Mechanics and Foundations Division, 98(7): 667-692.
Martin P. P., Seed H. B., 1982. One-dimensional dynamic ground response analyses[J]. Journal of the Geotechnical Engineering Division, 108(7): 935-952.
Seed H. B., Idriss I. M., 1970. Soil moduli and damping factors for dynamic response analyses. Report No EERC70-10. Berkeley:Earthquake Engineering Research Center, University of California Berkeley.


Analysis on Ellipse Fitting of Hysteresis Curve in Dynamic Rriaxial Test
Mao Yuanfeng*1), Ye Qingdong1), Shen Yupeng2)
1) The First Monitoring and Application Center, China Earthquake Administration, Tianjin 300180, China;
2) School of Civil Engineering, Beijing Jiaotong University, Beijing 100044, China
Abstract

There exist two difficulties in the data processing of dynamic tri-axial tests when we use sinusoidal cyclic loading to obtain dynamic shear modulus and damping ratio. One difficulty is that the elliptic hysteresis curve is not clear enough due to various noises, and another one is that the procedure of fitting ellipse is divergent or the results of fitting with larger uncertainty because selected ellipse fitting method is improper. In this study, in order to overcome these difficulties to some extents, we processed time series of dynamical stress and stain with filtering technology. To fit the elliptic hysteresis curve well, we calculated the fitted ellipse by determining the two focuses and the length of the long axis of ellipse which we obtained by combining principal component analysis with ellipse geometrical fitting method. The result shows that this method is easy to linearize the nonlinear ellipse fitted problem and relatively stable, and it is worth popularizing in the dynamical tri-axial test data processing.



主办单位:中国地震台网中心
版权所有:《震灾防御技术》编辑部
地址:北京西城区三里河南横街5号,   邮编:100045
邮箱:zzfy2006@126.com   电话:59959251;59959462;59959408
访问人数:1624320
动三轴试验滞回曲线的椭圆拟合分析
毛远凤1*) 叶庆东1) 沈宇鹏2)
《震灾防御技术》, DOI:10.11899/zzfy20170402