数值方法精选(九篇)

前言：一篇好文章的诞生，需要你不断地搜集资料、整理思路，本站小编为你收集了丰富的数值方法主题范文，仅供参考，欢迎阅读并收藏。

第1篇：数值方法范文

关键词： Matlab 数值方法； 指纹； 曲线拟合； 插值

中图分类号： TN911?34 文献标识码： A 文章编号： 1004?373X（2015）24?0027?04

Comparison and optimization of numerical fitting methods for fingerprint curves

LI Qing， ZHANG Wenjie， CHEN Jing

（School of Science， Beijing Forestry University， Beijing， 100083， China）

Abstract： Fingerprint identification， as a popular identification method， has been widely used in data encryption. The different curve fitting methods such as the least?squares fitting， polynomial curve fitting， monomial interpolation， Lagrange interpolation， Newton interpolation， polynomial interpolation in subsection， spline interpolation and so on are utilized in this paper to make the fingerprint curve fitting by means of Matlab software and on the basis of numerical methods. And then the fitting results and algorithm are compared to select an optimal method. The generalized extension approximation method is used in fingerprint curve fitting to deal with the bad fingerprint curves to further improve the fingerprint curve fitting. The summary， comparison and expansion of various fingerprint curve fitting methods were achieved.

Keywords： Matlab numerical method； fingerprint； curve fitting； interpolation

0 引 言

指纹识别作为近年来比较流行的识别方法，在数据加密、交通、门禁、司法等各个领域得到了越来越广泛的应用。在进行指纹曲线拟合时，可能会遇到指纹曲线类型多，拟合难度大，拟合效果不好的情况[1]。当前市面上的指纹识别，在断裂指纹的拼接这一方面，基本上是采用Gabor算子滤波的方式，虽然使用Gabor算子滤波效果不错，但是这种滤波要牵涉到一系列的卷积运算，通常是对整个指纹图像滤波，程序开销比较大，滤波时间也较长。然而，另一方面，如果用数值方法进行断裂曲线的拼接，它只是提取那些断裂的指纹曲线，用数值方法进行拟合，其他没有断裂的曲线，则保持原样。这样做的好处是，由于指纹中断裂的曲线一般不多，所以拟合比较快捷，程序开销较小。指纹断裂曲线的拼接属于指纹图像预处理环节，当把断裂指纹拼接上，之后的指纹识别就方便了。

曲线拟合是用连续曲线近似地刻画或比拟平面上离散点组所表示的坐标之间的函数关系的一种数据处理方法，其应用性非常强。本文基于数值方法，借助Matlab软件，利用最小二乘拟合、多项式曲线拟合、单项式基本插值、拉格朗日基本插值、牛顿基本插值、分段多项式插值和样条插值等不同曲线拟合方法对指纹曲线分别进行拟合，并对拟合效果和算法进行比较。在此基础上，针对劣性情况的指纹曲线，将广义延拓逼近法应用于指纹曲线拟合中，进一步完善指纹曲线拟合，实现了各种拟合方法的总结和比较以及扩展。

1 指纹数据的最小二乘拟合[2]

1.1 正规方程解指纹曲线的最小二乘拟合问题[3]

图1表明，残差为0.992 0时拟合质量良好，图2显示的是一种劣性情况，拟合质量比较差。

1.2 用QR分解法求最小二乘拟合

除了正规方程，最小二乘问题也可以由QR分解来解决。通过一些实现最小二乘法的Matlab内置命令，使用QR分解法来解超定方程组。正规方程法和QR分解法在解决最小二乘问题的差别不大，主要不同点在于计算所花费的浮点操作次数和数值的稳定性[4]。

图1 最小二乘拟合良性情况

图2 最小二乘拟合劣性情况

由图3，图4可知，对于良性问题，解正规方程和用QR分解法结果几乎一样，但对于劣性情况，QR分解法在使得残差最小这方面要优于解正规方程法[5]。

图3 QR分解拟合良性情况

2 指纹曲线的多项式拟合

由图5，图6可知，对于指纹曲线的劣性情况，二次多项式和三次多项式的拟合效果均不好，高次多项式拟合效果比低次拟合效果稍好，但并不等于次数越高效果越好。而简单的多项式拟合可以反映出曲线的大致弯曲程度，但对于劣性情况拟合的效果并不好。

图4 QR分解拟合劣性情况

图5 多项式拟合良性情况

图6 多项式拟合劣性情况

3 基本插值公式

3.1 用单项式基本插值公式进行插值拟合[6]

图7结果可见，警告信息指出了程序运行结果A=1×1016，A是病态的，RCOND = 3.170 723×1032是估计条件数的倒数，它的值很小，说明矩阵的条件数非常大。如图7所示，此时的插值函数出现严重的问题[7]。插值函数的曲线不是预料中的平滑曲线，而是不规则曲线。矩阵A的条件数非常大，因为其中的元素值的大小差异很大，它们中的最小元素是0，最大的元素是1×1016×4.045 4。在求解过程中的消去阶段，很小的舍入误差会导致很大的不确定性。在后面的代入过程中，会有更多的舍入误差，并对Ci的真值产生扰动，引起图中所示的振荡。舍入误差的问题是使用单项式基本插值公式进行多项式插值时固有的，这时因为Vandermonde矩阵通常是病态的。由于多项式求解中的舍入误差，插值函数对输入数据的逼近会变得不精确。

图7 单项式插值拟合

3.2 用拉格朗日基本插值公式对指纹离散坐标进行

插值拟合

如图8所示，没有类似于图7的警告信息，插值函数中也不会出现类似于图7的乱真振荡，这是因为拉格朗日多项式不需要线性方程组的解，也就不会有破坏Vandermonde方程组解的精度的病态性和舍入误差问题。

图8 拉格朗日插值拟合

3.3 离散点的牛顿插值公式对指纹离散坐标进行

插值拟合

如图9所示，牛顿形式插值采用均差形式比较方便，n阶插值多项式可通过向n-1阶牛顿多项式添加一个更高次项获得。它计算比较高效，并且随着多项式阶数的增加，它仍然能够保持比较可行的操作次数。牛顿形式插值的结果取决于支点的选取和插值多项式的阶数。

3.4 分段多项式插值对指纹离散坐标进行插值拟合

如图10所示，斜率不连续性表现得比较明显。在进行分段线性插值的实现中，比较重要的是在构造插值式时选用合适的支点对。在一个对任意数据集进行分段线性插值的Matlab函数中，支点的查找和插值函数计算都可以用Matlab程序自动完成。

3.5 折半搜索法进行分段线性插值对指纹离散坐标进行插值拟合

如图11，图12所示，此算法显示出一个弊端就是在节点处虽然是连续的，但曲线并不光滑。节点处导数的连续性没有受到约束，所以还有待于进一步改进。

图9 牛顿插值拟合

图10 分段式插值拟合

图11 折半搜索法分段线性插值（一）

图12 折半搜索法分段线性插值（二）

3.6 三阶样条插值

三阶倦条插值结果如图13所示，离散指纹数据对整个固定端点斜率进行三阶样条插值，插值的效果并不好。在中间空白区域，插值后得到的曲线受力断裂部分最近的两个端点的影响较大。这两个端点的斜率的大小会影响到断裂部分插值的效果。但整个数据对集的边界端点的斜率大小，对于断裂部分插值的效果影响不大。当整个数据对集的边界端点为不同值时，整个曲线的变化并不明显，需要进一步改进。

图13 三个阶样条插值拟合

4 广义延拓逼近法Matlab实现

首先将指纹曲线两个端点处的片段拟合起来，因为端点处的片段不方便进行单元域的延拓。程序中采用了Matlab内置的interp1来实现，用三阶样条同样可以实现，效果差别不大。除两端点外的其他片段都采用广义延拓逼近法实现。将每个片段的拟合函数的系数矩阵求出来，即可求得每个片段的逼近函数，再拼接起来，即可构造相应的指纹曲线[8]。由图14和图15可见，拟合的结果比较理想。可以通过延拓域的节点来约束当前段的拟合，可操作性较好，效率较高[9]。

5 对比分析和结论

