高精度交错网格有限差分算法的改进.docxVIP

高精度交错网格有限差分算法的改进 0 u3000交错网格有限差分算法 随着我国石油天然气勘探和开发的不断发展，研究人员面临的勘探对象和开发条件变得越来越复杂和困难。寻找复杂构造油气藏、岩性油气藏和裂缝油气藏等“剩余油”的任务艰巨［１－２］。为了解决这些复杂油气藏的勘探、开发问题,需要对地下复杂介质的地震响应进行高精度数值模拟研究。对地震弹性波方程进行数值模拟的方法主要包括有限差分法、有限元法和伪谱法［３］,而有限差分法由于计算速度快、精度高等特点使用最为广泛［４－５］。 自从Alterman和Karal首次将有限差分法应用到各向同性介质弹性波的模拟中后,由于其自身的优点,很快被用于各种地震勘探学的数值问题上,且在应用中不断发展［６］,先后出现了变网格有限差分［７－８］、非连续网格有限差分［９］、不规则网格有限差分［１０］、交错网格有限差分［１１－１３］、旋转交错网格有限差分［１４－１５］、可变时间步长有限差分［１６］、自适应可变空间步长网格有限差分以及隐式有限差分［１７－１８］。其中,交错网格由Madariage［１１］最早提出,Virieux［１２］首先将其使用到一阶速度-应力方程中。 为了避免对弹性常数进行空间微分,在弹性波正演模拟时,采用一阶速度-应力弹性波方程。这样可以在不同的时间层上使用不同的网格,分别进行应力和速度的计算和传播。Virieux发展的交错网格精度为O(Δt２+Δx２)(时间2阶、空间2阶),与常规网格相比,在没有增加计算量和存储空间的情况下,局部精度提高了4 倍,收敛速度也有所提高。Levander［１９］又将交错网格有限差分的精度提高到O(Δt２+Δx４)。随后董良国等［２０－２１］发展了更高精度的交错网格有限方法,精度达到了O(Δt４+Δx２u3000Ｎ)。为了使用较少的时间层,不增加计算存储空间,董良国等［２０－２１］将速度(应力)对时间的奇数次高阶导数转化为应力(速度)对空间的导数,把交错网格和高阶差分法成功地结合在一起。然而在高阶交错网格有限差分算法中,当变量对空间求三次导数时,采用的是低阶的差分格式和差分系数［２２－２５］,没有充分考虑“空间2 N阶精度”这一事实,这样做虽然有助于提高有限差分的计算效率,但要想获得更高精度的正演结果就有必要推导出准确的时间4 阶、空间2 N阶精度的交错网格有限差分格式及相应的差分系数。 差分系数是决定交错网格有限差分算法精度的关键,Taylor公式展开和最优化方法是求取差分系数的主要方式［１，１７，２６，２８］。本文以Taylor公式展开为基础,推导了可导函数任意次导数的任意偶数阶精度的差分近似式以及差分系数,完善O(Δt４+Δx２u3000Ｎ)一阶速度-应力方程组差分格式。 1 阶速度-应力弹性波方程 在二维各向同性介质xoz平面内,假定体力为0,一阶速度-应力弹性波方程为 其中:vｘ,vｚ为速度分量;τｘｘ,τｚｚ,τｘｚ为应力分量;ρ为密度;c１１,c１３,c３３,c４４为介质的弹性常数。 1.1 m阶差分近似 在使用交错网格有限差分算法时,速度和应力分别在t+Δt/2和t时刻进行计算。将和用Taylor公式展开,可得到2 M阶精度的时间差分近似,即 式中,Δt为时间步长。令M=2,式(2)就是常规的时间4阶精度差分近似。 为了减少计算内存,利用速度和应力的耦合关系,得到方程组式(1)的时间4阶精度差分近似,以式(1-1)为例: 1.2 空间一次导数 在常规交错网格算法中,速度(应力)分量的空间一次导数是由相错半网格点的应力(速度)分量计算的。设函数f(x)具有2 N-1阶导数,令x=x０±[(2n-1)/2]Δx,则由f(x)在x处的2 N -1 阶Taylor展开可以得到式(3)中空间一次导数的差分近似式: 式中:Δx为空间步长;Cｎ（Ｎ）为待定差分系数。通过求解下面的方程组就可以求出空间一次导数的差分近似式的差分系数: 1.3 应力的离散值 取x=iΔx,z=jΔz,t=kΔt,i、j和k分别表示空间和时间网格点,U、V分别代表速度分量vｘ、vｚ的离散值,R、T、H分别代表应力τｘｘ、τｚｚ、τｘｚ的离散值。则方程(3)的精度为O(Δt４+Δx２u3000Ｎ)的常规差分格式如下(Δx=Δz,下同): 可以明显地看出,式(6)中P１１、P１２、P１３、P１４及P１５的差分格式不是空间2 N阶精度的,而是用低阶差分格式代替。 2 速度-应力波动方程的差分格式及差分系数 为了得到交错网格有限差分的时间4阶、空间2 N阶精度准确的差分格式及差分系数,首先需要推导出函数对空间二次、三次导数的差分格式及差分系数,然后将其代入到一阶速度-应力波动方程组中,从而推导出该方程组精度为时间4阶、空间2 N阶的准确差分格式。 2.1 改进差分近似式 函数对空间任意