中国农业大学数值分析

2024-05-14

中国农业大学数值分析（精选8篇）

篇1：中国农业大学数值分析

2017中国农业大学数值分析考点总结

1、考试时间：1月14日上午9点（2.5-3个小时）

2、考试地点：一教401或者一教圆厅，考前一到两天会在邮箱发布每个人所在考场

3、考试工具：计算器、几张干净的草稿纸、证件

4、答疑时间：1月12日上午9点-下午2点；答疑地点：理学院407

一、题型：（1）填空题6-7个，每题两分，共计12-14分（2）计算题、证明题、公式推导题、简述（可能考比较）

二、具体内容第一章

1.2（1）五个概念：绝对误差（限）、相对误差（限）、有效数字；

（2）数值运算的误差估计：估计函数的误差（一元、多元）

1.3（1）算法稳定性概念；（2）设计算法五个原则注意：本章习题很重要第二章：分值 15分

（1）牢记拉格朗日、牛顿公式及余项公式（2）差商的计算（牛顿）（3）关注差商的性质3（4）埃尔米特差值，用两种方法构造（基函数、承袭性）（5）掌握为什么要分段低次插值；计算分段线性插值、分段三次插值

（6）简答题：简述样条插值过程（没看PPT2-6）第三章（1）最佳平方逼近和最小二乘法会做题第四章分值15分

（1）机械求积公式的定义；代数精度的求法（2）插值型求积公式定义和判定方法（2条）（3）会求求积公式的余项

（4）牢记梯形公式、辛普森公式、柯特思3个公式，会推导余项

（5）会计算复合求积公式

（6）高斯求积公式的定义；推导一个点两个点的高斯-勒让德公式，会计算（不）。

（7）定理证明：

1、插值型求积公式充要定理；

2、高斯点充要条件定理（不）第五章

（1）会用五种办法解方程组：高斯消去法、列主元素法、三角分解法、平方根法（包括改进的平方根法）、追赶法（2）会计算矩阵向量常用范数，2分，填空题（3）会计算矩阵条件数，判断病态

第六章（1）三种迭代法的构造和矩阵形式推导；（2）判断是否收敛（无定理）（PPT6-3 27疑问）第七章（1）不动点迭代法收敛判断：共四个定理（2）判断收敛的速度（3）牛顿公式推导（4）定理：牛顿局部收敛2阶证明第九章

（1）会推导欧拉五个公式（用差商代替求导来推导，用数值积分办法推导）

（2）会求局部截断误差（所有、特别作业7、8、11、12 例题7、8）

（3）一阶微分方程租和高阶微分方程会计算

总结：

1、3+3+5共11个公式（1）牢记拉格朗日、牛顿公式及余项公式（2）牢记梯形公式、辛普森公式、柯特思3个公式（3）欧拉五个公式

2、定理：第二章余项；第四章两个冲要；第七章一个收敛（牛顿）定理

最后一个题目：三选一

（1）总结本学期学过的加速方法1松弛法：待定系数法、松弛法（解方程组的）答三个满分写出一般规律；2基于误差分析埃特金、史蒂芬逊龙贝格（忘记记下来的一种）

（2）讨论两种非线性方程组转化为一元二次方程（3）简述样条插值的构造思路

篇2：中国农业大学数值分析

第一章误差有效数字算法设计若干准则

第二章拉格朗日插值牛顿插值埃尔米特插值分段插值三次样条插值插值余项插值基函数

第三章函数空间正交多项式最佳平方逼近最佳一致逼近（用切比雪夫）曲线拟合的最小二乘法

第四章代数精度牛顿-柯特斯公式复合求积龙贝格算法高斯求积数值微分

第五章高斯列主元消元 LU分解敏度分析矩阵条件数第六章雅可比迭代 G-S迭代 SOR迭代收敛定理

第七章二分法不动点迭代收敛定理收敛阶牛顿法弦截法

非线性方程组牛顿迭代法第八章规范化幂法反幂法

篇3：草莓冻结过程的数值分析

新鲜草莓是一种表面较易粘结的物料,采用普通固定床风冷速冻不利于其热质的均匀传递,影响速冻效率和品质。传统振动流化速冻对于高湿或有结团倾向的物料,虽然可以通过振动的作用使物料处于良好的流化状态, 但是柔嫩多汁的草莓果形易受到损坏。本文采用波形振动半流化型速冻机来进行草莓的速冻特性和速冻工艺的研究。

食品冻结时间和冻结速率是设计与评价冻结装置性能优劣的重要指标[1]。由于冻结方式与冻结装置结构的不同,会造成物料冻结时边界条件的变化,对于众多的冻结方式与冻结装置,考虑现有的开发产品,选择波形振动半流态化型冻结装置(如图1所示)进行分析。

流化床为长方体,形状为片状、条状或球状的物料紧密堆积在一起。冻结过程的特点及相变过程为固—固相变,即在冻结过程中始终保持固态,形状不变,为数值计算带来了很大方便。模型中物料的几何尺度相对网带宽度很小(L>>d),流化床内物料紧密堆积,各物料之间形成的空隙基本是一样的,对流场的影响是均匀的[2]。在垂直于流动的方向上,各点风速差别不大,所以在冷空气中物料边界一致,冷空气通过物料层温升为1～3℃。因此,可以用平均温度值来代替表面上各点的温度,冷空气温度沿网带长度方向近似不变。

1 草莓冻结过程传热模型

为了便于分析,可以将冻结过程简化,并做如下假设:

1) 物料内部的温度分布只与定性尺寸有关,物料之间的相互影响不作考虑。

2) 流体区域的传热为一维对流换热问题;物料与冷空气之间的对流换热系数在区域内恒定;在物料内部,热交换只考虑热传导,不考虑对流换热。

3) 冻结在温度T1*～T2*范围内进行,即存在相变区,在冻结过程中,不断后移的相变界面温度不变。

4) 物料的初始温度为开始进速冻机时物料的温度,冻结体初始温度为T0。

5) 相变冻结阶段的传热过程是热量控制过程,冻结速度由传导热量的速度决定。

6) 由于物料处于流化状态,所以在冻结过程中食品的热中心就是食品的几何中心。所谓热中心,是指降温过程中食品内部温度最高的点[3]。

7) 草莓为类球形,可视为球体的传热,在相变界面上进行的是一维热流传递,垂直于相变界面,物料不可压缩,且冻结过程中与外界无质交换。

以上假设使草莓流化床冻结过程的数学模型大为简化,同时为采用一维模型提供了充分的保证。

1.1 球状草莓传热模型

球状草莓由于几何尺度小,形状类球体,故采用一维球体导热微分方程。

