于 昊，喻 勇
(1. 西南交通大学 力学与航空航天学院，四川 成都 611756;2. 西安交通大学 应用力学与结构安全重点实验室，四川 成都 611756)
混沌现象是20世纪一项重大的科学发现，至今仍然具有很高的研究意义.混沌即在确定性系统中表现出来的随机的不规则运动，一个系统虽然可以用确定性理论来描述，但实际行为却表现出不确定性、不可重复、不可预测等性质，这就是混沌现象[1].非线性混沌科学不仅在理论上具有重大意义，其在生态、医疗、金融以及决策等问题上也具有很大的价值，混沌对于非常复杂、系统性疾病的研究很有帮助，也可以进行气象特征的分析等[2].此外，混沌还可用于信息加密，通过混沌系统得到的密码，具有了混沌的特征，使人们一时不能破解，也无法预测出密码的信息[3].混沌摆对于非线性系统的力学行为演示更为直接，能把混沌中所蕴含的确定性和不确定性展示出来，对于混沌理论的理解很有帮助[4].混沌摆还是力学教学中的热点.目前来说，虽然在力学教学过程当中并没有规定进行混沌摆理论的介绍，但是向学生适当进行相关知识的介绍不仅会打破学生确定性的概念，也会增加其对力学学习的兴趣[5].
唐有绮等[5]建立四自由度无阻尼混沌摆实验模型，进而建立其动力学方程，通过求解微分方程获得系统的摆角时程图，揭示了混沌现象.朱桂萍等[6]分别研究了混沌摆在保守系统和耗散系统下的动力学性质，得到了两种系统情形下的相图和角速度时序图，并分析了系统参数对系统的影响.孟勇[7]采用Maple软件对大角度单摆、双摆、傅科摆的摆动进行了模拟，分析了每个摆的运动情况和物理规律.李明达等[8]建立混沌摆模型，得到系统的相图和分叉图，分析系统对外界参数的敏感性特征.
石根华[9-10]提出的非连续变形分析方法(Discontinuous Deformation Analysis，DDA)平行于有限单元法，其在分析非连续变形问题中，可以说是最有效的方法之一.该方法以块体系统的位移作为未知量，从而建立系统的平衡方程，进行求解.由于其计算比较高效，目前已经被广泛应用到模拟落石、滑坡、隧道等方面.
目前大部分对于混沌摆运动的揭示都是通过建立系统的拉格朗日方程，然后对方程求解，得到角度的解析解，然后得到系统的相图、频谱图或者时程图来进行分析.虽然也可采用数值软件实现方程结果的可视化，但是前期对于方程的求解略微枯燥，对于不同的模型，可能求解出的结果非常复杂.本文采用非连续变形分析方法，利用Matlab编写DDA的计算程序进行两种摆运动的模拟，不需要人工求解复杂的微分方程，并且可以得到运动的动画，使显示更加直观，更容易理解.
但是在DDA程序运行过程中，需要对程序的具体参数进行合理取值，如弹簧刚度和时间步长等，如果取值不当，将会导致程序不能运行或运行出错.本文对于程序中具体参数按照文献[11-13]所讨论进行取值.
本节介绍DDA方程的构造方法以及两种模型的建立.
1.1 DDA方程的建立
DDA方程的建立关键是确定方程的子矩阵以及力矩阵的组成及其表达式，本节介绍DDA方程的建立依据以及子矩阵的确定.
1.1.1 位移模式
在本文中，只考虑二维块体，并且块体都为刚体，因此可以将原始DDA中的位移模式进行较大简化.
在原始DDA中块体的位移矩阵以及子矩阵都为6阶矩阵.在考虑二维刚体情况下，矩阵都将变为3阶矩阵，这将使得形式更加简单，计算更加方便.
块体内任一点(x,y)的位移U=(u,v)T可通过位移函数确定，即
U=T·D
(1)
1.1.2 平衡方程组
多个块体之间是依靠约束相互连接的，因此块体间形成一联立方程组，以本文中的四块体系统为例
(2)
对于二维刚性块体，系数矩阵kij(i,j=1,2,3,4)为3×3矩阵，Di,Fi(i=1,2,3,4)为3×1矩阵.系数矩阵kij的具体数值由块体和块体间的连接方式确定.
1.1.3 方程组子矩阵的确定
首先由于系统中每个块体受到重力影响，由此形成体积荷载：
(3)
其中S为块体面积，fy为体积荷载.该子矩阵加入到总体方程中的[Fi]中去.其次，块体在运动过程中会受到惯性力的影响，由此产生两个子矩阵：
(4)
(5)
式(4)加入到总体方程[kii]中去，式(5)加入到总体方程[Fi]中，其中，M为块体质量，Δ为时间步长，Vi(0)为该时间步的初始速度，为上一时间步末的速度，具体计算公式为：
(6)
最后，需要确定主摆和副摆之间的连接方式.在原始DDA中，对块体可以添加点位移、点荷载和锚杆连接等方式.主摆与复摆间的连接方式会对模拟结果产生重要的影响，不同的连接方式会产生不同子矩阵，通过多次实验，本文继续采用原始DDA中所用到的锚杆连接，锚杆连接方式会产生4个子矩阵.
采用锚杆将主摆和复摆短边中点相连接，由此会产生4个子矩阵：
(7)
(8)
(9)
(10)
式(7)加入到总体方程的子矩阵[Kii]中去，式(8)加入到总体方程的子矩阵[Kij]中去，式(9)加入到总体方程的子矩阵[Kji]中去，式(10)加入到总体方程的子矩阵[Kjj]中去.其中s是锚杆的刚度，l是锚杆的长度，i、j表示两个块体.[Ei]和[Gj]具体计算表达式：
(11)
(12)
1.2 建立模型
分别建立2自由度和4自由度混沌摆模型进行研究，本节简要介绍模型构成以及具体参数.
1.2.1 两自由度混沌摆模型的建立
建立如图1所示的2自由度无阻尼混沌摆模型，由一个T型主摆以及与其相连的一个副摆组成.转动时，通过给T型主摆一角位移，使其绕中间一固定点O旋转.T型主摆可看作由3根等大的摆组成，其中每个摆的质量为M，副摆的质量为M，每根杆的长度均为L.
图1 两自由度混沌摆模型
1.2.2 四自由度混沌摆模型的建立
建立4自由度无阻尼混沌摆，如图2所示.该装置由一个T型主摆以及所连接的3个副摆所组成.该T型主摆绕中间一固定点O旋转，并通过锚杆与3个副摆相连，这样四根杆组成四自由度无阻尼自由振动系统.当进行实验时，通过给T型主摆一角位移，使整个系统进行运动.T型主摆可看作由3根等大的摆组成，其每个摆的质量为M，副摆的质量为M，每根杆的长度均为L.
图2 四自由度混沌摆模型
1.3 理论分析
本节介绍2自由度和4自由度混沌摆理论结果的推导以及计算.
1.3.1 两自由度混沌摆理论结果
根据文献[2]，对两自由度混沌摆运动微分方程进行如下推导.
系统动能为
(13)
系统势能为
(14)
其中，m1、m2分别为主摆和副摆的质量，l1、l2分别为主摆、副摆的长度.
在本文中，不考虑阻力的影响，因此将动能、势能代入到保守系拉格朗日方程中，对拉格朗日方程采用4阶Runge-Kutta法进行求解，利用Matlab编写了相关计算程序，其中使用了ode45函数，得到摆的运动微分方程组：
4gm2sinθ1-3gm2cos(θ1-θ2)sinθ2]/
{l1[4m1+4m2-3m2cos2(θ1-θ2)]}，
2gm1sinθ2+2gm2cos(θ1-θ2)sinθ1+
{l2[4m1+4m2-3m2cos2(θ1-θ2)]}
(15)
1.3.2 四自由度混沌摆理论结果
根据文献[2]，因文献[2]公式输入错误，对其进行更改，得到正确的4自由度混沌摆运动微分方程推导如下.
系统动能为
(16)
系统势能为
m2l2(cosθ2+cosθ3+cosθ4)]
(17)
其中，m1、m2分别为主摆和副摆的质量，l1、l2分别为主摆、副摆的长度.
得到摆的运动微分方程组：
3m2sin(θ1-2θ2)-(4m1+5m2)sinθ1]-
3l1m2[sin 2(θ1-θ2)-sin 2(θ1-θ3)-
{8l1m1+3l1m2[5-cos 2(θ1-θ2)+
cos 2(θ1-θ3)+cos 2(θ1-θ4)]}，
(18)
在本文中，系统无阻尼且块体之间的连接都为光滑连接.块体具体参数如下，质量M=2×103kg，长度L=10 m，重力加速度g=9.8 m/s2.DDA程序参数如下：锚杆的刚度s取为5×1013N/m，锚杆的长度l取为固定值0.01 m ,时间步长为0.01 s.通过给T型摆一初始角位移使整个系统运动，初始条件如下：
θ′1=0，θ′2=0，θ′3=0，θ′4=0.
2.1 两自由度混沌摆仿真
各杆无初速度释放，可以得到两自由度混沌摆各杆角度随时间变化的具体数值，并将该数值与理论解进行对比，得到如图3、图4所示的角度对比图，可发现DDA计算的结果与理论结果相差不大，吻合性较好.
图3 在初始条件θ1下运动10 s的时程图
图4 在初始条件θ1下运动20 s的时程图
在运行20 s情况下，选取最下方块体左下角点，可以画出其轨迹曲线如图5所示，可以发现其运动为无规律的混沌现象.
图5 两自由度混沌摆一点轨迹图
2.2 四自由度混沌摆仿真
各杆无初速度释放，可以得到4自由度混沌摆各杆角度随时间变化的具体数值，并将该数值与理论解进行对比，得到如图6、图7所示的角度对比图，可发现DDA计算的结果与理论微分方程计算的结果相差不大，吻合性较好.
图6 在初始条件θ1下运动10 s的时程图
图7 在初始条件θ1下运动20 s的时程图
在运行20 s情况下，选取最下方块体左下角点，可以画出其轨迹曲线如图8所示，可以发现其运动为无规律的混沌现象.
图8 四自由度混沌摆一点轨迹图
从图8中可以发现，当初始条件较大时，系统即呈现出非周期运动的特点，杆的运动将变得不规律，数值模拟出的结果与理论解差别很小，说明了该DDA程序的适用性，也说明了对于块体间连接方式选择锚杆连接的正确性.
本文不同于以往用微分方程直接研究混沌摆的运动，采用非连续变形分析方法，对2自由度和4自由度混沌摆分别进行模拟，得到以下结论：
1) 运用Matlab软件编写DDA程序，分别将模型运行10 s和20 s，并与理论结果进行对比，发现角度变化曲线吻合较好，说明了该DDA程序的适用性；
2) 运用DDA方法所给出的锚杆连接方式，可以近似代替混沌摆中主摆与复摆之间的铆钉连接，并且结果相差不大；
3) 随着时间的延长，DDA计算的曲线与理论解得到的曲线仍会存在一定误差，是因为在每一时间步内加速度为常数，因此随着计算时间的增加，引起误差的积累，从而导致曲线出现偏差. 此外，对于连接块体的锚杆长度可能也会对结果产生影响.
猜你喜欢时程块体锚杆喷淋装置在锚杆钢剪切生产中的应用山东冶金(2022年1期)2022-04-19一种新型单层人工块体Crablock 的工程应用水运工程(2020年11期)2020-11-27隧洞块体破坏过程及稳定评价的数值方法研究湖南大学学报（自然科学版）(2020年5期)2020-06-03考虑增量时程贡献趋向和误差排序的多阻尼目标反应谱拟合*地震研究(2020年1期)2020-05-01模拟汶川地震动持时的空间分布规律研究地震研究(2019年4期)2019-12-19锚杆钢筋质量提升生产实践山东冶金(2019年1期)2019-03-30剂量水平与给药时程对豆腐果苷大鼠体内药代动力学的影响天然产物研究与开发(2018年10期)2018-11-06结构面对硐室稳定性的影响资源环境与工程(2016年3期)2016-06-09建筑施工中的锚杆静压桩技术中国房地产业(2016年2期)2016-03-01复合盾构在纵向锚杆区的掘进分析及实践工程建设与设计(2016年4期)2016-02-27
