（应用数学专业论文）反应精馏过程的一种数值模拟方法.pdf

摘要 反应精馏过程的一种数值模拟方法 摘要 反应精馏( r e a c t i v ed i s t i l l a t i o n ，简称r d ) 是化工过程强化的重 要方法之一，其数学模型一般为一具有强非线性的微分代数方程组， 用传统方法进行数值模拟较为困难。这篇论文主要目的就是建立反 应精馏过程的动态模型以及相应的数值求解方法，并对不同的求解方 法进行比较，以确定适合于所研究体系的数值方法。 r d 过程的平衡级模型是半显示2 指标的微分代数方程组。指标 越高，刚性就越强，寻找其数值解也就越困难。为了避免数值计算的 困难，本文用降指标的方法减少了微分代数方程组的刚性，从而为求 其数值解带来了方便；同时，采用变阶变步长的向后差分法( 一种单 支方法) ，既减少了计算时间，又能保证有高阶的精度。此外，由于模 型具有很强的非线性性，必须考虑其稳定性，因此文章还将单支方法 的曰一理论推广到微分代数方程组中，并得到了相应的b 一收敛结果。 应用本算法对m t b e 反应精馏系统进行模拟，将实验数据与上 述算法的模拟结果进行比较。结果表明：模拟结果与已有文献中的数 据吻合良好，这说明此算法是可靠的，并且计算结果能具有足够的 精度。另外，该算法占用内存少，收敛速度快，稳定性好，可方便地 在微机上进行r d 过程的模拟和分析。 关键词：反应精馏，指标，微分代数方程组，向后差分法 摘要 am 匝t h o d0 fn i m e r i c a ls i m u l a n o nf o r r e a c t i v ed i s t i l l a t i o n a b s t r a c t r e a c t i v ed i s t i l l a t i o ni so n eo f 山ei m p o r t a n tm e t h o d so fp r o c e s s s t r e n g 也e n i n gi nc h e m i c a le n g i n e e r i n g g e n e r a l l y ，m a t h e m a t i c a lm o d e lo f r e a c t i v ed i s t i l l a t i o ni sa d i e r e m i a l a l g e b r a i cs y s t e m w i t h s t r o n g n o n l i n e a r i t y i ti sd i m c u l tt on u m e r i c a l l ys i m u l a t et l l ep r o c e s sb yt h e 舰d i t i o n a la l g o r i t h ms t h i sp a p e ri sc o n c e m e dw i ma l g o r i t h mf o rt h e n u m e r i c a l i m e g r a t i o n o fs t i f f d i f f b r e m i a l - a l g c b r a i cs y s t e m 7 n l e u n s t e a d y - s t a t ee q u i l i b r i u ms t a g em o d e l i sb u i l ta n dm e t h o d sf o rs o l v i n gt l l e m o d e la r ei n v e s t i g a t e da n dd e v e l o p e d t h e s em e t h o d sa f ec o m p a r e dt of i l l d b e t t e rs u i t a b l em e t h o d sf o rm es y s t e m t h e e q u i l i b r i u ms t a g e m o d e lf o m sas y s t e mo f s t i 行 d i 腩r e n t i a l - a l g e b r a i ce q u a t i o n s ( d a e s ) o f i n d e x2 ni sw e l lk n o w nt h a tt h e s y s t e mw i t hh i g h e ri n d e xi sa ni i l - p o s e do n ea n dt h ei n d e xo fas y s t e mo f d a e si sr e g a r d e da sa ni n d i c a t i o no fc e n a i nd i 珩c u l t i e sf o rn u m e r i c a l m e t h o d s s ow ef i r s tp e r f o r ma na i l a l ”i c a l i n d e xr e d u c t i o nt oa v o i dt l l e s e 摘要 d i f f i c u l t i e s t h e nw eu s ean u m e r i c a lm e t h o db a s e do nt h eb a c k w 莉 d i 仃e r e n t i a t i o nf o m u l a e ( o n eo fo n e - 1 e gm e t h o d s ) w i mc h a l l g e a b l eo r d e r a n ds t e ps i z e ，w h i c hc a l ln o to n i yr e d u c ec o m p u t e rt i m e ，b u ta l s oo b t a i n e n o u g hp r c c i s i o n ，t o s o l v et h e s y s t e m i na d d i t i o n ，b e c a u s eo fs t r o n g n o n l i n e a r i t yo ft h es y s t e m ，i t sn o n l i n e a rs t a b i l i t y 枷s tb ec o n s i d e r e d s o ， s p e c i a l l y i nm i sp 印e r ，口一t h e o r yo f0 n e l e gm e t h o d s f o rs t i f fi n i t i a lv a l u e p r o b l 锄si se x t e n d c dt o s t i f rd i 脑e n t i a l - a l g c b r a i ce q u a t i o n sa 1 1 ds o m e c o r r e s p o n d i n gr e s u l t so f 曰一c o n v e r g e n c e a r ep r e s e n t e d am t b er e a c t i v ed i s t i l l a t i o nc o l u m nm o d c li ss i 脚l a t e da s a n u m e r i c a le x a l p l e t h en u m 舐c a lf e s u l t sm a t c hw e l lt h ed a t ai nm e r c f _ e r e n c e t 1 1 ee x 锄p l ei 1 1 u s t m t e st h a tt h ea l g o d t l l i i li sq u i t ee m c i e n tf o f s t i f fi n i t i a lv a l u ep r o b l e m sa n dc o m p u t a t i o n a lr e s u l t s 诵t l le n o u g hp r e c i s i o n c a nb eo b t a i n e db yi t m o r e o v e r ，m ea l g o 删1 mn o to n l yo c c u p yl e s se m s m e m o b u ta l s oo w n m o r er a 【p i dm t eo fc o n v e 唱e n c ea n d b e t t e rs t a b i l i 够i t i sq u i t ec o n v e n i e n tf o ru st os i m u l a t ea i l da n a l y z er dp r o c e s sb yt h e a l g o r i t h mw i t hp e r s o n a lc o m p u t e r s k e yw o r d s ：r e a c t i v ed i s t i l l a t i o n ，i n d e x ， d i 虢r e m i a l - a l g e b r a i c e q u a t i o n s ，b a c k 、7 l ，a r dd i 恐r e 而a t i o nf 0 咖u l a e 符号说叫 符号说明 组分数 塔板数 液相摩尔分数 汽相摩尔分数 相平衡系数 温度，k 汽相流率，m 0 1 s - i 液相流率，。1 s - l 进料量，。1 s _ l 进料比值 滞留量，m 0 1 s - 1 化学计量系数 反应体积 回流比 焓值，c a l m o l - f 第n 级塔板的流液量，。0 1 s _ l 反应速率，m 0 1 s - 1k g _ 1 活度 偏心因子 理想气体热容，c a l m o r i f 一 临界温度，k v i c n b 弓 巧 0 c 钆 q 0 r 彤 叮 甜 q 嚷 i p “蒸汽压，p a f时间，s ，彳 侧线流量与级问流量的比值 上角标 矿汽相 工 液相 f进料 o标准态 下角标 ，塔板号 f组分号 m第m 个反应 北京化工大学 学位论文原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师的指 导下，独立进行研究工作所取得的成果。除文中已经注明引用 的内容外，本论文不含任何其他个人或集体已经发表或撰写过 的作品成果。对本文的研究做出重要贡献的个人和集体，均已 在文中以明确方式标明。本人完全意识到本声明的法律结果由 本人承担。 学位论文作者签名：王玉花 2 0 0 6 年6 月2 日 北京化工大学硕士学位论文 第一荜概述 1 1 引言 第一章概述 反应精馏( r e a c t i v ed i s t i l l a t i o n ，简称r d ) 是化工生产中的重要分离手段之一。 这种将反应和精馏设备结合在一起的生产过程，彻底改变了长期以来人们对反应 和分离过程的传统认识，它将化学反应过程和精馏分离的物理过程结合在一起， 是伴有化学反应的新型特殊精馏过程。r d 概念自1 9 2 1 年由b a c l d l a u s 提出以来， 已经历了从3 0 年代到6 0 年代对特定体系的工艺探索，相对于传统的先反应后分 离过程具有节省资金、提高转化率、易操作、提高选择性、避免形成共沸物、避 免形成热点和飞温等一系列优点。但与此同时，r d 过程也受到挥发性限制、停留 时间限制、过程条件不匹配等困难。因此，对r d 塔的硬件设计要有严格的要求： 要使液体和催化剂颗粒有效接触、要使反应域汽液接触良好、要使反应段有充足 的持液量等条件。反应精馏过程因其内部各因素之问的强非线性耦合，导致其设 计与操作问题要比传统的反应器或精馏塔复杂得多，目前仍未形成完整的理论体 系。所以，当前有关r d 的理论研究和应用研究已成为各国专家们的研究热点。 由于化学反应和精馏之间存在着非常复杂的相互影响，物系之间往往呈现非 理想性，这使得操作条件和设备参数的微小变化对反应精馏的操作规律都有很大 影响。所以，要发挥反应精馏的优势，必须选择合适的设备和工艺条件，以达到 反应和精馏的最佳配合。由于反应精馏特有的复杂性，在反应精馏过程的设计、 放大、操作和控制方案的研究等方面存在一定的困难，人们对反应精馏规律的认 识，尚有待进一步的研究。因此利用数学模拟的方法，对这种装置进行模拟分析， 深入探讨反应精馏过程的生产规律，优化操作条件和设备参数等都是非常必要的。 自2 0 世纪7 0 年代以来，有关反应精馏的研究重点已从工艺转向数学模拟，目前 对反应精馏过程的数值模拟，其模型大多都是描述精馏过程的m e s h 方程组结合 反应速率方程r 而建立的。根据求解这些非线性方程组的方法的特点，可大致分 为泡点法、松弛法、牛顿拉夫森法、简捷法、分块法、解离法、同伦连续法等等。 北京化t 大学顾。l ：学位论文 第一章概述 前两种方法的模型简单，计算方便，占内存少，对初值要求不高。但在求解反应 精馏的高度非线性方程组时，常常难以收敛。所以该领域的大量近期文献广泛地 采用n c w t o n 法( 或它的变形) 对所有的独立方程进行同时求解，读者可以参考文 献”以获得更多的细节以及说明性的实例。n e 州o n 法收敛速度快，但占内存较大， 计算j a c o b i 矩阵耗时长，对初值也有严格的限制。因此寻找更先进的数值解法具 有更重要的意义。 1 2 相关领域的历史、现状和前沿发展情况 对r d 过程的模拟研究一直是反应精馏这一领域的热点问题之一，其可分为定 念模拟和动态模拟。数学模型则是对该过程进行模拟和理论分析的基础。目前主 要有两大类模型平衡级模型和非平衡级模型，以下将简要回顾文献资料中所 出现的数学模型及其模拟方法。 1 2 1 平衡( e q u i l i u m ，简称e q ) 级模型 r d 过程的平衡级模型是在传统精馏过程平衡级模型“。的基础上发展而来的， 其仅需在后者的液相物料衡算方程中添加用以表示反应的项即可。对r d 过程的模 拟尽管存在多种计算方法，但所采用的平衡级模型基本上是相同的，即假设各塔 板上汽液两相理想混合，离开塔板的汽液两相达到平衡，且反应仅发生t 液相。 与一般精馏类似，r d 的平衡级模型”用m e s h 方程组模拟平衡级塔板。系统设备 如图所示 图1 1 第层塔扳 f i g 1 一l1 h ee q u i l i 撕ms t a g c 2 北京化工大学硕士学位论文第一章概述 h 方程代表总能量衡算方程 警= h h 。吣f j h ：舢似h ：一q + 批h ( 1 5 ) 上式中的u ，是第，块塔板上的滞留量，通常认为u ，是液相的滞留量，仅在高 压的情况下才考虑汽相的滞留量；玉，是液相摩尔分数；是汽相摩尔分数；t ，是 相平衡系数；l 是汽相流率；l 是液相流率；c 是进料量；弓，是进料比值；h m 代表第f 个组分在第m 个反应中的化学计量系数；占代表第_ ，块塔板上的反应体 积；r ，是反应系数；彤，彰，彤分别代表汽相，液相和进料的焓值，其中下标 _ ，表示第，级塔板；f 表示第f 个组分；m 表示第m 个反应。 上式中左边时间导数项中的焓值代表塔板上总的焓值，一般仅取液相中的焓 值即可。有些作者在能量衡算方程中还添加了用于表示反应热的项。如果焓值是 以各组分的元素状态为基准的话，则反应热就被自动计入焓值中而无需用单独的 项来表示。 在稳态条件下，上述方程组中所有时间导数项均等于零。 另外模型中还包含反应动力学方程，对于缺少动力学方程的快速液相可逆反 应可用化学平衡方程代替动力学方程，对于慢反应，可用化学平衡来估算反应可 能进行的最大程度。除此之外，模型中还包括汽液相焓和相平衡常数关系式。 1 2 2 非平衡( n o n o q 试1 姗砌，简称n e q ) 级模型 为了更加完整地描述r d 过程，有必要详尽地考虑传递过程与化学反应之间的 相互作用。r d 过程的非平衡级模型是在传统精馏过程速率级模型的基础上发展而 来的。然而建立该过程的非平衡级模型并不像为之建立平衡级模型那么直截了当， 首先必须确定反应是均相反应还是非均相反应。整个非平衡级模型包括汽液相物 料衡算和能量衡算方程，汽液相传质与传热通量方程，相界面平衡关系式，归一 化方程，汽液相焓及相平衡常数关系式，反应速率和汽液相传质传热系数的特征 方程。而且在模型中还应包含r d 塔的水力学关系式。 4 北京化工大学硕十学位论文 第一章概述 1 2 3 髓过程的传统的数值模拟 对于稳态的平衡级模型，绝大多数有关r d 模型化的早期文献均致力于开发 用于求解稳态e o 级模型的方法。火部分这些方法都是或多或少地由求解传统精 馏的方法直接推广而来的，我们把用于求解e q 级模型方程组的基于计算机的计 算方法分成几类： 简捷法( s h o n c u lm e m o d ) 需要作大量的简化假设，以导出便于快速计算的 形式简单的近似方程组。由于在许多方法中过程都受到化学反应的影响，因而很 难得到适用于r d 过程的通用简捷法程序。有关文献。1 表明，在简捷法中常用的 化学平衡假设对于大多数r d 过程是不恰当的。 该领域的大量近期文献广泛地采用n e w t o n 法”1 ( 或它的变形) 对所有的独 立方程进行同时求解，尽管收敛速度快，但对初值却有严格的限制。 松弛法( r e l a 】【a t i o nm 劬o d ) 将m e s h 方程组写成非稳态的形式并对其进行数 值积分直至达到定态。k o m a t s u ”1 应用了这种方法并将平衡级模型的计算结果与 他的一个实验结果进行了比较，结果表明平衡级模型的组成分布从定性上来讲是 正确的。由于松弛法需要大量的计算时间，因而在实际中并不被经常使用。 同伦连续法( h o m o t o p y - c o n t i n u a t i o nm e t l o d ) 经常被用于求解用其它方法( 如 n e 、t o n 法) 难以解决的问题。读者可参考w j y b u m 和s e a d e r “1 撰写的文献以获 取该方法的相关信息。c h a l l g 和s e a d e r 。1 首次将这种方法应用于r d 操作。此外， 同伦连续法对参数敏感性研究以及定位多定态也是十分有效的。 a l e j s k i ，s z y m a n o w s h 和b o g a c 】【i 应用p o w e l l 开发的极小值法( m i l l i l n i z a t i o n m e l o d ) 求解了m e s h 方程组，但这个方法收敛的速度很慢。 另外，有些作者”1 1 引入变量的概念，将发生平衡反应的反应精馏过程模型变 换为普通精馏过程模型，减少了迭代变量，避免了反应量的计算，并易于采用普 通精馏的计算方法。并针对新模型的非线性强的特点，以松弛法得到初值，以n r 法为算法主体，对反应精馏过程进行模拟计算。 此外，还有些作者“甜采用了一种改进的平衡级模型，其中通过在相平衡方程 中引入板效率因子来补偿实际偏离平衡状态的误差。然而研究表明，板效率因子 北京化工大学硕士学位论文第一章概述 是十分复杂的，且化学反应的存在将对其产生较大影响。因而无法准确地估算板 效率因子，这就使平衡级模拟的可靠性受到很大的影响。 对于非平衡级模型，由于在非平衡级模型中考虑了汽液相传质与传热阻力， 方程个数大大增加，同时模型方程组较平衡级模型有较强的非线性，因而给求解 工作带来了很大的困难。用松弛法解非线性方程组难以收敛；用n e 叭o n r a p h s o n 法解非平衡级模型，对初值的要求又很高。实际计算表明，在r d 过程的计算中初 值很难确定，即使对简单反应物系的计算都难以收敛。把松弛法和n 喇o n r a p h s o n 法结合起来是求解非平衡级模型的一条有效途径。此外，还可采用具有极强的全 局收敛特性的同伦连续法求解非平衡级模型，通过适当地选取初值，可用该法找 到系统的多解乃至全部解。虽然该法的计算效率不如n e 叭o n 法，但它的可靠性极 好，并且同伦连续法对系统的参数敏感性研究也是十分有效的。目前国外建立了 以同伦连续法为基础的非线性动态系统，用于研究r d 过程中的传质现象，并借助 a u t o 数学软件简化非平衡级的模拟计算”“。 1 3 本课题研究的内容 本课题建立反应精馏过程的动态模型以及相应的数值求解方法，并对不同模 型以及求解方法进行比较，以确定适合于所研究体系的数学模型和求解方法，使 其计算效率较高，并且做到在满足条件较少的情况下，算法有较强的收敛性。 6 北京化工人学硕士学位论文第二章刚性问题和指标问题 许多文章m - 1 6 1 中的算法都是将代数方程和微分方程分开，先利用组分物料衡 算方程组计算液相摩尔分数再使用泡点法由加和方程计算各级塔板温度巧， 最后再将t i 代入总物料衡算方程和能量守恒方程计算汽液相流率巧，a 本 文从降低方程组的刚性着手，建立一种新的算法。 2 2 指标问题 2 2 1 微分代数方程组与指标 上述方程组( 2 1 ) 一( 2 4 ) 既含有微分方程，又包含代数方程。为了更好 的求其数值解，我们先来看下面的定义m 1 ： 对于方程组 f ( f ，y ，y 7 ) = o ，f c ( r 2 肿1 ，r 4 ) ，f o ，丁】 ( 2 5 ) ( 1 ) 当兰是非奇异的，则称式( 2 5 ) 是常微分方程组( o r d i n a r y d i f r 玎e n t i a l 卯 e q u a t i o n s ，简称o d e s ) ； ( 2 ) 当兰是奇异的，则称式( 2 5 ) 是微分代数方程组( d i 自f e r e n t i 小a l g i b r a i c c e q u a t i o n s ，简称d a e s ) ； ( 3 ) 特别地，当兰= o 时，则称式( 2 5 ) 是代数方程组( a l g e b r a i ce q u a t i o n s ) 。 咖 d a e s 与o d e s 之间有密切的联系，但两者之间也有很大的差别。首先来看 一下d a e s 的指标：伴随一个解，( f ) 的指标，是指表示微分代数方程组分离程 度的一个诈整数。这个数是使得能从超定方程细 f ( f ，y ，y ) = o 霉( ，y ，y 一) ：o 以( 2 6 ) 警y “一1 ) ) = 。 8 北京化工人学硕上学位论文 第二章刚性问题和指标问题 中唯一地解出_ ) 以至y 能直接依赖于y 和f 的最小的正整数p 。很显然，d a e s 的 指标都大于等于1 ，也可视0 d e s 的指标为0 。通常地，来源于实际问题的这些 d a e s 具有半显示的形式，即 ( 1 ) 半显示l 指标的h e s s e n b e r g d a e s 形式为 ；竺象等2 嚣l 。；善幂：j 吲。，卅 c z ， 1 0 = g ( x p ) ，( f ) ，f )g c ( r “州，r 2 ) 。 其中孥是非奇异的。 砂 ( 2 ) 半显示2 指标的h e s s e n b e r gd a e s 形式为 僻裟搿o x o ：i 三幂并吲。用s ，lo = g ( x ( r ) ，f )g c 。( 足”十1 ，r ) 其中挚学是非奇异的。 ( 3 ) 半显示3 指标的h e s s e n b e r g d a e s 形式为 i 工o ) = ，( 工( f ) ，y ( f ) ，z ，f )，c 1 ( r “+ 七“+ 。，r ”) y ( f ) = g ( 工( f ) ，y ( f ) ，f )g c 2 ( r “+ “1 ，r ) f e o ，丁】 ( 2 9 ) io = ( y ( f ) ，f ) c 3 ( r ”州，尺7 ) 其中罢冬芒是非奇异的。 很显然我们的问题( 2 1 ) 一( 2 4 ) 就是一个半显示2 指标的h s e n b e r g d a e s 。 d a e s 的指标不仅能够衡量系统的奇异程度，而且还能够暗示其数值解的难易程 度。一般地，指标越高，刚性就越强，问题也就越病态，寻找其数值解也就越困 难。为了避免数值计算的困难，可以考虑采用降低微分代数方程组的指标来减少 其刚性的算法求解，通过降指标的方法将其化为1 指标的微分代数方程组。 2 2 2 降方程组的指标 用能量守恒方程求出2 至n 一1 级塔板温度导数，并将加和方程变为 c ( 毫。一) = o ，对其求导变形，即求得汽液相流率的代数方程和第l 级第n 级 f 刊 9 北京化t 大学硕士学位论义第二章刚性问题和指标问题 塔板温度的导数方程。这样m e s h 方程组就化为半显示l 指标的微分代数方程组 即 ! = ( _ + ，t 。+ + + 川_ 川十一z 。一( 1 + ) _ i 。气， 一( 1 + 哆) ，x 。，+ v ，r 。，s ，) u ， f = 1 ，2 ，cj = 1 ，2 ，n 盟：熟! ：! 翌 出 眠，t ， 刍d 正 。 峨w 。h k 儿。h k + f | h ：一y i h j j h u i 一毛瓮i j d 警 出 型 鲁a l ， ，22 ，3 ，n l 盟：堑! ：! ! 兰 d t 母垫n ， 台d 氏“ o = 0 + 。+ 上川+ 一一( 1 + ) 0 一( 1 + 哆) 工，+ q ，。只。s ， 一“一 ，= 1 ，2 ，n o = 上一只矿 。= 私川鲁+ 砉等鲁 ( 2 1 0 ) 这样式( 2 1 0 ) 就具有式( 2 7 ) 的形式，由于系统具有很强的非线性性， 再将其化为o d e s 比较困难，所以我们直接对( 2 1 0 ) 式进行离散化，在离 散化之前要考虑问题的刚性，再选择合适的数值算法进行求解。 2 3 刚性问题 在航天、航空、热核反应、自动控制、电子网络及化学反应的过程等许多科 学技术领域及实际问题中，经常遇到可用常微分方程描述的物理和化学过程，这 些过程往往包含着多个相互作用但变化速度相差十分悬殊的子过程，这就导致相 应的常微分方程的解中既包含衰减十分迅速的分量，又包含有相对来说变化缓慢 的分量。当人们试图在解的慢变区间上数值求解这类问题时，尽管此时快变分量 1 0 北京化t 大学硕士学位论义 第二章刚性问题和指标问题 已衰减到微不足道，但这种快速变化的干扰仍严重影响数值解的稳定性和精度， 给整个计算带来很大实质困难，人们称上述这类过程具有“刚性”，描述这类过 程的常微分方程初值问题称为“刚性问题”。 定义2 1 2 叫初值问题僻竺o ：第 ( 2 【y t 町= 刁，叩“ 称为在区间k ，6 】上是刚性的，如果 ( 1 ) 该问题的解y ( f ) 在区间- ，6 】上是变化缓慢的( 这里所谓“变化缓慢”是 相对于善( f ，y ( 功的实部是模很大的负数的那些特征值而言的) ： 砂 ( 2 ) 对于解曲线上的每一点( f ，y o ) ) ，口f 6 ，j a c o b i 矩阵譬( f ，y ( f ) ) 至少有一 叫 个特征值丑其实部是模很大的负数，且所有特征值的实部都不取大的正值。 实际积分刚性系统时，人们感兴趣的是稳态解，因此至少要积分到暂态解衰 减到可以忽略不计的时刻。若试图在整个计算过程中使用显式方法，具体的说， 比如使用定步长显式e u l e r 法，则数值稳定性要求步长 o 是积分步长。s 级 = i r 1 1 1 1 9 e - k i l t t a 方法具有单步方法的固有优点：( 1 ) 计算是自开始的，可直接利用原 问题的初值作起始值，不需要计算附加起始值；( 2 ) 易于变步长计算，不需要像 线性多步法那样化成n o r d s i e c k 形式或缩小步长时进行插值；( 3 ) 稳定性恒等于1 。 许多r u n g e _ k u t t a 方法具有处理刚性问题所需要的高稳定性和高精度，这方面远远 优于线性多步法。然而不幸的是一般s 级隐式r 1 m g c k 嘣a 方法由于每步需要求解 一个由s m 个( 通常是非线性的) 方程联立的方程组，像通常采用n e w t o n r a p h s o n 迭代求解，其计算量近似的是线性多步法每步计算量的s 3 倍，当s 和m 较大时，计 算量大得惊人。另外在求解非线性刚性问题时，其收敛阶并不高，因此这类方法 目前有较少的实用价值，除非能设法大幅度减少其每步的计算量或加快其每步计 算的速度。 线性多步法 i q y 。= 层，( 。只。) f = lf = l 其中吒，屈是实系数，其突出优点是计算格式简单，每步计算量相对较少，但高 阶线性多步法的数值稳定性不好。高阶线性多步法稳定域的扩大是以损失稳定程 度和计算精度为代价的。强迫其稳定域极度扩大殊不可取，只有根据实际情况恰 当选择合适的稳定参数使在稳定域大小，形状，稳定程度及误差参数诸因素间取 得恰当平衡，才有可能收到较好的计算效果。所以从根本上改善高阶方法稳定性 的唯一出路是跳出线性多步法的狭窄圈子去构造高效的新型算法。 线性多步法的孪生姊妹单支多步法 北京化t 大学顶士学位论文第二章刚性问题和指标问题 q 只。= 矿( 层。，屈m 。) f = li = li = i 其中哆，层是实系数，与相应的线性多步法具有完全相同的线性稳定性，然而在 非线性稳定性方面，单支方法比线性多步法具有明显优势。我们的问题恰是刚性 的非线性问题，所以选取单支法求解，其中向后差分公式可视为这类方法的典型 代表。 北京化】= 大学硕十学位论文第三章数值模拟 3 1 数值求解 第三章数值模拟 如上所述，我们采用向后差分公式( 一种单支方法) 束求解d a e s ，为了提高 精度和计算速度，我们采用变阶变步长的算法求解。下面我们给出算法以及实现 算法的具体步骤。 3 1 1 向后差分公式 初值问题的向后差分公式 2 1 】( b a c k w a r dd i t i a lf o m m l a ，简记为b d f ) 矽( 。只+ 。) 其截断误差为 = 击矿1 y m i ) ( 伊) 将公式( 3 1 ) 写成通常形式 ( 3 1 ) 一l 虬+ 。= q 儿。+ 展，( 。只+ 。) ( 3 2 ) i ；0 其中 是步长，系数哆，羼可参见附录。 b d f 具有以下优点： ( 1 ) 每步计算仅需求解一个含有m 个未知数的非线性代数方程组，每前进一 步长解隐式方程组所需要的工作量比较小； ( 2 ) 在。点具有极端稳定性，在原点的稳定程度及方法的误差常数处于可接 受范围； ( 3 ) 向后差分法又是单支方法，因而在非线性稳定性和b 一收敛性方面优于 其他任何k 阶k 步隐式线性多步法； ( 4 ) 容易改变阶和步长。 1 4 丛 以一、 j v 一 。爿 型生型生2 兰望! 兰三! 丝堕兰 笙三童塑篁塑型 用边界轨迹法画出方法的稳定域，其中1 2 阶公式一一稳定，3 6 阶公式 爿( 口) 一稳定且s t i 行稳定，高于6 阶的公式因零稳定性被破坏而不能用。 3 1 2 离散化 显然，式( 2 1 0 ) 是个自治系统，为了简化，可以用下面的公式代替 j z = g ( j ，y )g c ( 月啦+ 3 ，月。+ ) “) io = 日( x ，y )h c 1 ( “圳，r 2 “1 其中x 是( c + 1 ) n 维向量，y 是2 n 维向量，g ，灯是非线性的，豢是非奇异 的，将其简写成 m 警一其中z 2 m 2 m “哉。 s ， 为了有效地求解式( 3 3 ) ，采取g e a r 提出的变阶变步长的向后差分法，变阶 变步长的向后差分法能节省计算时间，解刚性的微分代数方程组具有很好的效果。 对于方程组( 3 3 ) ，其离散化后成为 一i m ( 乙+ 。一啊互+ 。) = 展e + 。 ( 3 4 ) i = o 其中只+ t2f ( 乙+ t ) ，式( 3 4 ) 是一个隐式格式，需要用一个显式格式预估乙。， 可用 m ( 乙。一西乙。) = 雕，e + 。 ( 3 5 ) 来预估求出无+ t 的初值乏2 ，再用式( 3 4 ) 进行校正求出乙+ 。其中西，群是 常数，可参阅附录。由于此系统j 。c 。b i 矩阵的模| 筹l 很大，刚性较强，因此，在使 用迭代法求解式( 3 4 ) 时，采用n e 叭o n r a p h s o n 方法，迭代格式为 一l m 乏嚣”一乏乏+ 篡( 墨嚣一q 乙。) 】= 醇冀 孱删 ( 3 6 ) 扛0 她础- ( 1 瑚警) - i 。 北京化t 大学顾十学位论文 第三章数值模拟 为了提高计算速度和精度，我们采用变阶变步长的向后差分法求解，变阶变 步长的算法比较复杂，我们对其变量的存储方式提出以下算法。下面还是考虑标 准的自治系统的半显示1 指标的微分代数方程组 牌赫黜 ；詈茹暑一陬明 慨， i o = 占( z o ) ，y ( f ) ) ，g c ( r ，尺2 ) 初值为并( o ) = 刁，y ( o ) = 彳，满足相容条件0 = g ( 叩，彳) 。 为了变步长和变阶的方便，将方法表示成n o r d s i e c k 2 2 】形式。n o r d s i e c k 提出 储存，蟛，丢，及只，蜕，鲁杉，的近似值来代替储存x ( f ) ，y ( f ) 及其导数在 各个节点上的近似值，它的出发点是使改变步长所需要的附带计算简单一点a 令 以= k ，城，鲁，箬蛰) 】t 和= 眈，蜕，簧一，等) t 。 为了方便，下文中令c 乙，= ( 量 。 首先考虑变阶的方法，对于j j 阶方法，在计算过程中，若解向量的每个分量乙 的最后一个分量为等掣，则变阶为+ 1 时，可用话( 眷一墚：) 代替 ! 生：，这样就解决了阶数增高的问题。 他十1 11 1 下面给出解固定步长和阶次的微分代数方程组的计算步骤： ( 1 ) 使用原方程和初值利用上述变阶方法来估计阶次为七、步长为 的五和 k 的元素 五= 【，城，等，箬搿) 】t 纠腻孙舡) 】t 。 ( 2 ) 利用p a s c a l 三角矩阵d 、凰和k 计算一和i ，墨= 蹦。，i = d k 北京化工大学硕士学位论文第三章数值模拟 其中p a s c a l 三角矩阵d 为一上三角矩阵，它的第( f ，- ，) 个元素是 = 匀 地尘世_ 业，狐 m ! 。 ( 3 ) 若乞= t ，停止计算，输出以和l ；否则转入第4 步，计算以+ 。和。 ( 4 ) ( i ) 令只+ 。= j d 峨，e + 。= d ( i i ) 用n 州o n - r a p h s o n 法求下式得出6 和d j o = 矽( t + 。+ 应。6 或。+ 应。d ) 一( 茸+ 。+ 6 ) 【o = g ( 元+ 。+ 应。6 ，或+ l + 应。d ) ( i i i ) 令爿j + 。= j ：+ + 6 f ，t + ，= e + 。+ 讲；n = 行+ 1 转( 3 ) 。 其中，参数，的值可参阅附录，应。相当于，0 ，6 是m 维向量，d 是七维向量。 为了叙述方便， 令可：件 l d 步长的选取是一个比较复杂的问题，我们希望选取合适的步长以便控制每个 分量的误差，所以步长j l l 应使得不等式 g + ，| j s( 3 8 ) 对于方程组的所有元素都成立其中c + 。= i 鲁。 积分一步并检验( 3 8 ) 式是否成立，若成立则接受这一步；否则拒绝这一步。 对于下一步或重新计算抛弃的步，所用的步长估计为p ，其中p 由等式 q + 。七! “ 来确定。如果采用这个步长，并且误差正好与旷“成比例，则下次检验( 3 8 ) 式 时，( 3 8 ) 式正好满足。但是由于等( 露一础) 并不是不变的，所以为了保证不等 北京化工大学硕士学位论文第三章数值模拟 式( 3 8 ) 能够满足，常采用稍小一点的步长，其中p 由 l 岛2 瓦 来确定。为了变阶，还必须检验在其它阶的公式中所用的步长，即在七一l 和忌+ l 阶 公式中p 分别为 l 成，两 1 n + - 2 丙 占 q 女! 占 g + 2 七! l + 2 一( i + 2 后! “ 盯 在积分过程中，每一步都计算n ，以_ l ，n + i 选取其中最大的p ，以确定下 一步采用合理的步长和阶数。但在具体程序中，每步都改变步长和阶数不一定好， 应作以下考虑： ( 1 ) 如果这次失败，则估计矿； ( 2 ) 在上一次改变步长或阶数以后的七+ l 步，估计p ； ( 3 ) 如果在上次估计p 时步长不放大，则在其l o 步后估计p 。 这样，我们就实现了变阶变步长的向后差分算法。 1 r t哺llll1111，。j kri。llllllij 北京化丁大学硕士学位论文第三章数值模拟 3 2 稳定性分析 对于非线性系统( 2 1 0 ) 的数值稳定性准则，最简单的方法是采用局部线性 化准则，即假设j a c o b i 矩阵孚( f ，y o ) ) 在真解_ y ( f ) 的某个适当领域中是变化缓慢的， 缈 因而可将非线性系统近似等同于常系数线性系统。事实上，实际计算中一般用局 部精度和数值稳定性来控制积分步长，即是原积分步长如此之小，使得局部截断 误差不超过事先给定的误差容限，同时使得 箬( f ，y o ) ) 的每一个特征值不超出方 叫 法的绝对稳定区域。前者保证每一计算步所产生的误差足够小，后者保证各种数 值误差在计算过程中稳定地传播，且逐渐趋向于零。当然，如果使用隐式方法且 每步需要用某种迭代法求解一个非线性方程组时，积分步长还应满足迭代收敛性 的要求。但是，作为问题的另一面，我们必须看到上述实用判稳准则是将非线性 问题近似看作常系数线性系统处理而得到的，这种近似处理在某些情况下可能引 出错误的结果，所以不能作为非线性问题稳定性的理论基础。应寻找更好的理论 作为判别稳定性的基础。 当用高稳定的高阶隐式方法求解非线性刚性微分方程组时，人们指望“高阶” 会带来应有的更高精度，但事与愿违，实际计算精度往往并不像预想的那样理想， 且实际观测阶远远低于经典收敛阶，这种阶降低现象尖锐地表明经典收敛理论不 适合于非线性刚性微分方程组。因此，1 9 8 1 年，f r a n k ，s d m c i d 和u e b 础l u b c r 提 出了曰一收敛理论。 首先，给出单边“口s c h i t z 条件 r e ( ，。一虬，厂( f ，) 一厂o ，y 2 ) ) 圳y l y 2 旷，、吵l ，北r ”，f 疗，6 其中( ，- ) 为有限欧式空间标准内积；为相应的内积范数；，是任一单边l i p s c h i t z 常数。特别地，可取它为最小单边“p s d l i t z 常数f = f 0 以达到可能的最佳估计。由 于任何经典l i p s c h i t z 常数l 均可作为单边“p s c h i t z 常数，故恒有，0 茎。 下面给出口阶最佳日一收敛和b 一稳定的概念。 定义3 1 锄1 方法称为p 阶最佳b 一收敛的，如果以该方法从精确起始值y ( 0 ) 出 9 北京化工大学硕十学位论文 第三章数值模拟 发按定煳求解初值问题髅篙x 二第时，黼懒近删 的整体误差有估计峨一y ( 厶+ 砌) 1 | 兰c ( 屯) 9 ，o o ，非负有界函数，：( o 甜】一r 及实对称正定r ，矩阵g = 【晶 ，这里口，n g 仅依赖于方法，使得对于任给的初值 问题僻篙l 二第及任意两个平衡计蝴厕邯，研瑚和 o 一| j l ，i ) o ，z ，z ，f ) ，有j i y zj l 。( 1 + z y ( z ) ) 2 l l 歹一手0 。，o ，口。其中卜l b 定 义为。= ( u ，g u ) ：。f 主岛( 吩，“。 f ，产i b 一理论摒弃经典收敛理论而采用基于单边l i p s c h i t z 条件的条件估计 | 眵( f ) 一y ( f ) 0 i 旧一叩e x p ( ，( f d ) ) + 等【e x p ( ，( f 一4 ) ) 一1 ，口f 6 其中疗是其扰动问题歹( f ) = 于( f ，歹( f ) ) ，口蔓f 6 的初值，万= s u p i 旷( f ，歹) 一厂( f ，y ) l i 。 r e f 口6 1 ” “e r 。 由于厶三，所以最小单边l i p s c l l i t z 常数f 0 提供的条件估计优于任何经典l i p s c h i t z 常数提供的条件估计。对于许多刚性问题，最小单边l i p s c h i t z 常数仅有适度大小， 甚至可能为负的，此情形下给出的估计是逼真的。因而以它为基础的口一稳定及 b 一收敛理论能反映逼真的误差传播规律及逼真的误差估计。对于非线性刚性问题 的数值解，大量的数值试验实际观测结果表明，在步长的实际范围内经典理论严 重失真而b 一理论是f 确的。 引理3 1 踟单支方法是丑一稳定的，当且仅当它是一稳定的。 引理3 2 删若单支方法是口一稳定的，且经典相容阶为p ( l ，2 ) ，则该方法 是p 阶最佳b 一收敛的，且整体误差有估计 北京化t 大学顾十学位论文 第三章数值模拟 一l i i y 也) 一只忙q ( 乙) 忧) 一川+