最小二乘拟合能很好的跟踪离散数据点趋势。求单项式基本插值公式的系数需要解Vandermonde矩阵，而矩阵可能是病态的，这样会导致单项式系数不确定。另外单项式中的各项可能在大小上有很大差异，会导致多项式计算中的舍入误差。这些数值上的复杂性可在构造Vandermonde方程组之前通过对数据x进行转换和缩放来改进。拉格朗日基本插值公式插值在理论分析中很有用，由拉格朗日基本插值公式形成的多项式插值式的舍入误差比较小，但是计算比较繁琐。牛顿多项式基本插值法不但舍入误差比较小，而且计算公式也很高效.换句话说，对数值计算来讲，使用牛顿基本插值公式进行插值比使用单项式或拉格朗日基本插值公式更好。牛顿插值多项式的系数是输入数据集的均差，均差在推导三阶样条时很有用。需要注意的是，对均匀分布的支点数据进行任意阶的多项式插值时，必须防止由多项式摆动带来的误差。多项式摆动会影响任何基本插值公式上定义的多项式。分段多项式插值使用相对低阶的多项式，它们定义在输入数据的子集上。对一维数据来说，插值式是由很多插值式在某些点上连接而成，所以需要选择好合适的局部插值式。若采用三阶样条插值，插值式的总体平滑比较好，样条插值式的精确度要受端点条件的影响很大，但通过本文的实验，三阶样条在指纹曲线的拟合方面，效果并不好[5]。

图14 广义延拓逼近法（一）

图15 广义延拓逼近法（二）

对于劣性情况，本文将广义延拓逼近算法引进到指纹曲线的拟合。所谓广义延拓逼近算法是通过在延拓域上进行拟合逼近，在边界上进行插值约束处理，将插值法和拟合法有机地结合，从而形成了广义延拓算法自己的特点。广义延拓外推模型基本上继承了最小二乘拟合的长处，可以充分地运用较多的实验数据点，同时可以把最新的数据点锁住，充分发挥最新数据的作用与影响，也就是用最新的数据进行插值约束处理。通过本文的实验，广义延拓逼近算法在指纹曲线拟合方面的效果比较好。

6 结 语

本文利用Matlab工具，通过数值方法，对指纹曲线使用多种方法分别进行拟合，选择将广义延拓逼近法应用其中，进一步完善了指纹曲线拟合，实现了各种拟合方法的总结和比较以及扩展。通过数值方法对指纹曲线进行拟合，拟合效果较好，但其缺点在于需要对指纹曲线一根一根拟合，所需工作时间可能要较长。但是，这也是它的优点所在，正是由于它是一根一根的拟合，对于没有断裂的曲线，则保持原样，不需要所有曲线（包含未断裂的）都进行拟合，减少了工作量。本文所提到的这些拟合插值方法，不仅可使用在指纹曲线的拟合中，在其他需要曲线拟合的研究领域也有广泛的应用。

参考文献

[1] 李昊，傅曦.指纹模式识别系统算法及实现[M].北京：人民邮电出版社，2008.

[2] 解同信.最小二乘法求作拟合直线[J].北京工业职业技术学院学报，2006（3）：5?7.

[3]吕喜明，李明远.最小二乘曲线拟合的Matlab实现[J].内蒙古民族大学学报：自然科学版，2009（2）：125?127.

[4] 刘霞，王运锋.基于最小二乘法的自动分段多项式曲线拟合方法研究[J].科学技术与工程，2014（3）：55?58.

[5] 李庆扬，王能超，易大义.数值分析[M].北京：清华大学出版社，2001.

[6] 陈玉洁.多元多项式插值[D].长沙：国防科学技术大学，2003.

[7] 韩艳丽，李凤萍.插值与曲线拟合[J].中国西部科技，2009（11）：91?92.

第2篇：数值方法范文

关键词： 函数 定义域 值域 值域的求解方法

设 是非空的数集，如果按照某种确定的对应关系 ，使对于集合 中的任意一个数 ，在集合 中都有唯一确定的数 和它对应，那么就称 为从集合 到集合 的一个函数，记作 ，其中 叫做自变量。 的取值范围 叫做函数的定义域；与 的值相对应的 的值叫做函数值，函数值的集合 叫做函数的值域

由函数的定义可知，一个函数的构成要素为：定义域、对应关系和值域。其中函数的值域是一个较复杂的问题，又是高中数学中的一个难点。总体来讲，求函数的至于要注意以下几点：（1）值域的概念，即与 的值相对应的函数值的集合 ；（2）函数的定义域。当题目中未明确给出函数的定义域时，应先求出函数的定义域，在定义域的范围内求函数的值域；（3）函数的单调性。求函数的值域时，常常借助函数的最值来求解，而求函数的最值时，对函数的单调性的讨论往往是必不可少的；（4）函数的解析式。在求函数的值域时，往往要根据所给函数的解析式的形式，使用相应的方法。具体常用的求函数值域的方法如下：

（1）观察法

对于一些简单的常见的函数，通过观察就可以求出其值域。例如我们熟悉的一次函数的定义域是 ，值域也是 ；反比例函数 的定义域为 ，值域为 。

（2）配方法（或公式法）

（3）换元法

（4）分离常数法

（5）利用函数的单调性求值域

例5. 求函数 的值域

解：由题可知函数的定义域为 ，因为 和 在 上均为增函数，故原函数为 上的增函数.所以 ，故原函数的值域为

（6）利用函数的最值求值域

对于区间上的连续函数，利用求函数最大值和最小值来求函数的值域。

总之，同学们在学习的过程中应多注意积累，善于总结，从而在求解函数值域的问题中，才能迅速找到求解此类问题的比较简单且合适的方法。

第3篇：数值方法范文

【关键词】离心泵性能曲线；数据处理；正交拟合

离心泵特性曲线一方面用于设计制造中检测离心泵性能是否达到设计要求，同时也作为正确选择、使用离心泵的主要依据。离心泵的特性测试是为了获取离心泵的性能曲线。了解并掌握水泵性能对我院轮机、热动、建筑环境专业的学生来说是非常重要的，性能参数对于水泵设计、实验和使用都是十分重要的基本参数，又是学生不易理解的参数。我院动力实验室考虑实验经费问题，在原有基础之上，建立闭式系统水泵综合性能实验台，并且主要从实验手段上拓展学生的知识面和增强了动手能力，激发学生的学习积极性。

1、实验测试系统及原理

试验台装置示意如图1所示。

① 流量Q（m3/h）由流量计直接测量。

② 扬程H（mH2O）进口断面和出水断面之间的总水压头，即这两个过水断面之间的静压差和动压差。

③ 水泵轴功率N（kW）泵和电机采用联轴器传动，故电机输出功率为泵的轴功率（1）

其中N电―电机输入功率，由功率表测出kW；η电―电机效率；η传―联轴器传动效率，可取0.98。

④ 泵的效率η泵的有效功率与轴功率之比：

2、正交多项式数据处理方法

在绘制性能曲线时，以往学生只是将测试的数据在坐标纸上进行描点连线，单调枯燥，往往实验做完，不知道自己做了什么；而且手工方式无法达到精度要求，计算机编程处理虽然能满足精度要求，但需繁杂的编程过程，对选型设计人员要求较高。为了适应选型的要求，需要找到一个既快速又简便的方法，求出离心泵性能曲线的函数表达式，从而将离心泵性能曲线转换为计算机图形。

以H―Q曲线为例，提出利用正交多项式法对离心泵的性能曲线进行快速拟合。通过直观观察离心泵实际H―Q曲线的特点以及对大量数据的拟合和比较，认为采用多项式拟合法来拟合离心泵 H―Q曲线为较好。设拟合曲线为m次多项式，可将数学模型形式表示为：

式中Hi、Qi分别为扬程（m）、流量（m3/s）；bj为拟合系数（j=0，1，2，…，m）

3、实验计算实例及分析

通过实验，并经过公式测算，得到表1所示数据。对H―Q曲线进行拟合如图2，并得到其相关系数R2，同样的方法对其它曲线进行拟合，得到方程及相关系数如下：

为检验其可靠性，本人对三条曲线分别进行了正交多项式拟合得到如表2结果，可见拟合可以保证精度要求。

4、结论

