数值分析范文

导语：如何才能写好一篇数值分析，这就需要搜集整理更多的资料和文献，欢迎阅读由公务员之家整理的十篇范文，供你借鉴。

数值分析

篇1

关键词：地下水，渗流场，稳定性

1 引言

拟建基坑工程坑深19.5m，场地范围内地基土为软弱土，地下水丰富，水位埋深约为3m。本文基于对该基坑工程的初步设计，应用浸润线理论采用数值模拟分析的方法对基坑进行分部施工模拟，在开挖及降水的综合影响作用下计算分析基坑的应力、位移状态。

2 工程概况

基坑工程建设规模158127，拟建工程地下5层，地上42层。基坑周长315m，坑深19.5m。

3 工程地质、水文地质条件

3.1 地形地貌

拟建工程，建设场地地形平坦，相对高差不超过0.5m。归属于山前冲积平原地貌单元。

3.2 地层岩性

根据勘察钻探资料，建设场地表部为较厚的第四系覆盖层，详述如下：

3.3 水文地质条件

根据勘察资料，建设场地内地下水水位埋深2.5m～3.0m，主要为①杂填土层中的上层滞水与2、3粉质粘土层中的孔隙潜水。1杂填土层中的上层滞水水量较少，地层透水性较强，渗透性系数K=1m/d，2、3粉质粘土层中的孔隙潜水水量丰富，渗透性相对较弱，渗透性系数K=0.1m/d。

4 基坑支护结构设计

4.1 支护方案

根据建设基坑工程特征，工程场地地质条件。支护方案采用排桩+锚索，降水方案选用坑内积水明排。

4.2 计算参数的选取及设计结果

5 数值分析

数值分析中通过模型的建立及分步施工过程（开挖、降水）的模拟充分考虑基坑工程在开挖、降水综合作用下的应力与位移特征。

5.1 分步施工

将该基坑工程分解为12个分步工程：

1支护桩施工2基坑开挖至第1道锚索下0.5m，水位降至坑下0.5m3第1道锚索施工4基坑开挖至第2道锚索下0.5m，水位降至坑下0.5m5第2道锚索施工6基坑开挖至第3道锚索下0.5m，水位降至坑下0.5m7第3道锚索施工8基坑开挖至第4道锚索下0.5m，水位降至坑下0.5m9第4道锚索施工10基坑开挖至第5道锚索下0.5m，水位降至坑下0.5m11第5道锚索施工12基坑开挖至预定深度（19.5m），水位降至坑下0.5m。几何模型建立如图2。

5.2 地下水渗流分析

基坑场地范围内地层主要为粉质粘土，透水性相对较弱，基坑开挖过程中拟采用积水明排的降水方法，分步施工阶段基坑外侧的地下水不能完全的从坑壁排出，势必形成一定的水头差，在基坑及基坑外一定范围形成地下水渗流场。本文采用水位变化岸坡地下水位浸润的计算方法确定分步工程中降水后的地下水位线。然后将其赋入数值分析模型中，进行地下水渗流计算。

地下水位浸润线计算：，施工中降水后的坑内水位，施工中降水前的坑内水位，、降水前后计算点至坑壁的距离，降水后计算点的地下水位，降水前计算点的地下水位。

采用浸润线理论模拟基坑开挖降水过程计算孔隙水压力。降水后的水位线改变了基坑一定范围内原有的渗流场，孔隙水压力环境也随之改变（图3―5），桩后形成三角形的孔隙水压力（图6），最大孔隙水压力约为178kN/。

5.3 基坑应力位移分析

用数值分析基坑开挖降水后，在主被动土压力及孔隙水压力的综合作用下，计算坑壁的应力及变形状态（图7―10），坑壁应力呈三角形，最大应力约为417kN/，基坑外侧2倍的基坑深度范围内土体皆有不同程度的塑性变形，最大位移量约为130mm。

6 结论与建议

（1）基坑开挖降水作用改变了初始地下水环境，形成新的渗流场，最大孔压可达178kN/；

（2）开挖完成后，基坑坑壁应力呈三角形，最大应力约为417kN/，基坑外侧2倍的基坑深度范围内土体皆有不同程度的塑性变形，最大位移量约为130mm；

（3）建议对基坑外3倍深度范围的建构筑物进行位移监测。

参考文献

（1）张力霆，土力学与地基基础，高等教育出版社，2007年7月；

篇2

【关键词】数值分析科学计算信息与计算科学教学改革

【中图分类号】G642.0 【文献标识码】A 【文章编号】1009-9646(2008)09(a)-0044-02

科学计算是伴随计算机的出现而迅速发展并获得广泛应用的一门新兴交叉学科,是数学及计算机实现其在高科技领域应用的必不可少的纽带和工具。科学计算的方法和理论为科学研究与技术创新提供了新的重要手段和理论基础,科学计算现已成为与实验和理论并重的第三种科学研究方法。数值分析又称数值计算方法,它作为介绍科学计算的基础理论和数值方法的课程,对培养学生的科学计算能力,解决实际问题的能力以及创新能力起着非常重要的作用。数值分析现不仅是综合性大学信息与计算科学专业、应用数学专业的必修课程,也是信息类专业、工科类专业(如土木工程,航天航空,机械工程等)的必修课程,因此在进入21世纪后,随着计算机技术的飞速发展,社会对科学计算人才需求的变化,对数值分析课程进行教学改革,进一步提高数值分析的教学质量和教学效果,以适应新形势下社会对人才的需求,变得尤为迫切。

1 课程的内容和特点

1.1 课程的内容

数值分析课程是以微积分、线性代数、常微分方程等课程的基本内容为基础,以算法设计为手段,以计算机为工具实现算法,全面地介绍了科学研究及生产实践中各领域所遇到的普遍的数学问题的数值方法和理论,因此它是所有与科学计算密切相关的专业的一门重要的基础课程。

数值分析课程的基本概念有:误差、收敛性、稳定性、插值、计算量、存储量等,这些概念是用来刻画数值方法的准确性、可靠性、效率以及计算复杂度;核心内容包括:代数问题数值计算、函数的数值逼近论、计算几何、微分方程数值解、最优化方法、统计计算等方面[1]。

1.2 课程的特点

数值分析作为大学本科阶段的一门数学基础课程,它具有下面的一些主要特点:

1.2.1 实践性强

上文已提到数值分析的核心内容是研究应用计算机对科学研究和实际应用中常见的数学问题进行数值求解的各种方法。这些方法除了理论上要正确可行外,还要求通过数值试验证明是有效的、实用的。因此,与其它数学类基础课程有所不同,数值分析课程很强调实践,要求学生在学习了每个算法后用学过的计算机语言编程或借助于一些数学软件,如Mathematica、Matlab,完成算法的计算机实现,这样学生不仅知道理论上是如何计算,而且还能把结果真正地计算出来。有些算法尽管在理论上不够严谨,但通过实际计算、对比分析等手段证明是行之有效的,在实际应用中也是可以采用的。

1.2.2 应用性强

数值分析中的理论与算法并不是都来源于纯粹的微积分、线性代数、常微分方程等课程,很大程度也是由于科学研究和实际工程计算的需要,可以说它来源于实际又应用于实际。一方面,目前很多算法已经在物理力学、生物信息、计算机应用、土木工程、航天航空、石油勘探、环境工程、材料等领域得到了广泛的应用;另一方面,数值方法在这些领域的应用过程中会产生一些新问题,为解决这些新问题,研究人员和技术人员深入研究、不断攻关,从而又会促进数值分析理论和算法的发展。

1.2.3 计算公式多且复杂,不好记忆

数值分析课程的许多算法或是采用“构造性”的方法,把计算公式具体构造出来,如Lagrange插值,Hermite插值等;或是采用“递推”方法,将复杂的计算过程化简为简单计算过程的多次重复,例如:求解线性方程组的Jacobi迭代法,Gauss-Seidel迭代法等,从而使得算法里出现的计算公式多且复杂,学生不好记忆。

2 数值分析课程教学中存在的问题

目前,在数值分析课程教学中普遍存在下面几个问题:

2.1 内容多,学时少

数值分析的内容非常丰富,包含代数问题数值计算、函数的数值逼近论、计算几何、微分方程数值解等,但学时数却不多,一般只有64/54,教师使用传统的教学方法讲授这些知识时,由于计算公式多且复杂,推导繁琐,不得不采用填鸭式教学,每节课都在讲新课,拼命赶进度,这种做法既忽视了学生这个主体的作用,使学生处于被动学习的状态,时间久了就会失去学习数值分析的兴趣,产生厌学情绪;同时也扼杀了学生学习的主动性、积极性和创造性。

2.2 重理论,轻实践

前面已提到数值分析是一门与计算机关系密切的实践性很强的课程,但在传统的教学模式中,教师通常只注重讲授算法原理、误差分析和收敛性证明等理论内容,忽略了上机实践环节。学生动手实践机会少,就会使学生不能彻底理解算法的原理与运用,只知道大概是“怎么算”,但不会真正把结果算出来,造成理论与实践的脱节;同时也会使学生感到数值分析的学习很抽象、枯燥、难学,以致学习兴趣不高。

2.3 缺乏带有上机实验内容的教材

数值分析课程是20世纪50年代末在我国一些综合性大学开始开设的,由于当时计算机的应用在我国刚开始起步,所以当时使用的数值分析教材基本上都没有上机实验内容的。人类社会进入21世纪后,计算机技术发展迅猛,许多复杂的计算能以惊人的速度在计算机上得到满意的计算结果,但作为介绍科学计算理论和数值方法的数值分析课程,其大部分教材没能适应时代的变化,还是沿用旧教材的内容结构,没有增加上机实验内容,这给教学带来很大的困难,直接影响了人才的培养。

3 课程教学改革

为全面贯彻落实科学发展观,切实把高等教育重点放在提高质量上,教育部和财政部于2007年先后颁布了《教育部财政部关于实施高等学校本科教学质量与教学改革工程的意见》和《教育部关于进一步深化本科教学改革全面提高教学质量的若干意见》。在这些纲领性文件指导下,全国高校掀起了一轮新的教学改革活动,数值分析课程的教学改革也开展得如火如荼并已取得一定成效(见文[2,3])。很多专业开设了数值分析课程,不同专业有不同的要求,有些要求相差很大,下面谈谈我校信息与计算科学专业的数值分析课程教学改革的一些措施和想法。

我校信息与计算科学专业是2006年开始招生。2006级学生的数值分析课程安排在第四学期。为了做好数值分析课程的建设工作,我们专门成立了数值分析课程建设小组。以下就是我们在数值分析课程教学改革的一些措施和想法:

3.1 选择合适的教材

对于信息与计算科学专业的学生来说数值分析是一门专业基础课,要求比其它专业高,既要掌握算法的应用,也要掌握算法的原理、误差分析、收敛性分析等理论知识,因此选取的教材应兼顾这两方面。我们采用了林成森新编的《数值分析》4. 该教材内容全面,阐述严谨、详细、深入浅出,且每个常用的算法都给出了伪程序。我们知道很多学生都怕编程,选择了该教材将有助于学生动手编程上机实现算法。

3.2 内容选取合理,突出重点

信息与计算科学专业的数值分析课程的课时是64学时,在这有限的学时内要详细地讲授所有的内容是不现实的,也是不可能的。我们必须对教学内容重新优化设计,要有侧重,以主带次。在课堂上有限的时间里力求讲明白一个主题,重点讲授典型的、具有代表性的并能体现其思想方法的常用算法和理论,而对那些原理相近的内容则可留给学生课后自学,使课堂讲授的内容真正做到“少而精”。例如:在讲授解线性方程组的Gauss消去法时,我们在课堂上重点详细地介绍了Gauss消去法和Gauss列主元消去法,学生在理解了这些方法的原理之后,再结合高等代数知识就很容易自学Gauss按比例列主元消去法和Gauss-Jordan消去法。

3.3 加强实践环节,提高科学计算能力

数值分析有别于其它数学类基础课程,它是一门与计算机结合密切的实践性很强的课程。由于课时紧张,在2006级的培养计划里数值分析课程没有安排上机实验课,但我们任课老师在讲授完每一章后都会布置一些上机实验题让学生自己编程上机运算(注:学生在第二学期已经学了C语言),并且要求把程序和实验结果交给老师批改。通过上机实验,使学生加强了对课堂内容的理解,同时又和计算机语言结合起来,自己编程,动手计算,使学生较好地掌握常用的科学计算方法和技巧,积累了解决实际问题的能力,而且也培养了学生质疑问题的能力和实践创新能力。这种做法很受学生欢迎。为了进一步提高上机实验的质量,我们打算在新的培养计划中增加9～12课时的上机实验课。

3.4 采用现代化教学手段

由于数值分析内容丰富,计算公式多且复杂,推导过程繁琐,所以我们采用了多媒体教学。这样可节省板书、画图的大量时间,教师有充足的时间把基本概念,算法构造的原理等重点内容讲透彻,而且在多媒体教学过程中还可以结合运用Mathematics,Matlab等数学软件演示算法的计算过程和结果,使学生加深对抽象知识的直观理解,从而提高了学生的学习兴趣和学习积极性。

3.5 积极编写新教材

