时间：2023-05-31 09:32:42
开篇：写作不仅是一种记录，更是一种创造，它让我们能够捕捉那些稍纵即逝的灵感，将它们永久地定格在纸上。下面是小编精心整理的12篇数值计算，希望这些内容能成为您创作过程中的良师益友，陪伴您不断探索和进步。
英文名称：Journal on Numerical Methods and Computer Applications
主管单位：中国科学院
主办单位：中国科学院计算数学与科学工程计算所
出版周期：季刊
出版地址：北京市
语
种：中文
开
本：16开
国际刊号：1000-3266
国内刊号：11-2124/TP
邮发代号：2-413
发行范围：国内外统一发行
创刊时间：1980
期刊收录：
中国科学引文数据库(CSCD―2008)
核心期刊：
中文核心期刊(2008)
中文核心期刊(2000)
中文核心期刊(1996)
中文核心期刊(1992)
期刊荣誉：
联系方式
Abstract： The determination of the vertical bearing capacity of prestressed concrete pipe pile is a concern in the engineering field. There are many academic researches on the vertical bearing capacity of prestressed concrete pipe pile. At present， there are few abroad empirical formulas about the vertical bearing capacity of prestressed concrete pipe pile， and there is no accurate calculation method in China. If the vertical load capacity of the single pile is determined according to the empirical formula of the current code， the result is often much lower than that obtained by the static load test， which results in a large increase in the cost of the project. In this paper， combined with the project construction in Panjin， Liaoning， the static load test， theoretical calculation， numerical simulation and other methods are used to analyze the stress characteristics and the unique bearing mechanism of prestressed concrete pipe pile.
P键词： 预应力管桩；单桩竖向承载力；静载荷试验；数值分析
Key words： prestressed pipe pile；vertical bearing capacity of single pile；static load test；numerical analysis
中图分类号：TU473.1 文献标识码：A 文章编号：1006-4311（2017）06-0159-03
0 引言
预应力混凝土管桩具有成桩质量易控制、施工简便，单桩承载力高，工程造价低，节能、环保等诸多优点。目前在许多地区已得到越来越广泛的应用。但由于开口管桩结构的特殊性，其沉桩和承载机理较为复杂，国内外许多学者对此做了大量的研究。鉴于此，本文针对盘锦地区的一个工程实例，依据沉桩施工资料，对预应力管桩静载荷试验结果进行了具体分析，并讨论了产生此结果的原因，为其他类似工程的设计施工提供技术参考。
1 工程概况
本施工场地位于辽宁盘锦，大地构造位置处于新华夏第二沉降带，堆积了厚达6000～8000m的新生界陆相地层。地貌单元属辽河河口三角洲，地势平坦，地貌单一，地层主要由第四系全新统海陆交互相沉积物组成。各土层统计物理力学指标、桩基参数指标如表 1所列。
2 静载荷试验
①1#、2#桩：当荷载分别加到1800kN时，累计下沉量分别为12.47mm、12.46mm。当荷载分别加到1890kN和1880kN时，桩急速下沉，荷载已经加不上去，千斤顶已自动卸载，桩已丧失承载力，终止加载。此时1#、2#桩的累计总沉降量分别为45.75mm和46.26mm，如表2-1。单桩竖向抗压极限承载力，取Q-S曲线陡降段前一级荷载为1800kN，如图1和图2。
②3#桩：当荷载加到1800kN时，累计下沉量为12.13mm。当荷载加到1990kN时，桩急速下沉，荷载已经加不上去，千斤顶已自动卸载，桩已丧失承载力，终止加载。此时3#桩的累计总沉降量为47.41mm。单桩竖向抗压极限承载力，取Q-S曲线陡降段前一级荷载为1900kN，如图3。
3 单桩竖向承载力计算
按规范中给出的参考数据合理取值。通过计算得到的计算结果均低于静载荷试验值。计算结果如表2。
4 数值模拟分析
4.1 地基特性
除了管桩桩身混凝士材料外，桩周有多层土体，所以在定义材料属性性窗口中，定义多个土层属性。在该对话框的窗口中，定义各种材料的弹性模量、饱和容重、泊松比、粘聚力、内摩擦角等。
4.2 材料特性
不论是二维还是三维计算模型，都需要有一个合理的网格划分方法和网格密度把握，桩土结构涉及到的模型几何形状较规则，因而采用合理的网格划分方式可以使计算来得方便。桩土材料力学性能见表3。
4.3 荷载和边界条件
根据现场预应力管桩实测资料管桩可以承受1800kN竖向荷载，本模型对管桩桩顶施加1800kN压力，并对桩顶面所有节点进行耦合，使桩顶面集中力转化为均布荷载，荷载分13级进行加载，每级加载140kN。对模型边界进行约束Ux=0，Uz=0。
4.4 理论计算、有限元分析与现场载荷结果分析
为验证数值模拟的结果是否能够正确反映管桩桩身荷载传递规律，选取现场的静载荷试验数据，并用有限元对现场情况进行数值模拟。图5为l#、2#、3#桩的静载荷数据曲线和有限元模拟静载荷试验曲线的对比情况。
在加荷的初始阶段，沉降值与实际值较相近。实测情况的最后阶段，荷载达到一定值时，沉降值会有一突然增大的现象，数值会变得非常大，这表明桩土在这一时刻的平衡关系被打破，桩体承载力达到极限。由此可以看出，有限元数值分析在实际工程中有着较好的实用性。
5 结果对比分析
静载荷试验过程与勘察报告中所反应的土层的力学性状基本一致。根据根据土的物理指标与承载力参数之间的经验确定预应力混凝土管桩的单桩承载力时，计算值和现场实测值较接近，且偏于安全。预应力管桩属于端承摩擦桩，桩身承载力较多的依靠侧摩阻力提供。在试桩施工过程中，因沉桩时间很短，桩侧阻力发挥作用较小，静载荷试验反映出来的压力值主要来自于桩端阻力，其侧阻力的发挥较少，桩的极限承载力还没有完全发挥出来。
参考文献：
[1]JGJ94―2008，建筑桩基技术规范[S].
[2]JGJl06―2003，建筑基桩检测技术规范[S].
[3]施峰.PHC管桩荷载传递的试验研究[J].岩土力学，2004，26（1）.
引言
微处理器的工作过程是大量数据的输入--运算--输出的过程，其中相当数量的数据使用十进制形式表达。使用者希望微处理器的输入数据和输出结果能使用十进制形式表达，而在微处理器内采用二进制表示和处理数据更方便，所以在二者之间的数制转换是必要的。通常采用两种方式解决这一问题。
方法1：十--二进制转换电路将输入的十进制数据转换为相应的二进制数据，微处理器内部算术逻辑单元仍然执行二进制数据运算微操作，运算结果再进行二--十进制转换，将结果以十进制形式输出。
方法2：算术逻辑单元对二进制数据处理能力的前提下，增加少量硬件线路，使之对某种二进制编码形式表示的十进制数据具有直接处理能力，该算术逻辑单元能够接收特定二进制编码构成的十进制数据，可以产生相同编码组成的计算结果，在数据处理过程中该单元执行十进制数据运算微操作。
微处理器使用中涉及大量的数据输入输出操作，显然方法1不是理想的选择，因而从提高机器的运行效率，简化机器结构和保证系统时序结构的规整性考虑，方法2更有实用价值。 所以本文讲述了方法2为算法依据的BCD加减电路。
校验原理
在计算机得数值计算中，数值经常是以BCD码表示的十进制进行运算的。即一位BCD码用4位二进制位表示。但是BCD的加法需要两个加法器来完成，如果分析一下BCD数的加法过程，原因就很清楚。请看下面：
令A=1000,B=0111,这两个数都是正确的BCD码，如果两个操作数直接相加，结果不是一个BCD码：
1000
＋ 0111
1111
正确的BCD码加法运算应为1000＋0111＝（1）0101即8＋7＝15。其它BCD码操作数运算的结果也能产生不正确的BCD码结果。实际上当结果大于9或者有进位时，就要进行BCD的校验，以确保结果的正确性。
对于产生进位得情况，加法器直接提供了二进制的进位输出，即BCD修正信号Y=C.而对于结果大于9，需要修正的数为1010－1111。
把它们作为四变量布尔表达式的最小项，就能化简逻辑。即Y=E3E2+E3E1.其中E3 、E2、E1、E0是加法器的和的输出。综合以上结果可得BCD修正信号Y=E3E2+E3E1+C.修正电路如图一所示
下面就已四位并行加法器和一位串行加法器两种电路形式来讨论BCD码的验证。
图二 4位并行加法器BCD加法电路
图二所示为4位并行的BCD加法器电路。其中上面加法器的输入来自低一级的BCD数字。下面加法器BCD的输出E3、E2、E1、E0和COUT至高一级BCD数字,其A3和A1位接地,即当BCD校验信号为真时Y=1,A3A2A1A0=0110,以实现加6的调整.当不需要BCD调整时Y=0,此时A3A2A1A0=0000,从而使输出结果无变化.
虽然4位并行加法器运算速度较快，但是所用逻辑门较多。图三所示为一位串行BCD加法器。它是以牺牲速度以达到减少硬件逻辑门的目的，这种电路在对频率要求不高的系统中非常之适用。其中ADDER1、ADDER2均为一位全加器。ADDER1做主运算器,ADDER2做BCD校验运算器，不管是否做BCD校验，ADDER2的初始进位、借位始终为“1”。
图三中Z型门为延时电路，延时一个时钟周期，这样在外部电路控制下，经过四个时钟周期，得到一位十进制BCD结果E3E2E1E0.由电路图所以当C+(E3E2+E3E1)逻辑值为‘1’时，控制多路选择器选择A通路（A通路为序列1001），当C+(E3E2+E3E1)为‘0’时，选择B通路（B通路序列为1111），即需要校验时，多路选择器输出序列1001；不需要校验时，输出序列1111，与Z型门的输出对应相加，并且ADDER2的初始进位始终为‘1’，由此可完成BCD的校验工作。
图三 一位串行BCD加法器电路
下面是基于4位并行BCD加法器算法的一种快速BCD的加法器VERILOG硬件描述语言程序及其仿真结果。
module bcd_check (data_i,data_o,cy_i,cy_o,en,z_i,z_o);
input data_i;
input cy_i;
input z_i;
input en; //insructure
output cy_o;
output data_o;
output z_o;
wire [3:0] data_i;
wire cy_i;
wire en;
reg z_o;
reg cy_o;
reg [3:0] data_o;
//}} End of automatically maintained section
reg [4:0] TEMP_RESULT;
always @(data_i or cy_i or en )
if(en == 0)
begin
cy_o=cy_i;
data_o=data_i;
z_o=z_i;
end
else
begin
if(data_i[3]&&data_i[1] ||(data_i[3]&&data_i[2]) || cy_i==1)
TEMP_RESULT = {1'b0,data_i } + {1'b0,4'b0110 } + cy_i;
else
begin
TEMP_RESULT[3:0]=data_i;
TEMP_RESULT[4]=cy_i;
end
data_o = TEMP_RESULT[3:0];
z_o = | TEMP_RESULT[3:0];
cy_o = TEMP_RESULT[4];
end
例如：两个十进制数2189＋8075的正确结果应为11064，可是，相加运算后的结果为FEH,为此应进行BCD调整。将为经校验的相加结果0010，1001，1000，1001（十进制2989）＋1000，0000，0111，0101（十进制8075）＝1010，1001，1111，1110代人上述BCD校验模块，可得仿真结果如图四。
图四 仿真结果
有图可知data_o为1064 且进位输出为1，即总的结果为11064，这与2989＋8075＝11064的结果是一致的。
关键词：基坑；Blum法；数值模拟
中图分类号：TU47 文献标识码：A
1 前言
基坑工程是一项古老的综合性工程，支护方案的设计计算方法影响到支护方案的工程进度、工程质量和工程成本，在整个基坑施工工程中占有相当重要的地位。随着我国经济的不断发展，城市土地资源越来越紧张，地下空间的发展利用变得日益迫切这就对基坑工程的支护设计提出了更为严格的要求。同时由于基坑工程的地域性强、发生事故后损失大、补救困难且造价高等特点，就要求基坑工程的支护方案必须稳定可靠、经济合理选择基坑支护计算方案显得尤为重要。
2 基坑概况及场地周边环境
基坑概况。长春农贸集团股份有限公司商场大厦占地面积约为10000m2，基坑开挖深度5.8m，商场底部为一地下车库，占地面积5500m2。
场地周边环境。长春农贸集团股份有限公司商场大厦位于黑水路批发市场西侧10m处；其北侧为距其10m的长春市绵麻土特产品总公司（6层）和长春市农贸集团股份有限公司办公用楼（6层）；西侧为距其26m的东三条街；南侧为距其13m的黑水路。该建筑基坑深度为5.8m。场地的四层土层的物理力学性质指标见表1。
3 采用Blum法计算参数进行软件数值模拟分析
3.1 几何模型。根据本基坑的实际尺寸，取土层边界为基坑宽度的一倍，深度为基坑开挖深度的二倍。有限元计算中，用梁单元模拟支护桩；几何模型如图1所示，围护结构参数同前所述。
3.2 网格的划分。Midas/GTS软件提供有很多种单元形式，本文对于基坑周围土体采用的是四边形平面单元形式，维护结构采用的是直线1D单元形式，支护方案的网格的划分结果如图2所示，划分9351个单元，9562个节点。
3.3 边界条件。根据本基坑工程的特点，对位移边界条件做以下假定：模型的竖直方向允许发生变形，水平方向的左、右边界位移为零；下边界任意方向的变形都为零，边界条件如图3所示。桩的参数：桩的长度为10.03 m，直径为0.6m，桩的弹性模量为2.5×105N/mm2，泊松比0.2，容重为25 N/mm3。土层的参数如表2所示模量为2.5×105N/mm2，泊松比0.2，容重为25 N/mm3。土层的参数如表2所示。
3.4计算结果
根据上述方法的数值模拟表明Blum法安全系数2.9203mm，数值模拟结果最大水平位移1.66mm和最大竖向位移1.07mm，说明此方法可靠度较高，所以本基坑实际工程采用Blum法计算的桩长。
结论
随着越来越多高层、超高层及重大型工业建筑的发展，基坑开挖的深度日益加深，排桩支护结构的设计、施工等方面都面临很多问题，只要我们能认真进行方案的优选、方案的论证，设计理论不断完善、不断创新，施工工艺不断改进，基坑工程会得到更快更好的发展。
参考文献
[1]余志诚，施文华.深基坑支护设计与施工[M].北京：中国工业出版社.1997.
[2]徐杨青，王永宁，程杰林.模拟深基坑开挖和支护全过程的有限元数值分析[J].岩土力学，2002（S1）.
[3]谢猛，侯克鹏，傅鹤林，等.值梁法在深基坑支护设计中的应用[J].土工基础，2008.
随着科技的迅猛发展，越来越多的人意识到机器的计算能力在现代科学发展过程中举足轻重的影响力。现如今各领域的工程师也不遗余力地优化各自的计算模型和工具，从而让项目的计算性能得到进一步提升。而近年来，随着图形处理单元（GPU，Graph Processing Unit）的发展日益成熟，其应用的范围已从最初的计算机图形学领域，逐步扩展到数值计算领域。与传统的CPU相比，GPU在计算能力方面有着巨大的优势。本书主要关注GPU中数值方法的实现，着重介绍了目前GPU在数值计算领域的发展以及相关原理和方法，
全书共18章，分为4个部分，每部分都包含了一些典型的数值计算方法。第1部分 线性代数问题求解，包括第1-6章：1.基于GPU结构的稠密线性代数计算；2.基于GPU结构的三对角计算；3.矩阵计算，介绍LAPACK、GEMM以及MKL等高性能的线性代数计算库；4.LU和QR分解，介绍GPU编程中批量分解算法；5.线性系统下基于LU分解的CUDA算法，着重阐述小矩阵计算中批量线性求解方法的CUDA实现；6.矩阵向量乘积，着重介绍如何用GPU高效地实现科学计算中经典的稀疏矩阵与向量乘积（SpMV）算法。第2部分 介绍微分方程的空间离散化，包括第7-11章：7.GPU常微分方程的求解；8.GPU并行集成，介绍GPU中线程级并行算法的实现；9.介绍谱元法在非结构网格流动和波动问题中的应用；10. 利用局部修正SOR方法求解对流扩散问题；11.CUDA及OpenCL编程中的有限差问题。第3部分 随机数和蒙特卡洛方法，包括第12-15章：12.GPU伪随机数生成，介绍蒙特卡洛仿真算法在OpenCL下的编程实现；13.蒙特卡洛自动化积分，介绍CUDA编程中动态的并行计算问题；14. GPU加速计算实例：量子轨迹法；15.GPU动态系统中的蒙特卡洛仿真，介绍数值积分并给出蒙特卡洛方法的应用实例。第4部分 快速傅里叶变换以及N体问题，包括第16-18章：16.快速傅里叶变换（FFT）在GPU下的编程实现；17.介绍一种共享内存复用方法，在此基础上实现高效快速的傅里叶变换；18.N体问题仿真算法在GPU中的实现，探讨在增加并行性的同时如何减少对线程等系统资源的抢占问题。
本书内容丰富，而且特别注重实际应用，可作为GPU以及CUDA的学习和编程参考书。对从事高性能计算相关学科学习的高年级大学生、研究生和相关领域的研究人员，本书极具参考价值。本书要求读者有一定的线性代数、微分方程等数学功底且具有扎实的C或C++语言编程基础。
关键词：大跨度桥梁 模态分析 数值计算步骤
1、引言
寸滩长江大桥为主跨880m的双塔钢箱梁悬索桥，边缆跨度250m，北塔高199.5m，南塔高194.5m，桥塔梁上部分高度为117m，钢箱梁宽39m，塔宽39m，其他资料详见图纸说明。对此桥进行模态分析，为描述结构振动特性及减小振动对结构的不利影响提供重要参数，如模态频率、模态振型、模态阻尼等。
2、模态分析数值计算的一般步骤
①做好资料准备工作，了解工程概况，研读图纸和设计资料，记下桥梁的各项与模态分析相关的结构尺寸和材料性能；②选定合适软件进行数值分析，一般常用ANSYS等有限元软件；③根据所需结构尺寸和材料性能，通过有限元软件进行数值建模；建模过程一般先进行各部分节段划分，然后确定主梁、桥塔、主墩、悬索或拉索、锚固点等坐标，再计算各划分截面的实常数（需指出，采用不同的模拟单元，实常数的形式不一样），然后通过实常数建立单元，最后再将没有建入模型的部分如横隔板、风嘴、二期等可以通过加质量点的方式考虑进去，以便模拟更精确；④对模型施加约束；塔底、锚固点一般是固结约束，塔梁结合处一般采用CP命令进行耦合，耦合自由度依据设计说明；⑤最后再依次进行恒载静力分析和模态分析，提取各阶振动频率和模态。
3、寸滩长江大桥动力特性分析实例
3.1 动力特性分析资料及采用的模拟单元类型
塔柱、主梁和桥墩等均采用梁单元模拟；主缆和吊索均采用杆单元模拟；二期恒载采用质量单元模拟。所采用的边界约束条件、单元类型如表1、2所示。
表中：x为纵桥向，y为竖向，z为横桥向。0表示自由，1表示主从，d表示固结约束。
3.2 建立有限元模型
3.2.1 节段划分和坐标确定
依据设计和图纸资料，主梁每5m划分一段，桥塔在特殊截面位置（如与横梁、大缆连接等位置）需要单独划分，在一般位置同样每5m划分一段，大缆在吊杆位置划分，然后定出各部分的坐标，坐标原点可以任意选取，本例中坐标原点定在主跨跨中主梁截面底部位置，在建立各部分坐标的时候尤其要考虑主梁的纵向坡度1.5%。
3.2.2 实常数计算
主梁为钢箱梁截面，采用BEAM44单元模拟，需要简化截面，将风嘴、横隔板、二期等部分删除，然后在CAD中建立面域，通过面域massprop查询截面特性找到形心位置，将面域移动到形心位置，然后保存为*.sat文件导入ANSYS中，通过网格划分后，计算截面特性，CAD中查询的截面特性和ANSYS中计算的截面特性主要差别在于ANSYS中能计算出扭转惯性矩这个重要的实常数参数之一。桥塔为空心变截面，桥塔横梁为空心等截面，采用BEAM44单元模拟，在ANSYS中可以通过循环命令来实现截面特性的计算。大缆和吊索采用LINK8单元模拟，实常数形式不一样，其中初应变的确定需要调试，一般使跨中受力最小的初应变较为合适。质量点的实常数计算较为麻烦，需要计算质量惯性矩Im。
3.2.3 模型建立
通过实常数将各坐标点依次连接起来，依次连接北面桥塔、南面桥塔、桥塔横梁、主梁、大缆、吊索等。然后将简化后的二期、横隔板、风嘴等以质量的形式加载到各结点上。再镜像单元，对塔底和锚固点加约束、对塔梁进行耦合，形成完整的有限元模型。
3.2.4 提取模态分析计算结果
最后再依次进行恒载静力分析和模态分析，提取各阶振动频率和模态，提取前10阶结果见表3，并列出第1阶和第2阶频率对应的振型图见图2、3。
图1成桥状态第一阶振型 图2成桥状态第二阶振型
4、结语
通过寸滩长江大桥模态分析的实例可以清楚的了解到模态分析数值计算的一般步骤，掌握了一般步骤，对其他类型的桥梁进行动力特性分析时思路更明确，计算更快捷。同时需要了解模态分析的基本概念，有限元中模态分析的本质是求矩阵的特征值问题，所以“阶数”就是指特征值的个数。将特征值从小到大排列就是阶次。实际的分析对象是无限维的，所以其模态具有无穷阶。但是对于运动起主导作用的只是前面的几阶模态，所以计算时根据需要计算前几阶就能达到要求。
参考文献：
关键词 地层热阻 热响应测试地埋管换热器
中图分类号：TK172 文献标识码： A
Thermal resistance calculation on vertical ground heat exchanger
Li Jinghui1Wang Jiankui2 Lu Lin2 Fang Xugen1
1 Zhejiang Construction Division Building Energy Technology Co., Ltd.
2 Zhejiang Academy of Building Research & Design.ltd
Abstract: Using the data of grock-soil thermal response test, andaccordingto , this paper presents the numericalcalculationofthe vertical ground heat exchanger's thermal resistance calculation formula,andalso analyzeshow heat exchanger resistance affect the performance of heat exchanger. The calculation and analysis have a certain reference valueinguidingthe design and construction of ground heat exchanger.
Key word: ground heat resistance，heat responsetest，ground heat exchanger
0 引言
地层热阻是决定土壤源热泵合理设计和科学应用的核心因素，是影响地下埋管换热器传热性能及土壤源热泵系统节能与经济性的重要原因。换热器的热阻工程实际计算及热阻的影响对于设计人员设计及方案的优化有着重要的影响。
1 地埋管换热器的热阻计算
地埋管地源热泵系统地埋管换热器的设计计算主要有采用单位井深换热量、专业软件和依据《地源热泵系统工程技术规范》[1]（以下简称《规范》）作为目前指导地源热泵工程设计和施工最重要的规范，在其中第4.3.5条明确指出“竖直地埋管换热器的设计也可按本规范附录B的方法进行计算。
但是《规范》附录B中对地层热阻及短期连续脉冲负荷引起的附加热阻的计算公式的表述存在一定的歧义，也不完整，在实际应用中会产生较大的误差甚至错误[2]。本文以供热工况下的换热器的热阻为例各热阻值对进行计算分析。
1.1岩土热响应测试数据
表1 岩土热响应测试基础数据表
项目 测试孔 项目 测试孔
钻孔深度（m） 100 钻孔直径（mm） 135
埋管形式 单U型 埋管材质 PE管
埋管内径（mm） 26 埋管外径（mm） 32
钻孔回填材料 原浆 细沙 主要地质结构 粉质粘土
岩土热响应测试测试报告结果：埋管区域的平均综合导热系数为1.667 W/ m·℃，平均容积比热为2.034×106 J/m3·℃。岩土体初始温度20.1℃。
1.2值计算
公式1
式中：
——传热介质与U形管内壁的对流换热热阻，（m·℃/W）;
——U形管的内径（m）；
——传热介质与U形管内壁的对流换热系数[W/ m2·℃]。
换热器温差＜20℃，根据热物性测试报告定性温度按6℃，查表得=56.3×10-2W/ m·℃，ν=1.553×10-6m2/s ，=11.60 =1547×10-6N·s /m2
d：管内径=26mm管内媒质流速按0.9 m/s计算；=0.9×0.026/1.553×10-6=15067＞104；=0.023×15067×11.6=4020；K==4020×0.563÷0.026=87045 W/ m2·℃；=1÷(3.14×0.026×87045)=1.41×10-4 m·℃/W
1.3 值计算
公式2
公式3
式中：
——U形管的管壁热阻，（m·℃/W）;
——U形管导热系数[w/(m·℃)]；
——U形管的外径(m)；
——U形管的当量直径(m)；对单u形管，n=2；对双U形管，n=4。
=0.42~0.45 w/(m·℃)；=0.032 m； ——0.0452m
=0.0539 m·℃/W
1.4计算
公式4
式中：
——钻孔灌浆回填材料的热阻(m·℃/w)；
——灌浆材料导热系数[w/(m·℃)]；
——钻孔的直径(m)。
=1.47w/(m·℃)， =0.135m，= 0.0809m·℃/w
1.5计算
对应于单个钻孔：
公式5
公式6
式中：
——地层热阻(m·℃/w)；
——指数积分公式；
——岩土体的平均导热系数[w/(m·℃)]；
——岩土体的热扩散率m2/s；
——钻孔的半径（m）；
——运行时间（S）；
根据文献[2-4]对进行修正简化为下式：
公式7
= 公式8
公式9
=1.667w/(m·℃)，=/cρ，cρ=2.034×106J/m3℃，则=0.8196×10-6m2/s，钻孔的半径=0.0675（m）,运行时间=68×24×3600=5.8752×106（S）；z=0.0006
==0.3712 m·℃/w
1.6计算
公式10
式中：
——短期连续脉冲负荷引起的附加热阻(m·℃/W)；
——短期脉冲负荷连续运行的时间（s）。
根据文献[2-4]修正简化为下式：
公式11
=10×3600s
=0.1297 m·℃/w
2 热阻影响分析
根据《规范》中公式B.0.2-4中供热运行份额Fn取0.2。根据以上计算得各热阻值及比例见下表，竖直埋管换热器分项热阻计算分析。
表2 竖直埋管换热器分项热阻计算分析表
项目 ×Fn ×(1-Fn)
热阻值
（m·℃/w） 0.001 0.0539 0.0742 0.1038 0.1184
份额（%） 0.04 15.39 21.18 29.60 33.79
根据以上数据可知，单孔U型管回填材料的热阻、地层热阻和短期连续脉冲负荷引起的附加热阻为构成换热器热阻的主要部分，U型管内壁的对流换热热阻很小，对换热器热阻影响几乎可以忽略。
3 结论
3.1对竖直地埋管换热器各热阻值进行了详细的计算，结合相关文献对地层热阻及短期连续脉冲负荷引起的附加热阻进行了简化计算。
3.2 计算结果表明，单孔U型管回填材料的热阻、地层热阻和短期连续脉冲负荷引起的附加热阻为构成换热器热阻的主要部分，对于优化换热器设计方案具有一定的参考价值。
参考文献
[1] 中国建筑科学研究院.GB 50366-2005 地源热泵系统工程技术规范 [S].北京：中国建筑工业出版社,2009版
[2] 雷建平,於仲义. 关于《地源热泵系统工程技术规范》地层热阻计算式的商榷及应用分析[J], 暖通空调,2009,39(6):27-30
[3] 余延顺. 土壤耦合热泵系统地下埋管换热器传热模型的研究[J].暖通空调,2005,35(1):26-31
关键词：暂态稳定性；数值积分方法；RungeKutta方法；RungeKuttNystrm方法；RKNd方法
中图分类号：TM744 文献标识码：A
Fast Numerical Simulation of Power System Transient Stability by RKNd Methods
ZHANG Lei，WANG Fangzong，HU jiayi
（Electrical Engineering & Renewable Energy School，China Three Gorges University，Yichang443002，China）
Abstract：The RKNd method is a new kind of numerical integration methods， of which the order is higher than that of the traditional RungeKutta methods and RungeKuttaNystrm methods for the same stage. In this paper， the RKNd method is introduced to the numerical simulation of power system transient stability， and then a fast numerical simulation method has been proposed. The proposed method has been compared to both the traditional numerical integration method and the symplectic Gauss method using IEEE145bus power system， and the tested results show that the implicit RKNd method has the advantages both in calculation accuracy and in computational efficiency respectively over the symplectic Gauss method and the implicit trapezoidal rule. Therefore the proposed methods should be more suitable to numerical analysis of transient stability and other likewise problems.
Key words：transient stability；numerical integration method；RungeKutta method；RungeKuttNystrm method；RKNd method
1引言
数值积分方法是电力系统暂态稳定性分析计算的基本方法。最常用的数值积分方法大致包括隐式梯形积分法以及RungeKutta方法（RK方法），前者是隐式积分类方法，后者是显式积分类方法。
近年来，研究人员又提出了不少新的数值积分算法。文献[1]和文献[2]将辛Runge—Kutta算法（辛RK方法）、文献[3]将可分Hamiltonian系统的显辛算法、文献[4]将辛代数动力学算法用于暂态稳定性的计算，并对这几种新的数值积分方法进行了测试和对比分析。文献[5]和文献[6]分别将多级高阶辛RK算法以及多级高阶辛RungeKuttaNystrm算法用于暂态稳定性的并行计算。
文献[7]利用一阶常微分方程导出的二阶方程，借鉴Nystrm方法，提出了一类新的数值积分方法，即RKNd方法。RKNd方法的最大优点是：在相同级数情况下，RKNd方法可达到的最高代数阶比传统的RK方法高。在传统的RK系列方法中，s级的显式RK方法可达到的最高阶数是s阶；s级的隐式RK方法可达到的最高阶数是2s阶。但RKNd方法不同，2级的显式RKNd方法可以达到4阶；2级的隐式RKNd方法可以达到5阶。因此，与同级的RK方法相比，RKNd方法具有更高的计算精度；与同阶或略低阶的RK方法相比，RKNd方法具有更高的计算效率。
本文将RKNd方法引入电力系统暂态稳定性的数值计算。以IEEE 145节点系统为例，分别将2级4阶显式RKNd方法与传统的4级4阶显式RK方法、2级5阶隐式RKNd方法与2级4阶隐式辛RK方法进行了对比测试。测试结果验证了RKNd方法在计算效率方面具有明显的优势，因而可以推广应用于电力系统暂态稳定性及其它领域的数值计算。
2RKNd方法简介
RKNd方法既不同于传统的RK方法，也与RKN方法有所不同。对给定的2阶常微分方程初值问题
=f（t，x）=g（t，x），x（t0）=x0 （1）
关键词：地下水数值模拟系统、矿井涌水量计算、实际应用
中图分类号：TU991.11文献标识码： A
正文：
矿井涌水量的准确预测对于防止矿井突水、淹井等恶性突发事故有着重要的意义，同时也能大大降低生产成本，保障矿山的安全生产。【1】
一、地下水数值模拟系统
（一）概况
地下水数值模拟系统主要用来解决各类水文地质问题，并起到一定的预测作用，其模拟任务主要有四种，分别为地下水运移模拟、水流模拟、反应模拟以及反应运移模拟。在模拟模型的简历过程中，需要针对其中某一个目标，模型的建立步骤一般有这样几点：1、建立概念模型；2、选择数学模型；3将数学模型进行数值化；4、模型校正；5、校正灵敏度分析；6、模型验证；7、预测；8、预测灵敏度分析；9、给出模拟设计与结果；10、后续检查；11、模型再设计。
地下水数值模拟系统是随着计算机出现而发展起来的，以有限单元法为基本计算方法，在分割近似原理的指导下，将复杂的非线性问题简化为线性问题，从而避开了解析法求解微分方程时各种严格理想化的要求，使数值模拟系统更能灵活地适用于各种矿井涌水量的计算。【2】现目前常用的地下水数值模拟软件有美国开发的GMS（其中MODFLOW是世界上使用最广泛的三维地下水水流模型）、Visual MODFLOW、Visual Groundwater、PHREEQC、TNTmips等是运用较为广泛的软件。
二、实际应用
（一）郑煤盛源煤业有限公司概况
宝丰盛源煤业位于宝丰县大营镇宋坪村西南方，其由宝丰县大营镇宋坪村办煤矿和大营镇双鱼山二矿于2007年被郑煤集团整合而成,并于2010年加入中国有色金属工业集团。地处平顶山市宝丰县内，紧邻207国道，交通十分便利。矿井设计生产能力为30万吨每年，煤种为1/3焦，是优质的炼焦用煤。
1、可采煤层。
主要开采山西组下部的二1煤层以及一4煤层。
2、煤层标高。
二1煤层深埋280m~338m，煤层开采深度标高为-140m~0m;一4煤层深埋261m~400m，煤层开采深度地板标高为-160m~-30m。
（二）矿井的水文地质
该矿区的主要含水层有四个系层：寒武系上统崮山组，二又叠系下统山西组、下石盒子组，石炭系上统本溪组和太原组，第四系。
第一系层：寒武系上统崮山组含水层。白云质灰岩，厚59~131m，无泉水出露。
第二系层：二又叠系下统山西组、下石盒子组含水层。其中山西组是由二1煤层上部大占香碳砂岩段中的粗粒岩、砂岩段构成；下石盒子组是由下部中、粗粒砂岩组成，该层含砂岩裂痕承压水，富水性较差，不威胁煤矿开采。
第三系层：石炭系上统本溪组和太原组含水层，岩性为灰至深灰色结晶灰岩，含水层由L1~L8（L1灰岩厚为0.33~15.76m）薄层状灰岩和中粒砂岩组成。该岩层承压性较好，存储水量较多，水压较高，但其分布并不均匀。
第四系层：第四系含水层。由冲击岩、沙石岩组成，直接覆盖在下伏地层上，对矿井开采有一定的影响。
（三）模型的建立
概念模型的建立是一个极为复杂的过程，需要我们充分了解模拟地区的地质构造、水文地质、岩石矿物、气象、地形地貌、工农业利用等一切与地下水相关的关系点。【3】
在划分水文地质单元，确定模拟边界和范围之前，我们应该准备：地形地貌图、第四纪地质图、水文地质图、地下水等水位线及埋深图、模拟区遥感影像数据、有关的区域地下水方面的调查勘察研究报告及成果。
（四）模型的结构
首先应该对模拟区域的地质结构和水文地质条件加以概化，建立水文地质概念模型，然后构建相关数值模型。
1、有剖面线的位置的水文地质剖面图；
2、详细的钻孔、深孔资料（附有名称、坐标、孔标高、终孔深度、分层信息以及岩性描述等）；
3、以完整的水文地质单元作为模拟的区域范围，考虑到有些边界和范围过于偏远，应该考虑扩大模拟区的范围，适时可采用模型嵌套技术。
（五）模型的参数
1、潜水、承压含水层和弱透水层平、根据岩性和抽水试验分区的垂向渗透系数的分区图与数值；
2、承压水含水层存水率的分区和相关数值；
3、潜水含水层的导水程度的分区和数值；
4、弱透水层的存储率的分区图和数值；
5、各类渗水实验的资料和研究成果；
6、各层的有效钻孔的隙度。
（六）实际案例
1、确定一4煤层为计算区域，对计算区域进行三角分区和参数分区；
2、对一4煤层的含水系统内部进行概化，对其边界形态进行简化，对边界进水类型进行划分以及地下水的相关运动状态；
3、明确各分区的参数值，要详细了解一4煤层每个节点、每个分单元的信息，以及观察钻孔的水位信息和抽水孔水量信息，计算时段信息等。
在这个水文地质模型中将一4煤层共划分为20个分区，650个节点，分单元5512个，有效单元5082个，观测孔号5个。
表1 地下水位观测值拟合统计
其中我们选取了1号孔水位观测值进行研究，可得出：
表2 各分区参数计算
最后通过系统模拟可以计算出一4煤层的正常涌水量为222m²/h，最大涌水量为289m²/h。
（四）结语
近年来，随着科技的不断进步和发展，水文地质工作者们拥有了大量的科学决策和科学管理信息的方法。地下水数值模拟系统可以量化地下水的动态变化与人类开采活动的关系，可以比较不同开采方案并预测其对环境造成的影响。由此可见，地下水数值模拟系统在实践中是可以得到广泛运用的，它将在国家制定区域水政策以及各企事业单位开矿采矿事业中做出重大的贡献。由此可以预测，地下水数值模拟系统在矿井涌水量研究应用中的前景是无限广阔的，应该引起相关工作者的重视。
参考文献：
【1】魏军，《矿井涌水量的数值模拟研究》，2006年12月
关键词：大气科学；非线性计算不稳定；隐式格式
中图分类号：G642.0 文献标志码：A 文章编号：1674-9324（2013）42-0160-02
一、引言
根据物理定律，比如牛顿第二定律、质量守恒定律、能量守恒定律、气体试验定律等，可以得到支配大气运动的基本方程组。但由于大气运动基本方程组是一组高度非线性的偏微分方程组，很难求得其解析解，人们可通过数值的方法求其近似解（数值解），这就是数值天气预报。数值天气预报是一门实用性很强的应用基础学科[1]。通过差分方法求解大气运动基本方程组时，人们发现数值积分过程中会产生计算不稳定问题，这就需要采用恰当的差分格式或积分格式。在《数值天气预报》课程关于时间积分格式的讲解中[2]，重点介绍了“显示”格式，即差分方程的右端项全部为当前（或和过去）时刻的变量值，通过积数值分可求出方程左端的未来时刻变量值。但对于“隐式”及“半隐式”格式只是简单提及其概念和特点，比如隐式格式为用未来时刻变量值求出未来时刻变量值，它具有计算稳定、但计算复杂的特点。
在《数值天气预报》课程中，需讲解大量的公式推导和讲解，若不能配合简单而又形象的举例和图形，学生尤其是本科生作为授课对象，将很难理解和接受本课程中的相关内容，讲课的效果也将大打折扣。在非线性不稳定计算的举例中，教材中[2]虽给出了采用不同的初值和不同的差分方案对计算稳定性的影响，但并没有清楚地列出其求解过程，因此学生很难了解隐式格式差分的具体求解过程，对教材中列出的“显示”和“隐式”格式的各自优缺点更是难以理解。因此，需要对两种格式的计算过程进行相应的讲解，尤其是隐式格式。此外，教材中[2]对同一个微分方程构造的两个不同的差分方程中，除隐式格式和显式格式的差异外，还存在着对■采用了不同的差分格式，即显式格式采用中央差格式，隐式格式采用前差格式。本课程[2]已清楚的讲解到中央差格式虽具有较高的计算精度，但在时间差分计算时存在计算解的问题，若初值取得不当，则计算解会有较大的振幅。因此，从逻辑上讲，教材中给出的不同的差分方案的影响，实际上不仅仅来源于显式格式和隐式格式的差异，还来源于对时间微分采用不同差分格式的差异。这又加大了学生对显式格式和隐式格式特点的理解难度。
针对上述问题，本文将以简单的一维非线性平流方程为例，给出隐式格式差分方程的具体求解过程，重新探讨非线性计算不稳定现象，目的是使学生更好地了解显式格式和隐式格式差分方程的求解过程，深刻理解两种格式各自的优缺点。
二、非线性计算不稳定的计算实例
以大气科学中极具代表性的一维平流方程为例：■+u■=0，0≤x
■+■（■）=0，0≤x
或■+■（u■+■），0≤x
以上两式与教材[2]基本一致，不同的是这里的x取值范围并不到1。在大气科学中，方程或模式的计算可在全球或某一纬圈上进行。在该情况下，没有纬向侧边界条件。对于上式而言，可认为u在x=1的取值等于u在x=0的取值，也即循环边界条件。在构造上述微分方程相应的差分方程过程中，对（1）式和（2）式分别采用显式格式和隐式格式：
uin+1=uin-■[（ui+1n+uin）2-（uin+ui-+1n）2] （3）
uin+1=uin-■[（■i+1+■i+■i-1）（■i+1-■i-1）] （4）
其中上标n为第n步，下标i为第i个格点，■i=（uin+1+uin）/2。可见，与教材中不同的是，（1）式和（2）式中■均取了前差格式，这样可避免由于三个时间层计算而出现的计算解问题，有利于问题的讨论更加集中。
同样给定两种不同的初值，两者仅相差一个常数：
ui0=sin2πiΔx （5）
ui0=1.5+sin2πiΔx （6）
计算中，Δx取=1/3，Δt=0.004，则|u■|≤umax=|u■|=2.5×0.004×3=0.003
三、隐式格式的求解
显式差分方程（3）的求解过程即是将已知的n时刻u值代入等式右端算出等式左端未知的n+1时刻u值，可见，求解过程简单。至于隐式差分方程（4），其求解过程，较复杂。首先将（4）式写在[0，1）的x0=0、x1=1/3和x2=2/3三个格点上，并令m=-■，u0n+1+u0n=X，u1n+1+u1n=Y，u2n+1+u2n=Z，u0n=a，u1n=b，u2n=c。可见，a、b、c均为已知的第n步值。采用循环边界条件可得三元二次方程组：
X=2a+m（X+Y+Z）（Y-Z）Y=2b+m（X+Y+Z）（Z-X）Z=2c+m（X+Y+Z）（X-Y） （7）
将（7）式中的三式相加可得：X+Y+Z=2a+2b+2c，再令m（2a+2b+2c）=d，该d值也是已知的第n步值，（7）式可化为三元一次方程组：
X=2a+d（Y-Z）Y=2b+d（Z-X）Z=2c+d（X-Y） （8）
最终可利用已知的a、b、c和d值分别求得n+1步未知的Z、Y、X值：
Z=■Y=■X=2a+d（Y-Z）（9）
再分别将其减去c、b和a值，可得n+1步的u2n+1、u1n+1和u0n+1。由此可见，隐式差分方程的求解过程较为复杂。需指出的是，本文在[0，1）仅选取了3个格点，若选取教材中的10个点（Δx=0.1），则需在10个格点上写出10个差分方程，并进行联立，求解十元一次方程组，其求解过程更为复杂。
四、计算结果及分析
图1a和1b分别给出了两个初值、两种计算方案的计算结果。初值取（5）式用显式方案（3）式的计算结果表明（图1a实线），动能逐步增大，在500步以后突然急剧增加，出现按指数增加的趋势；但若给初值加上一个常数后（图1b中实线），总动能在4m2/s2左右变化，表明计算结果稳定。至于隐式方案（4）式，无论取哪种初值，结果均稳定。
至于产生如图1a中的不稳定现象，仍可利用混淆误差理论进行解释，即网格系不能正确分辨短波长的波动而导致不稳定。本文例子在[0，1）的一个周期范围内仅取三个格点，采用循环边界条件，即将[0，1）进行I=3等分。可见，该网格系只能正确识别平均值0波、波长为的3/2波和波长2Δx为3Δx的1波波动。若波数k1=k2=3/2的两个波动相互作用，则可产生0波和3波的波动。其中3波波动超出该网格系的识别能力，将会被错误的识别为0波。该0波的能量将不断的积累，从而可导致不稳定现象。该过程也可通过（10）式得到验证：
sin■=sin■=0=sin■cos■=cos■=1=sin■ （10）
因此，虽然本文的计算结果与教材中基本一致，但举例十分简单，这有利于学生的理解和接受。
五、结束语
本文通过简单的计算实例重新探讨了差分格式对非线性计算稳定性的影响。这里的“简单”，主要指将一维非线性平流方程的时间偏导项统一地取成前差格式，同时，差分方程仅写在三个格点上。从而，隐式格式差分方程的求解过程便成为三元一次方程组的求解过程。该求解过程比显式格式差分方程复杂，但计算结果稳定，充分体现出隐式格式和显式格式的优缺点。虽最终的计算结果与教材中[2]基本一致，但本文举例更为简单、易懂，且给出了详细的求解过程，有助于学生自己动手推导和计算求解，以加深其对显式和隐式格式的理解，并深刻体会各自的优缺点。
参考文献：
关键词：激光 相变硬化 温度场 模拟计算
中图分类号：TG665 文献标识码：A 文章编号：1672-3791（2013）07（b）-0080-03
激光是20世纪60年代产生的重大科学技术成果之一，随着现代化工业的发展，对产品的热物理性质提出了各种各样的、新式的和特殊的要求。激光加工对传统工业的改造发挥着很重大的作用，目前在很多方面，如汽车、电子、电器、航空、冶金、机械制造等重要部门占据了重要的地位。在新世纪，激光技术已经形成了自身的产业，在现代信息社会光电技术产业贡献颇大，激光产业已经是受到科学界高度重视的骨干产业。
凭借激光自身的特殊性质，加工精度高、适应性强、加工效率高、灵活性高等优良性质，能够进行激光相变硬化、激光焊接、激光合金化、激光切割、激光打孔、激光快速加工以及对材料的微型加工等作用。
激光相变硬化技术是激光加工技术中，起步比较早，另外发展比较成熟的一门高科技新型技术。目前国内外都有相关的研究，并取得了一批重要的成果，有些也已经应用于各种机械零件的表面硬化。但是在这些研究成果和应用中仍然可以发现很多不足之处，所以对激光相变硬化的进一步研究是十分重要的。
1 研究方案及意义
本课题主要对激光扫描加热准稳定温度场数值模拟计算设计进行研究。对于确定的材料，怎样选定合适的激光加工工艺参数，因为激光加热过的金属表层区域的金相组织、物理、化学性能会发生变化，有表面局部硬化的表现。其结果一般都采用试验或凭经验加以确定，这样就要耗费大量的财力和物力，也给激光技术的更广泛应用带来不便。采用计算机数据仿真模型，模拟激光扫描加热操作，得到相关数据，从而降低实验的研究成本，带来更大的市场商机。
我们通过这次的研究课题，制作出一个简单明了的，具有友好型界面的软件，可以计算出我们想要的一些相关数据。激光扫描加热处理材料，在材料中的每个点的温度值是不同的。首先我们所要讨论的是激光加工工艺参数与材料参数之间的关系，建立数据仿真模型，计算出合适的值。确定激光扫描加热区域温度场的数值仿真模型，然后对该模型进行求解，得到激光加热工程的温度场分布，之后对激光加热辅助切割的主要影响参数做进一步的仿真及相关实验研究，分析这些激光加工工艺参数对加工质量的影响，优化加工用的激光参数，如光斑中心离刀尖距离、激光的功率密度、工件的旋转速度、工件表面光斑直径、激光工作的温度场等。在我们对与激光加工工艺参数进行优化的同时，对于数据的计算也进行一定的优化，使计算公式能够更加的清楚。
2 数学模型的建立
2.1 问题的简化
在激光扫描加热过程中，由于问题往往比较复杂，材料各点的温度多变，难以进行分析，所以在求解过程中，我们做如下的设定，根据假设条件来简化方程。
（1）材料的热物理性质不随温度而变化。
（2）除了激光束开始辐射和即将结束两个阶段外，相对于运动坐标系的温度分布为准稳定态。
（3）相变潜热相对于激光束热流相当小，可以忽略不计，及H=0。
（4）在y=0，=0，符合绝热边界条件，及y方向的温度场分布以y=0面对称。
（5）激光束具有均匀的能力分布，q=Q/2a×2b，式中Q为激光束的输出功率，q为热流，2a，2b，为激光束的尺寸。
（6）材料对激光的吸收系数不随温度变化而变化，但与扫描速度有关，以F（V）来加以修正。
对问题简化了以后，为之后公式的推导提供了前提。激光对物体加工轨迹如图1所示。
2.2 导热偏微分方程
在三维直角坐标系中，导热偏微分方程为：
其中Kii（i=x，y，z）为x，y，z三个方向的导热系数，H为固态相变潜热，P为材料密度，C为材料比热。
对于各向同性材料，导热系数为常数，令×则上式可变成：
对准稳定，则变成：
2.3 激光热处理瞬态温度场快速算法
通常，激光热处理是一个短暂的过程，热影响局限于材料的表层，同时，由于材料在激光扫描后基体对表层热能的迅速扩散是使表层完成淬火过程的原因，所以只要研究光束在临近区域的温度变化，就可以有效的预计热处理的结果。因此在光束临近区域并与光束共同运动的动坐标中对问题进行讨论。
在半无限大均匀的介质上建立三维直角坐标系，坐标指向材料的内部，如果材料表面存在速度v沿着x方向运动的单位强度点热源，并在与光源共同运动的动坐标中温度分布与时间无关以及光源处于动坐标原点时，温度场计算公式是：
式中，To为光源与材料相互作用前材料基体的温度，在观察时刻t运动坐标系与固定的坐标系的关系为：
按照温度场的线性叠加性质，功率密度为的面热源在半无限大介质内激起的温度场就可以表示为与（1）式右端的卷积为：
沿用以上的推导公式，在t=0时开始在原点处引入一单位强度点热源，此后热源以速度v沿着x轴正向运动，在介质的热物性参数与温度无关，表面满足绝热边界条件的情况下，时刻t介质的温度长Ti（x，y，z）可以写为：
（上式中，为介质的热扩散系数；k为介质的导热系数；t为热作用时间；To为t=0时刻介质的初始温度分布（K）；是函数。
关键词:数值分析;Matlab;教学实践
Abstract:ThispaperpointsoutthefeaturesofnumericalanalysiscourseanditspresentshortcomingsandintroducesMatlab,anadvancedsoftwareofnumericalanalysis.Besides,withregardtotheteachingpractice,thepaperstudiestheteachingreformofMatlab-platform-basednumericalanalysisforpostgraduatesofengineeringmajorfromthefouraspectsofvisualizingtheabstract,simplifyingthecomplexcomputing,enhancingthebaseofnumericalexperiment,andimprovingabilityofcaseanalysis.
Keywords:numericalanalysis;Matlab;teachingpractice
随着计算机技术的迅猛发展,科学计算已经与实验研究、理论分析并称为科学研究的三大方法,借助于计算机和数值计算理论,人们能对自然科学、工程技术、经济管理及至人文社会科学领域中的数值模型求出数值解。因此,要求改革数学课程教学的呼声日益增长,而数值分析课程以其独特的特点首当其中,提高数值分析教学质量,培养学生数学素养,掌握实用算法并能熟练利用计算机求解成为当务之急。
一、数值分析课程的特点与现状
(一)课程内容多杂而教学时数少
本课程包括了数值逼近(插值、函数逼近,数值积分与微分),数值代数(线性方程组求解、矩阵特征值与特征向量),方程求根(非线性方程(组)求解,解微分方程)等内容。在运用传统教学方法讲授这些知识时,由于公式多,推导过程繁琐,加上教学时数少(54学时),很多内容难以说清楚,致使学生产生厌学情绪,教学效果大打折扣。
(二)计算复杂
解题时,一般都要进行大量的计算,不是一支笔、一张纸,外加一台计算器能顺利解决的,因此学生只能做一些简单的数值模拟题,难以领会和理解方法的计算要领和步骤,体会问题的条件和限制范围,理解一般问题和特殊问题的区别。
(三)重理论轻实践
传统课程只注重讲授数值方法的原理,课堂教学占去整个教学过程的绝大部分时间,学生没有实验、实践的深刻体会,不能全面理解和运用教材中的算法。
(四)直观性差
课堂教学中难以对一些重要概念、重要现象进行直观展示,学生只是被动记忆一些结论,并不真正理解。
二、Matlab——优秀的数值计算软件
美国Mathwork公司于1967年推出了适用于不同规格计算机和各种操作系统的数学软件包-Matlab[1],它集数值计算、符号分析、图形可视化、文字处理于一体,语法简单、操作方便、界面友好,只要有点Windows操作经验,在短时间内就可学会它的操作和使用方法,而且其编程效率远优于BASIC、FORTRAN、C等软件。它具有很好的开放性,以它为基础开发的二十多个工具箱,可用于解决诸多学科专业中的数值计算问题。该软件已成为发达国家高等院校理工科学生必须掌握的基本软件,也是科学研究和工程设计部门解决具体问题的一种标准软件。用过（中国整理）它的人都感叹:用Matlab处理矩阵-容易;用Matlab实现可视化-轻松;用Matlab编程-简洁!
三、数值分析课程教学改革实践
针对数值分析课程的特点和传统教学中的不足,笔者在多年的教学实践中,将Matlab软件应用于数值分析教学,取得了良好的教学效果。主要做法如下。
(一)抽象内容直观化
心理学研究表明,与抽象的内容相比,学生易于识记生动、形象、有趣的学习材料。如果能将抽象的数学知识直观的呈现在学生们面前,无疑将会极大地激起学生的学习兴趣,Matlab强大的可视化功能正好能做到这一点。一个典型的例子是在引入分段低次插值时,为了让学生更好地理解Runge现象,利用屏幕动态地显示f(x)=11+x2插值函数图象[2]31。
在[-5,5]上取等距节点(给定n,共取n+1个点)
xk=-5+k·10n,k=0,1,2,…,n。
构造拉格朗日插值多项式Ln(x),随着n的增大,在区间端点附近Ln(x)与f(x)接近程度越差,形象的说明了高次插值函数近似f(x)的效果并不好,为讲授分段低次插值做了个很好的引子。
借助于Matlab平台,函数逼近中的误差分布,数值积分方法的改进,迭代过程等等均能很容易直观地呈现在学生面前。
(二)复杂计算简单化
数值分析难,主要难在运算过程(公式)复杂,大部分时候难以靠一支粉笔在黑板上一步一步的演算,学生课后练习也只能做一些简单的数值模拟。而Matlab强大的数值计算功能,帮助我们解决了这个问题。
例如,用雅可比迭代、高斯-赛德尔迭代和超松弛迭代求解线性方程组的数值计算,计算量非常大,按照教材讲解,只能简单地告知结果,学生积极性不高,利用Matlab编写三个简单的小程序,便可将计算结果(甚至每一步的迭代结果)直观地展示给学生。
例1分别用雅可比迭代、高斯-赛德尔迭代和超松弛迭代解线性方程组[2]213
并比较收敛速度,取精度=10-5。其精确解为x*=(-1,-1,-1,-1)T。
教学中,用自编的程序,很快得出了下列结果,并比较了方法的优劣(见表1和表2)。
数值分析教学中,这样的例子比比皆是,只要教师认真设计,不仅能大大激起学生对学习内容及过程强烈的兴趣,而且还对初步培养学生科学计算能力起到了重要作用。
(三)数值实验强基础
数值分析是一门实验性较强的学科,上机实验不仅能加深学生对算法稳定性,理论可靠性及计算复杂性的理解,培养学生的编程能力,还能培养学生质疑问题的能力和创新精神。因此每章结束后,都提供一两个问题要求学生利用MATLAB软件,自编程序或利用其库函数求解问题,分析结果。如学习线性方程组迭代法后,提供一个高阶病态的线性方程组,要求学生用各种方法上机求解,并对结果进行分析,找出收敛较快的迭代法,寻求最佳松弛因子。
(四)案例分析长才干
实践性是数值分析课程区别于其它数学课程的一个重要特征。每个部分内容结束后,安排一个案例分析,帮助学生从“算”数学过渡到“用”数学。如海底测量(插值),估计水塔的水流量(插值、数值微分与积分),投入产出分析(线性方程组),商品的产量与价格(方程求根),导弹系统的改进(微分方程)等。通过这些案例分析,既使学生认识到数值分析的实用性,又让学生领略了Matlab的强大功能,积累了用数学软件解决实际问题的经验。由于教学时数的限制,案例分析一般由学生课后完成,教师负责指导。
四、结束语
数值分析是研究如何用计算机解决实际问题的课程,将Matlab与数值分析课程结合起来,开阔了学生思路,拓展了解决问题的方法,取得了较好的教学效果,学生做毕业论文,甚至在实际工作中遇到有关难题时,经常当面或通过邮件与我进行探讨。学生普遍反映数值分析课程内容多,实践性强,应用广泛,费时费力,但学下来最有用,收益终生。
参考文献:
[1]]周品,赵新芬.Matlab数学建模与仿真[M].北京:国防出版社,2009.
[2]李庆杨,王能超,易大义.数值分析(第四版)[M].武汉:华中科技大学出版社,2006.
[3]刘彩云,.《数值分析》课程教学改革的若干思考[J].长江大学学报:自然科学版,2009,6(1):358-359.