在数值处理的方法上，还有最小二乘法，插值法等，其精度和实用性各有差别，但正交多项式法测试表明：1）由正交多项式法得到的拟合曲线在工作区内与实测性能曲线基本重合，误差较小，因此，可以代替手工对离心泵性能数据进行处理并作性能曲线图。拓展学生的知识面和动手能力，激发学生学习兴趣。2）从返回的R2来看，置信度在0.99以上，说明该方法可信、可靠。3）该方法处理过程快速直观，无须进行烦琐的编程。4）本方法便于应用计算机进行离心泵的选型和计算，可大大提高工作效率其拟合精度均较高，因此，本方法具有广泛的实用价值。

参考文献

[1]孙杰.自动回归离心泵性能曲线的正交多项式法[J].水泵技术，1991（3），60～64

[2]郁飞.试验设计与数据处理[M].北京：中国标准出版社，1998：38～49.

[3]李丽.水泵综合性能实验台及其应用[J].实验技术与管理，2005（4），44～47

第4篇：数值方法范文

关键词：数值计算方法；创新意识；计算平台

作者简介：张俊丽（1980-），女，山东菏泽人，内蒙古民族大学数学学院，讲师。（内蒙古?通辽?028000）

中图分类号：G642.0?????文献标识码：A?????文章编号：1007-0079（2012）28-0087-01

随着科技的飞速发展和计算机技术的广泛应用，数值计算方法已成为重要的桥梁和工具深入到航天航空、地质勘探、汽车制造、桥梁设计、天气预报等各个领域，成为每一位科研人员和工程技术人员所必备的知识。为了满足社会需求，数值计算方法现已成为高等院校理工类学生的一门专业必修课程，其目的是让学生掌握设计数值算法的基本方法，培养学生分析问题与解决问题的能力，为以后用计算机解决科学计算问题打下坚实的基础。

一、“数值计算方法”课程的特点与教学现状

数值计算方法，简称计算方法，又叫数值分析，是一门研究数学问题的近似解并利用计算机进行数值实现的学科，是数学分析、高等数学、高等代数、概率统计等数学基础课的后续课程，它既有数学理论上的抽象性与严谨性，又有实验性与应用性的数值特征。计算方法课程的内容包括插值和拟合、数值微分和数值积分、求解线性方程组的数值方法（直接法和迭代法）、非线性方程数值解、矩阵特征值计算及常微分方程初值问题数值解法等；[2]它的计算对象是数学中的微积分、线性代数、常微分方程，只是它不像别的数学课程那样只是研究纯粹的数学理论，而是把数学理论与计算相结合，重点探讨数学问题的数值解法及应用；它的课程要求是在掌握算法原理的前提下设计算法编程实现。

“数值计算方法”是一门介绍科学计算的核心理论和基本方法的数学课程，它对培养学生的科学计算能力和解决实际问题的能力具有不可替代的作用。从20世纪80年代起，“数值计算方法”相继成为各高等院校数学及其他理工科（如物理、计算机等）专业本科生的一门专业基础课。但内蒙古民族大学（以下简称“我校”）的数值计算方法课程只在应用数学、信息与计算科学两个专业开设必修课，一般开设在第三或第四学期，理论课48学时，上机实验16学时，在别的学院（如物理、计算机等）没有开设该课程。该课程普遍存在的教学现状是：理论课内容多，学时少，各部分内容不连贯，公式繁多，枯燥乏味，使得学生产生厌学情绪；上机课时间紧，且一般集中上机，与理论课内容脱节，失去了上机实验操作的意义；很多时候这门课程的学习都结束了，学生还不清楚这门课程与原来的课程有什么联系，学习这门课有什么用，更无从谈起培养学生的创新能力；而且“数值计算方法”课程教学过程中还存在着教学内容陈旧、教学方式落后及考试形式单一等问题。针对该课程目前的教学现状，如何对该课程教学进行教学改革，是值得深入思考的问题。

二、关于“数值计算方法”课程改革的若干建议

根据前文分析可知，目前“数值计算方法”课程教学中存在着一些不容忽视的问题。那么如何进行教学改革，培养学生的实际应用能力，体现该课程在工程科学中的价值和意义，是值得数学界思考的问题。根据近年来我校师生在该课程教学中出现的问题，本文对“数值计算方法”课程教学改革提出以下几点建议：

1.优化教学内容，选择合适教材

“数值计算方法”课程讲授时既要强调它的理论结构与使用价值，又要注重提升它与计算机使用密切结合的实用性特点，所以该门课程对教材的要求很高。然而现行教材有的理论偏深，不适合普通本科生使用；有的内容陈旧，与实际联系缺乏；有的实用性强，但与实践结合的算例较少；[3]再加上该课程内容抽象，知识连贯性不强，定理和公式较多，推导过程烦琐，从而导致学生对该课程的学习没有兴趣，只是为了应付考试机械性地记忆公式。按照教育部关于“数值计算方法”课程在教学过程中应把握“重概念、重方法、重应用、重能力”的培养要求，对该课程的教学内容应灵活把握，知识点讲解应详略得当，不同专业的学生对该课程的要求不同，讲解的侧重点也应有所不同，最好选用的教材也不同。对数学类的学生来说，理论与实践应并重，而对于别的理工科的学生来说，不在于理论的论证与推导，而应侧重算法原理与实际应用。当选定教材后，在实际教学过程中还需要对教学内容灵活整合，对于一些复杂且后继课程将会深入学习的内容（例如微分方程的数值解法等），[4]可以略讲甚至不讲。不同地区的高校对该课程的教学要求也略有不同，例如我校处少数民族地区，学生的基础知识相对较差，在该课程授课时更应减少烦琐公式的推导，重在加强学生对知识点的掌握与实际应用能力的培养。鉴于该课程对以后学习和工作的重要性，我校建议除了数学与信息类的学生以外，别的理工科（如物理，计算机、信息工程等）的学生也应开设数值计算方法课程的选修课。我院本专业教师在包玉兰教授的带领下，根据我校学生的状况及多年积累的教学经验，编写了比较适合少数民族地区学生特点的数值计算方法教材，现已经出版在我校试用。该教材内容较浅，并配备一定量的习题和上机实验题，要求理论学时50～60学时（包含习题课），上机实验16～20学时，并且标注了一些选讲的内容，不同专业的学生可以针对性地学习，[5]基本上满足了我校学生对该课程教材的要求。

第5篇：数值方法范文

数学是科学之母。一门学科之是否成为科学，决定于该学科的问题描述是否能化归为数学。工程技术属于应用科学范畴，工程技术问题通过建立数学模型与数学产生直接联系。数学问题的分析解通常是极难得到的，因此必须归结为数值计算问题。例如：人造飞船的轨道研究、汽车耐撞性问题研究、大型桥梁设计、天气预报等都必须数值求解。数值计算方法作为研究数学问题的近似求解方法的课程，既有一般类数学课程理论上的抽象性和严谨性，又有工科类课程的实用性和实验性特征，是一门理论性和实践性都很强的学科。该课程理论涉及面广、方法应用性强、内容丰富，再加上随着计算机技术的飞速发展，优秀数学软件层出不穷，数值计算方法更能与计算机相结合，适应科学发展的需要，现已成为各高校大部分理工科专业的必修课程。在数值计算方法的教学过程中，笔者发现了很多问题。本文对其中的部分问题进行了分析，并提出了几点教学改革建议。

二、教学过程中存在的问题

以笔者所在的机械工程专业为例，起初该课程为学科选修课，选课学生少，且其中大部分是为了凑学分而来的，学习兴趣不高在所难免。后来学科培养计划改变，该课程归入专业必修课，选课学生数量增加了，但是学习热情还是不高。究其原因，主要有以下几点：

1.课程对数学基础要求较高。本课程主要解决以下几大类问题：非线性方程求根、线性代数方程组求解、矩阵特征值与特征向量的数值解法、插值与拟合、函数最佳逼近、数值微分与积分、常微分方程初值问题的求解等。需要先修的数学课程包括高等数学、线性代数等。学生只有掌握这些课程中的基本内容，才能学好数值计算方法课程。而这几门课程均是难度较大的数学课程，学生的掌握程度本来就不好，甚至学过后已经忘记。由于同时要学习其他机械专业课程，学生不愿再花大量的时间和精力去学习或复习相关的数学知识，特别是本来就对数学不感兴趣的学生。所以在课程学习中，学生就会陷入“听不懂，听不懂就没兴趣，没兴趣就不想听课，不听课就不懂”这样一个死循环。