由于在新的形势下数值分析课程的教育观念、教学要求以及教学内容已经发生了变化,所以教材的编写也应跟上时展的步伐。现虽然已有一些数值分析的新教材尝试了与计算机应用相结合,融进了Basic语言、C语言或Matlab数学软件编程实现算法的新内容,但大家认同的精品教材不多,因此仍需从事数值分析教学的同行们继续努力、不断改革、勇于创新,编写适合不同专业需求的新教材。由于合适的上机实验辅导教材奇缺,所以我们的课程建设小组已开始编写上机实验辅导教材,以解燃眉之急。

3.6 改革考核方法

考核是评估教学质量和学习效果的重要手段。数值分析传统的考试方法为笔试闭卷考试,主要是考察学生的识记能力,忽略了对学生编程上机实践的要求。在新的形式下,为提高大学生的应用能力和创新能力,我们对数值分析的考核方法进行了改革,把上机实验列入考核内容,即学生的期评成绩由笔试成绩、上机实验和平时成绩按比例确定。

4结语

我校信息与计算科学专业的数值分析课程建设起步晚,课程的教学改革仍需我们不断的探索与实践,我们也将继续努力,并借鉴兄弟院校该课程的教学改革经验,朝着精品课程目标建设数值分析课程。

参考文献

[1] 伍渝江,尤传华,丁方允.数值分析课程的继承与改革 [J].高等理科教育,2000,8(1): 46-49.

[2] 杜廷松.关于数值分析课程教学改革研究的综述和思考[J].大学数学,2007,23(2):8-15.

篇3

关键词：MATLAB；数值分析；教学改革

中图分类号：G421 文献标识码：A文章编号：1007-9599 (2011) 13-0000-01

Attempts to Strengthen the Teaching of Numerical Analysis by Matlab

Huang Cheng

(Foreign College Of Central South University of Forestry and Technology,Changsha410201,China)

Abstract:The combination of Matlab and Numerical Methods is a viable reform of teaching about Numerical Methods.This paper introduces the new teaching’s system and new curriculum.In order to solve the defects of curriculum,provides new ideas and approaches.

Keywords:MATLAB;Numerical Analysis;Teaching Reform

历年本科学生对数值分析的反映可表达为：抽象、冗繁、枯燥，只是为了考试而学。这种状况必须改变。

MATLAB在上世纪80年代推出后，已发展为一种功能全面的软件。它问世后，一个“用软件工具增强线性代数教学”的项目ATLAST在美国国家科学基金立项，并在90年代获得成功。其后欧洲日本等国家快速跟进，并在相邻学科中进一步拓展。

09年中南林业科技大学将“数值分析”改为“数值分析与MATLAB”，要求在教学领域有所创造和前进。我受教研究室的支持和帮助，在教学中进行了初步尝试。本文讨论该课程教学体系的改革和教学内容改变两个问题。

一、教学体系的变革

国内“数值分析”的多数是沿着方法原理、计算步骤、计算框图、计算速度、误差、收敛和实例这样一个体系展开A。这套体系成型于手工计算时代。在低级语言编程的时代作了一定改进，但仍显繁琐。目前在MATLAB下，各种方法已打包为指令，不再需要对各种方法进行编程。使用者只需知道方法原理、调用方法和指令即可。导致教学的立足点己经不同，这是一个根本性的改变。面对这样情况的出现，课时与内容矛盾不断加剧，国内计算器时代建立起来的教学体系到了必须改变的时候，也具备了改变的条件。

现在我们在“数值分析与MATLAB”课中讲授数值分析是按：问题提法、解题理念与原理、解题指令与多指令的配合、计算结果表示、实例和实验这样一套教学体系进行的。

这样有什么理由和好处呢？一个课程的教学体系与该学科体系、教学目的、学习者层次等因素有关。从教学目的、学习者层次角度而论，本科生学习“数值分析”课的目的是应用它来解工程和科学问题；从学科体系看，数值分析的学科发展与计算工具有着密切的关系[1]。可惜，目前教学领域未适应这种发展和变革。所以改革是很有必要的。MATLAB只是软件发展的一个阶段，但也引起了很多变化。它的影响表现为两个方面：一可以把传统分析方法、设计程式处理得更简捷；二是推动新方法、新程式面世。例如在MATLAB中，用4阶Runge-Kutta数值积分法解常微分方程，核心部份是[2]

〔t，Y〕=ode45（odefun，tspan，y0）

从应用角度看，详细的讲ode45内Runge-Kutta算法是不必要的。只需知道输入量的意义、输入方法，能看懂输出量t，Y的意义，就可以了。换句话说，知道解题指令与多指令的配合、计算结果表示等就可以用了。但为保证正确使用，就必须了解其适用范围、优缺点等。讲授时就不必追求数学上的公式推导，而应着力阐述概念。其中包括：问题的数学提法，解题方法的实质理念、几何概念和计算公式。而有关计算速度，因仅用计算次数来分析计算速度，显得不适用。而在MATLAB中有计时器可测出指令计算问题的时间，从而得出该方法是否接受；也可对比不同方法的速度，从而给出评价。所以分三条途径学习：一是在解线性方程组的章节中重点介绍用计时器进行实测和评价；二是以做习题的形式，要求用计时器对比各种方法的计算速度；三是在章节总结中，从原理和实算二方面进行归纳对比。有关误差分为二条线：一是把误差概念集中一章介绍；二是在讲迭代法解微分方程和数值积分的指令时，介绍误差控制和自适应性。这样的好处是：立足应用，原理概念清晰，联系实际。一定程度上克服了抽象、枯燥的问题，提高了学生实际使用能力。

二、教学内容改变

在新形势下如何规范精选教学内容是个急待解决的问题。我们的原则是：摒弃部分原有的分析方法和设计程式，传授以MATLAB为基础的分析方法和设计程式。将课程内容设定为三部份：一是MATLAB基础，含MATLAB语言基础和数据可视化基础；二是多项式运算、数值微积、插值逼近、误差分析、解微分方程等传统数值分析的内容；三是符号运算、线性系统分析和Simulink建模仿真介绍。较传统意义上的数值分析内容更为广泛。

课程中MATLAB语言基础是最基本的。同时为加强形象思维和图示计算结果，安排了数据可视化基础。鉴于多项式运算在线性系统分析中的作用；插值逼近、解微分方程等为传统数值分析中最基本内容，于以保留。数值微积在解微分方程中要用到，也选取了[3]。符号运算、线性系统分析和Simulink等应用面最广泛，作为高端应用代表入选。

上述教学内容前两部份数学基础主要为线性代数、高等数学。线性系统分析和Simulink建模仿真等，则涉及运算微积、控制论、仿真等专业知识。各专业学生对于此类基础概念参差不齐是个突出难题。此部分课内适当补充基础，教师掌握好课程深浅是成败的关键。因此对数值分析实行分类分级教学，不同专业和不同层次的学生采取不同的教学方式，从总体上提高教学质量是特别必要的。

通过教学实践，我们形成了新教学体系。虽然这套体系尚不完整，但为解决“数值分析”课的蔽端，提供了一种新思路。但所选方法是否抓住了最基本、最有用的方法？能否反映学科前沿？颇值得讨论。改革是个长期的过程，成效多大，更要通过社会与教学的长期考验才会有结论。写了上面的文字，希望能得到大家的指正。

参考文献：

[1]马昌风,林伟川.现代数值计算方法[M].北京:科学出版社,2008

[2]张志涌,杨祖樱.MATLAB教程(R2010a)[M].北京:北京肮空航天大学出版社,2010

篇4

关键词:数值分析;课程改革;开放式教学

中图分类号:G64文献标志码:A文章编号:1673-291X(2010)31-0293-02

引言

大学毕业生是国家宝贵的人才资源,其就业问题关系到国家经济建设、社会稳定和人民群众的根本利益,关系到高等教育的持续健康协调发展。高等学校作为国家和社会培养高素质人力资源的主要阵地,要以市场需求为导向设置专业,变革教学方法,提高教育质量,培养能够适应社会需要的复合型人才。

数值分析是一门理论与实践相结合的课程,它是信息与计算科学专业学生的主干课程,我们在强调它的理论结构时,更注重它的实用价值。同时它也是各种计算性科学的共性基础与联系纽带,实践证明数值分析作为科学计算的基础与核心,已被广泛应用于科学技术和国民经济的各个领域。如计算物理学、计算经济学、天气预报、大型水利工程的建设与设计等。如何根据市场需求改革数值分析课程的教学方法?如何提高教育质量,培养具有一定层次信息技术素养的大学生?如何提高学生的实际应用能力和创新能力?这些都是我们值得探讨的问题。

一、传统教学方法的不足

传统教法是以教师的讲授为主,课程记忆性内容太多。而课堂例题、课后作业和实验题目又都是固定的、单一知识为主的题目,程序式的问题较多,问题的模仿性太大,高层次、开放性和思维性问题较少。传统教法中主要以教材为主线,授课中多以演示法及互动法为主,并配备固定的实验和作业。这样教学虽然也能达到一定的效果,但缺点也很显然。首先,教学以理论为主线,实践性问题过少,学生只知道理论,不知道如何应用。其次,传统教法不能够体现出特优生、一般生和差生特点,更不能满足到特优生想学多、学好和学深的要求。传统教学方法忽略了学生的个性发展和对学生创造性思维的培养,不利于培养高层次的、可就业的IT人才。

二、课程改革的具体方案

为了弥补这些不足,在授课中应少讲理论,多讲实践应用,让学生在动手编制程序的过程中理会理论知识点,通过改正程序中出现的错误,细化知识,达到强化理论知识、精通程序编制的目的。其中有两个重要的环节,一是作业题目、实验题目的合理设计。让有能力的学生做一些有深度、有广度、有拓展的综合性题目,让差生做一些可独立完成的单一知识点题目,分层次、分需求教学,分层次、分需求管理,这样才能充分体现学生的自主学习能力和创新能力。二是教师要按学生的反馈信息、学生的需求合理修改、设置教学实践例题,引导学生解决复杂问题,让学生感受到问题的不同解法、程序编制的不同途径、算法的优劣等方面给编制程序带来的不同效果,影响学生形成好的编程习惯。我们的具体方案如下。

1.开放式作业

数值分析是一门集理论、抽象和设计于一身的计算机科学与技术专业的重要基础课程,它主要研究用计算机解决数学问题的数值方法及其理论。数值分析既有纯数学高度抽象性与严密科学性的特点,又有应用的广泛性与实际实验的高度技术性的特点,是一门与计算机使用密切结合的实用性很强的数学课程。通过本课程的学习,能使学生熟练掌握各种常用的数值算法的构造原理和过程分析,提高算法设计和理论分析能力,并且能够根据实际问题建立数学模型,然后提出相应的数值计算方法,并能编出程序在计算机上算出结果。因此,在实际教学中我们采取开放式作业。

开放作业包含两个方面的内容:一方面是作业内容的开放。教师依据教学要求,设计10道左右的作业题目,分为综合性题目、单一知识点难题和单一知识点基础题等类型,其中综合性题目主要以工程上的实际问题为背景,要求学生建立相应的数学模型并求解。这些题目由学生自主选择组合完成其中的2～3道题,每题有相应的难度系数,作为衡量作业质量的标准。另一方面是作业提交形式的开放。学生的作业可以写在作业本上上交,也可以提交电子作业。电子作业需要学生把编制的程序运行无错后提交,这种形式锻炼了学生的动手能力、提高了学生解决实际问题的能力和纠错的能力。作业成绩由选题的难度系统和选题的个数决定。

2.开放式实验

数值分析是数学的一个分支,但它不像纯数学那样只研究数学本身的理论,而是把理论与计算紧密结合,着重研究数学问题的数值方法及理论,并将这一理论在实践中应用,以达到解决实际问题的目的。为了让学生更好地理解他们所学的知识在实际生活中的应用,我们必须重视实验教学环节,在讲述基本原理和基本操作方法的同时,应注重培养学生的创新能力和工程实践能力,提高他们的程序设计能力和上机操作那个能力。因此在实际教学中我们采取开放式实验。

开放实验包含两个方面的内容:一方面是题型的开放。实验题型分为程序纠错题、程序填空题、程序验证题和程序编制题四类。另一方面是选题的开放。实验题目有15道左右,难度、深度、广度各不相同,学生自主选题组合,题目不固定,个数不固定。选什么程度的题目,选几道题完成由学生全面决定。老师唯一约束的是学生要绝对独立完成,教师可根据学生提交的实验程序文件验证。学生的实验成绩由选题的难度系统和选题的个数决定。

3.开放式教学

数值分析课介绍的数学方法大多数都有工程背景,对离散数据插值和对函数的数值逼近的方法来源于对实验数据处理、产品外形设计等工程实际问题,样条函数方法正是图像处理技术等方面的关键部分,FFT 技术能够快速处理数据,在机械、电子、信息、自动化工程中的实时信号处理中占有重要位置。因此,在实际教学过程中我们可以根据工程背景采取开放式教学。

