动力学计算

动力学计算（共12篇）

动力学计算 篇1

0 引言

高速飞行器在运动过程中受外界摩擦阻力和周围的温度变化影响极大, 所以在飞行器的结构设计时应将这些因素考虑进去, 热流固耦合分析是解决这种问题的关键。本文对飞行器高速运动时进行热流固耦合分析, 热分析是指对飞行器周围的温度场进行分析, 流分析是指对飞行器周围的流动场进行分析, 固分析是指对飞行器的结构场进行分析, 热流固耦合分析是量化不同物理场之间相互的影响。温度场和流体场需要针对飞行器周围流体场进行分析, 需要进行流体动力学计算, 计算的结果导入结构场中进行耦合分析。计算流体动力学需对飞行器周围流体场划分网格, 计算的结果储存在节点内, 而分析结构场中需对飞行器表面结构场划分网格, 导致在数据传递过程中, 不同物理场的网格节点无法统一。前人阐述对热流固耦合分析数值传递过程中, 只描述对两种物理场信息求解软件之间接口问题, 局限于两种物理场网格划分相对一致的近似计算情况, 精度有限。本文利用搜索函数和权函数思想, 对流体动力学计算出的数据进行处理, 将数据传递给结构场中进行分析, 给出处理此传递过程的数据传递公式, 精确顺序热流固耦合场仿真问题模拟。

1 流体动力学仿真计算结果传递

通过流体动力学计算得到某一时间点的数据信息。在飞行器运动过程, 利用计算流体动力学得到某飞行器周围流场网格节点内的数据信息, 包括节点压力与节点温度值, 并对飞行器表面的结构场进行网格划分。流体场划分网格节点数量比结构场中使用单元节点多, 所以无法直接运用流体动力学数值进行耦合分析。如图1所示为某飞行器外表面和周围场网格划分参考示意图, 实线表示飞行器外面表结构场中单元划分情况, 虚线表示流体场中单元划分情况;a、b、c、d为结构场分析使用单元节点, I、J、K、L为流体动力学分析使用单元节点;节点围成四边形为一个单元体。

用数学方程表示计算流体动力学得到的飞行器周围流场外节点载荷信息:

式中:xi、yi与zi为第i个流场单元节点坐标值;Pi是第i个流场单元节点压强值;Ti为是第i个流场节点温度值;n为此流场单元节点个数。

将求得的信息存储在分析过程中, 并记录在数组。建立如下函数, 并记录结构场分析飞行器外表面节点的载荷信息:

式中:pj为第j个结构场节点的压强值;Tj为第j个结构场节点的温度值;m为结构场单元节点个数。

设定搜索函数, 此搜索函数可记录节点内对应信息:

式中:ej为飞行器结构外表面标注的单元号;mI、mJ、mK、mL为对应周围流场的网格节点号。

由于存在对应关系, 流场网格比结构场计算网格稠密, 因此存在不同场网格的对应关系。该函数可找到各自结构场单元节点对应的最近流体场节点, 并记录于数组, 利用最近的4个流体场节点储存的数值可近似估计结构场节点对应的数值。

用几何思想, 求出4个流体场节点围成的几何形心, 并把此几何形心的数值赋予结构场ej上, 方程如下:

式中:wi为结构场单元4个节点数值的权重值;Φ (nj) 是构成流体场网格4个节点各自压强载荷值;Ψ (nj) 为流体场网格4个节点各自对应温度载荷值。

把各自流体场节点对结构场节点的空间距离当作衡量权重标准, 规定方程:

最后得到结构场单元的单元面载荷数组:

式中:lij (i=a、b、c、d;j=I、J、K、L) 为两个节点间空间距离;ej为第j个结构场单元的载荷信息;Pj为第j个结构场的压强值;Tj为第j个结构场的温度值。

2 流体动力学数据传递在飞行器运动中热流固耦合分析实例

在飞行器飞行中, 利用FLUENT和ANSYS软件分析热流固耦合场。在FLUENT软件得到流体场的温度和压强值, 温度值和压强值在ANSYS结构场中进行耦合分析。利用前两章节公式处理数据, 便可以将FLUENT计算的结果仿真出飞行器外表面的热载荷和力载荷, 这两种载荷同时施加在结构场单元上。如图2与图3所示, 为施加到飞行器的力载荷和热载荷信息的转换。

3 结论

物体在高速运动中, 需考虑热和力同时对结构体的影响, 单纯地把流体动力学网格叠加到结构场的网格这种方法适用于计算要求误差不高的情况, 而对数据处理, 再导入结构场分析, 这种方法使热流固耦合分析过程更加精确。本文完整热流固耦合数据处理工过程, 建立了FLUENT与ANSYS联合仿真平台, 填补了热流固耦合分析数据传递这一领域的空白。

摘要：某高速飞行器在超音速飞行中, 周围温度和压力会对其结构产生影响, 热流固耦合分析可以量化这一影响。由于热流固耦合仿真对象为多物理场, 所以流体动力学计算出的物理量导入结构场中需要进行数据传递。文中给出多物理场之间数据传递的推导过程和推导方法, 最终实现对流体动力学计算出的数据在结构场中仿真模拟, 量化温度场和压力场对结构场的影响。

关键词：流体动力学,热流固耦合,数据传递

参考文献

[1]赵善彪, 张天孝, 李晓钟.红外导引头综述[J].飞航导弹, 2006 (8) :42-45.

[2]孟庆超, 段萌, 张运强, 等.红外空空导弹整流罩技术的新进展[J].航空兵器, 2008 (2) :24-27.

[3]臧国才, 李树常.弹箭空气动力学[M].北京:兵器工业出版社, 1989:80-81.

[4]左克罗M J, 霍夫曼J D.气体动力学[M].北京:国防工业出版社, 1984:39-45.

[5]马铁犹.计算流体力学[M].北京:北京航空学院出版社, 1989:80-82.

[6]王淑华.高速飞行器头罩及侧窗热流计算[M].北京:航天科技集团701所, 1997:43-48.

[7]纪楚群.导弹空气动力学[M].北京:宇航出版社, 1996:11-133.

[8]GILBERT K G.Over view of aero-optics[G]//Progress in Astronautics and Aeronautics Series, 80, 1982:1-9.

[9]何友金, 张鹏, 彭军, 等.高速红外空空导弹整流罩气动加热及应力分析[J].红外技术, 2007, 29 (7) :373-376.

动力学计算 篇2

2.如图12-32所示，在铅垂平面内有重为G，长为ι的均质直杆OA，由油缸推力（或拉力）驱使绕O点摆动，G = 300N，ι= 1.2m。在图示位置，已知直杆OA的角加速度4rad/s2,求液压缸的拉力。

3.如图11-57所示，柴油机的曲柄OAr80cm，连杆ABL160cm，转速n116r/min。求当分别为0o和90o时，连杆AB的角速度和活塞B的速度。

4.图10-14所示重物G=5kN，从H=40mm处自由落下。已知直杆长l2m，截面积 A = 30×30mm2，弹性模量 E=200GPa。求冲击时杆件内的最大正应力。5.试校核图9-16所示千斤顶丝杠的稳定性。已知其最大承载Fp=150kN，有效直径d152mm，长度l0.5m，材料为Q235钢，s235MPa，稳定安全因数[n]st1.8。丝杠的下端可视为固定端约束，上端可视为自由端。

6.如图8-13所示，转轴传递的功率P = 2 kW，转速n100r/min，带轮直径D=250mm，带的拉力FT2Ft，轴材料的许用应力[]80MPa，轴的直径d=45mm。试按第三强度理论校核轴的强度。

7.空心管梁受载如图7-36所示，已知[]150MPa，管外径D=60mm。在保证安全的条件下，求内径d的最大值。

8.空心钢轴的外径D = 100mm，内径d = 50mm，要求轴在2.4 m内的最大扭转角不超过1.8o，材料G=80GPa。求：（1）轴内最大切应力。（2）当转速n100r/min时轴所能传递的功率。9.如图5-13所示，设钢板与铆钉的材料相同，许用拉应力[]160MPa，许用切应力[]100MPa，许用挤压应力[jy]320MPa，钢板的厚度t =10mm，宽度b = 90 mm，铆钉直径d = 18mm，拉力F = 80 kN，试校核该联接的强度（假设各铆钉受力相同）。

10.一钢制阶梯杆如图4-38所示。已知AD段横截面面积为AAD400mm2，DB段的横截面面积为ADB=250 mm2，材料的弹性模量E=200GPa。试求：（1）各段杆的纵向变形。（2）杆的总变形lAB。（3）杆内的最大纵向线应变。

11.如图3-24所示，一转轴上装有齿轮与带轮。已知皮带拉力Ft11300N,Ft2700N，皮带轮半径R=50cm，齿轮圆周力Ft1000N，径向力Fr364N，齿轮节圆半径r = 30cm，a0.5m，试求轴承A、B处的约束力。

《计算结构力学》教学改革 篇3

关键词：计算结构力学 实践教学 科研项目

中图分类号：G642文献标识码：A文章编号：1673-9795（2012）10（b）-0015-01

1 《计算结构力学》课程建设的意义

《计算结构力学》在各个特别是开设土木工程专业的高校中都受到重视。在《计算结构力学》课程的教学中，往往存在一种重复机械记忆学习的指导观念，即学生被动的接受知识，学生的主动性与能动性被忽略。这种机械的教与学的方式导致学生的独立性和能动性降低，从而严重阻碍学生个性与能动性的发挥。另一方面，随着社会经济的发展以及工程内容日益复杂，土木工程专业研究生的教育和教学都必须进行相应的改革，以满足土木工程专业对学生的更高要求，特别是《计算结构力学》课程的教学中提出了更高的标准。要求学生在了解常规结构分析方法的基础之上，能够进行复杂的结构分析计算；并且能够通过计算结果分析总结结构分析方法的一般规律。在工程结构日益复杂和计算精度要求不断提高的前提下，力学分析中必须借助计算机等高效辅助工具才能到达这一标准。

目前土木工程学科研究生教育还存在着培养出来的学生难以满足企事业单位的需求的问题。

主要表现为：

（1）理论与实际相脱离，对所需解决的工程项目不知如何下手，特别是动手解决实际工程问题的能力差，墨守理论，缺乏创新能力。

（2）缺乏科学思维和科学态度，愿意相信书本，缺乏足够的独立思考和决断能力。

（3）综合能力薄弱，知识面狭窄，难以成为高质量的复合型人才，不适应现代化经济建设和工程技术发展的要求。

2 《计算结构力学》教学改革

2.1 科研项目于教学相结合

研究生培养是培养青年高层次人才，需要做到“科学研究与高等教育的有机结合”。但是目前的状况是高等教育普遍存在理论和实践脱钩，不能将所学的知识应用到实际工程中，一流的学术并不能代表一流的人才培养质量。现在国内很多高水平教师忙于搞科研、做学术，但他们没有把那些高深的学术理论、先进的学科项目和自己所教授的课程结合起来。由于国内教育体制和管理制度的原因，大学生无法接触到大学教授们进行的那些学术活动，使这些最优质的教学内容失去了其宝贵的价值。为了改变这种状况，我们在《计算结构力学》的教学中，为了培养学生科学分析解决问题的能力，直接针对一些实际做过的课题，作为教学案例加入到课程内容当中，对学生进行讲解。例如某体外预应力钢-轻骨料混凝土结构研究，目的是建立一套完整的用于指导预应力钢-轻骨料混凝土组合梁设计承载能力的计算公式。对此课题研究方法和计算结构力学方法应用的讲解将提高学生分析解决问题的能力。

2.2 理论与实践相结合

当前高等教育存在大学与社会需求脱节，造成就业市场上同时出现“人才荒”和“就业难”。一方面，伴随着我国经济社会从粗犷型向节约型的转变以及经济结构的优化调整，新兴专业的高素质、高专业技术人才严重不足；另一方面，大学里一些专业设置和课程内容陈旧过时，没有及时跟上“转方式，调结构”的人才要求。我课题组在《计算结构力学》课程教学的基础上，组织学生参与一些正在进行的实际工程项目，例如课题组曾组织研究生参与“某市快速轨道交通轻轨双线（22+55+22）m预应力混凝土拱-连续组合箱梁吊杆张拉施工监控”的工程实际项目，此工程的目的是对此桥吊杆张拉阶段和吊杆索力调整阶段的施工监控，以确保该阶段施工过程中各主要受力构件不发生损害，保证施工顺利进行。学生通过这样实际工程的参与，使得研究生在掌握《计算结构力学》课程内容的基础上，提高了用《计算结构力学》分析解决实际工程问题的能力。

3 软件应用

计算力学已经引发一个令人振奋的新观点：理论、实验和计算成为现代科学的三大支撑；产生了一个新的领域“计算科学”。计算力学的延伸造就了CAE软件和产业，而CAE产业产生了巨大的社会和经济效益，在2005年美国总统信息技术顾问委员会给总统的报告“计算科学：确保美国竞争力”中指出，“计算科学采用先进的计算能力理解和求解复杂问题，已经成为对美国科技领导地位、经济竞争力和国家安全的关键，计算科学是21世纪最重要的技术领域之一”。就国内而言，近年来，CAE软件在土木工程中的应用越来越重要。因此，在掌握计算结构力学这门理论课程的基础上，必须要通过合适的软件实现计算的目的。这要求我们在课堂上就要给学生讲清楚其关系，这样，将来从事具体工作能够得心应手。

4 课程内容的完善与提高

首先在教学规划上，查阅国内外相关教材资料，制定课程建设计划和实施方案。根据计划制定符合专业培养目标、内容新颖、结构体系完整、重点难点突出、学时分配合理的教学大纲。其次教学内容上，如前所述我们引入一些科研内容，激励学生敢于思考，勇于实践，善于总结，培养学生创新意识。在教学观念上要树立“授之以渔”的教学观念。在教学过程中，通过一些实际工程项目的教学案例教学使得知识的学习不再是唯一的目的，而是认识科学本质、训练思维能力、掌握学习方法的手段；在教学手段上，多媒体、互联网、电子教案和计算机辅助教学等现代化教学手段和方法，同时配合传统的教学方式实现“组合式”教学。工程图片的形象和生动，能弥补课本知识的不足，提高学生的学习兴趣，能有效改善学生的学习效果，提高教学质量。

5 结语