2.教学课时的限制。该课程的内容覆盖很广，如“插值方法”这一章，就算法而言就有Lagrange插值、Aitken插值、Nevile插值、差分与差商形式Newton插值、Hermite插值、分段低次插值、三次样条插值、B-样条插值等。然而，总学时设置仅为32学时。即使不面面俱到，挑选几种典型的插值方法讲解，也需要花费不少时间。因此，教师的课堂时间主要用来讲解各问题相关算法的理论推导和算法设计，几乎没有帮学生回顾相关数学知识的时间，且在课堂内也没有时间及时将理论运用于具体问题。学生觉得这是一门纯数学课，枯燥无味又难懂，没有学习兴趣。

3.没有与计算机很好结合。数值计算方法的一大特点是面向计算机。一个好的数值算法要通过程序设计在计算机上实现，要求用最简练的语言、最快的速度、最少的存储空间来实现某种计算结果。要达到上述要求，要求教师和学生既要掌握数值计算算法，又要能熟悉并熟练使用计算机语言。而现在的课堂教学重点大都侧重在理论讲解上，没有很好地结合计算机编程，学生把这门课当成了数学课来上；且学生在课外也没有将课堂上学到的算法付诸于计算机上实现。导致该门课程理论与实践严重脱节，达不到预期的教学效果和教学目的。

三、如何提高课程的趣味性

综合上述教学中出现的问题，要想教好这门课、使学生学到知识，最重要的是要提高课程本身的趣味性。“兴趣是最好的老师”，学生有了兴趣，才会有学习的热情，才会把精力付诸于课程学习上。那么关键问题是如何提高该课程的趣味性，主要从下面几点出发。

1.结合专业特点，从实际出发，合理安排课时和教学内容。由于课时有限，而且授课对象主要是机械这样的工科类专业的本科生，课程的主要目的是培养学生具有机械工程工作所需的数学能力，培养学生用数学的思想方法分析问题、解决问题的意识和能力，并为后续的工作和学习打下良好基础。那么教师在安排课时要懂得取舍，选择与机械专业紧密相关的内容讲解，课程主要浓缩为如下几个主要内容：非线性方程的求解、线性方程组的求解、插值与拟合、数值微分与积分、常微分方程数值求解。而每个内容应该针对其中的经典算法进行讲解。如非线性方程的求解，只需就二分法、简单迭代法、Newton迭代法进行详细讲解，其他算法如弦割法、简单Newton法等只需简单提及即可；常微分方程的数值解法，只需对Euler法和Runge-Kutta方法详细讲解，其他内容略讲即可。例如非线性方程求解中，判断迭代法收敛的充分条件，复杂的证明过程可以略去不讲。这样一来，教学课时和内容安排合理，整堂课就不会全在枯燥无味的数学定理推导中度过，即使数学基础不是很好的学生也能掌握并运用。而且学生能运用定理判断一种迭代法是否收敛，本身也会获得一定程度的满足感和自信心，而学习兴趣也可以在这之上建立起来。

2.对学生的计算机编程能力要求。该课程研究的是几大类数学问题的数值算法，懂得算法之后，一定要结合计算机进行编程实现。但本门课程又不是专门的计算机编程课程，且由于学时限制，课堂上不可能有多余的时间教授学生编程知识，因此该课程的先修课程还需要掌握一门计算机编程语言。现今的主流商用数学软件主要有如下几种：Matlab、Mathematica、Maple、MathCAD、C/C++、Fortran等，应选用一种熟悉或较易掌握的软件将各种算法进行计算机实现。另外，也可选用如Mathematica这类商用软件进行编程，该类软件界面简洁，语言简单，且功能也比较强大，自学便能很容易上手。

3.将数学理论与计算机相结合。在课堂上利用数学软件，绘制出直观的可视化图片，这样可以把课程中涉及的抽象原理、方法以及复杂的计算过程直观地呈现出来，使学生对相关算法有更直观和清楚的认识，更容易理解，同时也可加强学生运用数学软件进行编程计算的能力。如对非线性方程求根之前，先要找出有根区间，这时可以运用数学软件先画出函数曲线图，找出有根区间的大概位置，然后选择某一算法编程，观察根在迭代过程中的收敛性特征；又例如讲解最小二乘法拟合曲线时，可以运用数学软件将拟合出来的函数图与原函数表对比，可更加直观地理解插值和拟合函数中存在的误差。

4.课程中穿插实践环节，让学生参与到课堂中来。某一算法或某类问题讲解完后，应举出一些算例，让学生在课堂上分组讨论解决的办法，选择怎样的算法合适，怎么编程实现等。对于一些相对较简单的问题，可以请学生直接在课堂上编程求解并运行结果，然后一起讨论该结果的可靠性，或者对编程和运算过程中出现的问题怎么改正等。让所有的学生都参与到课堂中来，以提高学生学习的兴趣，而且同时能提高学生当场解决问题的能力、语言表达能力、计算机编程能力等。

5.课堂教学生动多样化。教学时应充分利用多媒体提高教学效果。如在PPT中增加声音、图像、动画等多种形式，在教学过程中形成多种感观刺激，使原来学生误解的枯燥、抽象的数学课直观化、形象化、生动化，充分激发学生的学习热情，从而大大地提高学生汲取知识的效率。另外，还可以将教学方式多样化，避免教师“满堂灌”、“唱独角戏”的尴尬局面出现。除教师讲解外，还可让学生一起参与到课堂中来，如分成小组讨论某个算法的优缺点，某个具体问题的解法，或采用小组竞赛模式，针对某一问题看谁的算法简洁、效率高、结果可靠等。

6.选择学生感兴趣的算例。算例的选择应有特点，或与学生专业相关，或与学生感兴趣的事物相关，而不应该是单纯的数学习题，应联系相关的背景或出处。如对于车辆专业的学生，讲述曲线拟合这部分内容时，可以计算汽车车身外形曲线轮廓线为例讲述曲线拟合的过程，那么可先给出一些典型车型的外形轮廓图，然后针对某款车型，给出其轮廓线上某些型值点的数据表。学生在看到丰富多彩的汽车图时，首先会感到眼前一亮，兴趣马上会提高，课堂气氛也会得到活跃，而曲线拟合的知识也能很容易领会。

四、总结

第6篇：数值方法范文

关键词：ABAQUS；位移约束；海底管道

中图分类号：P752 文献标识码：A 文章编号：1007-9599 （2012） 17-0000-02

1 工程概述

海底管道铺设是海洋油气工程建设的一项重要内容。海底管道铺设的方法基本可以分为两类：铺管船法[1，2]和拖管法[3]，其中铺管船法包括S型铺管船法、J 型铺管船法、卷管式铺管船法；根据管道所处位置不同，拖管法分为水面拖、水下拖、近底拖和底拖。对于登陆段海底管道采用底拖法施工更具有可行性。对于底拖法施工可以在陆地焊接后，由陆至海利用绞车、绞盘、拖轮等设备牵引铺设；也可以在铺管船上焊接，由海至陆铺设，其牵引方法有：岸上设置绞车牵引，利用铺管船上的绞车反向牵引。

某登陆段管道采用底拖法铺设，在铺管船上焊接管道，岸上设置定滑轮，由铺管船上的绞车带动管道铺向岸边，见图1。

该方案中，铺管船到海床段的海底管道形成S型，管道受到拖管力、张紧器张力、自重、自身浮力、浮筒浮力、海床支撑力、海床摩擦力等载荷作用，为了底拖施工的安全进行，进行管道强度校核是十分必要的。以下给出了管道强度分析的关键参数：

管材为X65钢，钢管外径为813mm，壁厚22.2mm，钢管外敷防腐涂层，厚度2.8mm，防腐涂层外为混凝土层，厚度为80mm，管道长度总长575m，海床上管道长度为375m；海床摩擦系数为1.0；水深14m；张紧器张力为100kN；拖管力为350kN；由于绑缚浮筒，管道水下重量为540.7N/m。