开放教学包含两个方面的内容:一方面是教学实例的设计。传统教学中都以教课书中实例为主,教师在上面讲,学生在下面边听边看书上的程序,对照起多媒体投影中的程序和书本的程序一模一样,学生对教师的演示也就提不起兴趣了,感觉既枯燥又乏味。鉴于此,我们可以开放教学实例内容,由教师或学生提供实例题目,课上教师和同学共同分析问题、提供实例的数值解法,并给出相应的实现程序。教师可通过这种实例分析,引导学生解决多知识点综合性题目。例如,在介绍三次样条插值方法时,可以给出工程上飞机机翼型值线的实际问题,通过列举机翼上缘轮线的数据,要求同学们自己编程用三次样条函数画出机翼曲线。这样,不仅锻炼了学生的动手能力,而且还能进一步加深学生对数值分析这门课的学习兴趣。另一方面是基础知识的实践性教学。授课时教师精讲基本知识点,做适当的演示,然后留有一定的时间给出一些基础问题,让学生扮演教师角色,讲解具体的数值算法、编制程序,其他的学生可对程序进行优化、纠错。学生由被动的听讲变成了主动的讲解,角色的转换激发了学生的热情,吸引了学生的注意力。学生在讲的过程中,也对所讲内容有了更深层的认识。设计开放教学着重强化学生自我获取知识的能力、组织和交流能力。

4.充分利用网络资源,开放师生之间的交流

为了减少课程教学难度,提高教学质量,加强对学生自主学习的辅导力度,扩大辅导面,我们充分利用网络资源,开展了多种形式的导学、网上讨论咨询答疑、学习辅导活动,具体方式有以下几种:(1)网上讨论、咨询答疑。这种类型的咨询答疑活动是利用网络环境,由专业教师主持,通过BBS展开师生间、同学间的讨论、信息交流,适宜于进行学习方法讨论,教学信息的交流,简单课程学习问题的咨询答疑。(2)电话咨询答疑。电话咨询答疑,一般来说仅是单一的声音媒体交流,因此多用于简短信息的交流和简单概念性问题的咨询答疑,例如,确定考试或辅导的时间、地点等等。(3)E-mail信箱答疑辅导。E-mail信箱答疑辅导时纯文字媒体的交流,常常用于解题示范,例题展示,信息交流等类型的问题答疑辅导。(4)运用现代信息技术。运用现代信息技术编写并制作CAI教学课件和电子教案,在课件中的每一章明确本章的学习要求,课后有复习思考题。数值分析课程的教学大纲、试卷、试卷分析、标准答案也一律做成电子版的,以方便学生的查阅。

结语

数值分析是一门理论与实践相结合的课程,它是信息与计算科学专业学生的主干课程。目前,用计算机进行科学与工程问题的科学计算,已成为与理论分析,科学实验同样重要的科学研究方法。在数值分析课程教学中,从内容编排到课堂教学,我们都本着提高学生应用能力的思想,简化理论推导,重视实践,极大增强了学生利用数值计算方法解决实际问题的信心,同时也提高了学生的市场竞争力,为学生更好的就业提供了保障。当然,数值分析课程的教学改革任务是长期的、艰巨的,以上的工作只是我们的一点尝试。在今后的教学中,我们要时刻总结良好的教学经验,为数值分析课程的教学改革多作贡献,我们将广采百家之长,不断地改进我们的教学方式方法,争取取得更好的教学效果,培养更多的合格人才。

参考文献:

[1]王能超.数值分析简明教程[M].北京:高等教育出版社,1984.

[2]姚传义.面向应用提高数值分析课程教学效果[J].化工高等教育,2007,(2): 39-41.

篇5

关键词分解槽； Fluent；搅拌桨叶；数值分析

中图分类号TQ13 文献标识码A 文章编号 1674-6708（2014）110-0048-02

Numerical Simulation of optimum designing for Stirring blade of Precipitator Tank

Wang You

Guiyang Aluminum-Magnesium Design & Research Institute Co.Ltd.， Guiyang，Guizhou， China 550081

AbstractBased on the principle of hydromechanics similarity， this paper gives a numerical simulation analysis on the precipitator’s stirring blade （MIG）relevant design modification， and combined with the fluid analysis software Fluent. The paper competitively analyzes four aspects as the three dimensional flow field velocity distribution， solid content difference analysis， stirring power and the maximum shear stress， provides reference basis for design of stirring blade.

KeywordsPrecipitator Tank， Fluent， Stirring Blade， Numerical Simulation

0 引言

随着氧化铝生产大型化的发展，传统的Φ14m分解槽已不能满足生产要求，需要开发更大直径型的分解槽。分解槽大型化设计的主要难点是搅拌装置的设计，其搅拌在生产过程中既要满足料浆充分的混合悬浮又不破坏晶种的长大，因而其搅拌有一定的特殊性。搅拌装置设计的重点在于桨叶的选型，目前由于搅拌过程种类繁多，介质情况差异很大，实际使用的搅拌桨叶形式多种多样。目前的选型方法多是根据实践经验，选择习惯应用的桨型，再在常用范围内决定搅拌器的各种参数。也有通过小型试验，再进行放大的设计方法。随着计算流体力学的发展，运用流体分析软件对搅拌过程进行数值模拟技术已日趋成熟，本文就是在现有的氧化铝生产上通用的MIG型搅拌器的基础上，运用相似原理和Fluent软件提供的稳态多重参考系法(MFR)对设计的三种搅拌器进行数值模拟，并与原有Φ14m分解槽的MIG型搅拌器进行对比分析，得到适合大型分解槽搅拌使用要求的桨叶形式，为设备改进优化提供设计参考依据。

1 研究对象及模型建立

1.1 物理模型

分解槽整体模型如图1，槽体直径16m，高42m，内设置6层搅拌桨、1组挡板、在挡板对面设置提料管，建模中忽略提料管内部流场，忽略搅拌桨厚度。三种不同桨叶结构形式见图2，其中模型A是传统MIG桨叶形式，桨叶与轴的夹角为60°，模型B是将模型A中内桨叶分为上内桨及下内桨两部分，模型C是将模型A中桨叶与轴的夹角由60°增加到70°，具体结构尺寸见表1。

图1 整体搅拌分解槽模型

图2 桨叶模型图

模型A 模型B 模型C

分解槽内径（m） 16 16 16

液面高度（m） 38 38 38

桨叶层数 6 6 6

桨叶直径（m）底层 11.6 11.6 11.6

其它层 10 10 10

每层桨叶之间高度（m） 6 6 6

轴径规格（mm） Φ610X26 Φ610X26 Φ610X26

桨叶与轴夹角（°） 60 60 70

内浆分段 1段 2段 1段

挡板数量（含出料管） 2 2 2

转速（rpm） 4.4 4.4 4.4

表1 分解槽不同搅拌桨叶形式结构尺寸

1.2计算方法

本文选用Realizable k-ξ湍流模型，欧拉-欧拉多相流模型对分解槽内固液体系进行数值模拟。在模型中考虑相间作用力、虚拟质量力及升力对固体颗粒的影响，其中固－液两相间阻力系数的理论计算采用相间相互碰撞的Gidaspow 模型。

采用稳态多重参考系法（MRF），将各个计算区域分成两个或多个互不重叠的圆筒状区域，整个分解槽分为旋转区域和静止区域两部分，旋转区域的几何结构只有搅拌桨，静止区域的几何结构包括整个槽壁、挡板与提料管，旋转区域创建旋转坐标系，静止区域创建静止坐标系，搅拌桨相对内部子区域静止，实现搅拌桨的旋转。

1.3 工艺条件

表2是分解槽实际生产中的一组常用物性参数。

项目料液的密度

kg/m3 料液的粘度

Pa.s 颗粒密度

kg/m3 固含

g/l 颗粒

大小

μm 含量

（质量分率）%

数值 1753 0.0038 2424 1000 70 10.25/31

表 2 物性参数

2 模拟结果分析

2.1 分解槽内物料的三维速度矢量

图3为穿过搅拌桨叶中心X-Y平面的三维速度分布图，模型A、B形成的流场相似，都只在每层桨叶之间形成了非常明显的流体循环，流体在槽内基本是在每层桨叶之间流动，没有形成桨叶之间的两层流体循环，而模型C在每层浆之间形成明显的桨间循环，内外桨叶有明显地流体向上运动之后分别向内外桨叶的流场位置循环从而形成了明显的两层流体循环，导致颗粒在槽内提留时间要比模型A、B长，从而有利于颗粒的结晶长大，也同实际设计MIG型搅拌器的预期效果吻合。

图3 流场(X-Y平面)三维速度分布图

2.2流场内的均匀度

分解槽搅拌的主要目的之一还要保持溶液浓度均匀，保证晶种与溶液有良好的接触以利于析出晶体。通过模拟可以得到颗粒相在整个流场中的分布状况，以及确定颗粒相的高浓度区域。

图4给出了70μm颗粒的体积相分布情况，从图中可以看出，在分解槽底面上有比较明显的沉积，说明底层桨附近区域是沉积高危区，且易沉积的区基本可以分为两块，就是搅拌轴附近区域以及槽底边缘的区域。模型A、B的沉积区域明显多于模型C。

图4 70μm颗粒体积相分布图（X-Y平面）

流场内的最大固含差，可以在一定程度上反映出整个搅拌的颗粒相分布的均匀程度，本文根据固体颗粒体积分数换算为固含量，进而得到固含差，表 3给出了三种桨叶形式的最大固含差的计算分析值。

从表中可以看出，模型C的最大固含差最小，模型A最大，工业生产要求固含差控制在5%∽8%以内，从计算结果看，模型B和C可以满足。

模型颗粒

直径体积分数/% 固含/g/l 最大固含差/%

最大最小差值最大最小差值

A 70μm 31.38 28.93 2.53 1033.11 918.21 114.9 11.12

125μm 11.24 8.95 0.60

B 70μm 31.63 30.5 1.13 1046.2 968.4 77.8 7.45

125μm 11.53 9.45 2.08

C 70μm 31.39 30.43 0.68 1035 977.11 57.89 5.59

125μm 11.31 9.88 0.44

表3 三种桨叶形式的最大固含差

2.3搅拌功率

搅拌功率是搅拌中重要的参数，一定程度影响了生产成本和工业生产的现实可能性。

图5给出了运用Fluent计算的三种桨叶形式各层桨叶消耗功率分布情况。模型A消耗的总功率为106.4 KW，模型B消耗的总功率为137.1 KW，模型C消耗的总功率为115.6 KW，通过比较分析，在满足使用要求和经济性方面综合考虑，模型C的综合性能最好。

图5 功率分布图

3 结论

1）本文建立了大型分解槽搅拌桨叶的三种计算模型，并采用稳态多重参考系法对三种桨叶的搅拌过程进行了数值模拟计算，结论是模型C相较于模型A和模型B，搅拌流动效果较好，沉积区最少，均匀度最好，综合性能经济指标亦能满足生产需要；

2）通过与现有工业上使用的分解槽及其搅拌结构进行对比分析，运用Fluent计算所得的分解槽搅拌模型能满足实际生产对分解槽搅拌结构和工艺性能的要求，能为分解槽的大型化工业生产提供可靠的理论设计依据。

参考文献

[1]王凯，虞军，等.搅拌设备[M].北京：化学工业出版社，2003.

[2]钟丽，黄雄斌，贾志刚.用CFD研究搅拌器的功率曲线[J]. 北京化工大学学报，2003，30(5):5-8.

[3]李振花，何珊珊，万茂荣，谈遒.搅拌槽中的流体力学模型[J].高校化学工程学报，1996，10(1):22-29.

篇6

[关键词]雷电；接闪；电场数值分析

中图分类号：P427.32 文献标识码：A 文章编号：1009-914X（2014）47-0260-01

一、雷电

雷电是发生在大气层中的声、光、电等物理现象，通常认为是由于热空气上升，冷空气下降过程中的热交换，产生带有正负电荷的小水滴积聚形成积雨云，在积雨云（雷云）形成过程中，由于大气电场以及温差起电效应、破碎起电效应的同时作用，正负电荷分别在云的不同部位积聚。当这些电荷积聚到一定程度时，就会在云与云之间或云与地之间发生放电，也就是人们平常所说的“闪电”。其中云与地之间产生的放电现象又称为地闪，极容易对人类造成不可挽救的危害，也是我们进行雷电防护研究的主要对象。

雷电放电瞬间可产生数十千安，甚至数百千安的放电电流。雷电流能产生巨大的破坏力和很强的电磁干扰，给人类的生活、工作带来很大的影响，它引起的灾害是自然界十大灾害之一。

随着科学技术的不断进步，各类电子信息产品得到广泛应用，特别是电子信息系统的应用，极大的方便了人们的生活。但是，这些电子设备普遍存在着绝缘强度低、过电压和过电流耐受能力差、对电磁干扰敏感等弱点，一旦建筑物受到直接雷击或其附近区域发生雷击，雷电过电压、过电流和脉冲电磁场会通过供电线、通信线、接收天线、金属管道和空间辐射等途径侵入建筑物内，威胁室内电子设备的正常工作和安全运行。如防护不当，这些雷害轻则使电子设备误动作，重则造成电子设备永久性损坏，严重时还可能造成人员伤亡。随着全社会现代化水平的不断提高，雷电对电子和通讯等设施的破坏，而造成的经济损失及人员伤亡，远远超过了雷击火灾的损失，雷电灾害已成为“电子化时代的一大公害”