课题组通过《计算结构力学》对这门课程进行了实践教学改革研究，在完成基本的课程内容改革基础上，着重就培养研究生的创新能力，解决实际问题的能力进行了实践研究。首先，提出在教学内容上，科研项目应于教学相结合，从而提高研究生分析问题解决问题的能力。其次，在重视理论教学的同时，也要重视实践教学环节，并要加强理论与实践教学的结合。通过组织学生参与实际工程项目来提高学生这方面的能力。希望通过这些有益的教学改革能够培养出既能动脑又能动手，既不缺乏理论知识又具备较强实践能力的复合型人才。

参考文献

动力学计算 篇4

《计算流体动力学及其应用》是能源动力方向硕士和博士研究生的一门学科基础课程。自上世纪90年代开设以来, 前后经历了三次大的变革。在1995—2000年期间, 课程名为《叶轮机械的数值模拟计算》, 只针对能源动力学院流体机械专业方向的研究生开设;自2001年起, 课程改名为《计算流体动力学及其应用》, 面向能源动力一级学科及其下属各二级学科的研究生, 并成为能源动力学科方向研究生的一门基础课程;2006年以前, 课程授课内容以计算流体动力学方面的原理和方法为主, 选课学生主要为能源动力学科方向的研究生;从2006年开始, 为适应广大研究生的选课需要, 我们对教学内容进行了适当调整, 辅以CFD商用软件的实例和应用, 以期实现理论性与实践性并重的教学理念, 并将课程面向全校研究生开放。随着计算机技术的飞速发展, 计算流体动力学的应用日益广泛。众所周知, 计算机硬件水平的提升, 将相应地促进CFD商用软件功能更加强大, 应用更加广泛, 最终使得CFD商用软件得到了前所未有的发展。同时, 随着研究生招生规模的扩大, 使得选修《计算流体动力学及其应用》课程的研究生人数大增, 从上个世纪90年代的十几个学生, 到现在的一百多个学生, 而且涉及众多学科, 比如船海、化工、建筑、电气、交通、材料、光电等。《计算流体动力学及其应用》课程的历史与现状在一定程度上给我们将要进行的教学改革提出了新的要求, 同时也为我们指明了课程建设的新方向, 值得我们深入思考, 并付之于实践。

二、课程定位

《计算流体动力学及其应用》作为一门研究生的学科基础课程, 我们在进行改革之前, 应该首先考虑它的定位。华中科技大学作为一所教育部的“985”和“211”的高校, 一直以“研究型”大学著称。学校对于研究生的培养非常重视, 导师为每一位研究生制定了详细的培养计划, 课程的选修均有所考量。基于选修《计算流体动力学及其应用》课程的研究生人数众多, 涉及的院系广泛, 经任课教师讨论, 申请学校研究生院同意, 决定将该课程定位为高水平研究生课程。所谓高水平研究生课程, 初步确立的含义为, 高水平的学者, 采用高水平的教材, 以先进灵活的形式教授课程, 旨在培养学生坚实宽广的理论基础和系统深入的专门知识。高水平课程在内容上应该具有基础性、专业性和前沿性, 前沿性可以体现在任课教师结合自己的科研实践, 在讲授中融入一些与课程相关的前沿内容。

三、教改实践

基于高水平研究生课程这样一个定位, 我们开始着手进行课程的教学改革, 具体内容包括:组建教学团队、改革教学内容、建设实践教学平台。首先, 组建一支高水平、高素质的教学团队。教学团队由三位教师组成, 他们均具有博士学位, 高级职称。其中, 课程负责人张师帅副教授, 长期从事计算流体动力学及其应用方面的教学及研究工作, 自2006年起, 一直担任该课程的主讲教师;任课教师郭照立教授, 是目前国内计算流体动力学方面的顶尖学者, 国家杰出青年基金获得者, 并具有较高的国际知名度。郭教授团队在国内外权威学术刊物和会议上发表科学论文100余篇, SCI收录90余篇, SCI引用1200余次;任课教师陈胜副教授是一位青年学者, 在格子Boltzmann算法研究方面颇有建树。将《计算流体动力学及其应用》课程建设成一门高水平研究生课程, 得到了教学团队中每一位教师的支持, 大家一致赞同经常开展教学交流, 学习先进的教学方法和教学手段, 进一步提高教学效果。其次, 我们对教学内容进行了改革。教学团队根据选课研究生人数众多, 涉及的学科方向广泛, 重新制定了《计算流体动力学及其应用》课程的教学大纲, 确保讲授内容的基础性、专业性和前沿性。课程主要内容包括:控制方程的离散化方法、流场的求解计算方法、湍流模型及其应用、网格生成与计算技术、复杂流动的介观模型和数值方法、格子Boltzmann算法及其应用、经典CFD软件的基本用法等。而对于控制方程的离散化方法, 将重点介绍有限差分法和有限体积法;对于流场的求解计算方法, 将重点介绍SIMPLE及其系列算法;对于湍流模型及其应用, 将重点介绍k-ε模型及其应用;对于网格生成与计算技术, 将重点介绍结构网格和非结构网格的生成方法以及并行计算方法。同时, 还将邀请国内外计算流体动力学方面的专家学者前来开展专题讲座。对《计算流体动力学及其应用》课程进行教学改革是全体任课教师的共同愿望, 大家积极讨论, 并提出在现有教材的基础上, 编写具有自己特色的教材等建议。在改革教学内容的同时, 教学团队还利用学院现有的条件, 建立“计算流体动力学”软件平台, 该平台拥有高性能的计算工作站, 可以开展并行计算、直接数值模拟等大型计算研究, 为研究生开展离散方法、网格生成方法、计算方法以及复杂流动模型等研究工作创造了良好条件, 同时也为对计算流体动力学方面的前沿研究课题感兴趣的大学本科生开展创新性研究工作提供了良好条件。与此同时, 该平台还拥有多种商用CFD软件, 比如FLUENT、CFX、STAR-CD、PHOENICS、Flo-EFD等, 成为广大研究生开展自主学习、自主实践、相互交流的优良场所。还可以根据研究生的需求, 安排任课教师不定期地通过软件平台为学生解惑答疑, 引导研究生探索创新, 提高学术水平。

众所周知, 研究生学术水平的高低是一所大学学术水平的反映, 更是一个国家科技创新能力的反映。研究生不仅需要扎实掌握专业基础知识, 更需要具有较强的创新意识和创新能力。目前, 高等学校在培养研究生创新能力、提高研究生学术水平方面还有待加强。为此, 本文提出了通过对《计算流体动力学及其应用》课程进行教学改革, 并将之建设成具有基础性、专业性和前沿性的高水平课程, 进一步培养研究生的创新能力, 提高研究生的学术水平。同时, 本文对实践过程中的一些具体措施和经验进行了探讨。

参考文献

[1]吴新颖.论研究生教育的课程设置[J].科技导报, 2004, (2) :40-42.

[2]李学艺, 钟佩思.机械类全日制专业学位研究生教育探索[J].中国电力教育, 2013, (1) :22-23.

[3]李静波, 柴育玲.研究生教育创新模式研究[J].教育理论与实践, 2010, (5) :10-11.

[4]别敦荣, 万卫.学术性学位与专业学位研究生培养模式的特性比较[J].研究生教育研究, 2011, (4) :77-80.

计算地球流体力学总结 篇5

对于任意函数通过傅氏分解都可以看做无穷个波的叠加，波长为2L/n（0~2L）。但是对于数值计算，区间内离散点确定，离散傅氏级数得到分波的波长为2Ndx/n（n≤N），小于该波长的波离散后必然将与可以分辨的波混淆。由此可见离散型问题不可能分辨λ<2dx的分波。

非线性耦合短波

计算的不稳定性主要来自于计算中存在相对网格过小的短波，一旦产生波长小于2dx的分量，网格系统不能正确分辨，必将产生混淆而将其折叠刀大于2dx的波上。对于非线性发展方程，这样的耦合产生短波又折叠混淆的过程不断重复，即构成短波能量的虚假增长而导致计算的不稳定。克服了初始误差短波的增长，也就获得了计算的稳定。

uuu（00x1）当一原始usin(kx)平流耦合，txuu1uksin(kx)cos(kx)ksin(2kx)只要有k波，经都会出现新的小x22xx于格距的短波。

计算的不稳定主要来源

相对网格过小的波，1、对于线性问题，克服了初始误差短波的增长就获得了稳定性

2、对于非线性问题，即使初始条件中不含短波，由于非线性耦合作用也会不断产生短波，由此产生不稳定。

3、初值的选择同样会引起不稳定。这些都是由空间离散化造成的，即使步长减小也不能克服。差分格式抑制短波

差分离散后，原波长为n2Nx/(2mNl)2x的分波，就表现为n2Nx/l2x故离散型问题不可能分辨出波长2x的分波，（与差分分辨率一致），相对网格过小的短波是计算不稳定的主要来源，抑制短波的发展也就获得了计算的稳定。

有限元数值模型、不规则网格有限差分数值模型的原理、步骤、差异是什么？

有限元技术的思想：既然整体试函数难于选择和确定，就不如将求解域分成若干小区，在每个元内都用一个简单的（如空间线性的、二次的等）函数作为元内这一小局部试函数，并以某种光滑性进行联结，以构成最终的全域试函数，再依变分法或权余法求得所有小元内的试函数的待定系数，整体函数也就自然确定了。虽然一般说来这样确定的整体函数可能不如经典形式的函数具有整体的各阶光滑性，但可具有整体的连续性和低阶光滑性。（如果采取计算量较大的“谱元法”，则高阶光滑性也是可以达到的）。

不规则网格差分法，其网格分布可沿用有限元的单元分布，原来的“元”，即是差分法的网格或剖分网格。所不同的是差分网格的边长比求解域的尺度必须是个小量，而有限元并非一定要如此。

不规则网格差分法比有限元法的计算量可以节省很多，数值格式比较方便，当然要满足相应的稳定性条件。在有限元法中，水位在三角形内的差商表达式是由线 性插值函数得到的，且元内角顶三点以逆时针方向顺序计数。而在不规则网格差分法中，是利用边长比计算域水平尺度为小量时，三顶点上的函数值，在其内任一点P（x，y）展开Taylor展式，略去二阶小量项之后联立得到的。

有限元计算方法的三个优点：①很好地弥和岸界②网格疏密可依所计算的流场的几何与物

理特性而灵活确定③可以将自然边界条件融入内点的计算方程之中，不必单独列出。

隐式格式

显式：对两时间层格式，知道n时刻各空间层上函数值而推n+1时刻的值。隐式：包括n+1时间层上二个或多个节点处的未知值，用n时刻各空间点值不能直接得出n+1时刻各点值，必须联立求解一个与网格点数目相同的方程组。

与显示格式不同，隐式差分格式，使用隐格式求解，每个时刻包含较多的计算量。但是在稳定性上往往优于显式，因而对时空步长放宽要求，可以减少计算量。

有限体积法（FVM）又称为控制体积法。将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解的微分方程对每一个控制体积积分，便得出一组离散方程。其中的未知数是网格点上的因变量的数值。为了求出控制体积的积分，必须假定值在网格点之间的变化规律，即假设值的分段的分布的分布剖面。

有限体积法可视作有限单元法和有限差分法的中间物。有限单元法必须假定值在网格点之间的变化规律（既插值函数），并将其作为近似解。有限差分法只考虑网格点上的数值而不考虑值在网格点之间如何变化。有限体积法只寻求的结点值，这与有限差分法相类似；但有限体积法在寻求控制体积的积分时，必须假定值在网格点之间的分布，这又与有限单元法相类似。在有限体积法中，插值函数只用于计算控制体积的积分，得出离散方程之后，便可忘掉插值函数；如果需要的话，可以对微分方程中不同的项采取不同的插值函数。

开边界条件

再用数值方法对某一海域的某种海洋环境问题做计算模拟或预测时，只要不是对全封闭海计算，必然会有一人为划定的水中边界，该边界处各种状态与相邻的内域并无本质差异，成为开边界。开边界给出的优劣，往往对域内的计算响应很大。开边界条件类别：①确定（第一类）边界条件；②辐射条件；③强迫波开边界条件；a.辐射强迫条件；b.强迫波与自由波分别计算；c.时间分段条件；④海绵（Sponge）条件和无穷远边界条件；⑤按水深分段条件；⑥统一格式

代数坐标变换与微分坐标变换的差异是什么？

答：这两类坐标变换都将实际物理空间中的问题变化到变换空间去进行计算。代数坐标变换是两空间中的坐标关系以代数式相关，微分坐标变换是两空间的坐标关系以一种微分方程相联系。代数坐标变换对弥和一段弯曲度不甚严重的陆架岸界是简便易行的，但对弯曲度较大的海岸，特别是对封闭或半封闭海湾则不易实现，而微分坐标变换可以弥补这个缺点，能够很好的弥和弯曲度较大的岸界。它利用的是调和函数的性质：①解是唯一的或解只差一常数；②极值不可能出现在域内，即最大值和最小值均在边界上，选用调和函数为坐标函数。

干湿网格方法、坐标变换模型怎么使用？使用步骤是怎么样的？每步怎么做？做些什么内容？它们的优缺点是什么？

对于岸边地形较为平缓的潮滩低地，对岸边流场特别是漫水位置干湿网格方法特别关注，固定岸边界模式已不适应。干湿网格方法核心是建立一套判别准则，再任一瞬时对方程组进行数值计算前，判断哪网格干湿，从而有了瞬时交界线作为岸界，并在此给定边界条件。由于每一时刻判断出的岸界随时间变化，因此可以模拟出任一瞬时的淹水位置和整个过程的最大淹水范围： 干湿网格优点：如果已经具有了固定岸界的某个数值程序及其计算程序，那么只需在程序中加入一段干－湿判别（从而确定瞬时边界位置）的子程序，下一步的流场函数的计算则可以完全采用原有的计算程序，而岸边界的边界条件仍是法向流速为零。缺点：岸边界条件不符合流体运动学原理。它采用的是与固定岸边界相同的边界条件－法向流速为零。虽然它的边界依判别准则是在移动着，但那都是出于某种技术上的考虑而并非是从根本原理上得到的。

Johns 代数坐标变换模型：优点：岸边界条件具有很好的物理意义，符合流体运动学原理。缺点：对岸界的弯曲程度有一定的限制。