对球状草莓颗粒的相变导热焓方程式[4,5]为

$ρ_{p} \frac{\partial Η_{p}}{\partial t} = λ_{p} \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial Τ_{p}}{\partial r})$ (1)

式中 ρp—草莓颗粒的密度 (kg/m3);

Hp—草莓颗粒的焓值(kJ);

λp—草莓颗粒的导热系数(kJ/m·℃);

r—草莓颗粒的半径(m);

T—草莓颗粒的温度(℃)。

用控制容积法离散得到

$\begin{array}{l} r_{p}^{2} ρ_{p} (Η_{Ρ} - Η_{p}^{0}) / △ t = λ \frac{r_{e}^{2} (\frac{\partial Τ}{\partial r})_{e} - r_{w}^{2} (\frac{\partial Τ}{\partial r})_{w}}{2 △ r} = \\ \frac{λ}{2} \frac{(\frac{r_{E} + r_{p}}{2})^{2} \frac{Τ_{E} - Τ_{p}}{△ r} - (\frac{r_{p} + r_{w}}{2})^{2} \frac{Τ_{Ρ} - Τ_{W}}{△ r}}{△ r} (2) \end{array}$

1.2 模型初始条件和边界条件

初始条件 TP|t=0=T0 (3)

边界条件 $\frac{\partial Τ_{p}}{\partial r} |_{r = 0} = 0$ (4)

$- λ_{p} \frac{\partial Τ_{p}}{\partial r} |_{r = 0} = α (Τ_{p} |_{r = R} - Τ_{f})$ (5)

边界条件式(4)离散为

$\frac{- 3 Τ_{1} + 4 Τ_{2} - Τ_{3}}{2 △ r} = 0$ (6)

边界条件式(5)离散为

$(1 + \frac{3 λ}{2 △ r α}) Τ_{L 1} = \frac{2 λ}{△ r α} Τ_{L 2} - \frac{λ}{2 △ r α} Τ_{L 3} + Τ_{F}$ (7)

1.3 球状草莓的数值求解

数值求解必须先确定食品的物性,草莓的物性参数[6]如下:

水的质量分数ω/%:89.3

初始冻结温度Tm/℃:-0.89

比热容Cp (Siebel)/kJ·(kg·K) -1:

0.837+3.349ω 未冻相

0.837+1.256ω 冻结后

热导率/kJ·(m·℃)-1:λ=0.148+0.493ω

密度/ kg·m-3:ρ=(1-ωW)ρice+ωWρW

融化热/kJ·kg-1:302

用纯水的表面传热系数代表所有的食品表面传热系数,则

α=5.8+4.7ν0.84[1]

式中 α —草莓颗粒的表面传热系数(kJ/m2·℃);

ν —草莓颗粒的表面空气流速 (m/s)。

用三角分解和迭代法求解矩阵方程,根据上述公式编制了计算机程序。

2 数值模拟结果

依据草莓的物性参数,食品表面传热系数为α=5.8+4.7ν0.84。根据编制的球状草莓计算程序数值求解(初始条件为 20℃,中心-18℃,50层,△l=0.1),计算结果如表1所示。

3 结论

1) 通过对草莓在波形振动半流化食品速冻机中的冻结过程数值的模拟和比较分析,得到草莓不同风速和风温时的中心温度变化值。

2) 在数值模拟的基础上,分析了对冻结过程的影响因素。在较低风温时,冻结受风速影响较大,风速越高,冻结速度越快;而风速较低时冻结速度较慢。随冷风温度的升高,风速的影响逐渐减小。以上分析表明,风温和风速对冻结速度都有重要的影响,但其影响并非是线性的,在不同的风温与风速组合下,影响程度有所不同。提高风速会大大增加单位质量产品的能耗。从降低单位产品能耗的角度出发,采用较低的风速,通过提高振动强度,实现了单体的快速冻结,减少了运行费用。

3) 根据计算结果绘出了冷冻过程等时间温度曲线,使冷冻工艺的制定建立在更科学的基础上。

参考文献

[1]谢晶,徐世琼.水产品冻结时间的试验研究[J].制冷,1998(3):6-10.

[2]刘安源,刘中良,段钰锋.振动流化床干燥装置干燥特性计算方法研究[J].化学工程,1999,27(6):24-26.

[3]康景隆.快速冻结[M].北京:轻工业出版社,1996:117-125.

[4]陶文铨.数值传热学(2版)[M].西安:西安交通大学出版社,2001.

[5]钱壬章,俞昌铭,林文贵.传热分析与计算[M].北京:高等教育出版社,1987.

篇4：噪声扩散过程数值分析

关键词：噪声扩散数值模型

1 噪声扩散方式

关于噪声的扩散方式，可以借由信息传播进行描述，因为信息和噪声的上层概念都是消息，二者的传播方式非常相近。

基于非确定信息网络结构的研究，非确定性网络是指个体间通过随机碰撞传递信息。非确定性信息网络刻画了金融系统信息传播的非同时性、异质性。Cont和Bouchaud （2000）提出了基于完全随机碰撞交互的金融系统信息传播模型。该模型的优点是和真实金融系统的情况非常接近，个体无需任何关于整个系统的统计信息就可以进行决策。模型较好地刻画了群体性从众行为从形成、演化、消失的动态规律。应尚军等（2001，2003）[1][2]、杨春霞等（2005）[3]改进了个体间交互的规则，通过多种信息传播模型的构建，刻画了信息传播的途径、速度、范围等动态特征，还模拟了价格泡沫的形成和破灭。

基于特定复杂网络的研究。这类研究的思路是由研究者根据其对金融市场的观察与理解，指定金融系统遵循小世界网络或无标度网络模式，进而探索复杂网络演化特征及其对资产价格的影响。Hein和Schwind （2008）等研究了基于小世界网络的股市信息交互结构，发现了个体数量和信息交互强度的变化将影响价格形成过程。陈彦锟（2010）[4]基于无标度网络建立的信息传播结构，研究了泡沫的形成与崩馈的条件。

此外，林俊波（2005）[5]借用Shannon的无线电信号传递原理模型，界定证券市场的信源、信道和信宿，构建了证券市场信息传递模型，比较全面而抽象地概括了证券市场信息扩散过程。

邓忆瑞（2008）[6]借用物理学中“场”的概念，建立了信息扩散场，基于场论，采用逻辑推理和数学分析相结合的方法，构建并求解信息扩散场的扩散状态模型，将信息扩散机理用“场”语言描述出来，并利用马氏漂移链原理建立信息时空扩散模型，描述了信息扩散的时间扩展规律与空间分布特征，并利用计算机对模型进行了模拟。