雷电灾害造成的损失的大小是牵扯到社会许多方面的十分复杂的问题。雷电能造成人员伤亡，即使得建筑物起火、击毁、能对电力电话、计算机及其网络等设备造成破坏，雷电又是年年重复发生的自然现象，在每年的七八月份是雷暴的高发期，尤其是在热带地区，雷电次数就更多。但是现在对于雷电灾害及防雷知识，全社会缺少一定的认识，加上侥幸心理，所以极易造成雷击事故。因此，必须进行深入的研究和采取必要的措施。

就减轻雷电造成的损害而言通常考虑的措施有：一是加强预报工作，在雷暴来之前就能像预报天气一样，让人们做好准备，如拔掉电源插座等。第二就是加强新建建筑物的防雷击的能力，在建筑物建成之前做好建筑物的防雷设计和施工。目前国家气象部门很重视这个方面的工作，各个地区的气象部门也都开展了这样的业务。三是加强雷击时的救援工作。所以目前普遍采用的是加强防雷工程的做法。

当雷电击中建筑物时，由于雷电是具有高电压、大电流，作用时间极短的瞬变过程，通常在瞬间释放出巨大的能量，把被击中金属熔化，使物体水份受热膨胀，产生强大的机械力，或分解成氢气和氧气，产生爆炸，使建筑物遭到破坏。雷击产生的高温引起建筑物燃烧构成火灾，产生的高压引起触电。根据目前的防雷理论，无论采取哪种保护方法，都需要使用接闪器接闪，通过引下线将雷电流引下至接地装置，由接地装置散入大地中。在此过程中存在以下雷击安全隐患：

1）雷电流沿引下线传导过程中，在其周围存在很强的电磁场，可能引起感应过电压和过电流。

2）雷电流由散流装置入地过程中形成的电位梯度过大会导致行人因跨步电压而发生人身伤亡事故。

3）直接雷击时，雷电流在泄放和散流过程中因电阻压降和电感压降导致高电位通过静电感应在水平布设的信号线路和电源线路上产生的过电压损坏设备接口，并有可能导致反击及人身触电伤亡事故。

二、感应雷的灾害分析

1）散流时引起的过流（压）损坏：当雷电击中建筑物散流时，分流到配电系统、信号线路、其它金属管道中的雷电流引起设备过压（流）损坏或人身触电导致伤亡事故。

2）发生直接雷击，雷电流泄放时，建筑物内部分布着暂态电磁场，尤其以引下线周围最为强烈。此电磁场将会对建筑物内各个系统产生作用，引起设备误动作或损坏。

3）室内暂态磁场作用在信息系统环路上，将会产生感应过电压（流），导致设备接口或设备本身损坏。

4）雷雨云（积雨云）引起的感应雷击而发生损坏。当有雷雨云经过沿线上空或附近时，由于静电感应会在电源线路、信号线路、控制线路上感应出极性相反的静电荷，当雷云放电后，这些静电荷由于不能及时入地会产生过电压（流）损坏设备。

5）云内闪和云际闪对信息系统设备的影响。云内闪和云际闪产生的雷电电磁脉冲（LEMP）可引起内部设备因感应过电压（流）损坏。

三、接闪杆接闪瞬间电场数值分析

图1为本文研究对象：防护直击雷的接闪杆。图一

如图2所示，由于接闪杆为轴对称模型，所以本文只计算了其一个面。

图3为ANSYS环境下对接闪杆进行1：1建模后，对其施加kV电压，其X轴方向的电场强度分布图。由图3可知，其X轴方向的电场强度最大值及最小值均在接闪杆尖端。

篇7

关键词:数值分析；MC模型；HS模型

PLAXIS deformation in the pit in the application of numerical analysis

Ye Haibo

Jiading District, Shanghai Real Estate (Group) Co., Ltd.

Abstract: The Excavation of a numerical analysis of the key issues is the use of appropriate soil constitutive model. [1] PLAXIS geotechnical finite element analysis software is used to solve geotechnical deformation, stability and common issues such as groundwater seepage finite element family of software, which provides a mole - Coulomb (MC), soil hardening (HS ), soft soil creep (SSC) and other soil material model. This departure from the Project to discuss the numerical analysis of pit deformation PLAXIS for calculating the parameters of thinking, comparative analysis of the MC and HS model results and monitoring results were compared with the pit, can provide a reference for similar projects.

Keywords: numerical analysis; MC model; HS model

1工程概况

上海某地铁车站基坑工程为二地铁线十字相交处，后建南北向车站被已建的东西向车站分隔为南北两个区域，地质条件复杂，道路管线多，交通流量大，周围建筑物密集。本文对其拟建的北侧标准段区域进行分析，基坑南北长约65m、东西宽约25.2m～31.2m，基坑开挖深度约为24.0m，基坑保护等级定为一级。

1.1地质资料

基坑范围内主要涉及杂填土、素填土、褐黄～灰黄色粉质粘土、灰色砂质粉土、灰色淤泥质粘土、灰色粘土、灰色粉质粘土、暗绿～草黄色粉质粘土、草黄色粘质粉土、草黄色粉砂、灰色粉砂、灰色粘土层土。

1.2水文资料

本工程地下水主要有浅部土层中的潜水，及深部粉性土、砂土层中的承压水。上海年平均水位埋深在0.5～0.7m，低水位埋深1.50m。现场测得的地下水位埋深一般在1.15～1.25m之间。

1.3支护结构体系

1.3.1围护结构

围护结构采用1000mm厚地下连续墙，混凝土强度等级为水下C30。标准段地下连续墙深42米，入土比为0.74。

1.3.2支撑

基坑采用钢支撑和混凝土支撑，标准段设9道支撑，第2、4和7道分别为800×1000、1000×1000和1200×1000混凝土支撑，其余均为Φ609×16钢支撑，第3、5道支撑有移撑。

1.3.3地基加固

地基加固采用高压旋喷桩局部抽条加固，标准段加固范围为第六道钢支撑中心以下3米及坑底下3米，加固强度为qu≥1.2MPa。

2数值模拟分析

2.1计算模型

结合监测点位置计算断面取[M-7]轴×[14]轴处，此处基坑宽度为31.2m左右，最大开挖深度为24.00m。根据工程经验，基坑开挖宽度影响范围为开挖深度的3～4倍，基坑开挖深度的影响范围为开挖深度的2～4倍。计算模型可按此取水平方向上长度为75.6m，竖直方向上深度取为60m。

2.2单元模拟

数值模拟中土体采用平面应变15节点2-D等参土体(soil)单元；围护结构地下连续墙采用板桩墙(Plates)单元来模拟；支撑结构按照抗压构件(Anchors)来考虑；土与结构按特殊接触面(Interfaces)单元处理。

2.3边界条件

基坑开挖过程中，数值模拟中的位移边界条件：数值计算中地表为自由边界条件；模型左右两侧边界的侧向位移限制为零；模型底部边界的竖向位移和水平位移均限制为零。

2.4计算工况

①施工此范围内的地下墙结构及地基加固、相应灌注桩、格构柱。

②沿基坑深度从上至下挖土至相应支撑顶部，然后及时掏槽开挖浇筑混凝土支撑或架设钢支撑，直到挖至下一层板梁下，架设第三道钢支撑，浇筑下一层部分板带并架设钢筋砼支撑，中部预留施工阶段出土孔。

③待下一层板带及其钢筋砼支撑达到设计强度后，开挖至第四道钢支撑处。

④及时架设第四道钢支撑（将第三道钢支撑移至第四道钢支撑处），然后开挖至下二层板梁下，架设第五道钢支撑（备一道钢支撑），浇筑下二层部分板带并架设钢筋砼支撑，中部预留施工阶段出土孔。

⑤待下二层板带及其钢筋砼支撑达到设计强度后，开挖至第六道支撑处。

⑥及时架设第六道钢支撑（将第五道钢支撑的备撑移至第六道钢支撑处），然后开挖至坑底。

⑦做素砼垫层并浇筑砼底板，底板达到强度后不得拆除支撑。

2.5土体计算参数

由于本构模型的不同，土体主要计算参数的设置也不同。岩土工程分析计算中，MC本构模型经常使用，但其不能模拟土的一些特殊形为特性，而HS模型可以体现初次加载和卸载-再加载之间的刚度差别。本文分别采用以上两种土体本构模型进行模拟计算。

2.5.1 MC土体模型

在PLAXIS程序当中MC模型，共需要五个参数杨氏模量E、泊桑比υ、内聚力С、内摩擦角φ、剪胀角ψ。根据PLAXIS材料手册，粘土剪胀角ψ≈0，泊桑比С一般取0.15～0.25；内聚力С、内摩擦角ψ为土体强度参数，根据地质勘查报告取用；杨氏模量E为土体的刚度参数，可分为初始切线模量Ei、一点切线模量Etan、割线模量Esec和回弹模量Eur。其中，割线模量Esec代表土体平均刚度，多作为杨氏模量E。0

2.5.2 HS土体模型

Hardening-Soil模型采用Mohr-Coulomb破坏准则，主要涉及到刚度参数压缩模量、参考割线刚度、参考卸荷再加荷模量、卸载再加载泊松比Vur等。

压缩模量：根据地质勘察报告或相关技术标准结合经验选取；

参考割线刚度：按照正常固结粘土，取 (1～2)；

参考卸荷再加荷模量：根据相关工程经验并结合Plaxis中默认缺省关系，取。

卸载再加载泊松比Vur：对砂性土取Vur=0.15，粘性土取Vur=0.20。 [2][3]

2.5.3水泥土加固土体参数

考虑到本工程坑底加固为局部抽条加固，计算模型可从安全角度出发不考虑土体加固。如为满堂加固需考虑水泥土加固影响时加固后水泥土粘聚力c=0.2qu，内摩擦角φ取30°；初始割线模量Ei=E50=(60～154)qu。[2]0

2.6支护结构计算参数

2.6.1地下连续墙

因此HS模型的数据模拟结果更为真实。

3数值模拟计算结果评价

由2知，HS模型比MC模型可以给出更为符合实际情况的模拟结果，因此，对数值模拟计算结果的评价可仅比较HS模型计算得到结果与实际监测数据、规范经验形式的结果。

3.1.1规范控制值与经验预估0

基坑保护等级为一级时，（1）围护结构最大侧移控制指标为0.18%H，为43.2；（2）坑外地表最大沉降控制指标为0.15%H，为36.0。

3.1.2监测布置与监测数据

基坑部分监测项目与监测点开挖至基坑设计标高时，（1）坑内隆起（K1）监测结果为19.90；（2）坑外地表沉降监测数据结果最大值为19.90；（3）地下连续墙水平位移曲线最大值为50.61，已超过报警值。

3.1.3数值模拟计算（HS模型）结果的评价

（1）坑外地表沉降

根据HS模型数据模拟计算结果，由PLAXIS程序导出图8（b）坑外地表沉降Uy的数据表，并形成曲线图，叠加沉降监测点数据得到图12。由图可知，数值模拟坑外地表沉降曲线与监测数据点的分布比较吻合，且比较吻合规范经验公式沉降曲线的形态。

坑外地表沉降数值模拟曲线与沉降监测点比较图

（2）地下连续墙水平位移

采用HS模型数据模拟计算地下连续墙的水平位移小于实际监测值，原因：（1）数值模拟计算过程钢筋混凝土支撑刚度按照设计值取用，实际基坑开挖施工时并非在结构梁板达到100%强度后才开始下层土施工；（2）数值模拟钢支撑掏槽开挖较为理想状态，实际掏槽开挖局部土层还是会产生短时的“局部缺撑”状态。而地下连续墙的最大水平位移的位置，数值模拟（21.75m）与实际监测值（24.00m）略有差距，但其符合文献[2]上海地区地铁基坑统计结果，即“20m以上的基坑围护最大变形值一般位于开挖面以上，平均值0.9h深度处”。

（3）坑内隆起

很明显坑内隆起的实际监测值小于数值模拟结果（73.09m），主要原因是简化计算时未考虑结构梁板局部逆作先行浇筑后自重以及局部土体加固附加荷载的压土作用。

4结论

（1）以PLAXIS岩土工程软件采用HS土体模型对上海软土基坑工程进行数值模拟可以得到较为符合实际情况的结果，可以此评估基坑开挖的变形情况以及对周边环境的影响。

（2）采用钢筋混凝土支撑的板式支护体系围护结构的基坑工程，施工过程中需要确保支撑达到设计要求强度方可进行下层土体开挖施工，否则可能引起围护结构水平位移超限值。

（3）逆作法浇筑的结构梁板自重以及坑底加固作为附加荷载可以有效的抑制坑内隆起量，由于缺少具体先行浇筑结构范围，故本文未作相关的比较计算。

参考文献：

[1]上海市勘察设计行业协会，上海现代建筑设计（集团）有限公司，上海建工（集团）有限公司主编.DG/TJ08-61-2010.上海市工程建设规范《基坑工程技术规范》.上海，2010

[2]刘国彬，王卫东主编.基坑工程手册（第二版）.北京：中国建筑工业出版社，2009

[3]北京金土木软件技术有限公司.PLAXIS岩土工程软件使用指南.北京：人民交通出版社，2010

作者简历：

篇8

一、实验3.1

题目：

考虑线性方程组，，，编制一个能自动选取主元，又能手动选取主元的求解线性代数方程组的Gauss消去过程。