使用干湿网格方法时（Flather-Heaps模型），用连续方程计算出每个水位计算点的水位，从而组成瞬时水深H，在计算流速分量之前，要判断流速计算点的干湿情况。使用Johns 代数坐标变换模型时，所需增加的内容就是在每一时间层利用同样的边界条件计算出方程中的变量后，还需找到新的岸界和网格对应位置。

Heaps 谱方法的思想是怎样的？如果本征序列不好，计算量怎样？

答：思路：以谱方法将一个线性正压海洋的三维模式化为一系列二维方程进行数值求解。首先将两个水平流速分量在垂直方向上做谱展开，即令：

uArFr(z)vBrFr(z)，其中 ArAr(x,y,t)，BrBr(x,y,t)

r1r1然后解一个仅考虑潮的三维近海动力学的垂向本征问题，得到一系列本征值和对应的本征函数序列，且二者是正交的，可将其他量同u、v展开。接着以每个本征函数为权函数，对两个水平运动方程导出垂向 Galerkin 权余式，引用本征问题和流速的海面与海底边界条件，可将权余式化为二维方程组，最后利用本征函数归一化和本征函数正交性，得到流速u、v。

流速分解方法的提出是怎么提出的？它的物理意义是什么？

答：准平衡模式的流速分解：产生水平流速的两个作用力，一个是由水位梯度所表达的压强梯度力，一个是海面风应力，即梯度流和风海流。含时变项的流速分解：只取决于海面梯度而与深度无关的正压梯度流和一个随深度而变化的Ekman流，但两者皆考虑其时变效应。

σ 坐标

答：σ 坐标是一种适应地形的坐标变换，它既可以保证垂直分辨率一致，也给垂向分层计算带来很大方便。(z-ζ)/H 将海域在垂直方向变为无因次单位厚度，在这一变换海域中，垂向按该坐标分层可以方便在全海域得到同样的层数。

 坐标的一般形式：（－1，0）区间：

(ba)zbH，其中a，b为两个常数。

zhzH（0，1）区间：H

在陆地附近计算斜压部分的水平压强梯度力会出现问题，在线性模型中，Z坐标

Pzdz,取坐标，即,Zh在稳定的层花海洋下的x方向斜压梯度力为：0Xg中，这两项的永远是符号相反的。因此，在陆地处或海山附近，hp很大，此时（）xx也必然很大，从而一个较小的量由两个大量之间的差来决定。由此可知，在这一区域进行数值离散计算时，精度非常重要。

有限解析法 原理：将计算域划分为若干子区域，在每个子区域上将在整个求解域上不为常数的参数近似地认为是不变的，即为局地常数，从而可求出在各个子区间的解析解。（（比有限差分更具物理意义。如在对流项应用迎风格式，在有限解析法中自动满足。））

优越性主要表现在2各方面：

一、这样形成的方程，在每个子区间都保持了解析解的特征，而解析解能自动形成迎风性质，从而保证了物理上的输运性。

二、有限解析法求解定常对流扩散方程，该格式可将方程离散为对角占有型的代数方程组，从而很好的保证代数方程组的收敛性。ADI模型比全隐格式节省计算量，又比半隐半显格式提高了稳定性，从而加大了时间步长。

Taylor展式构造差分方程：将函数在一时空点上按时间做Taylor展开，按需要的计算精度截断，用微分方程将时间导数化为空间导数，再将空间导数用差分格式取代，即构造差分方程。

mu分裂算子原理：若微分发展方程有若干空间算子，即形如Ll(u)，则可以

t1在微小时间段内，分裂成若干发展方程求解，如

1uLl(u)(l1,2,3....m)mtxL分辨率：在微商化为差商的离散过程中，当x对于不同波长的分辨率能力不同，xL/2。各种差商都失去了分辨L波的能力，成为差分分辨率。

物理解：微分方程的差分方程解中与微分方程的解析解相同的解，符合物理规律的解。计算解：反之，不符合的解 待定系数法：线性齐次微分方程进行差分离散后，化为某时空离散点的函数值与周围若干点的函数值的代数方程，依次思路，可事先确定一含待定系数的代数方程，再依计算方程和Taylor展开需要的精度，最终确定这些系数。守恒性

总体海水质量守恒，二维：总机械能守恒，总体物质守恒，总体涡度守恒。二维运动边界：总体水量守恒，总机械能守恒，总体涡度守恒和总物质量守恒。三维：总机械能守恒，总物质量守恒 非线性计算不稳定的特点：

1计算不稳定是由于方程的非线性性质和非线性作用产生，一半不能用缩短时间步长的方法克服。

2、计算不稳定具有突变型。

3、不仅与差分格式有关还与初始条件有关

4、与方程类型及解得性质有关。

非线性计算不稳定的产生原因：

1、由于差分格式的混淆误差引起

2、由于能量守恒关系被破坏引起

3、由于函数网格变号引起

4、与问题的初始状态有关

5、由非线性解的非光滑性或激波解得性质引起

6、与非线性方程解的分岔、混淆性质有关。

2解决办法：

1、进行空间或时间的平滑

2、在物理变量F的方程中加入扩散项kF,即加入人工粘性。

3、构造具有隐式平滑或具有选择性衰减作用的差分格式

4、构造总能量守恒的差分格式。基础模型——地转平衡：水平压强梯度力与科氏力平衡。半地转平衡下的Rossby波模型：平面上，东西向地转平衡无运动，南北向惯性力与压强梯度力有运动。惯性平衡：水平惯性力与科氏力平衡Ekman平衡：湍粘性力和科氏力平衡。

非旋转模型——Laplace平衡：局部惯性力与压强梯度力（海面坡度引起）平衡。准平衡风增水模型（风暴潮）潮位压强梯度力与海面风应力平衡。

1，什么是斜压海洋？什么是正压海洋？一般情况下，模拟潮波不需要考虑斜压效应，但在什么情况下需要加入斜压效应？

答：正压海洋应理解为均匀海洋，即流场密度为常数；斜压海洋为非均匀海洋，即是层化海洋，密度不再是常数。在模拟内潮时需要考虑斜压项。

2，近海风海流模式应如何改进（讨论底边界条件）？若仍要在近海使用风海流模式，那么有什么缺点？

答：当在近海时，水深非常浅，底应力会非常强，必须考虑海底摩擦作用，这时方程为 fVVV()0，边界条件：z＝0（海面）a，z＝－h，： zzzVkV（底应力正比于底流速）。z5、对流项：依据海水微团中污染物的增加量取决于由海水带入微团的物质量而导出。扩散项：根据Fick定律导出。

6、lagrange余流： 对特定海水微团的流速进行跟踪平均得出的流速分量。Euler余流：是指在固定空间点上分解或潮平均出的流速分量。

7、近海二维动力学方程组中各项的物理意义是什么？

uuu1paaxbxuvfvgL2u txyxxHvvv1paaybyuvfugL2v txyyyHHuHv0 txy局部惯性力、平流项、科氏力、海平面倾斜造成的压强梯度力、海面大气压不均运造成的压强梯度力、风应力、水平湍黏性力

8、正确的运动岸边界条件是怎么提出的？

答：物质面上的流体质点总是沿着该面运动，它们的流速在界面法向的分量与界面几何点的法向移动速度相等，即 vnvsn（A式），v 为界面上流体质点的流速，vs 为界面几何点移动速度，n 为界面单位法线向量。如果运动界面可以用一个以时间变量t 为参数的曲面方程表示，即如 F(x,y,z,t)=C , C 为常数，则F的全导数必为零 dFFFdxFdyFdz0（B式）

dTtxdtydtzdt而其中的dydxdz=u , =v, =w, 即是该曲面上之几何点随时间的移动速度在dtdtdt（x,y,z）三个方向的分量。且由曲面方程还可知，该曲面之法线与F的梯度同向，即 nF。|F|FvsF0 t此时 B 式又可写成

依流体运动连续原理A式，则最终获得运动界面流体质点的边界条件 FvF0 t在最一般的情况下，运动岸边界的方程可取H＝0，则普遍形式的运动学边界条件为

HHHHHvn0 uv0 或者 tntxy9、在推导层内平均运动方程过程中，使用了哪些近似？ 答：1），界面处流速的导数以相邻两层间的差商取代，2），umzm

zm1udz 代替了zmzm1udz vm2zm1zmudz 代替了zm1zmuvdz

10、近海物质的输运机制及输运方程是什么？

答：如果某种物质在海水中没有源和汇存在，也不会因某种化学物质而产生或消亡，这种物质虽在海域内输运，但能量保持守恒，这种物质即称保守物质，它们的输运只通过对流、扩散完成。

海洋中的运移与变化可以归纳为三种机制，并可以统一在一个方程中表达出来，这就是对流（平流和垂直对流）、扩散和转化。在考虑垂直输运的问题中，将悬浮质的沉降也归于对流之中，转化则看作物质的源汇。

CiC(VHVM)HCi(wws)i tz＝C(Ki)KL2HCiQi(Cj)zzCi为某一状态变量的浓度，i,j1,2,......,N，即由N个状态变量组成一生态系统或水质系统；VH为水平流速；VM为生物群体水平回游（Migration）速度，ws为其沉降速度，K为垂直扩散系数；KL为水平扩散系数；Qi为源汇函数。

11、正压分层平均模式

自海底到海面分为若干层，认为确定层数和每层厚度，然后在每层内类似二维的做法，将三维方程组在层内平均为二维。由于水深不同，分层后海底近似阶梯状。若自海底到海面分为m层,第m层厚度为hm,下层坐标为zm上层为zm1,从而垂直流速不平均，只在边界处求值，hm=zm-zm1,设该层之层平均水平流速为Vm，Vm取第m层的中点。

12、垂直流速

动力学计算 篇6

关键词：回转扭矩；转速计算

1.概述

面对紧迫的形势，全液压钻机以其优良的性能，逐渐取代传统的立轴式岩心钻机。现已HYDX-5A型钻机为例，阐述动力头扭矩、转速的计算方法。

2.动力头的设计计算方法

2.1液压系统的计算

2.1.1 柴油机的型号及参数

HYDX-5A型钻机采用康明斯6CTA8，.3-C195型柴油机，额定功率/转速： 145KW/1900rpm，最大扭矩/转速：820Nm/1500rpm。

2.1.2 液压泵型号及参数

钻机采用萨澳90L75型闭式泵，最大排量：Vg泵=75ml/r，系统最高压力ΔP=320bar，泵输入转速 n泵=1500rpm。

2.1.3 液压马达型号及参数

钻机采用萨澳51L80型马达，最大排量：Vg马达max=80.7ml/r，最大排量：Vg马达max=，16.1ml/r。

2.1.4液压系统的计算

假设液压系统无泄漏的情况下，液压泵的输出流量等于液压马达的输入流量。

Q泵=Qe马达=Vg泵×n泵×1/1000=75ml/r×1500rpm×1/1000=112.5L/min

当Vg马达max=80.7ml/r时，获得最低输出转速nmin和最高输出扭矩Memax

nmin=Qe×1000×?v/Vg马达max=112.5L/min×1000×1/80.7ml/r=1394rpm

Memax=Vg馬达max×ΔP×?mh/20×3.14=80.7ml/r×320bar×1/62.8=411Nm

当Vg马达min=16.1ml/r时，获得最高输出转速nmax和最低输出扭矩Memin

nmax=Qe×1000×?v/Vg马达min=112.5L/min×1000×1/16.1ml/r=6988rpm

Memin=Vg马达min×ΔP×?mh/20×3.14=16.1ml/r×320bar×1/62.8=82Nm

2.2 机械系统的计算

2.2.1 变速箱高速，主轴箱低速

N1主轴=n马达输出转速× 齿轮1（?92mm）/齿轮2（?180mm）×齿轮3（?156mm）/齿轮4（?116mm）× 齿轮5（?80mm）/齿轮6（?235mm）×齿轮7（?174mm）/齿 轮8（?342mm）=0.119×n马达输出转速

M1主轴=M马达输出扭矩× 齿轮2/齿轮1×齿轮4/齿轮3× 齿轮6/齿轮5×齿轮8/齿 轮7=8.4Me马达输出转速

2.2.2 变速箱高速，主轴箱高速

N2主轴=n马达输出转速× 齿轮1/齿轮2×齿轮3/齿轮4× 齿轮9（?165mm）/齿轮10（?150mm）×齿轮7/齿 轮8=0.385×n马达输出转速

M2主轴=M马达输出扭矩× 齿轮2/齿轮1×齿轮4/齿轮3× 齿轮10齿轮9×齿轮8/齿 轮7=2.6Me马达输出转速

2.2.3 变速箱低速，主轴箱低速

N3主轴=n马达输出转速× 齿轮1/齿轮2×齿轮11/齿轮12× 齿轮5/齿轮6×齿轮7/齿 轮8=0.045×n马达输出转速

M3主轴=M马达输出扭矩× 齿轮2/齿轮1×齿轮12/齿轮11× 齿轮6/齿轮5×齿轮8/齿 轮7=22Me马达输出转速

2.2.4 变速箱低速，主轴箱高速

N4主轴=n马达输出转速× 齿轮1/齿轮2×齿轮11（?92mm）/齿轮12（?180mm）× 齿轮9/齿轮10×齿轮7/齿 轮8=0.146×n马达输出转速

M4主轴=M马达输出扭矩× 齿轮2/齿轮1×齿轮12/齿轮11× 齿轮10/齿轮9×齿轮8/齿 轮7=6.85Me马达输出转速

2.3主轴理论输出转速、扭矩

2.3.1 当nmin=1394rpm Memax=411Nm 时

N1主轴=166rpm M1主轴=3452Nm N2主轴=537rpm M2主轴=1069Nm

N3主轴=62rpm M1主轴=9042Nm N2主轴=203pm M2主轴=2815Nm

2.3.2 当nmin=6988rpm Memax=82Nm 时

N1主轴=832pm M1主轴=688Nm N2主轴=2690rpm M2主轴=213Nm

N3主轴=315rpm M1主轴=1804Nm N2主轴=1020pm M2主轴=561Nm

3.结论及建议

通过对动力头结构及计算方法的深入剖析和计算，加深了钻机设计人员对钻机设计理念的理解，同时也便于钻机用户在实际使用时，更好地掌握钻机的操纵规律，创造出更好的钻探效果。

参考文献：