铺管船各辊轴相对位置：为了考虑边界影响，张紧器前取2个辊轴。根据工程作业的铺管船情况，在张紧器后共8个辊轴，从船艏至托管架方向各个辊抽名称分别为：R1，R2，张紧器，R3，R4，R5，S1，S2，S3，S4，S5。在水平方向和竖直方向上每个辊轴距离辊轴R1的长度见表1，其中R1距离水面的高度为3.9m。

注：托管架上各个辊轴水平向至R1的距离考虑了拖管架角度。

2 ABAQUS数值模拟

以上工程施工中，铺管船到海底段管道形成S型，为大变形问题。ABAQUS[4]软件具有强大的非线性分析功能，在工程中有着广泛的应用。根据以上参数，采用软件ABAQUS模拟管道的底拖过程。铺管船和托管架上面的辊轴和管道的作用，以及管道和海床的相互作用都可以通过接触的方式处理。众所周知，接触为非线性问题，对管道、海床和辊轴的建模有一定的要求，如果处理不当则计算难以收敛。因此，本文通过位移约束的方式模拟了管道和辊轴的接触，通过位移约束和加载的方式模拟了管道和海床的相互作用。以下给出模拟过程及计算结果。

2.1 模拟过程

第一步：建立模型，考虑管道半径，管道竖直向坐标为4.3893m，管道单元B32，见图2。

第二步：根据各辊轴位置给出管道上相应的约束点。通过移动坐标系平面的方式建立新平面，新平面和管道的交点为约束点，见图3。虽然当管道大变形后约束点和相应辊轴位置不一致，但在本文的模拟中，这种不一致对结果的影响可以忽略。着泥点的位置可以根据经验确定，或者通过调试的方法得到：首先给出着泥点初始值，计算出着泥点的支反力，然后调整着泥点的位置，当支反力为零时，对应着泥点位置。

第三步：施加约束。根据各个辊轴相对R1在竖直向的长度得到管道约束点和海床段管道竖直向位移，施加位移约束。在海管铺设中，某些辊轴并不能起到支撑的作用，计算出各约束点的支反力，当其为拉力时，则放松该约束。张紧器的拉力通过简支约束管道端部体现，其他段管道在水平向可以自由移动。见图4。

第四步：施加重力载荷。水面以上和以下管道重力不同，水面与管道交点可以通过经验得到，也可通过迭代的方式求得。

第五步：施加海床摩擦力和拖管力。以均布载荷的形式施加摩擦力，根据管道水下重力和摩擦系数，可知摩擦力为540.7N/m。拖管力取350kN。

2.2 计算结果

按照以上步骤建立模型，计算得到管道应力场，见图5。其中上弯段最大应力为297MPa，下弯段最大应力为290MPa。

3 总结

某登陆段管道采用由海至陆的底拖法铺设，本文采用ABAQUS软件建立数值模型，计算了管道应力。上弯段最大Mises应力为297MPa，下弯段最大Mises应力为290MPa，为管道底拖强度分析奠定了基础。

参考文献：

[1]E Heerema. Recent Advancements and Present Trends in Deepwater Pipe- Lay Systems. OTC 17627， 2005.

[2]Braestrup M， Andersen J， Andersen L， et a1. Design and installation of marine pipelines.Blackwell Science Ltd.， 2005， 210-238.

[3]桑运水，韩清国.海底管道近岸浅水铺设的岸拖与海拖.石油工程建设.2006（4）.

[4]ABAQUS Version 6. 7 Documentation，ABAQUA，Inc.

第7篇：数值方法范文

关键词：波浪力；SPH；开孔沉箱；水动力；线性关系；规则波；光滑函数

中图分类号：TV139.26 文献标志码：A 文章编号：1672-1683（2014）06-0078-06

开孔沉箱式防波堤利用自身消浪室内外波动的相位差和消浪室内水体的紊动，能够有效消耗波能，减少波浪与港口结构物间的反射，降低反射系数（可从实体沉箱的1.0降低到0.3～0.8左右），也可以减少作用在沉箱结构上的波浪力（一般可减小至不开孔时的70%～90%左右），加上其施工简单，造价相对较低，因而成为低反射防波堤的典型型式。

自从Jarlan[1]提出开孔式沉箱防波堤的概念以后，陆续出现大量的相关研究成果：Goda等[2]提出了最实用的研究直墙压力分布规律的公式；Takahashi[3]对该公式进行修正，并应用到开口直墙结构中，提出了开孔板前后两面及消浪室后墙压力的分布都是基于三种峰值形式组合的理论；Hendrik Bergmann等人[4]分别对单层和双层前墙开孔板进行了实验分析，得出当相对宽度比为0.25时单层开孔板的总水平力比降低最显著的结论；Kyung等人[5]提出了基于微波幅理论研究不规则波与开孔沉箱相互作用的新型数学模型，并实验分析了波浪反射与消浪室宽度、波高和周期的关系；陈雪峰等人[6]通过大量实验，得出了反射系数和总水平力与相对宽度、相对波高等因素的近似经验关系；EI-Hafid等人[7]将影响系数χ项引入到Takahashi公式中，分别用2D物理模型实验和Dieppe 实测数据等方式对新公式进行了验证；陈雪峰等人[8]采用VOF方法研究了反射系数和总水平力；刘勇等人[9]基于线性势能理论，提出了以匹配特征函数和有限元法为基础的半解析方法，研究了不规则波作用在开口式沉箱上的总水平力和垂直力；姜俊杰等人[10]通过2D规则波波压力试验，提出了有顶板开孔沉箱所受波浪力的计算方法。

但是，由于波浪与开口式沉箱结构相互作用过程中涉及到波浪破碎、湍流等多种复杂的水动力现象，特别是消浪室内流体运动的复杂特性对开孔沉箱水平力的影响，因此要准确描述该过程，仅有实验分析和理论研究是不够的，还必需建立一个考虑多尺度结构和多物理效应的流固耦合数值模型。而光滑粒子流体动力学（SPH）[11]方法为这一工作提供了新的途径。

与传统的VOF数值方法相比，SPH是一种纯拉格朗日性质的无网格粒子自适应的方法，主要用于处理大变形、跟踪运动界面或自由表面，以及获取变量的时间历程等问题，非常适合处理波浪运动的非线性问题。

本文建立了基于光滑粒子流体动力学方法（SPH）的二维流固耦合数学模型，通过黎曼解和CSPM修正后，模拟开孔式沉箱的水动力条件。基本思路为：首先，将实体沉箱波压力计算结果与线性规则波理论进行对比，验证SPH方法的解决波浪与沉箱水动力问题的准确性；其次，分析不同消浪室宽度下消浪室内外水粒子的运动特性，得出消浪室内外波压力的分布规律，探讨总水平力与消浪室相对宽度、相对波高和相对水深的关系。

1 数值计算模型

1.1 水动力控制方程

本文采用Navier-Stokes方程组来描述流体运动的质量守恒和动量守恒方程，质量守恒方程为

式中：ρ为水的密度；u为速度矢量。

通过SPH方法的基本方程可将质量守恒方程离散为SPH粒子形式，即

式中：mj为粒子j的质量；uij为粒子i和j之间的相对速度矢量；Wij为粒子j对粒子i产生影响的光滑函数，与光滑长度h紧密相关Wij=W（Rij，h），Rij=|xi-xj|h。

由SPH原理得到粒子密度的求和法计算公式

描述无黏性流体运动的动量守恒方程可以表示为

式中：p为压力；F为体积力，通常为重力加速度。

可以得到SPH形式的描述流体运动的动量方程如下式：

式中：ij=hijuijxij|xij|2+2；Cij=12（Ci+Cj）；ρij=12（ρi+ρj）；hij=12（hi+hj）；uij=ui-uj；xij=xi-xj；u和x为粒子的速度和位置坐标；αΠ和βΠ为常数。对于自由表面流动问题，通常取αΠ=0.01，βΠ=0；ij=0.1hij；C为声速。最终SPH形式的动量方程为