（1）取矩阵，，则方程有解。取计算矩阵的条件数。分别用顺序Gauss消元、列主元Gauss消元和完全选主元Gauss消元方法求解，结果如何？

（2）现选择程序中手动选取主元的功能，每步消去过程都选取模最小或按模尽可能小的元素作为主元进行消元，观察并记录计算结果，若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。

（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。

（4）选取其他你感兴趣的问题或者随机生成的矩阵，计算其条件数，重复上述实验，观察记录并分析实验的结果。

算法介绍

首先，分析各种算法消去过程的计算公式，

顺序高斯消去法：

第k步消去中，设增广矩阵中的元素（若等于零则可以判定系数矩阵为奇异矩阵，停止计算），则对k行以下各行计算，分别用乘以增广矩阵的第行并加到第行，则可将增广矩阵中第列中以下的元素消为零；重复此方法，从第1步进行到第n-1步，则可以得到最终的增广矩阵，即；

列主元高斯消去法：

第k步消去中，在增广矩阵中的子方阵中，选取使得，当时，对中第行与第行交换，然后按照和顺序消去法相同的步骤进行。重复此方法，从第1步进行第n-1步，就可以得到最终的增广矩阵，即；

完全主元高斯消去法：

第k步消去中，在增广矩阵中对应的子方阵中，选取使得，若或，则对中第行与第行、第列与第列交换，然后按照和顺序消去法相同的步骤进行即可。重复此方法，从第1步进行到第n-1步，就可以得到最终的增广矩阵，即；

接下来，分析回代过程求解的公式，容易看出，对上述任一种消元法，均有以下计算公式：

实验程序的设计

一、输入实验要求及初始条件；

二、计算系数矩阵A的条件数及方程组的理论解；

三、对各不同方法编程计算，并输出最终计算结果。

计算结果及分析

（1）

先计算系数矩阵的条件数，结果如下，

可知系数矩阵的条件数较大，故此问题属于病态问题，

b或A的扰动都可能引起解的较大误差；

采用顺序高斯消去法，计算结果为：

最终解为x=(1.000000000000000,

1.000000000000000,

1.000000000000001,

0.999999999999998,

1.000000000000004,

0.999999999999993,

1.000000000000012,

0.999999999999979,

1.000000000000028)T

使用无穷范数衡量误差，得到=2.842170943040401e-14，可以发现，采用顺序高斯消元法求得的解与精确解之间误差较小。通过进一步观察，可以发现，按照顺序高斯消去法计算时，其选取的主元值和矩阵中其他元素大小相近，因此顺序高斯消去法方式并没有对结果造成特别大的影响。

若采用列主元高斯消元法，则结果为：

最终解为x=(1.000000000000000,

1.000000000000000,

1.000000000000000)T

同样使用无穷范数衡量误差，有=0；

若使用完全主元高斯消元法，则结果为

最终解x=(1.000000000000000,

1.000000000000000,

1.000000000000000)T

同样使用无穷范数衡量误差，有=0；

（2）

若每步都选取模最小或尽可能小的元素为主元，则计算结果为

最终解x=(1.000000000000000

1.000000000000000

1.000000000000001

0.999999999999998

1.000000000000004

0.999999999999993

1.000000000000012

0.999999999999979

1.000000000000028)T

使用无穷范数衡量误差，有为2.842170943040401e-14；而完全主元消去法的误差为=0。

从（1）和（2）的实验结果可以发现，列主元消去法和完全主元消去法都得到了精确解，而顺序高斯消去法和以模尽量小的元素为主元的消去法没有得到精确解。在后两种消去法中，由于程序计算时的舍入误差，对最终结果产生了一定的影响，但由于方程组的维度较低，并且元素之间相差不大，所以误差仍比较小。

为进一步分析，计算上述4种方法每步选取的主元数值，并列表进行比较，结果如下：

第n次消元

顺序

列主元

完全主元

模最小

6.000000000000000

4.666666666666667

4.285714285714286

4.133333333333333

4.064516129032258

4.031746031746032

4.015748031496063

4.007843137254902

4.003913894324853

4.001955034213099

0.015617370605469

4.001955034213099

从上表可以发现，对这个方程组而言，顺序高斯消去选取的主元恰好事模尽量小的元素，而由于列主元和完全主元选取的元素为8，与4在数量级上差别小，所以计算过程中的累积误差也较小，最终4种方法的输出结果均较为精确。

在这里，具体解释一下顺序法与模最小法的计算结果完全一致的原因。该矩阵在消元过程中，每次选取主元的一列只有两个非零元素，对角线上的元素为4左右，而其正下方的元素为8，该列其余位置的元素均为0。在这样的情况下，默认的主元也就是该列最小的主元，因此两种方法所得到的计算结果是一致的。

理论上说，完全高斯消去法的误差最小，其次是列主元高斯消去法，而选取模最小的元素作为主元时的误差最大，但是由于方程组的特殊性（元素相差不大并且维度不高），这个理论现象在这里并没有充分体现出来。

（3）

时，重复上述实验过程，各种方法的计算结果如下所示，在这里，仍采用无穷范数衡量绝对误差。

顺序高斯消去法

列主元高斯消去

完全主元高斯消去

选取模最小或尽可能小元素作为主元消去

1.000000000000000

1.000000000000001

0.999999999999998

1.000000000000004

0.999999999999993

1.000000000000014

0.999999999999972

1.000000000000057

0.999999999999886

1.000000000000227

0.999999999999547

1.000000000000902

0.999999999998209

1.000000000003524

0.999999999993179

1.000000000012732

0.999999999978173

1.000000000029102

1.000000000000000

1.000000000000001

0.999999999999998

1.000000000000004

0.999999999999993

1.000000000000014

0.999999999999972

1.000000000000057

0.999999999999886

1.000000000000227

0.999999999999547

1.000000000000902

0.999999999998209

1.000000000003524

0.999999999993179

1.000000000012732

0.999999999978173

1.000000000029102

2.910205409989430e-11

可以看出，此时列主元和完全主元的计算结果仍为精确值，而顺序高斯消去和模尽可能小方法仍然产生了一定的误差，并且两者的误差一致。与n=10时候的误差比相比，n=20时的误差增长了大约1000倍，这是由于计算过程中舍入误差的不断累积所致。所以，如果进一步增加矩阵的维数，应该可以看出更明显的现象。

（4）

不同矩阵维度下的误差如下，在这里，为方便起见，选取2-条件数对不同维度的系数矩阵进行比较。

维度

条件数

顺序消去

列主元

完全主元

模尽量小

1.7e+3

2.84e-14

1.8e+6

2.91e-11

5.7e+7

9.31e-10

1.8e+9

2.98e-08

1.9e+12

3.05e-05

3.8e+16

3.28e+04

3.88e-12

3.28e+04

8.5e+16

3.52e+13

4.2e-3

3.52e+13

从上表可以看出，随着维度的增加，不同方法对计算误差的影响逐渐体现，并且增长较快，这是由于舍入误差逐步累计而造成的。不过，方法二与方法三在维度小于40的情况下都得到了精确解，这两种方法的累计误差远比方法一和方法四慢；同样地，出于与前面相同的原因，方法一与方法四的计算结果保持一致，方法二与方法三的计算结果保持一致。

结论

本文矩阵中的元素差别不大，模最大和模最小的元素并没有数量级上的差异，因此，不同的主元选取方式对计算结果的影响在维度较低的情况下并不明显，四种方法都足够精确。

对比四种方法，可以发现采用列主元高斯消去或者完全主元高斯消去法，可以尽量抑制误差，算法最为精确。不过，对于低阶的矩阵来说，四种方法求解出来的结果误差均较小。

另外，由于完全选主元方法在选主元的过程中计算量较大，而且可以发现列主元法已经可以达到很高的精确程度，因而在实际计算中可以选用列主元法进行计算。

附录：程序代码

clear

clc;

format

long;

%方法选择

n=input(＇矩阵A阶数:n=＇);

disp(＇选取求解方式＇);

disp(＇1

顺序Gauss消元法，2

列主元Gauss消元法，3

完全选主元Gauss消元法，4

模最小或近可能小的元素作为主元＇);

a=input(＇求解方式序号：＇);

%赋值A和b

A=zeros(n,n);

b=zeros(n,1);

for

i=1:n

A(i,i)=6;

i>1

A(i,i-1)=8;

end

A(i,i+1)=1;

end

for

i=1:n

for

j=1:n

b(i)=b(i)+A(i,j);

end

disp(＇给定系数矩阵为：＇);

disp(＇右端向量为：＇);

%求条件数及理论解

disp(＇线性方程组的精确解：＇);

X=(A＼b)＇

fprintf(＇矩阵A的1-条件数:

＼n＇,cond(A,1));

fprintf(＇矩阵A的2-条件数:

＼n＇,cond(A));

fprintf(＇矩阵A的无穷-条件数:

＼n＇,cond(A,inf));

%顺序Gauss消元法

a==1

A1=A;b1=b;

for

k=1:n

A1(k,k)==0

disp(＇主元为零，顺序Gauss消元法无法进行＇);

break

end

fprintf(＇第%d次消元所选取的主元：%g＼n＇,k,A1(k,k))

%disp(＇此次消元后系数矩阵为：＇);

%A1

for

p=k+1:n

l=A1(p,k)/A1(k,k);

A1(p,k:n)=A1(p,k:n)-l*A1(k,k:n);

b1(p)=b1(p)-l*b1(k);

end

x1(n)=b1(n)/A1(n,n);

for

k=n-1:-1:1

for

w=k+1:n

b1(k)=b1(k)-A1(k,w)*x1(w);

end

x1(k)=b1(k)/A1(k,k);

end

disp(＇顺序Gauss消元法解为:＇);

disp(x1);

disp(＇所求解与精确解之差的无穷-范数为＇);

norm(x1-X,inf)

end

%列主元Gauss消元法

a==2

A2=A;b2=b;

for

k=1:n

[max_i,max_j]=find(A2(:,k)==max(abs(A2(k:n,k))));

max_i~=k

A2_change=A2(k,:);

A2(k,:)=A2(max_i,:);

A2(max_i,:)=A2_change;

b2_change=b2(k);

b2(k)=b2(max_i);

b2(max_i)=b2_change;

end

A2(k,k)==0

disp(＇主元为零，列主元Gauss消元法无法进行＇);

break

end

fprintf(＇第%d次消元所选取的主元：%g＼n＇,k,A2(k,k))

%disp(＇此次消元后系数矩阵为：＇);

%A2

for

p=k+1:n

l=A2(p,k)/A2(k,k);

A2(p,k:n)=A2(p,k:n)-l*A2(k,k:n);

b2(p)=b2(p)-l*b2(k);

end

x2(n)=b2(n)/A2(n,n);

for

k=n-1:-1:1

for

w=k+1:n

b2(k)=b2(k)-A2(k,w)*x2(w);

end

x2(k)=b2(k)/A2(k,k);

end

disp(＇列主元Gauss消元法解为:＇);

disp(x2);

disp(＇所求解与精确解之差的无穷-范数为＇);

norm(x2-X,inf)

end

%完全选主元Gauss消元法

a==3

A3=A;b3=b;

for

k=1:n

VV=eye(n);

[max_i,max_j]=find(A3(k:n,k:n)==max(max(abs(A3(k:n,k:n)))));

numel(max_i)==0

[max_i,max_j]=find(A3(k:n,k:n)==-max(max(abs(A3(k:n,k:n)))));

end

W=eye(n);

W(max_i(1)+k-1,max_i(1)+k-1)=0;

W(k,k)=0;

W(max_i(1)+k-1,k)=1;

W(k,max_i(1)+k-1)=1;

V=eye(n);

V(k,k)=0;

V(max_j(1)+k-1,max_j(1)+k-1)=0;

V(k,max_j(1)+k-1)=1;

V(max_j(1)+k-1,k)=1;

A3=W*A3*V;

b3=W*b3;

VV=VV*V;

A3(k,k)==0

disp(＇主元为零，完全选主元Gauss消元法无法进行＇);

break

end

fprintf(＇第%d次消元所选取的主元：%g＼n＇,k,A3(k,k))

%disp(＇此次消元后系数矩阵为：＇);

%A3

for

p=k+1:n

l=A3(p,k)/A3(k,k);

A3(p,k:n)=A3(p,k:n)-l*A3(k,k:n);

b3(p)=b3(p)-l*b3(k);

end

x3(n)=b3(n)/A3(n,n);

for

k=n-1:-1:1

for

w=k+1:n

b3(k)=b3(k)-A3(k,w)*x3(w);

end

x3(k)=b3(k)/A3(k,k);

end

disp(＇完全选主元Gauss消元法解为:＇);

disp(x3);

disp(＇所求解与精确解之差的无穷-范数为＇);

norm(x3-X,inf)

end

%模最小或近可能小的元素作为主元

a==4

A4=A;b4=b;

for

k=1:n

AA=A4;

AA(AA==0)=NaN;

[min_i,j]=find(AA(k:n,k)==min(abs(AA(k:n,k))));

numel(min_i)==0

[min_i,j]=find(AA(k:n,k)==-min(abs(AA(k:n,k:n))));

end

W=eye(n);