[1]成大先.机械设计手册.北京.化学工业出版社 2002.

动力学计算 篇7

人口是反映国情、国力基本情况的重要指标,是研究区域经济必须考虑的重要因素之一,是分析人口现状、制定发展规划首先要考虑的基本问题。科学的人口预测是制定和实现社会经济各项战略设想的出发点,是制定正确人口政策的科学依据。

长江三峡库区是中国西部大开发的组成部分。三峡工程自1993年开始兴建,到2009年建成,最终正常蓄水水位高达175米,淹没涉及湖北的的夷陵、兴山、秭归、巴东和重庆的万州、涪陵、渝北、巴南、巫溪、巫山、云阳、奉节、开县、忠县、石柱、丰都、武隆、长寿、江津以及重庆市区20个区县(自治县、市),总面积约为5.83万平方公里。三峡库区是全国18个集中连片贫困地区之一,重庆库区15个区县中,有国家级扶贫开发重点工作县10个,市扶贫开发重点县两个,绝对贫困人口64.4万,占全市贫困人口总量的70%,相对贫困人口193.9万,占全市总量的70%左右。

重庆市直辖以来,库区人口控制问题取得了很大成绩。据统计,2004年与1997年比较,人口出生率由12.2%下降为9.61%;人口自然增长率从5.29‰降到3.83‰;合法生育率由88.73%提高到92.39%;为全市经济社会全面协调和可持续发展创造了良好条件。然而,按照全面建设小康社会、建设长江上游经济中心、构建和谐社会的要求,库区人口和计划生育工作仍然存在一些突出问题,例如三峡库区人口稠密、经济不发达、城镇化程度低、生态环境脆弱,人口承载能力极差等问题尤为突出。

为了改善人口结构、提高人口素质、稳定低生育水平、引导人口合理分布,加强库区人口预测和控制,因此建立一个人口增长预测数学模型对三峡库区人口增长的短、中、长期发展趋势做出预测就显得十分重要了。

二、三峡库区人口预测动力学模型

英国统计学家Malthus曾得到经典的人口发展模型[1]:

设时刻t=t0时人口数为x0,时刻t的人口为x(t),x(t)可视为连续、可微函数,r是人口增长率(常数),于是x(t)满足微分方程的初值问题:

解初值问题(1)得:

由此模型可以看出,人口总数是按指数无限增长的。当t→∞时,x(t)→∞(r>0)。

显然,此人口模型只能用于短期预测,因为生存空间的环境所提供的条件只能供养一定数量人口,当时间趋近于无穷大的时候,人口总数不可能达到无穷大,所以此模型不能用于长期预测,必须进行修改。

由于中国实行计划生育政策,人口增长率r是时变的。为此,我们修正Malthus人口模型,即得到三峡库区人口发展预测动力学模型:

式中假设r=r(t)=a-b(t-t0)为时变的人口增长率,a和b为参数。求解得:

注意到(3)中的初始条件,则初值问题(3)的解为:

另外,如果已知t-t0时,r=r0及t=t时,r=rt,由假设可得如下方程组:

解方程组(5),可求得a和b。于是将a和b以及初始值P0代入(4)式中,即可求得三峡库区在t年的人口总数P(t)。

根据长江三峡库区生态、环境、土地、矿产等资源在一定时间内允许承载人口数量的极限,将三峡库区人口最大限度Pn=2 200万人[2],依据重庆市政府确定的人口自然增长率r0=3.83‰[3],以2010年为初始基年,t0=2010,t=2060,根据全国第六次人口普查公布的人口数据,P0=1 900万人(见表1)。则按(4)式计算,2060年人口增长率按线性规律变化时,应该下降到rt。于是,rt可如下求出:

解方程(6)得:

三、模型的Mathematica软件求解

将a,b,P0,t0和t代入(4),用Mathematica计算代码如下:

由上可知,从2010—2060年,全国人口增长率按线性规律减小时,则应该下降到2.034‰。据此可推算出各年的人口增长率和人口总数应该(如表2所示)。

用MATLAB做出三峡库区2010—2060年人口趋势图,程序如下:

图像(如下图所示):

从以上三峡库区人口预测表和图形可以看出,到2060年把三峡库区人口控制在最大环境人口容量2 200万以内,必须下大力气实行计划生育,采取有效措施,确保人口增长率逐步下降到2.03‰。为此,提出几点对库区人口发展工作的建议。

四、库区人口发展工作建议

人口多、底子薄、人均资源相对不足是直接影响库区经济发展和计划生育工作的瓶颈。采取切实有效措施解决人口问题将有利于建设长江上游经济中心和实现全面小康社会。为此笔者建议:

1. 坚持现行“城区稳、库区降”的生育政策。库区生态、环境、粮食、不可再生资源承载力有限,要实现库区人口、资源、环境和经济、社会协调发展和可持续发展,必须稳定现行生育政策不动摇。

2. 开展“关爱女孩行动”,采取综合措施降低出生人口性别比。加强对广大农民进行“生男生女都一样”的宣传教育,对其独生子女按照城市独生子女待遇同等对待;采取有效措施,提高女孩及其家庭的社会政治地位,提高女性受教育水平、科学文化素质,拓宽妇女就业门路,增加就业机会。多措并举,综合治理出生人口性别比失衡问题。

3. 大力发展库区生产,繁荣库区经济。推动城镇化进程,扩大就业,增加老百姓收入,繁荣库区经济,并增加对人口和计划生育事业经费投入力度。鼓励和提倡人口合理流动,繁荣库区经济,减小库区居民生存压力,推进城市化进程;加强生态移民,防止库区生态环境遭到破坏,增强库区承载和防灾能力。

4. 积极自主研制和吸取国外先进技术,引进长效避孕药物,对农村尽可能做到免费或成本价供应;大力推广农村合作医疗保险,并给予一定比例的津贴,让农民也享受与城市居民一样的国家公共福利,这有利于进一步消除一对农民夫妇只生一个独生子女的后顾之忧。

5. 加大教育和人才投入的力度,大力开发人才资源,充分发挥科教和人才对库区建设的主导作用,为库区经济社会发展营造良好的人才环境,把贫困地区的发展方向调整到以人力资本投资为主的发展道路上来,把扶贫开发转移到依靠科技进步提高人口素质的轨道上来。

摘要：以中国长江三峡库区人口作为研究对象,建立人口发展动力学模型,以2010年三峡库区人口数为初始值求解了所建立人口动力学模型。运用软件对数据进行计算处理,并利用软件做出了三峡库区2010—2060年人口发展趋势图,由此可预测出长江三峡未来人口发展状况。针对该地区未来人口预测,提出了人口发展控制的若干建议。

关键词：三峡库区,人口预测,动力学模型,人口控制

参考文献

[1]姜启源,谢金星,叶俊.数学模型[M].北京:中国统计出版社,2003.

[2]唐晓平,舒克盛.三峡库区人口承载力的动力学模型研究[J].南方人口,2007,(4).

[3]温永高,杨应敏,欧阳文.三峡工程重庆库区人口计生工作现状、问题及对策研究[J].人口与计划生育,2005,(5).

[4]重庆市统计局.重庆市2010年第六次全国人口普查主要数据公报[N].重庆日报,2011-05-03.

[5]宜昌市统计局.宜昌市2010年第六次全国人口普查主要数据公报[N].三峡日报,2011-05-30.

动力学计算 篇8

凸轮机构在机械工程领域被广泛应用。目前, 用于研究凸轮机构动力学的模型可以归纳为四种形式。第一种模型将凸轮机构看作单自由度的多刚体系统, 这种模型最为简单, 并被现有的一些机械设计教科书[1,2,3,4,5,6]所采用。第二种模型[7,8,9,10,11,12]将推杆看作可发生拉压变形的柔性体而将其他构件看作刚体的柔性多体系统。文献[7,8,9,10,11,12]采用第二种模型研究了凸轮机构的有关动力学特性。第三种模型[13,14,15,16]将凸轮轴看作可发生扭转变形的柔性体而将其他构件看作刚体的柔性多体系统。文献[13,14,15,16]应用第三种模型研究了计入凸轮轴柔性的凸轮机构的动力学行为。第四种模型[17]将推杆和凸轮轴都看作柔性体而将其他构件看作刚体的柔性多体系统模型。相对前三种模型而言, 第四种模型更符合实际, 但同时又是更为复杂的一种模型。

以上四种模型的共同缺陷是均没有考虑凸轮与推杆之间的摩擦。鉴于凸轮机构中运动副的摩擦主要表现为凸轮与推杆之间的摩擦, 而这种摩擦对机构系统的运动又具有一定的影响。因此, 在凸轮机构的动力学研究中, 计入凸轮与推杆之间的摩擦可有效提高其动力学研究的精确性, 使研究结果更加符合实际。基于此, 本文在计入凸轮轴柔性的基础上, 考虑凸轮与从推杆之间的摩擦, 建立了一种新的偏心凸轮机构的动力学模型。以该动力学模型为基础, 通过数值方法研究了偏心凸轮机构的运动响应, 分析了凸轮与推杆之间的摩擦以及凸轮轴的柔性对该机构运动的影响, 得到了相应的结论。

1 动力学建模

1.1 凸轮轴的扭转变形

图1所示为一偏心凸轮机构, 该机构由驱动盘、凸轮轴、偏心凸轮、推杆和机架组成。驱动盘、凸轮轴和偏心凸轮固连于一体, 作为主动件。机构运行特别是高速运行时, 凸轮轴的柔性将使凸轮轴不可避免地发生扭转变形或扭转振动, 这将导致推杆出现运动误差, 因此, 在凸轮机构的动力学研究中, 计入凸轮轴的扭转变形是完全必要的。要计入凸轮轴的扭转变形就需要首先确定用何种方法来描述凸轮轴的扭转变形。Cveticanin[17]在研究凸轮机构的运动稳定性时, 采用无质量的扭簧描述凸轮轴的扭转变形。这种方法虽然简单, 但不能有效描述凸轮轴各截面处的扭转变形, 也无法将凸轮轴的转动惯性计入动力学模型。为了克服这些缺陷, 本文应用有限元方法描述凸轮轴的扭转变形。为此, 将凸轮轴等分为n个轴单元, 并从毗邻驱动盘的轴单元开始依次编号 (编号为1, 2, …, n) 。图2所示的是其中的第i (i=1, 2, …, n) 个轴单元, 图2引入三套坐标系, 分别是固连于机架的惯性参考系O0X0Y0Z0、固连于驱动盘的动参考系O0XYZ和单元坐标系OiXiYiZi。上述三套坐标系的坐标轴O0Y0、O0Y和OiYi均重合于凸轮轴的中轴线。

采用线性插值函数模拟轴单元的扭转变形, 则第i个轴单元上任意横截面处的扭转角Θi (yi, t) 可表达为

Θi (yi, t) = (1-yi/li) θi-1 (t) +yiθi (t) /li (1)

式中, θi-1 (t) 、θi (t) 分别为第i个轴单元两端面的扭转角, 后文简写为θi-1、θi;yi为点Oi至第i个轴单元上任意横截面的距离;li为第i个轴单元的长度。

考虑到凸轮轴被等分为n个轴单元, 因此, 有

li=l/n (2)

式中, l为凸轮轴的长度。

1.2 偏心凸轮机构的动能

为了应用拉格朗日方法导出偏心凸轮机构的动力学方程, 需要写出该机构的动能关于广义坐标和广义速度的函数表达式。为此首先选取机构的广义坐标, 本文选取驱动盘的转角φ (t) (后文简写为φ) 和凸轮轴的各轴单元端面的扭转角θi为机构的广义坐标。

驱动盘的动能为

$E{}_{\text{k}\text{d}}=J{}_{\text{d}}\dot{\phi }{}^2/2$ (3)

式中, Jd为驱动盘对其转轴的转动惯量。

凸轮轴上的第i个轴单元的动能为

$E{}_{\text{k}i}=\frac{1}{2}{\displaystyle \int }{}_0^{l{}_i}\rho {\rm I}{}_{\text{p}}[\dot{\phi }+\dot{\Theta }{}_i(y{}_i,t)]{}^2\text{d}y{}_i$ (4)

式中, ρ为凸轮轴的密度;Ip为凸轮轴的横截面对其中心的极惯性矩。

将式 (1) 和 (2) 代入式 (4) , 运算后得到Eki关于广义速度的表达式:

$E{}_{\text{k}i}=\frac{\rho l{\rm I}{}_{\text{p}}}{6n}[\theta {}_{i-1}^2+\dot{\theta }{}_{i-1}\dot{\theta }{}_i+\theta {}_i^2+3\dot{\phi }(\dot{\phi }+\dot{\theta }{}_{i-1}+\dot{\theta }{}_i)](5)$

将各轴单元的动能相加, 便可得到整个凸轮轴的动能, 即

$\begin{array}{l}E{}_{\text{k}\text{c}\text{s}}={\displaystyle \sum _{i=1}^{n}E}{}_{\text{k}i}=\frac{\rho l{\rm I}{}_{\text{p}}}{6n}{\displaystyle \sum _{i=1}^n[}\dot{\theta }{}_{i-1}^2+\\ \dot{\theta }{}_{i-1}\dot{\theta }{}_i+\dot{\theta }{}_i^2+3\dot{\phi }(\dot{\phi }+\dot{\theta }{}_{i-1}+\dot{\theta }{}_i)](6)\end{array}$

凸轮的动能可表达为

$E{}_{\text{k}\text{c}}=J{}_{\text{c}}\dot{\text{ϕ}}{}^2/2=J{}_{\text{c}}(\dot{\phi }+\dot{\theta }{}_n){}^2/2$ (7)

式中, Jc为凸轮对其转轴的转动惯量;ϕ为凸轮的转角。

推杆的动能可以表达为

$E{}_{\text{k}\text{f}}=m{}_{\text{f}}\dot{h}{}^2/2$ (8)

其中, mf为推杆的质量;h为凸轮的轴心O到推杆底面的距离 (图1) , 显然h可以被表达为

h=R+esinϕ=R+esin (φ+θn) (9)

式中, R、e分别为凸轮的半径和偏心距。

将式 (9) 代入式 (8) 后得到推杆的动能关于广义坐标和广义速度的表达式:

$E{}_{\text{k}\text{f}}=m{}_{\text{f}}e{}^2(\dot{\phi }+\dot{\theta }{}_n){}^2\cos {}^2(\phi +\theta {}_n)/2$ (10)