宋逢明等（2002）[7]对中国股票市场的收益率与交易量进行大量实证研究的基础上，构造了中国股票市场区别于成熟资本市场的特殊信息传导模型，并验证了模型的适用性，还发现不同类型股票的投资者的构成不同，其信息传导结构也不同。

吴忠群（2011）描述了两种特殊的信息模式，泄漏式和爆炸式。信息以渐次传递方式被个体获知称为“泄漏式”，信息瞬间扩散到每一个系统中每个个体称为“爆炸式”。信息的传播途径也可分为通过媒体披露和通过个人交流。这样，不同传播速度和传播途径交织组合在一起，构成信息传播方式的多样性。

Kosfeld（2004）[8]用数理推导呈现了一个关于谣言对市场影响的模型，为其他关于谣言和资产价格的实证提供理论分析基础。模型中个体通过与周围邻居进行交流产生了谣言传播的可能。推导的结果显示，谣言最终消失，长期均衡价格等于谣言前价格；如果谣言仍存在，会造成与谣言有关资产的价格上升。

林春燕和朱东华（2005）[9]从股市内部信息传播如何影响股价的波动和交易量变化为研究角度，建立常微分方程组描述了在未知情者不具备学习能力的情况下，信息自身传播过程，并预测传播高峰，分析其过程中股票交易量的变化情况。该模型证明了两个定理：一是传播信息的投资者人数先单调增加，然后单调减少并趋于零；而是总有一部分人在信息公布前无法得知内部消息。

2 噪声的扩散过程建模

股市噪声从噪声源产生以后总是要流动的，不存在静止状态的噪声。所谓噪声流动也就是噪声的扩散，而噪声扩散又可称为噪声的传递，是指噪声从噪声源产生以后，经过传递渠道送达给信息的全部活动过程。

股市噪声扩散依据传播速度分类可分为两种特殊形式，泄漏式和爆炸式。泄漏式代表信息是以渐次传递方式被公众获知的，爆炸式代表信息是瞬间扩散到每一个受众的；依据传播途径分类可分为，通过媒体披露以及个人之间交流。不同传播速度和传播途径交织组合在一起，可以构成信息传播方式的多样性。

本章重点研究以“泄漏式”为传播方式和以“个人交流”为传播途径的噪声的扩散过程。此类噪声是以“人”为载体进行传播的。每个传播者只通过周围有限个人渐次将噪声扩散出去。本章将建立模型，研究噪声传播者的变化规律来说明此类股市噪声的扩散特征，可以为控制此类噪声扩散提供一定的理论支持。

2.1 模型基本假设

设投资者集合I，噪声产生时刻t。t时刻后，投资者开始分化成三个种群。第一类投资者在t时刻获知噪声并正在转告他人，称为噪声传播者，记为集合D；第二类投资者在t时刻获知噪声但不转告他人，称为传播终止者，记为集合K；第三类投资者不知晓该噪声，称为不知情者，记为集合U，并且有：

I=D+K+U （2-1）

设t时刻上述三类投资者占总投资者的比例分别用D（t）、K（t）、U（t）表示，所以有：

D（t）+K（t）+U（t）=1 （2-2）

对模型作出如下假设：

①噪声传播者通过有限个个体交流，逐渐将噪声传播出去；

②变量D（t）、K（t）、U（t）为连续可微变量。

③只有噪声传播者D传播给不知情者U才算有效传播。

假设③说明，噪声从人群D到D、D到K两种传播都是无效的，因为D、K两类人群均已知此噪声。

设噪声传播者平均传播率为常数λ，λ为大于等于1的常数，意味着单位时间内，平均每个噪声传播者有向λ个人传播噪声的传播能力，但只有向人群U传播才是有效传播，故噪声传播者在单位时间内传播噪声的有效能力为：

λU（t）（2-3）

设总体投资者的人数为N，单位时间内传播者传播噪声的有效传播人数为：

λNU（t）D（t）（2-4）

上式也表示噪声传播者D单位时间的增加量，即噪声传播者D增加的速度，可以用它来衡量噪声传播的速度。它是U（t）、D（t）的函数。

④传播者D会转换为传播终止者K，传播终止率正比于噪声传播者数量ND，比例系数为μ，μ为大于0且小于1的常数，所以噪声传播者D单位时间减少量为：

μND（t）（2-5）

A4：不知情者具有学习能力。当股价波动，成交量扩大时，会引起持有或对该股票感兴趣的投资者的关注。当有高于（或低于）均衡价的大笔订单成交时，将引起原本不知情的交易者改变他们对该股票的预期价格，有可能对该股进行买卖。噪声传播者越活跃，不知情者参与的概率越大。

设不知情者参与的概率为P（t），随时间的变化而变化，且与D（t）成正比，为不知情者学习能力参数，为大于0且小于1的常数，则有：

（2-6）

学习能力强的投资者会从交易量的变化中猜测噪声，并加入到传播人群中去，这部分人群的单位时间增加量为：

NP（t）U（t）（2-7）

2.2 噪声传播模型

2.2.1 噪声传播方程

当噪声出现时，系统内种群开始分化。不知情者接收到噪声，变成噪声传播者，将噪声扩散出去；之后噪声传播者不再传播，变成传播终止者。

投资者的角色转换过程是U?D?K，三类投资者数量随时间推移呈现“此消彼长”的特征。三类投资者所占比例的变化速率分别记为D′（t）、K′（t）、U′（t）。建立常微分方程组（2-8）如下：

U′（t）=-λU（t）D（t）-p（t）U（t）

D′（t）=λU（t）D（t）-μD（t）+p（t）U（t）

K′（t）=μD（t）

P（t）=D（t）

U（t）+D（t）+K（t）=1

U0>0，D0>0，K0>0（2-8）

2.2.2 种群的演化特征一

噪声是通过个体交流的方式逐渐在人群中扩散的，所以噪声扩散者应该是逐渐增加的。随着时间的推移，不断有新消息进入系统取代先前的噪声成为决策依据。噪声逐渐沦为过时的消息，传播者没有动力去传播噪声，表现为人群D数量的减少。

所以噪声有如下传播特征一：

噪声传播者的种群人数变化规律有两种可能：或单调减少，或先逐渐增加，达到峰值后逐渐减少。

2.2.3 种群的演化特征二

从现实情况来看，噪声的扩散不会至整个系统。由于信息不对称总是存在，交易者在有关交易信息之数量和质量的拥有上不相等，股市消息不会到达每一个人。也就是说，人群中总有人不知道噪声。

如下传播特征二：

不知情者不会随时间的推移而完全消失。

2.2.4 噪声扩散的峰值