W(min_i(1)+k-1,min_i(1)+k-1)=0;

W(k,k)=0;

W(min_i(1)+k-1,k)=1;

W(k,min_i(1)+k-1)=1;

A4=W*A4;

b4=W*b4;

A4(k,k)==0

disp(＇主元为零，模最小Gauss消元法无法进行＇);

break

end

fprintf(＇第%d次消元所选取的主元：%g＼n＇,k,A4(k,k))

%A4

for

p=k+1:n

l=A4(p,k)/A4(k,k);

A4(p,k:n)=A4(p,k:n)-l*A4(k,k:n);

b4(p)=b4(p)-l*b4(k);

end

x4(n)=b4(n)/A4(n,n);

for

k=n-1:-1:1

for

w=k+1:n

b4(k)=b4(k)-A4(k,w)*x4(w);

end

x4(k)=b4(k)/A4(k,k);

end

disp(＇模最小Gauss消元法解为:＇);

disp(x4);

disp(＇所求解与精确解之差的无穷-范数为＇);

norm(x4-X,inf)

end

二、实验3.3

题目：

考虑方程组的解，其中系数矩阵H为Hilbert矩阵：

这是一个著名的病态问题。通过首先给定解（例如取为各个分量均为1）再计算出右端的办法给出确定的问题。

（1）选择问题的维数为6，分别用Gauss消去法（即LU分解）、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何。

（2）逐步增大问题的维数，仍用上述的方法来解它们，计算的结果如何？计算的结果说明的什么？

（3）讨论病态问题求解的算法。

算法设计

对任意线性方程组，分析各种方法的计算公式如下，

（1）Gauss消去法：

首先对系数矩阵进行LU分解，有，则原方程转化为，令，则原方程可以分为两步回代求解：

具体方法这里不再赘述。

（2）J迭代法：

首先分解，再构造迭代矩阵，其中

，进行迭代计算，直到误差满足要求。

（3）GS迭代法：

首先分解，再构造迭代矩阵

，其中

，进行迭代计算，直到误差满足要求。

（4）SOR迭代法：

首先分解，再构造迭代矩阵

，其中，进行迭代计算，直到误差满足要求。

实验过程

一、根据维度n确定矩阵H的各个元素和b的各个分量值；

二、选择计算方法（

Gauss消去法，J迭代法，GS迭代法，SOR迭代法），对迭代法设定初值，此外SOR方法还需要设定松弛因子；

三、进行计算，直至满足误差要求（对迭代法，设定相邻两次迭代结果之差的无穷范数小于0.0001；

对SOR方法，设定为输出迭代100次之后的结果及误差值），输出实验结果。

计算结果及分析

（1）时，问题可以具体定义为

计算结果如下，

Gauss消去法

第1次消元所选取的主元是：1

第2次消元所选取的主元是：0.0833333

第3次消元所选取的主元是：0.00555556

第4次消元所选取的主元是：0.000357143

第5次消元所选取的主元是：2.26757e-05

第6次消元所选取的主元是：1.43155e-06

解得X=(0.999999999999228

1.000000000021937

0.999999999851792

1.000000000385369

0.999999999574584

1.000000000167680)T

使用无穷范数衡量误差，可得=4.254160357319847e-10；

J迭代法

设定迭代初值为零，计算得到

J法的迭代矩阵B的谱半径为4.30853＞1，所以J法不收敛；

GS迭代法

设定迭代初值为零，计算得到GS法的迭代矩阵G的谱半径为：0.999998＜1，故GS法收敛，经过541次迭代计算后，结果为X=(1.001178105812706

0.999144082651860

0.968929093984902

1.047045569989162

1.027323158370281

0.954352032784608)T

使用无穷范数衡量误差，有=0.047045569989162；

SOR迭代法

设定迭代初值为零向量，并设定，计算得到SOR法迭代矩阵谱半径为0.999999433815223，经过100次迭代后的计算结果为

X=(1.003380614145078

0.962420297458423

1.031857023134559

1.061814901289881

1.014037815827164

0.917673642493527)T；

使用无穷范数衡量误差，有=0.082326357506473；

对SOR方法，可变，改变值，计算结果可以列表如下

迭代次数

100

迭代矩阵的谱半径

0.999999433815223

0.999998867083155

0.999996830135013

0.999982309342386

1.003653917714694

0.974666041209353

1.011814573842440

1.042837929171827

1.017190220902681

0.945462001336268

1.014676015634604

0.896636864424096

1.090444578936265

1.107070542628148

1.006315452225331

0.873244842279255

1.028022215505147

0.790604920509843

1.267167365524072

1.061689730857891

0.990084054872602

0.846005956774467

1.051857392323966

0.653408758549156

1.486449891152510

0.783650360698119

1.349665420488270

0.664202350634588

0.054537998663732

0.126755157720745

0.267167365524072

0.486449891152510

可以发现，松弛因子的取值对迭代速度造成了不同的影响，上述四种方法中，松弛因子=0.5时，收敛相对较快。

综上，四种算法的结果列表如下：

算法

Gauss消去法

Jacobi法

GS法

SOR法（取）

迭代次数

不收敛

541

100

迭代矩阵的谱半径

4.30853

0.999998

0.999999433815223

0.999999999999228

1.000000000021937

0.999999999851792

1.000000000385369

0.999999999574584

1.000000000167680

1.001178105812706

0.999144082651860

0.968929093984902

1.047045569989162

1.027323158370281

0.954352032784608

1.003380614145078

0.962420297458423

1.031857023134559

1.061814901289881

1.014037815827164

0.917673642493527

4.254160357319847e-10

0.047045569989162

0.082326357506473

计算可得，矩阵H的条件数为>>1，所以这是一个病态问题。由上表可以看出，四种方法的求解都存在一定的误差。下面分析误差的来源：

LU分解方法的误差存在主要是由于Hilbert矩阵各元素由分数形式转换为小数形式时，不能除尽情况下会出现舍入误差，在进行LU分解时也存在这个问题，所以最后得到的结果不是方程的精确解

，但结果显示该方法的误差非常小；

Jacobi迭代矩阵的谱半径为4.30853，故此迭代法不收敛；

GS迭代法在迭代次数为541次时得到了方程的近似解，其误差约为0.05

，比较大。GS迭代矩阵的谱半径为0.999998，很接近1，所以GS迭代法收敛速度较慢；

SOR迭代法在迭代次数为100次时误差约为0.08，误差较大。SOR迭代矩阵的谱半径为0.999999，也很接近1，所以时SOR迭代法收敛速度不是很快，但是相比于GS法，在迭代速度方面已经有了明显的提高；另外，对不同的，SOR方法的迭代速度会相应有变化，如果选用最佳松弛因子，可以实现更快的收敛；

（2）

考虑不同维度的情况，时，

算法

Gauss消去

J法

GS法

SOR法（w=0.5）

计算结果

0.999999999966269

1.000000001809060

0.999999976372676

1.000000127868103

0.999999655764116

1.000000487042164

0.999999653427125

1.000000097774747

0.997829221945349

1.037526203106839

0.896973261976015

1.020345136375036

1.069071166932576

1.051179995036612

0.996814757185364

0.926343237325536

1.012938972275634

0.939713836855171

0.988261805073081

1.064637090535154

1.083633345093974

1.045060177115514

0.970603024778469

0.880212649657655

迭代次数

356

100

谱半径

6.04213

0.999999999208776

时，

算法

Gauss消去法

Jacobi法

GS法

SOR法（w=0.5）

计算结果

0.999999994751197

1.000000546746354

0.999985868343700

1.000157549468631

0.999063537004329

1.003286333127805

0.992855789229370

1.009726486881556

0.991930155925812

1.003729850349020

0.999263885025643

0.997442073306751

1.019069909358409

0.992278247786739

0.956441858313237

0.986420333361353

1.021301611956591

1.038701026806608

1.035942773498533

1.016693763149422

0.985716454946250

0.947181287500697

1.015776039786572

0.966429147064483

0.928674868157910

0.996931548482727

1.066737803913537

1.097792430596468

1.088030440855069

1.048110620811192

0.989919418572424

0.922840813704142

0.853252417221922

迭代次数

1019

100

谱半径

8.64964

0.999999999999966

时

算法

Gauss消去法

Jacobi法

GS法

SOR法（w=0.5）

计算结果

0.999999968723799

1.000002417094896

0.999994922439769

0.998640261957706

1.025668111139297

0.781933485305194

2.066840925345890

-2.279036697492128

7.532393125791018

-7.355047567109081

7.380667063930484

-1.129041418095142

0.425748747257065

1.733284233971601

0.817952344733362

不收敛

1.004385740641590

1.046346067877554

0.907178347707729

0.905763455949053

0.972521802788457

1.043731445367903

1.091535169448764

1.110090020703944

1.103129684679768

1.077168651146056

1.038514736265176

0.992259990832041

0.942151390478003

0.890785366684065

0.839876442493220

迭代次数

262

100

谱半径

6.04213

＞1

1.000000000000000

8.355047567109082

0.160123557506780

分析以上结果可以发现，随着n值的增加，Gauss消去法误差逐渐增大，而且误差增大的速度很快，在维数小于等于10情况下，Gauss消去法得到的结果误差较小；但当维数达到15时，计算结果误差已经达到精确解的很多倍；

J法迭代不收敛，无论n如何取值，其谱半径始终大于1，因而J法不收敛，所以J迭代法不能用于Hilbert矩阵的求解；

对于GS迭代法和SOR迭代法，两种方法均收敛，GS迭代法是SOR迭代法松弛因子取值为1的特例，SOR方法受到取值的影响，会有不同的收敛情况。可以得出GS迭代矩阵的谱半径小于1但是很接近1，收敛速度很慢。虽然随着维数的增大，所需迭代的次数逐渐减少，但是当维数达到15的时候，GS法已经不再收敛。因此可以得出结论，GS迭代方法在Hilbert矩阵维数较低时，能够在一定程度上满足迭代求解的需求，不过迭代的速度很慢。另外，随着矩阵维数的增加，

SOR法的误差水平基本稳定，而且误差在可以接受的范围之内。

经过比较可以得出结论，如果求解较低维度的Hibert矩阵问题，Gauss消去法、GS迭代法和SOR迭代法均可使用，且Gauss消去法的结果精确度较高；如果需要求解较高维度的Hibert矩阵问题，只有采用SOR迭代法。

（3）

系数矩阵的条件数较大时，为病态方程。由实验可知，Gauss法在解上述方程时，结果存在很大的误差。而对于收敛的迭代法，可以通过选取最优松弛因子的方法来求解，虽然迭代次数相对较多，但是结果较为精确。

总体来看，对于一般病态方程组的求解，可以采用以下方式：

低维度下采用Gauss消去法直接求解是可行的；

Jacobi迭代方法不适宜于求解病态问题；

GS迭代方法可以解决维数较低的病态问题，但其谱半径非常趋近于1，导致迭代算法收敛速度很慢，维数较大的时候，GS法也不再收敛；

SOR方法较适合于求解病态问题，特别是矩阵维数较高的时候，其优势更为明显。

采用高精度的运算，如选用双倍或更多倍字长的运算，可以提高收敛速度；

可以对原方程组作某些预处理，从而有效降低系数矩阵的条件数。

实验结论

（1）对Hibert矩阵问题，其条件数会随着维度的增加迅速增加，病态性会越来越明显；在维度较低的时候，Gauss消去法、GS迭代法和SOR迭代法均可使用，且可以优先使用Gauss消去法；如果需要求解较高维度的Hibert矩阵问题，只有SOR迭代法能够求解。

（2）SOR方法比较适合于求解病态问题，特别是矩阵维数较高的时候，其优点更为明显。从本次实验可以看出，随着矩阵维数的增大，SOR方法所需的迭代次数减少，而且误差基本稳定，是解决病态问题的适宜方法。

附录：程序代码

clear

all

clc;

format

long;

%矩阵赋值

n=input(＇矩阵H的阶数:n=＇);

for

i=1:n

for

j=1:n

H(i,j)=1/(i+j-1);

end

b=H*ones(n,1);

disp(＇H矩阵为：＇);

disp(＇向量b：＇);

%方法选择

disp(＇选取求解方式＇);

disp(＇1

Gauss消去法，2

J迭代法，3

GS迭代法，4

SOR迭代法＇);

a=input(＇求解方式序号：＇);

%Gauss消去法

a==1;

H1=H;b1=b;

for

k=1:n

H1(k,k)==0

disp(＇主元为零，Gauss消去法无法进行＇);

break

end

fprintf(＇第%d次消元所选取的主元是：%g＼n＇,k,H1(k,k))

for

p=k+1:n

m5=-H1(p,k)/H1(k,k);

H1(p,k:n)=H1(p,k:n)+m5*H1(k,k:n);

b1(p)=b1(p)+m5*b1(k);

end

x1(n)=b1(n)/H1(n,n);

for

k=n-1:-1:1

for

v=k+1:n

b1(k)=b1(k)-H1(k,v)*x1(v);

end

x1(k)=b1(k)/H1(k,k);

end

disp(＇Gauss消去法解为:＇);

disp(x1);

disp(＇解与精确解之差的无穷范数＇);

norm((x1-a),inf)

end

D=diag(diag(H));

L=-tril(H,-1);

