基于三阶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动力学仿真算法及参数设置分析
( 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 版求解器来完成。
- 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)。