根据噪声传播者的变化规律可知，传播者人数必然会在某一时刻达到峰值。到达峰值的时间与初始传播人数、传播终止率、噪声传播者的传播能力和不知情者的学习能力有关。一般地，初始传播人数越大，噪声传播者的传播能力越强，不知情交易者学习能力越强，峰值就会越大，到达峰值的时间就越早；传播终止率则与峰值成负相关。

3 数值模拟与分析

3.1 分析工具的选择

MATLAB是使用较为广泛的模拟工具，提供了7个求常微分方程数值解的函数：ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb。ode45是解决数值解问题的首选方法，将选择ode45进行方程求解。同时，将选择二维画图函数plot进行函数曲线绘制。

3.2 模拟分析

用MATLAB求解微分方程组，并画出D（t）、K（t）、U（t）三个种群人数演化曲线（演化曲线略）。模型中需要赋值的变量有：

①噪声传播者的初始比例D0在0～1内由随机数生成；

②噪声传播者的初始比例D0与未知情者的初始比例U0之和为1；

③传播终止者的初始比例K0为0；

④参数λ大于或等于1的常数；

⑤参数μ、的取值范围是（0，1）。

3.3 结果分析

①初始时刻，系统中噪声传播者为固定值，是整个系统的噪声传播源，其余人群则是不知情者。当噪声产生后，传播者开始传播噪声，不知情者人数逐渐减少。随着噪声的传播，知晓噪声但不传播的传播终止者人数从0开始逐渐增多。噪声传播者的变化规律与传播能力、传播终止率、未知情者学习能力和初始噪声传播人数都有关。

②其他条件不变，初始状态下，系统中知晓并传播噪声的人数越多，传播高峰来的越快，峰值越大。这符合逻辑推理和实际观察。

③其他条件不变，当传播者传播噪声的能力越强，传播高峰来的就越快，峰值就越大，符合逻辑推理和实际观察。

④其他条件不变，当传播终止率越小，峰值越大，但传播高峰时间变化不明显。

⑤其他条件不变，当未知情者学习能力越大，峰值越大，传播高峰时间变化也不明显。

⑥曲线均未涉及当时间足够长时，三类人群的数量特征。但上一章经过理论已经证明，噪声传播者最终会消失，整个系统只有未知情交易者和传播终止者。

参考文献：

[1]应尚军，魏一鸣，范英，等. 基于元胞自动机的股票市场复杂性研究——投资者心理与市场行为[J].系统工程理论与实践，2003，（12）：18-24+31.

[2]应尚军，魏一鸣，范英，等.基于元胞自动机的股票市场投资行为模拟[J].系统工程学报，2001（5）：382-388.

[3]杨春霞，王杰，周涛，等.基于自组织逾渗的金融市场模型[J].科学通报，2005（20）：127-131.

[4]陈彦锟.基于无标度网络的信用违约风险传染效应研究[J].统计与决策，2010（2）：20-23.

[5]林俊波.证券市场信息传导机制与信息披露制度研究[D].浙江大学，2005.

[6]邓忆瑞.基于网络维力的信息扩散研究[D].哈尔滨工程大学， 2008.

[7]宋逢明，唐俊.中国股票市场的信息传导与流动性需求[J].经济科学，2002（2）：46-57.

[8]Kosfeld M.Rumours and markets[J].Journal of Mathematical Economics，2005，V41（6）：646-664.

[9]林春燕，朱东华.证券市场信息传播的数学模型研究[J].数学的实践与认识，2005（11）：40-45.

作者简介：

刘芷彤（1992-），女，内蒙古通辽人，大学本科，学士，华北电力大学经济与管理学院。

篇5：中国农业大学数值分析

一、课程设计的目的与意义

本课程设计是配合《数值分析》课程而开设的一门实践课程，是培养学生分析和解决实际问题能力的重要实践教学环节，重点培养学生运用计算、Matlab软件对数值分析课程中的数值计算方法、实例进行编程实现，对巩固和深化学生所学课程的理论知识，增强学生的应用能力具有重要作用。

二、课程设计的内容

1、本课程设计的选题内容应属课程范围，能达到课程教学的目的和满足课程教学的要求，使学生得到较全面的综合训练。

2、题目的难度和工作量应适合学生掌握的知识和具备的能力，使学生在规定的时间内既工作量饱满，又能经过努力完成任务。

三、课程设计的方式

以集中学习为主；独立完成课程设计阶段规定的全部工作任务。

四、课程设计的基本要求

1、注重对学生运用所学知识，解决实际问题的能力和创新精神的培养。

2、完成一篇论文的撰写，不少于5000字。

3、课程设计的具体要求与格式应符合本学院规定。

4、严格遵守学习纪律，按时完成规定任务。

5、学生因特殊原因请假须履行手续，凡未请假或未获批准擅自离岗者，均按旷课处理。

五、课程设计成绩的评定

1、课程设计成绩采用五级分制：优、良、中、及格、不及格。

篇6：数值分析总结

2.数值分析的误差种类：1）截断误差：模型的准确解与数值方法准确解之间的误差。

2）舍入误差：实数形式的原始数据与有限字长计算机数据间的误差。

3.算法的数值稳定性与病态问题：1）若某算法受初始误差或运算过程中的舍入误差影响较小，则称为数值稳定。2）若微小的初始误差都会对最终结果产生极大的影响，则称之为病态问题。

二：1.Runge现象及其解决方法

Runge现象即高次插值的振荡现象，指增加节点固然能使插值函数 p(x)与被插值函数f（x）在更多的地方相等，但在两点之间p(x)不一定能很好地近似f（x），有时候误差非常大。

解决方法：分段低次插值（将插值区间分成若干小区间，在小区间内用低次插值）

2.样条插值思想：插值函数p(x)在插值区间[a,b]上有二阶光滑度，在分段的小区间

[xk,xk+1]上是低次多项式，同时满足p(xi)=yi.三：理解逼近问题与拟合问题：1）逼近问题：函数f(x)在区间[a,b]具有一阶光滑度，求多项式p(x)是f(x)-p(x)在某衡量标准下最小的问题。2）拟合问题：从理论上讲y=f(x)是客观存在的，但在实际中，仅仅从一些离散的数据（xi,yi）(i=1,2…)是不可能求出f(x)的准确表达式，只能求出其近似表达式φ（x）。

插值问题与逼近问题的特点和区别：1）相同点：它们都是求某点值的算法。

2）不同点：A，被插值函数是未知的，而被逼近函数是已知的。B，插值函数在节点处与被插值函数相等。而逼近函数的值只要满足很好的均匀逼近即可。C，求p(x)的方法不同。

四：Romberg求积法和Gauss求积法的基本思想：

