基于DDA的混沌摆分析_于昊.pdf

第 41 卷第 12 期大学物理Vol41 No122022 年 12 月COLLEGEPHYSICSDec 2022收稿日期:20220506;修回日期:20220617基金项目:西南交通大学 2021 年本科教育教学研究与改革项目(2103107)资助作者简介:于昊(1998)，男，河北承德人，西南交通大学力学与航空航天学院 2021 级硕士研究生通信作者:喻勇，Email:yuyong2000 126com基于 DDA 的混沌摆分析于昊1，2，喻勇1，2(1 西南交通大学 力学与航空航天学院，四川 成都611756;2 西安交通大学 应用力学与结构安全重点实验室，四川 成都611756)摘要:采用非连续变形分析方法(Discontinuous Deformation Analysis，DDA)对混沌摆进行模拟分别建立了两自由度和四自由度无阻尼非线性自由振动的混沌摆模型，并对两种模型的运动微分方程进行推导，利用 Matlab 求解得到两种自由度下的摆角的理论解对原始 DDA 方程进行简化，推导得到刚体 DDA 方程，并利用 Matlab 编写了 DDA 程序，对两种摆的运动进行数值模拟，分别得到了两种模型下不同摆杆的角度时程图，揭示了混沌现象，并通过与理论结果进行对比，说明了该程序的适用性关键词:混沌摆;非连续变形分析;数值仿真中图分类号:O 4 1文献标志码:A文章编号:1000-0712(2022)12-0017-05【DOI】1016854/jcnki1000-0712220233混沌现象是 20 世纪一项重大的科学发现，至今仍然具有很高的研究意义混沌即在确定性系统中表现出来的随机的不规则运动，一个系统虽然可以用确定性理论来描述，但实际行为却表现出不确定性、不可重复、不可预测等性质，这就是混沌现象1 非线性混沌科学不仅在理论上具有重大意义，其在生态、医疗、金融以及决策等问题上也具有很大的价值，混沌对于非常复杂、系统性疾病的研究很有帮助，也可以进行气象特征的分析等2 此外，混沌还可用于信息加密，通过混沌系统得到的密码，具有了混沌的特征，使人们一时不能破解，也无法预测出密码的信息3 混沌摆对于非线性系统的力学行为演示更为直接，能把混沌中所蕴含的确定性和不确定性展示出来，对于混沌理论的理解很有帮助4 混沌摆还是力学教学中的热点目前来说，虽然在力学教学过程当中并没有规定进行混沌摆理论的介绍，但是向学生适当进行相关知识的介绍不仅会打破学生确定性的概念，也会增加其对力学学习的兴趣5 唐有绮等5 建立四自由度无阻尼混沌摆实验模型，进而建立其动力学方程，通过求解微分方程获得系统的摆角时程图，揭示了混沌现象朱桂萍等6 分别研究了混沌摆在保守系统和耗散系统下的动力学性质，得到了两种系统情形下的相图和角速度时序图，并分析了系统参数对系统的影响孟勇7 采用Maple 软件对大角度单摆、双摆、傅科摆的摆动进行了模拟，分析了每个摆的运动情况和物理规律李明达等8 建立混沌摆模型，得到系统的相图和分叉图，分析系统对外界参数的敏感性特征石根华9-10 提出的非连续变形分析方法(Dis-continuous Deformation Analysis，DDA)平行于有限单元法，其在分析非连续变形问题中，可以说是最有效的方法之一该方法以块体系统的位移作为未知量，从而建立系统的平衡方程，进行求解由于其计算比较高效，目前已经被广泛应用到模拟落石、滑坡、隧道等方面目前大部分对于混沌摆运动的揭示都是通过建立系统的拉格朗日方程，然后对方程求解，得到角度的解析解，然后得到系统的相图、频谱图或者时程图来进行分析虽然也可采用数值软件实现方程结果的可视化，但是前期对于方程的求解略微枯燥，对于不同的模型，可能求解出的结果非常复杂本文采用非连续变形分析方法，利用 Matlab 编写 DDA 的计算程序进行两种摆运动的模拟，不需要人工求解复杂的微分方程，并且可以得到运动的动画，使显示更加直观，更容易理解但是在 DDA 程序运行过程中，需要对程序的具体参数进行合理取值，如弹簧刚度和时间步长等，如果取值不当，将会导致程序不能运行或运行出错本18大学物理第 41 卷文对于程序中具体参数按照文献 1113 所讨论进行取值1DDA 方程以及模型的建立本节介绍 DDA 方程的构造方法以及两种模型的建立11DDA 方程的建立DDA 方程的建立关键是确定方程的子矩阵以及力矩阵的组成及其表达式，本节介绍 DDA 方程的建立依据以及子矩阵的确定111位移模式在本文中，只考虑二维块体，并且块体都为刚体，因此可以将原始 DDA 中的位移模式进行较大简化在原始 DDA 中块体的位移矩阵以及子矩阵都为 6 阶矩阵在考虑二维刚体情况下，矩阵都将变为3 阶矩阵，这将使得形式更加简单，计算更加方便块体内任一点(x，y)的位移 U=(u，v)T可通过位移函数确定，即U=TD(1)其中，T=10(yy0)01(xx0)()，D=u0v0r0()T，(x0，y0)为块体的质心坐标，(u0，v0)为块体沿 x 和 y方向的位移，r0为块体绕质心(x0，y0)转动的角度112平衡方程组多个块体之间是依靠约束相互连接的，因此块体间形成一联立方程组，以本文中的四块体系统为例k11k12k13k14k21k22k23k24k31k32k33k34k41k42k43k44D1D2D3D4=F1F2F3F4(2)对于二维刚性块体，系数矩阵kij(i，j=1，2，3，4)为 33 矩阵，Di，Fi(i=1，2，3，4)为 31 矩阵系数矩阵kij的具体数值由块体和块体间的连接方式确定113方程组子矩阵的确定首先由于系统中每个块体受到重力影响，由此形成体积荷载:0fyS0()T(3)其中 S 为块体面积，fy为体积荷载该子矩阵加入到总体方程中的 Fi 中去其次，块体在运动过程中会受到惯性力的影响，由此产生两个子矩阵:2M2 Ti(x，y)TTi(x，y)dxdy(4)2M Ti(x，y)TTi(x，y)dxdy Vi(0)(5)式(4)加入到总体方程 kii 中去，式(5)加入到总体方程 Fi 中，其中，M 为块体质量，为时间步长，Vi(0)为该时间步的初始速度，为上一时间步末的速度，具体计算公式为:Vi()=2 Di Vi(0)(6)最后，需要确定主摆和副摆之间的连接方式在原始 DDA 中，对块体可以添加点位移、点荷载和锚杆连接等方式主摆与复摆间的连接方式会对模拟结果产生重要的影响，不同的连接方式会产生不同子矩阵，通过多次实验，本文继续采用原始 DDA 中所用到的锚杆连接，锚杆连接方式会产生 4 个子矩阵采用锚杆将主摆和复摆短边中点相连接，由此会产生 4 个子矩阵:slEiEiT(7)slEiGjT(8)slGjEiT(9)slGjGjT(10)式(7)加入到总体方程的子矩阵 Kii 中去，式(8)加入到总体方程的子矩阵 Kij 中去，式(9)加入到总体方程的子矩阵 Kji 中去，式(10)加入到总体方程的子矩阵 Kjj 中去其中 s 是锚杆的刚度，l是锚杆的长度，i、j 表示两个块体 Ei 和 Gj 具体计算表达式:Ei=TiTlxly()(11)Gj=TjTlxly()(12)其中 lx=1l(x1x2)，ly=1l(y1y2)，x1、y1、x2、y2为锚杆连接的 2 个点的坐标综上所述，块体系统的子矩阵由以上 3 部分组成，即体积荷载、惯性力和锚杆连接产生的影响第 12 期于昊，等:基于 DDA 的混沌摆分析1912建立模型分别建立 2 自由度和 4 自由度混沌摆模型进行研究，本节简要介绍模型构成以及具体参数121两自由度混沌摆模型的建立建立如图 1 所示的 2 自由度无阻尼混沌摆模型，由一个 T 型主摆以及与其相连的一个副摆组成转动时，通过给 T 型主摆一角位移，使其绕中间一固定点 O 旋转T 型主摆可看作由 3 根等大的摆组成，其中每个摆的质量为 M，副摆的质量为 M，每根杆的长度均为 L图 1两自由度混沌摆模型122四自由度混沌摆模型的建立建立 4 自由度无阻尼混沌摆，如图 2 所示该装置由一个 T 型主摆以及所连接的 3 个副摆所组成该 T 型主摆绕中间一固定点 O 旋转，并通过锚杆与3 个副摆相连，这样四根杆组成四自由度无阻尼自由振动系统当进行实验时，通过给 T 型主摆一角位移，使整个系统进行运动T 型主摆可看作由 3 根等大的摆组成，其每个摆的质量为 M，副摆的质量为M，每根杆的长度均为 L图 2四自由度混沌摆模型13理论分析本节介绍 2 自由度和 4 自由度混沌摆理论结果的推导以及计算131两自由度混沌摆理论结果根据文献 2，对两自由度混沌摆运动微分方程进行如下推导系统动能为T=12(m1+m2)l2121+16m2l2222+12m2l1l212cos(12)(13)系统势能为V=12g(m1+2m2)l1cos 1+m2l2cos 2(14)其中，m1、m2分别为主摆和副摆的质量，l1、l2分别为主摆、副摆的长度在本文中，不考虑阻力的影响，因此将动能、势能代入到保守系拉格朗日方程中，对拉格朗日方程采用 4 阶 ungeKutta 法进行求解，利用 Matlab 编写了相关计算程序，其中使用了 ode45 函数，得到摆的运动微分方程组:1=3l1m2cos(12)sin(12)21+2l2m2sin(12)22+2gm1sin1+4gm2sin 13gm2cos(12)sin2/l1 4m1+4m23m2cos2(12)，2=3 gm1cos(12)sin 12gm2sin 22gm1sin 2+2gm2cos(12)sin 1+2 21l1m1sin(12)+2 21l1m2sin(12)+22l2m2cos(12)sin(12)/l2 4m1+4m23m2cos2(12)(15)132四自由度混沌摆理论结果根据文献 2，因文献 2公式输入错误，对其进行更改，得到正确的 4 自由度混沌摆运动微分方程推导如下系统动能为T=12(m1+3m2)l2121+16m2l22(22+23+24)+12m2l1l212cos(12)+3sin(13)4sin(14)(16)系统势能为20大学物理第 41 卷V=12g(m1+2m2)l1cos 1+m2l2(cos 2+cos 3+cos 4)(17)其中，m1、m2分别为主摆和副摆的质量，l1、l2分别为主摆、副摆的长度得到摆的运动微分方程组:1=g 3m2cos(123)3m2cos(124)3m2sin(122)(4m1+5m2)sin 1 3l1m2 sin 2(12)sin 2(13)sin 2(14)214l2m2 sin(12)22cos(13)23+cos(14)24 /8l1m1+3l1m2 5cos 2(12)+cos 2(13)+cos 2(14)，2=3 gsin 2l1sin(12)21+l1cos(12)12l2，3=3 gsin 3+l1cos(13)21+l1sin(13)12l2，4=3 gsin 4l1cos(14)21l1sin(14)12l2(18)2数值仿真在本文中，系统无阻尼且块体之间的连接都为光滑连接块体具体参数如下，质量 M=2103kg，长度 L=10 m，重力加速度 g=98 m/s2DDA 程序参数如下:锚杆的刚度 s 取为 51013N/m，锚杆的长度 l取为固定值 001 m，时间步长为 001 s通过给 T 型摆一初始角位移使整个系统运动，初始条件如下:1=4，2=0，3=0，4=0，1=0，2=0，3=0，4=021两自由度混沌摆仿真各杆无初速度释放，可以得到两自由度混沌摆各杆角度随时间变化的具体数值，并将该数值与理论解进行对比，得到如图 3、图 4 所示的角度对比图，可发现 DDA 计算的结果与理论结果相差不大，吻合性较好在运行 20 s 情况下，选取最下方块体左下角点，可以画出其轨迹曲线如图 5 所示，可以发现其运动为无规律的混沌现象图 3在初始条件 1下运动 10 s 的时程图图 4在初始