（应用数学专业论文）弹性杆拟动力学模型和动力学模型的研究.pdf

摘要 弹性杆是一种重要的力学模型，许多工程构件和生物体如海底电缆、纤维、生 物和有机物大分子等，在一定条件下都可以模型化为弹性细杆讨论。近4 0 年来， 随着生物和遗传工程的发展，人们发现：d n a 是具有弹性特性的长链，可以利用 弹性细杆模型进行动力学分析。因此，人们将弹性力学的方法和分子生物学的知识 相结合对d n a 进行了研究，取得了大量成果。 本文研究了弹性杆的数学建模和数值仿真问题。内容包括弹性杆的曲面模型、 k i r c h h o f f 弹性杆的拟h a m i l t o n 方程及其辛算法、弹性杆动力学方程解的存在性及 数值模拟，具体工作为： ( 1 ) 基于弹性杆的k i r c h h o f f 假设并使用k i r c h h o f f 比拟，从新的角度提出了超细长 弹性杆曲面模型的微分方程组。 ( 2 ) 引入欧拉参数给出了超细长弹性杆的k i r c h h o f f 模型的拟动力学方程组。现有 的结果中，这种模型是用e u l e r 角y ，只缈) 为变量描述的。由于当0 = k g 时e u l e r 角 出现奇点，给数值计算造成困难。我们利用欧拉参数代替e u l e r 角，建立了以欧拉 参数为变量的描述超细长弹性杆的拟h a m i l t o n 方程，并采用辛算法建立了弹性杆长 距离保结构计算的数值计算模式并给出了数值结果。 ( 3 ) 根据l a x 方程解的存在性，证明了弹性杆动力学方程组解的存在性；在给出了 二次型形式的弹性势能函数的情况下，给出动态弹性杆在一段时间内的数值模拟。 关键词：k i r e h h o f f 弹性杆，欧拉参数，拟动力学，动力学方程，数值模拟 a b s t r a c t e l a s t i cr o d ，a l li m p o r t a n tc l a s s i c a lm e c h a n i c a lm o d e lw i d e l yu s e di nm a n y e n g i n e e r i n gs y s t e m ss u c ha s u n d e r s e ac a b l e s ，f i b e r sa n db i o l o g i c a lt i s s u e si n2 0 t h c e n t u r y , h a v eb e e na t t r a c t e dg r e a ta t t e n t i o nd u r i n gt h el a s t4 0y e a r sa sau s e f u lm o d e li n d e s c r i b i n gd n a a n dl o t so fa c h i e v e m e n t so nt h i st o p i ch a v eb e e no b t a i n e dr e c e n t l y i n t h i sp a p e r , m a t h e m a t i c a lm o d e l i n ga n dn u m e r i c a ls i m u l a t i o nm e t h o do ne l a s t i c r o d sa r es t u d i e d ，i n c l u d i n gt h es u r f a c em o d e lo fa ne l a s t i cr o d ，t h eq u a s i - h a m i l t o n e q u a t i o n sa n ds y m p l e c t i cm e t h o do fk i r c h h o f fr o d ，a n dt h ee x i s t e n c eo fs o l u t i o no ft h e e q u a t i o n sa n di t sn u m e r i c a ls t i m u l a t i o n t h em a i nr e s u l t si n c l u d e s ： ( 1 ) b a s e do nt h ea s s u m p t i o no fak i r c h h o f fr o da n du s i n gk i r c h h o f fa n a l o g u et e c h n i q u e ， ad i f f e r e n t i a le q u a t i o ns y s t e mi sr a i s e dt od e s c r i b et h es u r f a c eo fas u p e rs l e n d e re l a s t i c r o d ( 2 ) q u a t e r n i o ni si n t r o d u c e dt os e tu pt h eq u a s i - d y n a m i c a lm o d e lo fas u p e r - s l e n d e r e l a s t i cr o d p r e c e d i n gr e s u l t s ，t h o u g hr a i s e dl o n ga g o ，w e r ed e s c r i b e db ye u l e ra n g l e s ( 奶口，缈) t h es i n g u l a r i t yo fe u l e ra n g l e sa t 口= k n m a k e i t sn u m e r i c a lt r e a t m e n tv e r y d i f f i c u l t i nt h i sp a p e r , b yi n t r o d u c i n gq u a t e r n i o nt or e p l a c et h ee u l e ra n g l e si nt h e q u a s i d y n a m i c a lm o d e l s ，w ed e r i v e dt h eq u a s i h a mil t o ne q u a t i o n so fk i r c h h o f fr o d s o s y m p l e c t i cm e t h o dc a nb ei n t r o d u c e di ni t sn u m e r i c a lc o m p u t a t i o n ( 3 ) a c c o r d i n gt ot h ee x i s t e n c eo fs o l u t i o no fl a xe q u a t i o n s ，w ep r o v e dt h ee x i s t e n c eo f s o l u t i o no fd y n a m i cr o de q u a t i o n s u n d e rt h ec o n d i t i o nt h a tt h ee l a s t i ce n e r g yf u n c t i o n h o l d sq u a d r a t i cf o r m ，t h en u m e r i c a ls t i m u l a t i o no fd y n a m i cr o di sg i v e n k e y w o r d s ：k i r c h h o f fe l a s t i cr o d ；q u a t e r n i o n ；q u a s i d y n a m i c s ；d y n a m i ce q u a t i o n s ； n u m e r i c a is i m u l a t i o n 学位论文独创性声明、学位论文知识产权权属声明 学位论文独创性声明 本人声明，所呈交的学位论文系本人在导师指导下独立完成的研究成果。文中 依法引用他人的成果，均已做出明确标注或得到许可。论文内容未包含法律意义上 已属于他人的任何形式的研究成果，也不包含本人已用于其他学位申请的论文或成 果。 本人如违反上述声明，愿意承担由此引发的一切责任和后果。 做储躲酱缸 眺- p f 月日 学位论文知识产权权属声明 本人在导师指导下所完成的学位论文及相关的职务作品，知识产权归属学校。 学校享有以任何方式发表、复制、公开阅览、借阅以及申请专利等权利。本人离校 后发表或使用学位论文或与该论文直接相关的学术论文或成果时，署名单位仍然为 青岛大学。 本学位论文属于： 保密口，在年解密后适用于本声明。 不保密旷 ( 请在以上方框内打“”) 论文作者签名：鑫1 久 h 期： 导师签磊拘力 日期 ( 本声明的版权归青岛大学所有，未经许可，任何单位及任何个人不得擅自使用) 扩 年 啤 7峥 下可 引言 引言 本文的目的是研究弹性杆曲面方程、拟动力学模型和动力学学方程及其数值计 算和数值模拟的问题。 弹性杆是一种重要的力学模型，工程上的许多研究对象如纤维、海底电缆【l 】和 钻杆【2 】等，在一定条件下都可以模型化为弹性细杆讨论。弹性杆的研究是一个经典 的问题。早在1 7 3 0 年，伯努利和欧拉就已经开始研究弹性细杆在外力和力矩下的 变形问题。弹性细杆的一个重要而简单的经典力学模型是1 0 0 多年前k i r c h h o f f 提 出的基于刚性截面假设得到的模型，称为k i r c h h o f f 弹性杆【3 1 1 4 l 【5 1 【6 1 。 遗传科学的发展给弹性细杆的研究注入了新的活力。d n a 是遗传信息的携带 者，对d n a 和d n a 缠绕组织蛋白所构成的高等细胞染色体结构和功能的研究是揭 开遗传之谜的关键环节之一。19 5 3 年，w a t s o n 和c r i c k 证明了d n a 是双螺旋链状 结构的人分子；1 9 世纪6 0 年代和7 0 年代的进一步研究表明：d n a 是直径约为2 纳米的双螺旋结构的聚合长链分子，在不同的条件下通过缠绕组织蛋白卷曲成1 0 纳米或3 0 纳米的超螺旋结构。d n a 的这些特殊结构在遗传特性的建立和调控上有 着重大的影响。实验表明，d n a 的杆状螺旋结构和对组织蛋白的缠绕呈现出类似 弹性杆的力学特征1 7 1 1 8 1 ；同时，其分子又表现出微观的运动规律，如：在液体中漂 流等。因此人们把弹性力学的方法和分子生物学的实验和统计分析技术相结合1 9 1 对 其进行了研究，取得了大量成果并对遗传科学的发展具有实际的指导意义。通过对 d n a 弹性动力学模型的研究求解，可以给出d n a 在液体中扭结、缠绕等运动的数 值模拟，由此体现出d n a 分子在细胞核液体中的各种运动状态。这项研究用以理 解在基因激活、r n a 转录、d n a 修复和d n a 复制的过程中与染色体结构相关的 结构变化，为生物工程科技工作者关于d n a 的研究提供物理上的理论支持。 近3 0 以来，弹性细杆的k i r c h h o f f 模型在生物大分予的研究中得到新的应用【10 1 。 实验表明，d n a 长链的长度为纳米级的小段不容易变彤，表现出较强的刚性i 】【1 2 】； 而整根d n a 长链的长度和直径的比达到1 0 7 量级，在细胞核中的白山运动表现出 柔软的一面，可以任意弯曲和扭曲变形以容纳在狭小的细胞核内。这些实验现象和 k i r c h h o f f 弹性杆的一些假设一致。据此，b e n h a m t l 3 1 和l eb r e t 1 4 1 等在k i r c h h o f f 弹 性杆理论基础上建立了d n a 的弹性杆模型；s h i 和h e a r s t 【1 5 1 利用曲率r 、挠率彳 和扭角z 描述k i m h h o f f 圆截面弹性杆，给出了描述d n a 超螺旋结构的s c h r o d i n g e r 青岛大学硕士学位论文 方程。经典的k i r c h h o f f 模型一搬利用e u l e r 角描述，不利于数值计算，因此， k e h r b a u m 【1 7 】和h u 1 8 】引入e u l e r 参数代替e u l e r 角给出弹性杆轴心曲线的数学描述。 k i r c h h o f f 弹性杆模型的一个重要特点是其平衡方程和刚体定点转动微分方程 的相似性。k i r c h h o f f 弹性杆的自变量是弧长s ，把弧长换成时问t ，则其平衡方程和 某个刚体定点转动的动力学方程一致【1 9 1 2 0 1 。因此，刚体运动的动力学性质和弹性杆 的结构和平衡特性相对应，称之为k i r c h h o f f 动力学比拟理论【2 。有关的比拟关系 在n i z z e t e 和g o r i e l y t 2 2 】的论文中得到详尽的描述。根据k i r c h h o f f l ： ：拟原理，k i r c h h o f f 弹性杆的- 甲衡方程可以表述成类似动力学方程的结构，我们称之为拟动力学方程， 如拟l a g r a n g e 方程等。 利用拟动力学方程研究弹性杆，可以充分利用动力学的已有结果【2 3 】，如对称性、 守恒性【2 4 】【2 5 1 和稳定性等。目前的动力学比拟是利用e u l e r 角( 口，纠描述的【2 6 1 ，由 于在p = k n 处出现奇点，给数值计算带来了困难。而d n a 弹性杆模型的超细长、 大变形对数值仿真提出新的要求。在动力学问题中引入时间变量后，弹性杆的动力 学模型成为以弧坐标s 和时间，为双变量的动力系统。s i m o 和m a r d s o n 2 7 对该动力 学系统有了初步的论述；刘延柱建立了k i r c h h o f f 弹性杆的动力学模型1 2 3 。 出于d n a 弹性杆模型的数值计算的需要，本文给出以下结果： ( 1 ) 由于k i r c h h o f f 弹性杆是刚性截面，可以利用截面的刚体运动形成弹性杆表 面方程。我们导出了更为简单的弹性杆表面微分方程组，来作为d n a 自接触问题 的约束方程组。 ( 2 ) 引入e u l e r 参数代替e u l e r 角作为状态变量，导出了描述k i r c h h o f f 弹性杆的 拟l a g r a n g e 方程和拟h a m i l t o n 方程。与e u l e r 角描述的拟动力学方程截然不同的是， 本文导出的方程组是奇异的微分代数方程组。文中讨论了方程组的指标和解的存在 唯一性问题。 ( 3 ) 由于利用e u l e r 参数描述的拟l a g r a n g e 方程和拟h a m i l t o n 方程与传统的 e u l e r 角模型比较更有利于数值计算，作为上述模型的一个应用，本文利用所导出 的模型给出了计算k i r c h h o f f 弹性杆的保结构算法。与利用e u l e r 角描述的相应算法 比较，计算不受奇点的影响。另外，由于方程组已经不是微分方程组，因此，我们 对计算过程作了相应的处理。 ( 4 ) 注意到辛算法保持动力系统的运动对称性，由k i r c h h o f f 比拟原理，k i r c h h o f f 弹性杆可以比拟成刚性截面沿弹性杆轴线作匀速运动和转动形成。因此，与运动对 2 引言 称性相对应的是弹性杆的沿轴线方向的结构对称性。因此，上述保结构算法的物理 意义是在长距离计算中保持弹性杆力学模型的结构对称性。 ( 5 ) 给出了弹性杆动力学方程解的存在性。并对一类动力学方程进行了数值求 解和数值模拟，并形成了完整的数值计算程序。 本文的内容共分为三章：第一章主要介绍本文使用的有关符号、概念和相关弹 性杆基础知识和理论。第二章给出了基于弹性杆的k i r c h h o f f 假设，运用k i r c h h o f f 动力学比拟，从一个新的角度，导出了描述刚性截面弹性杆的曲面模型及数值算例； 基于超细长弹性杆的计算考虑，利用e u l e r 参数代替传统的e u l e r 角的描述方法，建 立了描述弹性杆平衡的拟h a m i l t o n 方程；采用辛算法建立了弹性杆长距离保结构计 算的数值计算模式并给出了数值结果。第三章主要给出了弹性杆动力学方程解的存 在性，并给出了数值模拟。 第一章弹性杆的k i r c h h o f f 方程及其四元数的引入 第一章弹性杆的k i r c h h o f f 方程及其四元数的引入 弹性杆的一个重要而简单的经典力学模型是1 0 0 多年前k i r c h h o f f 提出的基于 刚性截面假设得到的模型，称为k i r c h h o f f 弹性杆。一百多年来，人们对这种模型 进行了大量的研究，我们可以在许多教科书中看到这些工作的总结【2 3 】。 1 1 弹性杆的截面主轴坐标系 以空问中的一点o y g n a ，建立惯性坐标系( o - 勿f ) 。讨论空间中的一条光 滑曲线，( s ) ，其中s 为弧坐标。在曲线上任取一点p ，则曲线在点尸处的单位切向 量为 r ( s ) ：掣 1 - ( 1 ) 曲率为 小，= 警l 单位法向量和单位副法向量可以表示为 呻) = 南掣叫萨那) 单位向量n ，b ，t 组成以p 为原点，依附在曲线，( s ) 的右手坐标系， f r e n e t 坐标系。定义曲线在点p 处的挠率为 彳( s ) - 1 - ( 2 ) 称为曲线的 引入矢量 婢( s ) = r ( s ) b ( s ) + 彳( s ) 丁( s ) l 一( 3 ) 称为曲线的d a r b o u x 矢量，其物理意义是：当p 点沿曲线r ( s ) 以单位速度向弧坐标 s 的j 下向运动时，f r e n e t 坐标系相对惯性坐标系( d 一勃f ) 的转动速度。 单位向量，b ，r 关于j 的变化规律取决于f ( s ) ，由以下常微分方程组确定： 4 青岛大学硕士学位论文 _ d n = t o p x ：r b g f l s _ d b = t o p x b ：呵n g s i d t = 畔r = r n l - ( 4 ) 称为f r e n e t - s e r r e t 方程组。给定曲率r ( s ) 、挠率f ( s ) 和n ，b ，t 的适当初值条件时， 可从f r e n e t - s e r r e t 方程组解出矢量，刀，r 的变化规律。根据l - ( 1 ) 式，空间曲线，( j ) 可通过以下积分得出： ，( j ) = 【：r ( 仃) d 仃+ ，( o ) l 一( 5 ) 其中，( 0 ) 表示初始点的位置。 讨论一细杆，杆截面的几何中心连成的空间曲线，( s ) 称为曲杆的中心线。假设 杆的横截面为刚性截面；忽略弯曲引起的剪切变形，横截面与中心线正交；忽略中 心线的拉仲变形，任意两截面沿中心线的距离不变。 以中心线上的任意一点p 为原点建立截面主轴坐标系( p x y z ) ，各坐标轴的基 矢量分别为面，畋，以，其中z 轴与曲线，( s ) 的切线轴丁重合，即以= t 。设z 为x 轴 与轴的夹角，即截面相对f r e n e t 坐标系扭转的角度如图1 1 所示，则截面主轴坐 标系( 尸一舻) 和f r e n e t 坐标系( p n b t ) 之间的方向余弦如表1 1 所示 图1 1 截面主轴坐标系与f r e n e t 坐标系 5 第一章弹性杆的k i r c h h o f f 方程及其四元数的引入 b7 工 c o s 瓮 s i l l z 0 y - - s i n z c o sz o z 00 l 爱1 1 设当p 点沿曲线r ( s ) 以单位速度向弧坐标s 的正向运动时，杆的横截面以角速 度相对( d 一勿f ) 转动。为横截面相对惯性坐标系( o 一勃f ) 的绝对角速度，即 横截面相对f r e n e t 坐标系( p n b t ) 的相对角速度( d z d s ) 嵋- 与( p - n b t ) 相对惯性 坐标系( 0 一勿f ) 的牵连角速度坼( s ) 之和，即 国= ( 警) 以+ 嘶刮m 时坞以 - ( 6 ) 其中 q = r e s i nz ，锡= c c o s , 3 y ，妈= f + 譬1 ( 7 ) 国是与杆的弯扭变形和扭转变形相关的欠量，我们称之为弯扭度。 1 2k i r c h h o f f 方程1 2 8 1 1 2 9 l 1 2 1 弹性杆平衡方程 考虑一细长弹性杆，杆的起点为只。假设该弹性杆为等截面，且截面内两主轴 方向有相同的几何尺度；杆为均匀且各向同性，弹性系数为常数；忽略杆的体积力 以及杆与杆之间的分布力。设尸为杆中心线上任意一点，为的邻近点，尸和尸，点 相对惯性参考系原点o 的矢径分别为，和，+ ，相对只点的弧坐标分别为s 和 s + 厶。规定p 和尸，点处外法线矢量与弧坐标增大方向一致的截面为正截面。分析 即7 微元弧段内杆的平衡情况，设尸点的负截面受邻近截面作用的内力主矢和内力 主矩分别为一f 和一膨，尸点的正截面受邻近截面作用的内力主矢和内力主矩分别 为，+ 心和m + ，。在平衡状态下，作用力对点p 的主矢之和、主矩之和都为零， 即有 a f = 0 ，a m + a r x f = 0 6 l - ( 8 ) 青岛大学硕士学位论文 将l - ( 8 ) 式的各项郡除以血，并且令厶一0 ，司得惯性坐标系中的弹性杆的平衡方 程 _ d f ：口，掣+ r f ：口 1 ( 9 ) ds 利用惯性坐标系和截面主轴坐标系的关系，将平衡方程1 ( 9 ) 改写成在截面主轴坐标 系( p - x y 2 ) 的方程，即有 _ d f + 国f ：0 ，掣+ m + 以f ：01 - ( 1 0 ) d sd s 1 2 2 弹性杆的k i r c h h o f f 方程 设内力主矢，和内力主矩m 在截面主轴坐标系中的投影式分别为 f = 巧矾+ e d 2 + e d 3 ，m = m i d l + 如畋+ 毛d 3b ( t 1 ) 利用1 - ( 6 ) 和1 - ( 1 1 ) ，得到方程l - ( 1 0 ) 组在惯性坐标系( p 一舻) 的投影式为 i d e , + 鲍e 一鸭e = o 鲁+ 皑只一弼e = o 警+ 锡e 一鸭互= o 盟+哆坞d j + 缎朋， s 咝+皑m，d l s 31 咝+弼鸩ds 鸭鸩一最= 0 坞+ 互= 0 a ) 2 m j = 0 1 - ( 1 2 ) l - ( 1 3 ) 由于忽略了杆的体积力，各截面的内力主矢，为常矢量。将，的方向定为惯性 坐标轴f ，f 轴相对截面主轴坐标系各轴的方向余弦记为口，7 ，则有 互= f o r ，e = f p ，e = f 厂l ( 1 4 ) 其中f - - i f i 。将上式代入方程组l 一( 1 2 ) ，则有 7 第一章弹性杆的k i r c h h o f f 方程及其四元数的引入 _ d t t + 鸱y 一坞声：0_ + 鸱y 一坞= ( i s 警+ 皑口一鳓y ：o 孟+ 皑萨鳓y _ 0 誓+ q 一q 口= o 当弹性杆具有原始曲率和扭率时，设砰和趟为弹性杆的原始曲率， 杆的原始扭率，则截而内力主矩可表示为 m = 彳( q 一衅) ，鸠= 曰( 哆一缱) ，鸭= c ( 皑一霹) 将上式代入方程l 一( 13 ) ，可得到弹性杆的平衡方程 彳亟+ ( c b ) ds 、7 b do)2+(彳一c)ds 、7 c 堕+ ( b 一彳) ds 、7 鲍一趟 皑一霹 q 一钟 1 - ( 1 5 ) 衅为弹性 l - ( 1 6 ) 其中a ，b 分别为截面关于x 轴和y 轴的抗弯刚度，c 为截面关于z 轴的抗扭刚度。 常微分方程组1 - ( 1 5 ) 和1 - ( 1 7 ) 为描述弹性杆平衡的k i r c h h o f f 方程。若给出适当的初 值条件，则求解k i r c h h o f f 方程可得到口( j ) ，( s ) ，y ( s ) 及q ( s ) ，哆( j ) ，坞( s ) 从而确 定刚性截面在惯性坐标系巾的姿态。 1 3 欧拉参数的引入1 3 0 1 1 3 l i 1 3 2 i 1 3 1 方向余弦矩阵及其欧拉参数表示 将惯性坐标系( d 一勃f ) 平移到曲线r ( s ) 上任意一点p ，记平移后的坐标系为 ( p 一勿f ) ；以p 为原点建立弹性杆的截面主轴坐标系( 尸一彤) ；根据刚体的有限转 动定理“1 ，将坐标系( p 一勿f ) 绕过p 点的某个瞬时轴h 一次旋转矿角后就可以与截 面坐标系( 户一班) 完全重合。记瞬时轴h 的基矢量日相对于坐标系( 尸一舻) 的方向 余弦为扛，红，吃。设有一个过点p 的向量口，记口相对于坐标系( 尸一勿f ) 和截面主轴 坐标系( p - x y z ) 的坐标列向量分别为a 。和a 。那么存在坐标系( 尸一勿f ) 相对于坐 标系( 尸一彬) 的方向余弦矩阵a 。，使得 r 卜 1 0 o = = 盯 f f 0 一 + = 、l_，、l-，、l_， 。一p，p 衅 衅 谜 一 一 一 鸭 q 缝 青岛大学硕士学位论文 口o = a o l a l f 砰( 1 - c o s 矽) + c o s 痧红扛( 1 - c o s 矽) - h 3s i n e 4 i = l 红矗( 1 一c o s 矿) + 呜s i n e嘭( 1 一c o s # ) + c o s # 【呜啊( 1 一c o s t , ) 一吃s i n e 鬼吃( 1 一c o s 妒) + s i n 将a o 。的各元素用半角公式化作以s i n ( g ) 2 ) 署q c o s ( 矿2 ) 示，定义以下四个符号 g ，= c 。s 善，吼+ 。= 吃s i n 罢( 七= l ，2 ，3 ) l 一( 1 8 ) 称叮= ( g t ，q 2 ，q 3 ，q 4 ) 7 为e u l e r 参数。显然彳+ z + 钙2 十吼2 一l = o 。将a 。的各元素用e u l e r 参数表示，并记用e u l e r 参数表示后的4 0 。为q ( 5 ) ，则 q ( s ) = 2 ( q 卜g ；) 一l 2 ( q 2 q 3 + g i q 4 ) 2 ( q ：q 4 - q i 3 ) 2 ( q 2 吼- q l q 4 ) 2 ( 彳+ 菇) 一l 2 ( q 3 q 4 + g i q 2 ) 1 3 2 弹性杆中心线的欧拉参数表示 2 ( q 2 q 4 + g l q 3 ) 2 ( q 3 q 4 一吼g ：) 2 ( 彳+ 刃) 一l 1 一( 1 9 ) 设弹性杆中心线的端点r ( o ) 与惯性坐标系的原点d 重合，即r ( o ) = 0 。由卜( 5 ) 式得 ，( s ) = r 以( 仃矽仃 1 一( 2 0 ) 将d 3 ( s ) 相对于( d 一勿f ) 的方向余弦代入上式，得到q 。( k = l ，2 ，3 ，4 ) 表示的尸点的 惯性坐标善( s ) ，刁( s ) ，f ( s ) ，以确定弹性杆的中心线，即 孝( s ) = 2 f ( g ：( 仃) 吼( 盯) + g l ( 盯) 9 3 ( 盯) ) d 万 刁( s ) = 2 r ( 吼( 仃) 吼( 仃) 一q ，( o r ) q ：( 仃) ) d 万 f ( s ) = r ( 2 ( g ? ( 盯) + 幺( 仃) ) 一1 ) d 将上式两边求关于s 的导数，得到 善小) = 2 ( 9 2 ( s ) q 。( j ) + 吼( s ) q 3 ( s ) ) r f ( s ) = 2 ( g ，( s ) g 。( s ) - q 。( s ) q ：( j ) ) 卜( 2 1 ) f 7 ( s ) = 2 ( g ? ( j ) + 口：( s ) ) 一l 9 ll_ 痧 n n 西r 吼 吼 5 窖红啊 一 + 力力力 瞄 螂 | 宝 一 一 一 l l l ，【，i 啊红砖 吃吗 第一章弹性杆的k i r c h h o f f 方程及其四元数的引入 1 3 3 弯扭度矢量及内力主矢的欧拉参数表示 把q ( s ) 分解成两个3 4 矩阵r 与r 的乘积，即q = r + r 7 f ，- 9 2 q l 吼一吼、1 - q 2q l q 4 q 3 、 足2l二耄-吼q4一q吼l吼j足=、72一q吼4吼qlqlq 4：2 jl 一吼 9 3一9 2 、一吼 9 2 g l 刚体绕0 点转动的角度极小，可视为无限小量时，称为刚体的转动为无限小转 动。无限小转动a o 可用欧拉参数的无限小增量a q ，( _ ，= l ，2 ，3 ，4 ) 表示。设在 ( d 一勃f ) 中，a 0 = a 0 ( 。) = ( 硝叭，绣0 1 ，谚。) 7 ，刚体相对( d 一勿f ) 的四元数列向 量为“叭，刚体作无限转动后的四元数列向量为彳0 1 ，则“。与彳( 之问有如下关系： 彳= 础“o 卜( 2 2 ) 其中 础= q l 9 2 9 2g i 9 3吼 q 4 一q 3 一口3一吼 一q 4q 3 q l- q 2 吼g i 将式卜( 2 2 ) 两边减去“叭，得 鲥( 0 ) ：( 础一e ) 以。) - 去五叮a o ( o ) 卜( 2 3 ) 其巾 h a o = a o - 4 0 = ( g i ，匈2 ，却3 ，匀4 ) 7 将卜( 2 3 ) 式两边左乘盈，得到 a 0 ( o ) = 2 r + i ( o 卜( 2 4 ) 将q ( s ) 左乘卜( 2 4 ) ，使无限小转动a o 的投影列向量变换剑刚体坐标系，记作护， 则 a o 1 = q ( s ) a o o = 2 r r r r h a o 卜( 2 5 ) l o 青岛大学硕士学位论文 故无限小转动秒沿( 尸一垆) 各轴的投影，略去上标，得到 a o , = 2 ( - q 2 g l + g l 9 2 + 9 4 仍一吼吼) a 0 2 = 2 ( - q ，g 。一吼幻2 + g 9 3 + q 2 9 4 ) 卜( 2 6 ) 岛= 2 ( - q 4 a q i + q 3 a q 2 一9 2 吼+ g l 吼) 将卜( 2 6 ) 各项与厶相除，并令厶- - - - ) 0 ，得到用欧拉参数变化率表示的弯扭度各分 量 饵= 2 ( - q ：g ：+ g l q ：- i - q 4 q 3 一q 3 q 4 ) c 0 2 = 2 ( 一9 3 9 ：一9 4 毹+ g i q 3 7 + q 2 q 4 7 ) 1 一( 2 7 ) 坞= 2 ( - q 4 q 1 7 + q 3 q 2 一q 2 q 3 + g i q 4 7 ) l 互= 2 f ( q ：q 4 - q i q 3 ) 五= 2 f ( q 3 q 4 一g l g ：) 卜( 2 8 ) i e = f ( g ? 一云一孝+ 云) 第二章弹性杆的拟h a m i l t o n 方程及其辛算法 第二章弹性杆的拟h a m i l t o n 方程及其辛算法 k i r e h h o f f 弹性杆及其相应的动力学比拟方法近年来在d n a 模型的结构和运动 的研究中得到应用。而d n a 的超细长、大变形等特征给k i r c h h o f f 弹性杆的建模和 计算提出新的问题。本章基于超细长弹性杆的计算考虑，利用e u l e r 参数代替传统 的e u l e r 角的描述方法，建立了描述弹性杆平衡的拟h a m i l t o n 方程；采用辛算法建 立了弹性杆长距离保结构计算的数值计算模式并给出了数值结果。 2 1 弹性杆的能量密度函数 2 1 1 刚性截面上任意点的应变分量和应力分量 考虑一长度为三的具有刚性截面的弹性杆，以截面的几何中心p 为原点建立截 面主轴坐标系( p - x y z ) 。设刚性截面上任意一点q ，则q 点的坐标p 在( p - x y z ) 中 可表示为 p = x d l + 22 - ( 1 ) 忽略弹性杆中心线的拉仲变形，截面上每一点的位移只随截面的转动发生变化。根 据刚体转动的原理【3 3 】，截面的无限小角转动a o 引起的q 点的无限小位移脚可表示 为 a u = a o x p2 - ( 2 ) 将上式两边都除以z ，并令止一0 ，则可得 娑：缈p2 一( 3 )_ = 缈。( 3 ， 把矢量h 在截而主轴坐标系巾的分量记为( f = l ，2 ，3 ) ，则可把2 一( 3 ) 式写成分量形 式 誓= 一皑y ，挚o z = 皑x ，誓= q y 一皑x 2 - ( 4 ) 0 = 一皑y ，- = 一= 皑x ，- = = q y 一皑x 2 。( 4 ) d zo z 对上式积分，并假设u ( o ) = 0 ，则可得 = - a 五y z ，u 2 = t 0 3 x z ，u 3 = ( 0 4 y a 是x ) z 2 - ( 5 ) 把坐标z 限制为小量，可利用c a u c h y 公式f 3 4 】，得到q 点的应变分量如下 毛= 粤= o ，乞= 誓= 0 ，岛= 誓= c q y 一咝x 2 - ( 6 ) o xd ，o z 1 2 青岛大学硕士学位论文 乃：= 携。= 粤+ 粤= o ，托，= 乃。= 粤+ 誓= 一缈，如= = 挚+ 誓：螂2 奶 砂 o x d z蹴 。 一d z咖 。 假设弹性杆为均匀且各向同性的，则q 点的应力分量可表示为 q = 吒= 0 ，吒= e ( q y q x )2 - ( 8 ) 互2 = 乞l = 0 ，互3 = 毛l = 一g a r y ，= 弓2 = g 坞x2 一( 9 ) 其中，e ，g 分别为弹性杆的杨氏模量和剪切模量。 2 1 2 弹性杆的能量密度函数 弹性杆的应变能是e h 多 l - 力功转化成的弹性变形势能，故可用应变分量和应力分 量表示为【3 5 】 巨= 了1r ( q + 吒乞+ 吒岛+ 气：兀：+ 如+ 互，托，) a s d s 2 ( 1 0 ) 刚性截面的作用力对尸点的内力主矩可t h 应力分量的积分式表示为 m = 吒声，鸩= 一吖d s ，鸩= ( 乞，x 一互，y ) d s 2 - ( 1 1 ) 其中s 为截面的积分域。把2 一( 8 ) 、2 - ( 9 ) 代入上式，并注意到在截面主轴坐标系下有 x d s = 0 ，上y d s = 0 ，点x y d s = 0 则截面作用力相对点尸的主矩可表示为 m i = 吒闪s = 么q ，鸩= 一j o 3 x d s = b o o z ，托= ( 砀x 一互，y ) d s = c a 毛 其中a ，b 分别为截面绕x 轴和y 轴的抗弯刚度，c 为截而绕z 轴的抗扭刚度，且 a = e y 2 d s ，b = e 上x 2 d s ，c = g 点( x 2 + 少2 p 故弹性杆的应变能可利用弯扭度表示为 巨= 寻r 么q 2 + b 够2 + c 皑2 凼 设外力主矢为f ，y 为孝轴相对z 轴的方向余弦，则外力势能 e p = 一毫f y a s 其中，f = i t i 。故弹性杆的总势能 1 3 第二章弹性杆的拟h a m i l t o n 方程及其辛算法 e = 巨+ 乓= 啦( 讹) 2 + ) 2 + c ( 酊) 一，十 2 - ( 1 2 ) 记 r = 丢( 彳q 2 + 刀皑2 + c 鸭z ) 一f 7 2 ( 1 3 ) 为总势能密度函数。 2 2 弹性杆的拟h a m i l t o n 方程 由变分学的e u l e r 公式，将q ，q ，皑和f y 的e u l e r 参数表示1 - ( 2 7 ) 、1 - ( 2 8 ) 代入 能量变分2 - ( 1 2 ) ，注意到的隐含约束 g ( g ) = g 。g 一1 = 02 - ( 1 4 ) 则在约束2 - ( 1 4 ) 下，能量变分2 - ( 1 2 ) 取极小的l a g , r a n g e 方程为 忙( 嚣) 一薏- - 2 ：， 【g ( 窖) = 0 其中三= i 1g 盯m ( q ) q 一f 7 一g ( q ) 2 ，五为拉格朗日乘子， 爿q ；+ b 霹+ c 留 一爿吼碍2 + ( b c ) q 3 q 4 一励l 7 3 + ( c a ) q 2 q 4 ( a b ) q 2 q 3 一c q i 吼 一4 q i q 2 + ( b c ) q 3 q 4 月g j + b 霹+ c 谚 ( a b ) q i q 4 一c q 2 q 3 ( c a ) q i 吼一b q 2 q 4 一砌l 吼+ ( c a ) q z q 4 ( a b ) q i 吼一c 臼2 吼 月云+ b q ；+ ( 勾； 一a q 3 q 4 - i - ( b c ) q , q 2 把方程2 ( 1 5 ) 展开，有 竺。肚八伽勺 其中 f ( q ，q ) = 4 a q ；a l + 4 b q ；媲+ 4 c “呜一2 圪g l - - - 4 a q ；o c + 4 b q 4 0 ) 2 + 4 c q ；鸭+ 2 e 9 2 - - - 4 a q 4 。c z 一4 却叙+ 4 c q 2 呜+ 2 f o q 3 4 a q ；o l - 4 b q 2 7 0 4 4 回恕一2 磊9 4 ( 彳一b ) q 2 q 3 一c q l q , ( c a ) q i q 3 一却2 吼 一a q 3 q 4 + ( b c ) q i q 2 a q ；+ b 旌+ c q ； 2 - ( 1 6 1 通过计算可知，矩阵l m ( g ) i = o ，但我们有 定理1 对于任意初始值，方程2 - ( 16 ) 存在唯一解。 证明：通过；j i a 变换1 - ( 1 8 ) ，弹性杆k i r c h h o f f 方程的解转化为方程2 - ( 1 6 ) 的解。从 而我们可知方程2 - ( 1 6 ) 的解的存在性。另外，经过计算可知，矩阵 1 4 青岛大学硕士学位论文 掣 2 羽乃 的5 个4 阶行列子式分别为0 和a b c q , ( f = l ，2 ，3 ，4 ) 。由2 - ( 1 4 ) 式可知，q 不可能是 零向量，即a b c q ，( i = l ，2 ，3 ，4 ) 至少有一项不为零。从而我们可知，对于任何四元数 p ：a - 觇7 ：m ( g ) g ， 2 ( 1 8 ) o q 。 h ( p ，g ) = l p r m + p + u ( 窖) fp - - 露p ( p ，g ) p 7 = 一z l ( p ，g ) 一岛( g ) 7 l 2 - ( 1 9 ) 1 0 = g ( g ) f ( q ，口，) = 一芸m ( g ) 矿+ v 。t ( q ，n v 。u ( g ) c l f l 。 = 一鸭n 毒( y 批汀) - v q u 2 啦。， j 州们肛，2 - ( 2 1 ) 【q r x = 0 u 一7 加一品( 护眦帅) 一m q q 2 - v q u 1 5 第二章弹性杆的拟h a m i l t o n 方程及其辛算法 = 一见( 办g ) 一m q q 召 2 - ( 2 2 ) 即定理得证。 口 因为方程组2 - ( 1 6 ) 和2 - ( 1 9 ) 没有实际的动力学意义，我们把2 - ( 1 6 ) 和2 - ( 1 9 ) 分别 称为弹性杆的拟l a g r a n g e 方程和拟h a m i l t o n 方程。 2 3 弹性杆拟l a g r a n g e 方程的应用 我们从一个新的角度，即采用邻界面的思想，推导出弹性杆曲面方程。借助于 k i r c h h o f f l 拘动力学比拟方法，我们把超细长弹性杆的表面看作一条刚性封闭的平面 曲线点集x 。沿着空间曲线r ( s ) 移动和旋转得到( 图2 1 ) 。 图2 1 弹性杆不意图 设x 。为初始截面点集，x ( j ) 可以理解成x 。在时间s 处的位置，则我们有 x ( s + a s ) = a ( s + a s ，s ) ( x ( j ) 一，( s ) ) + r ( s + a s ) 2 - ( 2 3 ) 其中a ( s + a s ，s ) 是x ( s ) 到x ( s + a s ) 的旋转映射。记q ( s ) = ( 以( j ) ，d ：( s ) ，d ，( j ) ) 是坐标 系( p - x y z ) 到( d 一勿f ) 的方向余弦矩阵；那么截面x ( s ) 到x ( s + a s ) 的旋转变换可 通过两次坐标变换米实现，即有 a ( s + a s ，s ) = q ( s + 厶) q ( s ) 。 那么显然有a ( s ，s ) = q ( s ) q ( s ) t = l ，则 1 6 青岛大学硕士学位论文 三砌a ( s + a s , s ) - a ( s , s ) ：三f 聊垡g 竺! 垡! 兰! ：二望! ! ! 垒! 兰z 厶| _ o a sh , - - - o ， ：壬唿业些巡r ：辈出) 2 ( 2 4 ) ，_ + o a s 、7 d s 。、7 、7 在方程2 一( 2 3 ) 两边都减去x ( j ) ，可得 x ( s + a s ) 一? r ( j ) = 4 ( s + a s ，5 ) ( x ( s )