将式 (3) 、式 (6) 、式 (7) 和式 (10) 相加, 得到整个机构的动能关于广义坐标和广义速度的表达式, 即

$\begin{array}{l}E{}_{\text{k}}=\frac{1}{2}J{}_{\text{d}}\dot{\phi }{}^2+\frac{1}{2}(\dot{\phi }+\dot{\theta }{}_n){}^2[J{}_{\text{c}}+m{}_{\text{f}}e{}^2\cos {}^2(\phi +\theta {}_n)]+\\ \frac{\rho l{\rm I}{}_{\text{p}}}{6n}{\displaystyle \sum _{i=1}^n[}\dot{\theta }{}_{i-1}^2+\dot{\theta }{}_{i-1}\dot{\theta }{}_i+\dot{\theta }{}_i^2+3\dot{\phi }(\dot{\phi }+\dot{\theta }{}_{i-1}+\dot{\theta }{}_i)](11)\end{array}$

1.3 偏心凸轮机构的势能

偏心凸轮机构的势能包括弹性势能和重力势能两部分。弹性势能表现为凸轮轴的弹性势能和回复弹簧的弹性势能。凸轮轴上第i个轴单元的弹性势能为

$E{}_{\text{p}i}=\frac{1}{2}{\displaystyle \int }{}_0^{l{}_i}G{\rm I}{}_{\text{p}}(\frac{\partial \Theta {}_i(y{}_i,t)}{\partial y{}_i}){}^2\text{d}y{}_i$ (12)

式中, G为凸轮轴的剪切弹性模量。

将式 (1) 、式 (2) 代入式 (12) , 运算后得到

$E{}_{\text{p}i}=\frac{nG{\rm I}{}_{\text{p}}}{2l}(\theta {}_i-\theta {}_{i-1}){}^2$ (13)

将各轴单元的弹性势能相加, 便可得到整个凸轮轴的弹性势能, 即

$E{}_{\text{p}\text{c}\text{s}}={\displaystyle \sum _{i=1}^{n}E}{}_{\text{p}i}=\frac{nG{\rm I}{}_{\text{p}}}{2l}{\displaystyle \sum _{i=1}^n(\theta {}_i-\theta {}_{i-1})}{}^2$ (14)

回复弹簧的的弹性势能可写为

Eps=kλ2/2 (15)

式中, k、λ分别为回复弹簧的刚度和变形量。

设推杆运动至最低位置时, 回复弹簧的变形量为零, 这样当推杆运动至任意位置时, 回复弹簧的变形量为

λ=e (1+sinϕ) =e[1+sin (φ+θn) ] (16)

将式 (16) 代入式 (15) 后, 得到回复弹簧弹性势能的表达式:

Eps=ke2[1+sin (φ+θn) ]2/2 (17)

如果规定凸轮轴的中轴线所在的水平面为零重力势能面, 则整个机构的重力势能为

Epg=g[ (mc+mf) esin (φ+θn) +mf (R+a) ] (18)

式中, mc为凸轮的质量;g为重力加速度;a为推杆的质心到推杆底面的距离。

将式 (14) 、式 (17) 和式 (18) 相加后, 得到整个机构的势能关于广义坐标的表达式:

$\begin{array}{l}E{}_{\text{p}}=\frac{nG{\rm I}{}_{\text{p}}}{2l}{\displaystyle \sum _{i=1}^n(\theta {}_i-\theta {}_{i-1})}{}^2+\frac{1}{2}ke{}^2[1+\sin (\phi +\\ \theta {}_n)]{}^2+g[(m{}_{\text{c}}+m{}_{\text{f}})e\sin (\phi +\theta {}_n)+m{}_{\text{f}}(R+a)](19)\end{array}$

1.4 广义非有势力

作用在机构上的非有势力包括作用在驱动盘上的驱动力矩M、作用在推杆上的工作载荷P, 以及凸轮与推杆之间的摩擦力 (在凸轮机构运行中, 运动副的摩擦主要体现为凸轮与推杆之间的摩擦, 因此, 忽略其他各处的摩擦) 。力矩M和载荷P的虚功可分别表达为

δWM=Mδφ (20)

δWP=-Pδh=-Pecos (φ+θn) (δφ+δθn) (21)

为了研究凸轮和推杆之间的摩擦力, 画出推杆的受力图 (图3) 。图3中, NB和 NE为滑道作用于推杆的约束力, F为回复弹簧作用于推杆的压力, NH 和 FH分别为凸轮作用于推杆的正压力和摩擦力 (这里设凸轮与推杆之间无润滑油膜, 即凸轮与推杆之间的摩擦视为干摩擦) 。

根据牛顿第二定律, 推杆的运动微分方程可写为

$m{}_{\text{f}}\ddot{h}={\rm N}{}_{{\rm H}}-F-{\rm P}$ (22)

根据虎克定律, 回复弹簧作用于推杆上的压力可以表达为

F=kλ=ke[1+sin (φ+θn) ] (23)

将式 (9) 和式 (23) 代入式 (22) 后, 可得NH的表达式:

$\begin{array}{l}{\rm N}{}_{{\rm H}}=m{}_{\text{f}}e[(\ddot{{\displaystyle \phi }}+\ddot{\theta }{}_n)\cos (\phi +\theta {}_n)-(\dot{\phi }+\dot{\theta }{}_n){}^2\sin (\phi +\\ \theta {}_n)]+ke[1+\sin (\phi +\theta {}_n)]+{\rm P}(24)\end{array}$

根据库仑摩擦定律, 凸轮作用于推杆上的摩擦力为

FH=μNH (25)

式中, μ为凸轮与推杆之间的摩擦因数。

根据牛顿第三定律, 进一步可以得到推杆作用于凸轮上的摩擦力F′H (图4) :

F′H=FH (26)

由图4易写出摩擦力F′H的虚功表达式:

$\begin{array}{l}\text{δ}W{}_{{\rm H}}=-F^{\prime }{}_{{\rm H}}h\text{δ}\text{ϕ}=-\mu \{m{}_{\text{f}}e[(\ddot{{\displaystyle \phi }}+\ddot{\theta }{}_n)\cos (\phi +\theta {}_n)-\\ (\dot{\phi }+\dot{\theta }{}_n){}^2\sin (\phi +\theta {}_n)]+ke[1+\sin (\phi +\theta {}_n)]+\\ {\rm P}\}[R+e\sin (\phi +\theta {}_n)](\text{δ}\phi +\text{δ}\theta {}_n)(27)\end{array}$

将式 (20) 、式 (21) 和式 (27) 相加, 可得到作用在整个机构上的所有非有势力的虚功之和的表达式, 即

$\text{δ}W=\text{δ}W{}_{{\rm M}}+\text{δ}W{}_{{\rm P}}+\text{δ}W{}_{{\rm H}}=Q{}_{\phi }\text{δ}\phi +{\displaystyle \sum _{i=1}^{n}Q}{}_{\theta {}_i}\text{δ}\theta {}_i(28)$

$\begin{array}{l}Q{}_{\phi }={\rm M}-{\rm P}e\cos (\phi +\theta {}_n)-\mu \{m{}_{\text{f}}e[(\ddot{{\displaystyle \phi }}+\\ \ddot{\theta }{}_n)\cos (\phi +\theta {}_n)-(\dot{\phi }+\dot{\theta }{}_n){}^2\sin (\phi +\theta {}_n)]+\\ ke[1+\sin (\phi +\theta {}_n)]+{\rm P}\}[R+e\sin (\phi +\theta {}_n)](29)\end{array}$

Qθi=0 i=1, 2, …, n-1 (30)

$\begin{array}{l}Q{}_{\theta {}_n}=-{\rm P}e\cos (\phi +\theta {}_n)-\mu \{m{}_{\text{f}}e[(\ddot{{\displaystyle \phi }}+\ddot{\theta }{}_n)\cos (\phi +\\ \theta {}_n)-(\dot{\phi }+\dot{\theta }{}_n){}^2\sin (\phi +\theta {}_n)]+\\ ke[1+\sin (\phi +\theta {}_n)]+{\rm P}\}[R+e\sin (\phi +\theta {}_n)](31)\end{array}$

式中, Qφ、Qθi分别为对应于广义坐标φ和θi的广义非有势力。

2 偏心凸轮机构的动力学方程

偏心凸轮机构的动力学方程的拉格朗日方程如下:

将式 (11) 、式 (19) 、式 (29) ～式 (31) 代入式 (32) , 再经符号运算后, 得到

式 (33) ～式 (35) 就是偏心凸轮机构的动力学方程 (动力学模型) , 需要指出该模型同时计入了凸轮与推杆之间的摩擦和凸轮轴柔性的影响。

3 算例

下面针对某偏心凸轮机构, 进行动力学的数值模拟, 设该机构的具体参数如下:Jd=0.6431kg·m2, ρ=7800kg/m3, l=0.3m, Ip=4.9701×10-9m4, G=80GPa, Jc=0.4652kg·m2, mc=27.5675kg, R=0.15m, e=0.075m, mf=9.8083kg, k=1500N/m, μ=0.14, M=46N·m, P=20N。机构的初始状态为φ (0) =θ1 (0) =θ2 (0) =θ3 (0) =0和$\dot{\phi }(0)=\dot{\theta }{}_1(0)=\dot{\theta }{}_2(0)=\dot{\theta }{}_3(0)=0$, 模拟中凸轮轴被等分为3个轴单元。

采用四阶Runge–Kutta法[18]求解式 (33) ～式 (35) , 得到该机构的运动响应。图5所示为驱动盘转角φ的时间历程;图6所示为凸轮相对驱动盘扭振的时间历程 (凸轮相对驱动盘的扭转角θ3的时间历程) ;图7所示为推杆位置h的时间历程。

由图5～图7可以看出:凸轮与推杆之间的摩擦对机构的运动具有一定的影响。相对不考虑摩擦的情况而言, 这种影响可以概括为凸轮与推杆之间的摩擦使得机构的运动变慢 (图5) , 使得凸轮相对驱动盘的扭振幅度变小 (图6) 并使得推杆往复运动的频率变低 (图7) 。

凸轮轴的柔性是导致凸轮相对驱动盘发生扭振 (图6) 的原因, 这种扭振又必然会对推杆的运动精度产生影响。为了显示这种影响, 特计算出推杆的位置误差Δh (推杆的实际位置与名义位置之差) 的时间历程。图8给出了推杆位置误差Δh的时间历程。从图8可以看出:推杆的实际运动偏离了名义运动, 或者说推杆的实际运动是在其名义运动的基础上叠加一了微小振动Δh。Δh在一定程度上影响了推杆工作的精确性。

4 结论

(1) 凸轮轴的柔性是导致凸轮相对驱动盘发生扭振的原因。这种扭振又进一步对推杆的运动精度产生影响, 使推杆的实际运动在一定程度上偏离其名义运动。

(2) 凸轮与推杆之间的摩擦对机构的运动具有一定的影响。相对于不考虑摩擦的情况而言, 这种影响表现为:凸轮与推杆之间的摩擦使机构的运动变慢, 使凸轮相对驱动盘的扭振幅度变小并使推杆往复运动的频率变低。

摘要：为了更加客观地研究偏心凸轮机构的动力学行为, 在考虑凸轮与推杆之间的摩擦以及凸轮轴柔性的基础上, 应用有限元法和拉格朗日方程建立了一种偏心凸轮机构的动力学方程。以该动力学方程为数学模型, 通过数值方法研究了该机构的运动响应。依据所获得到的运动响应, 分析了凸轮与推杆之间的摩擦以及凸轮轴的柔性对机构运动的影响。得出的结论说明, 在研究偏心凸轮机构的动力学行为时, 计入凸轮与推杆之间的摩擦以及凸轮轴柔性是必要的。

动力学计算 篇9

传染病长期危害着人类健康,它们的不断的流行给人类生存和发展带来了巨大的威胁。2003年出现的SARS病毒,传染性极强,具有相当的毒力,击倒了很多人,令全世界都感到恐慌,但是病毒及其活动规律同样受数学定律的支配[1]。长期以来,人类都在与传染病进行斗争,并建立了大量的数学模型用于研究和分析各种各样的传染病问题。人们通过对传染病病毒的长期接触和研究,对病毒已经相对了解,也有了一定的预防措施。

今天,计算机网络已经成为人们工作和生活中不可缺少的一部分,但是,计算机网络的普及和发展也给计算机病毒带来了前所未有的发展和进化,给人类的生活和发展带来了很大的负面影响和损失,网络蠕虫对计算机系统安全和网络安全的威胁日益增加。各色的蠕虫病毒的传播途径和复杂的应用环境使网络蠕虫的发生频率只增不降,潜伏性并没有减弱,而是越来越强,覆盖面也不断扩大,不断增加由此带来的损失。自1988年11月第一个互联网蠕虫莫里斯(Morris)出现以来,蠕虫就向人们展示了它不同于传统计算机病毒的特点和巨大破坏力。进入21世纪后,蠕虫更是频频大规模爆发,其爆发频率、扩散速度、造成的危害超过了以往任何时期,引起了网络安全研究人员的广泛关注[7]。Spafford首次对Morris蠕虫的功能结构和工作机制展开研究。关于传播模型,IBM的Kephart,White和Chess在1991年—1993年研究了病毒传播模型,根据他们的研究,邹长春等人以CodeRed为例,对以微分方程为基础的蠕虫传播模型进行了分析[6]。

在此背景下,本文简要介绍了传染病动力学,阐述了蠕虫病毒的危害及建立蠕虫传播模型的意义,列举了几种常用的传播模型。

1 传染病动力学及模型简介

1.1 传染病动力学概述

传染病动力学[5]是开展理论性定量研究的重要方法之一,它以种群生长的特性,疾病的发生及在种群内的传播、发展规律,以相关的社会因素等为依据,建立一种数学模型,这种模型能反映传染病动力学特性,我们要揭示疾病流行规律,预测它的变化趋势,找出其流行的重要原因,找到最合理的预防和控制策略,就要对模型动力学性态进行定性、定量分析和数值模拟,分析疾病的是如何发生的,为防制决策提供理论依据。动力学方法在从疾病的传播机理方面反映流行规律方面是个不错的方法,考虑到了流行过程中的全局性态。将传染病动力学与生物统计学以及计算机仿真等方法联系起来,进一步认识传染病流行规律,建立更加符合实际的模型与防制策略。传染病动力学模型不仅在传染病研究方面影响深远,而且也为生物种群分布、植病、网络病毒的传播、谣言的传播等自然和社会科学问题的研究提供可靠依据。