1）复化求积公式精度较高，但需要事先确定步长，欠灵活性，在计算过程中将步长逐次减半得到一个新的序列，用此新序列逼近I的算法为Romberg求积法。

2）对插值型求积公式，若能选取适当的xk.Ak使其具有2n+1阶代数精度，则称此类求积公式为Gauss型。

五.Runge-Kutta方法的基本思想：

篇7：《数值分析》课程实验报告

学号：

学院：

机电学院日期：

2015 年 X 月X 日目录实验一函数插值方法 1 实验二函数逼近与曲线拟合 5 实验三数值积分与数值微分 7 实验四线方程组的直接解法 9 实验五解线性方程组的迭代法 15 实验六非线性方程求根 19 实验七矩阵特征值问题计算 21 实验八常微分方程初值问题数值解法 24 实验一函数插值方法一、问题提出对于给定的一元函数的n+1个节点值。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。

数据如下：

（1）0.4 0.55 0.65 0.80 0.95 1.05 0.41075 0.57815 0.69675 0.90 1.00 1.25382 求五次Lagrange多项式，和分段三次插值多项式，计算, 的值。（提示：结果为,）（2）1 2 3 4 5 6 7 0.368 0.135 0.050 0.018 0.007 0.002 0.001 试构造Lagrange多项式，计算的,值。（提示：结果为,）二、要求 1、利用Lagrange插值公式编写出插值多项式程序；

2、给出插值多项式或分段三次插值多项式的表达式；

3、根据节点选取原则，对问题（2）用三点插值或二点插值，其结果如何；

4、对此插值问题用Newton插值多项式其结果如何。Newton插值多项式如下：

其中：

三、目的和意义 1、学会常用的插值方法，求函数的近似表达式，以解决其它实际问题；

2、明确插值多项式和分段插值多项式各自的优缺点；

3、熟悉插值方法的程序编制；

4、如果绘出插值函数的曲线，观察其光滑性。

四、实验步骤（1）0.4 0.55 0.65 0.80 0.95 1.05 0.41075 0.57815 0.69675 0.90 1.00 1.25382 求五次Lagrange多项式，和分段三次插值多项式，计算, 的值。（提示：结果为,）第一步：先在matlab中定义lagran的M文件为拉格朗日函数代码为：

function[c,l]=lagran(x,y)w=length(x);n=w-1;l=zeros(w,w);for k=1:n+1 v=1;for j=1:n+1 if(k~=j)v=conv(v,poly(x(j)))/(x(k)-x(j));end end l(k,:)=v;end c=y*l;end 第二步：然后在matlab命令窗口输入：

>>>> x=[0.4 0.55 0.65 0.80,0.95 1.05];y=[0.41075 0.57815 0.69675 0.90 1.00 1.25382];>> lagran(x,y)回车得到：

ans =121.6264-422.7503 572.5667-377.2549 121.9718-15.0845 由此得出所求拉格朗日多项式为 p（x）=121.6264x5-422.7503x4+572.5667x3-377.2549x2+121.9718x-15.0845 第三步：在编辑窗口输入如下命令：

>> x=[0.4 0.55 0.65 0.80,0.95 1.05];>> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.0845;>> plot(x,y)命令执行后得到如下图所示图形，然后 >> x=0.596;>> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.084 y =0.6262 得到f（0.596）=0.6262 同理得到f（0.99）=1.0547（2）1 2 3 4 5 6 7 0.368 0.135 0.050 0.018 0.007 0.002 0.001 试构造Lagrange多项式，和分段三次插值多项式，计算的,值。（提示：结果为,）实验步骤：

第一步定义 function[c,l]=lagran(x,y)w=length(x);n=w-1;l=zeros(w,w);for k=1:n+1 v=1;for j=1:n+1 if(k~=j)v=conv(v,poly(x(j)))/(x(k)-x(j));end end l(k,:)=v;end c=y*l;end 定义完拉格朗日M文件第二步：然后在matlab命令窗口输入：

>>>> x=[1 2 3 4 5 6 7];y=[0.368 0.135 0.050 0.018 0.007 0.002 0.001];>> lagran(x,y)回车得到：

ans =0.0001-0.0016 0.0186-0.1175 0.4419-0.9683 0.9950 由此得出所求拉格朗日多项式为 p（x）=0.0001x6-0.0016x5+0.0186x4-0.1175x3+0.4419x2-0.9683x+0.9950 第三步：在编辑窗口输入如下命令：

>> x=[1 2 3 4 5 6 7];>> y=0.0001*x.^6-0.0016*x.^5+0.0186*x.^4-0.1175*x.^3+0.4419*x.^2-0.9683*x+0.9950;>> plot(x,y)命令执行后得到如下图所示图形，然后 >> x=1.8;>> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.084 y =0.1650 得到f（0.596）=0.6262 同理得到f（6.15）=2.3644 五、实验结论插值是在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点，它是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。

实验二函数逼近与曲线拟合一、问题提出从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。

在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。

t(分)0 5 10 15 20 25 30 35 40 45 50 55 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64 二、要求 1、用最小二乘法进行曲线拟合；

2、近似解析表达式为；

3、打印出拟合函数，并打印出与的误差，；

4、另外选取一个近似表达式，尝试拟合效果的比较；

5、* 绘制出曲线拟合图。

三、目的和意义 1、掌握曲线拟合的最小二乘法；

2、最小二乘法亦可用于解超定线代数方程组；

3、探索拟合函数的选择与拟合精度间的关系四、实验步骤：

第一步先写出线性最小二乘法的M文件 function c=lspoly(x,y,m)n=length(x);b=zeros(1:m+1);f=zeros(n,m+1);for k=1:m+1 f(:,k)=x.^(k-1);end a=f＇*f;b=f＇*y＇;c=a＼b;c=flipud(c);第二步在命令窗口输入：

>>lspoly([0,5,10,15,20,25,30,35,40,45,50,55],[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64],2)回车得到：

ans =-0.0024 0.2037 0.2305 即所求的拟合曲线为y=-0.0024x2+0.2037x+0.2305 在编辑窗口输入如下命令：

>> x=[0,5,10,15,20,25,30,35,40,45,50,55];>> y=-0.0024*x.^2+0.2037*x+0.2305;>> plot(x,y)命令执行得到如下图五、实验结论分析复杂实验数据时,常采用分段曲线拟合方法。利用此方法在段内可以实现最佳逼近,但在段边界上却可能不满足连续性和可导性。分段函数的光滑算法,给出了相应的误差分析.给出了该方法在分段曲线拟合中的应用方法以及凸轮实验数据自动分段拟合。

实验三数值积分与数值微分一、问题提出选用复合梯形公式，复合Simpson公式，Romberg算法，计算（1）（2）（3）（4）二、要求 1、编制数值积分算法的程序；

