基于三阶Adams格式的求解声波方程的多步算法
第三周 Adams操作讲解
创建体6
选择，并单击坐标点spring、rotate_center即可 重命名body6
2、添加约束
弹簧阻尼器
创建作用点(body6上与坐标点rotate_center重合的点)
选择，单击body6,并单击坐标点rotate_center即可 重命名spring2
半径(Radius)，深度(Depth)
要求圆角半径变化时，需选择 End Radius 项，输入末端半径值， 此时圆角半径为起始半径
半径和深度的设置，不能使几何 体分成两块
凸圆 Boss
半Байду номын сангаас(Radius)，高度(Height)
挖空 Hollow 外围添 加材料
壳体厚度 t（Thickness）
添加部分厚度 t（Thickness） 删除内部的形体
此功能可以将经过合并、相交和切割等处
Split
理形成的组合形体，还原为原先的基本形 体。
注：实体 1—操作过程中的第一个实体，实体 2—操作过程中的第二个实体。
几何模型细节
功能 倒角 Chamfer
图标
参数设置 斜面的宽度(Width)
说明
圆角 Fillet
开孔 Hole
圆角半径(Radius) 末端半径值(End Radius)
1、建模(2)
创建体3：选择，并单
击点point5、point5_2, 则体3创建，重命名 body3
1、建模(3)
创建体4
复制体2（在edit菜单下，copy命令），使复制体绕坐标点 rotate_center 的Z轴旋转-120度
利用ADAMS进行动态仿真分析的一般方法和过程111
基于ADAMS软件进行动态仿真分析的一般方法和过程摘要：本文通过对相关资料的总结归纳，介绍了虚拟样机的发展现况、ADAMS软件、特点以及利用其进行动态仿真的一般方法和过程。
并结合多功能开沟机液压系统进行了建模与仿真分析。
关键词：仿真 ADAMS 优化虚拟样机1、前言随着近代科学技术的发展，工程设计的理论、方法和手段都发生了很大的变化。
从计算机辅助工程(CAE)的广泛应用，到并行工程(CE)思想的提出与推行，从根本上改变了传统的设计方法，极大地促进了制造业的发展和革命。
但与此同时，人们已清楚地认识到：即使系统中的每个零部件都是经过优化的，也不能保证整个系统的性能是良好的，即系统级的优化绝不是系统中各部件优化的简单叠加。
于是，由CAX/DFX等技术发展而来，以系统建模、仿真技术为核心的虚拟样机技术(Virtual Prototyping)得到了迅速发展，并正成为各国纷纷研究的新的热点。
虚拟样机技术(Virtual Prototyping Technology)是当前设计制造领域的一项新技术，其应用涉及到汽车制造、工程机械、航空航天、造船、航海、机械电子、通用机械等众多领域。
它利用计算机软件建立机械系统的三维实体模型和运动学及动力学模型，分析和评估机械系统的性能，从而为机械产品的设计和制造提供依据。
虚拟样机技术可使产品设计人员在各种虚拟环境中真实地模拟产品整体的运动及受力情况，快速分析多种设计方案，进行物理样机而言难以进行或根本无法进行的试验，直到获得系统的最佳设计方案为止。
虚拟样机技术的应用贯穿着整个设计过程中，它可以用在概念设计和方案论证中，设计者可以把自己的经验与想象结合在虚拟样机里，让想象力和创造力得到充分地发挥。
用虚拟样机替代物理样机，不但可以缩短开发周期而且设计效率也得到了很大的提高。
本文以ADAMS为平台，简单说明一下进行虚拟样机的动态仿真分析的一般方法和过程。
2、ADAMS软件简介及特点ADAMS(Automatic Dynamic Analysis of Mechanical System)软件，是由美国机械动力公司(Mechanical Dynamics Inc，现已经并入美国MSC公司)开发的最优秀的机械系统动态仿真软件，是目前世界上最具权威性的，使用范围最广的机械系统动力学分析软件，在全球占有率最高。
工程案例—机器人Adams虚拟实验详细步骤
一.ADAMS软件简介虚拟样机仿真分析软件ADAMS（Automatic Dynamic Analysis of Mechanical Systems）是对机械系统的运动学与动力学进行仿真的商用软件，由美国MDI （Mechnical Dynamics Inc.）开发，在经历了12个版本后，被美国MSC公司收购。
ADAMS集建模、计算和后处理于一体，ADAMS有许多个模块组成，基本模块是View模块和Postprocess模块，通常的机械系统都可以用这两个模块来完成，另外在ADAMS中还针对专业领域而单独开发的一些专用模块和嵌入模块，例如专业模块包括汽车模块ADAMS/Car、发动机模块ADAMS/Engine、火车模块ADAMS/Rail、飞机模块ADAMS/Aircraft等；嵌入模块如振动模块ADAMS/Vibration、耐久性模块ADAMS/Durability、液压模块ADAMS/Hydraulic、控制模块ADAMS/Control和柔性体模块ADAMS/AutoFlex等[3]。
1.1ADAMS软件概述ADAMS是以计算多体系统动力学（Computational Dynamics of Multibody Systems）为基础，包含多个专业模块和专业领域的虚拟样机开发系统软件，利用它可以建立复杂机械系统的运动学和动力学模型，其模型可以是刚体的，也可以是柔性体，以及刚柔混合体模型。
如果在产品的概念设计阶段就采取ADAMS 进行辅助分析，就可以在建造真实的物理样机之前，对产品进行各种性能测试，达到缩短开发周期、降低开发成本的目的。
ADAMS，即机械系统动力学自动分析（Automatic Dynamic Analysis of Mechanical Systems）该软件是美国MDI公司（Mechnical Dynamics Inc.）开发的虚拟样机分析软件。
目前，ADAMS已经被全世界各行各业的数百家主要制造商采用。
adams_solver(ADAMS求解器培训资料)
ADAMS/Solver中的四种文件类型ADAMS/Solver与ADAMS/View间的关系.bin/.cmd/.adm三种文件的异同点ADAMS/Solver 中的dataset file (.adm)练习: 读取.adm文件模型文件：shell_contact /shell contactGearGenerator2003_demo_shell.admmo shell.admADAMS/Solver command file练习: 读懂.acf文件模型文件：shell_contact /shell contact GearGenerator2003_de mo_shell.acfmo shell.acfSimulating a model in ADAMS/Solveradams05r2 ru‐sbelt_2d_pt_circle_hht.acf exhht acf练习: 读懂.acf文件运行模型模型文件：shell_contact /shell contactGearGenerator2003_demo_shell.acfmo shell.acf练习：批处理运行set MDI_ADAMS_CPU_COLUMN="Y"call adams05r2 ru‐s belt_2d_pt_circle_hht.acf excall adams05r2 ru‐s belt_2d_pt_circle_si2.acf exbelt2d pt circle si2acf模型文件：HHT / run_models.bat模型文件run models batADAMS/Solver Subroutines用户子程序流程总结P1. 编写程序代码user.c；P2. 在VC中编译user.c，生成user.obj；P3. 使用adams命令行，由user.obj生成user.dll；d命令行由bj dllP4. 将user.dll加入到adams库中；P5. 在Function Builder中即可使用自定义的函数。
第4章ADAMS软件算法基本理论-(陈立平)机械系统动力学分析及ADAMS应用
第4章ADAMS软件基本算法本章主要介绍ADAMS软件的基本算法，包括ADAMS建模中的一些基本概念、运动学分析算法、动力学分析算法、静力学分析及线性化分析算法以及ADAMS软件积分器介绍。
通过本章的学习可以对ADAMS软件的基本算法有较深入的了解，为今后合理选择积分器进行仿真分析提供理论基础，为更好地使用ADAMS打下良好的理论基础。
4.1 ADAMS建模基础ADAMS利用带拉格朗日乘子的第一类拉格朗日方程导出――最大数量坐标的微分－代数方程（DAE）。
它选取系统内每个刚体质心在惯性参考系中的三个直角坐标和确定刚体方位的三个欧拉角作为笛卡尔广义坐标，用带乘子的拉格朗日第一类方程处理具有多余坐标的完整约束系统或非完整约束系统，导出以笛卡尔广义坐标为变量的动力学方程。
4.1.1 参考标架在计算系统中构件的速度和加速度时，需要指定参考标架，作为该构件速度和加速度的参考坐标系。
在机械系统的运动分析过程中，有两种类型的参考标架——地面参考标架和构件参考标架。
地面参考标架是一个惯性参考系，它固定在一个“绝对静止”的空间中。
通过地面参考标架建立机械系统的“绝对静止”参考体系，属于地面标架上的任何一点的速度和加速度均为零。
对于大多数问题，可以将地球近似为惯性参考标架，虽然地球是绕着太阳旋转而且地球还有自转。
对于每一个刚性体都有一个与之固定的参考标架，称为构件参考标架，刚性体上的各点相对于该构件参考标架是静止的。
4.1.2 坐标系的选择机械系统的坐标系广泛采用直角坐标系，常用的笛卡尔坐标系就是一个采用右手规则的直角坐标系。
运动学和动力学的所有矢量均可以用沿3个单位坐标矢量的分量来表示。
坐标系可以固定在一个参考标架上，也可以相对于参考框架而运动。
合理地设置坐标系可以简化机械系统的运动分析。
在机械系统运动分析过程中，经常使用3种坐标系：（1）地面坐标系（Ground Coordinate System）。
地面坐标系又称为静坐标系，是固定在地面标架上的坐标系。
ADAMS多体系统建模与动力学方程迭代求解
动力学子程序，计算出分割面上的作用力和作用力 矩，然后将这些作用力和作用力矩做为反作用施加 到主子系统上，计算主子系统在同一积分时刻的各 分割面上的位移、速度和加速度，如此反复迭代， 直到满足收敛条件。此辅子系统的动力学分析子程 序 是 做 为 ADAMS 主 模 型 中 用 于 计 算 作 用 力 的 子 程 序出现的，因此辅子系统的动力学模型就被嵌入到 ADAMS 主 模 型 中 ， 体 现 了 系 统 的 整 体 求 解 。
& + ∑ c i& &i Fx = M 0& x τ
i=1 N
(3) (4)
用上述方法对图 1 所示的刚柔耦合系统进行迭 M＝ 374.448kg ， M0 ＝ 78kg ， 代求解， 各参数值如下： L ＝ 2.6m ， 悬 臂 梁 断 面 为 圆 环 ， 内 径 d ＝ 0.0355m ， 外 径 D ＝ 0.0785m ，弹 性 模 量 E ＝ 2.1× 10 11Pa ，密 度 ρ ＝ 7800kg／ m3， 推 力 F ＝ 370N 。 用 有 限 元 法 计 算 悬臂梁的固有模态时， 将 梁 划 分 为 20 个 等 长 三 维 梁 单 元 ， 截 取 前 10 阶 正 则 模 态 计 算 。ADAMS 求 解 器 为 WSTIFF 积 分 器 ， 最 大 积 分 步 长 为 0.0002s ， 最 小 积 分 步 长 为 0.00019s 。 图 4 是 ADAMS 迭 代 计 算 出 的 刚 体 A 的 加 速 度 曲 线 ，图 5 是 用 自 编 程 序 计 算 出的刚体 A 的加速度曲线。由图可见，计算结果一 样 。 在 ADAMS 迭 代 中 ， 求 解 器 预 估 校 正 （ 迭 代 ） 次 数 最 多 为 23 次 ， 最 少 为 8 次 。
ADAMS动力学仿真算法及参数设置分析
3 ADAMS 动力学仿真的参数设置
通过对 ADAMS 动力学仿真算法的分析可以
— 28 —
・ φ ( q , u , t) = 0 φ( q , t) = 0 ・ ・ η =λζ 其中 : , =μ 解决指数 3 公式计算问题的方法是所谓的指数 下降法 ( Index reductio n met hods) : 将高指数的方程 化为低指数的形式 。ADAMS12. 0 提供了三种系统 公式 : I3 就是前面分析所用的公式 , 它是系统默认 的公式 。SI2 ( Stabilized2Index 2) 是所谓的指数 2 公 式 ,公式的形式如式 ( 12b ) 所示 。SI1 ( Stabilized2In2 dex1) 是所谓的指数 1 公式 , 公式的形式如式 ( 12c ) 所示 。 SI2 公式 ( 12b ) 不但可以避免 J aco bi 矩阵的病 态 ,而且状态方程中出现了约束方程的导数项 ,所以 在求解中可以对系统状态的速度误差进行控制 , 求 解更为精确 。SI1 公式也有类似的结论 , 但它计算 量太大 ,一般不采用 。 2. 刚性微分方程的求解 所谓刚性是指系统各分量之间变化速度相差很 大 ,例如各部件由于质量差别大而导致的运动变化 在速度上的差别 。求解这样的系统时如果积分步长 总是用与变化快的分量所对应的小步长求解的话 , 会导致计算量过于庞大而使问题不可解 。所以 , 当 变化快的分量衰减后希望能够使用较大的步长来求 解系统方程 。机械系统一般都具有刚性 。
然后 ,对 ADAMS 动力学仿真参数的设置进行了总结 ; 最后 ,指出了本文的意义和进一步研究的方向 。
[ Abstract] Firstly , give an int roductio n mechanical dynamics equatio ns ; Seco ndly , give an analysis o n AD2 AMS stiff integrator algorit hm2 GIS TIFF ; Thirdly , give a summary o n ADAMS dynamic simulatio n parame2 ter s selectio n. Finally , point o ut t he meaning of t his article and t he way to f urt her st udy.
利用ADAMS与MATLAB对平面三级倒立摆的仿真
利用ADAMS与MATLAB对平面三级倒立摆的仿真孔凡国;相飞;钟廷志【摘要】利用ADAMS建立了平面三级倒立摆的虚拟样机模型,并采用ADAMS与MATLAB联合控制的方式,首次完成了平面三级倒立摆的虚拟样机控制仿真.仿真结果表明,尽管X向和Y向系统采用了相同的控制器且初始参数设置完全一样,但仿真结果二者存在一定的区别,这证明了X向和Y向系统存在耦合现象.此外,利用ADAMS与MATLAB联合进行系统的控制仿真,能即时反映控制算法的效果,具有很大的推广价值.%Built in ADAMS planar triple inverted pendulum virtual prototyping model, and uses the ADAMS and MATLAB joint control simulation ways, first completed virtual prototype control simulation of planar triple inverted pendulum. The simulation results show that although X and Y direction system to use the same controller and initial parameters settings exactly the same, but the simulation results to X and Y direction system are different, which proves to X and Y direction system exists coupling phenomenon. In addition, ADAMS and MATLAB joint system control simulation, can reflect the effect of real-time control algorithm, with a great promotion value.【期刊名称】《工业仪表与自动化装置》【年(卷),期】2012(000)003【总页数】3页(P3-4,30)【关键词】三级倒立摆;ADMS软件;虚拟样机;仿真【作者】孔凡国;相飞;钟廷志【作者单位】五邑大学,广东江门529020;五邑大学,广东江门529020;五邑大学,广东江门529020【正文语种】中文【中图分类】TP2730 引言倒立摆系统作为一个高阶次、不稳定、多变量、非线性、强耦合的复杂被控系统［1－2］，在对其控制过程中能有效地反映控制中的许多关键问题，如镇定问题、非线性问题、鲁棒性问题、随动问题以及跟踪问题等，是进行控制理论验证的理想实验平台。
基于ADAMS的机器人动力学分析及轨迹规划
基于ADAMS的机器⼈动⼒学分析及轨迹规划2.1 串联机器⼈在ADAMS中⽤连杆模拟机械臂，对两⾃由度的机械臂分别进⾏运动学分析、动⼒学分析及机械臂的轨迹规划。
2.1.1 运动学分析下⾯是建⽴模型并对模型进⾏设置分析的详细过程。
（1） 启动ADAMS/View，在欢迎对话框中选择新建模型，模型取名为Robot_arm，并将单位设置为MMKS，然后单击OK。
（2） 打开坐标系窗⼝。
按下F4键，或者单击菜单【View】→【Coordinate Window】后，打开坐标系窗⼝。
当⿏标在图形区移动时，在坐标窗⼝中显⽰了当前⿏标所在位置的坐标值。
（3） 创建机械臂关节1（连杆）。
单击连杆按钮，勾选连杆的长、宽、深选项，分别将其设置为300mm、40mm、10mm，如图2.1所⽰。
在图形区单击⿏标左键，然后将连杆拖⾄⽔平位置时，在单击⿏标左键。
（4） 在连杆的右端打孔。
在⼏何建模⼯具栏单击打孔按钮，将半径Radius设置为10mm，深度设置为10mm，如图2.2所⽰。
然后在图形区模型附近单击⿏标左键，在与XY平⾯垂直的表⾯上单击⿏标左键。
然后修改孔的位置，在孔附近单击⿏标右键，选择【HOLE_1】→【Modify】，在弹出的对话框中，将Center的坐标值设置成（300，0.0，5.0），如图2.3所⽰。
（5） ⽤（3）的⽅法在关节1右端孔中⼼处创建关节2，如图2.4所⽰。
然后再将关节2向内侧平移10mm。
2.1 创建连杆设置（6）添加约束。
在关节1的左端与⼤地之间添加转动副，在关节1与关节2结合处添加转动副。
单击⼯具栏中的旋转副按钮，并将创建旋转副的选项设置为2Bod-1Loc和Normal Grid，然后在图形区单击关节1和⼤地，之后需要选择⼀个作⽤点，将⿏标移动到关节1的Marker1处出现center信息时，按下⿏标左键后就可以创建旋转副，旋转副的轴垂直于⼯作栅格。
然后⽤同样的⽅法创建关节1与关节2之间的旋转副。
2013版ADAMS与Matlab联合仿真（绝对正确版）
声明：这一份matlab与adams联合仿真的例子是在原有例子的基础上改进的，因为原有的例子很坑爹的，把PID的控制线路图给标错了。
在这里我重新把图给改正回来了。
这样大家就可以看到一份正确的实例了我用2013版adams仿真算了，好人做到底，我把例子完全改为adams2013版的算了。
大家也容易看懂一点再讲一下PID的控制算法为什么要这么搭，这样大家就明白了PID（比例（proportion）、积分（integral）、微分（derivative）控制器）比例就是对误差乘以一个系数积分就是对误差积分然后再乘以一个系数微分是对误差求导这里我们的误差就是当前角度（angle）去减目标位置的角度（0度）=angle 即误差就是angle点击adams_sub后可以看到，上面才是角度angle，下面是速度（Velocity）=angle的微分（就是导数的意思）所以整个控制图是这个样子的。
其中表示的是对angle乘以一个系数，我乘了10 。
这就是所谓的比例调节，即P调节就是对angle积分，（就是累加的意思），即I调节，然后乘以了系数10因为下面输出的是角速度，角速度就是angle的微分的嘛，所以不用做什么操作然后乘以了系数10然后这三个相加起来，因为是反馈调节嘛，所以使用了，这个累减表示。
大家可以了解一下PID的公式，其实没有好难的。
现在大家来看看当时很坑爹的那个线路图吧，大家应该能看出来哪儿不一样了吧这个就是那个坑爹的图了，怎么会对速度积分嘛。
图7-35 变量Velocity随时间的变化然后居然还有角度与速度的图。
还都能控制的了。
分明是假的图嘛。
大家严重鄙视！现在奉上我修改过后的哈我是宾朋来自西南科技大学机器人小组谢谢！ADAMS/Controls使用实例本实例以MATLAB作为外部控制程序，以偏心连杆模型为例，讲解ADAMS与MATLAB的联合仿真过程。
主要包括创建机械系统模型、模型参数设置、建立MATLAB控制模型以及结果后处理四个步骤。
基于ADAMS的三缸柱塞泵动力学仿真分析_张子胜
零件名称 大齿轮
质量 /kg 1291 09
Ixx / ( kg# m2 )
51 380
Iyy / ( kg# m2 )
51 380
Izz / ( kg# m2 )
101670
小齿轮
71 76
01 022
01 022
01 036
曲轴
2501 10
21 456
1061 240 1061220
连杆
331 50
由美国 MDI公司开发的机械系统动力学仿真分 析软件 ( ADAM S) 是目前应用最广泛的动力学分析 软件, 能够有效地进行建模、仿真及后处理。由于 ADAM S仿真软件的建模功能目前不够强大, 因此笔 者使用三维建模软件 SOL IDEDGE 对三缸柱塞泵建 模, 建模完成后转换为 Pro /E 的文件格式, 再通过 P ro /E与 ADAM S 的接 口软件 M echpro 导入 到 ADAM S仿真软件中进行动力学分析[ 1, 2] 及后处理。
质的部分流体特性, 忽略阀堵的滞后动作, 认为液
缸内压强保持 35 M Pa无波动。建立以速度为条件
的分段函数来模拟阻力大小,
0
v柱塞 < 0
运动阻力 F = 0
v柱塞 = 0。
350
v柱塞 > 0
图 5 曲轴与连杆 1连接 处 X 向和 Y 向受力曲线
曲柄转角每 120b的受力周期相同, Y 向最大受
收稿日期: 2006- 10- 16 (本文编辑 赵连禄 )
( 上接第 27页 )
结
论
( 1) 在钻柱受力 分析的基础上, 根据疲劳强 度理论, 建立了钻柱在不对称循环应力状态下的疲 劳强度条件。
求解声波方程的辛RKN格式
求解声波方程的辛RKN格式刘少林;李小凡;刘有山;陈世仲【摘要】将声波方程变换至Hamiltion体系,构造了适用于高效声波模拟的二阶显式辛Runge-Kutta-Nystr(o)m(RKN)格式,运用根数理论得到此格式的阶条件方程组.针对两个自由度的辛条件方程组,根据三次项截断误差最小原理得到—种误差最小辛格式;通过分析声波的时间演进方程的稳定性,选择不同的辛系数使演进方程更稳定,并得到了另一种更为稳定辛格式;在频散关系分析中,选择使数值频散最小的辛系数,得到第三种最小频散辛格式.在理论分析中,这组辛RKN格式相比常见格式在精度控制、数值频散压制以及稳定性提升等方面均具有明显优势;在数值实验中,通过具体算例验证了理论分析的正确性.【期刊名称】《地球物理学报》【年(卷),期】2013(056)012【总页数】9页(P4197-4205)【关键词】声波波场;辛RKN格式;截断误差;稳定性;频散关系【作者】刘少林;李小凡;刘有山;陈世仲【作者单位】中国科学院地质与地球物理研究所,中国科学院地球深部研究重点实验室,北京 100029;中国科学院大学,北京100049;中国科学院地质与地球物理研究所,中国科学院地球深部研究重点实验室,北京 100029;中国科学院地质与地球物理研究所,中国科学院地球深部研究重点实验室,北京 100029;中国科学院大学,北京100049;中国科学院地质与地球物理研究所,中国科学院地球深部研究重点实验室,北京 100029;中国科学院大学,北京100049【正文语种】中文【中图分类】P6311 引言地震波正演模拟技术一直是地震学研究的热点，无论是中小尺度的地震波逆时偏移［1－2］、全波形反演［3－4］，还是区域性地震动研究［5－6］以及全球地震波模拟［7－9］和地球介质非均匀性反演等方面，正演模拟都起着关键作用.经过几十年的发展，地震波模拟方法已逐渐完善，总体而言，这些方法可以分为三大类，即射线法、积分方程法以及波动方程直接解法，Carcione等［10］和 Yang等［11］对这些方法做了详细的回顾与讨论.本文主要讨论声波方程的直接解法.无论是经典的有限差分法［12－14］，还是伪谱法［15－16］、有限元法［4，6］、谱元法［17－18］等众多方法，在时间离散上常常使用二阶中心差分或Newmark格式，虽然计算效率较高，但低阶的时间近似与高阶的空间离散格式往往不匹配，以致影响了地震波模拟的精度.针对时间精度较低等问题，Dablain ［19］将Lax等［20］提出的Lax－Wendroff方法运用至地震波模拟之中，其思想是运用高次空间微分修正时间精度.为了计算方便，一般得到时间四阶精度.但是直接运用有限长度的网格点波场值修正时间高次项，精度往往不够，针对此问题，Tong等［21］利用Yang等［22］发展的 Nearly Anayltic Discrete Method （NAD）近似空间高阶微分，得到了时间四阶、空间八阶的两步Stereo－Modeling Method（STEM）地震波模拟方法，两步STEM法在实际地震波模拟中精度的提高、数值频散压制和数值方向各向异性控制等方面均明显优于Lax－Wendroff方法.Chen在研究Lax－wendroff方法之后发现该方法在实际运用中误差增长较快［23］，针对此缺点，Chen将高阶的辛 Runge－Kutta－Nyström（RKN）与辛Partitioned－Runge－Kutta（PRK）格式运用至声波模拟之中［24－25］，因辛格式严格保持 Hamiltion系统微分二形式不变的特性［26］，有效地解决了地震波长时间计算中的误差累积等问题.Li等［27－28］在时间离散上采用三级三阶的辛PRK格式、空间离散上采用褶积微分算子法，形成了一套较为高效的地震波模拟方法，这些格式都是显式辛格式.与显式辛格式不同的隐式辛格式，其优点为无条件稳定性，Luo等［29－30］主要讨论了二阶隐式辛格式，但隐式辛格式涉及到微分算子的求逆运算，直接LU分解消耗大量的时间，谱因式分解虽计算速度快，但面临着精度损失、边界处误差过大等弊端，改进的混合算法虽然精度和计算效率上有所提高，但仍无法满足高精度与高效率地震波模拟的要求.Ma等［31］在分析总结前人工作基础上，提出在时间离散上采用二阶的辛PRK格式［32］、空间离散上采用NAD算子，得到了一种高效、高精度的地震波模拟方法.本文将伪谱法离散后的半离散声波方程变换至Hamilton系统，在分析不同时间积分格式面临的诸多问题之后，提出在时间离散上使用二级的辛RKN格式，以保证较高的计算效率，运用根数理论得到辛条件方程组.针对两个自由度的方程组，为了实际模拟需要，从精度、频散、稳定性三方面优化系数，最后得到了三组优化系数，理论上分析了新得到的优化格式在求解声波方程时的优良特性.在实际运用之中，可以根据不同的实际需要选择不同的格式以得到更好的模拟结果.2 辛RKN格式2.1 声波方程的伪谱法离散在地震波模拟与成像中有限差分算子得到了广泛的应用［19］，但有限长度的有限差分算子精度较低、数值频散明显［11］.伪谱法在精度和数值频散压制方面明显优于有限差分法［33］，随着计算机硬件的发展，伪谱法有望运用至三维高精度地震波模拟与成像之中［16］.考虑如下形式的地震声波方程：其中，u（x，y，z，t）为声波波场值，c（x，y，z）为声速.在空间上采用均匀网格对（1）式离散，即uijk≈u（iΔx，jΔy，kΔz），i＝1，2，…，Ni，j＝1，2，…，Nj，k＝1，2，…，Nk，Δx，Δy，Δz为空间网格间距.若采用伪谱法计算（1）式中的空间微分，得到的半离散方程如下：表示空间节点的离散波数值，记f（u）＝c2 F－1［w2·F（u）］.其中，u为离散后的波场向量，F，F－1分别表示Fourier正变换和逆变换，w ＝［w1，1，1，…，wNxNyNz］，2.2 辛RKN格式及条件方程对于（2）式可以采用多种时间离散格式，如三阶的 Runge－Kutta方法（RK）［34］、Newmark 格式［5，7，8］、二阶辛 PRK 格式［31－32］、三阶辛PRK 格式［23，27－28，34］、四阶的辛RKN格式［24－26］以及最近Yang 等发展的显式化后的 RK格式与 Adams格式［35－36］.这些方法中二阶的辛PRK格式效率最高，每一个时间步只需两次正反Fourier变换，高阶格式虽然有较高的精度，但每个时间步计算量过大以致严重影响计算效率.在大尺度以及全球尺度地震波模拟时，由于空间微分计算量过大，每个时间步无法承受过多中间步，所以本文使用二阶辛RKN格式对（2）式进行离散.引入变元v＝u·，则方程（2）可降阶为考虑（3）式为 Hamilton系统情形，即（3）式可表示为如下形式：其中 H（u，v）为一标量函数，由（3）、（4）式知 H 满足如下等式：因此H满足：当f（u）为某一标量函数V的梯度时才可考虑（4）式的辛算法.对于一般的s级RKN格式可以写为［26］：其中，h为时间步长，ci，aij，¯bj，bj为辛系数.可以证明，如果（6）式是显式辛的，（6）式应满足如下方程［26，37］：在地震波模拟时，希望格式（6）的计算效率与二阶的PRK格式相近［31－32］，当s＝2时即满足高效地震波模拟的要求.在精度上希望格式（6）精度尽可能提高，对于非约化的辛RKN格式［26，37］（即格式中无冗余级，bj≠0），满足（7）式的格式（6）是p 阶的，可以用图形表示微分关系的根数理论［37－38］，即任何s－树sτ，都有一个有根的s－树ρsτ由sτ的一胖节点提升得到，ρsτ的权与密度满足如下等式：其中φ为ρsτ的权，γ为ρsτ的密度.如果二级的辛RKN格式是三阶的，有4个非多余有根s－树需满足（8）式，将辛条件（7）带入（8）式化简得到如下等式：2.3 一种误差最小辛RKN格式求解（9）式发现，方程无实数解，即二级的辛RKN格式无法使精度达到三阶，我们可以根据截断误差最小原理得到精度接近三阶的误差最小辛RKN格式，即（9）式中的前两式满足的同时，构造目标函数，使目标函数尽量最小，运用非线性优化方法，可以得到如下一种误差最小辛RKN格式（记为M1）：2.4 一种强稳定的辛RKN格式运用二级辛RKN格式进行地震波模拟时，时间步长与空间步长应满足一定的比例关系，否则计算失稳以至无法进行.为了使分析过程中量纲一致，引入变量~v＝vh，考虑如下的简谐解，由（6）式得到的时间演进方程为其中，θ＝ckh，e1 ＝b1c21＋b1c22，e2 ＝b1b2（c2－c1）.由（13）式知时间演进方程的特征值为其中为了保证＝1，由（14）式可以得到如下不等式：选择不同的c1，c2，满足（9）、（15）式的同时力图使θ取到最大值.最后得到的一种强稳定的辛RKN格式（记为 M2）：2.5 一种最小频散辛RKN格式由于用时间和空间的离散点值近似替代时间空间连续变化的函数，必然导致真实信号中的部分信号无法拾取而导致数值频散.本文用伪谱法计算空间微分，在Nyquist波数范围内伪谱法对波数完全覆盖，用格式（6）进行地震波模拟时，数值频散将主要来源于时间离散.由（13）、（14）式可知，（13）式左端的exp（iωh）应与λ相等，由此可以得到相速度c与真实速度c之比满足如下等式：其中，θ又可以表示为θ＝rkΔ，库朗数选择适当的库朗数（如r＝1），当kΔ由0变化至π，选择不同的系数使得相速度与真实速度最接近，为此构造如下目标函数：用非线性最优化方法得到一种最小频散辛RKN格式（记为 M3），3 三种辛RKN格式特性分析第2节从精度、稳定性和数值频散优化三方面得到了3种辛RKN格式，本节就这三方面与常见格式对比，以凸显本文格式的特性.对比选用的格式包括：二阶辛PRK 格式［31－32］（记为 M4）、二阶优化PRK格式［39］（记为 M5）和三级三阶PRK格式［27－28］（记为 M6）.3.1 精度分析将（10）式中的E1视为截断误差，对比M1—M6 6种格式的E1.对比结果如表1所示.就精度而言，三阶格式M6精度最高；M1精度其次，其精度逼近三阶格式M6；M2与M5精度接近，M3和M4精度最差.表1 6种辛格式截断误差对比Table 1 Error truncation of six symplectic schemes方法误差M1 M2 M3 M4 M5 M6截断误差0.000292 0.00217 0.0179 0.0087 0.00204 03.2 稳定性分析按照2.4节分析方法，可以得到M1—M6的稳定性条件.表2给出了6种格式从一维到三维的稳定性极限（即库朗数r所取的极大值）.从表2可以看出，理论上M3稳定性最好.在常见格式中三阶格式M6稳定性最好，二阶的M4格式稳定性最差，优化的二阶格式M5介于两者之间，整体而言常见格式稳定性要逊色于本文3种格式.表2 6种格式一维到三维稳定性比较Table 2 Stability limits of six symplectic schemes for 1D，2D，and 3Dscalar wave simulations维数方法1D 2D 3D M1 0.8127 0.5746 0.4692 M2 1.2732 0.9003 0.7351 M3 1.2436 0.87930.7180 M4 0.6366 0.4502 0.3676 M5 0.7208 0.5097 0.4162 M6 0.84860.6000 0.48993.3 频散分析按照2.3节的分析方法，可以得到M1—M6的频散关系，图1（a、b）为r＝0.5和r＝0.8时6种格式1D情形下的频散曲线.由图1a可知，r取较小值时M6频散最小，M3频散其次，M4与M5频散较大；在图1b中当r＝0.8时M4与M5已经失稳，故未在图中显示，当kΔx取较小值时，M6频散最小，但kΔx取较大值时，M6频散将超过M3.由图对比可知，M3频散一直维持在较低范围.当r＝0.8时，按照Basabe和Sen［40］频散限定标准，即频散小于0.01，1D情况下，M3一个波长的最小采样点数为2.38；同理2D情形下，M3一个波场内的最小采样点数为3.36.图1 6种格式r＝0.5和r＝0.8时的频散关系对比Fig.1 Comparison of numerical dispersion of six symplectic schemes when r＝0.5（a）and r＝0.8（b）4 数值试验4.1 半无限均匀介质模型为了测试这组辛RKN格式对声波的模拟精度，设计了2D均匀介质模型，模型大小为2550m×2550m，震源位于模型中心，接收点在震源左侧400m处.震源由主频为30Hz的Ricker子波激发，空间网格间隔为10m，时间步长为1ms.声波波速为3000m／s.定义接收点的总体误差其中，ui为t＝i h时刻的数值解，ur为解析解.模拟结果如图2，计算效率与总体误差如表3所示.由图2可知，6种方法模拟结果与解析解高度一致，M1、M2和 M6与解析解最靠近，M3、M4和M5偏离解析解较远，由于压制高波数处的数值频散，M3对低波数模拟的精度影响较大，表现在图2中为相位较为超前.由表3可以看出，M6误差最小，M1与M2误差其次，其中M1与M2精度足够接近；M3—M5误差较大；就计算效率而言M6的计算时间为M1—M5的1.5倍，在计算效率相同的情况下，M1、M2比M4、M5有较大的精度提升.4.2 分层介质模型为了检验本文方法的数值频散压制能力，设计如图3的分层介质模型，模型大小为3825m×3825m，分界面位于1920m深度处，上层介质的波速为2000m／s，下层介质波速为3000m／s，震源位于（1920m，－1710m）处，由主频为60Hz的Ricker子波激发产生，时间间隔为1ms，空间间隔为15m.模拟结果如图4所示.图2 6种格式计算得到的声波波形与解析解对比Fig.2 Comparison of waveforms generated by six sympletic methods and analytic solution图3 分层均匀介质模型以及模型参数Fig.3 Two layer homogeneous mediaand model parameters图4（a—d）分别为M4—M6和M3四种方法计算得到的0.5s时刻的波场快照.从图4a可以看到除了直达波、反射波、透射波以及首波外，还在直达波波前、透射波波前附近出现了强烈的数值频散，图图4 4种方法模拟分层声波介质中地震波传播0.5s时的波场快照（a—c）M4—M6方法；（d）M3方法.Fig.4 The snapshots of scalar wave propagation in two layered medium at time t＝0.5s表3 6种辛格式计算效率和总体误差对比Table 3 Computational efficiencyand total error comparison of six symplectic methods方法参数M1 M2 M3M4 M5 M6内存消耗／kB 3656 3656 3656 3656 3656 3656计算时间／s 827.450 828.280 829.500 826.390 830.700 1253.330总体误差E 2.070 2.071 3.262 2.984 2.292 1.5914b在透射波波前也出现了强烈的数值频散；图4（c，d）中只在透射波前出现了轻微的频散，M3和M6计算结果较为一致.M3—M5计算效率基本相同，M6计算耗时为M3—M5的1.5倍，在高频地震波模拟时使用M3，在保证计算效率的同时，能较好地压制数值频散.4.3 非均匀介质模型——SEG／EAGE盐丘模型为了测试本文方法在非均匀介质中波场计算的有效性和稳定性，选择SEG／EAGE 模型做测试，模型速度结构如图5所示.模型中波速变化范围为1524～4480m／s，除了若干个起伏分层界面外，还存在盐丘高阻体，使得介质横向速度变化极为强烈.选择震源位于模型中心，其为主频是40Hz的Ricker子波，空间网格间距为20m，时间间隔为1ms时M2、M5和M6三种方法的模拟结果如图6所示，时间间隔为2ms时的模拟结果如图7所示.当时间间隔为1ms时，M2、M5和M6可以得到几乎相同的波场快照，从图6中可以看到，由于介质的非均匀性，在速度突变的界面上产生了强烈的反射和多次反射，以及在速度突变的角点处产生了明显的绕射和散射，三种方法均是稳定的，并得到了可靠的结果.当时间间隔增加一倍，即为2ms时，M5已经失稳，而M2和M6依然稳定，从图7中可以看到两种方法得到的快照和图6中的快照无论是整体还是细节上都极其一致，说明本文M2使用大时间步长的模拟结果依然是可靠的.图5 SEG／EAGE盐丘模型Fig.5 The velocity profile of SEG／EAGE salt model图6 三种方法模拟得到的非均匀介质中0.5s时的波场快照（a）M2；（b）M5；（c）M6，时间间隔为1ms.Fig.6 The snapshots of wavefields in heterogeneous media when t＝0.5s（a—c）are calculated by M2，M5，and M6，respectively.Time interval is 1ms.图7 两种方法模拟得到的非均匀介质中0.5s时的波场快照（a）M2；（b）M6，时间间隔为2ms.Fig.7 The snapshots of wavefields in heterogeneous media when t＝0.5s（a—b）are calculated by M2，and M6，respectively.Time interval is 2ms.5 讨论本文对地震声波方程空间上使用伪谱法离散，时间离散上选用高效的二阶辛RKN格式，通过以误差最小、稳定域最大和频散最小为依据构造了三种辛RKN格式.在理论分析中，和常见格式对比，理论上论证了本文方法在精度、稳定性和数值频散等方面的优势；在数值实验中，通过与解析解、分层介质和非均匀介质中不同方法波场模拟结果对比，数值结果进一步佐证了本文方法在保证计算效率的同时，在精度提高、稳定域增加和数值频散压制等方面均具有明显改进.可以根据不同实际需要选择不同方法，如在高精度地震波模拟时选择M1，在高频地震波模拟时选择M3，在强烈非均匀介质中地震波模拟时选择M2.本文方法为高效地震波模拟、成像提供了一种可靠的选择.参考文献（References）［1］ Whitmore N D.Iterative depth migration by backward time propagation.53rd Annual International meeting，SEG，Expanded Abstracts，1983：827－830.［2］刘红伟，李博，刘洪等.地震叠前逆时偏移高阶有限差分算法及GPU实现.地球物理学报，2010，53（7）：1725－1733.Liu H W，Li B，Liu H，et al.The algorithm of high order finite difference pre－stack reverse time migration and GPU implementation.Chinese J.Geophys.（in Chinese），2010，53（7）：1725－1733.［3］ Song Z M，Williamson P R，Pratt R G.Frequency－domain acoustic －wave modelling and inversion of cross－hole data，PartⅡ：Inversion method，synthetic experiments and realdata results.Geophysics，1995，60（3）：796－809.［4］张美根，王妙月，李小凡等.时间域全波场各向异性弹性参数反演.地球物理学报，2003，46（1）：94－100.Zhang M G，Wang M Y，Li X F，et al.Full wavefield inversion of anisotropic elastic parameters in the time domain.Chinese J.Geophys.（in Chinese），2004，46（1）：94－100. ［5］ Komatitsch D，Liu Q，Tromp J.Simulations of ground motion in the Los Angeles basin based upon the spectralelementmethod.Bull.Seismol.Am.Soc.，2004，94（1）：187－206.［6］张怀，周元泽，吴忠良等.福州盆地强地面运动特征的有限元数值模拟.地球物理学报，2009，52（5）：1270－1279.Zhang H，Zhou Y Z，Wu Z L，et al.Finite element analysis of seismic wave propagation characteristics in Fuzhou basin.Chinese J.Geophys.（in Chinese），2009，52（5）：1270－1279.［7］ Komatitsch D，Tromp J.Spectral－element simulations of global seismic wave propagation—I.Validation.Geophys.J.Int.，2002，149（2）：390－412.［8］ Komatitsch D，Tromp J.Spectral－element simulations of globalseismic wave propagation—Ⅱ.Three－dimensional models，oceans，rotation and self－gravitation.Geophys.J.Int.，2002，150（1）：303－318. ［9］ Yan Z Z，Zhang H，Yang C C，et al.Spectral element analysis on the characteristics of seismic wave propagation triggered by WenchuanMs8.0earthquake.Science in China Series D：Earth Sciences，2009，52（6）：764－773.［10］ Carcione J M，Herman G C，ten Kroode A P E.Seismic modeling.Geophysics，2002，76（4）：1304－1325.［11］ Yang D H，Tong P，Deng X Y.A central difference method with low numerical dispersion for solving the scalar wave equation.Geophysical Prospecting，2012，60（5）：885－905.［12］ Virieux J.SH－wave propagation in heterogeneous media：Velocity －stress finite－difference method.Geophysics，1984，49（11）：1933－1942.［13］ Virieux J.P－SV wave propagation in heterogeneous media：Velocity－stress finite－difference method.Geophysics，1986，51（4）：889－901.［14］ Yang D H，Liu E，Zhang Z J，et al.Finite－difference modelling in two－dimensional anisotropic media using a fluxcorrected transport technique.Geophys.J.Int.，2002，148（2）：320－328.［15］ Kolsoff D D，Baysal E.Forward modeling by a Fouriermethod.Geophysics，47（10）：1402－1412.［16］龙桂华，李小凡，江东辉.基于交错网格Fourier伪谱微分矩阵算子的地震波场模拟GPU加速方案.地球物理学报，2010，53（12）：2964－2971.Long GH，Li X F，Jiang D H.Accelerating seismic modeling with staggered－grid Fourier Pseudo－spectral differentiation matrix operator method on graphics processing unit.Chinese J.Geophys.（in Chinese），2010，53（12）：2964－2971.［17］ Komatitsch D，Barnes C，Tromp J.Wave propagation near a fluid－solid interface：A spectral－element approach.Geophysics，2000，65（2）：623－631.［18］ Komatitsch D，Barnes C，Tromp J.Simulation of anisotropic wave propagation based upon a spectral element method.Geophysics，2000，65（4）：1251－1260.［19］ Dablain M A.The application of high－order differencing to the scalar wave equation.Geophysics，1986，51（1）：54－66.［20］ Lax P D，Wendroff B.Difference schemes for hyperbolic equations with high order of munications on Pure and Applied Mathematics，1964，17（3）：381－398.［21］ Tong P，Yang D H，Wang M X.A high－order stereomodeling method for solving wave equations.Bull.Seism.Am.Soc.，2013，103（2A）：811－833.［22］ Yang D H，Teng J W，Zhang J F，et al.A nearly analytic discrete method for acoustic and elastic wave equations in anisotropicmedia.Bull.Seismol.Soc.Am.，2003，93（2）：882－890.［23］ Chen J B. High－order time discretizations in seismicmodeling.Geophysics，2007，72（5）：115－122.［24］Chen J B.Modeling the scalar wave equation with Nyströmmethods.Geophysics，2006，71（5）：151－158.［25］ Chen J x－Wendroff and Nyström methods for seismic modelling.Geophysical Prospecting，2009，57（6）：931－941.［26］ Okunbor D，Skeel R D.Explicit canonical methods for Hamiltonian p.，1992，59（200）：439－455.［27］ Li X F，Li Y Q，Zhang M G，et al.Scalar seismic－wave equation modeling by a multisymplectic discrete singular convolution differentiator method.Bull.Seismol.Am.Soc.，2011，101（4）：1710－1718.［28］ Li X F，Wang W S，Lu M W，et al.Structure－preserving modelling of elastic waves：a symplectic discrete singular convolution differentiator method.Geophys.J.Int.，188（3）：1382－1392.［29］ Luo M Q，Liu H，Li Y M.LU decomposition with spectral factorization in seismic imaging.Chinese J.Geophys.，2003，46（3）：602－612.［30］ Luo M Q，Zhu G T，Liu H，et al.A hybrid matrix inversion method for 3－D implicit prestack depth migration.Chinese J.Geophys.，2003，46（5）：978－987.［31］ Ma X，Yang D H，Liu F Q.A nearly analytic symplectically partitioned Runge－Kutta method for 2－D seismic waveequations.Geophys.J.Int.，2011，187（1）：480－496.［32］孙耿.波动方程的一类显式辛格式.计算数学，1997，（1）：1－10.Sun G.A class of explicitly symplectic schemes for waveput.Math.（in Chinese），1997，（1）：1－10.［33］ Fornberg B.High－order finite differences and the pseudospectralmethod on staggered grids.SIAM J.Numer.Anal.，1990，27（4）：904－918.［34］汪文帅，李小凡，鲁明文等.基于多辛结构谱元法的保结构地震波场模拟.地球物理学报，2012，55（10）：3427－3439.Wang W S，Li X F，Lu M W，et al.Structure－preserving modeling for seismic wavefields based upon a multisymplectic spectral element method.Chinese J.Geophys.（in Chinese），2012，55（10）：3427－3439.［35］ Yang D H，Wang N，Chen S，et al.An explicit method based onthe implicit Runge－Kutta algorithm for solving waveequations.Bull.Seismol.Soc.Am.，2009，99（6）：3340－3354.［36］ Yang D H，Wang Lei，Deng X Y.An explicit split－step algorithm of the implicit Adams method for solving 2－D acoustic and elastic wave equations.Geophys.J.Int.，2010，180（1）：291－310.［37］冯康，秦孟兆.哈密尔顿系统的辛几何算法.杭州：浙江科学出版社，2003.Feng K，Qin M Z.Symplectic Geometric Algorithms for Hamiltionian Systems （in Chinese）.Hangzhou：Zhejiang Science ＆ Technology Press，2003.［38］ Hairer E.Geometric Numerical Intergration I （2nd ed）.Berlin and New York：Springer－Verlag.［39］ MaLachlan R I， Atela P. The accuracy of symplecticintegrators.Nonlinearity，1992，5（2）：541－562.［40］ Basabe J D D，Sen K M.Grid dispersion and stability criteria of some common finite－element methods for acoustic and elastic waveequations.Geophysics，2007，72（6）：T81－T95.。
一种使用隐式Adams方法求解三维
波场模拟
图(a)空间步长是40米，时间步长是0.004秒，频率是15Hz 图(b)空间步长是30米，时间步长是0.003秒，频率是15Hz 图(a)上下层速度分别为：3.5km/sec, 4.5km/sec 图(b)上下层速度分别为：2km/sec, 4km/sec
波场模拟
两图上下层速度都分别为：3km/sec, 5km/sec 图(a)空间步长是40米，时间步长是0.004秒，频率是15Hz 图(b)空间步长是30米，时间步长是0.003秒，频率是15Hz
2 3
4
V in, j ,k
− pL + rpL + r 2 L + r 3 pL
2 3
4
V in,−j1k ,
0 I 3×3 0 I 3×3 0 I 3×3 0 I 3×3 , , , L = Diag 1 1 1 1 D 0 D 0 D 0 D 0 ρ ρ ρ ρ
与精度是O(∆t 4 + ∆x 4 + ∆y 4 + ∆z 4 )的LWC
和精度是O(∆t 2 + ∆x 2 + ∆y 2 + ∆z 2 )的FDM比较
误差分析（数值误差）
误差分析（数值误差）
计算效率
∂ 2u ∂ 2u ∂ 2 u ∂ 2u = c2 2 + 2 + 2 + f ∂ 2t ∂ x ∂ y ∂ z
3D SSA方法
r= 5∆t 8∆t ∆t ,s = ,p= 12 12 12
(I
V
24×24
−r⋅L V
)
n +1 i, j ,k
基于ADAMS的超声波探头架仿真与研究
分 析 ， 绘 出弹 簧 的受 力 情 况 以及 探 头 的运 动 轨迹 ， 描 验证 了 机 构 设 计 的 合 理 性 ， 析 结 果 为 设 计 提 供 了科 学 依 分
据 ，对 设 计 工 作 具 有 很 重 要 的指 导 意 义 。
关 键 词 ：Ａ ＤＡＭＳ 动 力 学 ；超 声 波 ；探 头 架 ； 中 图分 类 号 ：Ｔ ７ ２ Ｈ８ ８ ． 文献 标 识 码 ：Ａ
收 稿 日期 ：２ ０ —６０ ；修 回 日期 ：２ ０— ９１ ０ ６０ —２ ０ ６０ — ７ 作 者 简 介 ：胡 金 萍 （ ８ 一 ， ， 苏 东 台人 ， 士研 究 生 。 １ １） 男 江 ９ 硕
图 １
探 头 架 立 体 图
２ 添 加 约 束
在 进 行 模 型 仿 真 分 析 前 我 们 要 对 探 头 架 添 加 约
基 于 ＡＤ ＡＭＳ的超 声 波探 头架仿 真 与研 究
胡金 萍 ，任 建 平
（ 北 大 学 机械 工 程 与 自动 化 学 院 ， 山西 太 原 中 ００５） ３ ０ １
摘 要 ：运 用 Ｓ ｌ Ｗｏ ｋ ｏｉ ｄ ｒ ｓ软 件 建 立 了探 头 架 的立 体 模 型 ，将 其 成 功 导 入 ＡＤ ＡＭＳ，对 此 机 构 进 行 了动 力 学 仿 真
０ 引 言
不变 。
手 推式超 声波 钢 轨探伤 仪是 一种 能够对 铁 路钢 轨 进 行连续 探 伤 的仪 器 ，具有 多路 发射 接收 系统 和报警 闸 门，一次 推行 中可 检 出钢 轨上 的各 种有 害缺 陷 ，其 操 作简 单 ，重量轻 、功耗低 。当探测 到钢 轨 内部 的损 伤时 ，具有 自动报 警 功能 ，探 伤人员 就可 以及 时对 损 伤进行 处理 。超声 波 传感 器在 手推式 钢 轨探伤 仪 中的 作用非 常重 要 ，超声 波直接 与 钢轨接 触 ，并发 射超 声 波进行 钢 轨 内部 探伤 ，同时接 受钢 轨返 回波 。如果 探 头架设 置 不合理 ， 会导致 超声 波传 感 器的意 外损 坏 ， 将 严重 时影 响它 的使 用寿命 ，同时探 头架 也是受 磨损 和 冲击 比较 多 的器件 之一 ，因 此有必 要对 其进 行运 动学 的分 析 ，验证其 设计 的合 理性 ，使 手推 式超声 波 探伤 仪能 够长 久耐 用 。
ADAMS方法和一般线性多步法
(5.3)
ji
x xk j = x xk + xk xk j = th + jh xk i xk j = xk i xk + xk xk j = ih+ jh
并用 yk+1, yk 分别近似 y(xk+1), y(xk ), 得 l 步显式Adams公式：
l1
yk+1 = yk + h
f l ,i k i
01
22
l i
= =
02÷÷
2,1 =
1 0
t dt = 1
1, 2
l i
= =
12÷÷
由 (5.6)得
C p+1 = C3 =
1 2
1
t(t + 1)dt =
5
.
0
12
k +1
=
5 12
y(3)(
k )h3.
由(5.2) 得
yk +1
=
yk
+
h 2 (3 fk
fk 1 ).
(5.6)
Cl+1 =
1 l!
yk+1
=
yk
+
h (23 12
fk
16fk 1 + 5 fk 2)
yk +1
=
yk
+
h (55 24
fk
59fk 1 + 37fk 2
9 fk 3)
C p+2
1 12 5 24
19 720
3 160
Cp+1
1 2
5 12 3 8 251 720
Created with SmartPrinter trail version
Adams公式的新证明
解得 :A = 入公式 :
251 646 264 106 19 ,B = ,C = ,D = ,E = 。代 720 720 720 720 720
其中 h 为步长 , x i + 1 = x i + h , y i = y ( x i ) , f i = f ( x i , y i ) 。 由于式 (8) 中有 4 个未知系数 , 只需 4 个方程 。由定 理 1 ,取 f ( x , y ) 分别为 1 , x - x i , ( x - x i ) ( x - x i + h) , ( x
x i ) dx = h
2
h (0 A - hB - 2 hC - 3 hD) = h ( 0 A + 0 B + 2 h C + 6 h D)
2 2
∫(x x i
x i +1
2
参考文献
5h 6
3
=
∫
x i
x
[ 1 ] 阑小林 ,蒋耀林 . 现代数值分析 [ M] . 北京 : 国防工业出版社 , 2004. 9 , 321 - 323. [ 2 ] 于寅 . 高等工程数学 [ M] . 武汉 : 华中科技大学出版社 ,2001. 10 ,482 - 486.
i +1
( x - x i ) ( x - x i + h) dx =
3
h ( 0 A + 0 B + 0 C - 6 h D)
=
∫
x i
x
i +1
( x - x i ) ( x - x i + h) ( x - x i + 2 h) dx =
9h 4
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
创新项目论文一种基于三阶Adams 格式的求解声波方程的多步算法China University of Mining & Technology-Beijing摘要一个准确、高效、低数值频散的正演算法能够提高反演精度、加快反演收敛速度，因此研究地震波场正演模拟技术具有重要意义。
区别于传统的空间离散方法，利用空间插值, 用网格点处的函数值及其梯度共同逼近空间高阶偏导数的方法称为近似解析离散化方法。
声波方程通过变换，并采用近似解析离散化方法进行空间离散，从而转变成为一个半离散化的常微分方程组，再利用三阶显式Adams格式进行时间推进，求解半离散化的常微分方程组，从而得到了一个新的求解声波方程的有限差分方法（AD-STEM）。
对AD-STEM进行了理论误差和数值误差分析、计算效率比较和数值波场模拟。
研究表明，与传统方法AD-LWC比较，AD-STEM方法数值精度更高，数值频散更低，更高效，且与解析解匹配更好。
AD-STEM方法能够通过压制数值频散而提高计算效率。
在无可见数值频散的条件下，AD-STEM的计算速度是AD-LWC的1.88倍，而存储量只有其72%，更适合在粗网格下进行大规模地震波场数值模拟。
关键词：近似解析离散化方法；三阶Adams格式；数值频散；有限差分目录1 绪论 (1)1.1选题背景和研究意义1.2粘弹性介质国内外研究现状1.3有限差分国内外研究现状1.4本文主要研究内容2 粘弹性介质的基本模型 (6)3方法介绍....................................................................................................................... 错误！未定义书签。
3.1 Stereo-modeling方法简介 (10)3.2 Lax-Wendroff correction方法简介 ...................................................................... 错误！未定义书签。
4 粘弹性介质中的波场数值模拟..................................................................................... 错误！未定义书签。
4.1 波场快照 (11)4.2 波形图.................................................................................................................. 错误！未定义书签。
4.3 SEG模型的地表地震记录 (14)5 结论 (18)6 参考文献 (20)1 绪论1.1 选题背景和研究意义随着我国石油和天然气工业的迅速发展，油气勘探的精度和难度也在不断地加大。
这就要求我们必须深入了解地震波在地下介质中的传播规律，以便我们能更加准确地进行地下油气的勘探。
地震波数值模拟技术是目前地震勘探领域人们研究地震传播规律的一个重要手段。
传统地震波数值模拟技术，一般假设地层介质是均匀和完全弹性的，但是我们发现：在这种假设情况下，数值模拟的结果和我们在野外观察到的实际记录有很大差别。
造成这种差别的主要原因是因为实际的地层介质并不是均匀和完全弹性的，我们应该用粘弹性介质来代替，这对于我们研究地震波波场的传播规律具有了非常重要的理论和实际意义。
1.2 粘弹性介质国内外研究现状为了对实际地震资料的解释以及地震波的偏移反演等处理提供可靠依据，我们需要建立不同的介质模型来模拟地下介质，发展相应的理论和勘探方法。
从以往的模拟结果来看，以经典弹性波理论为基础得到的理论记录实际地震波在大地中传播时所接收到的记录之间存在很大差别，为了消除这种差别，我们通常要对球面扩散、界面的透射损失等进行补偿，但是在补偿过后，理论记录上来自深层的地震波振幅仍然较小、低频成分较丰富，这说明地下介质对地震波有吸收衰减的作用，特别是对高频成分，吸收更为严重。
地震波衰减受多种因素影响，事实上现在没有一种机理可以描述所有环境条件下产生的损耗，介质的粘滞性就是其中的一个主要原因，为了弄清楚地震波在地层中的衰减机制，很多研究人员开始对粘弹性介质的理论和应用方法进行研究。
（一）国外研究现状早在1845年Stokes就开始研究粘弹性介质中的地震波，他提出要考虑具有耗损的固体就必须使这种固体具有类似于粘滞液体的性质，因此他给固体的剪切模量附加上一个剪切粘滞，并建立了由于粘滞型内摩擦引起的能量耗损的Stokes 粘弹性波动方程。
但是粘弹介质理论的发展却是非常缓慢的，直到上个世纪40年代，美国地球物理学家N.K.Ricker首次在粘弹介质理论中完整详细的描述了地震波在地下介质传播中时的衰减问题。
其后，人们纷纷开始研究粘弹介质中地震波的传播理论，提出了一系列的理论和研究方法。
Aki和Richards(1980)建立了一种新的粘弹性介质模型——标准线性固体模型，相对于开尔芬—佛格特模型和麦克斯韦尔模型，标准线性固体模型更加适合描述地震波衰减的物理机制。
粘弹介质地震波数值模拟技术的研究是从上世纪八十年代开始的，Day和Minster (1984)第一次成功的在粘弹性介质中利用Pade近似法进行2-D时间域数值模拟。
EmmeriCh和Korn (1987)发现这种方法存在明显的缺陷即质量低劣和计算效率低，于是他们提出了一种新的近似称之为“广义标准线性固体”(GSLS)，同时对粘弹性模量的有理式进行推导，并发展二维有限差分算法使之适合标量波的传播，提高了计算效率。
CarCione(1988)等对粘滞声波在地层中传播的进行了模拟并推导出了模拟方程。
Robertsson(1994)等设计了一种速度一应力有限交错网法，并利用它研究了二维粘弹性介质下纵波和横波的传播规律，随后他们又将该算法推广到三维介质情况。
(二)国内研究现状跟国外关于粘弹性介质理论及其数值模拟的迅速发展，国内对这方面研究起步较晚。
近年来国内地球物理学家也开始对地震波在粘弹性介质中传播进行研究。
张剑锋、李幼铭(1994)把地层介质假设为水平层状介质，利用混合变量粘弹性波方程，直接逐层求解位移、应力的方法进行数值模拟四。
毕玉英、杨宝俊(1995)给出一种实现二维粘滞介质完全波场模型计算的方法。
该方法的独到之处在于将传播时间分解成了传播水平时间和传播垂直时间两部分，平面源人射改为线源人射，无需繁杂的射线追踪，只考虑入射角随偏移距的变化情况便可获得包括反射纵波、转换波、多次波、折射波及面波等在内的多种波场的时间一空间道集的正演模拟记录，而且还能灵活地模拟井间和VSP地震剖面。
宋守根(1996)提出了一种提高地震剖面纵横向分辨率的粘弹性介质波场外推方法，利用该方法进行反演时，无需对积分方程作线性化近似，适应性较强，并通过实例对该理论推断予以验正。
张建贵(1999)等针对塔里木盆地的沙漠覆盖面大，地层埋藏深，地质构造幅度小的特点以大地介质为粘弹性介质为前提，利用粘弹性介质波动方程聚焦成像技术得到了一套高分辨率、高信噪比和高保真的地震处理剖面。
孟凡顺(2000)等人根据粘滞弹性理论，推导出了粘滞弹性波的有限差分运动方程并对任意复杂地质体粘滞弹性波进行了正演模拟；程昌钧(2001)等系统研究了粘弹性介质中波的逆散射问题及其求解过程；崔建军和何继善 (2001)以二维粘弹性波动方程为基础，研究了粘弹性波动方程F-K域的正演模拟和偏移方法的思路和理论基础，并阐明了推导过程；杨午阳(2003)提出了一种利用粘弹性声波波动方程进行偏移的新方法；奚先和姚姚(2004)通过波动方程的交错网格有限差分数值模拟方法，对地震波在二维粘弹性随机介质中的传播进行了模拟并研究了其波场特征；朱慧卿(2004)利用新定义的各向异性网格一广义拟一致网格对粘弹性方程各向异性有限元方法的超收敛性进行了分析[22]；刘财、张智(2005)以积分本构方程为基础，应用对应原理建立波动方程，对线性粘弹体中地震波场进行伪谱法模拟；宋常瑜(2006)等采用二维粘弹性波方程交错网格高阶有限差分数值模拟方法进行了井间地震粘弹性波场数值模拟，并对地震波的传播规律和衰减特征进行了研究；邵志刚(2007)将以往两种粘弹介质中地震波模拟方法的优点结合起来，以模型理论和积分本构方程为基础，从理论上分析了模型对地震波场的影响并采用交错网格有限差分法对粘弹介质中的地震波进行数值模拟；孙成禹(2007)针对以往粘弹性模型在描述介质品质因子对频率的依赖关系方面存在不足，与实际观测结果不符的问题，提出了一种新的粘弹介质模型的数值方法，该方法较好地描述了品质因子对频率的非依赖性，可以从理论上较准确地对实际介质的粘弹性对地震波场的影响进行研究；唐启军(2009)将Von Karman型的随机各向同性背景引入粘弹性单斜各向异性波动方程，并应用交错网格技术进行模拟；Yong Yun-dong (2010)等人推导出了交错网格有限差分的并行算法，并将其用于三维粘弹性随机溶洞介质数值模拟中，取得了良好的效果。
1.3有限差分国内外研究现状地震数值模拟作为地震勘探基础性研究，是我们认识地震波传播规律，检验各种处理方法正确性的重要工具，地震波的数值模拟是是研究地震波传播规律的最有效工具和手段，在地震勘探和地震学各工作阶段都有重要的作用，贯穿于地震资料的采集、处理、解释的整个过程中，同时也是地震反演的基础。
有限差分法是最常用的一种数值模拟方法，基本原理就是在解偏微分方程时用有限差分算子代替微分，将微分方程化为相关的线性代数方程，通过求解代数方程，得到偏微分方程的数值解。
有限差分法数值模拟技术相对于射线方法具有更高的精度，同时比有限元方法计算量小，因此在实际应用中占很重要的地位。
(一)国外研究现状有限差分法数值模拟技术的发展只有几十年的历史，早期对于地震波传播有限差分方法的研究都是基于二阶位移方程。
有限差分数值模拟方法最早由Alterman和Karal(1968)提出，他们首次将有限差分数值模拟方法应用于层状介质的模拟，并分析弹性波在层状均匀介质中的传播规律。
他们为利用差分方法解决各种勘探地震学的实际问题奠定了基础，随后差分方法被广泛应用于地震勘探中并在应用中不断得到发展。
Boore(1972)进行了非均匀介质地震波有限差分数值模拟并研究了地震波非均匀介质中的传播规律。
Alford (1974)对有限差分模拟的影响因素(网格大小和地震波传播方向)进行分析，网格间距控制地震波数值模拟的计算精度，相同精度下高阶差分与低阶差分对网格间距的要求不同型;Kelly(1976)利用有限差分法进行地震记录的合成，还对合成的地震记录进行一些常规的数据处理分析，进一步拓宽了数值模拟在实际地震处理中的应用；前面都是基于笛卡尔坐标系下常规网格的离散差分，Madariaga(1976)首次提出交错网格方法，并将此方法应用于弹性介质地震波的模拟，其模拟的差分精度为O(Δt2+Δx2)。