并联多通道内核热耦合流动不稳定性的研究毕业论文.doc

目录 并联多通道内核热耦合流动不稳定性的研究毕业论文并联多通道内核热耦合流动不稳定性的研究毕业论文 目录 1绪论 1 1 1流动不稳定性及其危害 1 1 2研究现状 2 1 3研究内容 4 2点堆模型和核热耦合模型 6 3并联多通道理论模型 9 3 1单相水热工水力微分方程 10 3 2两相水热工水力微分方程 11 3 3压降模型 14 3 3 1入口段 14 3 3 2加热段 14 3 3 3上升段 15 3 4无量纲化 16 3 5并联通道模型 18 4数值方法 20 4 1系统的刚性 20 4 2 GEAR 方法 24 4 2 1简单迭代法 24 4 2 2Newton 迭代方法 25 4 2 3Gear 方法 26 4 2 4Gear 算法的阶与步长的控制 28 5程序计算结果及分析 33 5 1不考虑核热耦合时的计算结果 34 5 2核热耦合时的计算结果 37 5 3计算结果对比及分析 40 6结论与展望 42 参考文献 43 附录 45 致谢 78 绪论 1绪论 1 1流动不稳定性及其危害 在加热的流动系统中 如果流体发生相变即出现两相流时 流体以非均匀 形态所出现的大的体积变化可能导致流动的不稳定性 这里所说的流动不稳定 性 是指在一个质量流密度 压力降和空泡之间存在着耦合的两相系统中 流 体受到一个微小的扰动后所产生的流量漂移或者以某一频率的恒定振幅或变振 幅进行的流量振荡 1 这种现象与机械系统中的振动很相似 质量流密度 压 降和空泡可以看作是机械系统中的质量 激发力和弹簧 在这中间 质量流密 度和压降之间的关系起着重要的作用 流动不稳定性不仅在热源有变动的情况 下会发生 而且在热源保持恒定的情况下也会发生 流动不稳定性在反应堆领域广泛存在 流动不稳定性会严重影响反应堆的 安全运行 这是由于 1部件可能遭受有害的强迫机械振动 机械振动和局部热应力周期性变化会 导致疲劳破坏 2会引起控制问题 对于流体冷却的反应堆 当冷却剂兼作慢化剂时尤为 重要 还会引起核反应性变化耦合反馈效应 3 会影响局部传热特性 或者引起沸腾危机提前出现 及临界热流密度 降低 两相流不稳定性大致可以分为两大类 静力学不稳定性和动力学不稳定性 静力学不稳定性是非周期性的改变系统的稳态工作运行点 它的基本特征是 系统在经受一个微小扰动后 会从原来的稳态工作点转变到另一个不相同的 稳态工作点运行 这类不稳定性是由于系统的流量和压降之间的变化 流型 转换或传热机理的变化所引起的 2 动力学不稳定性是周期性地 改变系统的稳态工作状况 这里惯性和反馈 效应是制约流动过程的主要因素 它的基本特征是当系统经受某一瞬间的扰 动时 在以声速传播的压力扰动和以流动速度传播的流量扰动之间的滞后和 反馈作用下 流动发生周期性振荡 这类不稳定性的产生主要是由于系统的 流量 密度 压降之间的延迟与反馈效应 热力学不平衡性以及流型的转换 等原因引起的 本文研究的流动不稳定性为动力学不稳定性 1 2研究现状 从上可见两相流动不稳定性的研究对于反应堆的安全运行以及事故安全 分析都具有十分重要的意义 关于这方面的 研究已经产生了很多成果 然而 大部分这方面的研究 都只是集中于单一的流体动力学不稳定性 而没有考虑 到核热耦合的相互作用 这种作用通过中子 即反应性 把热流密度与空泡份 额联系在一起 结果也就与流量联系起来了 1988 年发生在美国的 LaSalle 核 电站紧急停堆事故即是由于核热耦合产生不断的振荡产生的 LaSalle 事故被认 为是由耦合的核热水力不稳定性产生的 但它的机制还未被完全了解 2 目前关于核热耦合流动不稳定性的研究主要有时域法跟频域法两种 这两 种方法各有利弊 时域法得出的结果直接随时间而变 可以很直观地看出是否 发生了流动不稳定 但是时域法需要占用大量的计算时间 而且还容易产生由 于数值计算方法引入的计算振荡问题 频域法不产生各种各样的变量在振荡中 表现的瞬时 结果因此不会产生大量的中间计算结果 缺点是不直观 在实际的动力系统中 很多装置具有上下两个联箱间并联一组管道的结构 可以通过平均管分析得到关于整体的不稳定性 但是有一种情况却是系统总流 量保持不变 而管道之间已经发生了脉动 这是一种管道之间的相互作用 非 常隐蔽 更加危险 在近年来成为不稳定性研究中的一个热点和难点问题 一 般认为其求解的边界条件应该是总流量不变和两端压降不变 但是尚未见到实 验证明这一边界条件的合理性 在很多的模型求解中 一般仅仅将两通道之间 的压降保持相等 而这一压降是变化的 在文献 3 中提到的方法尽管保证了压 降不变 但在同一个时间步长内两个通道的压降是不等的 这一矛盾在解方程 绪论 组时无法避免 在这一点上需要更多的实验验证和理论研究 本课题的研究基 于前一种方法 即是求解时保持各个通道之间的压降相等 阮养强 4 基于多通量频域法 针对多个通道相互并联的情况 考虑相邻通 道之间以及通道与外部回路之间存在的多重反馈和耦合作用 研究了系统参数 变化以及各通道之间的工况偏差对多通道耦合系统密度波不稳定性的影响 作 者认为无工况偏差 无耦合作用的并联闭式多通道系统可以简化成单通道来分 析 耦合作用以及工况偏差的 存在对系统稳定边界有重要影响 入口节流系数的增加 相变数的减小以及雷 诺数的增加改善了两通道耦合系统的稳定性 但是入口节流系数比率的增加对 系统稳定性有非单调性的影响 这些同本研究中的一些结论可以说是不谋而合 的 Lee 和 Pan 5 对强迫循环的并联通道不稳定性进行了研究 研究针对沸水堆 参数 中高压 多通道 分析了 2 6 个通道并联的情况 同单通道相比 多通 道并联不一定会使系统稳定性增强 文中还采用了类似于 Lahey 的相空间分析 得到了并联通道的吸引子 在同一个脉动工况中 各个通道的吸引子并不相同 周云龙等 6 提出了一个简化 忽略过冷沸腾 热负荷均匀 焓值线性变化等 的并联通道密度波型不稳定性的非线性数学模型 认为系统的稳定性与初始扰 动的幅度无关 只是取决于运行条件 其判别条件是以振荡发散与否来考察不 稳定 许圣华等 7 也进行了高压系统平行通道自然循环不稳定性的研究 Munoz Cobo 8 等对沸水堆并联通道的不稳定性进行了研究 在热工方面采 用了类似于 Clausse 等的 Galerkin nodal 方法 在物理方面考虑了功率分布和燃 料动力学以及多种反馈作用 在并联通道求解过程中采用交替性的满足外部压 降作为边界条件和求解方法 但是这种方法是否合理还有待于进一步研究 黄彦平 9 等针对由七根双层套管组成的多管平行通道流动不稳定性实验段 进行了流动不稳定性实验 结果表明 实验段的管间脉动主要表现为密度波脉 动 不规则脉动和热力型脉动 3 类 与单通道流动不稳定性 两管平行通道管 间脉动 均有明显区别 4 Wu 等 10 在水力学直径 158 8 和 82 8um 的微通道中通过可视化的实验装置 观察到了周期性沸腾现象并记录下流量 温度 压力的脉动曲线 陈听宽等 11 多年来致力于高压汽水回路上的垂直并联管中汽液两相流不稳 定性的研究 通过大量实验确定了压力 质量流速 入口过冷度 热负荷以及 热负荷不对称分布 进口及出口节流和可压缩容积等对不稳定性的影响 得出 了压力降型和密度波型不稳定性的界限 给出并联管中计算不稳定性起始条件 的无因次方程 为大型直流锅炉和蒸汽发生器的设计提供依据 文中认为管间 脉动是在系统上游无可压缩容积时 热负荷达到一定程度出现的一种密度波脉 动 总流量和出口压力基本不变 其稳定性比单管要差 李会雄等 12 在其文章中报告了在高压汽水两相流实验台上对垂直上升并联 多通道中的汽水两相流密度波型不稳定性进行的系统的实验研究 发现了并联 多通道 三通道 中汽液两相流密度波型不稳定性的主要特征 确定了系统压力 质量流速 入口过冷度 热负荷 进出口节流和可压缩容积等对该类不稳定性 的影响 并与单通道内的密度波不稳定性进行了对比分析 文中认为多通道比 单通道更加不稳定 这一点同陈听宽等的结论是一致的 更有趣的结论是通道 数目为奇数的系统稳定性要好于通道数目为偶数的系统 这一点可能还需要更 多的实验数据来进行支持 李虹波 13 等对平行双通道系统的管间脉动实验进行了报告 其特色在于采 用了矩形流道的实验段 目前关于核热耦合流动不稳定性的研究主要有时域法跟频域法两种 这两种方 法各有利弊 时域法得出的结果直接随时间而变 可以很直观地看出是否发生 了流动不稳定 但是时域法需要占用大量的计算时间 而且还容易产生由于数 值计算方法引入的计算振荡问题 频域法不产生各种各样的变量在振荡中表现 的瞬时 结果因此不会产生大量的中间计算结果 缺点是不直观 1 3研究内容 在广泛调研并联多通道以及核热耦合资料的基础上 本文从时域法研究 在西安交通大学核科学与技术学院开发的 船用核动力装置多通道流动不稳定 绪论 性的分析程序 的基础上 构建具有六组缓发中子的点堆功率模型 嵌入核热 耦合模块 由于原程序的算法是单精度类型 嵌入核热耦合模块后 程序发生 了振荡 计算速度很慢 然后把程序数据类型都改成了双精度 但是由于不断 的耦合迭代 计算还是很慢 不断改变系统的压力 进口过冷度以及进出口阻 力系数并把结果与原程序运行结果对比后 分析研究 得出以下结论 1增大通道进口流体的过冷度 有利于多通道系统的稳定 2增大进口阻力和减小出口阻力有利于系统的稳定 3增加系统的压力 可提高系统的稳定性 4核热耦合的相互作用 反应堆的负反馈会严重系统的稳定性 西安交通大学本科毕业设计 论文 6 2点堆模型和核热耦合模型 根据核反应堆物理的分析结果 具有六组缓发中子的点堆方程可以如下表示 14 2 1 1 eff ii i ieff i ii k dn tn t n tc tq dtll k dc t n tc t dtl 式中 中子密度 n t 第 组缓发中子先驱核浓度 i c ti 中子寿命 s l 反应堆有效增殖因子 eff k 核反应堆外中子源每秒产生的中子数 q 第 组缓发中子先驱核衰变常数 i i 1 s 第 组缓发中子份额 i i 总缓发中子份额 i 在推导该方程时并没有考虑核反应堆中存在的各种反馈 即认为反应堆的 功率水平比较低 反应性稳定系数可以被忽略 由于中子代时间为 反 eff l k 应性为 并且外中子源 故方程 4 1 可以简化成 1 eff eff k k 0q 2 2 ii i ii ii dn t n tc t dt dc t n tc t dt 点堆模型和核热耦合模型 注意这里的 16 i 反应堆的反应性相对于反应堆的某一个参数的变化率称为该参数的反应性 系数 如反应性温度系数 空泡系数 功率系数等等 对于我们研究的这个堆 属于沸水堆 因此真正能对反应性起到反馈作用的是反应性温度系数和空泡系 数 应该指出 反应堆内的温度是随空间变化的 堆芯中各种成分 燃料 慢化 剂等等 的温度以及温度系数都是不同的 因此反应堆中的总的温度系数等于各 成分的温度系数之和 即 2 3 Ti i i T 在实际的反应堆中 一般温度系数只需考虑燃料温度系数和慢化剂温度 F T 系数即可 M T 由于多普勒效应 以低富集铀为燃料的反应堆中 燃料温度系数一般都是 负的 而当温度增加时 慢化剂温度系数却可能出现正值 尤其在寿期初当慢 化剂中硼的浓度比较大时更可能出现这种情况 空泡系数是指在反应堆中 冷却剂的空泡份额变化百分之一所引起的反应 性变化 以表示 即 M V 2 4 M V x 式中 空泡份额 x 一般来说 对于轻水堆来说 当出现空泡或者空泡份额增加时 是负效应 值得注意的是在本文中 由于直接求出空泡份额比较困难 因此用功率相对值 的变化来表征空泡份额的变化 因此空泡系数也用功率相对值的变化来表征 上述三种系数在几种典型反应堆中的反应性系数如表 2 1 所示 15 西安交通大学本科毕业设计 论文 8 表 2 1 几种堆型的反应性系数 反应性系数沸水堆压水堆重水堆 燃料温度系数 5 10 K 4 1 4 1 2 1 慢化剂温度系数 5 10 K 50 8 50 8 7 3 空泡系数 200100 00 在本文中 由于引入的扰动比较小 为了研究的方便 把上述三种系数都 设为常数 即取 5 2 5 10 F T 4 2 9 10 M T 4 5 0 10 M V 并联多通道理论模型 3并联多通道理论模型 对于广泛使用的压水堆动力装置 主要包括反应堆 一回路系统 二回路 系统三部分组成 一回路系统 又称反应堆冷却剂系统 其主要功能是在正常 运行时将堆芯产生的热量传给蒸汽发生器 使二回路工质变为蒸汽 二回路系 统再将产生的蒸汽转换为电能或者其他形式的能量 压水堆动力装置的原理流 程如图 3 1 所示 图 3 1 一体化反应堆流程简图 本研究所关心的即是这一系统中的具有并联多通道的堆芯装置 可将其简 化为由上下两个联箱和中间多个相连通道所组成的系统 如图 3 2 所示 针对该系统建立并联多通道分析模型 该模型建立基于 Clausse 和 Lahey 的工作 在其基础上加入了入口段和上升段的模型 加热段总长度为 划分 H L 为单相区和两相区 单相区长度为 为单相区节点数 入口段长度为 N L s N E L 上升段长度为 为了更好地分析问题 先做如下假设 16 R L 西安交通大学本科毕业设计 论文 10 1 在给定系统压力下单相流体物性为常数 2 均相流 3 轴向均匀加热 4 轴向均匀加热 5 入口流体过冷 6 两相流体处于热力学平衡态 7 通道分为入口段 加热段和上升段 Upper Plenum Lower Plenum heater section riser section entrance section throttle valve LR LH LE single phase zone two phase zone 图 3 2 多通道并联系统 3 1单相水热工水力微分方程 质量守恒方程 0 ll l u zt 3 1 式中 l 单相水的密度 t 时间 l u 速度 并联多通道理论模型 z 竖直方向坐标 动量守恒方程 0 2 22 z P gu D f u zt u lll e l ll ll 3 2 式中 l f 单相水的摩擦系数 e D 通道当量直径 能量守恒方程 A q hu zt h L lll ll 3 3 3 2两相水热工水力微分方程 基于均相模型的两相流体控制方程 质量守恒方程 0 HH H u zt 3 4 动量守恒方程 0 2 22 z P gu D f u zt u HHH e TP HH HH 3 5 能量守恒方程 A q hu zt u L HHH HH 3 6 上式中 下标H均表示位于加热段的两相水的均值 可见对于均相流模型 两相区和单相区的控制方程均可以写成如下形式 质量守恒方程 0 u zt 3 7 能量守恒方程 A pq hu z h t H 3 8 西安交通大学本科毕业设计 论文 12 动量守恒方程 z P g u zkzk D f u z u t ei H 2 2 2 3 9 式中 h 水焓值 q 热流密度 H P 加热周长 i k e k 水进 出口阻力系数 A 通道面积 密度表达式 f f hh 3 10 1 f fg fg f hh h v v f hh 3 11 对于加热段单相区 根据假设 节点间焓线性变化 可得到焓在单相区各节点 的表达式为 if s in hh N n hh 3 12 在加热段单相区对能量方程 3 8 结合式 3 10 3 12 在两个节点之间积分可得 到单相区各节点边界的变化 dt dL LL hhA Pq Nu dt dL n nn iff H si n1 1 22 3 13 上式中当 s Nn 时即是沸腾边界的表达式 通道加热段出口速度可以由下式求得 NHie LLuu 3 14 其中 fg fg H h v A p q 3 15 3 15 中 fg 表示饱和蒸汽和液态水的比体积之差 fg h 饱和蒸汽和液态水 的比焓之差 并联多通道理论模型 N L是加热段单相区长度 也就是随时间变化的沸腾边界 对于每一个通道其质量变化可以写成 jjejejif jH Auu dt dM 3 16 1 jM 对质量守恒方程式 3 7 积分可得两相区质量为 1 ln 2 ef ef fNH LLAM 3 17 故加热段总质量 2 MALM NfH 3 18 上式中考虑到在单相区水的密度变化很小 为简化积分 故在单相区采用饱和 水密度代替过冷水密度 dt d je 可利用式 3 17 3 18 联立求解得到 Mj LL dt dL uu dt d f e f e NHf fe N ef ef feeif je 2 1 ln1 1 ln 1 2 2 3 19 对于上升段 其流体质量变化为 d dM 1 r rreRu A t 3 20 上升段两节点之间流体质量可以通过下式计算 1 1 1 1 1 ln M rr R RR rr rr R RR r N LA N LA 3 21 3 3压降模型 3 3 1入口段 根据假设 1 入口段内的流体质量是不变的 并且入口段流体流速等于加 热段入口的流体流速 故入口段的压降可以通过下式计算得到 西安交通大学本科毕业设计 论文 14 2 E E1E E P 2 fi f u L fgL D 3 22 式中 E L 入口段长度 1 f 根据表 3 1 17 计算 表 3 1 摩擦系数关系式 区域关系式 层流 Re 1000Darcy 过渡区 2000Re1000 0 048单相区 紊流 2000Re Balasius 1 f 在1400Re800 范围内采用插值计算 25 0 E2 tp GDcf 3 23 3 23 中 c 根据单相流动摩擦计算式的系数选定 1 1 lvtp xx 3 24 式中 tp 两相平均粘性系数 l v 液态水和蒸汽的粘性系数 3 3 2 加热段 加热段总压降可以写成下面四项压降和的形式 afgI PPPPP H 3 25 具体计算式如下 加速压降 22 0 2 H ife L ea uudz z u p 3 26 摩擦压降 包含局部压降 并联多通道理论模型 f P 222 0 2 1 2 1 s s iii L L H L H ukdzu D f dzu D fH N N 2 2 1 eee uk 3 27 s 21 2 1 2 Nif H Lu D f efi fe NH H u LL D f ln 1 12 2 2 s AML uuu chHf ef iei 1 2 2 3 1 s 2 2 NH eff ef ie LL uu 1 ln s ef efNH LL 2 2 1 ifi uk 2 2 1 eee uk 重位压降 A gM gdzp L g H 0 H 3 28 惯性压降 1 H H 00 HH ef chfie i LL I AMLuu u A M dt d dzu dt d dzu t p 3 29 3 3 3 上升段 通过对上升段长度积分可得到总压降 22 RRR R RRR 22 LLL P eeNeRe e eeN RR R eRRRRRe L L R uug dt d u dt du u AL M uAgMAMu dt d dz z P R R R 3 30 3 4 无量纲化 令 H LLL f Hs uLtt 西安交通大学本科毕业设计 论文 16 s uuu H Lzz ssxfof uAQhhh 2 s uPP f 2 sH gLuFr jinjs uu 其中下标s表示标准值 j 表示第 j 个通道 根据上面的定义 可以将守恒方程组写成下面的形式 0 u zt 3 31 z P Eu F zzkuu zt h r N i ii 1 2 2 1 2 2 3 32 1 uh zt h 3 33 式中 1 当0 h 3 34 1 1 hNpch 当0 h 3 35 1 单相区 2 两相区 单相区节点微分方程 dt dL LL N N Nu dt dL jn jnjn jsub jpch sji jn 1 1 22 3 36 通道加热段流体质量 jejeji j uu dt dM H 3 37 加热段出口的密度方程 ln 1 1 1 1 ln 1 2 jejej je jijeje j je jejeje uu dt d dt d 3 38 并联多通道理论模型 各类压降计算式 2 E 2 j 1jE P js j ij u gL u 3 39 1 1 1 1 H Hj je jjjsub jijI MN uM dt d p 3 40 jr j jg F M p H 3 41 2 1 2 1 3 1 1 1 1 1 1 1 1 2 1 ln 1 1 1 2 2 H 2 H 2 2 2 2 1 jejejejijijj j je je jjpch je jjjpchji je ji je j jijjf ukukM N MNu uuP 3 42 2 jijejeja uup 3 43 2 j j j j jH j 2 j j Hj j j j jH j j jR P eeN R R e R R RRRRe u A AM u L L LFrMAAMu dt d R 3 44 其中 Ad dM Re R e H R u A t 3 45 dt d dt du dt due i 3 46 3 5 并联通道模型 从图 3 2 中可以看出 由于通道均和上下联箱相连 所以各个通道压降相 等 因此 M PPP 21 3 47 将各个压降计算式带入上式并通过化简写成下面得格式 MjB dt du A dt du j i j ji 3 2 1 3 48 西安交通大学本科毕业设计 论文 18 式中 jj MMA H1 H 3 49 2 1 3 1 1 1 1 1 1 1 1 2 1 ln 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 2 1 3 1 1 1 1 1 1 1 1 2 1 ln 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 j j j j jH j 2 j j Hj j j j jH j j H 2 H 2 2 2 1 2 2 j 1 2 2 H H H 2 H H 1 1 1 1 H 1 1 2 1 1 H 1 1 1 1H 1 1 1 11 H 1 1 2 1 11 1 1 H11 1 1 2 1 1 1 2 2 1 11 2 1 1 2 1 1 1 1 H 2 1 1 2 1 1 1 1 1 11 1 1 H 1 H 2 1 11 1 1 H 1 1 1 H eeN R R e R R RRRRe jj j jeje jjpch je jjjpchji ji je j jij js je ij ji ji je je jeje j je jjpch ji j j je jjpchje j je jpchj j eeN R R e R R RRRRe ee pch e pchi e i e i s e i i i e e ee e pch i e pche e pch j j u A AM u L L LFrMAAMu dt d M NMNu j uu u gL u u k uk Fr MN u dt dM M N dt d M N dt d M u A AM u L L LFrMAAMu dt d M NMNu uu u gL u Fr M u k uk N u dt dM M N dt d M N dt d M B R R 3 50 2 1 3 1 1 1 1 1 1 1 1 2 1 ln 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 2 1 3 1 1 1 1 1 1 1 1 2 1 ln 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1 1 j j j j jH j 2 j j Hj j j j jH j j H 2 H 2 2 2 1 2 2 j 1 2 2 H H H 2 H H 1 1 1 1 H 1 1 2 1 1 H 1 1 1 1H 1 1 1 11 H 1 1 2 1 11 1 1 H11 1 1 2 1 1 1 2 2 1 11 2 1 1 2 1 1 1 1 H 2 1 1 2 1 1 1 1 1 11 1 1 H 1 H 2 1 11 1 1 H 1 1 1 H eeN R R e R R RRRRe jj j jeje jjpch je jjjpchji ji je j jij js je ij ji ji je je jeje j je jjpch ji j j je jjpchje j je jpchj j eeN R R e R R RRRRe ee pch e pchi e i e i s e i i i e e ee e pch i e pche e pch j j u A AM u L L LFrMAAMu dt d M NMNu j uu u gL u u k uk Fr MN u dt dM M N dt d M N dt d M u A AM u L L LFrMAAMu dt d M NMNu uu u gL u Fr M u k uk N u dt dM M N dt d M N dt d M B R R 各个通道的流量变化和等于总流量的变化 根据这个关系可以得到下面的式子 1 2 H 2 H 1 M j jj M j jj tot i AA BA dt dW dt du 3 51 数值方法 4数值方法 整个核动力系统的动态仿真 在数学上最终可归结为求解如下形式的以时间 t 为 基本变量参数的变系数非线性常微分方程组初值问题 4 1 yytf dt yd 4 2 00 yty 微分方程组所表示的系统物理特性可通过其雅可比矩阵的特征值表 f J y t 示 特征值的实部 j e R 对应于振幅的增减 虚部 j e I 对应于振动的角频率 当特征值 的实部 mi ss RR 均 0 时 系统是稳定的 及任何初始扰动都随时间 t 的增长而衰退 如果令 m jR j e j 1 1 j 为时间常数 若系统是由若干个 成分 组合而成 不同的 成分 可以有不同的时间常数 其中最大的时间常数 j max max 则表达了全过程的活跃时间 而最小的时间常数 j min min 则表达了系统最 敏感 环节的反应速度 即过渡时间 如果最大时间 常数与最小时间常数相差悬殊 即 max min 1 这样的系统或微分方程组便为刚性 的或称病态的 4 1 系统的刚性 对于一般线性常系数系统 4 3 tyAy 若系统 4 3 为刚性 则矩阵满足下列条件 A 1 系数矩阵 A 的所有特征值的实部不是很大的正值 2 系数矩阵 A 至少有一个特征值的实部模为很大的负值 3 对应于具有大负值实部的特征值的解分量的变化是缓慢的 西安交通大学本科毕业设计 论文 20 对于非线性系统 可以用 Jacobi 矩阵的特征值 ytfy 0 Tt y f 进行分析 2 1 mi i 下面是我们常用来分析热工水力问题的方程组 xtyfy 4 4 1 11 1 n dfal ii vff w wp nnnn nnnn nn n LdW PPPP Adt qqdT dtC W HWH AAHq U tLA 刚性是系统本身的性质 也就是说如果一个系统是刚性的 那么描述这一系统的 方程组也是刚性的 其物理意义就是 求解慢变化的解 但是系统中存在快速衰减的 扰动 这种扰动使得慢变化解的计算复杂化 当数值积分这样一个系统时 一旦快变 分量消失时 希望选取的步长 h 适合于积分慢变分量 考虑求解方程 4 5 0 0 0 Ttyy tFtFtyty 其中 F t 是 t 的慢变函数 精确解 y t 为0 0 0 FyetFty t 或者 tFtyehtFhtY t 其中 h 为步长 因为 显然在一个非常短的求解时间后 即暂态阶段后 0 刚性解分量就不在解中出现了 这就意味着慢变化函数 F t 在大部分积 0 0 Fye t 分区间 0 T 上支配着要计算的解 解的第二个表达式说明在任何时刻 t 对慢变解 ty F t 的扰动将迅速衰减 微分方程 4 5 本身的特点是 当时 对各种初值 问题都是不稳定的 0 0 y 即系统时不稳定的 若 即的值非常小的时候 解的曲线几乎是平行线 如0 果初始值有误差 那么到终点几乎还是那个误差 这类问题传统的 RK 方法和线性多 数值方法 步法都能够很好的解决 若即具有绝对值大的负值 解的曲线迅速收敛 不0 管是什么值 经过一个暂态阶段后 解的曲线迅速称为特解 这种问题使用传 0 y F t 统的数值积分方法却会遇到困难 考虑显式 Euler 方法 1nnn hyyy 应用于式 4 5 得到 1 1nnnnn thFtFtFyhy 其中的可以看成是在时刻对解 F t 的一个扰动 当时 这 nn tFy n t11 h 个扰动将被放大 计算将出现数值不稳定 只有当 即11 h02 h 这个扰动才衰减 当时 即问题是刚性的 显示 Euler 法要求整个积分区0 间都必须采用小步长 许多传统的显示方法一般都按照最小时间常数来选取步长 即 使这些时间常数所对应的解是可以忽略的 这是因为多数显示方法求解方程 时 仅当充分小时才是稳定的 故对刚性系统进行数值积分时 ytfy yfh 为了保证数值方法的稳定性 显式方法使用的步长常常与所要求精度无关 如果考虑隐式 Euler 方法 11 nnn hyyy 同样代入式 4 5 得到 4 6 1 1 11 11 1 nnnnnn tFhthFtFhtFyhy 当时 这个数值积分方法可以迅速的使扰动衰减 由于假设0 h nn tFy F t 是光滑缓慢变化的函数 所以有可以用近似 于是有 1 nn thFtF 1 n tF 因而对于任意的 当时 1 111 1 nnnn tFtFhthFtFh 0 h 这时的数值积分过程可以根据计算精度选取步长 用绝对稳定域有限 11 nn tFy 的数值方法求解刚性问题的数值解 由于对步长的限制 导致计算步数很大 为减h 小计算步数 应对步长无限制 h 下面我们来看对于刚性方程组 4 6 的求解 4 7 1 0 404040 0 0 202119 1 0 201921 33213 23212 13211 yyyyy yyyyy yyyyy 西安交通大学本科毕业设计 论文 22 其真实解如下 4 8 40sin40 cos 40sin40 cos 2 1 2 1 40sin40 cos 2 1 2 1 40 3 402 2 402 1 ttety tteety tteety t tt tt 从图 4 1 可以看出 其他参数不变 提高精度以后 Gear 算法的计算间距减小 以适应更高精度的要求 图 4 2 是反应在 Gear 算法可以很好适应这样的刚性方程 既 能反应解向量中快速变化的部分 也可以改变步长适应缓变量 以提高计算速度 051015202530 0 1 0 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 计算间距 输出点数 计算精度 0 01 0 001 0 0001 图 4 1 计算精度对解的影响 0 020 000 020 040 060 080 100 120 140 160 18 1 0 0 5 0 0 0 5 1 0 计算值 计算距离 B C D 图 4 2 计算中瞬变量和缓变量对步长的不同要求 数值方法 4 2Gear 方法 Gear 方法是满足刚性稳定性的多步法 十分适合与刚性系统的仿真 Gear 法的 一般形式为 18 4 9 knk k i inikn fhyay 1 0 其中和是待定的个常数 其系数见表 7 1 其中的用 1 2 1 kiai k 1 kC 于计算阶段误差 Gear 一般是非线性隐式方程 要采用迭代法求解 表 4 1 1 6 阶 Gear 方法公式 4 2 1 简单迭代法 求解隐式方程通常采用预估 校正法 一般可以将这种预估公式取为 4 10 1 0 1 1 k i knkinikn fhyay 选取其中的系数和使得式 4 10 的精度阶也是 基本公式是用式 4 10 i a 1 k pk 求出的预估值 再代入 4 9 求出满足非线性方程组的解 kn y 0 kn y kn y 4 11 1 0 knknknk k i inikn ytfhyay 如果求解的问题是非刚性系统 可以使用简单迭代法求解式 4 11 4 12 1 0 1m knknknk k i ini m kn ytfhyay 其中 由式 4 10 给出 为保证收敛 要求 0 kn y 4 13 1 yfh k 如果求解问题是刚性系统 Jacobi 矩阵的模很大 因此 积分步长 h 除了yf 受到积分方法的计算稳定性的限制以外 还要受到不等式 4 13 的制约 这个不等式 西安交通大学本科毕业设计 论文 24 表明 步长 h 要限制到问题最小时间常数的数量级 这与用绝对稳定区域为有限的方 法求解刚性问题给的对步长的限制条件是一致的 因此 简单迭代法是不适用的 为 了克服这个缺陷 一般采用 Newton 迭代方法求解非线性方程组 4 11 4 2 2 Newton 迭代方法 方程 4 11 可以简写为 4 14 gytfhy knknkkn 其中 在当前的时刻中 这是一个常数 所以 有如下的基本迭代 1 0 k i ini yag 1 kn t 式 4 15 gytfhy m knknknk m kn 11 Newton 迭代法的基本思想是 把的非线性函数在的附近线性化 取其近似的 kn y m kn y 表达式 11m kn m kn m knkn m knkn m knkn yyyt y f ytfytf 代入 4 15 得到基本的 Newton 迭代格式 4 16 1 1 gyffhy y f hIyy m knknk m knk m kn m kn 将式 4 10 和 4 16 联合 可以得到 Newton 迭代法的预估 校正公式 4 2 3 Gear 方法 为了变阶 变步长方便 Gear 将方法表示成 Nordsieck 的向量形式 Nordsieck 提出存 储 的近似值来代替存储前几步的函数及其导数的近似值 目 n y n hy 2 n y h ty 的是使改变步长时所需的附加计算简单一些 将初值问题 4 17 0 0 yy ytfy 的解及其各阶导数展成 Taylor 级数 并乘以因子 则有 ty i hi 3 2 1 32 ty h ty h ty h tyhty 数值方法 4 18 3 3 3 2 2 1 1 2 ty h ty h ty h hty h 4 4 3 3 2 2 432 ty h ty h ty h hty h 记 4 19 Tk n k nnnn y k h y h y h yZ 2 1 2 其中 是的近似值 如果公式中所用的前几步导数 nnn yyy nnn tytyty 值较多 那么表达式 4 19 中就有更高阶的导数的近似值 由式 4 19 可以由 1 n Z 来计算 n Z 4 20 nnZPZ 1 其中是 Pascal 上三角矩阵 它的第 元素为 Pij 当 时 0 ij iij j j i 1 1 1 1 31 1 321 11 1111 k k kk P 当时为 0 ij 由于式 4 20 没有用到微分方程 按数值计算的观点来看 它是不稳定的 为了克服 这一缺陷 加一个校正公式 而式 4 18 只做为零次近似 Gear 给出了计算刚性问题 的公式 预测 1 0 knknZPZ 校正 4 21 1m kn m kn m knZFuLZZ Mm 2 1 0 终值 M knnZZ 其中为迭代矩阵 1 01 Y F hlIlL Z F u 西安交通大学本科毕业设计 论文 26 其中 依赖与步长 h 和 Jacobi 矩阵 同时通过也依赖于方法1 10 ll k uyf k 的阶数 如果变化慢 那么对于一步或者其中步长和阶不变的若干步 在迭代yf 过程中 矩阵无需重新计算 这样可以大大减少计算量 u 从 4 21 中可以看到 公式仅仅用到前一步的及其导数的信息 而没y nnn yyy 有用到多个点上的的信息 因此 式 4 21 是一种单步多值算法 利用这种算法可y 以自启动 而且可以克服多步法在实现变阶 变步长方面的困难 表 4 2 给出了不同 k 对于向量的分量值 L 表 4 2 k 对于向量的分量值L 阶数123456 0 l1 3 2 11 6 50 24 274 120 1764 720 1 l111111 2 l 3 1 11 6 50 35 274 225 1764 1624 3 l 11 1 50 10 274 85 1764 735 4 l 50 1 274 15 1764 175 数值方法 5 l 274 1 1764 21 6 l 1764 1 其他的多步法亦可用上述方法得出型如 4 21 式的预估校正公式 由于使用了 Nordsieck 向量 矩阵与 4 21 的预估式已与算法无关 即不论 Gear 还是 Adams 法 P 都可以用 4 21 的预估式作为预报公式 而且矩阵都是一样的 各种算法的差别仅P 仅在于向量的不同 L 4 2 4 Gear 算法的阶与步长的控制 由于刚性问题的解存在快变分量和慢变分量 在数值求解过程中 当解迅速变化 时 需要采用较小的步长 而当解趋于稳定 就应该采用大步长 此外 在考虑方法 的自启动时 开始的时候应该采用低阶公式和小步长进行计算 随后 在计算中 阶 数和步长要不时的做调整 以获得最优的步长和阶数 所谓的最优是在满足计算精度 的同时计算量最小 设求解时间区域 为整个区间的最大误差要求 为单位时间的最大容许 ba bae 误差 二者的关系如下 ab e 所以每一步的最大容许误差为 这个误差可以与局部截断误差联系起来 K 阶eh Gear 方法的局部截断误差有形式 4 22 2 1 1 1 kkk kk hOtyhCE 其中是依赖于方法的一个常数 表 2 给出了 1 6 阶 Gear 算法的值 如果忽 1 k C 1 k C 略 4 22 中的高阶小量 用代替则有 eh k E 4 23 kk k htyCe 1 1 根据这一表达式 利用 Nordsieck 向量实现程序的变步长和变阶比较容易 因为利 n Z 用 Nordsieck 向量容易进行误差估计和公式起步计算 下面是 Gear 方法的处理 n Z 对于 k 阶方法 局部截断误差为 4 24 2 1 1 1 kk n k kn hOtyhCE 西安交通大学本科毕业设计 论文 28 略去 4 24 式中的高阶小量 则第 n 步的截断误差为 4 25 1 1 1 tyhCE k n k kn 为了估计 先要估计 已知 n E 1 k n y Tk nnn Tk n k nnnn zzzy k h y h y h yZ 2 1 2 现在用差分来估计 即 k n y 1 k n y h yy y k n k nk n 1 1 对上式两端同乘可得 1 k h k 4 26 1 1 1 1 1 k n k n k n kk n k n k k n zzz k h h yy k h y 将 4 26 代入 4 25 式得 4 27 1 k nkn zkCE 由上式可知 要将第步计算所得的与前一步的这两个向量的最后一个分量相n n Z 1 n Z 减即可得到这一步的截断误差 因此 若计算过程中要求每一步的截断误差小于预先指定的误差量 必须选取步长 0 使得不等式成立 h 4 28 0 1 k nk zkC 在微分方程组的情况下 希望控制每一个分量的误差 对于选定的权向量 其分量W 是加在第 个变量的截断误差上的权分量 目的是为了对不同分量采用不同的误 i wi i y 差控制 一般可以取积分到这一步的已出现过的的绝对值的最大者 要求选取步 i y 长 h 使得不等式成立 其中方程组的每一个元都有一个和权分量 为范 k n z 2 L 数 4 29 0 2 1 W z kC k n k 程序中基本步长的控制是积分一步并且检验式 4 29 是否成立 如果 4 29 成立 则接 受这一步 否就抛弃这一步 对于下一步或者重新抛弃的步 所用的步长为 ah 其中 a 由等式 数值方法 来确定 0 2 1 1 W z akC k nk k 如果采用这个步长 并且误差又正好与成比例 即不变 则下一此检验式 1 k h k n z 4 27 正好满足 但是不总是不变 所以为了保证不等式 4 29 能够成立 常采 k n z 用稍微小一点的步长 在程序中 a 由来确定 4 30 k k n k k Wz kC a 1 2 1 0 1 2 1 1 为了变阶 还必须检查在其他阶的公式中所用的步长 对于降阶 有 4 31 k k n k k Wz kC a 1 2 1 0 1 1 3 1 1 其中 是向量中的最后一个分量 k n k k n y k h z nZ 对于升阶 需用二次差分来近似 即 2 k n y h yy y k n k nk n 1 1 1 2 由于 1 1 k k n k n h k zy 所以有 2 2 1 1 2 k n k k n k n k k n z h k h zz h k y 从而可得 4 32 2 1 2 2 2 0 1 1 4 1 1 k k n k k Wz kC a 在式 4 28 和 4 29 中因子 1 3 和 1 4 的选取是考虑在变阶变步长的情况下 误差的估 计更加间接 而且变阶还要增加额外的计算 所以系数取的更加严格 同时希望需要 变阶时 尽量利用降阶 因为降阶的计算量稍少一些 变阶变步长以后 Nordsieck 的向量的变化如下 nZ 西安交通大学本科毕业设计 论文 30 1 不变阶只变步长时 有 4 33 n k k k k T k n kk k nkn n Z a a a y k ha yhayZ 1 2 2 降阶变步长时 有 4 34 T k n kk k nkn ny k ha yhayZ 1 1 11 1 1 3 升阶变步长时 有 4 35 T k n kk k nkn ny k ha yhayZ 1 1 11 1 1 其中 要用向后差分获得 1 k n y 4 36 h yy y k n k nk n 1 1 在积分过程中 每一步计算 选取其中最大的 以便确定下一步要采用 11 kkk aaaa 的合理的步长和阶 但是 在程序中每一步都改变步长和阶不一定有利于加速计算 所以 1如果这一步失败 则重新估计 a 2在上一次变阶或者步长以后的 k 1 步估计 a 如果上次估计时步长不放大 则计算 10 步以后估计 aa 程序计算结果及分析 5程序计算结果及分析 本研究采用较为直观的时域法对船用核动力装置并联多通道流动不稳定性进行了 分析 根据上文提出的船用核动力装置并联多通道流动