2、分别用两种算法计算同一个积分，并比较其结果；

3、分别取不同步长，试比较计算结果（如n = 10, 20等）；

4、给定精度要求ε，试用变步长算法，确定最佳步长。

三、目的和意义 1、深刻认识数值积分法的意义；

2、明确数值积分精度与步长的关系；

3、根据定积分的计算方法，可以考虑二重积分的计算问题。

四、实验步骤第一步：编写各种积分的程序复合梯形程序如下：

function I=TX(x,y)n=length(x);m=length(y);if n~=m error(＇The lengths of X and Y must be equal＇);return;end h=(x(n)-x(1))/(n-1);a=[1 2*ones(1,n-2)1];I=h/2*sum(a.*y);复合Simpson程序如下：

function s = simpr1(f,a,b,n)h=(b-a)/(2*n);s1=0;s2=0;for k=1:10 x=a+h*(2*k-1);s1=s1+feval(f,x);end for k=1:(10-1)x=a+h*2*k;s2=s2+feval(f,x);end s=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3;end Romberg程序如下：

function I = Romber_yang(fun,a,b,ep)if nargin<4 ep=1e-5;end;m=1;h=b-a;I=h/2*(feval(fun,a)+feval(fun,b));T(1,1)=I;while 1 N=2^(m-1);h=h/2;I=I/2;for i=1:N I=I+h*feval(fun,a+(2*i-1)*h);end T(m+1,1)=I;M=2*N;k=1;while M>1;T(m+1,k+1)=(4^k*T(m+1,k)-T(m,k))/(4^k-1);M=M/2;k=k+1;end if abs(T(k,k)-T(k-1,k-1))

2、对于积分Ι=01sin⁡XXdx，f(0)=1,梯形积分T=0.94607307，辛普森积分S=0.94607308，Romberg积分R=0.94607307。

3、对于积分Ι=01eX4+X2dx，梯形积分T=0.39081248，辛普森积分S=0.39081185，Romberg积分R=0.39081885。

4、对于积分Ι=01ln1+X1+X2dx，梯形积分T=0.27218912，辛普森积分S=0.27219844，Romberg积分R=0.27219827。

五、实验结论，通过本实验学会复合梯形公式，复合Simpson公式，Romberg公式的编程与应用，掌握MATLAB提供的计算积分的各种函数的使用方法。

实验四线方程组的直接解法一、问题提出给出下列几个不同类型的线性方程组，请用适当算法计算其解。

1、设线性方程组  2、设对称正定阵系数阵线方程组  3、三对角形线性方程组二、要求 1、对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法；

平方根法与改进平方根法；

追赶法求解（选择其一）；

2、应用结构程序设计编出通用程序；

3、比较计算结果，分析数值解误差的原因；

4、尽可能利用相应模块输出系数矩阵的三角分解式。

三、目的和意义 1、通过该课题的实验，体会模块化结构程序设计方法的优点；

2、运用所学的计算方法，解决各类线性方程组的直接算法；

3、提高分析和解决问题的能力，做到学以致用；

4、通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。

四、实验步骤：

列主元高斯消去法的matlab的M文件程序 function [x,det,index]=Gauss(A,b)% 求线形方程组的列主元Gauss消去法，其中，% A为方程组的系数矩阵；

% b为方程组的右端项；

% x为方程组的解；

% det为系数矩阵A的行列式的值；

% index为指标变量，index=0表示计算失败，index=1表示计算成功。

[n,m]=size(A);nb=length(b);% 当方程组行与列的维数不相等时，停止计算，并输出出错信息。

if n~=m error(＇The rows and columns of matrix A must be equal!＇);return;end % 当方程组与右端项的维数不匹配时，停止计算，并输出出错信息 if m~=nb error(＇The columns of A must be equal the length of b!＇);return;end % 开始计算，先赋初值 index=1;det=1;x=zeros(n,1);for k=1:n-1 % 选主元 a_max=0;for i=k:n if abs(A(i,k))>a_max a_max=abs(A(i,k));r=i;end end if a_max<1e-10 index=0;return;end % 交换两行 if r>k for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;end z=b(k);b(k)=b(r);b(r)=z;det=-det;end % 消元过程 for i=k+1:n m=A(i,k)/A(k,k);for j=k+1:n A(i,j)=A(i,j)-m*A(k,j);end b(i)=b(i)-m*b(k);end det=det*A(k,k);end det=det*A(n,n);% 回代过程 if abs(A(n,n))<1e-10 index=0;return;end for k=n:-1:1 for j=k+1:n b(k)=b(k)-A(k,j)*x(j);end x(k)=b(k)/A(k,k);end 然后在命令窗口输入 >> A=[4 2-3-1 2 1 0 0 0 0;8 6-5-3 6 5 0 1 0 0;4 2-2-1 3 2-1 0 3 1;0-2 1 5-1 3-1 1 9 4;-4 2 6-1 6 7-3 3 2 3;8 6-8 5 7 17 2 6-3 5;0 2-1 3-4 2 5 3 0 1;16 10-11-9 17 34 2-1 2 2;4 6 2-7 13 9 2 0 12 4;0 0-1 8-3-24-8 6 3-1];>> b=[5 12 3 2 3 46 13 38 19-21];>> gauss(A,b)ans = 1.0000-1.0000 0.0000 1.0000 2.0000 0.0000 3.0000 1.0000-1.0000 2.0000 高斯-约当消去法maltab的M文件程序 function [x,flag]=Gau_Jor(A,b)% 求线形方程组的列主元Gauss-约当法消去法，其中，% A为方程组的系数矩阵；

% b为方程组的右端项；

% x为方程组的解；

[n,m]=size(A);nb=length(b);% 当方程组行与列的维数不相等时，停止计算，并输出出错信息。

>> A=[4 2-4 0 2 4 0 0;2 2-1-2 1 3 2 0;-4-1 14 1-8-3 5 6;0-2 1 6-1-4-3 3;2 1-8-1 22 4-10-3;4 3-3-4 4 11 1-4;0 2 5-3-10 1 14 2;0 0 6 3-3-4 2 19];>> b=[0-6 20 23 9-22-15 45];>> Gau_Jor(A,b)ans = 121.1481-140.1127 29.7515-60.1528 10.9120-26.7963 5.4259-2.0185 五、实验结论用LU法，调用matlab中的函数lu中，L往往不是一个下三角，但可以直接计算不用它的结果来计算，不用进行行变换。如果进行行变b也要变，这样会很麻烦。

