数值实验系统

2024-06-26

数值实验系统（精选八篇）

数值实验系统篇1

关键词：流体力学,源程序,CFD商业软件,数值实验系统

本科生和研究生开设的流体力学课程知识点容量大，概念与原理抽象复杂[1,2,3,4]，且与工程实际现象紧密结合，很多学生在学习该门专业课程时，普遍反映理论知识比较抽象，流体的运动规律难于理解，且不能与具体流动现象相结合。[5]因此，实验在流体力学教学中有着非常重要的作用，但由于条件限制，传统的物理实验常存在一些弊端和局限性。随着远程教育的不断发展，物理实验的局限性更加突现。数值实验或数值模拟具有低成本、不受时空限制等优点可以弥补传统物理实验的缺陷，作为传统实验教学的扩展和补充。

构建数值实验仿真系统，建立相关数值实验算例数据库，有利于提高学生的动手能力，数值实验有利于学生进行自主探究和独立学习，学生可根据自己的兴趣自行选择实验内容、自行制定实验时间和进度，自己发现问题、解决问题，学生由被动的知识受体转变为知识的主动探究者和学习的主体。在流体力学课程教学实践中，采用以“学生为主体、教师为主导”的任务式实验教学模式，把通过数值模拟演示能够更好地了解和掌握的内容抽取出来，基于数值实验系统平台给学生布置若干数值模拟任务，通过在计算机上的自主实践，在教师的指导下学生在完成相应任务的同时，总结学习相关知识。

可以通过对接科研成果将研究的数值模拟程序应用于数值仿真平台建设，促进研究流体力学数值实验内容体系，确定数值实验系统的经典实验项目。教师可将课程教学内容与任务式实验结合，加深学生对教学内容的理解和掌握。通过流体力学课程数值实验系统的建设，可促进该课程实验教学内容重新组织，形成新的实验教学与实践授课计划。根据学历层次及专业的不同，形成适应不同专业的任务型实验和实践项目。可探索理论教学考核与实验教学考核结合模式下的学生成绩考核方法。基于数值实验系统形成的教学模式及方法也均有推广意义。鉴于此，构建流体数值实验系统，搭建数值实验软硬件平台十分有益于流体力学课程教学。

1 数值实验系统的结构和功能

根据流体力学教学要求和安排，该数值实验系统结合流体力学的基本知识点，设置典型算例问题，有针对性地选择包含作者研发及部分开源的CFD源程序和CFD商业软件构建该数值实验系统。构建该系统旨在将流体力学的重要知识点理论与实践有机地结合起来，通过该虚拟实验系统，加深学生对流体力学基本知识的理解，通过动画模拟过程增强其对实验现象的认知，让其内容更加生动形象化，更能调动学生的学习兴趣和积极性，从而充分提高学生的动手能力、实践能力和创新思维能力，进而培养学生的科研能力。数值实验仿真系统满足学生的不同需求数值实验具有无限量复制特性，学生可以根据自己的需求，编译运行程序，调用数值实验数据库中的文件，对经典算例进行数值实验，并对计算结果进行可视化图形图像输出，并做相应的动画输出，虚拟现实中的流动，建立对该种流动的直观认识，并对数值实验结果与物理实验结果对比验证。框架结构根据上述设计思想，该数值实验仿真平台分为两大系统，其主要结构功能如图1所示。

该流体数值仿真系统包括源程序数值仿真系统和软件数值仿真系统两大子系统。各子系统均包含前处理功能模块、数值计算功能模块、后处理功能模块。源程序数值仿真子系统为作者研发的源程序和国外部分开源代码组成。其中前处理模块包括BFC网格程序[6]、二维网格生成程序[7,8]、三维结构网格生成程序[9,10]。数值计算模块包括国外学者研发的源代码Dolfyn[11],CAFFA 2D&3D[8],CAFFA 3D.MB及并行版本[10]等以及笔者研发的DLRŸBFC程序[6,12]。这些程序可在POWERSTATION或INTEL VISUAL FORTRAN编译环境中运行计算，并行版本程序可在高性能计算集群的LINUX系统中编译运行。源程序数值仿真子系统中的源程序可供FORTRAN程序设计语言掌握很好的本科生或研究生学习及用其针对相关流场算例进行数值模拟实验研究。软件数值仿真子系统由当前流行的商业CFD软件FLUENT,CFX,PHOENICS,STAR-CD,TECPLOT,ORIGIN等软件组成，该子系统亦包括前处理模块、数值计算模块和后处理模块，三大模块结合完成整个流体数值仿真过程。

通过该流体数值仿真系统，能自动设置边界条件并划分网格，读取网格数据求解，很好地进行多算例工况数值实验，尤其在系列化设计方案对比中。还可将计算结果与实验结果进行对比，验证数值计算的准确性。利用系统的后处理模块，将计算结果数据后处理可得到各种图形和动画，展示计算结果，评估设计方案，让那些不熟悉软件的工程设计人员直观地看到计算结果。从而提高学生对流体力学问题的感性认识和激发学习兴趣的角度出发，满足教学和实验手段的灵活性和内容的全面性要求。

2 构建数值实验算例数据库

为了帮助学生加深对流体力学课程中涉及的一些常见流动现象的认识，根据课程教学内容体系，基于数值实验系统平台构建几种常见流动的数值模拟教学算例库，其中既有基于源程序的算例，也有基于商业软件的算例。建立数值实验前处理、数值模拟和后处理数据文件。笔者在本科生及研究生的流体力学课程教学实践中，尝试设置了二维及三维顶盖驱动方腔流动、圆管层流和紊流、锥形渐扩管路层流和紊流、后台阶流动、正弦波壁流动、圆柱绕流、单体建筑风场等简单且有代表性的算例，其中每个算例给出流体基本参数和边界条件设置情况。并对数值计算结果进行了可视化图形图像处理，形成相关格式文件，可用于学生数值实验时参考和在流体力学多媒体教学中应用。

3 流体数值实验系统在实际教学中的应用

笔者将构建的流体数值仿真实验系统在实际教学工作中加以利用，取得了一定的效果，主要应用体现以下3个方面。一是用于辅助流体力学课程教学。工程流体力学具有较强的理论性及涉及宽广的知识点，使这门课的教与学都存在较大的困难。笔者结合多年来的教学实践，基于流体数值实验系统对流体力学课程教学内容之圆管中的层流与紊流、锥形渐扩管路系统内部流及后台阶流等典型流动进行了数值模拟，并对计算结果后处理成图片和动画，在多媒体课件中展示，增强了学生对相关流动的感性认识和理性认识。二是用于指导工程力学专业本科毕业论文工作。笔者利用流体数值仿真实验系统设置顶盖驱动方腔流、突扩管路流、建筑风场并行计算等选题，连续5年指导工程力学专业本科毕业论文10余名，获得了较好的效果，毕业论文等级均在良好以上。三是基于高性能集群和流体数值仿真实验系统面向力学专业本科生及工程流变学方向研究生指导开展创新课题研究，其中指导大学生实践创新课题1项，研究生课题2项，相关课题研究顺利。

4 流体数值实验系统在教学中的应用建议

4.1 数值实验教学地点与时间的选择

数值实验教学地点最好选择在机房，与此同时，允许学生使用自带的笔记本电脑同步进行求解。也可将数值实验教学地点选择在教室，基于高性能计算工作站和CFD并行版源程序或CFD软件，通过互联网PBS作业提交系统在线提交作业完成数值模拟任务。在流场数值模拟演示教学时，依计算机性能不同和算例问题的难易程度，计算机求解运算的时间将有不同，因此如何掌握与合理应用好这段时间尤为重要，可采取背景知识介绍和对计算结果的预测这两种方法解决。具体教学实践中可由任课教师针对某一具体的流体力学问题，采用数值实验系统中的源程序或软件演示从网格生成和边界条件设置到数值模拟整个求解过程，讲解所采用的理论以及结果分析。

4.2 精选算例建设数值实验算例数据库

精选数值模拟仿真算例安排教学，将所选数值模拟算例系统整理分类，完善相关数据文件构成算例数据库，可为今后教学参考借鉴。流体力学数值实验课堂教学演示算例的选取可遵循简单与有代表性两个原则。算例简单既能够减少计算机的运行时间，有利于教学时间的紧凑性，又能让学生易操作，方便学生将书本理论与数值实验结果对比；而有代表性的算例密切书本知识，贴近生活或工程实际，则有利于提高教学趣味性，同时开阔学生的视野。

4.3 结合创新课题培养学生创新能力

由于课堂教学时间有限，因此应在简单的教学演示算例的基础上，精心布置较为复杂或略有提高的课外任务。例如变换几何模型、改变边界条件及各种参数等，可结合本科生的毕业设计，研究生的课题任务，根据学生的编程基础合理选择CFD源程序或软件，将教师的部分科研任务交给一些学有余力并对科研感兴趣的学生完成。布置数值实验研究算例，既可以巩固课堂所学，又利于培养学生独立创新科研能力，后者非常符合目前教育界提出的加强本科生参与教师科研的精神。

4.4 数值实验与物理实验的选择及验证