1.2 传染病动力学模型的基本形式

在传染病动力学中,主要沿用的由Kermack与McKendrick在1927年用动力学的方法建立了SIR传染病模型。SIR模型将总人口分成三类:易感者(susceptibles),其数量可用t时刻还没有染病但有可能被该类疾病传染的人数来表示,记为S(t),;染病者(infectives),其数量可用t时刻已被感染成为病人而且具有传染力的人数来表示,记为I(t);恢复者(recovered),其数量可用t时刻已不再是染病者的人数来表示,记为R(t)。

则有公式N(t)=S(t)+I(t)+R(t),其中N(t)为总人口。

建立SIR模型是以以下三个假设为基础的:

①忽略人口的出生、死亡、流动等种群动力因素。人口始终为常数,即N(t)≡K。

②当一个病人与易感者有过接触后,就一定有了一些传染力。假设t时刻单位时间内,一个病人能传染的易感者数目与此环境内易感者总数S(t)成正比,比例系数为β,那么我们可以用βS(t)I(t)表示在t时刻单位时间内被所有病人传染的人数。

③t时刻,单位时间内不再是染病者的人数与病人数量成正比,比例系数为γ,单位时间内不再是染病者数量为γI(t)。

在以上三个基本假设条件下,易感者从患病到恢复的过程框图表示如下:

SIR基础模型用微分方程组表示如下:

传染病的流行规律可以对SIR模型的分析和解的渐近性态为基础进行研究。如果忽略出生与死亡等种群动力学因素,传染病若无潜伏期,动力学模型为:SI模型,患病后不容易治愈;SIS模型,患病后能恢复,恢复者没有免疫力;SIR模型,患病者回复后获得终身免疫力;SIRS模型,病人康复后只在短时期内有免疫力,单位时间内将有一些治愈者因免疫力缺失而可能又被感染。如果不忽略传染病的潜伏期,在三类人群中再增加感染而未发病者(Exposed)一类,可得到比SIR或SIRS模型更复杂的SEIR或SEIRS模型。若考虑种群动力学、疫苗接种、年龄结构等更复杂的因素,模型的参数和复杂程度则与以往不同,都将加大[5]。

2 蠕虫病毒概述

2.1 蠕虫的定义

蠕虫是一种病毒,是一种能够通过网络传播的恶性病毒,具有传播性、隐蔽性、破坏性等等,同时也具有自己的一些特征,自身的主动性和独立性更加突出。

Kienzle和Elder从以下4个方面对蠕虫做出了解释,即网络传播、破坏性、主动攻击和独立性。蠕虫实际上是一种恶意代码,它能通过网络传播,不用用户干预就能够独立地或者依赖文件共享发起主动攻击。根据传播策略,他们把网络蠕虫分为Email蠕虫、传统蠕虫和文件共享蠕虫三类[6]。

在2003年10月的世界蠕虫会议上,Schechter和Michael D.Smith提出了一类叫做Access For Sale蠕虫的新型网络蠕虫。它不仅具有上述特征,还有身份认证的特点[6]。

因此,本文认为:蠕虫病毒是在一种可以独立运行,利用计算机网络传播,不需要计算机使用者的干预展开攻击,是一种智能化、自动化的病毒技术。

2.2 蠕虫的特点

蠕虫病毒造成的损失是相当大的,也是比较难清除的。总得来说,蠕虫具有以下特点:

①利用操作系统和应用程序的漏洞进行攻击。此类病毒有“红色代码”、“尼姆达”、“求职信”、“冲击波”和“震荡波”等。

②传播方式多样。如“尼姆达”和“求职信”等,这类病毒可把文件、电子邮件、Web服务器和网络共享等作为传播途径。

③此病毒的制作技术与传统病毒有差别,其实现依靠的是最新的编程语言和技术,不难修改并产生变种,反病毒软件难以查杀。可以潜伏在网页中,在浏览网页时触发。

④与黑客技术结合,潜在的威胁和可能造成的威胁更大。

2.3 蠕虫的威胁

蠕虫病毒是利用计算机网络进行传播的病毒,其传播速度远远快于早期病毒。早期的档案型病毒,要征服全世界的话,要花很长时间,甚至需要几年的时间才能完成。到宏病毒出现的时候,病毒环球之旅的时间缩短了,甚至是以“星期”为周期。1999年出现了首只通过Email传播的Melissa病毒,打破了病毒绕行世界的纪录也再度,病毒环球之旅的时间大大缩短了,甚至是以“天”为周期。2000年,I Love You网络病毒传播的速度比Melissa宏病毒要高很多,而且具有极强的破坏力。2001年的Nimda蠕虫病毒在区区1天的时间内感染了不少于220万台的计算机。据ZDNET报导,2003年8月的Sobig蠕虫病毒,就在24小时之内,美国在线(AOL)就收到多达1150万封携带Sobig.F的电子邮件[9]。

病毒发作频率非常大,差不多每个月都会出现新的病毒疫情,有一百多种计算机病毒造成的影响较大。求职信、梅丽莎、FunLove、网络神偷、爱虫、新年贺卡、红色代码、尼姆达、DDOS、Slammer等病毒形成了比较严重的影响,给网络造成了不小得损害。这些病毒都不约而同的把网络作为主要传播途径,病毒的传播速度不断加快,给网络造成的损害越来越严重[9]。现在的蠕虫病毒,只需要15分钟,就可造成全世界的病毒疫情几乎同步爆发,其扩散全球的速度见表1。

2.4 建立蠕虫病毒传播模型的意义

随着网络发展与应用,蠕虫病毒也随之发展,其破坏性也日益增大,造成了巨大的损失。建立并研究蠕虫病毒的传播模型,可以对其进行有效的预测,为控制病毒提供有效的决策。目前,蠕虫病毒成为网络病毒的主流,研究蠕虫病毒的传播模型更具有现实意义。

蠕虫病毒传播模型是一种检验理论和评价定量猜想的实验工具,数学化的过程加深了对蠕虫病毒传播特点的认识。研究蠕虫病毒传播模型,可以用于预测蠕虫病毒未来的发展和传播趋势;为研究蠕虫病毒提供参考,通过模型参数搜集所需要的重要信息与数据;可以理论评估各种预防措施和解决方案的效果。建立蠕虫病毒传播模型,能更好地反映蠕虫病毒的传播规律和全局性,加深对蠕虫病毒的深入认识,使所建立防治策略更加符合实际。

3 蠕虫病毒传播模型

根据参考文献和资料,本文在此列举出几种常见的蠕虫传播模型:简单传播模型、Kermack-Mckendrick模型、SIS(Susceptible→Infectious→Susceptible)模型、双因子模型和BCM模型[6]。

3.1 简单传播模型

在简单传播模型(Simple Epidemic Model,简称SEM)中,每台主机保持易感染的和感染的两种状态。他们的状态转变过程可描述为:易感染→被感染。SEM模型的微分方程表达式为:

在公式中,I(t)为时刻t已被感染的主机数;N为网络中主机总数;β为主机感染率。当t=0时,I(0)为已感染的主机数,N-I(0)为易感染主机数。

令a(t)=I(t)/N,根据上式得方程:

3.2 Kermack-Mckendrick模型

Kermack-Mckendrick传播模型(简称KM模型)的主机保持易感染、被感染和免疫三种状态。KM模型的微分方程表达式为:

在公式中,I(t)表示时刻t依旧有感染性的主机数;R(t)表示时刻t已经从被感染的机器中免疫的主机数;J(t)表示到时刻t全部被感染过的主机数,是仍具有感染性的和已经从被感染的机器中免疫的主机数之和,即J(t)=I(t)+R(t);γ是主机从被感染的机器中移除的恢复率;β是感染率;N表示网络中全部节点主机;S(t)表示时刻t仍具有脆弱性的主机数。

KM模型是以SEM为根据考虑感染主机免疫的状态,更适合蠕虫传播的情况。但是,该模型依旧忽略了易感染主机和感染主机被补丁升级或人为对抗蠕虫传播的情况。另外,也不应该把感染率看做常量。

3.3 SIS模型

在SIS模型中,每台主机也只有两种状态:易感染状态与已感染状态,但与简单传播模型所不同的是,它考虑了一台主机被感染后会被修复的可能,但被修复的主机对该蠕虫仍然是易感染的。它的微分方程形式如下:

式中,J(t)表示至时刻t被感染的主机数;β表示主机感染率;N表示网络中主机总数;δ表示已感染主机的修复率。

显然,不应该假设β与δ为常量,这是与实际不相符的。另外,该模型假设一台主机被修复后状态又变回到易感染状态,也是与实际大多数情况所不同的。

3.4 Two-Factor模型(双因素传播模型)双因素传播

模型参考的外界影响因素和蠕虫对抗措施及各ISP节点或用户的对抗措施较多;网络蠕虫的快速传播导致一些路由器不能畅通,从而使网络蠕虫的传播速度不高。在这个模型中,β(t)、R(t)和Q(t)三个参数都是随着时间t的变化而变化的,双因素传播模型的微分方程表达式为:

在公式中,R(t)表示时刻t感染后被免疫的主机数;Q(t)表示时刻t被感染前就作了免疫处理的主机数;I(t)表示具有感染性的主机数;β(t)表示时刻t的感染率;S(t)表示时刻t易感染的主机数;J(t)表示时刻t已被感染的主机数,J(t)=R(t)+I(t);γ、μ和β0为常量。我们可从上式中得出I(t)和时间的关系式:

双因素传播模型弥补了SEM和KM模型的不足,具有他们所没有的优点,比他们更适合网络蠕虫的传播状态。但双因素传播模型忽略了大规模自动补丁或升级对抗网络蠕虫传播的情况。此外,采用蠕虫对抗蠕虫使网络中蠕虫传播变得更为不再简单。

3.5 Worm-Anti-Worm模型(WAW模型)

该模型考虑网络中有蠕虫A(恶意蠕虫)和蠕虫B(对抗蠕虫)两类蠕虫。蠕虫A的传播可分为两个阶段。在蠕虫B还没出现的时候,蠕虫A的传播行为遵循双因素模型。当蠕虫B出现以后,网络中蠕虫A的传播分为4种情况:蠕虫B查杀蠕虫A并为感染主机修补漏洞;蠕虫B对所有的易感主机修补漏洞,并查杀蠕虫A;蠕虫B只查杀蠕虫A;蠕虫B对所有的易感主机修补漏洞。下面我们以情况一为例来分析蠕虫A的传播模型。根据双因素模型,从时刻t到时刻t+Δt,易感主机数S(t)的表达式为:

在公式中,对于蠕虫B来说,S(t)是t时刻的所有易感主机,且网络中主机只存在易感染的和已感染的两种状态。蠕虫B的传播行为应与SEM模型一致。感染主机的微分方程数据表达式为:

其中,RB(t)是在t时刻蠕虫B修复的主机。根据双因素模型式,WAW模型的表达式为:

WAW模型考虑到了对抗性蠕虫的存在,使蠕虫后期的传播趋势能够预测的更精确。但这种模型也有一些问题,忽略了对抗性蠕虫本身的传播与其他限制因素之间的联系等。

4 结束语

蠕虫作为一种计算机互联网发展下的一种新型病毒,再加上计算机网络缺乏中心控制能力以及蠕虫的自主传播特性,导致目前蠕虫病毒频频爆发,对计算机网络带来严重威胁。从目前蠕虫病毒的发展情况看,蠕虫入侵行为的自动化以及更多的黑客技术用到蠕虫病毒的编写中,蠕虫技术也会不断的完善和发展,使得检测和防范蠕虫病毒具有一定的难度。因此,建立并研究蠕虫病毒的传播模型,掌握蠕虫病毒传播的特性,对预测蠕虫病毒未来的发展和传播趋势,对蠕虫病毒的防治具有至关重要的作用。

对于蠕虫病毒的防治,本文认为要做到一下几个方面:

动力学计算 篇10

近来,基于图形处理器(GPU)的通用计算技术(GPGPU)得到了飞速的发展,图形处理器被越来越广泛地应用到图像处理以外的计算领域[1]。虽然,GPGPU的概念在很早就被提出,但因受限于编程复杂度和精度问题,而一直没有得到广泛的应用。随着CPU性能瓶颈的到来,以及现代GPU在可编程性上的不断提高,GPGPU开始得到越来越多的重视,主要的GPU厂商NVIDIA和AMD都推出了自己的GPGPU编程模型,并不断地提供改进。统一计算设备架构CUDA(Compute Unified Device Architecture)是由NVIDIA公司提供的面向GPU的并行通用计算构架,提供支持使用C语言进行开发的编程模型,降低了开发难度,并在计算能力1.3以上的GPU设备中加入了对双精度的支持。使其可被广泛应用于分子动力学、气象预报、市场分析、天体物理学、机器学习等领域。

分子动力学模拟指使用经典力学,通过分析系统内每个粒子的受力状态,求解系统内粒子位置和速度,确定粒子的运动状态,进而推知系统的性质。分子动力学模拟效果很大程度上取决于分子力场,现在广泛使用的分子力场大多使用原子中心固定点电荷来计算静电相互作用势能,忽略了非原子中心电荷的作用和分子随外势变化而引起的静电极化和电荷转移。文献[2,3]根据密度泛函理论和电负性均衡原理提出了原子键电负性方法(ABEEM)[4]计算大分子的电荷分布情况,并基于此提出ABEEMσπ模型。该模型能够更好地反映非原子中心电荷的作用和分子外势的影响,更加精确地描述大分子体系的电荷分布和极化作用,进而提供更恰当的静电作用势能描述,最终取得良好的分子动力学模拟效果。

本文基于MPI+CUDA混合并行环境,实现了应用ABBEMσπ模型计算大分子电荷分布。通过将占据大部分程序运行时间的解线性方程组运算移植到集群GPU上进行,借助GPU的并行计算能力加速整个程序。

1 ABEEM σπ模型

