基于三阶Adams格式的求解声波方程的多步算法

合集下载

二维地震波场的组合型紧致有限差分数值模拟

二维地震波场的组合型紧致有限差分数值模拟汪勇;石好果;周成尧;桂志先【摘要】地震波场数值模拟在地球物理勘探和地震学中具有重要的支撑作用.本文将组合型紧致差分格式用于声波和弹性波方程的数值模拟中.根据泰勒级数展开和声波方程,建立了位移场时间四阶离散格式,并将组合型紧致差分格式用于位移场空间导数的求取,然后对该差分格式进行了精度分析、误差分析、频散分析和稳定性分析.理论研究结果表明:① 该差分格式为时间四阶、空间六阶精度,与常规七点六阶中心差分和五点六阶紧致差分相比,具有更小的截断误差和更高的模拟精度;② 每个波长仅需要5.6个采样点,且满足稳定性条件的库郎数为0.792,可以使用粗网格和较大时间步长进行计算.所以该方法具有占用内存少、计算效率高和低数值频散等优势.最后,本文进行了二维各向同性完全弹性介质的声波和弹性波方程的数值模拟,实验结果表明本文提出的方法具有更高的计算精度,能够大幅度的节约计算量和内存需求,对于三维大尺度模型问题具有更好的适应性.【期刊名称】《地球物理学报》【年(卷),期】2018(061)011【总页数】16页(P4568-4583)【关键词】组合型紧致差分;声波方程;弹性波方程;数值模拟;数值频散;稳定性条件【作者】汪勇;石好果;周成尧;桂志先【作者单位】油气资源与勘探技术教育部重点实验室(长江大学),武汉 430100;长江大学地球物理与石油资源学院,武汉 430100;胜利油田勘探开发研究院西部分院,山东东营 257001;油气资源与勘探技术教育部重点实验室(长江大学),武汉 430100;长江大学地球物理与石油资源学院,武汉 430100;油气资源与勘探技术教育部重点实验室(长江大学),武汉 430100;长江大学地球物理与石油资源学院,武汉 430100【正文语种】中文【中图分类】P6310 引言地震波场数值模拟是应用数值方法模拟地震波在地下介质中的传播过程，计算在地面或地下各观测点地震记录的一种数值模拟方法，是地球物理勘探和地震学的重要基础.目前，常用的数值模拟方法主要有射线追踪法和波动方程法，其中波动方程法有伪谱法、有限元法、边界元法、谱元法和有限差分法(Alterman and Karal，1968；Graves，1996；董良国等，2000a；Saenger et al.，2000，2004).随着数值模拟技术的发展和生产实践的要求，围绕着提高有限差分计算效率(黄超和董良国，2009；唐佳等，2016)、模拟精度(Liu，2013；李振春等，2016)、算法稳定性(Virieux, 1986；董良国等，2000；杜启振等，2015)、压制数值频散、处理复杂介质(Fletcher et al.，2008；王璐琛等，2016；程玖兵等，2013)和吸收边界条件(Berenger，1994；Pan et al.，2012)等方面，研究人员已经发展了多种方法.特别是在如何压制数值频散和提高计算效率方面，杨顶辉等在这方面做了大量相关工作(Yang et al.，2009，2012a,b，2014；Yang and Wang, 2010，Ma and Yang，2017；Tong et al.，2013；Zhou et al.，2015；He et al.，2015；Liu et al.，2017；Jing et al.，2017)，取得了许多有意义研究成果.提高差分格式数值计算精度最直接的方法就是在差分计算时增加网格节点数，但这也增加了计算量和所需的存储空间.紧致差分方法恰好能够较好地解决这个矛盾，同时紧致差分是一种隐式差分格式，具有较好的稳定性，这些优势也使得它成为目前研究较多的有限差分方法之一.1989年，Dennis和Hudson(1989)针对Navier-Stokes方程提出了空间四阶的紧致格式，能够使用粗网格获得较高的精度.1992年，Lele(1992)构造了求解一阶导数和二阶导数的紧致差分格式.Adams 和Shariff(1996)提出了紧致ENO格式，用于求解激波湍流相互作用问题.Chu和Fan(1998)提出了三点组合型紧致差分格式，并将其用于求解对流扩散方程.随后人们针对不同的问题和方程，提出了许多种不同的紧致差分格式，广泛用于空气动力学、流体力学和电磁波方程等方面.在地震波场数值模拟方面，王书强等(2002)和Zhou和Greenhalgh(2003)先后将Lele(1992)提出的紧致差分格式用于弹性波方程的数值模拟中.Yang等(2003)年给出各向异性情况下的紧致差分方法以及声波和弹性波数值模拟结果，并给出相应的吸收边界条件.Du等(2009a,b)利用紧致交错网格差分方法对横向各向同性介质地震波场进行了数值模拟，Liu和Sen(2009)研究了任意阶空间导数的紧致格式和一阶空间导数的交错紧致格式，分析了它们的精度和计算效率，并将其用于声波和弹性波数值模拟中.Kosloff等(2010)提出可以利用当前点两边各任意多个点计算导数值的一般紧致格式，并基于Remes方法求取差分系数.杨宽德等(2011)研究了Biot流动和喷射流动耦合作用下波传播的FCT紧致差分模拟方法.此外，Chu和Stoffa(2010)、Liu(2014)还研究了频率域地震波数值模拟中的紧致差分方法.而组合型紧致差分格式在地震波场模拟中的应用，尚未见到文献报道.本文在声波和弹性波方程时间四阶离散格式的基础上，将组合型紧致格式应用到位移场空间偏导数的求取，实现了二维各向同性介质地震波场的数值模拟.1 组合型紧致有限差分方法原理早期的紧致差分格式是基于Hermite多项式构造而来的，在Hermite多项式中，相邻三个节点的函数值、一阶和二阶导数值关系可以表示为其中fi表示节点函数值，和分别表示一阶和二阶导数值.1992年，Lele(1992)对Hermite公式进行了扩展，构造了求解一阶导数和二阶导数的紧致差分格式，表示为(2)其中a,b,c,α,β为差分系数，h表示网格间距.对上述Lele差分格式进行泰勒公式展开和求解差分系数，可以得到一阶导数五点中心六阶精度差分格式：(3)二阶导数五点中心六阶精度差分格式：(4)由公式(3)和(4)可以看出，紧致差分(Compact Difference，以下简称CD)格式的特点就是用相邻节点的函数值求解若干节点上的导数值.而常规中心差分(Central Finite Difference，以下简称CFD)仅求解中心节点的导数值.另一方面，也可以认为常规中心差分是紧致差分的特例，即公式(2)中的α=β=0，例如一阶和二阶导数的六阶精度差分格式表示为(5)其中b1=0.75, b2=-0.15, b3=0.0167, c0=-2.7222, c1=1.5, c2=-0.15,c3=0.0111.常规中心差分如果要得到六阶精度的差分格式，需要用到七个节点，而紧致差分只需要五个点，这使得紧致差分格式比常规差分处理边界节点上更为方便.且一般情况下，在相同网格间距时，紧致差分格式比常规差分具有更高的精度，具有更小的数值频散(王书强等，2002).Chu和Fan(1998)等构造了精度更高的三点六阶组合型紧致差分格式(Combined Compact Difference，以下简称CCD)：(6)与公式(3)和(4)对比可以看出，上述组合型紧致差分格式只需要相邻的三个节点就可以同时求得一阶和二阶导数的六阶精度近似值，比常规紧致差分的节点数更少.公式(6)中的一阶导数和二阶导数是耦合的，既可以同时求出也有利于波形的保真性.2 声波方程组合型紧致有限差分方法分析二维声波方程可以表示为(7)式中u(x,z)为地震波位移，v(x,z)为地震波速度.利用截断的泰勒公式表示n+1和n-1时刻的位移场可以得到：(9)两式相加，略去高次项，并代入声波方程，得到位移场时间四阶精度的差分格式：(10)公式(10)为位移场的三层显式差分格式，利用它就可以进行地震波场时间层的推进.在公式中含有位移对x和z的二阶和四阶导数，对这些导数将利用组合型紧致差分格式(6)进行求取.假设模型纵向和横向节点数为m和n，h为空间网格大小.利用公式(6)求公式(10)中的空间偏导数∂2u/∂z2和∂2u/∂x2，矩阵表示为AE=BU， FC=UD,(11)其中A和C为公式(6)左端的差分系数矩阵，大小分别为2m×2m和2n×2n；E和F为待求位移场空间一阶和二阶导数矩阵，大小分别为2m×n和m×2n；B和D为公式(6)右端的差分系数矩阵，大小分别为2m×m和n×2n；U为位移场矩阵，大小为m×n.这些矩阵分别表示为(12)(13)(15)(16)(17)(18)公式(11)可以用追赶法进行求解.由公式得出位移场ui,j可同时求得它的空间一阶和二阶导数∂ui,j/∂x、∂ui,j/∂z、∂2ui,j/∂x2和∂2ui,j/∂z2，表示为E=A-1BU， F=UDC-1,(19)差分公式(10)中的空间四阶导数可以在求得二阶导数后，将二阶导数值当作U，再次使用公式(19)进行求取.2.1 精度分析不论是利用CCD格式还是利用常规的七点六阶CFD或五点六阶CD格式，进行声波方程数值模拟时，它们在时间层推进方式上是相同的，即都是利用差分格式(10)，时间差分为四阶精度.CD、CFD和CCD不同之处在于它们分别利用公式(4)、(5)和(6)计算位移空间高阶导数.虽然在声波方程的差分格式中并没有空间一阶偏导数，但在下文中的第4节进行弹性波场数值模拟中，存在空间的混合偏导数，需要通过一阶导数的迭代求取，所以这里对这三种格式计算的一阶和二阶导数的截断误差进行对比，结果见表1，α、β、a、b和c为公式(2)中的差分系数.从表1可以看出，虽然三种方法都能达到空间六阶精度，但截断误差有较大的差别.CFD和CD差分格式计算一阶偏导数的截断误差约是CCD的40倍和4.4倍，计算二阶偏导数的截断误差约是CCD的36倍和8.5倍，这说明CCD方法的具有更高的差分精度.表1 一阶和二阶导数不同差分方法的截断误差比较Table 1 Truncation errors in various difference schemes for the first and second derivative calculationsαβabc截断误差uxCFD003/2-3/51/1036/(7!)×f(7)h6CD1/3014/91/904/(7!)×f(7)h6CCD/≈-0.9/(7!)×f(7)h62ux2CFD003/2-3/51/1072/(8!)×f(8)h6CD2/11012/113/110≈-17/(8!)×f(8)h6CCD/-2/(8!)×f(8)h62.2 误差分析通过模拟二维平面谐波初值问题并计算模拟误差，定量分析CCD方法的数值模拟精度.二维平面波初值问题可以表示为(20)其中v是平面波的波速，θ 0是初始时刻波阵面法线方向(即传播方向)与x轴的夹角，f0是平面简谐波的频率.上述初值问题的解析解为(21)在二维波场数值模拟中，假设f0=20 Hz，θ0=π/4，均匀介质的波速为3600 m·s-1，模型长度和深度均为2000 m，纵横向网格长度相同，采样时间为1 s.在不同的空间网格长度和时间步长条件下，计算数值模拟的相对误差.相对误差定义为(22)式中表示数值解，u(tn,xi,zj)表示解析解.研究表明(Yang et al.，2014)，优化的近似解析离散化方法(Optimal Nearly-Analytic Discrete，以下简称ONAD)也具有较高的数值模拟精度和计算效率，为了比较，在此计算CCD、CD、ONAD和CFD四种方法在不同模拟参数情况下的相对误差曲线见图1.从图1可以看出，随着空间网格长度和时间步长的增加，四种方法的相对误差均会逐渐增加，并且随模拟时间的增加而增加.便于比较，将实验结果列于表2.从表中可以看出，四种情况下，误差相对关系基本一致，例如当空间步长为20 m，时间步长1 ms时，CCD算法的相对误差仅为0.066%，CD次之，为0.089%，CFD相对误差为0.943%，这与精度分析的结论是一致的.需要说明的是，本次实验的ONAD算法相对误差比CFD方法大，这是由于ONAD只有空间四阶精度，而本文的CFD空间精度为六阶.当采用较细网格时，如空间步长10 m，时间步长0.5 ms时，四阶ONAD方法的误差要小于六阶CFD方法，这也说明了ONAD方法在提高模拟精度方面也有积极的意义.用较小的空间步长(10 m)和时间步长(0.5 ms)时，CCD算法的精度显著提高，相对误差仅有0.0008%，并且相对误差增长缓慢，这说明了采用CCD格式的声波模拟结果精度较高，能进行较长时间的地震波场数值模拟.图1 四种数值模拟算法的相对误差曲线(a) CCD; (b) CD; (c) ONAD; (d) CFD.Fig.1 Relative error-time curves of the numerical simulation using four different methods and models(a) CCD; (b) CD; (c) ONAD; (d) CFD.表2 四种方法在不同情况下的最大相对误差(%)比较Table 2 Relative errors in four different schemes and parameters模拟参数CCDCDONADCFDdx=10 m,Δt=0.5ms0.00080.00160.1190.596dx=15 m,Δt =0.5 ms0.0120.0150.690.62dx =15m,Δt =1 ms0.0080.0190.5560.629dx =20 m,Δt =1 ms0.0660.0892.160.943 2.3 频散分析在数值模拟时，如果空间网格长度使用过大，会造成较大的求解误差，并会产生所谓的数值频散现象，所以频散关系分析既是判断数值模拟方法优劣的重要方法，也是确定空间网格大小的重要依据.通过对声波方程的组合型紧致差分格式进行数值频散分析，进一步分析该方法的适用条件.首先考虑x方向，令简谐波解=Aexp×[I(ihkx+jhkz-kv′nΔt)]，且：(23)如果数值模拟时不存在误差，则A,B,C应该满足以下关系：B=(Ikx)A, C=-(kx)2A,(24)但在差分的实际计算中，所产生的数值波数(又称修正波数)与解析波数不同.令：(25)将公式(23)和(25)代入差分格式(6)，可以得到方程组：(26)其中：解上述方程组可得：(27)同理z方向可得：(28)式中φx=φcosθ，φz=φsinθ，且φ=kh.θ是波的传播方向与x轴的夹角，k为波数，h为空间步长.将公式(23)、(27)和(28)代入声波方程的差分格式(10)中，并令x=kv′Δt，α=vΔt/h.其中v′为数值波速，v是实际波速，Δt是时间步长，α称为Courant数，利用欧拉公式，化简可得CCD频散关系：]]/12.(29)通过上述频散关系确定φ后可以解得对应的x.定义数值波速与真速度的比值为(30)在压制数值频散方面，优化的近似解析离散化方法(ONAD)方法也具有较好的效果，所以选取CFD、CD和ONAD 三种方法与之进行比较.在理想情况下，如果不存在数值频散则速度比γ恒等于1.γ偏离1越大，则说明该方法的数值频散越严重，反之则说明该方法能更好地压制数值频散.图2显示了以上四种方法在不同的α和θ的时候的速度比曲线.取φ∈[0,π]作为横坐标，它是波数与空间步长的乘积，单位波长内采样点数N=2π/φ，所以横坐标也可以看作N由逐渐减小至2.从图中的速度比曲线可以看出：①随着空间采样点数的减小，四种方法的频散现象逐步加剧.CCD、CD和ONAD方法的数值频散比CFD方法要小，其频散曲线更趋近于1，这也说明了紧致差分(CCD和CD)和近似解析离散类(ONAD)方法在压制数值频散方面的都具有较大的优势，都适合使用较大的空间步长，进而提高计算效率；②假设数值速度在理论速度的99.9%～100.1%以内表示不存在频散，则CCD、CD、ONAD和CFD 对应的φ的极小值分别为1.12、0.98、1.06和0.91，所以它们在最小主波长内需要使用的网格点数分别为5.6、6.4、5.9和6.9个，即CCD和ONAD方法对于空间网格长度的要求最低，其次是CD方法，CFD方法最为严格；③CCD方法的空间网格长度过大时，速度比曲线上翘，数值速度大于真速度，这与另外三种方法相反.这一结论也能从图3c中可以看出，当φ=π，即一个周期内取2个样点，蓝色的速度比曲线大部分要大于1.图2 不同差分方法的速度比曲线(a) θ=0,α=0.25; (b) θ=π/8,α=0.25; (c)θ=π/4,α=0.25; (d) θ=0,α=0.45; (e) θ=π/8,α=0.45; (f) θ=π/4,α=0.45.Fig.2 Velocity ratio curves of the numerical simulation in different methods图3 不同差分方法的方位角-速度比曲线(a) CFD; (b) CD; (c) CCD; (d)ONAD.Fig.3 Azimuth-speed ratio curves of the numerical simulation in different methods图4 空间网格32 m，时间步长1 ms，420 ms时刻波场快照(a) CFD ; (b) CD; (c) CCD; (d) ONAD.Fig.4 Snapshots at 440ms: Δx=32 m，Δt=1 ms为研究该差分方法的数值各向异性特征，即速度比与传播方向的关系，绘制极坐标速度比曲线，如图3所示.图中的极角表示地震波传播方向与x轴的夹角，极径表示速度比.四种方法在空间采样(φ=π/3)满足频散要求时，都没有各向异性特征.在空间采样不足时(φ=7π/5，φ=9π/5)均会出现速度各向异性，0°(或90°)与45°方向上的速度差别最大.利用二维模型来说明CCD方法在压制数值频散方面的优势.模型长度和深度均为4 km，纵波速度3600 m·s-1，震源子波为f(t)=sin(2πf0t)exp×主频20 Hz.根据频散分析结论，CCD方法在最小主波长内只需要使用约5.6个点，所以设置空间网格间距32 m，时间步长1 ms.图4为四种方法数值模拟的波场快照，可以看出CCD和ONAD方法的结果清晰，没有频散现象，CD方法在0°和90°方向上存在微弱的频散现象，而CFD的频散非常严重，与前面理论分析结果一致.2.4 稳定性分析稳定性条件也是有限差分数值模拟中一个非常重要的问题，是影响差分方法计算效率的重要因素.采用Fourier方法(张文生，2006)对声波方程的组合型紧致差分格式进行稳定性分析.定义：(31)其中为空间网格大小，kx和kz为x和z方向的视波数，i、j和n分别为空间和时间网格坐标，代入公式(6)中可得：(32)令θ1=kxh，-π/2≤θ1≤π/2，并利用欧拉公式，将上式化简可得：(33)解方程可得：=Aξn/h2,(34)其中(35)同理，定义：=eIh(kxi+kzj), =eIh[kxi+kz(j-1)],=eIh[kxi+kz(j+1)], =eIh(kxi+kzj),=eIh[kxi+kz(j-1)],=eIh[kxi+kz(j+1)],且令θ2=kzh， -π/2≤θ2≤π/2，代入公式(6)可以得到：=B ξn/h2,(37)其中(38)公式(10)是一个三层显式差分格式，为了分析其增长矩阵，将其改写为：(39)式中=φneIh(kxi+kzj).将公式(34、35、37和38)代入(∂2u/∂x2=exp[Ih(kxi+kzj)]和(∂2u/∂z2=exp[Ih(kxi+kzj)]，并代入声波差分格式(39)中，化简并写成矩阵格式为(40)其中(41)为增长矩阵，其中α=vmaxΔt/h为该差分格式的Courant数.为使差分格式(40)满足稳定性条件，需要其增长矩阵的谱半径小于等于1.直接求解该增长矩阵的特征值较为复杂，用二分法求该问题的解(宋国杰，2008)，计算得到的近似稳定性条件为α=vmaxΔt/h<0.792.2.5 计算效率分析采用2.3节所用模型进行数值模拟，采样时间1 s，分析CCD、CFD、CD和ONAD四种方法的计算效率，各方法使用的数组大小、个数和算法时间复杂度等列于表3.表3 不同方法计算效率比较Table 3 Comparison of computational efficiency between different methodsCCDCFDCDONAD网格半径3753每个波长采样点5.66.96.45.9Δx=Δz/(m)32262831数组/个1281225Courant0.790.980.930.527Δt/ms7774时间复杂度O(n3)O(n3)O(n2)O(n2)时间/s0.6670.8681.041.79从表3中可以看出：CCD和ONAD差分方法的网格半径最少，仅需要3个点，这对于边界处理最有利，但ONAD空间精度仅为四阶，要低于CCD的六阶精度.CFD和CD差分方法虽然空间精度能达到六阶，但网格半径也随之增加.CCD和ONAD均为低数值频散的差分方法，满足无数值频散要求的空间网格长度大于CFD和CD方法.但是ONAD方法需要计算位移场的梯度，这就增加了数组的个数和计算时间.同样的声波方程，CCD仅需要12个数组，而ONAD需要25个，在网格长度相近的情况下，CCD方法占用内存最小.从稳定性来看，要达到同样的时间四阶精度，CFD和CD方法的Courant数最大，CCD次之，ONAD最小.考虑空间步长的不同，CCD、CFD和CD能允许较大的时间步长(7 ms).从算法的时间复杂度来看，CCD和CD算法的复杂度要大于CFD和ONAD，其中n表示网格点数.值得注意的是，在模型大小相同的时候，由于网格长度的差异，会造成四种算法的空间网格节点数n不同，也就造成了实际所需运算时间的不同，表中运算时间是在同一台计算机和相同编程环境下模拟得到的.综合模拟精度、空间网格长度、时间步长和占用内存四个因素来看，CCD方法进行声波方程数值模拟时，它具有高精度(时间四阶、空间六阶)、高计算效率(较大Courant数、占用内存少)和低数值频散(粗网格)的优势.需要说明的是，表3中四种方法各参数的选取时，均不考虑计算精度，空间网格和时间步长仅满足无数值频散和稳定性的要求.但在实际计算时，为了精度考虑，二者均需要按要求适当减小，这势必会进一步增加计算量和存储量.3 弹性波方程组合型紧致有限差分方法各向同性完全弹性介质中的二维波动方程可以表示为(43)其中，u(x,z)和w(x,z)表示x和z方向的质点位移，VP(x,z)和VS(x,z)表示纵横波波速.同声波方程一样，利用公式(8)和(9)，并代入波动方程，将u对时间的导数转换为u和w对空间的导数，可以得到时间四阶精度差分格式：(44)同理可以得到w的时间四阶精度差分格式为(45)弹性波差分格式与声波方程不同，公式(44)和(45)中除了位移u的空间二阶导数∂2ui,j/∂x2和∂2ui,j/∂z2外，还有二阶混合偏导数∂2ui,j/∂x∂z和其他四阶偏导数.假设偏导数∂2/∂x2，∂2/∂z2，∂/∂x和∂/∂z的差分算子为和则二阶混合偏导数∂2/∂x∂z的差分算子为四阶偏导数∂4/∂x4的差分算子表示为即对于二阶混合偏导数，可以在求出∂ui,j/∂x和∂ui,j/∂z之后，再次使用公式(19)，对另一个方向求导，从而求出空间二阶混合偏导数∂2ui,j/∂x∂z.同样，对差分公式中出现的u和w的其他高阶偏导数，也可以按照上述同样的方法进行迭代求取，这里不再赘述.4 模型试算为了验证方法的数值性能，我们采用组合型紧致方法对公式(43)和(44)表示的二维各向同性完全弹性介质的弹性波方程进行数值模拟.4.1 均匀介质模型模型长度和深度均为10 km，介质为泊松体，纵波速度3460 m·s-1，横波速度2000 m·s-1，在模型中心处激发10 Hz主频的雷克子波.由于最小速度为2000 m·s-1，按频散分析的结果，每个波长内需采样5.6个点，所以设置网格长度为35 m.同时，由于最大速度为3460 m·s-1，按稳定性条件，取时间步长为8 ms.长时程数值模拟是数值模拟中的一个重要内容(Chen,2006; Li et al.,2012; Gao et al.,2016).对于大尺度模型和长时程传播，即便是十分微弱的误差也会随着传播时间的增加而不断累积，最终可能导致不容忽视的数值假象.为了分析CCD方法的长时程数值模拟结果，数值模拟时没有使用吸收边界条件，图5(a—d)显示的是1 s、2 s、5 s和10 s时刻的位移水平分量波场快照，地震波场特征清晰，无可见的数值频散现象.为了比较，对该模型采用相同的参数，使用时间四阶空间六阶精度CFD格式进行数值模拟，图5(e—h)为相同时刻的波场快照.可以看出，在满足频散要求和稳定性条件的情况下，CCD方法具有计算效率较高、长时程稳定的特点，能适应较长时间的数值模拟，计算结果较为理想.在相同参数下，CFD结果存在很明显的频散现象，并且随时间的增加，频散更加明显，严重地影响了数值模拟的效率和精度.图5 均匀介质模型弹性波数值模拟波场快照(a) CCD(1 s); (b) CCD(2 s); (c)CCD(5 s); (d) CCD(10 s); (e) CFD(1 s); (f) CFD(2 s); (g) CFD(5 s); (h) CFD(10 s).Fig.5 Wave field snapshots using elastic wave equation in homogeneousmedium(a—d) are the snapshots at 1,2,5,10 s for the displacement component in X direction using CCD; (e—h) are the snapshots at 1,2,5,10 s for the displacement component in X direction using CFD.4.2 水平层状介质模型设置水平层状介质模型，共20层，每层厚度400 m，模型总厚度8 km，模型长度同样是8 km.第一层纵波速度为3500 m·s-1，往下各层纵波速度比上一层增加100 m·s-1.模型为泊松体，即纵横波波速比为1.732，所以模型中最小和最大波速为2020 m·s-1和5400 m·s-1.在地面中间处(X=4 km，Z=0 m)激发10 Hz雷克子波震源，人工边界采用PML边界条件(刘有山等，2012)处理.按前文分析结果，设置纵横向网格长度为36 m，时间步长为5 ms.图6显示的是地面接收到的地震记录，记录长度8 s，由于直达波能量明显强于反射波，所以图中对地震记录进行了增益处理.多层模型中地面模拟接收到的地震记录非常清晰，没有数值频散和不稳定数值结果，地震记录中的直达波、反射波和转换波显示清楚，说明组合型紧致有限差分算法可以有效地模拟弹性波在多层各向同性介质模型中的传播.为了检验本文算法长时程数值模拟的结果，图7给出在较粗网格条件下(网格长度36 m，时间步长5 ms)，CCD和CFD方法分别得到的检波点(X=4 km，Z=0 m)处的波形图，并与利用CFD方法在较细的网格(网格长度10 m，时间步长1 ms)情况下得到的近似理论解的参考波形进行对比.图7a显示三种情况下得到的反射波到达时间一致，从图7(b，c)的局部放大结果来看，在粗网格和长时程模拟条件下，CCD方法得到的地震记录平滑无数值频散，与之前的分析结果一致，且CCD方法和近似解析解的波形曲线基本重合，这说明了CCD方法的正确性.而CFD方法在相同条件下存在较严重的数值频散，且模拟时间越长，模拟结果与近似解偏差越大.这些结果和对比说明，CCD方法能够长时间计算稳定性，而且也能很好地压制数值频散.5 结论与建议有限差分法是地震波场数值模拟的重要方法之一，在目前常规有限差分数值模拟中为了获得更高的模拟精度，需要使用更多的节点，这样会增加计算量和所需的存储空间，还会增加人工边界处理的难度.本文从组合型紧致差分格式出发，建立了时间四阶、空间六阶的声波方程的差分格式，对该差分格式进行了模拟精度、频散曲线和稳定性分析，并与常用的几种格式进行了综合比较，该格式可同时计算空间一阶和二阶导数，且具有使用节点少(3个)，低数值频散(每个波长采样5.6个点)，较高稳定性(Courant数为0.792)和模拟精度(空间6阶)的特点，适用于大尺度和较长时程地震波场数值模拟.文中还建立了各向同性介质弹性波方程的差分格式，并进行了数值模拟，模拟结果显示地震波场特征清晰，说明了该方法的适用性.该方法不仅可以用于弹性介质的波场模拟，还可以进一步推广到二维或三维的各向异性介质、黏弹介质和双相介质等复杂介质的数值模拟中，为今后研究地震波传播规律、逆时偏移、全波形反演等工作提供了一种有效的方法.组合型紧致差分格式仅需使用三个节点，就能使一阶和二阶空间偏导数的差分精度达到传统差分方法的六阶精度，因而能够大幅度节约内存需求和计算量，对于三维大尺度模型的正演和反演等问题均具有重要意义.图6 水平层状介质模型弹性波方程数值模拟地面地震记录(a) X方向位移分量； (b) Z方向位移分量.Fig.6 Seismic records of the model′s surface using elastic wave equation in horizontal layers media(a) Displacement component in X direction； (b) Displacement component in Z direction.图7 接收点(4 km，0 km)处接收到的反射波水平位移波形图(a) 整道反射波记录对比； (b) 4000～4500 ms局部放大对比； (c) 6000～6500 ms局部放大对比.其中，红色加点实线是粗网格条件下的CCD数值结果，绿色实线是粗网格条件下。

微分方程和偏微分方程的数值解法

期权定价模型
描述金融衍生品的定价过程，如布莱克-舒尔斯模型就是一个偏微分方程。通过求解该方程，可以得到期权的理论价格以及相应的风险参数。
投资组合优化
在投资组合理论中，常使用微分方程来描述资产价格的动态变化和投资者的风险偏好。通过求解这些方程，可以得到最优的投资组合配置策略以实现风险与收益的平衡。
数值解法需要保证稳定性和收敛性，即当离散间隔趋近于零时，数值解应趋近于真实解。
02
常微分方程的数值解法
欧拉方法
基本思想
通过逐步逼近的方式，利用已知点的信息来推测下一个点的信息。
公式推导
基于泰勒级数展开，忽略高阶项得到近似公式。
优缺点
简单易懂，但精度较低，仅适用于简单问题。
改进方法
采用改进的欧拉方法或预估-校正法提高精度。
物理问题中的微分方程和偏微分方程
牛顿第二定律
描述物体运动的基本定律，可以表示为二阶常微分方程。通过求解该方程，可以得到物体的位移、速度和加速度等运动学量。
热传导方程
描述热量在物体内部传递的过程，是一个偏微分方程。通过求解该方程，可以得到物体内部的温度分布以及热量的传递速率。
波动方程
描述波动现象（如声波、光波等）的传播过程，是一个二阶偏微分方程。通过求解该方程，可以得到波的传播速度、振幅、频率等波动特性。
工程问题中的微分方程和偏微分方程
结构力学中的弹性力学方程
描述结构在受力作用下的变形和应力分布，是一个偏微分方程。通过求解该方程，可以得到结构的位移、应力和应变等力学量，为工程设计提供重要依据。
流体力学中的纳维-斯托克斯方程
描述流体运动的基本方程，是一个偏微分方程。通过求解该方程，可以得到流体的速度、压力和温度等流场特性，为流体机械设计和优化提供指导。

基于MATLAB／ADAMS的平面三自由度并联机构的运动学和动力学分析及控制的初步设计

度，且以并联方式驱动的一种闭环机构。并联机
构由于具有累积误差小、运动惯量低、负载能力
１动学分析运
１１．运动学数学模型的建立
强、刚度较大等特点，已成为一种潜在的高速度、高精度运动机构ｎ。反解法是一种已知机构工作部分的运动情况而逆向推导主动件运动情况的研究方法，并联机
Ａｂｓｒｃ：Ｔｈｎｍａｉｓａｄｙａｃａａｙｉｆａｐａａ－ｔａｔｅｋｉｅｔｎｄｎｍｉｓｎｌｓｓｏｌｎｒ３ＤＯＦａａｌｌｃｐｒｌｅｍｅｈｎｓｃａｉｍｉｏｄｃｅｔｒｕｈｓｃｎｕｔｄｈｏｇＭＡＴＬＡＢ，ａｄｈｎｔｅｃｌｕａｉｎｒｓｌｓａｅｓｏｄｉｕｖｏｍ．Ａｆｅｈａ，ｔｒｕｈａｓｍｕａｉｎｉａｃｌｔｏｅｕｔｒｈｗｅｎａｃｒｅｆｒｔｒｔｔｈｏｇｉｌｔｏｎＡＤＡＭＳｔｅｃｌｕａｉｎｒｓｌｓｖｌｄｔｄｗｉｈｉｓ，ｈａｃｌｔｏｅｕｔｉａｉａｅｔｔｓａｃｒｃ，ｎｅｋｎｅｔｑａｉｎｔｄｆｏｔｅｄｔａｃｕａｙａｄｔｉｍａｉｅｕｔｓｆｔｒｍｈａａｔｔＭＡＴＬｈｃｏｉｅｈＡＢａｃｌｔｄａｅｔａｓｌｎｅｎｏＡＤＡＭＳｔｉｌｔｈｏｔｏｆｃｌｕａｅｒｒｎｐａｔｄｉｔｏｓｍｕａｅｔｅｃｎｒｌｏｔｉｃａｉｍ．ｓｌｏｂｏｈｍｅｈｄａｅｌｉｈｏｎａｉｎｆｒｔｅｓｂｅｕｅｔｄｖｌｐｎｆｈｏｒｌｙｔｍ．ｈｓｍｅｈｎｓＲｅｕｔｆｍｔｔｏｓｈｖａｄｔｅｆｕｄｔｏｕｓｑｎｅｅｏｍｅｔｏｅｃｎｔｓｅｓｒｏｈｔｏｓＫｅｗｏｄ：ｐａａ３ｙｒｓｌｎｒ－ＤＯＦａａｌｌｐｒｌｅｍｅｈｎｉｍ；ＭＡＴｃａｓＬＡＢ／ＡＤＡＭＳ；ｋｎｍａｉｓｎｌｓｓ；ｄｎｍｉｓｎｌｓｓｉｅｔｃａａｙｉｙａｃａａｙｉ；ｍｏｅｉｇｎｄｌｎａｄｓｍｕａｉｎ；ｐｒｌｎｒｏｔｏｅｉｎｉｌｔｏｅｉａｙｃｎｒｌｓｇｍｉｄ

基于 ADAMS 的振动过程频率特性分析

基于ADAMS的振动过程频率特性分析李晓静1，杨丰翔2，刘保军3【摘要】摘要：机械系统动力学分析软件ADAMS的后处理模块，能够帮助实现模型调试、试验验证、设计方案改进和结果显示功能，从而便于从可视化角度深入研究设计的有效性。

以一个简单多体动力学模型为例进行振动分析，采用ADAMS PostProcessor进行数据的后处理，研究仿真分析过程。

利用FFT曲线图进行分析，发现动力学模型的加速度频率特性中谱密度幅值的峰值发生在前19～20 Hz处，据此可以研究模型系统的性能。

【期刊名称】新技术新工艺【年(卷),期】2015(000)009【总页数】3【关键词】模型；处理；频率特性虚拟样机仿真是以并行工程思想为指导，建模仿真理论为核心，以各领域计算机辅助仿真软件为工具，进行产品各种性能测试和评估的过程。

机械系统动力学分析(Automatic Dynamic Analysis of Mechanical System，ADAMS)软件有助于工程技术人员快速建立机械系统虚拟样机，分析其性能，更好地理解机械系统的运动[1-3]。

通过分析各种设计方案，精确表示出载荷的变化，计算其运动模型的轨迹、速度和加速度状况等。

虚拟样机在工程中的应用是通过虚拟样机软件实现的。

比较有影响的软件是美国MSC公司的ADAMS、比利时的DADS以及德国的SIMPACK。

其中，美国MSC公司开发的ADAMS软件占据市场50%以上份额。

ADAMS软件领先的“功能化虚拟样机”技术，已经开始应用于汽车、航空、铁道、兵器、核能及工程机械等领域。

ADAMS PostProcessor是ADAMS软件的后处理模块，可以用来处理仿真的结果，动态显示仿真过程，以及完成曲线编辑和数字信号的处理功能。

后处理模块以可视化形式深入分析设计方案的可靠性，实现树状结构的搜索，层次清晰，可用于分析检索模型对象[4-5]。

1 曲线图的处理ADAMS软件配置了若干工具包用于对曲线图的处理。

第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动力学仿真算法及参数设置分析

J aco bi矩阵 , 可算出 qk + 1 , uk + 1 ,λ k + 1 , qk + 1 , u k + 1 ,λ k+1 , 重复迭代校正步骤 , 直到状态向量满足收敛条件 , 使系统方程 g ( y n + 1 , y n + 1 , t n + 1 ) = 0 近似成立。 ( 3) 检查积分误差并优化积分步长和多项式阶
( 8)
得到图 1 。 1. 机械系统动力学方程的秩式 ( 1) ～ ( 3) 是代数微分方程组而不是常微分方程组[ 5 ] ,即方程中的导数项不能化为显式的形式 ( 常微分方程) 。由式 ( 11) 可以看出 , 当步长很小时 ,J a2 co bi 阵将出现很大的项 , 这在数学上称为病态。病态的 J aco bi 阵将使求解线性方程组时产生严重的舍入误差 ,这是由代数微分方程自身的特点决定的。式 ( 1) 可以转化为式 ( 12a ) 的形式 , 在数学上称这样的方程为指数 3 ( index - 3 ) 公式[ 5 ] 。所谓指数是指 : 将代数微分方程转化为常微分方程所需要的微分的次数。 T λ- A T F = 0 Mu + φ q
1
I I
( 9) ( 10)
0 0 ( c)
・ φ ( q , u , t) = 0 φ( q , t) = 0 T・ T Mu + φ qη - A F = 0
ζ= 0 u - q +φ q
T・
0
k
Δq - F Δu = - G ( 11) Δ λk - Φ k 其中 :方程左边的系数为系统的 J aco bi 矩阵 ; 5 F/ 5 q 系统刚度矩阵 ; 5 F/ 5 u 系统阻尼矩阵 ; 5 F/ 5 u 系统质量矩阵式 ( 11 ) 表明 , 求解非线性代数方程组的问题已转化为求解线性方程组的问题。通过分解系统

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

基于ANSYS_LS_DYNA的应力波反射法的数值模拟_张乐婷

专家专稿文章编号:1009-6825(2010)32-0001-02基于ANS YS /LS -DYNA 的应力波反射法的数值模拟收稿日期:2010-07-25作者简介:张乐婷(1985-),女,兰州交通大学土木工程学院岩土工程专业硕士研究生,甘肃兰州 730070余云燕(1968-),女,教授,兰州交通大学土木工程学院,甘肃兰州 730070张乐婷余云燕摘要:采用一种基于AN S Y S /L S -DYNA 的非线性动力有限元分析方法,对桩身进行了数值模拟,得到了桩的应力波反射特点,验证了基于AN S Y S /L S -DYNA 非线性动力有限元分析方法在桩身应力波分析问题中的可行性。

关键词:应力波,有限元,数值模拟,AN S Y S /LS -DYNA 中图分类号:TU 473.1文献标识码:A0 引言由于桩基础可以把荷载传至稳定层,并达到了安全可靠的效果,所以随着经济的迅速发展,桩基础已被广泛用于各类建筑物、桥梁、港口等结构。

但由于桩基础作为一种隐蔽工程,难免在施工和使用阶段出现问题而不被及时发现,为了保证工程质量,必须对桩基质量进行检验,以便准确判断出缺陷的类型,测出缺陷的位置和程度,采取经济合理、易操作的补救措施,防止事故的发生。

由此可见桩身完整性检测对桩基工程而言具有极为重要的意义。

所以桩基础的检测成为众多学者研究的课题,同时也出现了很多有效的研究方法如回传射线矩阵法[1-3]、以波动理论为基础的动测技术等。

AN S Y S 是一个融结构、热、流体、电磁和声学于一体的有限元分析软件,近年来广泛用于工程领域。

其有限单元法(或称有限元法)是当今工程分析中获得最广泛应用的数值计算方法。

由于它的通用性和有效性,受到工程技术界的高度重视。

伴随着计算机科学和技术的快速发展,现已成为计算机辅助设计和计算机辅助制造的重要组成部分。

在AN SYS /LS -DYNA 中,ANSYS 仅仅为LS -DYNA 提供前后处理,具体求解过程由LS -DYNA 版求解器来完成。

工程案例—机器人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和MATLAB自行车机器人联合仿真运动

编号：毕业设计(论文)说明书题目：基于ADAMS和MATLAB自行车机器人圆周运动仿真研究学院：机电工程学院专业：机械设计制造及其自动化学生姓名：学号：指导教师单位：机电工程学院姓名：职称：副教授题目类型：☐理论研究☐实验研究☐工程设计☑工程技术研究2014年5月28日摘要自行车是我们日常很常见的交通工具，自行车机器人是在自行车的基础上，加上自动控制系统来达到不需要人为控制的自行车，简单的说它就是把自行车和机器人结合起来的一种机器。

在驱动设备的驱动下，它可以通过它的控制器进行控制自主完成运动。

我们都知道，自行车主要通过它的车轮子和地面间的摩擦力来运动，地面与车轮的接触跟点接触差不多，所以他的的平衡动力学模型是相当的复杂的。

因为机械系统的物理样机和动力学模型建立起来太过麻烦，而且成本很高。

所以在文中我们要用到虚拟样机，虚拟样机技术的好处在于不用建立物理样机，只需要用三维建模软件建立好几何模型，再给模型加上合理的约束、驱动力矩、质量属性等必要因素，这样就能简单便捷地到了我们想要的模型，大大减少了我们的工作量，使得研究起来更简单。

另外，结合Matlab软件进行联合仿真，这样一来，整个联合仿真过程就显得更加的简单了。

得助于这些软件，这次的研究联合仿真运动控制的成本和时间都节约了很多。

本文主要研究自行车机器人基于虚拟样机的圆周运动的平衡控制。

具体步骤如下；首先用所学软件建立它的几何模型。

文中本人使用的是PRO/E建立自行车机器人的三维模型。

然后装配好的自行车机器人模型保存为T_X文件以方便它导入到ADAMS2012，紧接着在动力学软件ADAMS2012中加上必要的属性从而得到带着动力学模型的机械模型。

再然后把课本学到是控制工程的思想拿出来运用，并且结合数学软件Matlab的Simulink模块进行联合仿真，实现了自行车的圆周运动任务。

关键词：PRO/E；Simulink；ADAMS ; PID控制；圆周运动；自行车机器人AbstractBicycle is our daily common means of transportation,The bicycle robots on the basis of the bike and combineding with automatic control system to achieve not need manual controling.Simply say it is a machine combines with the bicycle and robot.By the driving of the driving device, it can be done by the controller to control the independent movement.As we all know, bicycle mainly move by the friction which between the ground and its wheels.The contact of the ground and wheels is similar to point contact.So the balance of his dynamic model is quite complicated.Since the mechanical system dynamic model of physical prototype and built up too much trouble And the cost is very high.So in this paper we use virtual prototype,The advantages of virtual prototyping technology is not set up the physical prototype,it Only need to create a geometric model with three dimensional modeling software.Then added to the model of reasonable constraint, driving moment, necessary factors such as quality attributes.So you can simple and easy to have the model waht we want , greatly reducing the workload of us,It is more easy to make research.In addition, because of the combined simulation of Matlab software ,so the simulation process is more simple.To aid in the software, this study combined simulation movement control cost and time saving a lot.This paper mainly studies bicycle robot based on virtual prototype of the balance of the circular motion control,Specific steps are as follows,First of all use the software we have learn to establish its geometric model.In this paper, I use PRO/E set up bicycle 3 d model of the robot.Then bicycle assembly robot model is saved as a T_X documents for import it into ADAMS2012.Then add the necessary property in the dynamics software ADAMS2012 with a mechanical model of dynamic model is obtained.Then take out your the thinking of the control engineering you have learn to application.And combined with the mathematical software Matlab Simulink module for combined simulation, realize the circular motion of bicycle task..Keywords:PRO/E ;Simulink ; ADAMS; PID control; circling motion;bicycle robot目录摘要 (1)Abstract (2)目录 (3)引言 (4)1 绪论 (5)1.1 自行车机器人概述 (5)1.2 三维软件PRO/E软件的简介 (8)1.3 虚拟样机技术概况 (9)1.4 动力学软件ADAMS简介 (10)1.5 控制工程中的PID控制概述 (10)1.6 数学软件Matlab/Simulink简介 (11)1.7 课题背景及设计任务 (12)2 自行车机器人的建模 (13)2.1 运用三维软件PRO/E对自行车机器人进行几何实体建模 (13)2.2 运用动力学软件ADAMS对自行车机器人进行虚拟样机建模 (20)3 ADAMS和MATLAB联合起来对自行车机器人仿真 (24)3.1 自行车机器人动态平衡原理 (24)3.2 控制策略及控制器的设计 (24)3.3 联合仿真调整PID参数 (27)4 自行车机器人回转运动特性分析 (30)4.1 无干扰力下自行车机器人的圆周运动特性 (30)4.2 有干扰力的自行车机器人的圆周运动特性 (33)5 结论 (35)谢辞 (36)参考文献 (37)引言自行车机器人是人们提出来的一种比较新型的智能交通工具或者说是运输工具，它运用了机器人的智能，还有自行车的简单，将这两者巧妙结合，在机器人研究领域是一种新的概念。

基于Adams插值方法的单种群Logistic模型的数值解法

基于Adams插值方法的单种群Logistic模型的数值解法孙礼俊【期刊名称】《《蚌埠学院学报》》【年(卷),期】2019(008)005【总页数】5页(P76-80)【关键词】阿达姆斯插值方法; 单种群; 数值解【作者】孙礼俊【作者单位】蚌埠学院理学院安徽蚌埠 233030【正文语种】中文【中图分类】O175; S931单种群模型是种群生态学的基础模型，用于抽象表示种群的密度(或个体数量)随时间变化而改变的规律的模型；它常通过逻辑斯谛(Logistic)方程来解决，这种模型需要物种满足生态系统中的环境条件：存在天敌，食物、空间等资源也不足；物种增长趋势呈S形，是在环境资源等有限条件下种群成长趋势的最佳数学模型，而且应用很广。

胡喜生，范海兰等[1]应用改进逻辑斯谛模型高精度预测人口数量增长问题；任玉珑，史乐峰等[2]采用逻辑斯谛模型预测了国内汽车行业今后成长发展的趋势；黄洋洋[3]建立了一个模型来讨论制造业企业的违约风险；罗泰晔[4]运用逻辑斯谛模型的系数研究网络舆情传播规律；唐旭，朱灵伟等[5]利用逻辑斯谛模型分析新一轮退耕还林工程四川省农户参与意愿的影响因素，为政府优化实施政策提供参考依据。

尽管逻辑斯谛模型分析解决实际问题有非常重要的价值，但是对于比较复杂的逻辑斯谛模型，微分方程的精确解是得不到的，所以探索它的近似解计算方法就显得非常必要而且意义非同一般。

Adams插值方法是有限差分多步法中的一种常用方法，是通过数值积分或数值微分建立差分模型，运用反复迭代和现代计算技术，通过离散点来逼近理论解，从而得到精确度较高的近似解的一种方法。

该方法容易掌握且应用广泛，王磊[6]研究了方法的理论问题，包括稳定性条件、误差以及计算效率等；廖翠萃[7]提出了更稳定优异的方法以求解刚性常微分方程；程小昊[8]研究了求解振荡问题的保结构线性多步法，解决了周期性振荡微分方程的求解问题；许波峰，唐植懿等[9]将多步法应用于离散涡线对流控制方程，准确地预估了风力机稳态气动性能。

ADAMS多体系统建模与动力学方程迭代求解

Modeling of Multi- System and Iteration Solving of Dynamic Equations Basing on ADAMS
MA Ji-sheng (No.1 Department, Ordnance Engineering College, Shijiazhuang 050003, China)
积分步长，再行预估校正，这时外界积分程序也必须倒退到上一积分时刻，重新进行积分。 (3) 外界程序的积分步长必须与 ADAMS 的积分步长完全一致。为了避免在外界程序中计算 ADAMS 的积分步长的麻烦，可在 ADAMS 的求解器设置中，将初始积分步长、最大积分步长和外界程序的积分步长设成一样， ADAMS 的最小积分步长可稍小，这样 ADAMS 基本上是按定步长积分。要注意的是，取相同值的三个积分步长要足够小，以保证外界程序和 ADAMS 的求解器均能收敛。 3.4 计算结果
1 a/m.s－ 2 0.9 0.8 0.7
j= 1
&j = − ω 2 && τ ∑ d ij& i τ i − c ix
N
力 Fx 由一个用户子程序计算，在该子程序中， & & 通过调用 ADAMS 的状态变量获取子程序给出， x & j ( j = 1L N ) ，代入 (3) 式得到 F x 。 Fx 代入 (4) 式得到 & τ
收稿日期： 2003 － 10 － 22 ；修回日期： 2004 － 01 － 15 作者简介：马吉胜（ 1967－），男，河北人，博士，军械工程学院教授，1994 年毕业于南京理工大学，从事多体系统动力学、振动信号处理、模式识别与神经网络研究。

adams多体动力学拓扑

Adams多体动力学拓扑概述Adams多体动力学是一种用于模拟多体系统运动的数学方法。

它基于牛顿力学和运动学原理，通过建立动力学方程来描述系统的运动行为。

动力学拓扑则是在Adams 多体动力学中的一个重要概念，用于描述系统中各个物体之间的联系和相互作用。

动力学拓扑的基本概念动力学拓扑是通过定义物体之间的连接关系和相互作用来描述多体系统的结构和行为的方法。

在Adams多体动力学中，动力学拓扑可以通过节点和边的方式表示。

节点表示系统中的物体，边表示物体之间的连接关系或相互作用。

动力学拓扑可以用于描述不同类型的多体系统，如机械系统、电子系统、生物系统等。

在机械系统中，节点可以表示机械零件，边可以表示零件之间的连接关系或相互作用。

在电子系统中，节点可以表示电子元件，边可以表示元件之间的电路连接或电磁相互作用。

在生物系统中，节点可以表示生物分子或细胞，边可以表示分子之间的化学反应或细胞之间的信号传递。

动力学拓扑的建立需要考虑系统的结构和物体之间的相互作用。

通过分析系统的结构和相互作用，可以确定系统中的节点和边的数量，以及它们之间的连接方式和特性。

这些信息可以用于建立动力学方程，并通过求解方程来模拟系统的运动行为。

Adams多体动力学模拟Adams多体动力学模拟是一种基于动力学方程的数值计算方法，用于模拟多体系统的运动行为。

在Adams多体动力学模拟中，系统的运动行为可以通过求解动力学方程来得到。

动力学方程描述了系统中各个物体的运动行为和相互作用。

它可以通过牛顿第二定律和运动学原理来建立。

牛顿第二定律指出，物体的加速度与作用在物体上的合外力成正比，与物体的质量成反比。

运动学原理则描述了物体的位置、速度和加速度之间的关系。

Adams多体动力学模拟可以通过数值计算方法来求解动力学方程。

常用的数值计算方法包括欧拉法、龙格-库塔法等。

这些方法可以将动力学方程离散化为一系列时间步长上的更新方程，通过迭代求解来模拟系统的运动行为。

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控制模型以及结果后处理四个步骤。

求解声波方程的辛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.。

微分方程数值解法实验报告-Runge-Kutta格式、三阶Admas

《微分方程数值解法》课程实验报告1 实验内容要求：用经典三级三阶R-K 格式求解微分方程初值问题，并给出误差分析。

2算法描述先用C 的方法写出一个算法动态库,里面封装龙格库塔算法函数采用迭代原理，用递归实现,生成并模块化导出dll,lib 文件,两个文件中均包含了该函数偏移地址,在在源文件中隐式链接该库里的龙格库塔函数,从而得出结果.3 实验数据与实验结果 ⎩⎨⎧=-+-==1|1'0x y x y y1.0)1,0(=∈h x4 程序代码清单：Win32动态库Algorithm.dll 代码: #include <stdio.h> #include <math.h>#define e 2.718281828459045double x[11]={0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0}; double y[12]={1.0,0}; int x1=0,x2=1; int count=1;void RungeKutta(double h){ double K1,K2,K3; double ExactSolution; double Error; if(count==12){ return; } K1=-y[count-1]+x[count-1]+1; K2=-(y[count-1]+h*K1/3)+(x[count-1]+h*K1/3)+1; K3=-(y[count-1]+2*h*K2/3)+(x[count-1]+2*h/3)+1; y[count]=y[count-1]+h*(K1+3*K3)/4; ExactSolution=x[count-1]+pow(e,-x[count-1]); Error=y[count-1]-ExactSolution; printf("%lf %lf %lf %lf\n",x[count-1],ExactSolution,y[count-1],Error);count++; RungeKutta(h); }模块化导出文件Algorithm.def代码:LIBRARY AlgorithmEXPORTSRungeKutta @1入口函数实现功能代码:#include <stdio.h>#pragma comment(lib,"../lib/Algorithm.lib")int main(){double h;printf("用三级三阶龙贝格库塔方法解微分方程y'=x-y+1,x属于区间(0,1)\n");printf("请输入系数h的值:\n");scanf("%lf",&h);printf("-----------------------------------------------------\nx 精确解 R-k解yn 误差:\n");RungeKutta(h);return 0;}5:运行结果:1 实验内容要求：用三阶Admas 预报修正格式求解差分方程初边值问题。

一种使用隐式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

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 结论 (17)6 参考文献 (19)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)。