数值实验也存在一些如无法获得操作技能、无法全面考虑物理实验中出现的各种影响因素等缺点，数值计算结果需要用实验结果验证。在教学中，怎样发挥它们各自的特点，如何取长补短，是值得我们深入研究的。在教学中我们要明确以下几点：一是物理实验有其不可替代性，学生具体操作技能的形成必须通过物理实验获得；二是虚实互补，将数值实验和物理实验有机地结合起来，数值实验结果可以用物理实验结果相互验证，充分发挥各自的优势；三是因材施教，教师应根据算例流动的特点、目的、学生已有的知识和能力水平以及两个实验的特点等因素来采用物理实验还是采用数值实验。

5 结束语

数值分析课程实验报告篇2

实验名称用二分法和迭代法求方程的根

成绩

一、实验目的

掌握利用二分法以及迭代法求方程近似根的方法，并学会运用 matlab 软件编写程序，求解出方程的根，对迭代法二分法进一步认识并灵活运用。

二、实验内容

比较求方程 5 0xx e   的根，要求精确到小数点后的第 4 位 1.在区间[0,1]内用二分法； 2.用迭代法1/5kxkx e，取初值00.25 x .三、算法描述

1、二分法：二分法是最简单的求根方法，它是利用连续函数的零点定理，将汗根区间逐次减半缩小，取区间的中点构造收敛点列{ }来逼近根 x.2、迭代法：迭代法是一种逐次逼近的方法，其步骤是首先给定一个粗糙的初始值，然后用一个迭代公式反复修正这个值，知道满足要求为止。

四、实验步骤1、二分法：

(1)计算 f(x)在区间[0,1]端点处的值 f(0)和 f(1)的值；

(2)计算 f(x)在区间【0,1】的中点(0+1)/2=1/2 处的值 f((a+b)/2)；

（3）如果函数值 f(1/2)=0，则 1/2 是 f（x）=0 的实根，输出根 x，终止；否则继续转（4）继续做检验。由于 f(1/2)≠0，所以继续做检验。

（4）如果函数值 f(0)* f(1/2)<0,则根在区间[0,1/2]内，这时以 1/2 代表 1；否则以 1/2 代表 0；，此时应该用 1/2 代表 1.（5）重复执行(2)(3)(4)步，直到满足题目所要求的精度，算法结束。2、迭代法