实验五解线性方程组的迭代法一、问题提出对实验四所列目的和意义的线性方程组，试分别选用Jacobi 迭代法，Gauss-Seidel迭代法和SOR方法计算其解。

二、要求 1、体会迭代法求解线性方程组，并能与消去法做以比较；

2、分别对不同精度要求，如由迭代次数体会该迭代法的收敛快慢；

3、对方程组2，3使用SOR方法时，选取松弛因子ω=0.8，0.9，1，1.1，1.2等，试看对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者；

4、给出各种算法的设计程序和计算结果。

三、目的和意义 1、通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较；

2、运用所学的迭代法算法，解决各类线性方程组，编出算法程序；

3、体会上机计算时，终止步骤或k >（给予的迭代次数），对迭代法敛散性的意义；

4、体会初始解，松弛因子的选取，对计算结果的影响。

四、实验步骤第一步编写实验所需的Jacobi迭代法，Gauss-Seidel迭代法，SOR迭代法的程序。

Jacobi迭代法：

function [x,k,index]=J(A,b,ep,itmax)if nargin<4 itmax=100;end if nargin<3 ep=1e-5;end n=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while 1 for i=1:n y(i)=b(i);for j=1:n if j~=i y(i)=y(i)-A(i,j)*x(j);end end if abs(A(i,i))<1e-10|k==itmax index=0;return;end y(i)=y(i)/A(i,i);end if norm(y-x,inf)

function [x,k,index]=G(A,b,ep,itmax)if nargin<4 itmax=100;end if nargin<3 ep=1e-5;end n=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while 1 y=x;for i=1:n z=b(i);for j=1:n if j~=i z=z-A(i,j)*x(j);end end if abs(A(i,i))<1e-10|k==itmax index=0;return;end z=z/A(i,i);x(i)=z;end if norm(y-x,inf)

function [x,k,index]=SOR(A,b,ep,w,itmax)if nargin<5 itmax=100;end if nargin<4 w=1;end if nargin<3 ep=1e-5;end n=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while 1 y=x;for i=1:n z=b(i);for j=1:n if j~=i z=z-A(i,j)*x(j);end end if abs(A(i,i))<1e-10|k==itmax index=0;return;end z=z/A(i,i);x(i)=(1-w)*x(i)+w*z;end if norm(y-x,inf)

1、设线性方程组  2、设对称正定阵系数阵线方程组 3、三对角形线性方程组五、实验结论迭代法是解线性方程组的一个重要的实用方法，特别适用于求解在实际中大量出现的，系数矩阵为稀疏阵的大型线性方程组。通过此次实验学会了Jacobi迭代法，Gauss-Seidel迭代法，SOR迭代法的程序编写，并掌握了它们各自的优缺点及其适用条件。

实验六非线性方程求根一、问题提出设方程有三个实根现采用下面六种不同计算格式，求 f(x)=0的根或 1、2、3、4、5、6、二、要求 1、编制一个程序进行运算，最后打印出每种迭代格式的敛散情况；

2、用事后误差估计来控制迭代次数，并且打印出迭代的次数；

3、初始值的选取对迭代收敛有何影响；

4、分析迭代收敛和发散的原因。

三、目的和意义 1、通过实验进一步了解方程求根的算法；

2、认识选择计算格式的重要性；

3、掌握迭代算法和精度控制；

4、明确迭代收敛性与初值选取的关系。

四、实验步骤第一步：编写实验所需的程序。

function [x_star,index,it]=DD(fun,x,ep,itmax)if nargin<4 itmax=100;end if nargin<3 ep=1e-5;end index=0;k=1;while k

1、，x1=0，x2=0。

2、，x1=无穷大，x2=-0.3473。

3、，x1=1.8794，x2=1.8794。

4、，x1=-0.3473，x2=-0.3473.。

5、，x1=1.8794，x2=1.8794。

6、，x1=1.8794，x2=-0.3473。

五、实验结论对于非线性方程，求它的解析解有时候是很困难的，但采用数值方法可以很容易地求它的近似解。此次实验就是采用迭代法求非线性方程的根。对于一个非线性方程，选用不同的迭代形式，因为其收敛程度不一样，造成其效率与精确度有很大的差别。

实验七矩阵特征值问题计算一、问题提出利用冪法或反冪法，求方阵的按模最大或按模最小特征值及其对应的特征向量。

设矩阵A的特征分布为：

且试求下列矩阵之一（1）求，及取结果（2）求及取结果：

（3）求及取结果(4)取这是一个收敛很慢的例子,迭代次才达到结果（5）有一个近似特征值,试用幂法求对应的特征向量，并改进特征值（原点平移法）。

取结果二、要求 1、掌握冪法或反冪法求矩阵部分特征值的算法与程序设计；

2、会用原点平移法改进算法，加速收敛；

对矩阵B=A-PI取不同的P值，试求其效果；

3、试取不同的初始向量，观察对结果的影响；

4、对矩阵特征值的其它分布，如且如何计算。

三、目的和意义 1、求矩阵的部分特征值问题具有重要实际意义，如求矩阵谱半径，稳定性问题往往归于求矩阵按模最小特征值；

2、进一步掌握冪法、反冪法及原点平移加速法的程序设计技巧；

3、问题中的题（5），反应了利用原点平移的反冪法可求矩阵的任何特征值及其特征向量。

四、实验步骤第一步：写出实验所需的幂法求最大特征值及反幂法求最小特征值的程序。

幂法程序：

function [m,u,index]=TZ(A,ep,itmax)if nargin<3 itmax=100;end if nargin<2 ep=1e-5;end n=length(A);u=ones(n,1);index=0;k=0;m1=0;while k<=itmax v=A*u;[vmax,i]=max(abs(v));m=v(i);u=v/m;if abs(m-m1)

function [m,u,index]=FTZ(A,ep,itmax)if nargin<3 itmax=100;end if nargin<2 ep=1e-5;end n=length(A);u=ones(n,1);index=0;k=0;m1=0;invA=inv(A);while k<=itmax v=invA*u;[vmax,i]=max(abs(v));m=v(i);u=v/m;if abs(m-m1)

λ3=3.4723，x3=（1.0000 0.5229 0.2422）T。，λ1=21.3053，X1=（0.8724 0.5401 0.9973 0.5644 0.4972 1.0000）T；

λ6=1.6214。

五、实验结论求n阶方阵A的特征值和特征向量，也是实际中常常碰到的问题。通过此次实验掌握了用幂法和反幂法求一个方阵的最大特征值和特征向量，绝对值最小的特征值和特征向量。

实验八常微分方程初值问题数值解法一、问题提出科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler法，Rung-Kutta方法求其数值解，诸如以下问题：

（1）分别取h=0.1,0.2,0.4时数值解。