式中：P*ij=piρjcj+pjρici+ρiciρjcj（uRj-uRi）ρjcj+ρici，；URij=uRjρjcj+uRiρici+（pj-pi）ρjcj+ρici；ρi、pi、ci、ρj、pj、cj分别为相互作用的两粒子的密度、压力、声速；uRi、uRj为两粒子速度ui、uj在其连线方向的投影值。

本文中的沉箱的边界采用虚粒子法，根据牛顿第三定律原理，可得出水粒子与沉箱粒子相互作用时，进行力和能量的传递，进而可以得出沉箱粒子受力。

数值波浪水槽造波板及如何消除波浪的二次反射相关理论内容参见高睿等人[16]所述。

2 数值计算域和水动力验证

本文采用图1所示2D数值计算域和开孔沉箱模型，坐标原点设置在静水面与开孔板相交处。波浪水槽尺寸为6 m×1.5 m，水槽右边界设置开口式沉箱，沉箱高h=1.0 m，沉箱宽度B1=0.45 m，开孔率ε=0.4（ε=S开孔沉箱开孔面积/S消浪室迎水面面积）。数值计算参数见表1，取值点分布见图2。流体计算域采用SPH方法模拟，粒子初始间距为0.01 m，约32 245个。

为了消除反射影响，在右端设置人工黏性消波层，用来吸收入射波。数值水槽周期T=1.4 s，波高H=0.12 s，共计算10 s。以距离造波板1 m的静水位位置作为观察点，观察其波压强历时规律，见图3。

从图3可知，粒子间距为0.02 m时，波周期发生明显不稳定的现象，且波谷位置出现毛刺；粒子间距为0.01 m时，初始波压强较粒子间距0.005 m偏大，但0.5 s后，两者基本一致。考虑到计算效率和成本，本文计算时采用粒子间距为0.01 m进行数值模拟。

为了检验SPH方法模拟波浪与沉箱相互作用时，沉箱结构压强计算的正确性，取开口率ε=0时极限情况的实体沉箱，将水槽粒子静水压强沿水深分布与水力学理论进行对比。计算中，实体沉箱宽度B1=0.45 m，波高H=0.10 m，波浪周期T分别为1 s和1.2 s。静水压强沿水深的分布见图4，同时，图5中给出了T为1.0 s、H为0.08 m，10 s时水粒子运动的形态，水粒子分布较均匀，初始情况下粒子压强值大小也比较准确。

同时，本文提取了实体沉箱动水压强沿水深分布的计算值，并与李玉成等人[17]所述基于线性波理论直墙上的波动压强沿水深分布的解析解进行对比；水粒子作用下实体沉箱波动压强沿水深的分布见图6，纵坐标自下至上为水底至静水位无量纲化后的z值，横坐标为无量纲化后的波动压强P。

由图6可知，在长波作用下，实体沉箱的波动压强分布较为均匀，而在短波的作用下，实体沉箱的波动压强分布随着水深的增加减小的趋势更为明显。SPH方法的波压强的计算值与线性波理论值是一致的，因此采用SPH方法数值模拟波浪与防波堤的相互作用是可行的。

SPH方法最大的优势是能够解决大变形问题，以及模拟波浪的非线性特性。图7给出了周期T为1.0 s，波高H为0.12 m，消浪室宽度B分别为0.3 m时，水粒子运动10 s时，实体和开孔沉箱消浪室内外水粒子压力分布和速度矢量。图7（a）中，水粒子在实体沉箱半水深位置形成漩涡流动；图7（b）中，消浪室内水粒子向下运动剧烈，反射波与入射波的叠加比较明显，在开孔沉箱前形成驻波。可见，SPH方法可以比较形象直观地模拟波浪与开孔沉箱相互作用这一复杂的水动力过程。

3 消浪室宽度对波压力分布的影响

作用在开孔沉箱上的波压力，沿着波浪传播的方向可分成三个部分：开孔板外侧；开孔板内侧；消浪室后墙。为了方便与现有理论进行对比，在此所讨论的波压力是特定位置的点，开孔板外侧取1、2、6点，开孔板内侧取7、8点，消浪室后墙取10、11点。通过取值计算，分别得出消浪室宽度B分别为0.15 m、0.2 m和0.3 m时，波浪周期T为1 s、波高H为0.08 m时波压力分布（图8）。为了便于比较，图8中还给出了实体沉箱波压力分布，以及开孔沉箱外侧压力分布规范解。

从图8可以看出：实体沉箱波压力分布是呈梯形分布的，在静水位位置的波压力最大；消浪室宽度为0.15 m时，开孔板外侧波压力值与实体沉箱波压力值基本相等；消浪室宽度为0.2 m时，波压力值有明显的减小；消浪室宽度增加到0.3 m时，开孔板外侧波压力值最小。由规范可知，零压力点应位于静水位上1.275倍波高位置，即0.102 m处。本文计算得到波压力零点为0.04 m。

上述分析依据的是有限的取值点，下面对能够与水粒子接触所有沉箱边界进行分析。从实体沉箱及开孔沉箱波压力分布图（图9）可以看出以下现象。

（1）实体沉箱波压力最大值较静水位偏下，并且分布形式大致呈梯形，但不完全是线性的。

（2）消浪室后墙波压力分布规律大致与实体沉箱一致，且消浪室宽度越大，波压力值越小。

（3）1号板外侧波压力值随水深减小而逐渐变大，基本呈线性分布，且消浪室宽度越大，非线性越明显。

（4）与板外侧相比，1号板内侧波压力分布，波压力随水深减小增加的更为明显。

（5）与板内侧相比，开孔沉箱2号板内外侧波压力分布的非线性更强。

4 总水平力影响因素的研究

4.1 开孔沉箱总水平力与影响因素关系的实验对比

根据刘勇等人[9]研究，开孔率对开口式沉箱垂直力影响较大，对水平力影响较小（约为4%）。因此，本文忽略开孔率对水平力的影响，仅采用开口率ε=0.4时进行分析。开孔沉箱模型设计为前墙开孔，开孔幅度为水下0.3 m至顶。规则波周期T分别取1.0 s、1.2 s、1.4 s；波高H分别取0.08 m、0.1 m、0.12 m；消浪室宽度B取0.15 m、0.2 m、0.3 m。为了便于比较，还模拟了前墙为实体的沉箱受力情况，其波要素与前墙开孔时相同。每种工况组合数值模拟计算三次，取其平均值进行分析。为了便于与已有试验进行对比，沉箱固定点水平压力取值点为图2中的第2点，波长的近似计算采用Katsardi[18]提出的经验结论，分别用Fk和Fs表示开孔沉箱和实体沉箱所受到的单位长度总水平力，以Fk/Fs为纵坐标作为分析开孔沉箱总水平力指标，反映开孔沉箱对总水平力的影响。对Fk和Fs进行无量纲化处理，无量纲因子为容重γ、波高H、水深d和沉箱Y方向宽度b。分析结果见图10。

图10（a）是H和B不变，而改变波长L的情况，图10（b）是B变化，而L和H不变的情况。很明显，两种情况下，相对宽度（B/L）与单个点总水平力Fk/Fs之间的关系近似呈线性关系。图10（c）显示，只改变H而B和L保持不变时，相对波高（H/L）与总水平力之间也是近似呈线性关系；图10（d）显示，B和H不变而改变L时，相对水深（d/L）与总水平力比之间同样仍近似呈线性关系。

陈雪峰等人[6]针对开孔沉箱进行了大量实验，研究认为B/L对水平力比的影响最大，H/L的影响次之，d/L的影响最小。陈雪峰等人[8]提出的开孔式沉箱总水平力比的实验公式为

图11给出基于SPH方法得到的数值解析解与上述实验公式解的对比结果，可以看出，SPH方法对开孔沉箱水平力的计算与实验公式的结果基本吻合，大多数点在y=x（1±10%）的包络线内。可见，在类似本文数值域和波浪要素的条件下，SPH数学方法的计算结果可以作为分析开孔沉箱所受总水平力大小的参考依据。

4.2 开孔沉箱总水平力相位差和非线性研究