ABEEMσπ模型是依据密度泛函理论和电负性均衡原理提出的有机和生物大分子体系的电荷分布计算方法,有较好的精确度。相比于使用以原子为中心的固定电荷模型[3],ABEEMσπ模型充分考虑到分子中化学键的存在,它将分子划分为原子区域和化学键区域,并进一步考虑到化学键中的双键结构,将化学键双键区域划分为一个σ键区域和四个π键区域(每个双键原子各有两个π键区域),σ键电荷置于构成双键的两个原子之间共价半径之比处,π键区域电荷中心垂直于双键所在平面,位于双键原子上下两侧[4]。 最终将大分子划分成六种不同的区域:单键原子、双键原子、孤对电子lp、单化学键区域、双键σ键区域和双键π键区域。ABEEMσπ分子体系总能量表达式如下:

$\begin{array}{l}E{}_{ABEE{\rm M}}={\displaystyle \sum _{i=1}^{{\rm N}}\{{\displaystyle \sum _a[}E{}_{ia}^{*}-u{}_{ia}^{*}q{}_{ia}+\eta {}_{ia}^{*}q{}_{ia}^2]+}\\ {\displaystyle \sum _{j=1}^{{\rm N}}}[\frac{1}{2}{\displaystyle \sum _{a\in i}{\displaystyle \sum _{b\in j}k}}{}_{{\rm H}-bond}(R{}_{lp,{\rm H}})\frac{q{}_{ia}q{}_{jb}}{R{}_{ia,jb}}+\\ \frac{1}{2}k{\displaystyle \sum _{a\in i}{\displaystyle \sum _{b\in j(a\ne b)}\frac{q{}_{ia}q{}_{jb}}{R{}_{ia,jb}}}}]\}(1)\end{array}$

其中N表示整个体系中分子的个数,u${}_{ia}^{*}$表示分子i中任意一个区域a的价态化学势,qia表示分子i中区域a的电荷,η${}_{ia}^{*}$表示分子i中区域a的价态硬度,lp和H表示形成氢键的孤对电子和氢原子,k为矫正系数。

由式(1),根据电负性均衡原理可以得到体系内第i个分子的任意一个区域a的有效电负性,表达式如下:

$\begin{array}{l}x{}_{ia}=x{}_{ia}^{*}+2\eta {}_{ia}^{*}q{}_{ia}+\\ {\displaystyle \sum _{b}k}{}_{{\rm H}-bond}(R{}_{lp,{\rm H}})\frac{q{}_{ib}}{R{}_{ia,ib}}+k{\displaystyle \sum _{b(b\ne a)}\frac{q{}_{ib}}{R{}_{ia,ib}}}+\\ {\displaystyle \sum _{j=1(j\ne i)}^{{\rm N}}\left[{\displaystyle \sum _{b\in j}k}{}_{{\rm H}-bond}(R{}_{ia,jb})\frac{q{}_{jb}}{R{}_{ia,jb}}+k{\displaystyle \sum _b\frac{q{}_{jb}}{R{}_{ia,jb}}}\right]}(2)\end{array}$

其中x*ia=-u*ia,是分子i在区域a的价态电负性。

ABEEMσπ模型中,单个中性分子的总电荷可以表示为:

$q{}_i={\displaystyle \sum _{a}q}{}_{ia}=0$ (3)

即,分子i的总电荷等于它的所有区域电荷之和,成中性。

根据电负性均衡原理,当分子处于平衡结构时,各区域的电负性相等,等于分子电负性。即$x{}_{ia}=\overline{x}‚\overline{x}$为分子电负性。如果一个分子可以划分m个区域,则根据式(2)我们可以得到m个方程,根据式(3)可以再得到一个方程。这样就可以形成以区域电荷qia和分子电负性$\overline{x}$为未知数,以η${}_{ia}^{*}$,u${}_{ia}^{*}$参数的m+1阶线性方程组,并且其系数为正定矩阵。通过求解此线性方程组就可以得到分子各区域电荷qia,从而了解分子内的电荷分布状况,为下一步分子力场计算提供精确的数据支持。关于方程组中矫正系数k和参数的确定方法以及各变量的具体意义请参考文献[5]。

2 算法GPU加速

实验发现在使用ABEEM σπ模型计算大分子的电荷分布时,串行程序95%以上的时间用于解大型线性方程组。因为该分子模型的电荷分布方程组的系数矩阵是正定矩阵,所以原程序使用Cholesky分解进行方程组求解。将程序中的Cholesky分解部分运算移植到GPU上进行能够极大地提高整个算法的效率,借助CUDA可以实现算法到NVIDIA GPU的映射。

2.1 CUDA概述

CUDA作为NVIDIA推出的GPU通用计算模型,包含了从开发工具到底层驱动一直到硬件构架的一整套综合体系。NVIDIA已经在最新的GPU内核构架Fermi上将通常的流处理器SP(Streaming Processor)更名为CUDA核心(CUDA Croe), GPU正走出专业图像处理领域,更多进入到通用计算领域。相比于CPU明显的性能优势使得GPU在应用方也得到越来越广泛的重视,G200核心的浮点运算能力达到1TFlops以上,Fermi (G300)核心的浮点运算能力达到2.4TFlops,超过同期CPU十倍以上速率。

CUDA编程模型可以采用C语言进行开发,不需要借助图形学API。一个CUDA程序通常包含两部分:CPU的部分和GPU部分。通过CPU和GPU的协作完成整个任务的运算,这也是其名字Compute Unified Device Architecture的由来。GPU内部只有很少的逻辑控制元件,需要CPU来进行任务调度和一些串行方面的工作,GPU则主要用于并行数据运算。CPU通过调用核心函数来启动GPU线程组。GPU的线程组分为三级结构:grid、block、thread,每一个block(线程块),可以最多包含一定数量(G80核心521个、G300核心1536个)的thread(线程),而多个block组成一个grid。Block中thread可以采用一维、二维或三维的形式,通过threadid进行线程识别。Block也以同样的形式被组织在grid中,并通过blockid识别。数据也同样采用不同的维度形式存储到GPU存储器中,合理的数据存储形式可以满足数据对齐[6]的要求,合并线程访问,提高带宽利用率,CUDA程序中带宽的优化往往比计算更重要[7]。

2.2 Cholesky分解GPU移植

Cholesky分解是线性代数计算领域重要的基础运算,是求解系数为对称正定矩阵的线性方程组的主要方法。N阶正定矩阵A,可以分解成A=LLT的形式,L是下三角矩阵,LT为其转置矩阵。为避免分解过程中对矩阵A对角线元素的开方运算,通过对L和LT矩阵的对角线归一化处理,得到改进后的cholesky分解A=LDLT,公式如下:

算法基于集群并行计算系统,集群的各节点机采用CPU+GPU混合构架,并提供CUDA运行环境支持,节点之间采用MPI消息传递接口进行通信。N阶电荷分布线性方程组,其系数矩阵A,集群节点数为p。将准备好的系数矩阵A以行块卷帘的形式分割存储在各节点中,即将矩阵A的每m行划分为一个块,将每个块以卷帘的形式分发存储在各节点中,这种存储方式可以兼顾负载平衡度和通信开销。在本算法中考虑对齐GPU线程束warp[8](32个线程),m的大小最好选取为16×2n。

集群的所有节点,进行i=0,1,2,…,N-1循环。对于第i次循环,具体算法描述如下:

① 若系数矩阵A的第i行存储在本节点,则将其以一维形式拷贝至GPU全局存储器,并对调用GPU计算aij/aii,j>i。将计算结果拷贝回节点然后广播给其他节点。若第i行并未存储在本节点,则进入循环后接收i行所在节点发送来的数据,并将其拷贝至GPU。

② 将节点机中存储的行块拷贝至GPU全局存储器,拷贝依靠以下三条规则进行:一,块中数据的行号必须大于i,i行所在的块只拷贝其中大于i的部分。二,块中的数据只截取一部分拷贝至GPU存储器中,假设块中数据是矩阵A的第j行(j>i)到第j+m-1行,则截取块中j列之后的数据进行拷贝,如图1的阴影部分所示。三,数据以二维的形式拷贝至GPU存储器。

③ 每当完成一个块拷贝,CPU调用内核函数,启动GPU线程组对该数据块进行计算。线程组采用每个block含有256个thread,以16×16二维形式组织在block中,block同样以(m/16)×[(n-j)/16](产生余数时向上取整)形式组织在grid中,这样线程和块数据做到了一一对应,每个线程只负责处理块中的一个数据。

④ 将处理完的块从GPU存储器拷贝至节点内存,从节点将下一个符合规则的块拷贝至GPU进行计算,直至所有块完毕,完成本次循环。

当所有循环完毕就完成了分解运算。

2.3 GPU移植算法优化

上述算法,数据在异构系统中的计算和拷贝是依次串行进行的,GPU设备在数据计算和数据传输之间相互等待,而CPU在这个过程中一直处于空闲状态。这使通信开销在算法中所占比重过大,浪费了部分计算资源,影响算法性能。为此,采用异步并发的方法,对算法进行优化。即通过对任务合理划分,使异构平台各设备,能在同一个时间段分别处理各自的任务,从而提高资源利用率,提高算法性能。任务能够并发的基本要求就是并发各任务不存在直接相关性,算法优化分为两个方面:

(1) 主机和设备之间的异步并发。其中包括主机和设备内核执行的并发、主机和数据传输的并发。CPU在调用设备内核函数或带有Async后缀的CUDA API存储器拷贝函数后,在完成任务之前,就会将线程控制权返回给CPU,此时可以使用CPU进行其他任务。在此,将原算法中小部分计算任务分配到CPU进行,来提高资源利用效率。我们将原算法中①步骤分配到CPU进行,即,对于i次循环,CPU负责计算系数矩阵A的第i+1行的运算a(i+1)j/a(i+1)(i+1), GPU负责其他部分的计算,CPU的计算将与GPU的计算以及数据传输并发进行。

(2) 设备内核函数的执行和数据传输的并发。CUDA支持数据传输和内核执行的并发(需要计算能力1.1以上的设备),即GPU在处于计算状态的同时,也可以进行与正在进行的计算不相关的数据在设备存储器和主机存储器之间的传输。将原算法中②、③、④步骤改为并发进行,具体优化方法为:在设备存储器中开辟同大小的缓存区1、缓存区2,在第i次循环,对节点机中的等待GPU处理的数据块,将一个数据块传输到设备存储器缓存区1中,传输完毕后,调用GPU内核函数对缓存1中的数据进行分解,此时GPU处于忙状态,控制权返回到CPU线程,调用Async后缀的拷贝函数将下一个待分解数据块发送至缓存2中,系统总线PCI-E进入忙状态,此时数据传输和内核执行重叠。当数据传输与GPU执行的完毕后,调用GPU内核函数对缓存2中的数据进行分解,并将缓存1中计算完毕的数据拷贝至节点机内存中,然后将下一个待分解数据块发送至缓存1中,此时再次进入并发状态。重复上面的流程,直到本次循环待分解数据处理完毕。

通过优化,寻找到不存在直接相关性的任务,将其同时并发到系统的各个设备,使CPU、GPU、PCI-E总线这三个设备在某些时间段同时处于工作状态,各自独立完成部分任务。提高利用效率,减少系统资源浪费。

3 性能测试

实验环境采用由5台通用PC构建的集群系统。配置为主频2.4GHz的Intel Xeon CPU,4GB DDR3内存,GeForce GT 260显卡,主机与显卡直接通过PCI-E2.0X16总线连接。集群采用Linux RH9系统,安装CUDA Toolkit2.1与MPI并行编译环境。

实验采用了分子系统电荷区域总数分别为1728、3956、8497的三个模型,对模型在不同节点数下,MPI程序MPI、MPI+CUDA程序M_CU、优化后的MPI+CUDA程序M_CU_O计算时间进行测试,结果如表1所示。

从表1中,可以发现GPU的加速效果随着模型的增大而越来越明显,当模型电荷区域总数很小时,几乎没有加速效果,当模型达到8497时单节点的加速比可以达到2倍以上,随着节点数目的增加,GPU的加速效果也逐步降低。说明GPU在进行大数据量单指令流计算时有很高的性能优势,当计算量减小后,通信开销对性能影响变大。优化后的算法相比于未优化的算法节省了25%-30%左右的时间。通过优化,将一部分任务交由CPU完成,提高了CPU的利用率,减少了由这部分任务引发的PCI-E通信开销。并且将GPU计算和数据传输重叠,提高了GPU和系统总线PCI-E的利用率,尽量减少了各部分资源的相互等待时间,提高了算法效率。

4 结 语

本文实现了基于MPI+CUDA异构并行计算平台的ABEEMσπ模型分子电荷分布算法,针对ABEEMσπ模型电荷分布计算的特点,通过将主要的计算部分移植到GPU上进行,利用GPU对计算加速。基于异构平台的特点,采用异步并发的方法对算法进行优化。通过测试对比证明了移植优化后的算法是高效可靠的,在大分子电荷分布计算方面有大的性能提升,为下一步分子动力学模拟中静电力计算提供了精确的电荷分布模型。

摘要：在集群与GPU组成的异构并行计算平台上,使用MPI+CUDA混合编程模型,实现基于ABEEMσπ模型的分子动力学模拟中电荷分布的计算。通过对电荷分布分布求解中的计算部分移植到GPU上进行,并针对算法中通信开销大和资源未充分利用的问题,通过异构平台的异步并发方法进行优化,提高了求解效率。性能测试结果表明,相比于单纯MPI并行算法,优化后GPU加速的异构并行算法,在化学大分子模型电荷分布计算上,有着明显的性能优势。

关键词：GPU通用计算,集群系统,分子动力学,电荷分布,ABEEMσπ模型

参考文献

[1]张浩,李利军,林岚.GPU的通用计算应用研究[J].计算机与数字工程,2005,33(12):60-62.

[2]Geerlings De Proft,Langenaeke.Conceptual density functional theory[J].Chem Rev,2003,103:1793-1873.

[3]Wilfried J Mortier,Swapan K Ghosh,Shankar S.Electronegativity equali-zation method for the calculation of atomic charges in molecules[J].J Am Chem Soc,1986,108:4315-4320.

[4]杨忠志,丛尧,王长生.原子-键电负性均衡方法中的σπ模型及应用[J].高等学校化学学报,1999,20(11):1781-1783.

[5]Yang Z Z,Zhang Qiang.Study of Peptide Conformation in Terms of the ABEEM/MM Method[J].J Comput Chem,2005,27:1-10.

