（应用数学专业论文）dna力学结构的建模数值仿真和图形后处理.pdf

摘要 以超细长弹性杆作为d n a 的双螺旋结构模型来研究d n a 的力学结构和动力学 性态，基于k i r c h h o f f 的动力学比拟方法，导出了没有外界约束的条件下描述静 态弹性杆表面的曲面微分积分方程组，在此模型的理论基础上，将弹性杆结构 模型视为刚体有限转动模型，推导出适合数值计算的四元数形式的微分方程。采 用相关处理，通过对弯扭弹性杆的数值实验，讨论d n a 分子的空间结构特性和数 值仿真及其图形后处理。文章用曲面回扫方法结合刚体有限转动的有关理论研究 d n a 弹性杆的空间结构及空阳：j 曲面图形的描绘，包括曲面的空间形状和弹性杆自 身的扭转，为分析弹性杆的自身接触问题和d n a 动力学性态的计算机仿真及其图 形后处理提供了模型和算法的支持，是研究d n a 的动力学模型的一个方便的工 具，通过实例发现这一方法也可以用于其它一些问题的图形处理。文章引入描绘 d n a 弹性杆的h a m i l t o n 函数，研究了包括d n a 弹性杆力学模型的h a m i i t o n 方程 的表示，相应的辛算法的推导以及计算结果的分析等问题；探讨了利用保结构的 数值方法计算超细长弹性杆的力学结构问题，针对d n a 弹性杆超细长特性，通过 数值实验发现辛算法在超长弹性细杆的长距离数值模拟中可以保持其结构特征， 从而给出了辛算法在结构问题中的应用。 关键词：超细长弹性杆；四元数；h a m i i t o n 函数；保结构算法；数 值仿真 a b s t r a c t i nt h i sp a p e r , as u p e r - l o n ge l a s t i cs l e n d e rr o di su s e da sas t r u c t u r a lm o d e lo f d n at o s t u d yt h ed y n a m i c a lp r o p e r t i e sa n ds t r u c t u r eo fd n a m o l e c u l e s b a s e do n k i r c h h o f f sa 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 l i n t e g r a le q u a t i o ns y s t e mi sg i v e nt o d e s c r i b et h es u r f a c eo fas t a t i cs u p e r - l o n ge l a s t i cs l e n d e rr o d b a s e do nt h i sm o d e l ， a n dr e g a r d i n gt h ee l a s t i cr o da saf i n i t er o t a t i o nm o d e lo fr i g i d b o d y ，e q u a t i o n s s u i t a b l ef o rn u m e r i c a lc o m p u t a t i o ni nq u a t e r n i o ns t y l ea r ed e r i v e d t h r o u g ht h e n u m e r i c a le x p e r i m e n t so fe l a s t i cr o dw i t hc u r v a t u r e ，t h ep r o p e r t yo ft h er o d ss p a t i a l s t r u c t u r ei sd i s c u s s e d n u m e r i c a ls i m u l a t i o na n df i g u r ep o s t - p r o c e s s i n ga l g o r i t h m sa r e d e s i g n e d c o m b i n i n gs w e p tv o l u m e st e c h n i q u e sw i t ht h et h e o r i e so fr i g i db o d y , a n e wm e t h o di si n t r o d u c e di ns t u d y i n gt h es p a t i a ls t r u c t u r eo f e l a s t i cr o d s ，i n c l u d i n gi t s s p a t i a ls h a p e s ，t h et o r s i o n sa n dt h es u r f a c ec o n t o u r s ，w h i c hp r o v i d e sac o n v e n i e n tt o o l f o rn u m e r i c a ls i m u l a t i o na n dp o s t - p r o c e s s i o na sw e l la sf o rs t u d y i n gt h ed y n a m i c a l p r o p e r t i e so f ad n a b yi n t r o d u c i n gt h eh a m i l t o nf u n c t i o no f t h ee l a s t i cr o dm o d e l ，a s t r u c t u r e - p r e s e r v i n g n u m e r i c a lm e t h o di si n t r o d u c e di n s i m u l a t i n g m e c h a n i c a l s t r u c t u r e so fs u p e r - l o n gt h i ne l a s t i cr o d s ，t h ep r o b l e m ss t u d i e d i n c l u d i n g t h e r e p r e s e n t i n go ft h em e c h a n i c sm o d e lo fe l a s t i cr o di nh a m i l t o ne q u a t i o nf o r m ，t h e d e r i v a t i o no fc o r r e s p o n d i n gs y m p l e c t i ca l g o r i t h ma n dt h ea n a l y s i so fc o m p u t a t i o n r e s u l t s a i m i n ga tt h es u p e r - l o n gp r o p e r t yo ft h ee l a s t i cr o d s ，n u m e r i c a le x p e r i m e n t s a r em a d e t h es y m p l e c t i ca l g o r i t h mi sf o u n dp r e s e r v i n gt h es t r u c t u r ei nl o n g - r u n n u m e r i c a ls i m u l a t i o n i ti st h ef i r s ta p p l i c a t i o no fs y m p l e c t i cm e t h o di ns t u c t u r e p r o b l e m s k e y w o r d s ：s u p e r - l o n ge l a s t i cs l e n d e rr o d ；q u a t e r n i o n ；h a m i l t o nf o r m ；s t r u c t u r e p r e s e r v i n g ；n u m e r i c a ls i m u l a t i o n 青岛大学硕士学位论文 学位论文独创性声明 本人声明，所呈交的学位论文系本人在导师指导下独立完成的研究成果。文中依 法引用他人的成果，均已做出明确标注或得到许可。论文内容未包含法律意义上已 属于他人的任何形式的研究成果，也不包含本人已用于其他学位申请的论文或成果。 论 校享有以任何方式发表、复制、公开阅览、借阋以及申请专利等权利。本人离校后 发表或使用学位论文或与该论文直接相关的学术论文或成果时，署名单位仍然为青 岛大学。 本学位论文属于： 保密口，在年解密后适用于本声明。 不保密影 论文作者签名 导师签名：赵冀戮胥 日期：m 滓6 月f 日 日期：刀口7 年f 自月j 臼 ( 本声明的版权归青岛大学所有，未经许可，任何单位及任何个人不得擅自使 用) 引言 弹性杆问题的研究历史 引言 早在1 7 3 0 年。伯努利( d b e r n o u l l i ) 和欧拉( l e u l e r ) 就已经开始研究弹性杆在 外力和力矩作用下的变形等力学问题。k i r c h h o f f 于1 8 5 9 年在刚性截面假定和无体 积力、无轴向变形的条件下，根据弹性杆的平衡微分方程和刚体定点转动微分方程 之问的相似性，提出了弹性细杆平衡的动力学比拟理论。弹性细杆模型有着广泛的 实际背景，电缆洲、绳索、钻杆、纤维乃至自然界中攀缘植物的细茎都可以以弹性 细杆作为其力学模型。k i r c h h o f f 模型为研究弹性细杆的静力学和动力学提供了新 的途径。弹性细杆的非线性力学因在分子生物学中成功地模拟蹦a 等一类大分子的 静态平衡而再度受到关注和重视。 在分子生物学领域中，自2 0 世纪中期w a t s o n 和c r i c k 提出了d n a 分子的双螺 旋三维结构模型以来，关于d n a 的基础理论研究不断突破，用具有原始扭率的弹性 细杆作为d n a 的宏观力学模型的理论研究得到了实验的肯定，d n a 弹性杆模型中的 重要物理参数如杨氏模量、泊松比、抗弯和抗扭刚度等已经借助生物学实验的手段 和统计物理的方法得到嘲f 】“”，促进了以遗传工程为首的生物技术的迅速发展。因 此，宏观的建模方法和经典的分析工具可以在d n a 的研究中得到充分的应用。 d n 双螺旋结构由两条螺旋形核苷酸链和联系二螺旋线的碱基组成，为了研究这 种三维结构的几何形态及其结构性质，采用具有原始扭率的弹性细杆作为d n a 的宏 观力学结构模型，许多经典的数值计算的基本原理和方法在d n a 的力学结构模型的 研究中得到充分的应用。于是分子生物学与力学、数学等产生了结合点。自7 0 年 代以来，关于蜊a 力学模型的研究逐渐形成一支交叉学科分支，讨论d n a 、r n 以及 细菌纤维等各种生物分子链的几何形态和结构性质。人体细胞的最大染色体所含d n a 分子的螺旋半径约为2 n m ，而长度可达7 c m 。秆长约为半径的3 5 1 0 7 倍。因d n 弹性杆模型的极端细长性和超大变形性，使得对d n a 弹性杆数值仿真的研究完全不 同于对传统弹性杆的研究，自7 0 年代以来，利用数学方法研究和讨论d n a 、h n a 以 及细菌纤维等各种生物大分子链的几何形态和稳定性等阀题成为一个具有重大实际 意义的课题。 研究弹性杆模型的目的和意义 研究此课题的目的是研究超细长弹性杆的数学建模、结构动力学分析如稳定性 分析等，以及数值仿真和图形处理问题。超细长弹性杆模型是以d n a 分子结构和 青岛大学硕士学位论文 运动为实际背景所建立在其他领域如海底电缆，工程钻杆等的研究中也有重要的 应用具体的意义如下； ld n a 分子的弹性力学性质对其生物学功能( 折叠成染色体超螺旋、作为基因转录、 复制的模板) 有极大的影响。建立数学模型来数值仿真d n a 分子的力学结构，得 到的结果对研究d n a 分子的力学性质，动态稳定性等过程有重要的实际意义。对 研究d n a 分子的自制、转录、重组、修复等生物学过程的研究具有重要的参考价 值。 2 对d n a 超螺旋结构进行数值仿真，通过对仿真结果的分析和研究，人类可以更加 科学的根据一定的目的，对d n a 分子进行体外加工操作。再引入受体生物，以改 变后者的某些遗传性状，从而培养生物新类型以此来避免或治疗遗传疾病。 3 超螺旋结构的d n a 分子的三维形状对d n a 分子所处的生理盐溶液中离子的浓度非 常敏感，所以d n a 的初曲率与离子浓度的影响有关；因此研究具有初曲率的d n a 分子构象具有较为普遍的意义，d n a 分子的不同构象，反映了生理盐浓度的大小变 化 4 具有初曲率的闭环d n a 结构在复制、转录、重组、修复中不断改变着自身的拓朴 结构，引起超螺旋参数如扭转数，缠绕数的变化。因此可以通过改变扭转数来获 取具有初曲率的d n a 分子稳定状态的不同构象，以满足不同生物学过程的需要。 国内外研究动态 d n a 分子结构模型的研究是近3 0 年来数学、生物学和力学研究的新课题，一 般情况下，d n a 分子被模拟成均质的、各向同性的、线性圆截面弹性细杆，关于d n a 弹性杆的研究已经有大量的文献发表在生物数学、物理化学、生物化学和高分子化 学以及数学物理等学科的刊物上。但与此紧密相关的弹性杆的研究却已经有近3 0 0 年的历史早在1 7 3 0 ，d b e r n o u l l i 和l e u l e r 就开始研究弹性杆在外力和力矩作用 下的变形问题。1 8 5 9 年k i r c h h o f f 提出了弹性杆平衡的动力学比拟理论。三年后 c l e b s h 在其弹性力学著作中论述了这一理论。k i r c h h o f f 嘲理论作为弹性杆静力学模 型的理论基础，l o v e 在其弹性杆著作中对此理论作了详尽的论述。 2 0 0 0 年以来，国内学者在这一研究领域傲了许多工作。国内最新讨论了受曲面 约束弹性细杆的平衡问题，从动力学的观点讨论了超细长弹性杆的平衡稳定性。 刘延柱等油船1 对d n a 弹性杆模型的力学性质和稳定性作了深刻系统的研究，他 们在弹性杆受到几何约束情况下建立弹性杆的平衡方程。、研究了弹性杆的动力学问 题，推导出了曲率和挠率的动力学方程及其哈密尔顿系统的静力学方程模型。 用数值方法对d n a 弹性杆的力学结构模型进行仿真，对研究d n a 分子生物学 性质极为重要但由于d n a 结构的超细长特性，d n a 结构模型的数值仿真也给我 2 引言 们提出了新的问题，即在数值模拟中存在长距离数值计算的精度和保结构问题。利 用传统的数值方法如r u n g e k u t t a 方法计算，当计算距离很长时，往往出现一团乱 麻的现象，失去了d n a 特有的缠绕结构。如何在计算中保持d n a 的特有结构的研 究目前刚刚开始，但有关保结构算法的研究已经取得了丰硕成果。 冯康于1 9 8 4 年提出了h a m i l t o n 系统的保结构辛算法，从理论上清晰的阐明了 传统的数值方法导致系统“能量”损耗的根本原因。并用大量数值实验结果证明对 于保结构动力系统的计算，尤其是对于长距离的数值模拟，辛算法具有传统数值方 法无法替代的优势保持系统的结构特征。冯康首次将h a m i l t o n 系统几何结构的 思想引入数值分析，开辟了一个新的研究领域，并在近十年的时间里在辛算法及其 在数值分析中的应用等方面进行了广泛深入的研究，取得了丰硕的成果。使人们越 来越意识到辛算法的重要性。并用大量数值实验结果证明辛算法已在很多领域包括 天体力学，量子化学，非线性波，不可压流体，大气物理和地物勘探数据处理等， 找到了很好的应用 主要研究内容 d n a 弹性杆的数学建模问题 由于d n a 超细长弹性杆的长久的历史背景和广泛的实际背景，目前已经有不同 类型的模型嘲，如由k i r c h h o f f 方程所导出的欧拉角形式的模型，由以曲率和挠率 为变量，建立的由s c h r o d i n g e r 方程嘲所描述的数学模型等，但由于上述模型在不 同程度上都存在缺陷，如欧拉角出现奇点问题，等本文将根据研究的需要，建立由 欧拉参数或四元数表示的d n a 弹性杆模型以及h a m i l t o n 形式的能量模型。 d 帆弹性杆数值仿真问题 k i r c h h o f f 的动力学比拟方法可以类似地应用于弹性杼数学模型的图形处理 本文将k i r c h h o f f 的动力学比拟方法和四元数方法相结合，得到弹性杆结构图形的 微分代数方程模型和相应的图形描述算法。这一方法不但可以解决弹性杆的结构仿 真问题，而且可以推广为计算机仿真的图形回扫算法n 町保结构算法在需要保持固 有结构的长距离数值模拟的模型中有着传统数值方法无法比拟的优势，如何利用这 两种方法对d n a 弹性细杆进行数值模拟也是文章解决的主要问题。 本文的内容安排 本文的引言为全文的绪论，主要介绍了弹性杆问题的研究历史，研究d n a 弹性 杆模型的目的和意义，d n a 弹性杆模型国内外研究动态，本文主要研究内容和解决 3 青岛大学硕士学位论文 的主要问题等。第一章主要建立了d n a 弹性杆力学结构的几种数学模型、介绍了各 种模型的相关特点。第二章对第一章的模型进行数值仿真，设计相应的数值算法， 并对仿真的数值结果进行分析。第三章主要讨论了弹性杆方程的h a m i l t o n 形式及其 相对应的辛算法，讨论了辛算法在结构问题计算方面中的优势。 4 第一章弹性杆的k i r c h h o f f 方程及其相关模型 第一章弹性杆的k ir o h h o f f 方程及其相关模型 一般情况下，d n a 分子被模拟成均质的、各向同性的、线性圆截面弹性细杆， 关于d n a 弹性杆的研究已经有大量的文献发表在物理和化学物理、生物化学和高 分子化学以及分子生物学等学科的刊物上，其中大多数工作是以鼬r c h h 0 丘理论为基 础的解析研究。 1 1 模型的提出与改进 1 1 1 基本假定 讨论一长度为l 的细杆，杆截面的几何中心连成的不自交的空间曲线称为曲杆 的中心线，假定： 1 ) 中心线在变形前后均为2 阶以上光滑曲线； 2 ) 杆的长度和曲率半径远大于横截面的尺寸； 3 ) 横截面为刚性截面； 4 ) 忽略弯曲引起的剪切变形，横截面与中心线正交： 5 ) 忽略中心线的拉伸变形，任意两截面沿中心线的距离不变： 6 ) 杼的弹性常数为常值，应力和应变满足线性本构关系； 7 ) 相邻截面可绕中心线作相对扭转，扭角为的弧长的连续函数。 1 1 2 弯扭度矢量的引入 极端细长且具有超大变形的弹性细杆不同于传统弹性力学的研究对象。在刚性 截面假定的基础上，曲杆的几何形态由截面沿中心线的移动和扭转所体现。中心线 的几何形态的改变由杆的弯曲变形引起，截面中心线的转动体现杆的扭转变形截 面的弯扭度矢量，即截面的无限小位移矢量对弧坐标的变化率是描述曲杆局部几何 形态的参量，由中心线的曲率和截面的扭率构成。 以杆中心线上任意点尸为原点，建立截面的f r e n e t 坐标系p t n b 设z 为x 轴与轴，y 轴与口轴的夹角，即截面相对f r e n e t 坐标系的扭转的角度， 如图( 1 1 ) ，在曲线几何学理论。1 上中，可以定义d a r b o u x 矢量“1 ，( j ) = k ( s ) b + r ( s ) r l - ( 1 ) 青岛大学硕士学位论文 圈1 1 弹性杆截面扭转示意圈 其中| i ( s ) ，r ( s ) 为曲线的曲率和挠率。 在此基础上，定义新的矢量 o j - - - - c o f + 馔她= ( o l e l + o ) 2 e 2 + 0 7 3 e 3 1 - ( 2 ) 其中 0 3 1 = k s i n x , 0 7 2 = k c o s x ，伤= f + 譬 则称1 - ( 2 ) 所定义的向量为弹性杆的弯扭度矢量。 1 1 3 弹性杆的平衡方程 在三维空间的惯性坐标系d 一勿f 中讨论一长度为l 的细长弹性杆，杆的始端和 终端份分别记为p o 和尼p 和p ，点相对固定参考点d 的矢径分别为，和，+ a r ，相 对p 和p 点的弧坐标分别为j 和s + a s 以杆中心线上任意点p 为原点，建立截面 的主轴坐标系p 一习瞄各坐标轴的基矢量分别为岛，t 2 ，岛，其中z 轴与曲线c 的 切线轴r 重合，即c ，= r 。考虑微元弧段杆，的平衡( 如图l 。2 所示) ， 图1 2 弹性杆徽元受力分析图 p 和p 点处截面的内力主矢和主矩为别为e 膨和f + a f , m + z l m ，利用微元体的 平衡可以导出弹性杆的平衡微分方程“1 为： 6 第一章弹性杆的k i r c h h o f f 方程及其相关模型 竺：o d a l q 3 ) d m + r f ：0 d s 改为相对截面主轴坐标系进行，基矢量r 改记为岛，方程1 部) 改写为 坚+ ，：o ，l 删 d m + x m + p ，f ：0 r = 圜 西= b 哆屿) l 锕 誓懈e 飞驴o - 峨f f + ( ，巧一。e ) = o 鲁地五吨耻o 7 l - ( 7 ) 青岛大学硕士学位论文 彳訾+ ( c b ) 仍鸭一c 砖吐一即= o 曰i d ( 0 2 + 似一c ) q 伤+ c 硝q + 鼬= o 1 - ( 8 ) c i d o ) 3 + ( 丑一铷：q = o 、 以上封闭的微分方程组1 - ( 7 ) ，1 - ( 8 ) 为描述弹性杆的k i r c h h o f f 方程。 1 2 弹性杆四元数形式的k i r c h h o f f 方程 1 2 1 欧拉四元数 根据刚体有限转动的欧拉定理定义欧拉四元数g ， 窖= q q ：q 3 q 4 ) ，四元数集合记为q 。 q 上的单位球面 墨= 矗l 口? + g ；+ 爵+ 订= 1 上的元素可以表述为 将四元数吲记为 1 - ( 9 ) = c o s 罢州= s i n 毋一：l ，2 ，3 )1-(10)ql q k s m - 三( k 2 = 三，“= = i ，) j 二二 其中 ( 七= l ，2 ，3 ) 为弹性杆截面的有限转动轴h 相对惯性参考系0 一善，嘭的方向余 弦妒为有限转动角。将弹性杆截面主轴坐标系p 一驴相对惯性坐标系0 一翻f 的 方向余弦用四元数表示，得到方向余弦矩阵m 彳+ g ；一爵一g ：2 ( q ：q 3 一9 1 吼)2 ( 口：q 4 + 铂9 3 ) 1 q ( j ) = 缸( s ) = l2 ( q 2 q 3 + q l q 4 ) 衍一钌+ 酊一磊2 ( q 3 q 4 一q l q ：) i l - ( 1 1 ) 【2 ( q 2 q 4 一吼9 3 )2 ( q 3 q 4 + q - q 2 ) 疗一旌一爵+ 玉j 将q i ( s ) ( 七= 1 , 2 ，3 ，4 ) 对j 的导数作为新的变t q 。o ) ( i = l ，2 ，3 ，4 ) ，即有 孕= g ( 七；1 ，2 , 3 ，4 )1 - ( 1 2 ) m 弹性杆的弯扭度蛾( f = l ，2 3 ) 可用欧拉参数及其导数表示为“1 ： q = 2 ( 一q 2q l + q l g + 吼g 一吼a ) - - - 2 ( - q 3q l g 伤+ 吼q 3 + 9 2 幺)1 - ( 1 3 ) 吃= 2 ( _ 9 4q l + 吼伤一窖2 幺+ g 。a ) 将1 - ( 1 3 ) 代入到1 ( 7 ) ，l - ( 8 ) 得到： 0 第一章弹性杆的k i r c h h o f f 方程及其相关模型 誓+ 珥( 一吼q q , q 2 + 吼q 3 + 叮2 q 。) e 一( 吨奶+ q 3 q 2 一9 2 q 3 + q l q 4 ) f 2 = o - 幔d + 2 ( - q , q 1 + 吼q 9 2 q 3 + q i q 4 ) f i 一( - q 2 蜴+ 吼q + 吼q 3 一吼q 。) e 】= o l - ( 1 4 ) 鲁+ 2 ( - q 2 q i + 9 1 q 2 + 五q 3 一吼q ) e 一( - 仍q l 一吼q 2 + 吼q 3 + g ：q 。) e 】= o - q - a + 孵+ 识瑙+ 2 ( q 砷( 硼一娼嘲+ 必) l 硼+ 鹏一鹕+ 孵i i 训l 州砷( 硼+ 鸱一鹕+ 心) 卜蝴一幔+ 鹏+ 必_ i 趟2 ) 1 4 r ：4 1 = o - - 蝴- q a + 弼+ 必+ 2 ( 卅功( 刮a + q , g 一必+ 岛g ) 硼+ 鹕+ 鹕一磁一( 州 1 - ( 1 5 ) 卅司( 硼+ 娼+ 鹏一哟) 【吲a + 鸱一必+ 啦- 回2 ) l + ( r , 4 l = 0 胡+ 码磁+ 循+ 2 ( 酮( 蝴+ 哂q + 必一心) i 硼一啦+ 鹕+ 必_ 巡2 ) l - 们( 嘲一鸱+ 鹏+ 码) 【硼+ 鹏+ 娼一必i 引i = o 由于无分力，截面的作用力的主矢，为常矢量，可将f 的方向作为惯性坐标系的f 轴，可确定f , ( s x i = l ，2 3 ) 为欧拉参数的以下函数 石= 2 磊( 9 2 q 4 一q l q ，) e = 2 最( 吼g + q l q 2 )1 7 1 6 ) e = 2 f o ( q j 一云一酊+ 云) 其中e o 为端部作用力的模。将卜( 1 6 ) 带入卜( 1 4 ) 后与1 - ( 1 2 ) ，1 ( 1 5 ) 构成封闭的方程 组。 1 3 考虑分布力的弹性杆的s c h r o d i n g e r 方程 以杆中心线上任意点p 为原点，建立截面的f r e n e t 坐标系p t n b 设z 为x 轴与轴，y 轴与君轴的夹角，即截面相对f r e n e t 坐标系的扭转的角度， 曲线的中心线确定以后，只需再确定扭角随弧坐标的变化规律，曲面的空间几何形 态便可以完全被确定，其中p - x y z 和p t n b 与之间的方向余弦如下表所示： 表1 1 方向余弦 丑r x c o s zs m z o ) ， 一s m zc o s z o z oo1 9 青岛大学硕士学位论文 根据d a r b o u x 矢量 街f ( j ) = 七( s ) 口+ r ( s ) t 和弯扭度矢量 ，毋= + ( 挚，= 纬毛+ 呸f 2 帕3 白 卜( 1 7 ) 的关系，截面的主轴坐标系e - x y z s 的变化规律烈j ) 由以下微分方程确定 棚： i d e t = x p ，( i = 1 ，2 ，3 ) 1 - ( i s ) 将1 - ( 1 8 ) 带入1 - ( 1 7 ) 后，导出。 要：( 升掣泌一心 ds 鲁咄+ ：+ r 咖炽 t 鲁= r ( c o s 妃一s 斑觞) 上面所述两个模型都未考虑弹性杆之间的分布力，实际上若考虑弹性杆有可能与其 他物体接触，则存在分布的接触力。因此更完善的平衡方程应将分布力和分布力偶 考虑在内。因此需在模型的微元段上加上单位长度的分布力和分布力偶埘。则相 应的弹性杆的平衡微分方程为： 坚+ ，；0 | d 。s ， 1 ( 2 0 ) 坐+ r f + 扰：o 将求导过程改为相对截面主轴坐标系进行，基矢量r 改记为岛，方程改写为 将分布力，在主轴坐标系p 一印譬中投影，分量记为z ( j = l ，2 3 ) 。若不计摩擦力和切 向的分力，则方程中，的切向分量为0 。若分布力全部略去，方程组转化为1 - ( 7 ) 。 若将在f r e n e t 标架下表示为 f = 五+ 五bl ( 2 2 ) 利用f r e n e t 标架和主轴标架p x y z 的关系，结合1 ( 2 2 ) ，在k i r c h h o f f 的圆截面弹 i o m 却协 印 吖 扣产 第一章弹性杆的k i r c h h o f f 方程及其相关模型 性杆的简化假设下，导出利用超细长圆截面杆的s c h r o d i n g e r 方程，得到曲率k 和挠 率t 的微分方程组 窘叫r 一争2 邯+ 譬= 鲁 r 韭d a + 2 塑d a ( f 一马2 = 鲁 、7 彳 青岛大学硕士学位论文 第二章弹性细杆的数值仿真 自从d n a 三维螺弦结构模型提出以来，关于d n a 弹性细杆模型的研究工作已 经取得了不少成果，其中大多数是以k i r c h h o f f 理论为基础的研究嘲。另外需要注 意的是在弹性杆模型的研究中所采用的数值方法，不管采用何种方式的研究，结果 都需要将d n a 的三维结构模拟出来，这样才能直观的表达分子的空间拓扑结构， 为了更加清晰的表达d n a 螺弦的状态，需要用d n a 弹性杆的扭转来实现，本章采 用m a t l a b 编制计算程序”1 对具有弯曲和扭转性能的弹性杆进行界面处理。 本章研究超细长弹性杆的表面图形描述问题，包括曲面的空间形状和弹性杆的 扭转。这里我们设弹性杆是由圆形的平面刚性区域，其中心沿着一条空间曲线运动 而形成，运动过程中圆形平面区域始终垂直于曲线。在这种前提下，我们导出了描 述弹性杆曲面的微分，积分方程组，给出了求解这一方程组的数值方法，并利用给出 的方程和算法分析了具体问题。 2 1 弯扭线的提出嘲 仿真中，为区分弹性杆的扰性线( 中心线) 和形象的表达弹性杆的扭转的概念， 本文提出了弯扭线的概念：即尚未弯曲与扭转情况下与中心线平行的直杆边缘“母 线”在扭转与弯曲后的形状。此“弯扭线”与弯曲后直杆的中心线形成的对比，充 分表达了围绕中心线缠绕与扭转的意思。根据文献“弹性杆的非线性力学”提出的 弯扭度的概念，本文采用弯扭线的概念，来描述d n a 弹性细杆在空间弯扭后的拓 扑结构。 现研究以长度工的细杆，将弹性杆横截面的几何中心连接成不自交的空间曲线 c l 称为中心线，以空间固定点为原点建立惯性坐标系o 一孝，7 f ，以细杆端点乩为原 点沿中心线建立弧坐标s ，并在曲线c l 上以任意点，为原点建立主轴坐标系 p x y z ，其中：轴为中心线的切线，也是杆的横截面的法线。当p 点沿中心线c 以 单位速度朝弧坐标j 的正向运动时，截面以角速度国相对o 一勿f 转动。显然，截面 内任意矢量的端点也随着p 点移动和绕切线转动而在空间划出一条轨迹，形成一条 缠绕的空间曲线，我们称为“弯扭线” 第二章弹性细杆的数值仿真 2 2 旋转矩阵 图2 1 弯扭线示意图 设0 - x t y t a t 和d 一以y i 孔是以d 为同一原点的不同坐标系，对应的基向量分别 为e ，和e k ，则同一矢量4 可以用两种不同的基表示出来 口= 4 “= 气t 童御。 其中口 “。口为向量口的坐标阵列。将上式后一等号的两边用e t 点乘，得到： 口 f ，= 以4 n 其中a a 为3 x 3 标量矩阵，定义为 4 。= f ，= 0 ，) p 一( 1 ) 由z - ( 1 ) n - - 以判断，相同基之间的旋转矩阵为三阶单位阵，1 5a n ( 或钆) = e 并且 a 2 = 山= 镌 与实际情况相符合。 我们用以下两个算例来实现旋转矩阵的应用。 f c o s t l 例1 ：螺线，( ，) = ls i n tl i r j 螺线在t = 0 点基本三棱形的三个基向量分别为 嘞= ( 。孚孚) ，风= t 一，。啪h = ( 。一孚孚 在任意点，的本三棱形的三个基向量分别为 删= ( 一譬如，譬瞄，孚 青岛大学硕士学位论文 旋转矩阵为 以= ( d ( ，) ) “= p ( o = ( - c o s t - s i n t0 ) = f 孚如，一孚删 压 一一s m f 2 c o s t 压 - - - 一s i n t 2 其中的值为螺线在r = o 点基本三棱形的第f 个基向量和螺线在任意点，的基本三 棱形的第_ ，个基向量的点积( i , j = l ，2 ，3 ) 而 s i n t 2 以叫岍i 孚伽， s i n t 2 互 一- 一c o s t 2 一s i n t j - - - 一c o s t 2 然后利用m a t l a b 作图，可以得到如下的图形。 s i n t 2 压 一- - - 一c o s t 2 s l n r 2 图2 2 带扭转螺旋杆示意图 例2 ：考虑直角坐标系口一习留的d x y 平面上的圆c ： ，c r ，= ( ：芋； 。r z 万 选择口一彬空间中的小圆环。 1 4 。卜2 m 。一2孚争等 。卜2 一2学参半 第二章弹性细杆的数值仿真 。 f 1 0 阳1 砸) = i o i ， 。 【s i n 0j 下面考虑两种情况： ( 1 ) 让小圆环c 的圆心始终在大圆c 上，且在任意时刻，大圆c 都垂直小圆c 所在平 面，小圆环在绕坐标系原点o 沿大圆转动的同时自身不发生任何转动，此过程用到 一个旋转矩阵( 图2 3 ( a ) 所示) ( 2 ) 让小圆环c 的圆心始终在大圆c 上，且在任意时刻，大圆c 都垂直小圆f 所在平 面，小圆环在绕坐标系原点o 沿大圆转动的任意时刻自身也以一定的速度绕大圆在 该时刻的切线发生转动，此过程用到了两个旋转矩阵( 图2 3 ( b ) 所示) 。 图2 3 ( a ) 图2 3 ( b ) 一一 图2 4 ( a ) 小豳绕大圆旋转运动表面示意图 2 3 刚体的有限转动的欧拉定理 图2 4 ( b ) 小圆环绕大圆旋转且自身扭转表面示意图 d n a 弹性杆模型以其极端细长性和超大变形而完全不同于传统弹性力学的研 究对象，也给d n a 结构的图形描绘和动力学性质的数值仿真也增添了困难。在理 论问题中，人们往往将弹性杆结构模型视为剐体有限转动模型。 青岛大学硕士学位论文 定理内容：刚体绕定点的任意有限转动可由绕过该点的某根轴一次有限动某个有限 角度来实现。 将刚体上的不动点记为0 ，如图所示过该点建立刚体的连体基和考察刚体运动 的参考基，分别记为 气= ，y ，矿) f ，= x r ) ，) 图2 5 刚体连体基和参考基示意图 刚体在参考基p 上的姿态与该刚体连体基p 6 相对于参考基p 的姿态是一致。它可以 用基p 6 相对于基e 7 的方向余弦阵4 。来描述， 小鞋批萎芬j | ； i a 咕一a ，1 = 一( 4 l + 彳2 2 + 彳站) 矛+ ( 4 i l + a 趋+ 4 ) 名一l = 0 由此可得到如下结论：对于任意两个基p 7 与基矿存在一个矢量，它在两基的 1 6 第二章弹性细杆的数值仿真 坐标阵相等。此结论也可理解为将矢量p 作为一个旋转轴，基p 。是基f 绕p 转过一 个有限角度后到达的新的方位。考虑到此矢量的存在性，可得定理结论。将单位矢 量p 称为由基到基矿一次转动矢量。转过的有限角称为一次转动角，记为，o 。从 上面分析不难看出，刚体相对于基f 7 的不同姿态均可绕相应的一次转动矢量和作相 应的一次转动角来实现，也就是说刚体的不同姿态对应不同的一次转动矢量和一次 转动角 2 4 数值离散与数值方法的应用 2 4 1a d d s 方法和递推方法 设0 一勿f 是惯性坐标系，弹性秆轴心曲线为空间中的2 阶以上的光滑空间曲 线 ，= ，( j )2 - ( 2 ) 其中是j 弧长参数。在曲线上任一点，建立曲线的f r e n e t 坐标系p n b t ，其中三 个基矢量n 、占分别是曲线2 - ( 2 ) 的切向，主法线方向和副法线方向的单位矢量。 借助于k i r c 砌o f f 的动力学比拟方法，我们把超细长弹性杆的表面看作一条刚性封闭 的平面曲线点集j 沿着空间曲线r ( s ) 移动和旋转得到( 图2 6 ) 。按照k i r c _ i l i ，o f f 弹性杆 的假定，我们把弧长j 看作拟时间坐标，即将弹性杆曲面沿着弧坐标j 变化看作一条 封阕的平面刚性曲线j 随。时闻”s 的匀速运动。设x o ) 是j 在“时刻”s 的位置， 它可以看作由点集r 运动得到。 ； 图2 6 ：弹性杆示意圈 1 7 青岛大学硕士学位论文 我们把集合与某元素的差定义为集合中每个元素与该元素差的集合。则 x ( s ) = o ( s x x ( o ) 一，( 0 ) ) + ，( o ) + ( r ( s ) 一，( 0 ) ) = q ( j ) ( x ( o ) 一，( o ) ) + ，( j )2 - ( 3 ) 其中 x ( o ) = f q ( s ) 是由f 一，( 0 ) 所在平面到x ( o ) 一，( o ) 所在平面的旋转矩阵n 8 1 。q ( 0 ) = j 。在2 ( 3 ) 式中对s 微分得到 i d x = 警( 猁，( 。) ) 川s ) 2 - ( 4 ) 注意至l jq ( s ) 是正交矩阵，利用q 7 左乘方程组2 - ( 3 ) 得到 x ( o ) 一r ( o ) = q ( j ) 7 ( x ( j ) - r ( s ) )2 - ( 5 ) 由于 一 ，( s ) 一，( o ) = 【r ( t ) d t 2 - ( 6 ) 将2 - ( 5 ) 式和2 ( 6 ) 式代入2 式得到求解曲面网线x ( s ) 的微分方程初值问题 警= 罢) r ( 耶) _ r w 瑚+ ，( o ) ) + 趴砂2 仰 ix ( o ) = f c 考虑弧坐标s 的等距节点= 地s 七= 0 , 1 ，2 ，。设鼍为x ( s d 的数值近似，并 记 f ( x ，。= 辈q ( s ) 7 ( x ( d f r ( ，) 出+ ，( o ) ) + r ( s ) 2。-(8)ls 由 利用a d a m e s 方法1 得到方程组2 - ( 7 ) 的数值离散公式 取= f ( x i ，) 一d 出q q ( ) 7 ( x ( ) 一只+ r ( o ) ) + t ( ) 2 - ( 1 0 ) 其中 五= f l r 西 2 - ( 1 1 ) 2 ( 1 1 ) 式中的积分项最每步需要利用数值积分计算，随着屯的增大计算量会大大增 加，这在微分方程组的边值问题中是较大的计算负担。为了减少计算量，可以考虑 利用以下迭代“o 的方法：由于 只“= r “r ( t ) d t = 忍+ e “r ( t ) d t 2 - ( 1 2 ) 1 8 锣 2 +， = i 吖b只 ，舯 厶 = 吖k吩 。川 中其 第二章弹性细杆的数值仿真 利用梯形公式近似2 - ( 1 2 ) 式右端积分项得到递推公式 最+ im 最+ 等( f ( 以) + r ( s ) )2 - 0 3 ) 利用2 - 0 3 ) 式计算，计算效率和精度都比直接利用梯形公式计算只好得多。 下面利用上文中的结果给出弹性杆的图形描述。为简单起见，考虑圆截面弹性杆。 将分布力，在在f r e n e t 标架下表示为 f = f 。n + f b b 在k i r c h h o f f 的圆截面弹性杆的简化假设下，利用超细长圆截面杆的s c h r o d i n g e r 方程得到曲率k 和挠率t 的微分方程组： 窘哪一) 2 _ c l t + 了t - 3 = 鲁 2 - 0 4 ) 誓+ 2 誓c r 一争= 鲁 2 - ( 1 5 ) 出出、2 7 彳 其中 hm 2 q 5 j 一了 玛h 是材料和运动有关的常数。 将方程2 - 0 4 ) 和方程2 - ( 1 5 ) 的数值解代) 、_ f r e n e t - s e r r e t 方程 盟：峦一玎 凼 坚：捌 出 塑：删 d s 2 - ( 1 6 、 解得正 丑。带入算法2 - 0 0 ) 2 - 0 1 ) 2 - ( 1 3 ) 即可解出弹性杆的曲面图形。 注意2 - ( 1 4 ) 2 - ( 1 5 ) 秉12 - ( 1 6 ) 式是边值问题，通常采用试射法或差分法计算。而2 - ( 1 0 ) 圣= 1 9 青岛大学硕士学位论文 圣( 力= 垂( o ) + f | 一f ( f ) o o t ) d t r 0 f ( f ) 一r ( f ) 、 o ( s t + 1 ) 却( 啪= 刊町00 卜( t ) d t * ( ii i + i 垂c ，+ ii 。庐c ) 1一尘f 鱼r 22 _ a s f 10 一- a s r 01 一i r f 川1 i 丑。i _ i 瓦+ j 1 a s f 一a s 茁 一尘r1o 了a s r ol 2 - ( 1 7 ) 2 4 1 8 ) f 肌1 2 q 9 ) 差分格式2 - 0 9 ) 是隐式二阶方法，计算稳定，计算量也很少。如果需要高阶精度 的算法，可以利用2 - ( 1 7 ) 式类似地导出r u n g e k u t t a 型方法。 2 4 2 四元数形式的微分代数方程表示的模型 超细长弹性杆动力学模型在d n a 分子结构模型的平衡、稳定性等问题的研究 中有重要的应用。为了便于其数值仿真，图形后处理以及研究表面接触等问题的处 理，需要建立弹性杆的表面模型。 四元数在刚体及其有限转动的研究中有重要的优势，文章利用k i r c h h o f f 弹性杆 模型的动力学比拟技巧，建立了由四元数的描述超长弹性杆曲面的微分积分方程 组。将惯性坐标系0 一手穆f 平移到曲线，( s ) 上任意点户为原点后的位置记为p x y z 。 设坐标系p x y z 绕过p 点的某个瞬时轴日转过矿角后与在点尸建立的f r e n e t 坐标 系p - n b t 完全重合。根据刚体有限转动定理，上述旋转过程是可以实现