（应用数学专业论文）超细长弹性杆的结构模型及数值仿真算法的研究.pdf

摘要 近年来，对作为d n a 宏观力学模型的超细长弹性杆的数学建模、数 值计算和数学性质的研究有重要的发展。本文在这些成果的基础上，针 对d n a 分子带有静电且静电斥力对其结构有较显著影响这一特点，在刚 性截面、线性本构关系等基本假定成立的条件下建立了存在静电斥力作 用的超细长弹性杆的k r i c h h o f f 方程和欧拉四元数的平衡方程。 由于所导出的数学模型是含有复杂积分项的微分积分方程组，而且 积分限是与未知量有关的函数，直接计算有困难且计算量较大。我们对 问题的数值计算进行了研究。对于初值问题，我们考虑利用迭代方法进 行求解，并建立递推公式简化了计算。对于边值问题，则利用有限元方 法进行了直接离散求解。 另外，对于弹性杆曲面方程，我们利用k r i c h h o f f 动力学模拟的技 巧，建立了弹性杆表面的曲面微分方程并给出了图形处理的计算实例和 结果分析。 作为对d n a 力学模型的弹性杆的模型的研究，我们期望我们的工作 可以在研究预测d n a 超螺旋性质的变化、静电斥力的影响等方面得到应 用。 关键词：超细长弹性杆；d n a ；静电斥力；数学建模；数值计算 a b s t r a e t i nr e c e n ty e a r s ，g r e a tp r o g r e s sh a sb e e nm a d ei nm o d e l i n g ，n u m e r i c a l s i m u l a t i n g a n da n a l y z i n ga s u p e r s l e n d e re l a s t i cr o da sa ni m p o r t a n t s t r u c t u r a lm o d e lf o rs t u d y i n gd n as t r u c t u r e b a s e do nt h ea c h i e v e m e n t sa n d t h ef a c tt h a tt h ee l e c t r o s t a t i cr e p u l s i o ni nd n ac a nn o t a b l ya f f e c t si t s s t r u c t u r e ，ak i r c h h o f fe q u a t i o na n di t se u l e rq u a t e r n i o nr e p r e s e n t a t i o na r e g i v e ni nt h i sp a p e r a ss t r u c t u r a lm o d e l sf o rs t u d y i n gt h es u p e rs l e n d e re l a s t i c r o d sw i t he l e c t r o s t a t i cr e p u l s i o n s i n c et h en e wm o d e l sa r ec o m p l i c a t e dd 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 n sw i t h u n k n o w ni n t e g r a lp a t h s ，c o m p u t a t i o nt e c h n i q u e sa r ed e s i g n e di nn u m e r i c a l s i m u l a t i o n t oi n i t i a lp r o b l e m s ，a ni t e r a t i v em e t h o di si n t r o d u c e dt od e a lw i t h u n k n o w ni n t e g r a lp a t h s ，a n dar e c u r s i v ef o r m u l ai sd e s i g n e dt or e d u c et h e a m o u n to fc o m p u t a t i o n t ob o u n d a r yp r o b l e m s ，a na l g o r i t h mb a s e do nf i n i t e e l e m e n tm e t h o di sg i v e n i na d d i t i o n ，b yu s i n gt h ek i r c h h o f f sd y n a m i c a la n a l o gt e c h n i q u e ，a d i f f e r e n t i a le q u a t i o ns y s t e mi ss e tu pt od e s c r i b et h es u r f a c eo ft h ee l a s t i c r o d a sam o d e ls t u d yf o rd n as t r u c t u r e s ，w ee x p e c to u rr e s u l t sc a nf i n d f u r t h e ra p p l i c a t i o n si nf o r e c a s t i n gt h es u p e r - h e l i xp r o p e r t i e sa n dt h ea f f e c t i o n o fe l e c t r o s t a t i cr e p u l s i o n k e yw o r d s ：s u p e rs l e n d e re l a s t i cr o d ；d n a ；e l e c t r o s t a t i cr e p u l s i o n ；n u m e r i c a l s i m u l a t i o n ；m a t h e m a t i c a lm o d e l l i n g 青岛大学硕士学位论文 主要符号表 彳杆绕x 轴的抗弯刚度 量中心线的副法线矢量 b 杆绕y 轴的抗弯刚度 c 中心线 c 杆的抗扭刚度 d 杆的刚度矩阵 e 单位阵 ，气，f r e n c t 坐标系( p - n b r ) 的基矢量 巳( j = l 2 ，3 ) ( p 叼功的基矢量 气，白，勺( d 一善印f ) 的基矢量 ，截面作用力主矢 刎9 截面作用力主矢的增量 ，复主矢，f = e + i e e ( _ ，= 1 , 2 , 3 ) f 对( p - x y z ) 的投影 ，单位长度接触力 ，复分布力，= z + 巩 石，磊，沿法线轴和副法线轴的投影 g 杆的剪切模量 ，。 h h = ( 2 h a ) 一加二 五，5 截面对工轴和y 轴的惯性矩 五截面的极惯性矩 ，s c h r 6 d i 啊方程的初积分常数 ，单位长度杆的惯量张量 k 挠性线相对弧坐标的角频率 三杆的长度 工7 有限元法第歹单元的长度 吖截面作用力主矩 a m 截面作用力主矩的增量 肘复弯矩，m = 肘+ i m 2 m 单位长度分布力 中心线的法线矢量 3 6 一截面边界线的外法线基矢量 d 固定参考点 ，中心线上任意点 p 户点的邻近点 凡杆的起始端 凡杆的终止端 只u = o l ，功有限元法的节点 层杆的杨氏模量 q x i = 1 , 2 , 3 a ) q i = d q i 出 q a i = 1 2 ，3 ，4 ) 广义坐标( 欧拉参数) ，p 点相对d 点的矢径 s 横截面面积 s 7 第歹单元的的截面积 s 弧坐标 r 中心线的切线矢量 x , y ，：截面主轴 葺有限元法第i 个变量 而j 变量毛在节点弓处的值 u = l 国f 沿x , y 轴的投影 口，厉，f 轴相对( p 叼功的方向余弦 出矢径，的增量 出弧坐标增量 口欧拉角中的章动角 六玎，f 定坐标轴 善复曲率，手= q + i 吐 刁= 声+ i ， = m z + i a 3 ，杆的泊松比 r 中心线的曲率 pq 相对p 的矢径 主要符号表 p 曲率半径 p 杆的密度 i 角位移矢量 ，有限转动角位移 z 截面的相对扭角 f 中心线的挠率 口弯扭度 ， 仉d a r b o u x 矢量 国，【，= l ，2 ，3 ) 毋对( p 搁的投影 执秆的扭率 国；u = 1 ，2 3 ) m ，的原始值 m ：杆的原始扭率 青岛大学硕士学位论文 学位论文独创性声明 本人声明，所呈交的学位论文系本人在导师指导下独立完成的研究成果文中依 法引用他人的成果，均已做出明确标注或得到许可论文内容未包含法律意义上已 属于他人的任何形式的研究成果，也不包含本人巴用于其他学位申请的论文或成果 本人如违反上述声明，愿意承担由此引发的一切责任和后果 论文作者签名： 孙楞 4 0 日期： 叼年r 月j 1 日 学位论文知识产权权属声明 学位论文知识产权权属声明 本人在导师指导下所完成的学位论文及相关的职务作品，知识产权归属学校。学 校享有以任何方式发表、复制、公开阅览、借阋以及申请专利等权利。本人离校后 发表或使用学位论文或与该论文直接相关的学术论文或成果时，署名单位仍然为青 岛大学。 本学位论文属于： 保密口，在年解密后适用于本声明。 不保密口。 ( 请在以上方框内打“ ” ) 论文作者签名：莩r 永纵日期：川年，月；i 日 导师签名：翘乡佳力应日期：) 口0 7 年f 月日 4 1 引言 引言 弹性细杆的动力学模型是一个古老的问题，有关研究可以追溯到2 0 0 多年前 d a n i e lb e r n o u l l i 和e u l e r 的工作1 8 5 9 年k i r c h h o f f 对不考虑轴向变形的圆形 截面的细杆，在刚性截面假设和忽略体积力的条件下，利用弹性细杆的平衡微分方 程与刚体定点转动微分方程之间的相似性，提出了弹性细杆平衡的动力学比拟理论 日，建立了独特的建模和分析方法。1 0 0 多年来，在k i r c h h o f f 的理论和方法基础 上发展出研究弹性细杆动力学问题的数学建模和动力学分析的系统的理论与方法， 并在许多实际问题如海底电缆、纤维等的数学建模、结构稳定性分析和数值计算等 方面得到广泛的应用啪嘲嗍而且，由于数学微分方程模型的相似性，这一研究 也与其他的领域发生联系，在卷浪的传播、太阳黑子的形成等的研究中得到应用n 1 。 近3 0 年来，k i r c h h o f f 弹性杆理论和动力学比拟的建模方法在分子生物学的 研究中得到应用和发展。如大肠杆菌，肺的纤毛和d m 等都可以借助上述方法进行 研究。以d n a 为例，自2 0 世纪中期w a t s o n 和c r i c k 提出了d n a 分子的双螺旋三维 结构以来，关于d n a 分子的结构理论研究不断被突破，用具有原始扭率的弹性细杆 作为d n 的宏观力学模型的理论研究得到了实验的肯定，d n a 弹性杆模型中的重要 物理参数如杨氏模量、泊松比、抗弯和抗扭刚度等已经借助生物学实验的手段和统 计物理的方法得到御刖。因此，宏观的建模方法和经典的分析工具可以在d n a 的研 究中得到充分的应用 、 d n a 分子的力学性质对其生物功能( 折叠成染色体超螺旋、作为基因转录、复 制的模板) 有重要的影响紧缩的d n r 超螺旋结构存储着生物学过程必需的能量和 信息，对复制、转录、修复等生物学过程的研究具有普遍的意义。因此，在近3 0 年 中，相关的研究受到重视，成为生物学、医学、物理学和应用数学学科交叉研究的 重要领域。 但是，研究d n a 的结构的弹性杆模型又有其特殊的内容。人体细胞的最大染色 体所含d n a 分子的螺旋直径约为2 n m ，而长度可达7 c m ，杆长为半径的3 5 1 0 ，倍 如此细长的分子链必须往复缠绕以保证能够被容纳在半径仅l o u m 的狭小的细胞核 空间中。因此，d n a 的弹性杆模型具有极端细长性和超大变形性，这是与经典的弹 性杆模型的研究对象完全不同的。这样的细长的大变形d n a 分子像毛线团一样往复 缠绕，给d n a 结构分析和数值仿真也提出了新的课题。为了区别这模型与传统的 的弹性杼模型，我们称之为超细长弹性杆 3 0 多年来，d n a 弹性杆的k i r c h h o f f 模型及其相应数学理论和方法的研究工作 取得了很多成果，w h i t e ( 1 9 6 9 ) n 町，f u l l e r ( 1 9 7 1 ) “”最早研究了d n a 双螺旋结构 的拓扑性质，引入了缠绕数作为三维空间曲线的几何不变量，并研究了相应的数学 青岛大学硕士学位论文 性质。c r i c k ( 1 9 7 6 ) “”讨论了描述d n a 双螺旋结构的连接数，扭转数、缠绕数及其 相互关系。b e n h a m ( 1 9 7 7 ，1 9 7 9 ) “”“1 和l eb r e t ( 1 9 7 8 ，1 9 7 9 ) “剐”1 等在k i r c h h o f f 弹性杆理论基础上建立了d n a 弹性杆模型。b e n h 卸m ( 1 9 8 3 ) “”和l eb r e t ( 1 9 8 4 ) “” 讨论了封闭曲杆连接数的不变性。 k i r c h h o f f 方程并不是弹性杆平衡的唯一表达方式，s h i “”和h e a r s t 啪1 等 ( 1 9 9 4 。1 9 9 5 ) 以弹性杆的曲率和挠率为变量，建立了s c h r o d i n g e r 方程形式的另 一种数学模型。w e s t c o t t “1 ，t o b i a s 蚴，s w i g o n 嘲。v a nd e rh e i j d e n ( 2 0 0 3 ) 珏玎， s c h u b e r t ( 2 0 0 3 ) 嘲等研究了杆间自相接触情形。k e h r b a u m m l ，k a ih u 啪1 等( 1 9 9 7 ， 1 9 9 9 ) 给出用欧拉参数代替欧拉角的数学表达。 一般情况下，d n a 分子被模拟成均质的、各向同性的、线性圆截面弹性细杆， 关于d n a 弹性杆的研究已经有大量的文献发表在生物数学、物理化学、生物化学和 高分子化学以及数学物理等学科的刊物上，考虑杆间自接触的文献有：d n a 超螺旋 的有限元分析；d n a 超螺旋自接触模型；国内最新讨论了受曲面约束弹性细杆的平 衡问题，从动力学的观点讨论了超细长弹性杆的平衡稳定性。 2 0 0 0 年以来，国内学者在这一研究领域做了许多工作。刘延柱教授等嘲咖 对d n a 弹性杆模型的分析力学性质和稳定性作了深刻系统的研究，他们在弹性杆受 到几何约束情况下建立弹性杆的平衡方程，研究了弹性杆的动力学问题，推导出了 曲率和挠率的动力学方程及其哈密尔顿系统的静力学方程模型。 、尽管有许多新的成果出现，超细长弹性杆的研究及其在生物学和遗传学等领域 的应用仍然有许多问题需要我们进一步探索。例如，( 1 ) 对于超细长弹性杆的数值 计算问题由于模型的超细长和大变形的特点，需要对长距离的非线性边值问题进 行计算，如何使得长距离仿真保持系统特有的缠绕结构，不出现一团乱麻的失真现 象，对算法有特殊的要求。( 2 ) 超细长弹性杆所描述的生物模型如d n a 等，是带有 生物电荷的。它们所受重力通常可以忽略，但许多生物学学者指出，d n a 分子带有 静电斥力，这一体积力是不可忽略的。由于超细长的d n a 分子像毛线团一样缠绕在 很小的细胞核中，同种电荷的相斥作用必然影响这一细长分子的形状和结构。而且， 带有体积力的超细长弹性杆模型是非线性微分积分方程组，其积分限是未知的，这 对于数学处理和数值计算都带来困难。 主要研究的内容： 本文针对前面的两个问题展开研究，探讨受到静电斥力影响的超细长弹性杆的 数学问题，包括数学模型和数值方法的建立和分析。所研究的内容包括： i 、建立带有静电斥力的d n a 弹性杆的静力学平衡方程模型。由于连续模型是存 在奇点的，不利于数值计算。但由于这些奇点是可去奇点，我们通过无穷小代换和 2 引言 分段处理的方法，消除了模型中豹奇点，为模型的实际应用创造了条件 z 、设计了相应的微分积分方程的数值计算方法。对于初值问题，我们采用递推 的方法简化计算，得到了效率较高的算法。对于长距离的边值问题，我们采用有限 元方法进行数值离散。且由于积分曲线是未知量，因此，各单元是无法预知的为 了解决这一问题，我们设计了迭代方法，利用这一方法计算，得到了收敛的结果。 3 、在研究过程中我们还发现，利用k i r c h h o f f 的动力学比拟方法同样可以解决 超细长弹性杆的图形绘制问题。一条反复缠绕的细杆曲面的图形描述是一个复杂的 问题。但我们利用k i r c h h o f f 的动力学比拟方法，可以建立细杆表面的图形微分方 程，并设计图形描绘的跟踪算法，使得图形的描绘清晰且简单。 3 青岛大学硕士学位论文 第一章超细长弹性杆的平衡方程 1 i 预备知识 考虑空间中长度为l 的光滑曲线c ，在其无伸缩变形的条件下，建立沿曲线c 的弧坐标s ，点昂为曲线原点，曲线上任意一点p 在曲线上的位置由弧坐标j 确定 以空间中的固定点d 为原点，建立定参考坐标系( o 一勿f ) ，尸点在空间中的位置 由相对固定点o 的矢径，确定，为弧坐标s 的单值连续可微函数。矢量r c s ) 完全确 定曲线c 的几何形状 r ( s ) = 石d r i - c i ) 价目 l - ( 2 ) i v 0 ) ： 坚( 删)i-(3)asj r r 为沿曲线c 在户点处切线方向的单位矢量，称为切线矢量。r ( s ) 为曲线c 在p 点 处的曲率。n ( s ) 为沿d r 凼方向的单位矢量称为曲线c 在p 点处的法线矢量。 、，d 图1 1 切线的无限小角位移 从上图i i 可以看出，曲率r 等于p 点移动到相邻的p 点时，切线矢量的无限小角 位移4 妒与弧坐标增量厶之比，当心专0 时的极限，可表示为对弧坐标s 的导数 扣阱。m i 。ma 厶五- - * o 型d s l 删 i 加l “厶 一 因此曲率r 为度量曲线弯曲程度的参数，其倒数称为曲线c 在p 点处的曲率半径， 记作p ( s ) 4 第一章超细长弹性杆的平衡方程 烈j ) ：三 g - 定义以下矢量为曲线c 在尸点处的副法线矢量，记作口( s ) s c s ) = r ( s ) o ) 图1 2f r e n e t 坐标系 l - ( 5 ) l - ( 6 ) 如图1 2 由依附于曲线的矢量，丑，r 组成以p 为原点的右手坐标系( p n b t ) ， 称为曲线的f r e n e t 坐标系。| i v ，b ，r 所在的轴分别称为法线轴、副法线轴和切线轴，。 其基矢量厶，e b ，e t 分别与矢量，曰，r 相等。切线轴和法线轴张成的平面( r ，) 称为 曲线c 在p 点处的密切平面，法线轴和副法线轴张成的平面( ，b ) 称为曲线c 在p 点处的法平面，副法线轴与切线轴张成的平面( 丑，n 称为曲线c 在p 点处的从切平 面。设想p 点沿曲线c 以单位速度朝弧坐标j 的正向运动。则法平面绕副法线轴转动 的角速度为r ( s ) 。由于副法线矢量曰与密切平面正交，密切平面绕切线轴转动角速 度的模为 f ( s ) ：目 l - ( 7 ) 1 0 5 l r c s ) 称为曲线c 在p 点处的挠率。挠率为零时密切平面在空间中无转动，因此平面 曲线上任意点处的挠率均等于零。上述法线和副法线的定义均要求曲率r 不得为零。 在以下讨论中，除特别指出的情形以外，均认为g - o 条件得到满足。 引入矢量嗥( d ，定义为 口o ) = ( j ) 四+ r ( s ) t1 ( 8 ) 略o ) 称为曲线的d a r b o u x 矢量，其物理意义可理解为当p 点沿曲线e 以单位速度 朝弧坐标s 的正向运动时，f r e n e t 坐标系相对( d 一勃f ) 的转动角速度。m f ( s ) 与运动 学中的角速度概念相似，区别在于时间自变量r 被弧坐标j 替代利用式1 彤) l ( 6 ) 和l - ( 8 ) 中与曲率有关的第一项也可表示为 ， 一个 g b = r 1 - ( 9 ) 5 青岛大学硕士学位论文 矢量，曰，r 随弧坐标j 的变化规律可根据f o f ( s ) 由以下微分方程确定 _ d n 一= ，n 竺= 脚b 1 - ( 1 0 ) 出 。 i d t = t 将式l - ( 8 ) 代入计算后，化作 d n ；f 口一灯 出 i d b = 一f 1 - ( 1 1 ) d t ：a n d j 此矢量方程组称为f r e n e t s e r r e t 方程嗍。给定曲率r ( s ) 和挠率f ( s ) 为弧坐标j 的已 知函数，可从f r e n e t s e r r e t 方程解出矢量，丑，r 的变化规律。则空间曲线c 可根 据式1 - ( 1 ) ，由以下积分得出： ，o ) = r r p ) d 盯+ ，( 0 ) 1 - ( 1 2 ) 其中矢径，( o ) 表示昂点的空间位置曲率r ( s ) 和挠率f o ) 是确定空间曲线c 的两个 独立变量。借用分析力学名词，将确定曲线或曲杆几何形态的独立变量数目称为曲 线的自由度，则空间曲线的自由度为2 。 1 2 超长弹- 眭杆的平衡方程 1 2 1 基本假定： 讨论长度为l 的超细长弹性杆的平衡问题，曲杆的中心线为杆截面的几何中心 连成的不自相交的空间曲线c 。作类似于文献【2 8 】的假定： ( 1 ) 杆的中心线是2 阶以上的光滑曲线； ( 2 ) 横截面为刚性平面； ? ( 3 ) 忽略弯曲引起的剪切变形，横截面与中心线正交； ( 4 ) 忽略中心线的拉伸变形，任意二截面沿中心线的距离不变； ( 5 ) 相邻截面可绕中心线作相对扭转，扭角为j 的连续函数； ( 6 ) 杆为等截面，且截面内二主轴方向有相同的几何尺度； ( 7 ) 杆为均匀各向同性，弹性常数为常值，应力和应变满足线性本构关系。 6 第一章超细长弹性杆的平衡方程 图1 3 固定标榘( p 一刁z ) 与f r e n e t 坐标系( ，n b t ) 建立截面的主轴坐标系( p 一磅譬) ，其各坐标轴的基矢量分别为q ，e 2 ，岛，其中z 轴 与曲线c 的切线轴t 重合，即岛= = t 。设z 为x 轴与轴，y 轴与曰轴的夹角， 即截面相对f r e n e t 坐标系扭转的角度( 图1 3 ) ，则( p - n b t ) 与( p 一移曙) 之间的关 零= 0 剩 s ， ls 峨c o s x0 l 【0 0 1 j 曲杆的中心线确定以后，只需再确定扭角z ( s ) 随弧坐标的变化规律，曲杆在空 对( 0 一手，嘭) 转动m 为截面相对定坐标系p 一善叩f ) 的绝对角速度，即截面相对 f 姗e t 坐标系( p 一砌r ) 的相对角速度( 哆乞) 岛与( p 一脑r ) 相对定坐标系( d 一善，7 f ) 的牵连角速度，之和。利用式1 ( 8 ) 和表i - ( 1 3 ) 可导出 = 绵+ ( 誓) 岛哪+ 呸吃+ q 巳 ，删 青岛大学硕士学位论文 q = | r s i n z ，d 9 2 - - - - 9 0 0 s z ，鸭= f + d 出z 1 ( 1 5 ) 上式中沿切线轴的分量峨为挠率f 与相对扭角变化率a z d s 之和，即截面绕切线 轴的绝对转角对弧坐标j 的变化率，称为杆的扭率，以区别于前面定义的挠率f 。有 以下几种不同情形： 过扭转：d x d s 与f 同号，o j 3 f ， 纯扭转；d z d s 为零，鸭_ - - f ，欠扭转：毗出与f 异号，鸭 0 ，当p 一卅 占时，引入投影矩阵 q ( s ) = j 一订4 知h 声= 跏，即 胁乞褊俺棚c 怯精 将2 - ( 3 3 ) ，2 - ( 2 5 ) 两式代x 2 - ( 3 1 ) 式得到 州* k p 2 k n + k p 2 q ( s ) 一怎挚 = 鲁咖彬q 上潞 2 - ( 3 3 ) 2 - ( 3 铆 2 - ( 3 价 2 - ( 3 7 ) 将式2 - ( 3 6 ) 作为分布力代入2 - ( 3 ) 式，并注意到分布力偶在静力学问题中可以忽略， 即得到带静电斥力的弹性杆模型。 第三章d n a 弹性杆的数值方法及其仿真 第三章d n a 弹性杆的数值方法及其仿真 d n a 作为生物学的重要部分，自从其三维螺旋结构模型提出以来，关于d n a 弹性细杆模型的研究工作已经取得了不少成果其中大多是在k i 化h h o f f 的假设条件 下对其的解析研究。但对任意条件下的弹性杆，如带静电斥力，则必须用数值方法 来计算其几何形态。 3 1 初值问题的数值方法 考虑方程组2 - ( 2 8 ) 坐一r l 一竺1 2 一巴r + ! 一：五 万叫r 刊1 “j 如彳 彭生+ 2 坚f f 一里1 ：五 其中体积力分量由式2 - ( 3 6 ) 确定其中 工= 扣乜器叽 石形。上揣叽 方程组是非线性微分积分方程，而且积分路径是未知的。我们分两步来处理这一 问题。第一步，设积分路径是已知的。通常的方法是首先离散积分项，化为普通的 非线性常微分方程组，然后利用r u n g e - k u t t a 方法或a d a m a ；方法等求解，这些方法 已经有成熟的算法和完善的软件包但对于我们的问题，由于，的值的本身与整个 积分区问的未知量有关，而且积分套积分，利用r u n g e - k u t t a 方法求解的每一步要 计算大量的函数值，对于积分区间很长的情况，计算量很大。我们考虑利用递推的 方法简化计算。记 邝峥埘i 矿莘罴奇。 - - - - qp 2 枷+ 助2 q ( s ) g ( 力 3 - ( 1 ) 其中 1 9 青岛大学硕士学位论文 q d 2 h 上静 r ( s ，= ir ( a ) 姐 对于小的步长d s ，由于 g 。+ 妒p 。l 嚣等 矿羔：五g ( s + ( t s , o 删 + d s ) d g ( s + d s 一”3 也) 矿靠：l ，+ 出) r 一 而 g ( s + d s ，o + d s ) = r _ t ( x ) d x t , a t 、= f r ( 名) + r “a ) 舨一f “r ( 五) 以 一 烈s ，盯) + 孚仃( 盯) + t ( a + d j ) 一r ( s ) 一r o + d j ) ) 3 - ( 3 ) 3 - ( 3 ) 式的离散利用了二阶精度的梯形公式，如果需要高精度的算法，也可以采用更 高阶的积分公式如s i m p s o n 公式等离散。 利用递推公式3 - ( 3 ) ，我们把每一步计算g ( s ，力大量的积分计算简化成简单的 递推计算，从而对g 0 + 凼) 的计算，可以利用计算a ( s ) 的已有数据g ( s ，盯) 和递推过 程3 一( 3 ) 得到计算g ( s + d s ) 的数据g ( s + d s ，o + d o ) ，因而不需要计算大量的积分， 使得计算量大大减少 利用类似的思想，也可以建立g ( s + 凼) 和g ( s ) 之间的递推关系。但是由于积分 项分母的非线性，需要多次t a y l o r 展开。难以得到高阶、简明的算法。 下面考虑积分限未知的情况。在我们的模型中，积分限与曲线的形状有关例 如，考虑方程组2 - ( 2 8 ) ，2 - ( 3 6 ) 进行数值分析的主要困难是曲线积分 ，刮，o ) = ，( o ) + l t ( s ) d s3 - ( 4 ) f 的积分路径r 是待求变量。它依赖于电荷分布，而电场的作用又决定积分路径。这 种耦合关系给问题的计算带来很大困难。为了解决这一问题，我们考虑通过对电荷 密度p 和路径，分别作迭代来处理这一问题。注意到微分积分方程组2 - ( 2 8 ) ，2 - ( 3 6 ) 的解是路径r ：，= ，0 ) 的函数，我们可以将微分，积分方程组记为 t o ( r ) = 0 这一方程组可以利用各种方法化为等价的形式 ，= 圣( ，) 3 - ( 5 ) 如n e w t o n 迭代形式等。实际上，我们的计算过程为由方程组2 - ( 2 8 ) 解得曲 第三章d n a 弹性杆的数值方法及其仿真 率r 和挠率f ，然后代入f r e n e t - s e r r e t 方程组 盟：虚一灯 山 坚；一州 出 d t ：槲 凼 解得r ，最后利用3 - ( 4 ) 式解得，。注意到3 小) 式已经是3 - ( 5 ) 的形式了。 注：积分3 ( 4 ) 要计算一组值r ( s ，) 。每一个值都利用数值积分计算计算 量较大对于等距节点的情况，我们可以利用公式 r ( s “i ) = ，0 1 ) + r “r ( 盯) d 盯 计算。比如，利用梯形法可以得到 = + 等【r o ，) + r ( t + 1 ) 】 j = 0 , 1 ，一l + 对3 ( 4 ) 式作迭代 ，：“= 咖( ) ( f = o ，l ，)3 - ( 6 ) 对于迭代的初始近似t o ( s ) ，我们通过对电荷密度p 进行迭代求解，迭代过程为： ( 1 ) 取p o = 0 ，即初始假定系统不带电荷，对不同形状的弹性杆，问题可得到解 析解或数值解，= r ( s ，岛) 。 ( 2 ) 取0 = 岛 岛 成= p 。若已知f i r - - - - ，( 8 ，屏，f l _ 1 ) ，其中参数r ，- l 表示在 计算过程中，2 - ( 3 0 ) 式中的r 是曲线l 。的单位切向量。将其代入方程组2 - ( 3 ) ，取 p = 岛+ l 并利用1 - ( 1 1 ) ，1 - ( 1 2 ) ，1 - 0 5 ) 式得到，p ，卢h ，r ，) ， 注：计算过程可能出现发散，这是由于圣的条件不好一个常用的方法是寻找 预处理矩阵4 。使得迭代 r t “= ( i - 4 ) 噍+ 4 i 圣( ，i ) 有好的收敛性。 3 2 边值问题的数值方法 由于我们讨论的超细长弹性杆模型是结构问题，其模型常常是边值问题边值 问题的数值计算难度远远大于初值问题，通常的方法有s h o o t i n g 方法和直接离散的 方法这里的问题是：由于计算路径是未知量，给直接离散带来困难。我们采用直 接离散的方法，并且对积分路径进行迭代 本节考虑更一般的模型2 - ( 1 4 ) 的数值离散问题。 2 1 青岛大学硕士学位论文 3 ，2 1kir c h h o f f 方程的欧拉参数形式 模型2 一( 1 4 ) 是利用欧拉角为变量建立的这种模型的一个突出问题是不利于数 值计算。由于欧拉角表示的方程在0 = 七石时出现奇点，造成计算过程的失真，因此 通常采用欧拉四元数把模型变换为适合计算的形式 根据刚体有限转动的欧拉定理雠耵定义欧拉四元数g ，导数为q g = 吼 吼 吼 g p ls i n ( ) p 2 洲约 见咖( 约 c o s ( 约 q = 吼 9 2 吼 霉 其中p=hp ：p s r 鼋满足：9 7 9 = l 截面主轴坐标系( p - x y z ) 相对( o 一孝叩f ) 的方向余弦矩阵： 4 = r l t ，4 = l r i 吼嘞吼咱ijg 吼- 9 2 嘞i r = lg ，9 4 - g i - q 2l ，工= i - q 3g g l - q 2i l - q 2g lq 4 嘞j l 吼- g l 吼吧j + 爵+ 云一旌一爵2 ( q l q 2 一吼吼) 2 ( q l q 3 + 9 2 9 ) 1 则4 = l2 ( q t q 2 + q 3 q 4 ) 霹+ 订一彳一面2 ( q 2 q ，一q l q ) i l2 ( q l q 3 - q 2 q 4 ) 2 ( q 2 q 3 + q i q 4 ) 爵+ 订一彳一爵i 即； 匮 = 2 毒三豸戛 对i - ( 1 2 ) 求导写成微分形式： ，= 毛= 置( g ) f 幻) 岛 窖l 9 2 吼 g 3 - ( 8 ) 3 - ( 9 ) 3 - ( 1 3 - ( 1 1 ) 3 - ( 1 2 ) 3 - ( 1 3 ) 3 - 0 4 ) 第三章d n a 弹性杆的数值方法及其仿真 由l ，l - ( 3 2 ) 叫8 厮= 号，y = 三，= 孚m = 等得到欧拉四元数表示 的触行方程的标准形式为 一= q ，q = 圭r ( q ) 田+ 圭工t ( 口，= 置( g ) 工t ( g ) 鸭 3 _ ( 1 5 ) m = 、，一1 【a ( 一。) + 毛彳2 f o 岛= f 1 00 1 j = l0 。0l l o o y 。j 3 2 2 带静电斥力的四元数模型及其有限元离散 由于d n a 弹性杆带有负电荷，这样杆与杆之间就存在排斥作用，来保证杆中心 线之间维持一定距离以避免互相穿越，这种静电斥力可离散化为带电单元之问的作 用 设 黾= 窖i ，心+ = q ( 七= 1 2 ，3 ，4 ) 黾+ i = e = l 2 3 ) 3 - ( 1 6 ) 把3 - ( 1 3 ) 代入2 - ( 1 4 ) 则有； 啊2 彳+ + + 一1 = o 3 - ( 1 7 1 ) 如= 警一而- - - 3 - ( 1 7 2 ) 岛= 鲁一= o 双1 7 3 ) = 警一而_ o 3 - ( 1 7 4 ) 鬼= 鲁一毛= 0 3 - ( 1 7 5 ) k = 一而鲁+ 五鲁+ _ 鲁一毛警+ 2 ( c 彳) ( 噶毛一_ + 五而+ 而黾) 叫q + 而毛一屯而+ 玛黾一( 曰2 ) 一( 纠彳) ( 一弓+ 黾一而而+ 五矗) 卜毛毛一_ 氏+ 而而+ 而毛】一( 迸2 ) 一( 而。4 ) = 0 3 - ( 1 7 6 ) 青岛大学硕士学位论文 呜= 喝鲁一警+ 五鲁+ 屯警+ 2 ( 曰) ( 一粕毛+ 而一而”五 ) 屯毛+ 毛+ _ 弓一毛 一( 群2 ) 一( c 彳) ( 而黾+ 而+ 毛而一而) - x x 5 + x s x 6 一而而+ 而黾卜( 霹2 ) 一( 而4 ) = o 3 - ( 1 7 7 ) 嚏；啊鲁+ 毛鲁一乇鲁+ 而鲁+ 2 ( b c ) ( - x , x , + 毛一而而+ 而黾) 而毛一_ + 五而+ 屯黾一( 趟2 ) 一( 们) ( 而毛一x 4 x 6 + x ，x 7 + x 2 x t ) 卜屯+ 而一_ 而+ 而黾。( 砰2 ) = 0 = 鲁+ 2 ( 一而毛一+ 而而+ 而吒) 葺- 一 ( 一x , x 5 + x 3 x 6 一x 2 x 7 + x t x s ) x t l + “2 ) = o 忙誓+ 2 ( 一_ 墨+ 而毛一楠+ 毛黾) 与一 ( b 而+ + 毛而一而黾) 薯- + ( 五2 ) = o 忙誓+ 2 ( 吨而+ 而+ 而一毛黾) 3 - ( 1 7 8 1 3 - ( 1 7 9 ) 卜( 1 7 1 0 ) ( - x s x j 一+ x _ l x 7 + x 2 x s ) x 9 + ( f 3 2 ) | - o 3 - ( 1 z 1 1 ) 以上列出的1 1 个方程危构成封闭的方程组，分布力的规律给出以后，可以解出1 1 个未知变量毛。 应用有限元方法将弹性杆延中心线用n + 1 个节点只= o ，l ，彩划分为撑个单 元马叮- - - 1 ，2 帕，其中的第7 单元岛介于第，一1 节点c 。与第_ ，节点之间。 将第七变量在第- ，节点另处的值记作吒( 七= l ，2 ，l l ；j = o j ，叻设杆 的变形前总长度为厶，第了单元d i 的变形前长度为工”若单元划分均匀，则 厶o = 厶加设单元d i 在受力状态下的轴向应变为s 7 ，则变形后的长度为 工，= 工_ o ( 1 + 岛) 3 - ( 1 8 ) 当单元分隔足够精细，即砧足够小时，可认为截面形状和物理参数在单元内为常值， 且可足够精确地令变量黾在单元d _ 内按线性规律插值。引入介于零与i 之间的无量 纲变量孝，将d 7 单元内的变量吒记作屯j ，令 x k a ( 善) = 【l 一善k j - l + 参t ( t = 1 ，2 ，l l ；j = 1 , 2 ，n ；o f s l ) 3 一( 1 9 ) 将上式在d i 单元内积分，得到的平均值等于两端节点值的平均，记作墨7 第三章d n a 弹性杆的数值方法及其仿真 耳j = k j g ) d f = 兰学o b l ，2 ，l l ；_ ，= l ，2 ，功 3 - ( 2 0 ) 以撤号表示变量对s 的导数， 在b 单元内的导数。7 也可利用两端节点值计算 户掣：手掣：毕。乩如，l l ；川加埘) 3 一( 2 1 ) 。d ， 三：d 善工： 、。一 ” 上式中的岛为变形后的单元长度，如式3 - ( 1 8 ) 所示。其中的轴向应变岛可用轴向力 e 在该单元内的平均值，即，7 的平均值表示为 。 勺：誓；毕 3 - ( 2 2 ) x ：2 置j 一 其中骛为岛单元的抗拉刚度 第f 单元受到的静电引力可以利用积分公式2 书6 ) 给出 e 5 j k 户2 村，+ 幼2 q c ，p ! 。再：i ! j j ；：；警f 拈 。q s ， 其中以为第i 个单元的中点。 最可以分解为 巧= 毛 。 3 - ( 2 4 ) 这里磊是第i 单元和第j 单元之间的静电斥力其离散解为 f i - - 6 弓- p 2 纽 l j 当s l sj i z 毛彬q 揣厶句当l - 一除z 硒) 注：静电斥力是以平方阶随距离增加而减小的。许多文献采用截断的方法 减少3 - ( 2 4 ) 式的项数来简化计算。但对于d n a 等超细长弹性杆，这样傲比较 困难。因为其毛线团似的缠绕，使得弧长参数的距离和实际单元的距离没有可 循的关系。在实际计算中这样做计算量减少不显著。但可以利用对称性简化计 算。 将( 3 2 0 ) ，( 3 2 1 ) 和( 3 2 5 ) 代入方程组( 3 - 1 7 ) 即得到离散方程组。 由于方法的积分区域是未知量，我们的单元实际上是未知的，我们的有限 元削分有一定困难解决这一困难也可以利用前面初值问题的迭代技巧 青岛大学硕士学位论文 第四章弹性杆的曲面方程及应用 超细长弹性杆模型如d n a 的结构模型，往往是象毛线团一样的缠绕结构， 其图形的描绘比较困难。我们发现，利用k i r c h h o f f 的比拟方法，可以建立弹 性杆的曲面方程，从而简化图形的表述问题。 类似于k i r c h h o f f 的比拟方法，我们把固定截面弹性杆曲面看作一个垂直 于弹性杆轴线的平面曲线沿轴线运动而得到。如图4 1 ： r 圈4 1 设 ， f x = x o ) ro=y=y(s)4-(1) 【：；：( s ) 是形成弹性秆表面的初始曲线，o = r ( o ) 。我们有 j ( s ) = 4 ( s l 一，( o ) ) + t o ) 4 - ( 2 ) 其中一p ) 是旋转矩阵，它描述，0 运动到r ( s ) 的旋转量。对方程组4 一( 2 ) 两端求导得 到 i d r = 警( ，o 一，( o ) ) + “j ) ，0 = j ( 0 ) 4 - ( 3 ) 注意到a ( s 1 是正交矩阵，从4 - ( 2 ) 中解得 ，o 一，( o ) = a o ) 7 ( r 0 ) 一t o ) ) 4 - ( 4 ) 代入4 一( 3 ) 得 d ：r = ：掣4 0 ) 7 ( j ( j ) 一，( 劝+ r ( j ) ，f o ；j ( o ) 4 - ( 5 ) a 5d s 又由于 第四章弹性杆的曲面方程及应用 a r = j 知 i d a 部) r + _ 鲁= 。 从而4 - ( 5 ) 还可表示成 芸啡) 鲁( 邢) - ，阱t o ) ，r o = f ( 0 ) 4 - ( 6 ) 4 - ( 7 ) 通过求解方程组4 一( 5 ) 或4 一( 7 ) 可以得到弹性杆曲面的网线图。更重要的是我 们得到了曲面数据，可以在弹性秆曲面接触问题中得到应用。 青岛大学硕士学位论文 第五章例子及分析 窘叫r 一罢) 2 三矿= 五 r 虿d 2 t + 2 警( f 一争= 五 其中q = 尝一了m 2 ，五，五为静电斥力在法线和副法线方向的分力，可写成 z = 形，a = f 形 5 - 0 ) “一c 争嘶2 上器砂席 工= 竽= 二掣竖竺l - s 钧 五：莽k2 二：掣2 墨= ，删 式子2 - ( 2 8 ) ，5 - ( 2 ) ，5 ( 3 ) 就是d n a 分子带有静电斥力的平衡方程的表达式。 为了便于分析，我们设圆截面弹性杆起始端有一较强的电荷碍，而杆的其他位 置有均匀弱电荷。因而只考虑9 与其他电荷的作用，而忽略其他电荷之问的作用。 t ，一嘲缈搬署盯 s 删 第五章例子及分析 ( a ) 种o o 0 0 2 卿0 0 0 0 3 青岛大学硕士学位论文 ( c ) 卿0 0 0 0 4 ( d ) q - - o 0 0 0 0 5 圈5 1 ：弹性杆受静电斥力，影响情况 第五章例子及分析 图( 5 1 ) 是弹性杆受静电斥力影响的情况( 实线代表没有力的情况，