（应用数学专业论文）弹性杆模型方程的数值谱方法.pdf

摘要 弹性杆是科学研究和工程开发应用的重要物理模型，近年来在生物学研究中也 发挥了重要作用。而数值仿真是弹性杆研究的重要方法之一。本文基于谱离散方法 给出了运动弹性杆的高精度数值算法，主要工作包括 ( 1 ) 弹性杆结构动力学模型的谱离散方法。算法对基于k i r c h h o f f 假设导出的描 述弹性杆结构的非线性s c h r o d i n g e r 方程采用复三角谱离散得到高精度的数值计算方 法。与实g a l e r k i n 方法比较，这样的离散在相同精度的前提下，大大节省了计算量。 ( 2 )基于k i r c h h o f f 假设的弹性杆动力学方程是非线性偏微分代数方程组。文献 中的相应的谱精度算法一般是用来求解边值问题的，通常只用于对空间变量，时间 变量的离散采用r u n g e k u t t a 方法等传统的方法离散，这样离散一般达不到谱精度。 本文将空问变量的谱离散方法和时问变量的谱延迟修正技术相结合，给出了新的数 值计算方法。这一算法不但精度达到了谱精度，而且是a 稳定和保辛的，对于长时 间的数值模拟提供了可靠算法工具。 由于不同初边值条件和弹性杆数值仿真的不同需要，本文分别给出了基于三角 谱和l a g e n d r e 正交多项式谱的数值离散格式。另外，由于模型非线性导致的数值方 法的复杂性，利用通常的解法计算量很大。本文采用预条件加速等技巧简化了相应 的计算。 关键词：弹性杆s c h r 6 d i n g e r 方程 动力学方程谱g a l e r k i n 方法谱延迟校 正 a b s t r a c t e l a s t i cr o d ，a ni m p o r t a n tp h y s i c a lm o d e lf o rs c i e n t i f i cp r o b l e m sa sw e l la se n g i n e e r i n g p r o j e c t s ，a t t r a c t e dg r e a ta t t e n t i o ni nb i o l o g i c a lr e s e a r c hi n r e c e n ty e a r s o n eo ft h em o s t i m p o r t a n ta p p r o a c h e si ns t u d y i n ge l a s t i cr o d si s n u m e r i c a lm e t h o d t oi m p r o v et h e p r e c i s i o no fn u m e r i c a ls i m u l a t i o n ，h i g hp r e c i s i o nm e t h o d so nd y n a m i cm o d e l so fa l l e l a s t i cr o da r es t u d i e di nt h i sp a p e r t h em a i nr e s u l t si n c l u d e i 、an u m e r i c a lm e t h o df o rs o l v i n gs t r u c t u r a lm o d e lo fa l le l a s t i cr o di sg i v e n b yu s i n g f o u r i e rs p e c t r a lb a s i si ng a l e r k i nm e t h o dt os o l v en o n l i n e a rs c h r o d i n g e re q u a t i o nb a s e d o nk i r c h h o f fh y p o t h e s i s ，n u m e r i c a ls o l u t i o no fs p e c t r a lp r e c i s i o nc a nb eo b t a i n e dw i t h r e l a t i v e l yl e s sc o m p u t a t i o n i i ) t h ed y n a m i c a le q u a t i o n so fa ne l a s t i cr o di s an o n l i n e a rp a r t i a ld i f f e r e n t i a l a l g e b r a i c e q u a t i o ns y s t e m s p e c t r a lm e t h o d sf o rs o l v i n gt h ee q u a t i o n so ft h i sk i n du s u a l l yu s e s p e c t r a l m e t h o d st od i s c r e t es p a t i a lv a r i a b l e sb u tr u n g e k u t t am e t h o dt o d i s c r e t e t e m p o r a lv a r i a b l e ，t h e r e f o r et h ep r e c i s i o nu s u a l l y c a nn o tr e a c hs p e c t r a lp r e c i s i o n an e w m e t h o dw h i c hc o m b i n e st h es p e c t r a lm e t h o do ns p a t i a lv a r i a b l e sa n ds p e c t r a ld e f e r r e d c o r r e c t i o nm e t h o do nt e m p o r a lv a r i a b l e si sg i v e ni nt h i sp a p e r t h em e t h o di so fs p e c t r a l p r e c i s i o n ，a s t a b l ea n ds y m p l e c t i c ，w h i c hm e a n sm o r er e l i a b l ef o rl o n gt i m en u m e r i c a l s i m u l a t i o n s p e c t r a lm e t h o d sb a s e do nf o u r i e rb a s i sa sw e l la st h a tb a s e do nl a g e n d r e o r t h o g o n a lp o l y n o m i a l sb a s i s ，a r ed i s c u s s e di nt h i sp a p e rf o rd i f f e r e n ti n i t i a lv a l u ea n d b o u n d a r yv a l u ec o n d i t i o n sa n dd i f f e r e n t n e e do fn u m e r i c a ls i m u l a t i o n i na d d i t i o n ， p r e - c o n d i t i o n i n ga c c e l e r a t i o ns k i l li su s e dt oo v e r c o m et h ec o m p u t a t i o n a lc o m p l e x k e yw o r d s ：e l a s t i cr o ds c h r 6 d i n g e re q u a t i o n s d y n a m i c se q u a t i o n ss p e c t r a l g a l e r k i nm e t h o ds p e c t r a ld e f e r r e dc o r r e c t i o nm e t h o d 青岛人学硕士学位论文 学位论文独创性声明 本人声明，所呈交的学位论文系本人在导师指导下独立完成的研究成果。文中 依法引用他人的成果，均己做出明确标注或得到许可。论文内容未包含法律意义上 已属于他人的任何形式的研究成果，也不包含本人已用于其他学位申请的论文或成 果。 本人如违反上述声明，愿意承担由此引发的一切责任和后果。 论文作者签名： 扣戡 醐越少狮加 学位论文知识产权权属声明 本人在导师指导下所完成的学位论文及相关的职务作品，知识产权归属学校。 学校享有以任何方式发表、复制、公开阅览、借阅以及申请专利等权利。本人离校 后发表或使用学位论文或与该论文直接相关的学术论文或成果时，署名单位仍然为 青岛大学。 本学位论文属于： 保密口，在年解密后适用于本声明。 不保密吃 ( 请在以上方框内打“”) 敝储戳。节三式 导师猛名：幻汾 日期z 勺户年聊日 1 日期叨6 r 卯 3 3 青岛大学硕士学位论文 引言 弹性杆的研究有广泛的工程背景，许多工程构件如钻杆、纤维和地下电缆等都 可以模型化为弹性杆研究。近年来，随着生物工程技术的发展，人们通过微观实验 发现许多大分子如大肠杆菌和d n a ，都呈现弹性杆的结构特征。因此，近3 0 年来， 弹性杆模型和理论在分子生物学，生物化学，物理化学等学科的科学研究中得到用 用和发展，成为新的研究热点。 弹性杆理论也是一个古老的理论。从1 8 3 0 年d b e r n o p u l l i 和e u l e r 的开创性工 作算起，其研究已经有1 8 0 年的历史。2 0 世纪5 0 年代以来，d n a 和分子生物学的 研究快速发展。人们从微观实验中发现，大分子的结构物理性质和弹性杆类似，如 一个小的片段表现出很强的刚性而较长的片段则表现出柔性。因此在近3 0 年中，弹 性杆成为研究d n a 和大分子的重要动力学模型。 弹性杆的理论包括静力学和动力学两部分。18 5 9 年k i r c h h o f f 的弹性杆动力学 比拟理论奠定了弹性杆静力学的理论基础。这一理论在c l e b s c h1 8 6 2 年的著作中得 到系统的论述n 。根据k i r c h h o f f 理论，在k i r c h h o f f 假设下建立的结构动力学模型， 和刚体定点转动的动力学方程一致，只是弧坐标s 变换成时间变量t 。这一比拟技 巧可以把对动力学问题研究的成熟结果应用于弹性杆的结构力学分析，成为1 0 0 多 年来弹性杆结构分析的重要工具。相应的结果在n i z z e t e 和g o r i e l y ( 1 9 9 9 ) 1 的论文中 得到总结。另外，基于对d n a 弹性杆模型的研究需要，近3 0 年来，人们发展出许多 新的模型和数学工具描述d n a 的结构特性，如w h i t e ，f u l l e r 嘲，c r i c k 吲和b e n h a m ” 等引入的研究超细长、大变形弹性杆的几何理论，s h i 和h e a r s t 建立的描述弹性杆 结构的s c h r o d i n g e r 方程，s t a r o t s i n 啪】，h e i j d e n 等给出的弹性杆自接触现象的几何 和力学描述等。k i r c h h o f f 弹性杆假设忽略拉伸正应变和弯曲剪应变的作用，对许多 问题的描述有较大的偏差。比如，实验表明，d n a 在运动中出现明显的拉伸和剪切 应变。作为对k i r c h h o f f 弹性杆模型的改进，c o s s e r a t 考虑杆的轴向线应变和弯曲剪 应变等冈素，建立了更精确的弹性杆平衡方程，自此之后，许多文献给出了c o s s e r a t 弹性杆的不同表述。 相比之下，弹性杆的动力学理论要复杂的多，近年米也更加受人关注。t a b o r ， k l a p p e “1 9 9 6 ) n 妇导出了杆截面的运动学关系式，以及用拓扑学参数描述杆几何形态 的运动学方程。g o r i e l y , s h i p m a n ( 2 0 0 0 ) u 厶1 3 1 在线性化动力学方程的基础上，用数值 方法研究了螺旋线平衡的一次近似稳定性。t s u r u ( 1 9 8 6 ，1 9 8 7 ) n t 巧 ， t o b i a s ，c o l e m a n ( 1 9 9 6 ) u 训等研究了杆的扭转转动，a n t m a n ( 1 9 7 9 ) n 7 d 羽，c o l e m a n ( 1 9 9 2 ) “训 研究了弹性波问题。o o l d s t e i n 幽1 ，p o w e r s 雎，w o l g e m u t h 等( 1 9 9 8 ，2 0 0 0 ) 瞄1 则研究了弹 引言 性杆在黏性介质中的运动。 弹性杆的数值仿真是近年来发展起来的重要研究方法。经典的弹性杆模型通常 使用欧拉角描述。虽然欧拉角在描述弹性杆复杂几何特性时具有优点，但由于欧拉 角存在奇点这样的缺陷，影响了数值计算的精度。同时方程的强非线性也给数值计 算带来困难。近年来，人们引入e u l e r 参数和一般四元数的方法来解决这一问题。 由于弹性杆模型是在空间曲线上建立的局部坐标描述的偏微分方程组或偏微分代数 方程组，因此，许多作者研究了相应的数值计算，如基于c o s s e r a t 瞄弹性杆的有限 元方法，保结构算法和s h i ，h e a r s t 泓捌等最近提出的弹性杆的数值精确算法等。 谱方法是用于微分方程字间离散的一种高精度数值方法。谱方法首先是和三角 函数和正交多项式联系在一起，包括f o u r i e r 级数，l e g e n d r e 多项式，c h e b y s h e v 多 项式，j a c o b i 多项式等。这些函数系本身具有正交性，因此在用作有限元方法的基 函数和检验函数时，得到的刚度矩阵是小带宽甚至对角型矩阵。同时它们又有很高 的逼近精度，在微分方程的积分中有最高的代数精度。根据检验函数分类的不同形 成谱g l e r k i n 方法，配置方法，t a u 方法等。t a u 方法是对于g a l e r k i n 方法在周期边 界条件问题上的改进，也被称为p e t r o v g a l e r k i n 方法。 2 0 世纪8 0 年代中期之后发展起来的较新的谱方法是将g a l e r k i n 方法和g a u s s i a n 四边形公式结合在一起，先运用满足边界条件的基函数近似方程的解，运用g a l e r k i n 离散变分得到离散的方程组，计算出方程解在四边形节点的值，得到方程的数值解， 这种方法会涉及到数值积分，称为g - n i 方法。g o t t l i e b 和o r s z a g 嘶1 出版的专著首先 将谱方法的一些理论与方法统一起来作为一个整体数学研究对象。谱方法被应用于 更多的问题，包括变系数和非线性方程等等。1 9 8 8 年以后，有关谱方法研究的文献 就越来越多。b o y d ( 1 9 8 9 ) 脚。详细阐述了谱方法，给出许多新的结果，特别是对 于无边界圆柱形和球形坐标系统的相关问题给出了较好的研究结果。f u n a r o ( 1 9 9 2 ) 瞄1 和g u o ( 1 9 9 8 ) 恻讨论了用多项式展开近似微分方程。t a d m o r ( 1 9 9 8 ) ，g p t t l i e d ，和 h e s t h a c e n 恻用谱方法计算了一阶双曲闻题，2 0 0 2 年c o h e n 用来求解波动方程等， 这些是谱方法最初被用来来解决一些特定的问题。8 0 年代后期之后关于经典谱方法 的研究主要是围绕复区域上的高阶问题。谨1 。 近年来，关于谱方法的研究更为广阔，研究角度也呈现出多样性。2 0 0 6 年， g r e e n g a r d o 列从对于微分方程误差校正的角度，提出的谱延迟算法给我们提出了新的 时问变量的离散思路。这种方法通过对时间变量把谱延迟技巧，用于求解常微分方 程可以达到谱精度。h u a n g 。h 龉蚓柏：在其基础上考虑将谱延迟校正和n e w t o n k r y l o v 子空问方法结合使用，并将其推广到微分代数方程 ：，用于求解常微分方程组和微 分代数方程组，可以达到谱精度，并使得方法更加快速。 2 青岛大学硕士学位论文 有关运用各种数值方法解决弹性杆模型方程的研究一直在继续，近几年有关成 果包括2 0 0 5 年刘延柱等人建立用欧拉参数表示的k i m h h o f f 方程，导出杆的总势能 变分。对于静电引力和几何约束两种不同计算方案，列出离散化的弹性杆平衡方程 2 0 0 6 年m i n g h a oc a id 刀等人运用时间可空间区域的守恒元和解元方法( c e s e ) 计算 了细长弹性杆的线性和非线性波动方程模型，2 0 0 7 年张光辉汹1 运用k i m h h o f f 弹性 杆模型的动力学比拟技巧，建立了由四元数的描述超细长弹性杆曲面的常微分，积 分方程组，利用a d a m e s 方法和递推方法设计了方程的数值解法。进两年研究较多 的是关于c o s s e r a t 杆数值方法，运用g a l e r k i n 方法和c o s s e r a t 点方法来研究拉伸变 形的弹性杆方程删，由此可以看出越来越多的数值方法别用于各类弹性杆模型方 程的求解。 本文将谱方法用于求解弹性杆模型方程，包括弹性杆平衡方程和弹性杆动力学 方程，运用谱方法的高精度性，能够得到较好的数值结果，从而能够描述弹性杆在 满足一定平衡条件形态以及动力学中受力条件下的运动状态。本文共分为四章，分 别包括： 第一章介绍谱方法的基本知识，包括谱基函数和微分代数方程的加速谱延迟校 正方法。谱基函数我们主要介绍了本文用到的f o u r i e r 基函数和特征l a g r a n g e 多项 式基函数，包括形式，性质等等。另外简要叙述了h u a n g 的求解微分代数方程的一 般步骤，以及方法中包含的一些其他方法的介绍。 第二章是弹性杆s c h r s d i n g e r 方程的谱方法，本章首先简要推导了弹性杆平衡方 程。针对圆截面螺旋杆的情形，方程满足周期为2 n 的边界条件，构成了常微分方程 组的两点边值问题，运用f o u d e r 函数系作为基函数和检验函数的g a l e r k i n 方法将方 程离散得到变分方程，为代数方程组。另外对于给定初值的s c h r s d i n g e r 方程，构成 常微分方程的初值问题，运用加速谱延迟校正方法求解。对于两类倩形，运用不同 方法，分别计算得到数值结果，描绘出了弹性杆的平衡状态。 第三章是弹性杆动力学方程的谱方法，本章首先导出基于刚性截面假设的 k i r c h h o f f 弹性杆动力学方程组的形式，为偏微分方程组。首先将方程在空问一卜离散， 根据基函数的不同特点，选择f o u r i e r 函数系和特征l a g r a n g e 特征多项式作为基函 数对方程离散，变分得到半离散的微分代数方程组。然后运用加速谱延迟校正方法 进行时间离散，得到方程的数值解，具有谱精度，进一步根据弹性杆的模型方程计 算模拟出弹性杆不同时间的运动状态。 第四章最为本文的总结叙述了本文的问题研究的主要结论，以及需要改进的地 方。 第一章预备知识 1 1 弹性杆预备知识 第一章预备知识 由于弹性杆中心线的描述涉及到曲线几何学，首先简单介绍相关知识m 川1 。 以空间中的固定点o 为原点建立定参考坐标系( d 一勿f ) ，讨论光滑空间曲线 c 。选定曲线上的异点为原点，建立弧坐标s 则曲线上的任一点p 可以由矢径，确 定。考虑曲线p 点的切矢量 r ( j ) ：i d r ( 1 1 ) r ( s ) 对s 的导数的模定义为曲线c 在p 点处的曲率，记作 = 吲 沿d r d s 方向的单位矢量称为曲线c 在p 点处的法线矢量， 理论知 o ) = i 1i d t ，( 誓o ) 定义曲线c 在尸点处的法线矢量 记作n ( s ) ，由极限近似 ( 1 2 ) b ( s ) = t ( s ) x ( s ) ( 1 - 3 ) 则矢量，曰，z 组成以p 为原点，依附于曲线的右手坐标系( p - n z t ) ，称 为曲线的f r e n e t 坐标系各坐标轴分别称为法线轴、副法线轴和切线轴。 记 邢，= 吲 f ( s ) 称为i i | l 线f 在p 点处的挠率。 由矢量，b ，随弧坐标的变化规律得到 d n ：仍一g d 5 d b ：一z n 出 d t ：槲 出 此矢量方程组称为f r e n e t - s e r r e t 方程。给定曲率砥s ) 和挠率f ( s ) 为弧坐标s 的已 青岛大学硕士学位论文 知函数，可从f r e n e t - s e r r e t 方程解出矢量n ，b ，t 的变化规律。则空间曲线6 - 司 根据式( 1 1 ) ，由以下积分得出 ，( s ) 2j ：丁( d 仃+ ，( o ) 其中矢径，( 0 ) 表示昂点的空间位置。 在曲杆几何学中，以中心线上的任意点尸为原点建立截面的主轴坐标系 ( p 幻囝，各坐标轴的基矢量分别为a ，晓，岛其中z 轴与曲线c 的切线轴t 重合， 即岛= q = t 。设z 为j 轴与轴，轴与曰轴的夹角，即截面相对f r e n e t 坐标系 扭转的角度。由各坐标轴的基矢量与曲杆的曲率鬈、挠率f 和相对扭角z 的关系 堕：r 甜亟、1 乞一肫。泓 (14)ds kd s 2 “ 生d s 一( 甜亟d s ) ”您i 呲 5 ， k，1 “ 鲁叫c 。一s i n x e 2 ) ( 1 6 ) 可知曲杆的几何形态可以由茁、f 和z 这3 个独立变量决定。令彩为截面相对惯性 坐标系的绝对角速度，可以得到关系式， 国= 娇+ ( 警弘= 劬q + 哆p ：+ 皑p ， 其中 劬硝s i n 石0 ) 2 = k c o s 石皑可+ 警 因此也可将c o , ( s ) ( i = 1 , 2 ，3 ) 作为确定曲杆几何形态的3 个独立变量。 设k = x ( o ) 是初始截面点集，它沿着中心线运动到x ( s ) 处得到弹性杆的曲面 方程 警= 詈酣( 耶) - ，) + p 3 ( 1 7 ) 【x ( o ) = x o 其中，q ( s ) = ( e l , e 2 ，巳) 是x ( o ) - r ( o ) 到x ( s ) - r ( s ) 的旋转矩阵。 由该方程可以确定曲杆在空问中的几何形态。 1 2 谱离散方法 第一章预备知识 谱方法嘲是用于微分方程空间离散的一种经典数值方法。谱方法的关键步骤在 于试探函数和检验函数的选取。试探函数是一些适当的基函数的组成的线性组合， 用来表示方程的近似解。而检验函数是用来保证截断函数展开能够满足微分方程和 相应的边界条件。这种基于试探基函数的谱方法我们称为经典谱方法，建立在全局 无限可微的区域，并且试探函数为正交函数集。而检验函数选取的不同形成了不同 的方法，包括g a l e r k i n 方法，配置法，t a u 方法等。我们以三角函数集和l e g e n d r e 多 项式作为谱基函数集来研究。 1 2 1f o u r i e r 展开 选取【o ，2 万】上的正交函数集合为 么( x ) - - e 服，k 为整数， 则有 “( 堋灿= 2 硪= ：笔， 对于【o ，2 万】上的复值函数材，有 玩= 石1f石“。)p船，七：0 x ，+ l ，控， 2 芴工“矿觑，七2 ，+ l ，控， 如果z ，在【o ，2 万】上r i e m a n n 可积，则上式积分是存在的。这样对于任一函数甜， 其f o u r i e r 展开可以表示为 甜= u k 办_ 一 定义【o ，2 万】上的内积( ，v ) = r 7 ”( x 声( x ) 血 在实际计算过程中，我们将函数u 近似表示为上式的截断形式 甜= u k 九，n 为正整数 ( 1 8 ) 这样就要将正交函数集定义存以下具有2 n + 1 个自由度的函数集里 s ：+ 。= e i h , k = o ，1 ，2 ，n ) 即测试基函数的阶为2 + 1 ，( 1 8 ) 式称为u 的2 + 1 阶截断f o u r i e r 展开式。 在许多实际应用中，特别是对于数值方法，不能应用f o u r i e r 展开的标准形式， 6 青岛大学硕士学位论文 因为对于要求的函数，我们不知道它的f o u r i e r 展开的系数，这就需要一些近1 以离散 形式的计算 对于【o ，2 万】上的函数甜，我们考虑在下面一些点 x ，：堕，歹：0 ，。1 ，2 nx = 二，= ， 2 n + l 。 。 则函数材的2 + 1 阶截断f o u r i e r 展开式在这些点的系数为 i k - 赤善甜( t ) p 一，七一跏， 1 9 并由正交关系式， 上2 n + l 争k 0 委仁翟巩舻0 + 1 匀l其他 可知 甜( _ ) = 玩卢，j = 0 ，2 n i 一 通过插值公式可以得到函数甜的离散f o u r i e r 展开式 u n ( x ) = u k e 触 k n 将( 1 9 ) 式代入( 1 1 0 ) 式可得 2 _ ( x ) = z ，( _ 鼽( x ) j = o 其中( x ) = 丽l 荟n p 忆) ( 1 1 0 ) 函数 ( 工) 即是三角函数集最构成的线性空间中的函数，满足 ( 而) = 磊，歹= o ，j j 1 2 n 若，h 为由正交函数集最组合得到的函数，我们可以定义其在此字间上的 内积， u n , ) = 熹薹z ，( _ ) 矿( 一) t h 以匕可知，在求解微分方程时，只要求得截断f o u r i e r 展开式在离散点的值，就可 第一章预备知识 以得到方程的近似解，而这些离散点的值具体求解过程在下一章结合弹性杆模型方 程分析。 1 2 2 特征l a g r a n g e 多项式 首先介绍l e g e n d r e 多项式，在区间【- l ，l 】上，对于任意整数七，其形式为 。cx，=专笠1c一-，(；)(2七k一27ixk-2ti=o 厶( x ) = 专( 一1 ) 扩“ 厶 ， 其中【k 2 】表示七2 的整数部分。于是有厶( x ) = l ，厶( x ) = x 。 对于卜l ，l 】上任意函数甜可以用l e g e n d r e 多项式展开 玩= ( 七+ 批小) 厶出 与以上f o u r i e r 展开类似，j 毫y - i ，i 】上l e g e n d r e 多项式展开内积 ( 州) m = f 。甜( z p ( x ) 出 对于一般的数值近似，对上式采用g a u s s 数值积分。我们考虑离散形式，因此需要 选取一些合适的节点。l e g e n d r e 多项式系列的节点选择有多种，包括 g a u s s - l e g e n d r e ( g l ) 节点，g a u s s l e g e n d r e - r a d a u ( g l r )节点和 g a u s s - l e g e n d r e - l o b a t t o ( g l l ) 节点等，同时运用不同的节点在求数值积分的时候它们 对应的权函数也不同，我们介绍一下g l l 节点，其形式为 x o = - l ，h = 1 ，x ，( _ ，= 1 ，2 ，n - i ) 为l n 的零点， 它所对应的权函数形式为 驴n ( n 初。0 + 1 ) “( x ，) 1 v 一” 以l e g e n d r e 多项式和g l l 节点为基础引入特征l a g r a n g e 多项式 小) = 丽1 丽( 1 - x 2 ) 揣牛0 一， 满足，( ) = 靠，j ，k = o ，n 。 3 、- 、 x ，l 己 甜 。脚 = 掰 青岛大学硕士学位论文 对于【一l ，l 】上任意函数材可以用n 阶特征l a g r a n g e 多项式展开为 | 材( x ) 名甜( x ) = u k 丸( z ) k ；0 并且可以得到 其中 塑盟。q甜i(t)d ，钍工i x 、j ， 岛= 掣= 揣m 蝴厶( _ ) ( 一一_ ) “一山 一业掣 i f i = j = o ； 4 1 3 微分方程的加速谱延迟方法嘶瓯“1 i f i = j = ； o t h e r w i s e 在数值计算时，对于偏微分方程我们，般将空问变量离散得到常微分方程和微 分代数方程，因此微分方程的谱方法主要体现在求解微分方程和微分代数方程。近 年米g r e e n g a r d 1 和h u a n g 等瓶矧提出的对时间变量的谱积分方法和k r y l o v 迭代 方法用于求解微分方程和微分代数方程可以达到很高的精度。 考虑以下形式的微分代数方程组 ，( “，旷f ) = 口 ( 1 1 1 ) 其中h c ， f 是光滑可微的，初值m o = u ( o ) 。 令盯= u ，则由p i c a r d 积分得 m ：m 0 + f u ( r ) d f ( 1 1 2 ) 选取时间的离散网格0 = ，o f i 0 丁，在【o ，f i 】上取p 个g a u s s l e g e n d r e l o b a t t 点0 = t 1 t 2 t p = “并记 ，= ，1t 2 ，p 7 若给定向量方程数值解痧= u ，致，) ，将( 1 1 2 ) 代入( 1 1 1 ) 运用g a u s s 9 半。 叭一 第一章预备知识 积分离散并表示成矩阵的形式为 户( 磊+ ( 心 l ) 矿，衫，f ) = o ( 1 1 3 ) 其中s 是谱积分矩阵，j 是n x n 的单位矩阵。 应用谱积分矩阵s 离散( 1 1 2 ) 得到 面= 碗+ ( a t s 0 1 ) o ( 1 1 4 ) 从而得到原方程的数值解。 将原方程组的求解转化为对( 1 1 3 ) 式的求解问题。首先运用低阶方法给定方 程的近似解矿= 玩，玩，厅， 7 ，玩= 彰，彰扪，彰) 是维行向量， = l ，2 ，p ，d ( f 表示在g a u s s 点f l ，t ，处的低阶近似解，活l ，2 ，n 。根据误差 占( ，) = u ( f ) 一厅( f ) 定义微分代数方程的误差方程 f ( u o + f ( 疗( f ) + 6 ( f ) ) 妣厅( f ) + 6 ( f ) ，f ) = o ( 1 1 5 ) 其中疗( ，) 是d 的多项式插值函数。 同样运用低阶方法得到近似误差否= 五，覆，瓦 2 ，乏= 彭，碍2 1 ，霉 是维 行向量，离散( 1 1 5 ) 式得到 户( 玩+ ( 缸譬o l ) 矿+ ( 腐 l ) 否，疗+ 否，7 ) = 口 ( 1 1 6 ) 将上式整理写成关于近似解扩的函数如下 否= 疗( 厅) ( 1 1 7 ) 求解原方程的解转换为求解衫使得否；0 ，于是考虑方程 膏( 驴) = o ( 1 1 8 ) 根据隐函数定理，膏的雅可比矩阵，时为 厶= 要= 一( 蓦+ 誓必 - l ( 蓦+ 誓螂 = 一n ( 蓦+ 誓旃 - i ( 誓( 西一心) = 一“c - 9 ， 其中a 尹a 厅是非奇- 异的，因为j 是s 的逼近，当a t 很小时，厶趋近于一，。 运用n e w t o n k r y l o v 方法求解方程( 1 1 8 ) 。先用n e w t o n 法得到方程的收敛解 u 。l = 玑一触 l o 青岛大学硕士学位论文 其中阳为线性方程，( 玑) 阳= 日( 玑) ，问题转化为对线性方程组的求解。在k r y l o v 子空间e ( ，( 玑) ，日( 玑) ) = 印册 日( 玑) ，( 虬) 日( 玩) ，( 玩) 日( 矿，) ) 中搜寻 方程，( 玑) 舾= 日( 玑) 的最优解， 运用g m r e s 方法使其满足 ，埘= ，( 玩) 刃一日( “) 与k r y l o v 子空间正交，或者使乇在厶中的范数最小。 从以上加速过程中求得数值解疗，带入方程( 1 1 4 ) ，得到原方程的数值解。 对于非线性常微分方程初值问题，形如 豁) = 厂( f ， ( ) ) ，【o ，t 】 ( 1 2 0 ) u ( o ) = u o 其中u o ，u ( t ) ec ，f ：r x c 一c 其离散过程与微分代数方程的思想类似，不同点在于运用残差思想得到解的误差的 迭代校正方程，同样运用n e w t o n - k r y l o v 方法求解，具体过程见 3 4 3 5 3 4 8 。 第二章弹性杆s c h r s d ln g e r 方程的谱方法 2 1 圆截面弹性杆s c h r 6 d i n g e r 方程 基于刚性截面假设的k i r c h h o f f 方程m 矧是描述弹性杆平衡方程的经典模型方 程，其形式为 竺+ 劬y 一锡= 0 一+ 劬y 一锡2 j d p _ + 伤口一q y ：0 一 = + 伤口一q y = d _ z 。+ q 夕一纰口= 0 爿皇兰i + ( c b ) 6 0 2 国3 一c 硝魄一，芦= o b 孥+ - 一c ) c 0 3 t o , + c 硝q + f a = o d j c d ( 。舀3 + ( b a ) c o , c 0 2 ：o 其中各截而内力的主矢，为常量，f 为其模。口，7 为f 方向坐标轴与( p - x y z ) 各轴 的方向余弦，a ，b 为截面绕x 轴和y 轴的抗弯刚度，c 为截面绕z 轴的抗扭刚度，衅 为弯扭度q ( i = 1 , 2 ，3 ) 在松弛状态下的原始值。该方程是以q ( j ) ( 净1 , 2 ，3 ) 作为确定 曲杆几何形态的3 个独立变量而建立的方程。 1 9 9 4 年s h i ，h e a r s t 阱_ 引等对于圆截面杆的特殊情形，将k i r c h h o f f 方程的变量 解耦，化作以曲率| r 、挠率z 和相对扭角z 为未知变量的非线性s c h r o d i n g e r 方程， 形式为 窘一t 詈) 2 - - c 。i ( + - - 鲁 汜- ， r 警+ 2 誓( f 一守鲁 晓2 ， 出出k2 ， 彳 譬= 锡。一矿 ( 2 其中系数c 定义为 1 2 青岛大学硕士学位论文 h 朋2 q 5 j 下 h = ( 2 h 彳) 一五盛，皑。= 坞( 对于圆截面杆) ，h 为哈密顿函数，肼= 名( 伤。一趟) ， 见= c 4 ，石，石为单位长度接触力厂沿法线轴和副法线轴的投影。j r ，f 分别为曲 杆中心线的曲率和挠率。 2 2 边值问题的g a l e r k i n 谱方法 运用第一章中谱离散，将方程( 2 1 ) ( 2 2 ) 离散。考虑圆截面螺旋杆s c h r 6 d i n g e r 方程，方程满足边界条件 r ( o ) = r ( 2 万) ，f ( o ) = r ( 2 x ) 于是方程方程( 2 1 ) ( 2 2 ) 组成一个非线性的微分方程组 生一t2)2-t,-+!-,-3ds2 e l k + 2 = 鲁一一气卜。芎 r 韭+ 2 坚f f 一竺1 ：五 (isd s 、2 a r ( o ) = r ( 2 万) ，f ( o ) = r ( 2 ，r ) 选取谱基函数 p ，b ，k = 0 ，l ，2 ，) ，运用g l e r k i n 方法n 6 4 7 馘侧将方程组离散 ， ( s ) = 盎扩， k = - 疋i2k ( s ) 一罢= n 蟊p 出，t 。= 瓦 厶t = 一n 将上式代入方程组( 2 1 ) ( 2 2 ) 并与 p ”，七= o ，l ，垃，j 作内积，得到方程弱形 式 订等一卟一手+ 3 卜= 。 巾警+ 2 警( 一一詈) 卜= 。 p = 0 ，l ，2 ， 由谱积分离散 r 疗等庐山卅删。荟。毫 r m 一册2 奄嘶 r 石p 秘出= 2 刀血 由 一 r 5 吾( ) 3 e 枷出= 万奄氏气 一 描心 r 石等秘出划。厄雹弓 r 万等( 扎詈户出= 摊。+ l + p - - o 蠢弓 将上述结果代入方程弱形式得 一2 ：x k 2 碡一2 万反毛z 1 2 2 耽i 碡+ 万丘氏气= o厶一 一o ，j一 二一 一一而。一 一 胪” + k + p l = o t + l z“舢z 等乜 2 m 反蟊+ 4 疵血毒= o一一“-一 通过以上代数方程的求解，得到r ，f 的数值解。将了的数值解代入方程( 2 3 ) 可 以求得相对扭角z 的数值解，再将运用1 c ，f 和z 的解代入( 1 4 ) 一( 1 6 ) 式可以 得到旋转矩阵q ( s ) 的值，【l i 方程组( 1 7 ) 计算得到圆截面杆的几何形态。 2 3 初值问题的加速谱延迟校正方法 考虑一般弹性杆的k i r c h h o f f 平衡方程的初值问题，将方程组( 2 1 ) ( 2 2 ) 的 加上边值条件r ( o ) = r o ，石( o ) = ，形式为 坐一r(f一2-cik-+ds2 吉矿= 砉一吖r j 2 着 r 韭+ 2 韭一一旦1 ：五 d sd sl 2 a r ( o ) = 岛，v ( o ) = ，f ( o ) = 令v ：_ d g ，方程组化为 生r(f一)2+clk-ds 21 2 矿+ 五a = r lf 一 一+ 二卫 f 1 4 青岛大学硕士学位论文 韭=孚(f一舟-4五-vds c a 一= 一l 。r 一一l - 一 誓l2 ， d x 一= 1 ， 出 r ( o ) = ，v ( o ) = x o ，f ( o ) = 方程转化为一阶常微分方程的初值问题。令口= ( e l l ，) ，u o = ( r o ，r o ，r o ) ，方程可写 成( 1 1 9 ) ( 1 2 0 ) 的形式。运用加速谱延迟校正方法，先给定方程的近似解”o ，记 矿( s ) 为其l a g r a n g e 插值函数，根据由p i c a r d 积分得到的误差方程 占( s ) = r 巾，l f o ( f ) + 占( f ) ) 一巾，矿( f ) ) 出+ ( s ) ( 2 - 4 ) 其中 ( s ) = + f ( 柚。( f ) ) 如一u o ( j ) 当8 ( s ) 很小时，( 2 4 ) 式近似为 占( s ) = r 无。( 删。) 山+ ( s ) ( 2 5 ) 其中j o 是厂( f ，矿( f ) + 占( ，) ) 的j a c o b i a n 矩阵。 运用g a u s s i a n 节点离散( 2 5 ) 式，并运用谱积分矩阵写成 ( ，一a t s j ) g = s 匕式可以应用g m r e s 方法加速，得到校正误差，从而得到方程的数值解。 2 4 数值计算结果 对于圆截面弹性杆s c h r s d i n g e r 方程给出的不同边值和初值条件，选择相同的物 理参数，分别运用g a l e r k i n 谱方法和加速谱延迟校正方法进行计算，得到较好的数 值结果。 2 4 1 边值问题数值计算结果 对于方程组( 2 1 ) ( 2 2 ) 中的物理变量选取m = l ，c i = 0 9 ，以= 五= 0 ，基 函数选择= 4 。求解方程组得到弹性杆中心线曲率r 和挠率f 的计算结果如图l 所 示，将他们代入微分方程方程得到方程的误差如图2 ，从图中可以看出其误差精度 达到1 0 一，图3 描述的是该方程组所决定的圆截面螺旋杆的几何形态。 1 5 图2 1 曲率r 和挠率f 的数值解 0 图2 2 原始方程代入数值解之后的误差 图2 3 对应的弹性杆 1 6 青岛大学硕士学位论文 从得到的数值结果可以知道，运用f o u r i e r 函数系作为试探函数和检验函数的 g a l e r k i n 方法作为微分方程的一种数值解法可以得到较好的精度，例中取基函数的 阶数为4 ，由于方法的稳定性，可以取更高的阶，得到更好的数值结果。因此基于 f o u r i e r 函数系的g a l e r k i n 方法用于解决弹性杆的平衡方程非常有效。 2 4 2 边值问题数值计算结果 对于初始方程，选择物理变量册= l ，q = o 9 ，五= 五= 0 ，并给定问题的初值 r ( o ) = 岛，v ( o ) = s o ，f ( o ) = f o ，弹性杆弧长j = 2 万。大步长选择厶= 2 石，g a u s s 节点选择p = 1 5 ，计算得到弹性杆中心线曲率彭和挠率f 如图2 4 所示，并与四阶龙 格库塔方法比较数值解的误差如图2 5 所示 图2 4 曲率r 和挠率f 的数值解 图2 5 运用龙格库塔法解的误差 图2 6