初值问题的精确解。

（2）用r=3的Adams显式和预-校式求解取步长h=0.1，用四阶标准R-K方法求值。

（3）用改进Euler法或四阶标准R-K方法求解取步长0.01,计算数值解，参考结果。

（4）利用四阶标准R-K方法求二阶方程初值问题的数值解（I）(II)(III)(IV) 二、要求 1、根据初值问题数值算法，分别选择二个初值问题编程计算；

2、试分别取不同步长，考察某节点处数值解的误差变化情况；

3、试用不同算法求解某初值问题，结果有何异常；

4、分析各个算法的优缺点。

三、目的和意义 1、熟悉各种初值问题的算法，编出算法程序；

2、明确各种算法的精度与所选步长有密切关系；

3、通过计算更加了解各种算法的优越性。

四、实验步骤 function [x,y]=euler(fun,x0,xfinal,y0,n);if nargin<5,n=50;end h=(xfinal-x0)/n;x(1)=x0;y(1)=y0;for i=1:n;x(i+1)=x(i)+h;y(i+1)=y(i)+h*feval(fun,x(i),y(i));end 实验程序及分析（Ⅰ）（1）、算法程序 function E =Euler_1(fun,x0,y0,xN,N)% Euler向前公式，其中 % fun为一阶微分方程的函数 % x0,y0为初始条件 % xN为取值范围的一个端点 % h为区间步长 % N为区间个数 % x为Xn构成的向量 % y为yn构成的向量 x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;for n=1:N x(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(fun,x(n),y(n));end T=[x＇,y＇] function z=f(x,y)z=4*x/y-x*y;(2)、运行程序 >> Euler_1(＇f＇,0,3,2,20)结果：

篇8：轮式离心风扇数值模拟分析

在小型离心风扇中, 叶轮主要分为轮式和板式。本文主要分析轮式离心风扇的内部流场。

流体机械内部流场一般比较复杂, 很难进行可视化研究。一般都是靠实验来测试风扇P-Q曲线性能来进行完成。但对于内部流场的认知不是很清楚。随着计算机技术的发展, 近些年来数值分析方法也慢慢应用在流体机械中, 这对于风扇内部以及流体机械内部流场的分析起到了至关重要的作用。

计算流体力学 (Computational Fluid Dynamics, 简称CFD) 是20世纪60年代起伴随计算机技术迅速崛起的学科。经过半个世纪的迅猛发展, 目前这门学科己相当成熟, 而成熟的一个重要标志就是近十几年来, 各种CFD通用性软件包陆续出现, 成为商品化软件, 为工业界广泛接受, 性能日趋完善, 其应用范围不断扩大。至今, CFD技术的应用早己超越传统的流体力学和流体工程的范畴。

1 轮式离心风扇的有限元模型

1.1 轮式离心风扇有限元模型的建立

在轮式离心风扇的结构设计中主要考虑轮式扇叶由于其结构强度不如板式结构, 所以在生产注塑过程中变形量要大一些。在分析了该叶轮的结构特点和材料的力学性能的基础上, 采用Pro-e建立了离心风机3D模型, 并用Gambit进行前置处理建立有限元模型 (如图1) 。

由于小型散热风扇主要应用在笔记本系统中, 整个结构不是很规程存在的棱角和不规则曲线部分比较多, 在进行网格化前, 需要对这部分设计进行改进, 否则网格生成非常困难, 甚至网格化不成功。在进行风扇的网格化前, 必须处理完毕。

处理后建立的的3D模型如图1所示。

轮式扇叶的有限元模型如图2所示。

离心风扇整体的有限元模型如图3所示。

1.2 边界条件和计算模型

由于小型离心风扇主要应用在笔记本散热系统, 所留给风扇的空间有限。设定风扇在额定工况下工作:由于小型风扇内部风速不大。内部流场认为是不可压缩稳定流场, 模型计算时忽略空气重力对于流场的影响。离心风扇空气入口边界条件采用压力入口 (PRESSURE_INLET) ;并严格AMCA 2210285标准测试规范建立了数值风洞模型, 离心风扇空气空气出口边界采用压力出口 (PRESSURE_OUTLET) 。

该离心风扇中, 其额定转速定义为3200RPM。在风扇模型中定义两种网格区域, 分别是流道静网格区和扇叶动网格区, 所以在流道和扇叶之间采用MRF多重参考系统进行耦合 (Multi-reference.Frame) , 压力速度离散采用Sample算法和k-ε湍流模型。

2 叶轮离心风扇的模拟结果分析

在作离心风扇模型数值分析之前, 首先需要决定迭代次数的选择。迭代次数选择的小, 模拟计算结果不精确;迭代次数选择过大, 浪费运算时间。所以关于迭代次数的选择, 要根据监测出风口和进风口速度分量的收敛性后才能决定, 在迭代次数至3500次以后, 其计算残值已经在10-4以下, 误差值已经很稳定。因此在迭代次数的选择上, 可以用计算残值小于10-4为标准, 一般的计算分析时间约为为20小时。

图4为出风口处速度切片图。出风口处的速度分布, 在出风口出上部空气流速最快, 而在下部稍微靠上一点的位置, 流速比较慢。

进一步分析出风口的速度流场分部。在图5中, 分析了X轴正向的风速, 可以明显的看到在此出风口中部靠近上进风口处, 沿X轴正向风速分布有空洞, 有部分回流现象。而且附近的流速也比较低。

具体分析此出风口的流线图 (如图6) 。可以明显看到在出风口中部靠近上进风口处有回流现象。主要是由于在出风口开口过大, 后部的空气又被旋转的扇叶吸入进来, 从而表现出在后部有回流现象, 而且附近的流速比较慢。

4 结论

本文通过以上对轮式离心风扇出风口流场分析可以得出以下结论:在采用轮式叶片时, 由于叶轮的曲率比较大, 在出风口处流场分布不太均匀, 其风速最小的地方会出现在出风口中部靠上部分。在散热系统设计中要考虑到此问题。建议可以缩小叶轮的径向尺寸, 来减少回流现象。另外也可以变更叶轮的曲率形状来改变出风口流线, 进而改进次问题。

参考文献

[1]王福军.计算流体动力学分析-CFD软件原理与应用[M].北京:清华大学出版社, 2004.

[2]韩占忠, 王敬, 兰小平.FLUENT流体工程仿真计算实例与应用[M].北京:北京大学出版社, 2004.

[3]吴子牛.计算流体动力学基本原理[M].北京:科学出版社, 2001.

本文来自 360文秘网(www.360wenmi.com)，转载请保留网址和出处

【中国农业大学数值分析】相关文章：

中国波兰分析05-24