[6]NVIDIA.NVIDIA CUDA Programming Guide,version4.0[R].2011.

[7]Garland M,Grand S L,Nickolls J,et al.Parallel Computing Experiences-With CUDA Micro[J].IEEE,2008,28(4):13-27.

动力学计算 篇11

【摘要】计算结构力学是在结构力学基础上为研究生开设的一门重要基础课程, 将为一系列后续课程的学习打下重要基础。本文从建立立体式课程教学体系、创新教学内容、创新教学方法、以科研促进教学等措施方面进行了教学改革尝试。

【关键词】研究生教学;教学改革;高等教育

【中图分类号】G420 【文献标识码】A 【文章编号】1001-4128(2009)10-0061-03

根据培养跨世纪人才的需要,高校的课程教学必须同时代的发展与科技的进步相适应,在内容、体系及观念上有所更新和突破[1-2]。目前,研究生课堂教学面临着众多的困难与挑战。从1999年持续至今的研究生扩招的背景下,研究生的招生规模和人数成倍增长,这给研究生的课堂教学提出了更严峻的考验。研究生教学所处的困境,以及由此所带来的研究生教学中存在的问题,一方面受传统研究生办学模式的影响,另一方面也由研究生教育发展过程中产生的诸多问题所致。改善研究生课堂教学是提高研究生培养质量的必然,也体现出研究生教学改革的紧迫性[3-4]。计算结构力学属于一门用于工程结构计算的技术基础课,是联系基础理论和专业知识的桥梁,在有关工程专业中占有重要地位。课程的这一性质决定了它具有理论性与实践性的特点,其教学不应只注重理论,偏离实践,而应二者兼顾。传统的计算结构力学课程有失偏颇,只侧重基本理论及其解题方法、技巧,比较忽视基本理论在工程实际中的应用,因而存在不完整性和缺乏先进性的问题,不利于学生能力的培养和提高,也不符合现代教育的要求。为此,根据当前形势发展的需要, 为了提高教学效果、教学水平,加快该课程的改革和建设步伐, 我们从2005年开始着手改革计算结构力学课程,借鉴了相关力学课程的改革经验,积极探索新方式,在计算结构力学的教学体系完善、教学模式、师资队伍建设等都尝试一些创新的方法,并取得了较好的效果。

1紧密结合专业建设,建设立体式课程教学体系

工程技术系统日趋增加的复杂性,更需要多学科的技术支撑,为此建立协调完善、融会贯通的课程体系十分必要。我们根据结构力学,计算数学,计算机应用学科在计算结构力学主要环节的内容交叉情况,从加强基础的角度,设置课程模块并进行程序设计,较好的促进了学科内容的相互融合渗透,知识结构系统完整,模块化教学循序渐进,建立了一个完整协调的立体化课程教学新体系[5]。计算结构力学立体化教学系统主要由三个模块组成,包括传统的纸质课本教材,一个教学系统和一个学习系统(见图1),它们分别对应的是客体,环境和主体三个教学要素。教学系统主要是教师授课时使用的多媒体教学课件模块,学习系统主要是用于学生课后复习以及网上提交作业的课件模块。三个模块之间协同作用,互相补充,共同完成教学任务,其中,纸质教材采用自主编写的面向船舶与海洋工程专业研究生教学的新教材。教学系统和学习系统又可相应分为几个子系统,各子系统包括若干功能模块,见图2。该课程体系还较好的体现了学科内容与交叉。近几年的教学实践表明,这一新的教学体系不仅精炼了教学内容,更可通过系统知识的模块化程序设计充分激发学生的主动进取精神和创新意识。极大的提高了学生的实践能力和计算机应用能力。

2创新教学内容,反映国际前沿学术动态

面对科学技术的飞速发展,现有的教学内容要不断更新,才能紧跟时代的前进步伐,使我们所传授的知识具有新鲜性。因此为了在教学中引进反映与计算结构力学相关领域的科技成果的新内容、新理论和新技术,提高教学内容的起点,我们在计算结构力学整个授课过程中开设两次反映当前与该领域相关的最新学术进展的科技讲座。讲座的内容是根据相关领域的国内、国际权威会议上的最新交流成果摘集而成。这些讲座不仅给学生开阔了视野,而且对其今后的研究方向也有着重要的影响,对培养高水平、高素质的人才有重要意义。

3创新教学方法、教学模式,提高对学生综合能力的培养

在体系改革的同时,教学方法的改革也要跟上去,这两者是相辅相成的。教为主导,学为主体是教学的原则与规律,但教师教什么,让学生学什么却必须认真研究。我们在计算结构力学课程的教学过程中本着使学生获得必要的知识的基础上,更注重使学生获得汲取知识的能力,发现问题、解决问题的能力以及独立思维和自我判断的能力的培养,因此我们在教学方法、教学模式上作了以下创新尝试:

3.1注重教师的主导作用,培养学生科学的思维方法

在讲课过程中注意围绕重点、难点精讲多练,对讲授内容注重讲概念,讲思路,讲方法,讲关键。例如在讲解平面问题和空间问题的有限元法中给学生理清有限元方法解决问题的基本思路,即结构离散化、形成结构刚度矩阵、计算等效节点荷载、列方程、解方程这一条主线。使学生在学习本课程的其他内容时,清楚这些内容具体是为了解决上述那一部分问题,从而培养学生科学的思维方法,系统的学习方法。

3.2引用多媒体、数值仿真实验室的动态教学方法,提高学生的理论联系实际能力

为了使学生适应新的需求,在不影响基本理论学习的情况下,更新计算手段,引用多媒体教学之后,可以适当延伸一些课程内容,通过直观、动态的多媒体课件,使学生有着很好的感性认识,之后具体讲解这些过程是应用了那些基本原理、基本方法得到的,调动同学的学习热情。此外,在结构运动学和动力学部分、结构的非线性部分、结构稳定性问题、断裂力学问题、流固耦合问题等章节中根据理论联系实际的原则,分别采用典型工程结构问题建立了数值仿真实验室系统,该系统主要包括典型薄板弯曲变形数值仿真、船舶结构强度分析数值仿真、冲击荷载作用下船舶动响应仿真、波浪荷载作用下海洋平台动响应仿真。

3.3采用双语教学模式,锻炼学生的外语使用能力

以教材内容为依据,制作适用于双语教学的多媒体课件。在实际授课过程中,设置模拟环境,按照一定的方式来进行的一种互动式授课及实训,使学生在亲自参与的过程中,提高对专业知识理解,并锻炼学生对语言的恰当使用的能力。这种互动式授课可以增强学生灵活运用英语的能力, 同时也加深了对专业知识的理解和感悟。

3.4建设网络教学系统,加强教师和学生信息交流

网络教学由于覆盖面广,可以实现更大范围的信息资源共享。因此,为了迎合现代教育观念的先进的教学和学习方式、增强学生与学生之间,教师与学生之间的协作交流,我们建设了计算结构力学教学网站。该网站以沟通、交流为宗旨,设有教学区和学习区。计算结构力学网站的建立极大的激发了同学对该课程的学习热情,根据反馈信息,采用网络这一新颖的辅助教学措施,不仅可以促进学生之间的协作交流,使学生对学习的内容会理解得更深刻,学习思路更开阔,而且可以使教师通过学习者提出的问题和要求中产生一些新的教学思路,并反馈到教学设计当中。

4结束语

学生在校期间获取一定量的知识是非常重要的,知识是一个人素质的一部份。但教师不应只满足于教学生会做几道题,尤其对于研究生而言,更应注重教学生学会学习,学会科学的思维。素质教育应着眼于学生的主动性,独立性,创造性的发展与提高。在研究生阶段,教师应多给学生提供参与活动,表现自我的机会,使学生在成功的体验中不断进取、开拓、创新,同时提高学生的自我意识,自身调控能力,鼓励学生大胆创新,求异思维,标新立异。以上是四年来对我们对计算结构力学课程进行改革的一些尝试。提供一些看法和做法与同行们交流探讨,以不断推动教学改革,完善计算结构力学研究生教学。

参考文献

1张爱华,高等教学质量观及高教改革论 [J]. 辽宁师专学报,2004(4):111-112.

2姚秋杰,关于高教改革的若干思考[J].北华大学学报,2004(2):2-5.

3张林,杨心德.浅析当前研究生教学存在的问题与对策[J]. 科技信息,2008(4):17-18.

4陶玉萍,纪军, 研究生教学模式探究[J]. 高层次人才培养技术与创新管理, 2007(3):64-66.

动力学计算 篇12

国家农机购置补贴政策的实施, 极大地促进了农业装备现代化, 提高了农田作业效率;同时对发展精准农业提出了要求, 以促进农田作业的精细化管理, 合理调配农业资源。目前, 国内精准农业机械化处于建立卫星定位基准站, 采集各类农田信息, 逐步实现精准播种、精准施肥、精准灌溉、精准动态管理和精准收获等阶段。这些“精准”字眼的背后显示发展旱地农业机械大型化、自动化、智能化和信息化的迫切要求。就自走式农业机械而言, 转向系统直接影响司机的劳动强度、整机行驶的速度与安全性, 进而影响着农业机械的作业效率与作业精度。随着国外重型拖拉机先进产品的进入和国内拖拉机市场与用户的成熟, 用户对拖拉机的驾驶舒适度、产品可靠性和操控性能的要求越来越高, 因而对轮式拖拉机进行整机性能的研究与开发提出了更高的要求。

1计算多体动力学概念

多体系统是指由大范围相对运动的多个物体构成的系统。计算多体系统动力学是一门涉及多体系统动力学、计算方法与软件工程的交叉学科, 是一般力学面向工程实践的新学科, 是机械系统虚拟样机技术的重要组成部分。其根本目的是最大限度地提高计算机对复杂系统运动学、动力学和控制系统的分析与综合能力。计算多体系统动力学的应用是指在机械领域对多体系统进行数学建模、仿真分析、检测与控制的统称。由于应用领域广泛和功能强大的商业软件的支持, 使得机械系统的虚拟样机技术的广泛应用成为可能。传统的拖拉机研发过程包括整机性能参数的选取、总体方案设计、论证与评审、图纸设计、物理样机实现、物理样机检测与参数对比、可靠性考核、改进与优化设计, 直至达到设计要求。在产品实现过程中存在大量的工作反复, 需要大量的时间、资金和人力。在研发过程中加入虚拟样机技术后, 增加了数学模型的修正环节, 虽不能减少研发流程, 但可以在图纸发布前利用合理模型解决后续可能出现的问题, 易于形成产品的系列化, 从而缩短了产品进入市场的周期, 减少了成本。

2计算多体动力学应用

1990年英国的拖拉机研究人员采用横向和纵向以Maxwell流变模型、垂向以Voigt-Kelvin流变模型模拟拖拉机轮胎作用, 分别建立11、12和14自由度共10种不同形式拖拉机模型, 以不同的行驶速度通过硬质标准路面, 利用数值仿真的方法研究拖拉机动态特性, 与物理实验结果对比, 振动均方根值的误差接近30%。在每一数值的仿真中, 选取轮胎侧向、纵向刚度与阻尼为常值, 未考虑变动的轮胎载荷、动态的接地印痕对轮胎的刚度与阻尼的影响。

2000年, 美国堪萨斯州立大学的研究人员采用负载变形、非滚动垂向自由振动、非滚动平衡负载变形、滚动垂向自由振动和滚动平衡负载变形五种方式, 利用特殊制造的铰接测试框架对不同充气压力下的轮胎的径向刚度进行测试, 并研究了轮胎振动的滞后现象, 得出了几点结论:

(1) 非滚动垂向自由振动测试的结果在所有充气压力下轮胎刚度最大。

(2) 负载变形、非滚动自由振动、非滚动平衡负载变形法得到的结果在所有充气压力下轮胎刚度近似。

(3) 滚动垂向自由振动和滚动平衡负载变形方法的结果在充气压力为83 k Pa和124 k Pa时接近。

(4) 当用并联弹簧与粘性阻尼表示轮胎简化模型时, 非滚动自由振动的结果对轮胎径向刚度的表达是不充分的。

(5) 负载变形、非滚动自由振动测试结果的差异是因为轮胎的滞后现象。

(6) 当用并联弹簧与粘性阻尼模型表示轮胎的简化时, 轮胎的特性应在期望的前进速度下测试。

(7) 在低充气压力下, 当用并联弹簧与粘性阻尼模型表示轮胎的简化时, 滚动自由振动对径向刚度的表达是不充分的。

2001年, 瑞典农业大学农业工程系的研究人员利用线性模型的方法仿真前后轴悬置的拖拉机在两种标准跑道上行驶, 研究了各线性、非线性悬架参数及拖拉机结构参数对悬架工作行程及拖拉机平顺性的影响。但模型只考虑了车身的垂向、俯仰、滚翻和前后轴的垂向和滚翻自由度。

2010年, 澳大利亚研究人员利用Cosmos Motion2005软件建立拖拉机与振动深松犁的仿真模型, 对农田作业中驾驶室座椅的振动进行了研究, 对比分析力仿真与实验结果, 验证了仿真模型, 为振动深松犁的结构改进与市场化提供了理论模型。

国内有关农业机械计算多体动力学技术的研究工作成果公开的还不多, 北农大研究人员用计算多体动力学方法建立了传统拖拉机牵引机组的22自由度动力学开环树形模型, 进行了基于农田不平度的拖拉机机组输入谱推导与振动仿真的研究, 未考虑特殊工况对拖拉机动态特性的影响及结构的要求。研究人员分别以装载机、集材拖拉机为例对铰接摆动桥形式的开环树形结构行走机械进行越障、模态和稳定性分析与研究, 采用车架与后桥2个质量刚体, 没有对前后车架的铰接约束及连接前后车架的转向油缸的动态特性进行描述。2009年山东理工大学的研究人员以非线性自适应随机接地印迹轮胎力学模型为计算前提, 用计算多体动力学方法建立了拖拉机牵引机组的22自由度动力学模型, 进行了基于农田不平度的拖拉机机组输入谱推导与振动仿真的研究。用RMS高度、自相关长度和自相关函数三个主要参数反映土壤表面不平度特征, 推导了拖拉机组的六轮输入谱矩阵, 研究了机组受农田不平度激励所产生振动的频域特性。

3总结