U=-triu(H,1);

%J迭代法

a==2;

%给定初始x0

ini=input(＇初始值设定:x0=＇);

x0(:,1)=ini*diag(ones(n));

disp(＇初始解向量为:＇);

xj(:,1)=x0(:,1);

B=(D^(-1))*(L+U);

f=(D^(-1))*b;

fprintf(＇(J法B矩阵谱半径为：%g＼n＇,vrho(B));

vrho(B)

for

m2=1:5000

xj(:,m2+1)=B*xj(:,m2)+fj;

norm((xj(:,m2+1)-xj(:,m2)),inf)

break

end

disp(＇J法计算结果为:＇);

xj(:,m2+1)

disp(＇解与精确解之差的无穷范数＇);

norm((xj(:,m2+1)-diag(ones(n))),inf)

disp(＇J迭代法迭代次数:＇);

else

disp(＇由于B矩阵谱半径大于1，因而J法不收敛＇);

end

%GS迭代法

a==3;

%给定初始x0

ini=input(＇初始值设定:x0=＇);

x0(:,1)=ini*diag(ones(n));

disp(＇初始解向量为:＇);

xG(:,1)=x0(:,1);

G=inv(D-L)*U;

fG=inv(D-L)*b;

fprintf(＇GS法G矩阵谱半径为：%g＼n＇,vrho(G));

vrho(G)

for

m3=1:5000

xG(:,m3+1)=G*xG(:,m3)+fG;

norm((xG(:,m3+1)-xG(:,m3)),inf)

break;

end

disp(＇GS迭代法计算结果:＇);

xG(:,m3+1)

disp(＇解与精确解之差的无穷范数＇);

norm((xG(:,m3+1)-diag(ones(n))),inf)

disp(＇GS迭代法迭代次数:＇);

else

disp(＇由于G矩阵谱半径大于1，因而GS法不收敛＇);

end

%SOR迭代法

a==4;

%给定初始x0

ini=input(＇初始值设定:x0=＇);

x0(:,1)=ini*diag(ones(n));

disp(＇初始解向量为:＇);

A=H;

for

i=1:n

b(i)=sum(A(i,:));

end

x_star=ones(n,1);

format

long

w=input(＇松弛因子:w=＇);

Lw=inv(D-w*L)*((1-w)*D+w*U);

f=w*inv(D-w*L)*b;

disp(＇迭代矩阵的谱半径：＇)

p=vrho(Lw)

time_max=100;%迭代次数

x=zeros(n,1);%迭代初值

for

i=1:time_max

x=Lw*x+f;

end

disp(＇SOR迭代法得到的解为＇);

disp(＇解与精确解之差的无穷范数＇);

norm((x_star-x),inf)

end

pause

三、实验4.1

题目：

对牛顿法和拟牛顿法。进行非线性方程组的数值求解

（1）用上述两种方法，分别计算下面的两个例子。在达到精度相同的前提下，比较其迭代次数、CPU时间等。

（2）取其他初值，结果又如何？反复选取不同的初值，比较其结果。

（3）总结归纳你的实验结果，试说明各种方法适用的问题。

算法设计

对需要求解的非线性方程组而言，牛顿法和拟牛顿法的迭代公式如下，

（1）牛顿法：

牛顿法为单步迭代法，需要取一个初值。

（2）拟牛顿法：（Broyden秩1法）

其中，

拟牛顿法不需要求解的导数，因此节省了大量的运算时间，但需要给定矩阵的初值，取为。

实验过程

一、输入初值；

二、根据误差要求，按公式进行迭代计算；

三、输出数据；

计算结果及分析

（1）首先求解方程组（1），在这里，设定精度要求为，

方法

牛顿法

拟牛顿法

初始值

计算结果X

0.905539609855914

0.905539493347151

1.085219168370031

1.085218882394940

0.672193668718306

0.672193293825304

迭代次数

CPU计算时间/s

3.777815

2.739349

可以看出，在初始值相同情况下，牛顿法和拟牛顿法在达到同样计算精度情况下得到的结果基本相同，但牛顿法的迭代次数明显要少一些，但是，由于每次迭代都需要求解矩阵的逆，所以牛顿法每次迭代的CPU计算时间更长。

之后求解方程组（2），同样设定精度要求为

方法

牛顿法

拟牛顿法

初始值

计算结果X

0.500000000009699

0.499999994673600

0.000000001063428

0.000000572701856

-0.523598775570483

-0.523598762908871

迭代次数

CPU计算时间/s

2.722437

3.920195

同样地，可以看出，在初始值相同情况下，牛顿法和拟牛顿法在达到同样计算精度情况下得到的结果是基本相同的，但牛顿法的迭代次数明显要少，但同样的，由于每次迭代中有求解矩阵的逆的运算，牛顿法每次迭代的CPU计算时间较长。

（2）对方程组（1），取其他初值，计算结果列表如下，同样设定精度要求为

初始值

方法

牛顿法

拟牛顿法

计算结果

0.905539609855914

1.085219168370031

0.672193668718305

9.211852562357894

-5.574005400255346

18.118173639381205

迭代次数

CPU计算时间/s

3.907164

4.818019

计算结果

0.905539609855914

1.085219168370031

0.672193668718305

9.211849682114591

-5.573999165383549

18.118182491302807

迭代次数

2735

CPU计算时间/s

8.127286

5.626023

计算结果

0.905539609855914

1.085219168370031

0.672193668718306

0.905539493347151

1.085218882394940

0.672193293825304

迭代次数

CPU计算时间/s

3.777815

2.739349

计算结果

0.905539609855914

1.085219168370031

0.672193668718306

0.905548384395773

1.085220084502458

0.672219278250136

迭代次数

188

CPU计算时间/s

3.835697

2.879070

计算结果

9.211852448563722

-5.574005155684773

18.118173976918605

Matlab警告矩阵接近奇异值，程序进入长期循环计算中

迭代次数

CPU计算时间/s

4.033868

计算结果

0.905539609857335

1.085219168371536

0.672193668734922

Matlab警告矩阵接近奇异值，程序进入长期循环计算中

迭代次数

CPU计算时间/s

12.243263

从上表可以发现，方程组（1）存在另一个在(9.2,

-5.6,

18.1)T附近的不动点，初值的选取会直接影响到牛顿法和拟牛顿法最后的收敛点。

总的来说，设定的初值离不动点越远，需要的迭代次数越多，因而初始值的选取非常重要，合适的初值可以更快地收敛，如果初始值偏离精确解较远，会出现迭代次数增加直至无法收敛的情况；

由于拟牛顿法是一种近似方法，拟牛顿法需要的的迭代次数明显更多，而且收敛情况不如牛顿法好（初值不够接近时，甚至会出现奇异矩阵的情况），但由于牛顿法的求解比较复杂，计算时间较长；

同样的，对方程组（2），取其他初值，计算结果列表如下，同样设定精度要求为

初始值

方法

牛顿法

拟牛顿法

计算结果

0.500000000009699

0.000000001063428

-0.523598775570483

0.499999994673600

0.000000572701856

-0.523598762908871

迭代次数

CPU计算时间/s

2.722437

3.920195

计算结果

0.500000000011085

0.000000001215427

-0.523598775566507

0.331099293590753

-0.260080189442266

76.532092226437129

迭代次数

CPU计算时间/s

5.047111

5.619752

计算结果

0.500000000000916

0.000000000100410

-0.523598775595672

1.0e+02

-0.001221250784775

-0.000149282572886

1.754185881622843

迭代次数

CPU计算时间/s

3.540668

3.387829

计算结果

0.500000000000152

0.000000000016711

-0.523598775597862

1.0e+04

0.000026556790770

-0.000020396841295

1.280853105748650

迭代次数

CPU计算时间/s

2.200571

2.640901

计算结果

0.500000000000005

0.000000000000503

-0.523598775598286

矩阵为奇异值，无法输出准确结果

迭代次数

CPU计算时间/s

1.719072

计算结果

0.500000000002022

0.000000000221686

-0.523598775592500

矩阵为奇异值，无法输出准确结果

迭代次数

149

CPU计算时间/s

2.797116

计算结果

矩阵为奇异值，无法输出准确结果

迭代次数

CPU计算时间/s

在这里，与前文类似的发现不再赘述。

从这里看出，牛顿法可以在更大的区间上实现压缩映射原理，可以在更大的范围上选取初值并最终收敛到精确解附近；

在初始值较接近于不动点时，牛顿法和拟牛顿法计算所得到的结果是基本相同的，虽然迭代次数有所差别，但计算总的所需时间相近。

（3）

牛顿法在迭代过程中用到了矩阵的求逆，其迭代收敛的充分条件是迭代满足区间上的映内性，对于矩阵的求逆过程比较简单，所以在较大区间内满足映内性的问题适合应用牛顿法进行计算。一般而言，对于函数单调或者具有单值特性的函数适合应用牛顿法，其对初始值敏感程度较低，算法具有很好的收敛性。

另外，需要说明的是，每次计算给出的CPU时间与计算机当时的运行状态有关，同时，不同代码的运行时间也不一定一致，所以这个数据并不具有很大的参考价值。

实验结论

对牛顿法和拟牛顿法，都存在初始值越接近精确解，所需的迭代次数越小的现象；

在应用上，牛顿法和拟牛顿法各有优势。就迭代次数来说，牛顿法由于更加精确，所需的迭代次数更少；但就单次迭代来说，牛顿法由于计算步骤更多，且计算更加复杂，因而每次迭代所需的时间更长，而拟牛顿法由于采用了简化的近似公式，其每次迭代更加迅速。当非线性方程组求逆过程比较简单时，如方程组1的情况时，拟牛顿法不具有明显的优势；而当非线性方程组求逆过程比较复杂时，如方程组2的情况，拟牛顿法就可以体现出优势，虽然循环次数有所增加，但是CPU耗时反而更少。

另外，就方程组压缩映射区间来说，一般而言，对于在区间内函数呈现单调或者具有单值特性的函数适合应用牛顿法，其对初始值敏感程度较低，使算法具有很好的收敛性；而拟牛顿法由于不需要在迭代过程中对矩阵求逆，而是利用差商替代了对矩阵的求导，所以即使初始误差较大时，其倒数矩阵与差商偏差也较小，所以对初始值的敏感程度较小。

附录：程序代码

%方程1，牛顿法

tic;

format

long;

%%初值

disp(＇请输入初值＇);

a=input(＇第1个分量为：＇);

b=input(＇第2个分量为：＇);

c=input(＇第3个分量为：＇);

disp(＇所选定初值为＇);

x=[a;b;c]

%%误差要求

E=0.0001;

%%迭代

i=0;

e=2*E;

while

e>E

F=[12*x(1)-x(2)^2-4*x(3)-7;x(1)^2+10*x(2)-x(3)-11;x(2)^3+10*x(3)-8];

f=[12,-2*x(2),-4;2*x(1),10,-1;0,3*x(2)^2,10];

det_x=((f)^(-1))*(-F);

x=x+det_x;

e=max(norm(det_x));

i=i+1;

end

disp(＇迭代次数＇);

toc;

%方程1，拟牛顿法

tic;

format

long;

%%初值

disp(＇请输入初值＇);

a=input(＇第1个分量为：＇);

b=input(＇第2个分量为：＇);

c=input(＇第3个分量为：＇);

disp(＇所选定初值为＇);

x0=[a;b;c]

%%误差要求

E=0.0001;

%%迭代

i=0;

e=2*E;

A0=eye(3);

while

e>E

F0=[12*x0(1)-x0(2)^2-4*x0(3)-7;x0(1)^2+10*x0(2)-x0(3)-11;x0(2)^3+10*x0(3)-8];

x1=x0-A0^(-1)*F0;

s=x1-x0;

F1=[12*x1(1)-x1(2)^2-4*x1(3)-7;x1(1)^2+10*x1(2)-x1(3)-11;x1(2)^3+10*x1(3)-8];

y=F1-F0;

A1=A0+(y-A0*s)*s＇/(s＇*s);

x0=x1;

A0=A1;

e=max(norm(s));

i=i+1;

end

disp(＇迭代次数＇);

toc;

%方程2，牛顿法

tic;

format

long;

%%初值

disp(＇请输入初值＇);

a=input(＇第1个分量为：＇);

b=input(＇第2个分量为：＇);

c=input(＇第3个分量为：＇);

disp(＇所选定初值为＇);

x=[a;b;c]

%%误差要求

E=0.0001;

%%迭代

i=0;

e=2*E;

while

e>E

F=[3*x(1)-cos(x(2)*x(3))-0.5;x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;exp(1)^(-x(1)*x(2))+20*x(3)+(10*pi-3)/3];

f=[3,x(3)*sin(x(2)*x(3)),x(2)*sin(x(2)*x(3));2*x(1),-162*x(2)-81/5,cos(x(3));-x(2)*exp(1)^(-x(1)*x(2)),-x(1)*exp(1)^(-x(1)*x(2)),20];

det_x=((f)^(-1))*(-F);

x=x+det_x;

e=max(norm(det_x));

i=i+1;

end

disp(＇迭代次数＇);

toc;

%方程2，拟牛顿法

tic;

format

long;

%%初值

disp(＇请输入初值＇);

a=input(＇第1个分量为：＇);

b=input(＇第2个分量为：＇);

c=input(＇第3个分量为：＇);

disp(＇所选定初值为＇);

x0=[a;b;c]

%%误差要求

E=0.0001;

%%迭代

i=0;

e=2*E;

A0=eye(3);

while

e>E

F0=[3*x0(1)-cos(x0(2)*x0(3))-0.5;x0(1)^2-81*(x0(2)+0.1)^2+sin(x0(3))+1.06;exp(1)^(-x0(1)*x0(2))+20*x0(3)+(10*pi-3)/3];

x1=x0-A0^(-1)*F0;

s=x1-x0;

F1=[3*x1(1)-cos(x1(2)*x1(3))-0.5;x1(1)^2-81*(x1(2)+0.1)^2+sin(x1(3))+1.06;exp(1)^(-x1(1)*x1(2))+20*x1(3)+(10*pi-3)/3];

y=F1-F0;

A1=A0+(y-A0*s)*s＇/(s＇*s);

x0=x1;

A0=A1;

e=max(norm(s));

i=i+1;

end

disp(＇迭代次数＇);

篇9

关键词：教学辅助；非线性方程求根；用户图形界面

中图分类号：TP319文献标识码：A文章编号：1009-3044(2012)13-3140-03

The Application of MATLAB GUI in Numerical Analysis Course Teaching

CHEN Li-hong, ZHOU Zhi-gang

(Department of Mathematics and Computer College, Wuhan Textile University, Wuhan 430073, China)

Abstract：According to the difficult in numerical analysis course teaching, it is discussed that the MATLAB GUI(Graphical User Interfaces，GUI) application in numerical analysis course teaching with solving nonlinear equations as example. the view is proposed that designing MATLAB GUI for numerical analysis will fully arouse teachers and students both aspects of the enthusiasm to improve quality of numerical analysis course combined with MATLAB GUI function.

Key words：auxiliary teaching；solving nonlinear equations；MATLAB GUI

《数值分析》是理工科院校数学、力学、物理、计算机等专业的教材，它在专业课程体系中占有重要地位。该课程的主要任务是研究用计算机解决数学问题的数值方法及其理论，它的内容包括函数的数值逼近、数值微分与数值积分、非线性方程数值解、数值线性代数、微分方程数值解等,这些均与计算机紧密结合在一起。如果采用传统的教学方式，一方面需要花大量的时间在黑板上绘图和计算，在有限的学时内无法进行内容的扩展；另一方面学生理解和接授知识时感觉枯燥、难度大。MATLAB软件的推出为该门课的教学提供了有利工具，利用MATLAB强大的绘图及仿真功能，可以将抽象的内容以形象的图形表示出来，既可演示复杂系统的未知结果，又可改变系统参数，演示系统随参数变化的变化结果或变化趋势，有助于学生对抽象理论的理解。然而在课堂上应用MATLAB演示所讲授内容，需要临时编程，这对于有限的课堂时间多有不便，而且界面也不直观，若能将MATLAB的可开发的GUI功能结合数值计算中的典型算法构造开放式的用户界面，既可充分发挥MATLAB的强大的计算功能，又可避免记忆繁琐的命令，不仅方便老师在课堂直观演示，而且便于学生课下自己设计系统，添加代码实现更多的演示功能，将会充分调动教师和学生两方面的积极性，全面提高课程的教学质量。

1图形用户界面设计简介

图形用户界面是由窗口、光标、按键、菜单、文字说明等对象（Objects）构成的一个用户界面，用户通过一定的方法（如鼠标或键盘）选择激活这些图形对象，使计算机产生某种动作或变化，如实现计算或绘图等。MATLAB的GUI编程可以用两种方式实现。一是GUI设计工具GUIDE，它的优点是非常容易入手，风格很像VB，相关控件可以随意拖动，GUI设计简单、省时，但GUIDE的一个严重缺点是无法直接创建核心对象；二是利用M函数构建GUI，即M文件界面设计，这种方法需要解决数据传递问题，如何正确实现回调函数中用户菜单或控件的句柄传送是M文件成功创建GUI的关键。事实上，不管采用哪种设计方法，事先都要分析界面所要求实现的主要功能，明确设计任务，并站在使用者的角度审查界面功能及界面的控件布局，然后进行代码编写，对功能进行逐项检查，调整完善界面功能。图形用户界面设计的一个基本原则要求具有简单性，即设计界面时应力求简洁、清晰地体现出界面的功能和特征，为此要尽量使用用户所熟悉的标志和符号，尽量删去可有可无的功能，尽量多采用图形结果，尽量减少窗口数目，力避在不同窗口之间进行来回切换。

2实例仿真及分析

非线性方程的迭代解法求根是数值分析课程的一个重要内容，初始迭代点及迭代函数的正确选取是求根的关键，为了使学生对迭代法求根有清醒的认识，下面以非线性方程迭代法求根的GUI实现说明MATLAB GUI对数值分析课程的辅助教学功能。

不动点迭代法求根中需要选取迭代公式，确定迭代初始点、精度，不动点迭代法求根的界面如图1。图1不动点迭代法求根的GUI界面

界面中设置了五个edit控件，分别用于输入方程f(x)、迭代公式、迭代初始点x0、精度tol和最大迭代次数；四个Push Button控件，分别用于绘图、求解、重设参数和退出界面；一个axes控件，用于显示函数f(x)的图像；为体现设计的简洁性，界面中只设置一个List? box控件，所有的结果都将在Listbox控件中显示，这样设计使界面更加合理化。系统能输入任意的方程，通过huatu_pushbutton3控件得到其图像，很容易判断该方程在零点的大致位置，即迭代初始点x0。输入方程后，单击画图控件，可以得到函数f(x)的图像，并显示在界面中。用编制好的GUI演示求解程f(x)=x3-x-1=0在x0=1.5附近的根x*，并用两种迭代公式求根，迭代公式分别为x= 3，初始点x0=1.5，精度tol=0.000001,最大迭代次数N=20，左键单击不动点求解控件，得到求根运行界面，如图2。在运行界面中得到运行结果，并且在函数图像中标出了通过运行得到的方程的根。

选取迭代公式x=x3-1，初始点x0=1.5，精度tol=0.000001,最大迭代次数N=50的运行界面，左键单击不动点求解控件，得到如图3的求根界面。图2,图3分别为同一方程取不同迭代函数求根的运行界面，由此可以让学生直观的看出不动点迭代法求根在选取不同迭代函数时，得到的收敛效果不同，直观的体现了迭代函数的重要性。

用Newton法来求方程f(x)=x3-x-1=0在x0=1.5附近的根，精度tol=0.000001,最大迭代次数N=20。编制的GUI演示结果可以让学生感受到Newton法求根的收敛速度比不动点迭代法求根的收敛速度快。

学生通过以上非线性方程求根的GUI，很容易体会到不动点迭代求根选取迭代函数的重要性及不动点迭代与Newton法求方程根的区别。同时设计的GUI具有开放性，可以让学生课后添加控件与代码，实现GUI更多的功能，这样不仅能够提高学生对数值算法的理解，而且极大提高学生学习数值分析课程的兴趣及编程解决实际问题的能力。

3结束语

将MATLAB GUI与数值分析课程结合起来,教师可以现场演示数值方法，开阔了学生学习数值分析课程的思路。若针对数值分析课程的所有教学重点内容编制一个辅助教学仿真软件，这对于数值分析课程的可视化教学、学生的数值实验更有意义。

参考文献：

[1]张志涌.精通MATLAB [M].北京:北京航空航天大学出版社,2000.

[2]陈垚光,王正林,毛涛涛.精通MATLAB GUI设计[M].北京:电子工业出版社,2008.

[3]尚涛,石端伟,安宁,等.工程计算可视化与MATLAB实现[M].武汉:武汉大学出版社,2002.

篇10

关键词: 改善; 模拟；气流组织; 分析; 比较;

中图分类号： F407.474文献标识码： A

1 机舱物理模型

机舱内设备及管线众多，结构复杂，建立模型时必须进行简化来使模型适合数值模拟计算，模型简化的基本原则是: 对空间气流场产生作用较小且温度影响不大的设备管线进行移除; 对结构形状复杂，简化后对模拟结果影响较小的设备规则化; 对厚度边界做无厚度壁面处理; 忽略风管，在风口段用立方体代替，底面为风口;抽风围阱对流场影响小，以风口面代替等;该船舶机舱模型空间大小为: 长 21m，宽 13m，高 10m，主要设备有汽轮机、锅炉、涡轮增压机组、高温除氧器、等离子过滤器、部分烟气管道等；该舱底层为铺板层，沿高度方向分别为第 1、2 格栅层，顶部为甲板层，风口均布置在第1、2 格栅层下方，为防止送风直吹发热表面，造成较大热应力，风口均采用矩形风口，风向垂直向下，尺寸大小为0.4m×0.2m 的风口55个，尺寸0.4m×0.28m 的风口 4 个 ( 侧 45°方向增加局部绕流) 。机舱简化后几何模型如图 1 所示。

图 1 机舱几何模型

2 流体数学模型及边界条件

2.1 流体数学模型

CFD 技术经过长时间的发展，在计算流体流动传热传质方面已经非常成熟，本文通过 Fluent软件CFD技术模拟机舱内气流组织分布，机舱内通风为湍流对流换热问题，本文选用图1模型，为获得较好的模拟结果，提高模拟效率，对模拟进行如下假设和简化:

1) 机舱内流体为不可压流体，密度符合假设;