(1)提供迭代初值25.00 x；(2)计算迭代值)(0 1x x  ；

(3)检查|0 1x x |，若   | |0 1x x，则以1x代替0x转(2)步继续迭代；当   | |0 1x x时

终止计算，取作为所求结果。

五、程序

（1）二分法程序：

function y=bisection(fx,xa,xb,n,delta)

x=xa;fa=5*x-exp(x);

x=xb;fb=5*x-exp(x);

disp(“[

]”);

for i=1:n

xc=(xa+xb)/2;x=xc;fc=5*x-exp(x);

X=[i,xa,xb,xc,fc];

disp(X),if fc==0,end

if fc*fa<0

xb=xc;

else xa=xc;

end

if(xb-xa)

end

（2）迭代法程序：

function y=diedai(fx,x0,n,delta)

disp(“[

]”);

for i=1:n

x1=(exp(x0))/5;

X=[i,x1];

disp(X);

if abs(x1-x0)

fprintf(“The procedure was successful”)

return

else

i=i+1;

x0=x1;

end

六、实验结果及分析

（1）二分法：

实验结果如下：

[

]

1.0000

0.5000

0.8513

2.0000

0.5000

0.2500

--0.0340

3.0000

0.2500

0.5000

0.3750

0.4200

4.0000

0.2500

0.3750

0.3125

0.1957

5.0000

0.2500

0.3125

0.2813

0.0815

6.0000

0.2500

0.2813

0.2656

0.0239

7.0000

0.2500

0.2656

0.2578

--0.0050

8.0000

0.2578

0.2656

0.2617

0.0094

9.0000

0.2578

0.2617

0.2598

0.0022

10.0000

0.2578

0.2598

0.2588

--0.0014

11.0000

0.2588

0.2598

0.2593

0.0004

12.0000

0.2588

0.2593

0.2590

--0.0005

13.0000

0.2590

0.2593

0.2592

--0.0001

14.0000

0.2592

0.2 593

0.2592

0.0002

15.0000

0.2592

0.0001

依据题目要求的精度，则需做二分十四次，由实验数据知 x=0.2592 即为所求的根

（2）迭代法：

实验结果如下：

软岩巷道锚杆托盘数值实验篇3

锚杆支护技术已在国内外得到普遍应用, 并成为巷道支护的主要技术手段[1]。回采巷道主要采用组合锚杆支护技术, 以锚杆为主体并辅以托盘、钢带等组合构件, 利用钢带或钢筋托梁等组合构件将各根单独作用的锚杆有效联结起来, 共同支护巷道围岩, 提高锚杆支护的整体效果[2]。特别在软岩回采巷道支护过程中, 锚杆托盘的对巷道支护起着十分重要的作用, 不同尺寸的托盘对软岩的控制效果大相径庭, 因此, 利用数值实验方法研究不同类型锚杆托盘的受力特征及其对围岩的加固机理, 将为现场工程实践锚杆托盘选型提供科学依据。

2 锚杆托盘数值模型

数值模拟效果清楚直观, 尤其是可以方便地调节锚杆支护参数而得到不同的试验结果, 易于分析锚杆加预应力后的作用机理。巷道围岩锚固后数值模拟的技术思路是:在巷道周边取一分离体, 在平面应变情况下, 研究预应力全长锚固和端头锚固锚杆的作用机理, 分析锚固体在不同的锚杆支护强度 (密度、预应力) 下力学参数的变化, 假设围岩是匀质体。相当于从巷道一侧或顶板“取出”的岩体, 其尺寸的大小应能代表锚固体的一般状态, 有关研究证实, 大部分岩石在 >0.5～1.0m时, 性状即稳定, 可认为该尺寸以上即代表岩体 (A-试件横截面积) 。选定模型的原型尺寸为3m×4m×3m, 锚固分离体所处位置、形状及尺寸如图1所示。

建立模型采用以上所提出的锚固分离体模型, 见图2, 模型坐标采用直角坐标系, 坐标原点在模型左侧面 (相当于巷道一侧自由面) 的中心。X轴沿巷道轴向, Y轴沿巷道断面的水平方向, Z轴沿巷道断面的垂直方向。计算模型取为4m×6m×4m的六面体。模型共360000个单元。锚杆长2.5m, 沿Y方向布置于模型的中心, 外锚头 (锚尾) 位于坐标原点。端锚锚固, 锚固段长1.0m, 自由段长1.5m。采用衬砌 (liner) 单元来模拟托盘, 托盘大小0.15m×0.15m, 0.2m×0.2m, 0.3m×0.3m, 厚度10mm、15mm。

通过对于不同大小尺寸托盘锚杆分别施加大小不同的预应力, 研究不同的锚固力作用下 (40kN、50kN及60kN ) 端头锚固锚杆轴力以及托盘应力分布规律。

3 物理力学参数确定

应用科学有效的方法, 较为合理地确定煤岩体的物理力学性质参数这一问题, 国内外许多工程技术人员已做了相关的研究, 许多研究都是集中于寻找在煤岩体性质参数和煤岩块性质参数之间建立一种比例关系。一般情况下, 煤岩体的性质参数都要小于煤岩块的性质参数, 所以它们的比值一般都要小于1, 这样就定义这个比值为减小因子。而对于泊松比, 煤岩体的取值要比煤岩块大, 它们之间的比值定义为增大因子。在参考、吸收国内外研究成果的基础上, 采用表1所示的煤岩体性质参数减小因子。

模型介质的力学参数和煤相符, 此处选取的回采巷道介质为中等强度的煤体, 模型介质的物理力学参数见表2, 锚杆选用直径为20mm的螺纹钢, 由于钻孔直径与锚杆直径之差为6～10mm为宜, 以7～8mm为最好, 所以钻孔直径选为28mm, 锚固剂为树脂锚固剂, 锚杆用cable单元模拟, 其物理力学参数见表3和表4。

托盘材料为A3钢, 屈服极限240MPa, 强度极限400MPa, 弹性模量200GPa, 泊松比为0.25。

4 锚杆托盘数值模拟分析

4.1 锚杆对岩体力学性态影响

(1) 由图3可以看出:预应力锚杆张拉力在锚杆托盘周边的煤岩体中形成一个压应力集中区, 可称为表层压缩区, 该应力集中区的范围沿锚杆径向约为1-1.2D左右 (D为托盘边长) , 沿锚杆轴向为0.3～0.4m;在内锚头周边煤岩体为拉压应力交汇区, 拉应力区似“杯形”, 除拉应力区外, 模型内全为压应力区所覆盖, 说明加锚后改善了煤岩体的应力状态, 锚固效果是明显的。

(2) 随着预应力值的增加, 煤岩体表层压缩区的延伸范围表现出增大的趋势, 增大范围不大;压应力量值影响的范围随锚固力的增加而增加。另外, 随着预应力值的增加, 煤岩体内部的拉应力区范围变化不大。由此可以说, 锚杆预应力值的提高将进一步影响托盘的作用, 使得托盘更为挤紧和压密表层煤岩体, 限制煤岩体的变形并改善煤岩体的力学性状, 初锚力的增大对改善表层煤岩体的力学性质作用较为显著, 但对深层煤岩体作用不大。

(3) 因为托盘的弹性模量远远大于模型介质的弹性模量, 所以托盘相对模型介质来说可以视为刚体介质, 使得托盘下的煤岩体的轴向应力等值线形状为“钟形”, 这也正是模拟托盘的一个体现。

(4) 在树脂锚固端岩体出现明显的拉伸现象, 且最大拉伸值出现在锚杆中部 (即树脂锚固段靠近巷道表面一侧) , 随着向岩体深部延伸呈减小趋势。同时, 随着预应力的增加树脂锚固端的应力值亦呈增加趋势, 但增幅相对较小。

4.2 锚杆轴力特征

图4～图6为同一托盘尺寸不同预应力下锚杆轴力分布特征曲线, 可以看出, 不同预应力大小作用下锚杆的轴力分布形态及特征基本一致, 随着预应力值的增加, 锚杆的轴力相应增加。由此可见, 锚杆的初锚力对锚杆的工作阻力影响较大。但是, 随着预应力值的增加预应力的损失也增大。这是因为岩体的变形随着锚杆预应力值的增大而增大, 那么初始张拉程度大的锚杆最终的收缩量也大, 都存在不同程度的预应力损失。这里的预应力损失主要原因是岩体的徐变 (受荷区内岩体内部结构各个组成单元在应力作用下, 将产生塑性压缩或相对变位, 且随时间变化而变化) 。

随预应力值增大图中锚固段的斜率也在同时增大, 但可以看出在深层岩体范围内的锚杆上的轴力分布形式衰减的很快, 有趋于一致的趋势, 也就说明, 预应力值的增大对岩体内部的受力状态的扰动和改善越来越小, 即增大预应力对改善深层岩体的力学性态作用较为有限。故在采用预应力锚杆对回采巷道进行加固支护时, 不能无休止的增大预应力, 可以在强度允许施加的范围内, 在设计时采用较小预应力, 但必须符合规范要求。因此, 针对回采巷道软弱岩层宜选用较适中预应力来锚固围岩。

4.3 托盘应力特征

(1) 图7为托盘应力特征, 其中左图为托盘正应力云图, 右图为托盘剪应力云图。从正应力图中可以看出, 托盘中部受到锚杆拉力作用, 应力呈正值, 该压应力范围尺寸较小, 仅处于托盘中心50mm半径范围内;在该压应力作用范围之外, 存在一个岩体挤压的负应力区域, 该区域为椭圆状;在该负应力区域外侧存在明显压应力环, 但应力值相对托盘中心压应力较小。剪应力分布云图表明, 托盘中心100mm半径内均呈现明显的剪切特征, 其中, 应力值较高区域呈现环状分布, 该应力环半径约为50mm, 其他剪切范围内应力值较小。

(2) 对图7分析可知, 在预应力值与托盘厚度不变时, 随着托盘长宽尺寸的增大, 托盘中心受压应力值逐渐减小, 该应力的影响范围较小, 正应力应力环面积亦呈增大趋势, 但应力环的应力值变化较小;剪应力区域影响范围并为发生变化, 但最大剪应力值有所增加, 幅度较小。可以看出, 随着托盘支护面积的增大, 锚杆对围岩的控制范围亦增加, 但其正应力值增加幅度并不大, 所以, 在巷道浅部围岩较为软弱松散时采用大托板支护之必要的, 即可增大锚杆的支护面积, 起到控制更大范围围岩的效果, 也可提高锚杆的锚固效果, 避免了小托盘内陷的缺点, 防止了锚杆失效。

(3) 同时, 在预应力值和托盘长宽尺寸不变时, 随着托盘厚度的增加, 正应力环面积呈增大趋势, 当托盘厚度为10mm时, 应力环仅分布于托盘内部, 当托盘厚度为15mm时, 应力环与托盘边界呈相切或相割状, 且面积增大近一倍多;剪切应力环应力值随着托盘厚度的增加而减小, 由此可看出, 适当增加托盘厚度对托盘受剪状况有改善作用。

(4) 在托盘尺寸相同, 随着预应力的增大, 托盘中心受压应力逐渐增大, 但影响范围变化较小, 同时, 在托盘尺寸为 mm和 mm时, 应力环面积随着预应力的增加而减小, 而 mm托盘应力环变化则不明显;剪应力特征表明, 剪应力环面积未发生变化, 但应力值随着预应力的增加而增大。

6 结论

(1) 预应力锚杆通过预应力张拉在锚杆托盘周边的煤岩体中形成一个压应力集中区, 称为表层压缩区, 在内锚头周边煤岩体表现为拉压应力交汇区, 除拉应力区外, 模型内全为压应力区所覆盖, 说明加锚后改善了煤岩体的应力状态, 锚固效果是明显的。

(2) 随着预应力的增大, 煤岩体表层压缩区的延伸范围表现出逐渐增大的趋势, 压应力量值影响的范围随锚固力的增加而增加。锚杆预应力的提高将进一步影响托盘的作用, 使得托盘更为挤紧和压密表层岩体, 限制煤岩体的变形并改善岩体的力学状态, 初锚力的增大对改善表层煤岩体的力学性质作用较为明显, 但对深层岩体作用不大。

(3) 随着托盘厚度的增加, 锚杆自由段的轴力有所增加, 但增幅不明显;在预应力和托盘厚度相同的情况下, 锚杆轴力并未随着锚杆托盘尺寸的增加而变化。随着预应力的增加, 锚杆自由段的轴力曲线变化出现不均衡, 但整体的变化趋势不变。

(4) 托盘中部受到锚杆拉力作用, 法向应力呈正值, 该压应力范围尺寸较小, 仅处于托盘中心50mm半径范围内;外侧存在明显压应力环, 但应力值相对托盘中心压应力较小。托盘中心100mm半径内均呈现明显的剪切特征, 应力值较高区域呈现环状分布。

摘要：利用数值实验方法, 对软岩巷道条件下不同类型托盘力学性能进行了模拟研究, 结果表明, 合理增大托盘尺寸, 既可增大锚杆的支护面积, 起到控制更大范围围岩变形的效果, 也可提高锚杆的锚固效率, 避免了小托盘内陷的缺点, 防止了锚杆失效。

关键词：锚杆托盘,锚固效率,数值实验,预应力

参考文献

[1]康红普.煤巷锚杆支护成套技术研究与实践[J].岩石力学与工程学报, 2005, 24 (21) .

[2]康红普, 王金华, 等.煤巷锚杆支护理论与成套技术[M].北京:煤炭工业出版社, 2007.

“偏微分方程数值解法”实验设计篇4

针对学习内容设计适当的程序, 不仅能够巩固基础理论知识, 还能够使学生了解到如何利用学科知识解决实际问题, 有利于提高学生学习的兴趣, 促使之将所学知识直接用于自己的研究领域. 下面主要从理论验证和学科具体应用两个不同的侧重点来设计实验的内容.

第一部分主要针对有限差分方法, 对三种不同类型的方程, 即双曲型、抛物型以及椭圆型方程整理课堂演示所需数值实验, 布置数值实验作业, 以巩固学习效果. 所依据的教材为清华大学出版社《偏微分方程数值解》 ( 第二版, 陆金甫、关治[2]) , 以及科学出版社《偏微分方程数值解法》 ( 第二版, 孙志忠[3]) , 部分代码参考《MATLAB语言常用算法程序集》[1].

第二部分主要为联系具体学科所涉及的偏微分方程求解问题进行计算. 由于时间以及本人知识有限, 仅选取了一个与水资源环境相关的实例[4].

一、理论验证

1. 双曲型方程的有限差分方法

( 1) 一阶双曲型方程不同格式计算结果对比. 对

比较迎风格式、Lax-Friedrichs格式、Lax-Wendroff格式以及Beaming Warming格式的计算效果, 可以观察到迎风格式改变方向解不收敛, 网格比不满足要求时解不收敛, 二阶格式会出现震荡.

布置练习题:

试用迎风格式、Lax-Friedrichs格式以及Lax-Wendroff格式取h = 0. 1, λ = τ/h = 0. 5, 计算到t = 0. 5, 1.

( 2) 二维问题分数步长法.

由于高维问题的计算稳定性限制的更加严格, 而分数步长法可以放宽稳定性条件. 因此下面对二维的对流方程

具体格式为:

计算结果如下图.

( 3) 二阶双曲型方程的显格式.

当网格比满足r =aτ/h= 0. 25 < 1, 数值解稳定, 空间步长和时间步长同时缩小到原来的1/2时, 最大误差缩小到原来的1/4.

2. 抛物型方程的有限差分方法

( 1) 扩散方程显隐格式以及CN格式计算结果对比. 比较显格式与加权隐格式以及全隐格式稳定性要求的差别, 使学生从直观上了解使用显式格式时一定要注意到网格比对稳定性的影响, 要求aτ/h2≤1/2, 否则解不收敛. 对加权隐格式选择θ =1/2的加权因子能够显著提高求解的收敛速度, 从而可以使学生理解以Taylor展开作为工具的误差分析的重要性.

显格式数值解在t = 0. 1时的误差| u ( xi, tk) - uki| .

( ⅲ) 显格式当时间步长减少到原来的1/4, 空间步长减少到原来的1/2, 误差减少到原来的1/4, 收敛速度为O ( τ + h2) .

CN格式当时间步长和空间步长同时减少到原来的1/2, 误差减少到原来的1/4, 收敛速度为O ( τ2+ h2) .

布置练习题:

3. 椭圆型问题的有限差分方法

( 1) 两点边值问题: 这是一维的二阶椭圆问题. 通过步长的逐步减半, 可以考察椭圆问题的收敛速度为二阶.

( 2) 二维Possion问题五点差分格式, 应用五点差分格式计算如下问题:

可以观察到:

当x与y方向步长减少到原来的1/2, 误差减少到原来的1/4, x方向与y方向收敛阶均为二阶.

二、为了使学生了解本学科知识如何运用到其专业所遇到的问题上, 我们以一位本校地下水专业老师的论文为例, 讨论了如何将模型概化为一个偏微分方程, 以及如何求解, 如何结合其专业对解作出合理的解释. 这个例子简单但是很具有代表性, 一个研究滨海不同距离的地下水位变化的情况的问题, 最终归纳为一个一维的抛物问题.

式中: H为以平均海平面为基准面的水位标高, m; S为承压含水层的储水系数, 无量纲; T为承压含水层的导水系数, m2/ h; x为距海岸的距离, m; t为时间, h; H0为海潮变幅的一半 ( 即振幅) , m; t0为海平面波动周期, h; L为内陆边界距海岸的距离, m; HL为距海岸L米处的地下水位, m.

用隐格式离散为:

线性代数方程组采用传统的求解三对角方程组的追赶法, 最终将解随时间的变化曲线画出来, 可以分析出滨海水位与海潮一样具有相同周期, 但时间上滞后, 其水位变化幅度与海岸线的距离具有衰减性, 越远水位变化幅度越小.

这个具有实际应用背景的数值算例的引入, 向学生展示了将实际问题概化为模型问题, 求数值解, 并将所求数值解结合原问题进行解释的过程, 为学生提供了利用学科知识解决具体问题的范例.

摘要：“偏微分方程数值解法”是计算数学基础课, 它的发展离不开计算机, 因此需要紧密结合程序来完成对算法的实现.本文针对不同类型的方程精心挑选了数值算例, 首先通过课堂演示来加深对理论知识的印象, 并提高学生兴趣;其次布置数值实验操作促进学生真正了解不同数值格式;最后引进一个建模实例, 为其知识未来的应用打好基础.

现代数值计算的实验教学探索篇5

随着计算机应用技术的发展, 人们开始意识到数值计算方法已经广泛地应用于物理学、力学、航空航天、土木工程、机械工程、风险投资、经济管理等众多领域。《数值计算》是大学教育开展的适合理工类学生利用数值技术求解各类工程问题的一门课程, 其理论抽象但实践性很强。然而传统的教学模式注重理论推导和定理、公式的证明, 很少会利用数学软件进行相应的实验训练, 因此, 学生在课程学习完成后只掌握了数值计算中的基本定义、定理和算法原理及推导过程, 而缺乏运用数值方法去解决实际问题的能力。近年来, 各高校逐渐重视对学生数值计算实操能力的培养, 在数值计算课程中开设了上机实验课时。在实验教学课堂上, 一方面, 由于老师们习惯了进行传统的理论方法教学, 往往抓不住实验训练的重点, 教学模式单一, 盲目地给学生布置上机训练题进行练习, 然后进行讲解;另一方面, 学生经常抱怨上机训练课时太少, 实验内容太多而且复杂, 面对着陌生的数学软件操作界面却迟迟没法动手操作。

数值计算实验是要利用计算机技术解决实际工程项目中抽象出来的数学问题, 其关键并不在于数学软件如何使用, 而是在于如何把数值计算方法转变成为计算机程序代码, 实现自动化计算过程, 这其中主要需要学习掌握的是理论方法转化为程序代码的技术, 也就是我们常说的算法设计。

二、算法设计教学的重要性

《数值计算》课程重点是培养学生运用数值方法解决实际问题的能力, 而解决实际问题必须结合计算机技术进行实验计算, 那么, 教学的关键在于教会学生如何将实际问题抽象成数学模型, 利用计算机技术完成数学建模并求解的过程, 也就是把数值方法转化成为程序代码的算法设计的过程。

很多实际问题抽象成为数学问题以后, 建立模型的求解过程需要用到各种数值计算方法。因此, 有必要将数学建模思想与数值计算实验进行有机融合。在《数值计算》的实验教学环节中, 应结合实际的数学建模案例进行数值方法的讲解, 并通过算法设计实现计算机编程, 课件算法设计在实际解决问题方面起着至关重要的作用。

现代数值计算实验教学的主要教学内容是如何把理论方法转化为计算机程序来实现运算, 也就是算法设计, 通过有目的、有重点的教学, 让学生感受到平时觉得很无用武之地的方法理论通过算法设计可以很好地用来解决实际问题, 使学生知道算法设计的重要性, 从而对这方面的学习产生浓厚的兴趣。

三、《数值计算》算法设计的实验教学

(一) 明确教学目的。

《数值计算》实验课程的授课目的应该是培养学生具备能够把实际问题抽象成为数学问题, 考虑采用数值计算方法解决问题的途径, 并能通过算法设计, 把相应的数值方法转化为计算机程序代码, 实现计算机自动运算得出结果。

(二) 突出教学内容。

《数值计算》实验课程是培养学生的算法设计能力, 实验课程的教学内容应该有相应的侧重。比如, 太理论化的定理、推论的证明过程可以简化为介绍性讲解, 而一些应用性较强的习题 (或小课题) 可以安排更多课时, 进行重点训练并着重讲解。目前的数值实验只是占据数值计算课程的小部分学时, 如果可能的话, 建议把它独立出来作为一门课程来讲授。

(三) 改善教学方法。

尝试把《数值计算》的教学与数学软件、数学建模、最优化理论等课程的教学相结合, 在讲授理论方法的同时也介绍实际问题的抽象化, 在算法介绍的过程中也可以通过实例讲解算法设计的过程。例如, 可以引用近年的数学建模竞赛题目进行讲解, 特别是数据量很大的题目, 重点说明数值计算方法的算法设计过程, 如果有条件, 可以在实验课程的后期, 安排一次数学建模设计竞赛, 使学生领会学有所用。

另一方面, 网络教学在教学改革中发挥着越来越大的作用。通过建设网络平台, 将课程录像、教学文件、课件、试题库等学生关心的信息以及便于教师间相互交流的信息放在网上, 为学生提供了强大的在线自主学习的环境。《数值计算》实验教学特别适合采用网络教学方式。

四、几种核心数值算法的实验教学设想

(一) 插值与拟合技术。

插值与拟合是从工程数据中寻找数学规律, 通过函数逼近或者数值逼近的方法找到离散点集的约束关系, 从而达到获取整体规律的目的。在教学过程中, 可以通过实例解析、实验引导、数学软件工具箱的操作演示, 帮助学生掌握插值与拟合的应用方法。

(二) 迭代逼近技术。

迭代逼近是一种典型的数值方法, 主要应用于方程求根、方程组求解和矩阵求特征值。其基本思想是逐次逼近, 先取一个粗糙的近似值 (即迭代初始值) , 然后利用同一个递推公式反复校正, 直至达到指定的精度要求为止。在教学过程中, 可以通过方程 (组) 校正举例, 重点说明如何赋予初始值, 如何进行逐次逼近求得方程 (组) 的解, 并要求迭代逼近方法满足逼近性和简单性。

(三) 简化压缩技术。

简化压缩技术是运用简单的计算方法解决复杂实际问题的过程, 把复杂的问题经过若干次压缩, 目的是通过“简单的重复”解决复杂问题, 减小解题的运算规模, 直到求出问题的解。简化压缩技术依据问题固有的特点, 采用循环、降维、去噪等简单的算法技术, 建立某种形式简单、思路明确、结构清晰的递归公式, 使复杂问题的计算过程简化, 层次得以压缩。在教学过程中, 可以通过层次解析、实例讲解、实验引导, 帮助学生理解简化压缩算法的精髓。

(四) 松弛调整技术。

松弛调整技术的核心是适应计算参数的可调性, 重复处理计算模型参数, 优化参数的调整、校正和自适应能力, 体现了把复杂问题简化处理的精髓。其中每一步松弛因子的选取, 会使用到数值迭代和加速收敛技术。教学过程中需要着重讲解如何选择合适的松弛技术和如何调整松弛因子, 使得学生能够比较容易地从理论学习过渡到实验训练。

五、结语

现代《数值计算》实验课程应该集中讲授如何实现算法设计, 采用以数值方法及其理论为基础, 以强化实验能力为引导, 执行理论和实验相结合的教学方法, 强调拟合技术、逼近技术、压缩技术、松弛技术的算法设计, 结合课堂讨论、案例教学、网络教学等途径, 进行数值计算实验教学。

参考文献

[1] .袁海燕, 安宇芳, 李敏静, 马艳芬.《数值分析》课程设计实践教学的几点探讨[J].学理论, 2013

[2] .黄陈思, 刘蓉, 王美清.浅谈数值实验在数值分析课程教学中的重要性[J].中国科教创新导刊, 2011

[3] .马淑芳, 张春蕊, 曲智林, 王文龙.数值分析课程教学改革研究与思考[J].林区教学, 2011

数值实验系统篇6

利用压电陶瓷材料逆压电效应制作的复合型薄膜压电振子 (压电陶瓷材料与金属基板复合) 已

被广泛应用到各个领域[1,2,3,4,5,6,7,8,9,10], 其常用结构有矩形和圆形两种, 矩型薄膜压电振子常简化为悬臂梁进行分析, 而圆形薄膜压电振子由于边界条件和结构的复杂性, 分析难度较大。对圆形薄膜压电振子, 目前多采用有限元分析法获得结构参数对变形量的影响关系, 如Christopher等[11]采用有限元法对固定、简支边界进行了静态驱动分析, 但分析中没有涉及驱动频率等方面的内容。Heyliger等[12]对压电圆片进行了自由振动分析, 该压电圆片结构中不包括非压电材料。阚君武等[13]利用瑞利法对圆形压电振子进行了弯曲振动分析, 从机电耦合系数最大化的角度分析了结构参数优化设计方法, 但没有考虑薄膜压电振子中黏结剂对其性能的影响。

为了获得弯曲变形薄膜压电振子的变形量, 本文从弹性薄板变形理论出发, 建立薄膜压电振子弯曲变形挠曲线方程, 得到压电陶瓷、黏结剂、金属基板材料及结构参数与振子变形量的关系, 利用MATLAB编程进行数值模拟计算和分析, 并与实验测试结果进行对比, 为提高薄膜压电振子工作性能提供了有效的方法。

1 圆形压电振子微元模型

圆形压电振子 (振子直径与厚度之比大于5) 常用黏结剂 (如环氧树脂等) 将轴向极化的压电陶瓷片和金属基板牢固地黏结在一起, 将金属基板周边固定支撑, 并给压电陶瓷施加交流信号, 则振子就会产生弯曲振动变形。在变形过程中认为各层之间不发生相对滑动且处于平面应力状态, 直法线假设成立, 因此可用薄板弯曲理论进行计算。

设压电陶瓷片半径为r1, 厚度为h1;黏结层 (胶层) 半径为r2 (r1=r2) , 厚度为h2;金属基板半径为r3, 厚度为h3。当压电振子驱动信号频率与其一阶模态频率一致时, 压电振子同心圆上各点位移相等[14], 因此, 取振子直径方向微元结构即可将其转化为梁结构, 设中性层位置为中性层与金属基板底面距离h, 压电陶瓷弯曲变形力矩为M1, 在金属基板上产生反作用力矩为M11, 由于金属基板周边固定, 故会在固定位置产生力矩M2, 振子弯曲变形微元模型如图1所示。

2 振子挠曲线方程

根据薄板弯曲理论, 中性层与金属基板底面距离[15]为

$\begin{array}{l} h = \frac{1}{2} \frac{1}{\frac{E_{3} h_{3}}{1 - γ_{3}^{2}} + \frac{E_{2} h_{2}}{1 - γ_{2}^{2}} + \frac{E_{1} h_{1}}{1 - γ_{1}^{2}}} \cdot \\ {\frac{E_{3} h_{3}^{2}}{1 - γ_{3}^{2}} + \frac{E_{2} [(h_{2} + h_{3})^{2} - h_{3}^{2}]}{1 - γ_{2}^{2}} + \\ \frac{E_{1} [(h_{1} + h_{2} + h_{3})^{2} - (h_{2} + h_{3})^{2}]}{1 - γ_{1}^{2}}} (1) \end{array}$

式中, 下标1、2、3分别表示压电陶瓷、胶层、金属基板;E为弹性模量;γ为泊松比。

压电振子直径比厚度大得多, 因此压电振子属于薄膜, 根据经典层合板理论[15], 压电振子各层具有相同的挠度, 即应变沿厚度方向呈线性分布, 应力呈非线性分布, 因此有

ε (z) =ε1=ε2=ε3=zk (2)

式中, k为应变斜率;z为厚度方向坐标。

设胶层与金属层交界面处应力为σi, 则有

$σ_{3} = z \frac{σ_{i}}{h_{3} - h}$

应变斜率为[15]

$k = \frac{σ_{i}}{h_{3} - h} \frac{1 - γ^{3}}{E_{3}} (3)$

胶层应力为

$σ_{2} = \frac{E_{2}}{1 - γ_{2}} ε_{2} = \frac{E_{2}}{1 - γ_{2}} z \frac{σ_{i}}{h_{3} - h} \frac{1 - γ_{3}}{E_{3}} (4)$

根据压电材料的d-型方程, 压电陶瓷在电压激励下产生应变为:ε′1=d31U/h1, 其中, d31为压电常数;U为驱动电压。则压电陶瓷应力为

$\begin{array}{l} σ_{1} = \frac{E_{1}}{1 - γ_{1}} (z k - ε^{'}_{1}) = \\ z \frac{E_{1}}{1 - γ_{1}} \frac{σ_{i}}{h_{3} - h} \frac{1 - γ_{3}}{E_{3}} - d_{31} \frac{E_{1}}{1 - γ_{1}} \frac{U}{h_{1}} (5) \end{array}$

设 $α = \frac{1 - γ_{3}}{1 - γ_{2}} \frac{E_{2}}{E_{3}}, β = \frac{1 - γ_{3}}{1 - γ_{1}} \frac{E_{1}}{E_{3}} ‚ η = \frac{E_{1}}{1 - γ_{1}}$ , 根据压电振子力矩平衡条件:

∫h3-h-hσ3zdz+∫h2+h3-hh3-hσ2zdz+∫h1+h2+h3-hh2+h3-hσ1zdz=0

可以求出:

$\begin{array}{l} σ_{i} = \frac{3}{2} η U d_{31} h_{1} (2 h_{3} - 2 h + 2 h_{2} + h_{1}) (h_{3} - h) \cdot \\ [h^{3} + (1 - α) (h_{3} - h)^{3} + (α - β) (h_{3} - h + h_{2})^{3} + \\ β (h_{3} - h + h_{2} + h_{1})^{3}]^{- 1} (6) \end{array}$

则

$\begin{array}{l} Μ_{1} = \int_{- h}^{h_{3} - h} σ_{3} z d z = \frac{1}{2} η U d_{31} h_{1} (2 h_{3} - 2 h + 2 h_{2} + h_{1}) \cdot \\ [(h_{3} - h)^{3} + h^{3}] [h^{3} + (1 - α) (h_{3} - h)^{3} + \\ (α - β) (h_{3} - h + h_{2})^{3} + β (h_{3} - h + h_{2} + h_{1})^{3}]^{- 1} (7) \end{array}$

由于压电振子属于叠层板结构, 且压电陶瓷和胶层直径要小于金属基板直径, 故根据经典叠层板变形理论[16], 需首先计算出压电陶瓷和胶层结构的等效弹性模量E′和等效泊松比γ′, 然后求出压电振子的等效弹性模量Ev和等效泊松比γv, 最后分别求出金属基板非压电陶瓷区域的挠曲线方程w3 (r) 和压电振子的挠曲线方程w (r) , 从而得到压电振子变形量计算方法。

压电陶瓷和胶层结构的等效弹性模量E′和等效泊松比γ′分别为

$\begin{array}{l} E^{'} = \frac{h_{1}}{h_{1} + h_{2}} E_{1} + \frac{h_{2}}{h_{1} + h_{2}} E_{2} + \\ \frac{h_{1} h_{2}}{h_{1} + h_{2}} \frac{E_{1} E_{2} (γ_{1} - γ_{2})^{2}}{h_{1} E_{1} (1 - γ_{2}^{2}) + h_{2} E_{2} (1 - γ_{1}^{2})} \\ γ^{'} = \frac{h_{1} E_{1} γ_{1} (1 - γ_{2}^{2}) + h_{2} E_{2} γ_{2} (1 - γ_{1}^{2})}{h_{1} E_{1} (1 - γ_{2}^{2}) + h_{2} E_{2} (1 - γ_{1}^{2})} \end{array}$

压电振子的等效弹性模量Ev和等效泊松比γv分别为

$\begin{array}{l} E_{v} = \frac{h_{3}}{h_{1} + h_{2} + h_{3}} E^{'} + \frac{h_{1} + h_{2}}{h_{1} + h_{2} + h_{3}} E_{3} + \\ \frac{h_{3} (h_{1} + h_{2})}{h_{1} + h_{2} + h_{3}} \frac{E^{'} E_{3} (γ_{3} - γ^{'})^{2}}{h_{3} E_{3} (1 - γ^{'}^{2}) + (h_{1} + h_{2}) E^{'} (1 - γ_{3}^{2})} \\ γ_{v} = \frac{h_{3} E_{3} γ_{3} (1 - γ^{'}^{2}) + (h_{1} + h_{2}) E^{'} γ^{'} (1 - γ_{3}^{2})}{h_{3} E_{3} (1 - γ^{'}^{2}) + (h_{1} + h_{2}) E^{'} (1 - γ_{3}^{2})} \end{array}$

金属基板弯曲刚度 $D_{3} = \frac{E_{3} h^{3}}{12 (1 - γ_{3}^{2})}$ , 压电振子弯曲刚度 $D_{v} = \frac{E_{v} h^{3}}{12 (1 - γ_{v}^{2})}$ , 压电振子变形分为两部分:一部分是压电陶瓷区域的变形w1, 另一部分为非压电陶瓷 (即金属基板部分) 区域变形w2, 压电振子挠曲线方程为

$w (r) = {\begin{cases} w_{1} (r) 0 \leq r \leq r_{1} \\ w_{2} (r) r_{1} \leq r \leq r_{3} \end{cases}$

$\begin{array}{l} w_{1} (r) = \frac{Μ_{11} r_{1}^{2}}{2 D_{3} [(1 - γ_{3}) r_{3}^{2} + (1 + γ_{3}) r_{1}^{2}}] \cdot \\ (- r_{1}^{2} + 2 r_{3}^{2} \lg \frac{r_{1}}{r_{3}} + r_{3}^{2}) + \frac{Μ_{1}}{2 D_{v} (1 + γ_{v})} (r_{1}^{2} - r^{2}) \\ 0 \leq r \leq r_{1} (8) \end{array}$

$\begin{array}{l} w_{2} (r) = \frac{Μ_{11} r_{1}^{2}}{2 D_{3} [(1 - γ_{3}) r_{3}^{2} + (1 + γ_{3}) r_{1}^{2}]} \cdot \\ (- r^{2} + 2 r_{3}^{2} \lg \frac{r}{r_{3}} + r_{3}^{2}) r_{1} \leq r \leq r_{3} (9) \end{array}$

根据连续性条件 $\frac{d w_{1} (r)}{d r} |_{r = r_{1}} = \frac{d w_{2} (r)}{d r} |_{r = r_{1}}$ 可得

$Μ_{11} = \frac{Μ_{1} E_{3} (1 - γ_{v}^{2}) [(1 - γ_{3}) r_{3}^{2} + (1 + γ_{3}) r_{1}^{2}]}{E_{v} (1 - γ_{3}^{2}) (1 + γ_{v}) (r_{3}^{2} - r_{1}^{2})} (10)$

将式 (10) 代入式 (8) 和式 (9) 得到压电振子弯曲变形挠曲线方程:

$\begin{array}{l} w_{1} (r) = \frac{Μ_{1}}{2 D_{v} (1 + γ_{v})} \frac{1}{r_{3}^{2} - r_{1}^{2}} \cdot \\ [r_{1}^{2} (r_{1}^{2} - 2 r_{3}^{2} \lg \frac{r_{1}}{r_{3}} - r_{3}^{2}) + (r_{3}^{2} - r_{1}^{2}) (r_{1}^{2} - r^{2})] \\ 0 \leq r \leq r_{1} (11) \end{array}$

$\begin{array}{l} w_{2} (r) = \frac{Μ_{1}}{2 D_{v} (1 + γ_{v})} \frac{r_{1}^{2}}{r_{3}^{2} - r_{1}^{2}} (r^{2} - \\ 2 r_{3}^{2} \lg \frac{r}{r_{3}} - r_{3}^{2}) r_{1} \leq r \leq r_{3} (12) \end{array}$

根据式 (11) 、式 (12) 可以计算出压电振子各点处的变形量, 该方程全面反映了压电陶瓷、金属基板和黏结剂材料参数、尺寸参数、驱动电压大小等与变形量之间的关系。

为了更加直观地表述各参数之间的影响关系, 下面通过数值计算对压电振子弯曲变形特性进行深入分析。

3 数值计算与分析

数值计算中压电振子所用材料及参数如表1所示, 利用MATLAB编程计算分析电压、金属基板材料、基板与压电陶瓷的直径、厚度以及胶层厚度等结构参数对压电振子变形量的影响, 以确定压电振子的最优参数。

3.1 驱动电压与振子变形量的关系

从式 (7) 、式 (11) 、式 (12) 可以看出, 其他参数确定不变的情况下, 输入电压越大, 压电振子变形量越大。图2是压电振子基板直径为35mm、厚度为0.3mm, 陶瓷直径为25mm、厚度为0.3mm, 不同电压下压电振子直径方向上各点最大变形量曲线, 从图2中可以看出, 压电振子同一半径圆上各点的变形量相同, 这与圆形薄膜压电振子一阶模态弯曲变形规律相符。当压电振子施加电压达到200V时, 中心点变形量可达到49μm, 电压为100V时, 中心变形量可达到24μm。图2中反映的是各点最大变形量, 计算中取电压最大值, 没有考虑频率的影响, 因此图2反映的是压电振子的静态变形量随电压的变化, 电压越大, 变形量越大。

3.2 金属基板材料与振子变形量的关系

金属基板是压电振子的重要构成部分, 选择不同的基板材料, 振子工作性能不同。图3反映了相同条件下 (输入电压为200V, 基板直径为35mm、厚度为0.3mm, 陶瓷直径为25mm、厚度为0.3mm) , 不同基板材料对压电振子变形量的影响。从图3中可以看出, 分别选用铜和钛合金作为基板材料的压电振子变形量接近, 钛合金生物兼容性好, 但成本高;铝合金基板的压电振子变形能力较弱。因此, 为获得较大的变形量, 应选择弹性模量较大的金属材料作为基板材料。

3.3 金属基板与压电陶瓷尺寸对振子变形量的影响

压电振子重要的结构参数是陶瓷片、基板直径以及厚度, 为了便于定量分析, 计算时, 铜基板直径为35mm、厚度为0.3mm。图4反映了压电陶瓷片与基板直径比对变形量的影响, 基板直径一定的情况下, 陶瓷片直径越大, 振子变形量越大;由于压电陶瓷硬且脆, 当陶瓷与基板直径相同时, 固定振子边缘会使压电陶瓷损坏而不能正常工作, 考虑到振子安装时固定金属基板边缘, 因此直径比选择在0.8～0.9之间在制作工艺上是可行的。

图5反映了陶瓷片与基板厚度比以及胶层与基板厚度比对振子变形量的影响, 计算中基板厚度选为0.3mm。从图5中可以看出, 陶瓷与基板厚度比为0.8时, 变形量达到最大;胶层黏结强度足够的情况下, 胶层厚度越小, 振子变形量越大。但在实际制作中, 根据制作要求常选陶瓷与基板厚度比在0.6～1.1之间, 黏结层厚度控制在0.05mm以内。文献[17]对黏结工艺及黏结层对振子阻抗的影响进行了深入研究, 为复合压电振子制作提供了工艺方法。

4 实验结果分析

为了验证上述分析结论, 定做了表2中不同尺寸的圆形压电振子 (制作工艺保证胶层厚度相同, 均为0.05mm) , 利用多普勒激光测振仪 (LV1610) 分别进行振幅测试, 测试时在基板上标定出对称位置的测试点, 图6反映了驱动电压为150V、频率为120Hz时4个样品的变形量曲线及样品1、3的计算曲线, 各样品的测试结果与数值模拟规律相符, 如样品1的陶瓷与基板直径比为0.75、厚度比为1, 样品3的陶瓷与基板直径比为0.71、厚度比为0.66;根据数值分析结论, 样品1的测量值比样品3的测量值大。在数值分析中未考虑频率变化对变形量的影响, 因此认为计算值为其静态变形量, 测试中施加直流电信号驱动。初步测试压电振子在各自谐振频率点振动时, 中心点变形量仅为其静态变形量的5%～10%, 而在80～150Hz频率段 (准静态频率) , 中心点变形量为静态变形量的50%～60%, 这主要受压电振子幅频特性及机电耦合特性的影响[18], 从图6中可以看出, 振子1、振子3的准静态变形量是其静态变形量的55%。图7反映了1号样品在准静态频率段振幅特性, 驱动频率在80～130Hz时, 振幅较大。

5 结论

利用弹性薄板变形理论推导了复合薄膜压电振子一阶弯曲变形挠曲线方程, 真实反映压电陶瓷、黏结剂、金属基板材料与结构参数对变形量的影响;通过编程计算, 得到驱动电压、基板材料、振子结构参数 (陶瓷与基板的直径比、厚度比) 等对变形量的影响;并通过实验测试加以验证。为了提高压电振子变形量, 可采用以下方法:

(1) 适当提高驱动电压, 采用弹性较好的金属材料作为基板;

(2) 在结构允许的情况下, 尽可能增大陶瓷与基板直径比;

(3) 陶瓷与基板厚度比为0.8时, 振子变形量最大, 实际制作中, 厚度比可选在0.6～1.1之间, 保证黏结强度的情况下, 胶层厚度越小, 对提高振子变形量越有利;

(4) 压电振子中心点的谐振振幅仅为静态变形量的5%～10%, 为了获得较大的振动变形, 薄膜压电振子应工作在准静态频率段, 通过测试振子的准静态幅频特性, 确定振子最佳工作频率。

数值实验系统篇7

1 瓦斯爆炸的实验研究

1.1 实验系统

瓦斯爆炸的实验系统如图1所示。系统主要由爆炸试验管道、配气系统、点火控制系统、压力采集系统、高速摄影仪系统等五个部分组成。爆炸实验采用内径100mm, 有效管长1 500mm的石英玻璃管道, 整条管道两端密封, 顶部设置泄压口, 管道一端设置点火电极, 利用可调高能点火器点燃预混气体, 如图2所示。可调高能点火器的型号KTD-A, 输入电压为220V, 点火频率3~20次/s, 如图3所示。实验压力数据通过SZSC数据采集器采集, 采集频率为10 000ftp/s。

1.2 实验条件和实验步骤

实验室温度在15~25℃, 相对湿度为40%~70%, 气压为大气压。由于压力采集装置对采光有一定要求, 因此实验均在天气晴朗、光线充足的条件下进行。实验分别在距离点火端500、1 000mm处设置压力监测, 爆炸发生时, 记录下该处的爆炸压力。实验操作过程如下:

(1) 用真空泵将实验管道抽成真空状态, 充入实验所需体积数的甲烷气体, 打开进气阀, 充入空气;

(2) 打开循环泵, 使甲烷和空气充分混合3min;

(3) 检查其他阀门是否关闭, 打开泄压阀;

(4) 按下点火器开关, 点燃甲烷预混气体, 实现爆炸;

(5) 压力采集系统记录下压力数据。

1.3 实验结果分析

1.3.1 甲烷爆炸的火焰传播的速度结果分析

李孝斌提出了基于MATLAB图像处理技术对实验照片进行分析, 拍摄照片与实际大小之间有一个比例系数, 据此可以得出甲烷爆炸实际可见光前沿传播位移和速率。计算如式 (1) 、式 (2) 所示。

式中:L前沿为可见光前沿的位移;d拍摄为拍摄有效管长;N为照片沿管长像素数;n前沿为可见光前沿位置像素值;n起始为可见光起始位置像素值;v前沿为可见光前沿的速率, m/s;v拍摄为拍摄速率, 帧/s;L起始i-1、L前沿i-1分别为第i和i-1计时时刻可见光前沿位移, m。

表1为甲烷反应可见光前沿的位移和速率分析。可以看出, 气体被点燃后, 浓度相对集中, 火焰前沿被拉伸, 使火焰加速传播到最大值128m/s, 相对应的压力也快速上升。随着泄爆膜破裂和燃烧物耗尽, 速率越来越小, 直至火焰熄灭。

1.3.2 甲烷爆炸的压力结果分析

图4为体积分数9.5%, 距离点火端500、1 000mm处2测点的压力采集曲线。可以看出, 在0~3s内, 500mm处的爆炸压力始终大于1 000mm处的爆炸压力, 说明距离点火端越近, 压力越大;500mm处的爆炸压力出现两个压力峰值, 分别在为7 000、6 000Pa。由于压缩波的存在, 爆炸波重复叠加造成震荡, 产生第一个高压。又因横向压缩波的存在或当泄压膜破裂之后释放出的未反应的预混气的燃烧, 出现第二个峰值。而1 000mm处会出现负值, 说明爆炸的瞬间, 压力全部释放。由于管道内剩余混气的燃烧压力上升到正值, 在这个阶段, 压力曲线出现一个峰值;3s以后, 两条曲线都趋于稳定。

2 瓦斯爆炸的数值模拟研究

2.1 网格划分及初始条件的设置

考虑到模拟三维模型时计算量庞大, 管内流体的模拟往往需要简化模型。且爆炸实验管道为轴对称的圆柱形, 对于轴对称元件, 将三维模型转化为二维模型是最为简便而常用的方法。因此, 分别将长1 500mm, 横截面积为100、1 000mm2的三维管道, 简化为长1 500mm, 宽为100mm的矩形管道, 忽略管道厚度, 利用流体动力学软件Fluent对体积分数为9.5%的甲烷/空气预混气体爆炸进行模拟。

2.1.1 边界及流场初始条件

模型的四个面边界类型均设为绝热壁面, 温度为常温300K;

对于二维封闭管道均匀混合瓦斯气体爆炸, 初始时刻为t0, T (t0) =300 K, P (t0) =1.01×105, u (t0) =0, v (t0) =0;

甲烷体积分数为9.5%时, 各组分的质量分数分别为:

2.1.2 点火初始条件

对于密闭管道内均匀瓦斯预混气体或局部瓦斯气体爆炸, 假设管道封闭端面的瓦斯预混气体遇到点火源发生爆炸。点火源产生能量加热火焰附近的局部瓦斯预混气体, 使其温度迅速升高达到着火温度而点燃;然后瓦斯气体借助火焰传播使整个瓦斯/空气预混气体着火燃烧。根据热点火理论, 在模拟计算中可将高温已燃气团设置为点火区, 该区域的温度和组分浓度为:T (t0) =2 000K,

2.1.3 网格划分

通过采用5、10、50、100 mm划分网格模拟, 得出其他3种网格模拟时对计算收敛性影响很大, 并且会降低计算的精度, 采用10mm划分一个网格能够达到理想的模拟效果, 网格总数为1 500。

2.2 模拟结果分析

2.2.1 小尺寸管道的模拟分析

(1) 甲烷爆炸的速率模拟结果。图5为甲烷爆炸的速率模拟和实验结果对比图。模拟结果显示, 速率随着时间先增大后减小, 模拟最大值达到148m/s, 实验最大值为128m/s。实验曲线从最大值开始下降, 且低于模拟值。最后两条曲线趋于稳定。

(2) 甲烷爆炸的压力模拟结果。图6为体积分数为9.5%, 距离点火端500、1 000mm处2测点的压力采集曲线, 模拟结果与实验相差不大。500mm处压力峰值为7 500、5 800Pa。1 000mm处首先出现负值, 然后压力上升趋于稳定。

通过采用实验方法与模拟方法对比的方法, 得出了实验与模拟的速率和压力参数的变化趋势基本一致。因此, 通过模拟方法模拟大尺寸管道内的压力特征及火焰发展规律, 为进一步研究实际大型巷道内爆炸的实验研究提供支撑。

2.2.2 大尺寸管道的模拟分析

其他条件不变, 将小尺寸管道放大数倍, 模拟长为15m, 内径为2m的管道内甲烷爆炸的火焰传播规律和压力变化。

(1) 甲烷爆炸的压力模拟结果。图7为体积分数9.5%, 距点火端200、500、1 000cm处3个测点的压力采集曲线。由图7可以看出, 在0~4s内, 离点火端最近的2m处的压力始终最大;从开始点火后不到1s的时间内发生爆炸, 压力急剧上升, 最大压力为8 000Pa。变化趋势与小尺寸管道基本一致。

(2) 甲烷爆炸的火焰传播及发展过程模拟。图8为甲烷气体从点火到爆炸火焰传播过程图。图9为对应的速度矢量变化图。

由图8、图9可以看出, 甲烷气体被点燃后, 首先形成球面火焰 (0.2s) , 火焰在管道中呈半球形沿各个方向传播。随着火焰的传播, 由于管壁的作用, 火焰前沿轴向伸展, 火焰被拉伸, 导致火焰的加速和压力的上升。由于在管壁处形成漩涡, 形成反向流, 然后反向流变得越来越强, 从而延缓了火焰的传播。在反转流场的拖拽和漩涡的推挤作用下, 火焰前沿在靠近管壁处向前传播比中心处快, 最后形成一个与小尺寸管道中相似的“郁金香”火焰, 但这种现象没有小尺寸管道中的明显。随着火焰的发展, 漩涡向后移动并逐渐衰减。随后火焰变成一个平面向前传播, 因为平面火焰没有凸面火焰的优势, 火焰的接触面积变大, 速度减小, 在火焰的边缘反应物质得不到补充, 当火焰的边缘接触到管壁时, 逐渐被管壁所冷却。在火焰开始熄灭前, 火焰的传播速度再次下降, 靠管壁处的火焰与管壁发生分离, 直至火焰熄灭。

4 结论

(1) 实验条件和模拟条件下的压力变化和速度变化趋势一致。随着时间压力先增大减小, 随后有一个小的浮动增加, 再次下降, 最大值为7 000Pa。速率曲线先上升后减小, 随着压力的增加对应着速度也增加到最大值148m/s。

(2) 压力曲线出现两个峰值, 一个是泄爆压力, 一个为剩余的混气的燃烧产生的压力。

(3) 大尺寸管道压力曲线与小尺寸的基本一致, 但大尺寸的爆炸压力大于小尺寸的爆炸压力。

(4) 由于漩涡的存在, 形成反向流, 因此在大尺寸的甲烷爆炸的火焰传播特征出现类似小尺寸管道的“郁金香”火焰。并且平面火焰的接触面积变大, 使得燃烧物质得不到补充, 促使燃烧熄灭。

参考文献

[1]Ibrahim S S, Masri A R.The effect of obstructions on overpressure resulting from premixed flame deflagration[J].Journal of Loss in the Process Industries, 2011, (3) :213-221.

[2]Bi Ming-shu.Numerical simulation of premixed methane-air deflagration in large L/D closed pipes[C].Applied Thermal Engineering, 2012.

[3]Zhu Chuan-jie.Numerical simulation of blast wave oscillation effects on a premixed methane/air explosion in closed-end ducts[J].Journal of Loss Prevention in the Process Industries, 2013:851-861.

[4]胡铁柱.瓦斯爆炸传播规律数值模拟研究[D].北京:中国矿业大学, 2005.

[5]李孝斌.矿井瓦斯爆炸感应期内光学特征研究[D].西安:西安科技大学, 2006.

[6]李孝斌.半封闭空间甲烷爆炸初期火焰传播时空分布研究[J].Natural science, 2011, 30 (6) :20-24.

[7]黄文祥, 李树刚, 李孝斌, 等.不同点火能量作用下管道内瓦斯爆炸火焰传播特征[J].煤矿安全, 2011, 42 (8) :7-10.

[8]VitalyBychkov.Flame acceleration in the early stages of burning in tubes[J].Combustion and Flame, 2007, 150 (4) :263-276.

[9]肖华华.管道中氢-空气预混火焰传播动力学试验及数值模拟研究[D].北京:中国科学技术大学, 2013.

[10]Jozef J arosinski.Combustion phenomena[M].US:CRC press, 2009.

数值实验系统篇8

关键词：线圈,导热系数,数值模拟

1 前言

线圈径向导热系数的测定,对于研究线圈的浇注、固化等过程有着重要的指导作用。基于稳态传热法基本原理[1],采用自设实验装置对模拟线圈径向导热系数进行了测定,获得了线圈内外壁和内部盛水的温度分布情况,按照热力学平衡方程推导线圈的导热系数。为了模拟和多次重复实验过程,采用国际著名的三维计算流体动力学软件FLOW-3D建立了全尺寸线圈基于控制体积的有限差分计算模型[2,3],考虑导热、对流和辐射三种传热方式以及密度扩散方程,计算了导热系数分别为λ=0.25W/(m·K)、λ=0.3W/(m·K)、λ=0.35W/(m·K)和λ=0.4W/(m·K)四种工况,与实验测得的模拟线圈内外壁温度分布进行对比,结果表明λ=0.3W/(m·K)这种工况更接近改进的实验结果。

2 实验测试

2.1 实验备件

准260mm×365mm模拟线圈(不含导线,见图1),加热棒,固态继电器,热电偶,测温仪,笔记本电脑,电度表,空气开关,隔热塑料桶等。

2.2 实验步骤

(1)用玻璃胶将脱模纸平整地粘贴在薄木板上,然后将线圈底部与粘有脱模纸的薄木板密封,以防桶内盛装的水泄露。

(2)在线圈同一环面内外壁相对应的位置粘贴热电偶,具体方案如图2所示。

(3)将60℃的温水倒入线圈中,盖上密封盖,防止热量从顶部散失。同时开启测温仪,每隔10s采集一次扫描点的温度。

(4)图2中的外表面测温点:A1、A2、A3、A4和B1、B2、B3、B4和C1、C2、C3、C4,对应的内表面测温点:D1、D2、D3、D4和E1、E2、E3、E4和F1、F2、F3、F4。

2.3 实验结果及分析

图3出示了同一环面外表面B1、B2、B3、B4和内部水温的实验测试曲线,图4出示了同一环面内表面E1、E2、E3、E4和内部水温的实验测试曲线,图5出示了线圈内水温在不同区域的分布和变化情况,图6出示了外表面B3和内表面E3的实验曲线。

图3可以看出,实验开始持续4000s,线圈外表面处于瞬态传热过程,温度在升高,随后达到稳态,温度缓慢下降并与水温下降保持同步,窗口凸台B4处温度较其它部位温度偏低,主要是由于凸台处较厚,热量传导慢所致。

从图4可以看出,内表面温度与水温有一定差别,主要是由于两种不同介质在界面存在热阻导致。

从图5可看出,上部水温比底部水温偏高约6℃,这主要是由于水的对流扩散引起,即水的温度越高,密度越小,因此线圈内水的温度从底部到顶部逐渐升高,存在明显的温度分层现象。

图6显示了线圈同一高度内外表面温度变化情况,实验开始持续4000s,线圈传热属瞬态过程,随后则达到了稳态(系统处于热平衡状态)。

径向导热系数λ的计算,从图6取4000s到6000s为一个稳态传热单元,线圈内表面E3点在4000s时的温度是55.5℃,在6000s时的温度是54.5℃,温度差△t=1℃,而线圈外表面B3点的温度是41℃。可近似认为E3点在4000s到6000s之间水所损耗的热量约等于线圈从E3到B3传导热量。其中线圈内径d1=210mm,线圈外径d2=260mm,高度h=365mm,t1=54.5℃,t2=41℃,水密度ρ=1.0×1.03kg/m3,水比热c=4.2×103kJ/(kg·K),一个稳态传热单元的时间t=2000s。

水损耗的热量:

线圈传导的热量:

同理,以6000s到8000s为一个稳态传热单元,λ=0.26W/(m·K);以8000s到10000s为一个稳态传热单元,λ=0.28W/(m·K);以10000s到12000s为一个稳态传热单元,λ=0.31W/(m·K)。

实验误差主要来自忽略了线圈外表面辐射换热和环境的对流换热,由于实验中使用了隔热塑料桶,保证线圈周围环境温度均匀且稳定,因此由线圈与环境的对流换热引起的误差很小。另外,水温的不均匀性也是导致误差的主要原因之一,为此改进了实验方法,重新进行了测试。

2.4 改进的实验

为了避免线圈内部水温分层现象,在前面实验将水倒入线圈中,同时加入少量食盐,利用重力使盐的浓度上低下高,因而比重也是上低下高,由此可以抑制因顶部散热而导致的内部上下环流,使温度趋于一致并更接近模型的假设。其它步骤完全相同。

图7出示了同一环面外表面A1、A2、A3、A4和内部水温的实验测试曲线,图8给出了同一环面内表面D1、D2、D3、D4和内部水温的实验测试曲线,图9给出了同一环面外表面A1和内表面D1的实验曲线。

根据图9推算导热系数λ,取5000s到10000s为一个稳态传热单元,线圈内表面D1点在5000s时的温度是52.8℃,在10000s时的温度是49.5℃,温度差△t=3.3℃,而线圈外表面A1点的温度是38.2℃。可近似认为D1点在5000s到10000s之间水所损耗的热量约等于A1点导热所消耗的热量,计算过程同上。

以5000s到10000s为一个稳态传热单元得出:λ=0.31W/(m·K);以10000s到15000s为一个稳态传热单元,λ=0.33W/(m·K);以15000s到20000s为一个稳态传热单元,λ=0.28W/(m·K)。

实验误差主要来自忽略了模拟线圈外表面辐射换热和环境的对流换热,由于实验环境温度基本保持恒定,因此这部分的误差很小,实验结果可靠。

3 模拟线圈导热系数的数值分析

3.1 CFD模型建立

采用Solidwork软件建立试验线圈的三维几何模型[4],如图1所示。采用FAVOR技术对其进行网格剖分,形成流体区域和固体区域,如图10所示。按照试验要求在线圈内部盛满初始温度为60℃的温水。选择传导、对流和辐射换热模型,牛顿流体、密度依赖温度变化模型。线圈初始温度设为27℃,线圈和水的热物性参数见表1。

分别对λ=0.25W/(m·K)、λ=0.3W/(m·K)、λ=0.35W(m·K)和λ=0.4W/(m·K)四种工况进行了模拟,获得线圈和水内部温度分布情况以及导热系数对其温度分布的影响。

3.2 计算模型的验证

为了以实验验证λ=0.31W/(m·K)线圈CFD模型的正确性,图11给出了同一环面外表面和水温度模拟曲线,与图7相比,无论趋势还是具体数值基本吻合。图12给出了同一环面内表面温度模拟曲线,与图8相比,无论趋势还是具体数值基本吻合。图13出示了同一环面内外表面温度模拟曲线,与图9相比,无论趋势还是具体数值基本吻合。图14给出了内部水温分布模拟结果,也存在温度分层现象。对比结果表明,线圈模拟实验的计算模型合理,结果可靠。

4 结论

通过线圈导热系数测定实验和数值模拟,得到以下结论:(1)模拟线圈(不含导线)的导热系数约为0.3W/(m·K)。(2)水中加入少量食盐可改善线圈内盛水的温度分布,有效地减弱了水温分层现象。

参考文献

[1]戴锅生.传热学(第2版)[M].北京:高等教育出版社,1999.

[2]曾攀.有限元方法(第5版)[M].北京:清华大学出版社,2008.

[3]庄茁,岑松.有限元方法(第5版)[M].北京:清华大学出版社,2006.

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

【数值实验系统】相关文章：

数值分析实验报告07-23

数值模拟实验报告08-14

河工大数值分析实验04-16

气膜冷却流场的实验研究和数值模拟的分析04-25