（一般力学与力学基础专业论文）作为dna力学模型的弹性细杆的数值仿真.pdf

作为 d n a力学模型的弹性细杆的数值仿真 摘 要 本文首先对作为 dna 力学模型的弹性细杆静力学的基本原理进行了研究 采用打靶法与能量法的结合对弹性细杆的作用机制进行深入探讨即将传统 的解析研究与数值研究结合起来 为了进一步研究弹性细杆的力学性质本文将采用 euler 四元数表示的 kirchhoff 方程来研究受力挤压作用下的弹性细杆的拓扑构形将得到的微分 方程与约束条件组成微分代数方程后再转化为微分方程规范形式以便求解 为满足边界条件应用数值打靶法求解边值条件并将弹性细杆在力作用下 的拉压过程用 matlab 仿真出来.同时对由于误差导致的违约现象进行处理并 针对欧拉参数的特征选取合适的修正系数以保持方程的稳定性 接着将 dna 的双螺旋结构模拟成具有初曲率的弹性细杆并考虑弹性细 杆的自接触采用有限元离散方法对弹性细杆进行分割运用变分法对细杆 的势能求变分导出离散化的非线性方程考虑到碱基平面姿态参数的守恒 增加带有惩罚因子的约束方程组成系统方程为求解此非线性方程采用 levenberg-marquart 算法实施数值计算非线性方程迭代的初始值由打靶法定 出逐步改变初始曲率得出不同的空间拓扑形状 最后对传统的 kirchhoff 解析解给出仿真通用算法及其程序对欧拉四 元数形式的数值解给出通用后处理仿真程序并采用飘带和弯曲管的形状来 直观显示 dna分子的形状并讨论了不同的初曲率对 dna 三维空间拓扑构 形的决定性影响作用从而说明处于不同浓度生理盐溶液中的 dna 分子会 产生不同的分子构像这满足了不同生物学过程的需要 总之根据作为 dna力学模型的弹性细杆的模型可预测 dna超螺旋性质 的变化这对研究 dna 在复制转录重组修复中发生的性质改变具有 一定的现实意义 关键词 dna弹性细杆 kirchhoff方程 euler四元数 levenberg-marquart matlab 超螺旋 自接触 拓扑构形 违约修正 弯扭线 the numerical simulation of a thin elastic rod as the model of dna abstract in this paper, we study the postulate of statics model of thin elastic rods as the model of dna. adopt the method of combining the shooting technique with energy approach to further probe into the mechanism of action on thin elastic rods, that is, connect traditional analytical research to numerical research. to further study mechanical property of thin elastic rod, this paper will apply the kichhoff equation which takes the form of euler quaternion and study the topological configurations of thin elastic rod which is under the extrusion of pressure. add the constraint condition to differential equation we can get differential-algebraic equation (dae). easily to be solved, dae will be transformed into criterion form. to satisfy boundary condition, we apply shooting technique to solve boundary value condition, and imitate the procedure of pulling and pressing of thin elastic rod suffering strength. simultaneously deal with the constraint stabilization phenomenon which is caused by error, and according to the euler quaternion character, we select proper corrector coefficient to retain the stability of pure differential equation. next, the double helical polymer is modeled as a thin elastic rod with intrinsic curvature, a finite element analysis based on principle of potential energy minimum is formulated so as to treat self-contact. after discretization of the thin elastic rod, apply variational principle to potential energy of rod and deduce the equation of nonlinear discretization. take the conservation of attitude parameter on the base plane surface into account, add the constraint equation with penalty factor so as to form into system equation. to solve it, the levenberg-marquart arithmetic is introduced to numerical computation. initial value of nonlinear equation iteration can be determined by previous shooting technique on euler-kirchhoff equation, alter the intrinsic curvature step by step, we can figure out the spatial topological configurations. finally, the universal algorithm and program is put forward in order to simulate the traditional kirchhoff analytical solution, and also outline universa l postprocessor to emulate the numerical solution in the form of euler quaternion. in the shape of ribbon and tube with curvature, intuitively display the shape of a dna molecule. we discuss that different intrinsic curvature can be have an crucial affection on the spatial topological configurations, thereby explain the phenomenon which dna polymer under physiological salt solution with different concentration come into being different molecular conformation, which satisfy different need of biological process. in conclusion, according to the statics model of thin elastic rod as the model of dna, people can predict the change of property of dna supercoiling, which have an actual significance during various biochemical processes, including replication, transcription, recombination and restoration. key words: thin elastic rod of dna, kirchhoff equation, euler quaternion, constraint violation correction, levenberg-marquart, matlab, supercoiling, self-contact, topological configurations, crankle line 上海交通大学 学位论文原创性声明 本人郑重声明所呈交的学位论文是本人在导师的指导下 独立进行研究工作所取得的成果除文中已经注明引用的内容外 本论文不包含任何其他个人或集体已经发表或撰写过的作品成果 对本文的研究做出重要贡献的个人和集体均已在文中以明确方式 标明本人完全意识到本声明的法律结果由本人承担 学位论文作者签名黄磊 日期 2005 年 2月 1日 上海交通大学 学位论文版权使用授权书 本学位论文作者完全了解学校有关保留使用学位论文的规 定同意学校保留并向国家有关部门或机构送交论文的复印件和电 子版允许论文被查阅和借阅本人授权上海交通大学可以将本学 位论文的全部或部分内容编入有关数据库进行检索可以采用影 印缩印或扫描等复制手段保存和汇编本学位论文 保密在 年解密后适用本授权书 本学位论文属于 不保密 请在以上方框内打 学位论文作者签名 黄磊 指导教师签名 包光伟 日期 2 0 0 5 年 2 月 1 日 日期2 0 0 5 年 2 月 1 日 主要符号表 主要符号表 s 弧坐标 e 杆的杨氏模量 g 杆的剪切模量 n 有限元划分的单元数 l 弹性细杆的总长度 a 弹性细杆的截面半径 a,b 截面绕主轴x和y的抗弯刚度 c 杆截面绕主轴z的抗扭刚度 wr,tw,lk缠绕数、扭转数、联系数 x i, y i 截面相对x轴和y轴的惯性矩 z i 截面相对z轴的极惯性矩 拉格朗日乘子 截面的弯扭度 0 杆截面初始曲率和挠率 f v 杆平衡状态下的截面作用主矢 m v 杆平衡状态下的截面作用主矩 rb a, br a 连体基与参考基的方向余弦阵 j 本构关系中的惯量张量坐标阵 r 1 e, r 2 e, r 3 e 空间固定的直角坐标系 b 1 e, b 2 e, b 3 e 截面转动的主轴坐标系 r 弹性细杆的中心线 c r 弹性细杆的弯扭线 q 表征截面姿态的欧拉四元数 q 欧拉四元数的一阶导数 l,r 欧拉四元数的相关矩阵 t,s 四元数组成的转换矩阵 el e 弹性细杆的弹性应变能 u 杆单元间纳米力所作总功 j x 弹性细杆的单元节点变量 x 节点变量组集的总体变量 1 , 2 杆端部约束能 1 , 2 端部约束的拉格朗日乘子 l r 细杆终端的位置向量 l 细杆端部角位移向量 k 组集的系统应变能刚度阵 p 组集的系统应变能力阵 b 组集的系统自接触刚度阵 0 l 有限元划分的长度矩阵 4 e 44 单位矩阵 dk 非线性方程的惩罚因子 ， 违约修正振荡参数 ， 线性插值的自然坐标 c 扭转数离散的转换矩阵 lm levenberg-marquart 算法 d 自接触势能积分面域 第一章 绪论 1 第一章 绪论 1 1 弹性细杆平衡问题的研究历史 早在 1730年daniel bernoulli和 euler就开展了对弹性细杆力学问题的 研究kirchhoff于 1859利用弹性细杆静力学与刚体动力学的相似性建立了描述 弹性圆截面细杆在无轴向变形无体积力杆间接触力情况下静力学的动力学 相似理论基础弹性细杆非线性力学因在分子生物学中成功地模拟 dna等一 类大分子的静态平衡而重新受到关注形成经典力学与分子生物学交叉的新的 研究领域 从结构力学到生物化学以及生物工程弹性细杆的平衡和稳定均有 着广泛的应用背景如卷曲的海底电缆与井壁有接触的钻杆动态断裂的绳 索细丝结构的生物分子细菌纤维螺旋状攀藤植物的细茎人体毛细血管 都可以看作为弹性细杆的力学模型kirchhoff模型为研究弹性细杆的静力学与 动力学提供了一个有效的途径 在分子生物学领域中自从 20世纪中期 watson和 crick提出 dna分子 的双螺旋三维结构模型以来关于 dna的基础理论研究不断取得突破以遗 传工程为首的生物技术得到迅速发展dna双螺旋结构由两条螺旋形核苷酸链 和联系二螺旋线的碱基组成为有可能研究这种三维结构的几何形态及其稳定 性可以简化为具有初始扭率的圆截面弹性细杆以这种力学模型为基础经 典力学的基本原理和方法都能得到充分应用所需的物理常数如杨氏模量 泊松比截面的抗弯刚度以及抗扭刚度均可由试验来确定于是分子生物学与力 学产生了一个结合点自 70年代以来关于 dna力学模型的研究逐渐形成一 支交叉学科分支即将传统的力学研究方法与分子生物学的试验研究方法结合 讨论 dnarna以及细菌纤维等各种生物大分子链的几何形状和稳定性 一般情况下dna分子被模拟成均质的各向同性的线性圆截面弹性 细杆关于 dna弹性杆的研究已经有大量文献发表在物1理和化学物理物理 化学生物化学和高分子化学以及分子生物学等学科的刊物上其中大多数工 第一章 绪论 2 作是以 kirchhoff理论为基础的解析研究如:大尺寸结构弹性模型,受约束圆截 面 dna的扭曲与螺旋, dna双螺旋的扭转刚度与超螺旋统计力学,高超螺旋 dna弹性模型,环形 dna的超螺旋转变, 环形 dna的三维形状, 基因力学模型 考虑到杆间自接触的如dna超螺旋的有限元分析,dna超螺旋自接触模型 有限长 kirchhoff弹性杆的边值问题国内最新讨论了受曲面约束弹性细杆的平 衡问题从动力学的观点讨论了弹性细杆的平衡稳定性 人体中的 dna可长达 12 m ,而分子的螺旋半径约为 2 nm ,因此弹性细 杆所描述的是不同于传统弹性力学的研究对象:一定刚性而又可以卷曲的 dna 分子 dna的分子构像依赖于碱基序列拓扑限制不同的外界环境及蛋白 质的相互作用这种构像多态性与碱基序列本身的多态性一样对于核酸的功能 有重要作用为此本文弹性细杆的研究模型是通过有限元来离散 dna分子 结构着重从模拟碱基平面姿态的欧拉四元数入手来建立近似的数值计算模型 另外碱基对构成的平面与螺旋轴近似垂直螺旋轴穿过碱基平面这也是将 dna分子描述成弹性细杆的一个重要物理基础 1 2 弹性细杆课题研究的意义 1 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在复制转录重组修复中不断改变着自身的拓 第一章 绪论 3 扑 构型引起超螺旋参数扭转数缠绕数 的变化 因此可以通过改变 扭转数来获取具有初曲率的 d n a分子稳定状态的不同构象以满足不同生 物学过程的需要 1 3 数值计算的现状与发展 关于 d n a弹性细杆的数值计算如 y a n g ( 1 9 9 3 ) 4 基于 k i r c h h o f f理论的有限 元方讨论自接触的问题 y a n g , z h a n g , t o b i a s , o l s o n ( 1 9 9 4 ) , s c h l i c k ( 1 9 9 2 , 2 0 0 0 ) 和 w e s t c o t t ( 1 9 9 7 ) 2 针对 d n a超螺旋的自接触模型a l e x a n d r e ( 2 0 0 3 ) 1 打靶法 求解受边界约束的 k i r c h h o f f 弹性细杆的平衡构形 关于 d n a在平衡状态附近的波动以及振动的模态问题有待于进一步的讨 论 1 4 本文主要研究目标及内容安排 为方便研究首先从经典的 k i r c h h o f f 方程推导出欧拉参数形式的微分方 程进一步深入探讨 dna分子对多重蛋白结构的缠绕性质建立了 dna弹性 细杆受力挤压情况下的一般模型对受到边界条件约束的生物分子建立的打靶 法不仅有助于能量法的分析研究而且对于利用光学捕获手段操作 dna单分 子有着重要的指导意义 打靶法的结论将能说明 是在其结构上引起扭曲的 dna 损伤是可以修复的d n a 链与蛋白质的缠绕对转录过程的影响 接着将 d n a的双螺旋结构模拟成具有初始曲率的弹性杆并考虑弹性杆的 应变和自接触采用有限元离散方法对弹性杆进行分割并用变分法对弹性杆 的势能求变分导出离散化的非线性代数方程利用序列决定 d n a结构的性质 导出研究对象的初始曲率根据计算的结果着重研究闭环 d n a的超螺旋性质 直观地用飘带和流管来图示 d n a分子的形状能量法的结论将试着解释核质体 的超螺旋化 dna局部性质的变化可能导致生物学过程规律的改变这与配合基或蛋白 质对 dna的包装作用相关包装过程将引起共价键或氢键的联结或断裂从 而导致分子局部弹性常数的改变由此产生分子构像的变化因此有限元的计 第一章 绪论 4 算方法为研究这种局部性质变化引起总体构像改变的特性提供了可行的途径 采取 m a t l a b 编写计算程序前处理输入欧拉角度等参数后处理产生三个 描述超螺旋参数值输出图形并与试验结果分析比较 1 数值计算方法最小势能原理 采用打靶法对变形后的 k i r r c h h o f f方程满足的边界条件实施计算采用 欧拉四元数作为变量建立势能最小化非线性的方程组并采用非线性最小二 乘法中的经典方法实施数值计算有限元离散方法为计算的常用工具 2 仿真程序( 开环闭环 d n a ) 通过分别对研究目标的数值结果的仿真编制一套针对开环与闭环情 况的数值仿真程序采用带有网格弯曲弹性杆和飘带分别模拟所需要的结果 并对传统的 k i r r c h h o f f 解析解给出通用的后处理算法 3 前处理输入描述碱基平面的方位角等参数由序列决定结构计算产生初 始 曲率同样也可以由打靶法的结果作为能量法迭代计算的初始值并在计算 结束后作相关的参数与物理性质的分析比较 4 后处理计算超螺旋几何参数( 联系数扭转数缠绕数) 输出 d n a 构型 分别对各参数的计算公式进行变化以适应计算模型的需要并进一步利 用超螺旋参数来讨论分子的平衡结构与相关物理现象相比较解释 5 与相关结果对照并分析比较 根据实际的情况将数值结果的结论与已经存在的关于 dna 分子的理 论知识对比解释在此基础上进一步根据分析的结论推断新的关于 dna 弹 性性质并预测下一步的研究方向与目标 以上便是本文研究的主要目标和内容全文各章节的具体内容安排如下 第一章为绪论主要介绍弹性细杆平衡问题的研究历史弹性细杆课题研 究 的意义数值计算的现状与发展并提出对弹性细杆研究的新方向与期望 第一章 绪论 5 第二章为弹性细杆弯扭的 kirchhoff方程的违约校正求解本章主要介绍欧 拉四元数形式的 kirchhoff 方程的推导对组集的系统微分代数方程进行变换 得到纯微分方程的形式并在此基础上违约修正利用边界条件打靶从而建 立了弹性细杆在拉压情况下的力学模型算例能从一定程度上较好地解释 dna 分子的包装过程针对弯扭的数值模拟讨论染色质的基本单位核小体的作用机 理 第三章为弹性细杆平衡的变分极值原理主要介绍弹性细杆势能变分原理 利用有限元划分弹性细杆采用欧拉四元数作为节点坐标并组集成总体变量 对弹性细杆的应变能作用于弹性细杆外力作的总功 van der waal 势能弹性 细杆端部的几何约束等等进行变分导出系统的非线性方程并将欧拉参数守 恒条件一同作为系统方程的一部分进行迭代求解数值计算过程中采用非线性 最小二乘的 levenberg-marquart 方法为加快方程的收敛速度在此基础上采 用正切调谐因子调整迭代的步长 第四章讨论了 dna 弹性细杆数值仿真问题首先介绍了数值仿真应用背 景 接着在弯扭度的基础上提出弯扭线的概念针对解析解的情况给出了基于向量 递推弯扭线通用算法并对混沌解进行了数值试验 第五章主要讨论了反映 dna弹性细杆超螺旋性质的重要参数扭转数联 系数缠绕数并对相关公式进行推导使得能够轻易地使用为计算和讨论 dna分子的拓扑构形提供了重要的参考依据 第六章为全文总结将指出本文的主要贡献及创新 第二章 弹性细杆弯扭的 kirchhoff 方程的违约校正求解 6 第二章 弹性细杆弯扭的 kirchhoff 方程的违约校正求解第二章 弹性细杆弯扭的 kirchhoff 方程的违约校正求解 为了进一步研究弹性细杆的力学性质，本章将采用 euler 四元数表示的 kirchhoff 方程来研究受力挤压作用下的弹性细杆的拓扑构形；将得到的微分方 程与约束条件组成微分代数方程后再转化为微分方程规范形式以便求解； 为满足 边界条件，应用数值打靶法求解边值条件，并将弹性细杆在力作用下的拉压过程 用 matlab 仿真出来.同时对由于误差导致的违约现象进行处理， 并针对欧拉参数 的特征，选取合适的修正系数以保持方程的稳定性。 dna 的复制与转录能力与 dna 分子局部的弹性力学性质有着密切的关系，细 长的弹性分子被挤压在核小体内，这将导致分子构形的变化；同样与蛋白质分子 的缠绕与挤压也将有同样的现象产生；为此，本文建立了 dna 弹性细杆受力挤压 情况下的一般模型。 首先建立欧拉参数形式的 kirchhoff 方程， 并与系统的约束方程组成微分代 数方程；为求解微分代数混合方程，采用增广法进一步将此混合方程转化为标准 的纯微分方程，为满足系统的边界条件，对纯微分方程进行打靶法，同时对由于 误差导致的违约进行修正；通过数值计算与仿真，希望能够进一步了解 dna 分子 的一些弹性力学性质。欧拉参数的无奇点性给方程的求解带来极大的方便；利用 欧拉参数表示的 kirchhoff 方程使得微分方程更加简洁与规范。 最后，在弹性细杆弯曲模型的基础上，将上述理论推广到弯扭的状况，并进 一步讨论弹性细杆的弯扭与 dna 分子包装作用之间的内在关系。 21 模型的提出与假设21 模型的提出与假设 现考虑一端部受力作用的简支简支细长弹性细杆的弯曲，如图 2-1 所示， 杆的截面几何中心连成的空间曲线c为杆的中心线， 以其中一个端点o为原点沿 中心线建立弧坐标s.则由曲线上任意一点p 相对于固定点o的矢径)(sr v 完全确 定曲线c的几何形状，杆的弧长为l，即为拉直状态下的长度。 第二章 弹性细杆弯扭的 kirchhoff 方程的违约校正求解 7 211 基本假定：211 基本假定： a) 杆为均匀各向异性，本构方程为线性关系； b) 杆等截面，截面沿主轴方向有相同几何尺寸； c) 忽略杆的体积力以及杆与杆之间的接触力； d) 杆存在着初始曲率与挠率； e) 杆的长度和曲率半径远远大于截面的尺寸； f) 杆无伸长或缩短，任意二截面沿中心线的 距离不变。 图 2-1 弹性细杆的空间姿态 fig 2-1 the spatial attitude of thin elastic rod 212 平衡方程212 平衡方程 如图 1 所示,以o点为原点建立空间固定的直角坐标系 )(o(参考基 r e v );以杆中心线上的任意点p为原点， 建立截面的主轴坐标系(xyzp)(连体基 b e v ),引入 v 为p点处刚性截面相对弧坐标s的角位移变化率； 设p点处截面上作 用的内力主矢和主矩分别为f v 和 m v ， 3 e v 为z轴的基矢量.利用微元体的平衡可以 导出弹性细杆的平衡方程为 31： 0 d d 0 d d 3 =+ = fe m f v v v v s s (1) 弹性细杆本构关系： )( 0 jm vv v = (2) j j 为惯量张量，其矩阵形式为： = c b a 00 00 00 j，且本构关系在截面连体基 b e v 上 的坐标阵为： = 3 03 2 02 1 01 00 00 00 m c b a b (3) 第二章 弹性细杆弯扭的 kirchhoff 方程的违约校正求解 8 a,b,c分别为截面绕x轴和y轴的抗弯刚度以及绕z轴的抗扭刚度， 可表示为： x eia=， y eib=， z gic =；其中e为材料的杨氏弹性模量，g为剪切弹性模 量，)1/(+=eg，为泊松比； x i， y i为截面相对x轴和y轴的惯性矩， z i为 截面相对z轴的极惯性矩。设圆截面的半径为a，则有： 4 4 a ii yx =， 2 4 a iz = (4) )( 00 s ii =) 3 , 2 , 1(=i为细杆在无力作用（包括力与力偶）状态下存在的弯扭度分 量 32，即杆存在着初始曲率和挠率。 由力平衡方程，得到截面作用力主矢为一常量： 0 ff vv = (5) 力矩平衡方程在截面连体基 b e v 上的表达式为： 0 d d 03 =+fem m v v v v v s (6) 即： 0)()( 0300 =+fe v vvvvvv jjjj (7) 其在截面连体基的矩阵式为： 0 )( )( 0300 =+fejj (8) 22 弹性细杆欧拉参数形式的 kirchhoff 方程22 弹性细杆欧拉参数形式的 kirchhoff 方程 根据刚体有限转动的欧拉定理 22定义欧拉四元数q ，导数为q，引入l与r () () () () = = 2/cos 2/sin 2/sin 2/sin 3 2 1 4 3 2 1 p p p q q q q q， = 4 3 2 1 q q q q 10= 则 levenberg-marquart 算法总体原则设计如下： )()()()( 1 1kkt k kktkk xfxjiuxjxjxx + += (62) = = 法：提高收敛速度接近 收敛接近最速下降法：保证 newtongaussnextuu xfxfif repeatuu xfxfif kk kk kk kk ;/ )()( ; )()( 1 1 在实际计算中，为了使得方程组的收敛进一步加快，在(62)中增加一正切松弛因 子以调谐方程迭代的步长大小，并获得满意的效果。 323 带惩罚因子的非线性方程的 lm 求解323 带惩罚因子的非线性方程的 lm 求解 本文在 lm 法的基础上加以改进，使得能快速收敛到一个存在解。考虑到非 线性方程与约束方程的耦合形式， 将约束方程乘以惩罚因子后与(59)合并成总体 的非线性方程组，如下所示： =+ = = +=+ 0 xxdk 0l 0rr qqpb)(kx i t i l t ) 1)(1( )( 2 2211 i tt (63) 此非线性方程组共有 5n+11 个方程， 可用于确定 4n+4 个节点变量、 6 个 lagrange 乘子 1 ， 2 和 n+1 个惩罚因子，方程组封闭。 3333 算例分析算例分析 第三章 弹性细杆平衡的变分极值原理 33 331 不考虑接触力的开环情况 331 不考虑接触力的开环情况 考虑弧长4=l(无量纲化)的弹性细杆， 细杆两端固定， 起点为o()t000， 利用打靶法得出的结果作为能量法迭代的初始值，结果如下： 类 型 端部受力 端部力矩 终端位移约束 0.40006= x f 1.0492= x m 0.68473= x r -0.40174= y f 0.53201= y m-0.68473= y r 打靶法 1=tw 2.5938= z f -3.1401= z m 3.2134= z r -0.063= x f -1.2885= x m 0.68473= x r 0.053634= y f0.84702= y m-0.68473= y r 能量法 -0.29419= z f -3.2212= z m 3.2134= z r 表 3-1 打靶法与能量法的参数比较 1 table 3-1 contrast of parameter between shooting technique and energy method （a） （b） 图 3-3 受边界条件限制的弹性细杆平衡构形 1 fig 3-3 the balance configuration of thin elastic rod which suffer boundary constraint 第三章 弹性细杆平衡的变分极值原理 34 从表 3-1 可以知道:能量法平衡时端部受力要远远小于打靶法平衡时端部受 力，势能最小状态下端部受力是打靶法的 15，13，11；端部力矩是打靶 法的 123，159，102. 类 型 端部受力 端部力矩 终端位移约束 2.5997= x f 0.020794= x m0.68593= x r -2.5877= y f 0.21537= y m-0.68399= y r 打靶法 2=tw 12.524= z f -6.2728= z m 3.3135= z r -0.043837= x f3.1787= x m 0.68593= x r 0.042455= y f-1.2339= y m-0.68399= y r 能量法 -0.23321= z f -6.2948= z m 3.3135= z r 表 3-2 打靶法与能量法的参数比较 2 table 3-2 contrast of parameter between shooting technique and energy method （a） （b） 图 3-4 受边界条件限制的弹性细杆平衡构形 2 fig 3-4 the balance configuration of thin elastic rod which suffer boundary constraint 第三章 弹性细杆平衡的变分极值原理 35 从表 3-2 可以知道:能量法平衡时端部受力要远远小于打靶法平衡时端部受 力，势能最小状态下端部受力是打靶法的 1.7，1.64，1.86；端部力矩是 打靶法的 15200，573，100. 332 开环结果分析与结论 332 开环结果分析与结论 从表 3-1 和表 3-2 中还可以知道，在能量最小状态下细杆端部受力与力矩 的方向可能也已经发生变化；在初始扭转量增加的情况下，端部受力减小更多； 而补偿的代价是扭矩的增加越多；另外，注意到由于端部固定不动，垂直端部的 轴向力矩几乎没有发生改变。 打靶法的物理意义是 dna 分子为实现某种生物学过程的需要而端部受到拉 压力和力矩作用下所达到的初始状态；能量法的物理意义则是在此初始状态下 dna 分子为实现总体势能局部最小而进一步到达的状态。从图 3-3 和 3-4 显然看 出：在打靶法的基础上得出的弹性细杆平衡构形并不稳定，经过能量法迭代计算 后得出的平衡构形比打靶法更为复杂：细杆端部受力减弱，力矩有所增加；螺旋 次数明显增多，弯曲与扭曲更为均匀。 文献 23中指出：当染色体复制时，细胞分裂成 2 个压缩的 dna 团块结构，称 为核质体，核质体最明显的特征之一是 dna 分子卷曲成螺旋形式，这是 dna 分子 超螺旋化的结果：即 dna 分子的末端以相反方向旋转，卷曲成螺旋形以减少张力 强度；显然从数值试验的结果与图形中可以验证这一点，即弹性细杆端部张力的 减弱与螺旋次数的增加的对比。 基因组的超螺旋化， 还意味着超螺旋产生的扭应力状态可以在一个结构内沿 着 dna 分子自由传递。图 3-4(a),(b)的对比分析可以看出：弹性细杆的螺旋结 构从局部扩展到整个结构中， 使得局部最小能量状态的弯曲与扭曲比初始状态更 加均匀；曲杆的中心线也向更加紧凑的方向靠拢，这对于说明长约 1mm 的 dna 分子如何折叠、凝聚以适应几个微米的细胞容量有一定的意义。 333 考虑接触力的闭环情况 333 考虑接触力的闭环情况 考虑弧长4=l(无量纲化)的弹性细杆， 细杆两端封闭， 起点为o()t000， 第三章 弹性细杆平衡的变分极值原理 36 同样利用打靶法得出的结果作为能量法迭代的初始值，结果如下： 类 型 端部受力 端部力矩 终端位移约束 -2.0786= x f -0.75415= x m 0= x r 1.256= y f -1.3797= y m 0= y r 打靶法 1.3279= z f -1.5704= z m 0= z r 0.086586= x f 0.00011959= x m0= x r -0.056443= y f 0.00077974= y m0= y r 能量法 5-0.0006699= z f-0.13007= z m 0= z r 表 3-3 打靶法与能量法的参数比较 3 table 3-3 contrast of parameter between shooting technique and energy method 图 3-5 受边界条件限制的弹性细杆平衡构形 3 fig 3-5 the balance configuration of thin elastic rod which suffer boundary constraint 从表 3-3 可以知道:能量法平衡时端部受力和力矩都要远远小于打靶法平衡 时端部受力，势能极小状态下端部受力是打靶法的 4.1，5，0.05；端部 力矩是打靶法的 0.0158，0.0493，8.283. 第三章 弹性细杆平衡的变分极值原理 37 334 闭环结果分析与结论 334 闭环结果分析与结论 从闭环的结果来看，能量法平衡时不仅端部受力减小很多，而且端部力矩 也减小很多； 弹性细杆的这种性质说明闭环平衡结构比开环平衡结构更容易使得 dna 分子的张力和力矩得到缓解， 而这些张力并不能释放到分子之外， 只能在 dna 内部采取分子扭曲的形式来补偿 43。从图 35 可以看出，在扭转数很小的情况 下，弹性细杆只需稍稍改变形状就可以使得内部的张力和力矩极小化。 【展望】 如果在计算中考虑初始曲率 0 的影响，可以通过非线性迭代逐步增加初始 弯曲曲率和初始扭转曲率的大小， 以获得具有初始曲率的 dna 弹性细杆的平衡构 形；在初始曲率增大到一定程度时候，整个闭环分子的构形将会产生麻花状的构 形，关于初曲率的工作由于时间的限制将在后续的研究工作中进一步讨论。 第四章 dna 弹性细杆数值仿真 38 第四章 dna 弹性细杆数值仿真第四章 dna 弹性细杆数值仿真 41 数值仿真应用背景41 数值仿真应用背景 自从 dna 三维螺旋结构模型提出以来，关于 dna 弹性细杆模型的研究工作 已经取得不少成果，其中大多数是以 kirchhoff 理论为基础的解析研究，另外是 考虑复杂情况的数值方法；不管采用何种方式的研究，结果都需要将 dna 的三维 状态模拟出来，这样才能直观地表达分子的空间拓扑构形，为了更加清楚表达 dna 螺旋的状态，需要用弹性细杆的扭曲来实现之，本文采用 matlab 编制计算 程序对具有弯曲和扭曲性能的弹性细杆进行界面处理。 仿真中，有使用 ribbon 带来表示细杆的弯曲与扭曲的变形 20，有采用带网 格的弯曲细杆来图示模型 18,19，但是至今为止尚未能给出这种后处理界面通用程 序和相应的算法，为区分挠性线（中心线） ，本文首先提出曲杆弯扭线的概念： 即尚未弯曲和挠曲情况下与中心线平行的直杆边缘直线在弯曲和挠曲后的形状， 充分表达了围绕中心线缠绕和扭转的意思。随后，为了能够编制较为通用的后处 理界面程序 fzdna，应用有限单元划分 dna 弹性细杆，采用向量单向递推的方法 来确定任一弯扭线的空间位置，该向量是确定弯扭线的关键因素；而通过对弯扭 线阵列后即可形成直观的界面处理程序，输出的图形为带有网格的弯曲杆。显然 ， 采用带有网格的弯曲弹性细杆来描述具有一定刚性而又可以卷曲的 dna 分子是 比较适合的，dna 分子的各种拓扑学状态能够在这里得到一定出程度上的反映。 本文所阐述的通用算法不仅适应于传统的 kirchhoff 弹性细杆，也适应于具 有初曲率以及考虑复杂杆间接触力的情况，甚至适应于更为复杂的情形；因为不 受任何假设条件的限制，因而具有较好的通用性能。 42 曲杆弯扭线的提出42 曲杆弯扭线的提出 根据“弹性细杆平衡的动态稳定性” 33一文提出的弯扭度概念，本文提出近 似的弯扭线的概念，以表征弹性细杆的空间拓扑构形。如图 4-1，现讨论一长度 第四章 dna 弹性细杆数值仿真 39 )(sr l p 0 p 弯扭线 2 c 中心线 1 c a z y x p a )(sr o 为l的弹性细杆， 将杆截面的几何中心连成的自不相交的空间曲线 1 c称为中心线 。以 0 p点为原点沿中心线建立弧坐标s，以空间固定点o建立惯性参考基 r e()o，并以 1 c上任意点 p为原点建立截面的主轴坐标 基 b e()xyzp，其中z轴为截 面的法线，与中心线的切线重 合。当p点沿中心线 1 c以单位 速度朝弧坐标s的正向运动时， 截面以角速度相对 r e转动。 图 4-1 弹性细杆弯扭线示意图 为截面的弯扭度，表示为： fig 4-1 the sketch of crankle curve of thin elastic rod bbb 332211 eee+= (1) 其中 3 称为杆的扭率，是截面相对于惯性参考基绕切线轴的变化率。显然，截 面内任意一矢量的端点也将随着p点移动和绕切线的转动而在空间划过一个轨 迹， 形成一条缠绕 1 c的空间曲线 2 c， 我们称之为弯扭线。 在截面内沿主轴y)( 2 b e 取一长度为a的微小线段，则可以由此确定弯扭线 2 c的空间位置坐标阵： brb ass 2 )()(earr+= (2） 其中 rb a为 b e 基相对于 r e 的方向余弦阵，)(sr为 1 c的坐标阵， b e2是主轴y的单位 向量坐标阵。 这样， 通过对弯扭线的阵列后， 划分网格并使用 surfsurf