2) 舱内流体流动及传热视为稳态过程;

3) 各壁面的辐射传热不计;

2.2 边界条件设置

利用开局让棋法进行几何建模及划分网格，因机舱内结构复杂，模型网格采用非结构化网格，对风口及重要发热面进行局部加密，网格数量为 88 万.

舱内通风是采用典型的置换通风与自然通风相结合的形式，边界条件设置如下:

1) 送、回风口边界条件: 送风及机械抽风口为速度入口边界，送风方向与出风口面垂直；

2) 自然排风口边界条件: 通风时保证舱内正压，因此自然排风口采用压力出口边界；

3) 固体壁面边界条件: 舱内发热固体壁面采用定热流密度，壁面厚度为0.

3 设计方案及模拟结果分析

3.1 设计方案

机舱内送风口位置相对固定，采用下送上回的通风方式，回风进入抽风围阱后通过管道与各风机相通；本文设计锅炉燃烧空气量占总送风量50% 左右，从舱内吸入会造成内部各区域气流压力梯度过大，因此采用独立系统从外界抽进燃烧空气。根据《ISO8861-1998造船--柴油机船舶机舱通风设计要求和计算基准》进行热负荷及通风量估算，得到该通风系统排除舱内余热所需风量为160000m3/h，舱内设计温度54℃，送风温度39℃，风机开度 100%时，送风口的设置送风速度为9.4m/s，考虑到送风量对气流组织及舱内温度的影响，设定送风开度来进行分析，分别为 100%，85%，70%.机舱环境的优化可通过改变送抽风口位置形状、面积及送风量等来实现，为获得最佳气流组织形式，该舱通过增大抽风围阱及隔屏来增加舱内空气扰流，防止大漩涡及通风死角，理论上，机械抽风口分布越多，抽风围阱尺寸越大，会有效改善舱内气流组织，但是机舱内空间狭小，因风口而增加的风管会减少舱内可用空间; 且过多地增加隔屏会给工作人员管理维护带来不便，经过分析给定 4 种调整方案来进行研究，见表 1

3.2 模拟结果分析

3.2.1 气流组织分析

计算收敛后，选取典型截面进行分析，由图 1可以看出，两锅炉间至主机右侧面区域散热面积最大，流场也最复杂，因此截取y=﹣0.5m 面进行气流组织分析，y=﹣0.5m 截面如图 2所示.

在假定机舱送风口位置固定情况下，要消除漩涡和通风死角，一般采用下列方法: 一是设置障碍物，障碍物对室内气流组织有显著影响; 二是布置空气射流通风系统，由于气流至底层后从下至上流动，障碍物迎风面不存在滞留区，同时考虑经济成本问题，对主机右侧安装隔屏作为障碍物来消除此漩涡，隔屏位置如图 2 所示"

4 结语

模拟结果表明，利用CFD技术对机舱气流组织进行模拟分析是一种高效简洁的方法，获得的结果是可信的，且通过模拟分析可知:

1) 合理的布置机械抽风口能有效地改进舱内流场，在舱内增加隔屏作为障碍物并非完全是负面影响，在不影响舱内设备运行和人员工作情况下，增加隔屏能有效改变舱内气流组织

2) 自上而下的通风方式能使各工作面温度得到较好的控制，人员活动区域温度符合设计要求

3) 风机开度直接影响舱内温度，通过改变风机开度能节省通风资源，机舱内环境的优化有很多方式，可通过改变风口数量大小、形式、位置等，这些还需要进一步的研究和补充

参考文献:

［1］向立平，王汉青空调客车内气流组织与污染物浓度场数值模拟［c］中南大学学报( 自然科学版) ，

篇1

篇2

篇3

篇4

篇5

篇6

篇7

篇8

篇9

篇10

热门标签

相关文章

相关期刊

数值计算与计算机应用

Numerical Mathematics

高等学校计算数学学报

护理与康复

精品范文