在实验中，仅根据有限点取值用来分析开孔沉箱总水平力的研究是不精确的，因此本文考虑了能与水粒子接触所有的沉箱边界，来研究的波浪力与相应影响因素的关系。一般情况下，沉箱总水平力F总由F前、F内、F后三个部分组成，图12给出了作用力方向和前述工况下总水平力的历时曲线，

由图12可知，F总的峰值总是出现在F前和F后最大值之间，即F总总是出现在波峰进入消浪室后，接触到消浪室后墙之前；F前和F后的相位差随着消浪室宽度的减小而减小。很明显，随着消浪室宽度的减小，四个力是趋于同相位的。此结论与刘勇[9]的实验结论一致。

通过进一步分析相对宽度B/L、相对波高H/L和相对水深d/L对总水平力的影响，以Fk/Fs作为分析指标，计算工况与5.1节相同，分别得出对应的四种结果（图13）。

由图13（a）可知，当B/L较小时，总水平力有小幅的增加；而B/L达到0.11左右时，总水平力有减小的趋势，即B/L对总水平力比的影响非线性，这一点在物理模型实验尚未发现。图13（b）反映同样规律。图13（c）中波高H变化时，相对波高H/L对开孔沉箱总水平的影响也是非线性的，但是当波高比较大、波高H为0.12 m时，H/L对总水平力比的影响却是呈现线性变化。图13（d）显示，相对水深d/L对总水平力是非线性的。

由于开孔消浪结构的消浪机制比较复杂，其尺度的选择与波浪力的计算只能依赖于理论分析、数值计算和实验研究相结合的办法来确定。从已发表的文献来看，理论分析和数值计算可以得出反射系数，而波浪力的分析至今较多是试验结果。对于开孔沉箱波浪力的计算，日本港湾研究所高桥（1996）推荐采用修正后的扩展合田公式，其中给出了静水位和墙底的压力值，认为静水面上压力零点位置为0.7H，三点间波压力呈线性分布，并且认为波压力与相对水深d/L呈二次关系。

本文假定总水平力与相对水深d/L呈二次关系，与相对波高H/L呈一次关系，采用最小二乘法拟合，将数值结果进行逼近，来描述各因素与总水平力比Fk/Fs的关系。拟合简化关系式为

拟合的相关系数R=0.921，符合拟合方程的相关性要求。与试验研究相比，式（13）给出了不同于式（12）形式的结果，说明总水平比与其相关影响因素之间的关系可能是非线性的。因此，以波压力特征点进行总水平力研究是不精确的，同时此结论也为类似条件下计算开孔沉箱总水平力的给出了一个新的参考。

5 结论

（1）SPH方法能够较真实地模拟复杂的流固耦合作用过程，特别是消浪室内部粒子的相互碰撞以及波浪的回流过程。经过对比，数值计算得出的波压力的结果与现有的理论结果比较吻合，因此，基于SPH方程的数学模型可以被用来进行波浪与开孔沉箱研究工作。

（2）在相同波浪要素条件下，消浪室后墙的波压力随消浪室相对宽度的增加而减小。在静水位以上，波压力分布基本是线性分布的，而在静水位以下，波压力的分布出现非线性现象。随着消浪室宽度的增加，总水平力与其组分之间存在着同相位到明显相位差的转变，总水平力峰值介于消浪室前墙和后墙峰值之间，并且消浪室前墙和后墙的总水平力增大的趋势相比消浪室内侧明显。说明总水平力并不是随着消浪室宽度的增加一直会减小。当消浪室相对宽度B/L为0.11时，开孔沉箱总水平力出现减小趋势。

（3）假定总水平力与相对水深d/L呈二次关系，与相对波高H/L呈一次关系，可以得出总水平力与其相关影响因素的非线性关系式，为开孔沉箱受力分析和工程设计提供了一个新的参考依据。

参考文献（References）：

[1] JARLAN GE.A Perforated Vertical wall Breakwater[M].Dock Harbour Auth.1961，XII（486）：394-398.

[2] GODA Y.Random seas and design of maritime structures[M].Press of University of Tokyo，1985：302-306.

[3] TAKAHASHI S，SHIMOSAKO K.Wave pressure on a perforated caisson[M].Port and Harbor Re-search Institute，Proc .Hydro-Port Yokosuka，1994，1：747-764.

[4] HENDRIK BERGMANN，HOCINE OUMERACI.Wave loads on perforated caisson rreakwaters[J].Coastal Engineering，2000：1622-1635

[5] KYUNG DOUG SUHA，JAE CHUN CHOIB.Reflection of irregular waves from perforated-wall caisson breakwaters[J].Coastal Engineering，2001：141-151.

[6] CHEN Xue-feng，LI Yu-cheng，SUN Da-peng.Regular waves acting on double-layered perforated caisson[C].International Offshore and Polar Engineering conference，2002：736-743.

[7] El-HAFID TABET-AOUL，ELOI LAMBERT.Tentative new formula for maximum horizontal wave forces acting on perforated caisson[J].Journal of Waterway，Port，Coastal and ocean Engineering / January/Feberary，2003：34-40.

[8] CHEN Xue-feng，LI Yu-cheng，Wang Yong-xue.Numerical simulation of wave interaction with perforated caisson breakwaters[J].China Ocean Engineering，2003，1（17）：33-43.

[9] LIU Yong，LI Yu-cheng，TENG Bin.Total horizontal and vertical forces of irregular waves on partially perforated caisson breakwaters[J].Coastal Engineering，2008：537-552.

[10] JIANG Jun-jie，LI Yu-cheng.Experimental study on vertical wave forces of roof perforated caisson under regular waves[J].Port and Waterway Engineering，2011：136-142.

[11] LIUG R，LIU M B.Smoothed particle hydrodynamics：A mesh free particle method[M].Singapore：World Scientific，2003.

[12] Colagrossi A，Landrini M.Numerical simulation of interfacial flows by smoothed particle hydrodynamics[J].Journal of Computational Physics，2003，191：448-475.

[13] MONAGHAN J J.Smoothed particle hydrodynamics[J].Annual Review of Astronomical and Astrophysics，1992，30：543-574.

[14] RANDLES P W，LIBERSKY L.Normalized SPH with stress points[J].Int.J.Number.Methods Eng.2000，48（10）：1445-1462.

[15] PARSHIKOV A N，MEDIN S A，LOUKASENKO I I.Improvements in SPH method by means of inter particle contact algorithm and analysis of perforation tests at moderate projectile velocities [J].Int.J.Impact Eng.2000，24（8）：779-796.

[16] 高睿.SPH强非线性水动力学数值模型的应用与改进[D].大连：大连理工大学，2011.（GAO Rui.Application and improvement of SPH strongly nonlinear hydrodynamic numerical model[D].Dalian：Dalian University of Technology.（in Chinese））

[17] 李玉成，滕斌.波浪对海上建筑物的作用[M].海洋出版社，2002.（LI Yu-cheng，TENG Bin.Function of the waves on the offshore structures[M].The Ocean Publishing Company，2002.（in Chinese））

[18] KATSARDI V，SWAN C，An experimental study of shallow water wave statistics on mild bed slopes[C].ASME Conference Proceedings，2011：711.

第8篇：数值方法范文

关键词：流变特性；锚喷支护；悬吊理论

1 工程概况

1.1 工程概况

高家梁井田位于鄂尔多斯万利矿区南部，鄂尔多斯市东南约8km处，设计年产量600万吨，是一个新开发的矿区。正在建设中的高家梁煤矿主井口底板标高1398.506m，坡度-14度，净断面积17.3m2、毛断面积19.2m2；副井筒井口底板标高1398.702m，坡度-5.5度，净断面积 17.8m2、毛断面积20.5m2；风井井口底板标高1398.318m，坡度-20度，净断面积 17.3m2、毛断面积18.9m2。顶设计成圆拱形，净断面半径为2500mm。

1.2 软岩斜井锚喷网支护数值模拟对比分析

高家梁煤矿斜井分主、副、风井，文章以副井为研究对象进行数值模拟的对比。高家梁副井的际支护为两帮及拱顶为锚喷网支护，底板每施工10m将泥化的砂岩清除，清除泥化砂岩的厚度为300mm，换为300mm厚片石砂浆基层，其上再浇300mm 厚C30混凝土面层。本节先拟用全断面锚喷网支护（底板拟喷射300mm混凝土，并设置4根锚杆）进行数值模拟，与实际采用的支护方式进行对比，分析两者的优缺点。

1.3 全断面锚喷网支护数值模拟

采用ADINA中的Native建模方式，计算区域左右两侧和上下两侧分别向外延伸斜井高度和跨度5倍，取50m×50m。支护结构包括喷射混凝土及锚杆。岩石以及喷射混凝土采用平面四节点等参单元，锚杆采用rebar单元。

在有限元计算中，边界条件分为应力边界条件和位移边界条件，应力边界条件在ADINA中由设置外载荷来实现，位移边界条件由设置模型边界约束来实现。本模型只有位移边界条件。在水平方向上，Y向的边界上关闭Y-Translation自由度（图中字母C表示），以模拟远离斜井左右的土体边界没有位移，垂直方向上，模型下表面边界关闭Z-Translation自由度（图中字母B表示），以模拟远离斜井的深层土体没有竖直方向的位移。

开挖分两次开挖，第一次挖去上半部分，然后加载支护。第二次挖去下半部分，然后加载支护；喷射混凝土厚度：0.2m，底板喷射混凝土厚度为0.3m。围岩、混凝土以及锚杆计算参数根据前文试验测定结果，整理如表1.1所示。

将制作好的试样放在10KN拉力材料试验机上进行加载，并绘制拉力―变形曲线；然后根据试验结果设计正交试验，确定最佳配比；得出结果后，重新制作试件；然后再将试样放在拉力材料试验机上进行加载并分析结果。

1.4 试验结果全断面锚喷网支护模拟后处理分析

利用ADINA里的云图，可以清楚的看到斜井位移及应力的变化情况，进而对支护前后斜井的变化情况进行分析。

开挖完成后会出现一定的应力集中现象，在锚杆支护完成后应力集中现象更加明显，通过ADINA后处理导出各点的应力数值，可知支护完成后最大有效应力为6.4×106Pa左右，最大剪应力为3.5×106Pa左右，两者均在合理变化范围之内。

通过支护前后对比分析可知，支护结构对围岩应变有显著的控制作用，尤其是底板4根锚杆对竖向应变的控制非常明显，斜井横向应变以井口为中心，左右对称，应变值相近，支护后横向最大应变为3.5×10-4，竖向最大应变5.4×10-4。

锚喷网支护对斜井围岩变形有显著的控制作用，尤其是竖向位移，开挖完成为打锚杆前，底鼓现象比较明显，底板进行锚喷网支护后有效的控制了底鼓，但顶板下沉控制不是很明显，横向位移两帮收敛值相近，以井口为中心左右对称。

2 数值模拟对比结果分析

软岩斜井采用锚喷网联合支护结构，改变了围岩内部应力，变被动支护为主动承载，原来需要一定强度支架支撑的岩体变成了支护结构的主体，充分利用了岩体的自身强度。锚喷网是一种主动支护形式，可以有效地保证斜井的使用断面。锚喷网支护形式具有机械化程度高、施工速度快、综合效益好、工人劳动强度低、能充分发挥围岩自承能力等优点。

通过全断面锚喷网支护与实际采用的锚喷网+混凝土支护的数值模拟对比分析可知，全断面锚喷网支护能有效控制围岩变形，且底板底鼓明显减小，斜井最大位移也明显减小，但通过现场实际检测，锚喷网+混凝土支护能够有效控制围岩变形，且此支护方式要比全断面锚喷网支护施工简单，且节约成本。综合考虑，采用锚喷网+混凝土的支护方法是比较理想的选择。

3 结语

文章通过数值模拟分析对比，可知拟采用的支护措施与实际采用的支护措施各有优缺点，拟采用的支护措施在控制围岩的变形上比较理想，这也是实施支护的目的所在；实际采用的支护措施会造成一定的底鼓现象，但经济效益和可推广性较强，在能够保证围岩变形在控制范围之内的基础上，实际采用的支护方法有较大的优势。

参考文献

第9篇：数值方法范文

关键词：水下航行体；变水域；数值模拟计算

0 引言

本文拟通过潜艇以相同的速度由无限水域进入不同阻塞比的有限水域和以不同的速度进入同一阻塞比的有限水域两种方法来分析在不同参数影响下潜艇进入变水域的阻力性能的规律（阻塞比β=潜艇面积/隧道面积）。通过几组参数计算结果的对比，揭示水下航行体进入变水域阻力性能与速度和阻塞比的关系，并通过此关系验证动网格技术对变域情况数值模拟的有效性，可作为后续导弹出筒等相似计算的方法依据。

1 基本理论与数值方法

1.1雷诺时均纳维叶-斯托克斯方程

雷诺时均纳维叶-斯托克斯方程和雷诺时均连续性方程：

―湍流动能生成项。其中的一些常数值如表1。

2 计算结果与分析

2.1 动网格模型

动网格区域与主体部分建模是分离的两个部分，而它们之间的数据交换是通过在动网格的外表面区域与主体部分与动网格交界内部区域之间完成的。

2.2 计算结果分析

经过数值模拟计算，可得到如图1数据。

从图1计算数据分析可以得出，潜艇表面压力在前体与后体变化比较剧烈，潜艇头部尖点处压力达到最大值，在前体与平行中体和后体与平行中体交接处变化最大，平行中体上压力值基本保持不变。在潜艇进入受限域的过程中，压力差值逐渐变大，在潜艇完全进入受限域也就是20s后压力差达到最大值并保持稳定。因此，后续计算数据主要以潜艇完全进入受限域后气动性能数据为主。

2.3 阻塞比的影响

阻塞比为潜艇横截面积与受限域横截面积的比值。当潜艇横截面积一定的时候，受限域横截面积越大，阻塞比越小。阻塞比对潜艇通过受限域时的水动力效应有很大的影响。为研究阻塞比影响，分别对阻塞比为0.08、0.12和0.2的数值模型进行了模拟计算，计算结果如图2。

在受限域内，随着阻塞比的增大，受限域内的压力值也增大。潜艇表面所受的压力也随之减小，在前体与后体不是很明显，在平行中体上受到阻塞比的影响相对较明显。在潜艇头部到达受限域至潜艇完全进入受限域的时候，随着阻塞比的增大，压力也增大。但是在潜艇完全进入受限域后情况却正好相反，阻塞比大的受限域潜艇表面最大压力值反而更小。随着阻塞比的增大，在进入受限域时潜艇阻力变化值也增大，完全进入受限域后，阻塞比对潜艇阻力的影响很小。

2.4 速度的影响

速度增大也会增大潜艇表面压力差值，压差值的增大会在潜艇表面产生比较大的扭矩，这对潜艇的结构安全性有很大的威胁。为研究速度对潜艇压力的影响，分别对速度为10m/s、15m/s和20m/s的潜艇进入受限域进行了数值模拟，得到结果如3所示。

潜艇速度对潜艇所受的压力最大值的影响特别大。随着潜艇的速度的增大，潜艇所受压力最大值的值也随之明显增大，并且潜艇受到的压力最大值的改变幅度也随之增大。也就是说，随着潜艇速度的增大，潜艇在受力最大值点的压力差值也随之增大。潜艇的速度值越大，潜艇在无限域中受到的阻力越大，在进入受限域过程中阻力差的值也越大。

3 结论

本文基于无粘、不可压缩流紊态流场的雷诺时均方程和k-？两方程紊流模型模拟了潜艇通过变水域时的水动力性能，建立了潜艇-变水域水动问题数值计算模型，应用有限体积法进行求解。对水下航行体在动网格计算法下得出的水动力性能进行了分析，得出了以下结论。

（1）利用CFD软件FLUENT建立了潜艇通过变水域的水动力学模型，并通过计算分析得出了潜艇由无限域进入受限域的水动力性能。

（2）随着阻塞比的增大，受限域内压力值也随之增大，但是潜艇表面压力值却是随之减小，在平行中体最明显；潜艇表面最大压力值的压力差值也增大，潜艇在进入受限域过程中总阻力值的最大值与阻力差值也随之增大，但是在无限域中和完全进入受限域后的总阻力值变化很小。