下载文档原格式
地震波波动方程数值模拟方法地震波波动方程数值模拟方法主要包括克希霍夫积分法、傅里叶变换法、有限元法和有限差分法等。
克希霍夫积分法引入射线追踪过程，本质上是波动方程积分解的一个数值计算，在某种程度上相当于绕射叠加。
该方法计算速度较快，但由于射线追踪中存在着诸如焦散、多重路径等问题，故其一般只能适合于较简单的模型，难以模拟复杂地层的波场信息。
傅里叶变换法是利用空间的全部信息对波场函数进行三角函数插值，能更加精确地模拟地震波的传播规律，同时，利用快速傅里叶变换(FFT)进行计算，还可以提高运算效率，其主要优点是精度高，占用内存小，但缺点是计算速度较慢，对模型的适用性差，尤其是不适应于速度横向变化剧烈的模型.波动方程有限元法的做法是：将变分法用于单元分析，得到单元矩阵，然后将单元矩阵总体求和得到总体矩阵，最后求解总体矩阵得到波动方程的数值解；其主要优点是理论上可适宜于任意地质体形态的模型，保证复杂地层形态模拟的逼真性，达到很高的计算精度，但有限元法的主要问题是占用内存和运算量均较大，不适用于大规模模拟，因此该方法在地震波勘探中尚未得到广泛地应用。
相对于上述几种方法，有限差分法是一种更为快速有效的方法。
虽然其精度比不上有限元法，但因其具有计算速度快，占用内存较小的优点，在地震学界受到广泛的重视与应用。
声波方程的有限差分法数值模拟对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：)()(2222222t S zu x u v t u +∂∂+∂∂=∂∂ （4-1） (,)v x z 是介质在点（x , z ）处的纵波速度，u 为描述速度位或者压力的波场，)(t s 为震源函数。
为求式（4-1）的数值解，必须将此式离散化，即用有限差分来逼近导数，用差商代替微商。
为此，先把空间模型网格化（如图4-1所示）。
设x 、z 方向的网格间隔长度为h ∆，t ∆为时间采样步长，则有：z∆,i j1,i j +2,i j+1,i j-h i x ∆= （i 为正整数）h j z ∆= (j 为正整数)t n t =∆ (n 为正整数)k j i u , 表示在(i,j)点，k 时刻的波场值。
波动问题中的三维时域粘弹性人工边界一、本文概述在波动问题研究中，粘弹性人工边界作为一种重要的数值模拟方法，被广泛应用于地震工程、岩土工程、结构动力学等领域。
本文将重点探讨三维时域粘弹性人工边界在波动问题中的应用。
我们将对粘弹性人工边界的基本理论进行介绍，包括其发展历程、基本原理以及在波动问题中的应用背景。
随后，我们将详细介绍三维时域粘弹性人工边界的建模方法、数值实现过程以及关键参数的选取。
我们还将分析三维时域粘弹性人工边界在波动问题中的优势和局限性，以及在实际应用中可能遇到的问题和解决方法。
我们将通过具体案例来展示三维时域粘弹性人工边界在波动问题中的实际应用效果，并总结其在实际工程中的应用前景。
本文旨在为从事波动问题研究的学者和工程师提供一种有效的数值模拟方法，以更好地理解和解决实际工程中的波动问题。
通过本文的介绍和分析，读者可以深入了解三维时域粘弹性人工边界的基本原理、数值实现方法以及实际应用效果，为相关研究提供有益的参考和借鉴。
二、波动问题基本理论波动问题，作为物理学和工程学中的核心领域，主要研究波在介质中的传播规律。
波的传播受介质特性、波的初始条件和边界条件等多种因素影响。
波动问题涉及弹性力学、动力学、波动方程等多个学科分支，其基本理论为理解和分析复杂波动现象提供了基础。
在波动问题中，波动方程是描述波传播行为的关键。
一维情况下，波动方程可以表示为 (\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2})，其中 (u) 是波的位移，(t) 是时间，(x) 是空间坐标，(c) 是波速。
这一方程描述了波在均匀、无阻尼介质中的传播行为。
对于三维情况，波动方程需要考虑三个空间维度，形式更为复杂。
同时，波动方程还需要结合具体的介质特性，如弹性模量、密度等，来求解特定问题的波动行为。
在波动问题中，边界条件对于波的传播具有重要影响。
FLAC3D动力分析中的人工透射边界和地震波施加方法从动力学的角度上看，动力响应是确定惯性（质量效应）和阻尼起着重要作用时质点或质点系动力学特性和响应的技术，它包括自振、冲击、谐振动、随机振动等分支。
动力学最早应用于结构抗震设计，自上世纪50年代逐步借鉴到岩土抗震设计中。
动力发展历程可总结为静力理论，反应谱理论和时程分析理论三个阶段。
我们知道，地震的三要素为振幅、频谱和持时。
静力理论只考虑了地震引起的最大振幅，属于拟静力法；反应谱理论考虑了振幅和频谱，但在设计中仍然把地震惯性力视为静力，只能算准动力法；时程分析理论考虑了振幅、频谱和持时，是严格意义上的动力分析法。
通常时程动力分析选用的地震波来自：（1）根据设计反应谱人工合成的场地波；（2）场地附近地震台记录的实测地震波。
由于实测地震波中掺杂了许多噪声和干扰信号，因此在使用前必须滤波去噪、频谱分析、积分变换和基线修正。
滤波去噪是为了消除噪声和高频波，频谱分析是为了检测地震波持时内所含的频率分量和振幅，积分变换可以转换地震加速度波为速度波或位移波，基线修正则是为了消除非平稳地震波中的弹性位移零线漂移、基线偏移等现象，大崎顺彦在其著作《地震动的谱分析入门》中做了详细而生动的说明，并附出了地震波处理的Fortran源程序。
鉴于FLAC3D软件是岩土领域广泛应用的时程动力分析软件，这里以著名的埃尔森特罗波（El Centro）为输入激励，研究基于FLAC3D软件的地震波处理和计算方法。
网站“http://www. /data.htm”提供了31秒的El Centro加速度波数据。
有兴趣者可按《地震动的谱分析入门》的方法选取了前8秒的地震加速度波（共401个记录），然后补零配成了512个记录的加速度波以采用快速傅里叶变换法，首先采用FLAC3D Fish函数库的filter函数进行滤波去噪，然后采用fft函数进行快速傅里叶变换，得到傅里叶加速度谱和功率谱，接着采用integrate函数积分两次求得速度波和位移波，并计算地震位移零线漂移值。
人工边界方法
韩厚德
【期刊名称】《数学建模及其应用》
【年(卷),期】2012(0)3
【摘要】人工边界方法是数值求解无界区域上偏微分方程的一类有效计算方法。
本文以二维Poisson方程的外问题和声波方程为例,介绍人工边界方法的基本思想和核心技术。
【总页数】10页(P1-10)
【作者】韩厚德
【作者单位】清华大学数学科学系
【正文语种】中文
【中图分类】O241.82
【相关文献】
1.动力分析中人工地基边界处理方法对比研究 [J], 高颖;郭庆林;杨永辉;卞艳山;张聪
2.黏弹性人工边界地震波动输入方法综述 [J], 王晓东
3.重力坝动力分析黏弹性人工边界及其地震动输入处理方法 [J], 谯雯
4.一种基于人工提取缺陷块的边界搜索方法 [J], 马天航;胡家铖;郑莉;刚蓓;刘思娇
5.基于人工边界方法的西藏旁多土石坝非线性动力分析 [J], 吴悦;郭永刚;胡锦因版权原因，仅展示原文概要，查看原文内容请购买。
动力分析中人工地基边界处理方法对比研究高颖;郭庆林;杨永辉;卞艳山;张聪【摘要】针对地基模型边界的近似处理问题,运用有限元方法分析了固定边界、等效粘弹性边界和无限元边界近似处理后的地基动力响应规律,借助已有的现场实测数据分析了各种边界处理方法的可操作性和精确性.结果表明:等效粘弹性边界和无限元边界均能反映边界辐射振动特性,从准确性角度来看,无限元方法具有更多优势.%For the foundation boundary approximation in the vibration analysis,dynamic responses of foundation treated by fixed boundary,equivalent viscoelastic boundary and infinite element boundary were obtained by using the fmite element method.Operability and accuracy of various boundary methods were analyzed according to the measured data.The results show that both the equivalent viscoelastic boundary and the infinite element boundary can better reflect boundary vibration characteristics,but infinite element method has a better advantage in terms of accuracy.【期刊名称】《河北工程大学学报（自然科学版）》【年(卷),期】2018(035)001【总页数】4页(P1-4)【关键词】土体动力响应;人工边界;无限元;等效粘弹性边界;有限元【作者】高颖;郭庆林;杨永辉;卞艳山;张聪【作者单位】河北工程大学土木工程学院,河北邯郸056038;河北工程大学土木工程学院,河北邯郸056038;河北工程大学土木工程学院,河北邯郸056038;河北工程大学土木工程学院,河北邯郸056038;河北工程大学土木工程学院,河北邯郸056038【正文语种】中文【中图分类】TU47运用有限元法研究结构-地基相互作用问题时，常常由于有限元模型难以反映无限域介质的辐射影响效应而使研究受挫。
地震波方程人工边界的两种处理方法
崔兴福;张关泉
【期刊名称】《石油物探》
【年(卷),期】2003(042)004
【摘要】在分析以往人工边界处理优缺点的基础上,提出了利用波动方程的频散关系式,得到可以吸收任何方向入射波的自适应校正吸收边界条件和旋转校正吸收边界条件.同时,采用波阵面法和能流密度法判别入射波方向,克服了Pade近似吸收边界只对正入射波具有较好吸收性,而对非正入射的波吸收性不好的缺点.数值试验结果表明了本方法的有效性.
【总页数】4页(P452-455)
【作者】崔兴福;张关泉
【作者单位】中国石油辽河油田分公司勘探开发研究院计算所,辽宁,盘锦,124010;中国科学院数学与系统科学研究院计算数学与科学工程计算研究所,北京,100080;中国科学院数学与系统科学研究院计算数学与科学工程计算研究所,北京,100080【正文语种】中文
【中图分类】P631.4
【相关文献】
1.适用于VTI介质地震波方程的PML吸收边界 [J], 杨佳佳;郭鹏
2.基于边界元法的地震波方程时间域数值模拟 [J], 刘志伟;王忠仁;裴海侠
3.准P波方程紧致交错网格井间地震波场模拟及边界条件 [J], 孟凡顺;张亮;李景岩;
李洋森
4.柱坐标声波方程正演两种PML边界加载方式的对比 [J], 郭锐;王尚旭
5.重力坝动力分析黏弹性人工边界及其地震动输入处理方法 [J], 谯雯
因版权原因，仅展示原文概要，查看原文内容请购买。
第42卷第4期2003年12月石 油 物 探GEOPHYSICAL PROSPECTIN G FOR PETROL EUMVol.42,No.4Dec.,2003文章编号:100021441(2003)0420452204地震波方程人工边界的两种处理方法崔兴福1,2,张关泉2(1.中国石油辽河油田分公司勘探开发研究院计算所,辽宁盘锦124010;2.中国科学院数学与系统科学研究院计算数学与科学工程计算研究所,北京100080)摘要:在分析以往人工边界处理优缺点的基础上,提出了利用波动方程的频散关系式,得到可以吸收任何方向入射波的自适应校正吸收边界条件和旋转校正吸收边界条件。
同时,采用波阵面法和能流密度法判别入射波方向,克服了Pade近似吸收边界只对正入射波具有较好吸收性,而对非正入射的波吸收性不好的缺点。
数值试验结果表明了本方法的有效性。
关键词:自适应校正;旋转校正;波阵面;能流密度中图分类号:P63114 文献标识码:ATwo processing methods for artif icial boundary of seismic w ave equationCui Xingfu1,2,Zhang Guanquan2(puting Center of Exploration&Development Research Institute of CNPC Liaohe Oilfield Com pany,Panjin124010,China;2.Institute of Com putational Mathematics and Scientific/Engineering Computing,Academy of Math2 ematics and System Sciences,Chinese Academy of Sciences,Beijing100080,China)Abstract:In this paper,two absorbing boundary conditions,adaptive correction condition and rotation correction condi2 tion,were derived to absorb incident waves coming from any directions by using dispersion relation,based on an analy2 sis of the advantages and disadvantages of existing boundary conditions.The determination of incident wave direction by wave front and energy flux density was also described.Numerical ex periments were performed and their results were presented,which indicated the efficiency of these methods.K ey w ords:adaptive correction;rotation correction;wave front;energy flux density 实际人工模拟地震勘探是在半无界空间中进行的,而在计算机上进行数值模拟地震波在介质中的传播,只能在有限的模型空间中实现。
地震波到达人工边界将产生虚假的反射波,干扰了实际地震波传播的机理,使仿真剖面变得模糊不清,不利于地下地层构造信息的解释。
自1969年Lysmer和Kuhlemeyer[1]首先提出人工边界处理问题,发展至今,已形成多种人工边界的处理方法,如阻尼边界[1～3]、吸收边界[4～6]、透射边界[7]、移动边界等,用Pade近似法得到的吸收边界条件[4,5]是目前边界处理效果比较好的一种。
我们正是在分析这种边界处理优缺点的基础上,从不同的侧面提出了这种边界的2种校正方法,即自适应校正法和旋转校正法。
1 声波方程自适应校正吸收边界条件1.1 二阶精度吸收边界的推导由二维声波方程u t t=c2(u x x+u zz)(1)进行三维Fourier变换得到频散关系式ω2=c2(k2x+k2z)=k2c2(2)式中,k2=k2x+k2z,k x=k cosθ,k z=k sinθ,θ是波前k和k z的夹角。
引进中间参量a(θ)=1-cosθ1-cosθ2(3)由图1表示的方向波数与波前关系可以推得cosθ2=(α-1)k+k xαk(4)利用k x=k cosθ,cosθ=2cos2θ2-1和2α2-4α+4=11+cosθ及方程(4)得到在频率波数域右行波的边界条件ωkx-1cω2+11+cosθk2z=0(5)收稿日期:20021113;改回日期:20030323。
作者简介:崔兴福(1965—),男,高级工程师,博士在读,现从事地震资料处理方法研究工作。
基金项目:本文得到国家973重点基础研究项目(G199903280)资助。
由Fourier 逆变换得到时间空间域中右行波的边界条件u x t +1c u t t -c1+cos θu zz =0(6)同理可以得到左行波、上行波、下行波的边界条件。
方程(6)中存在一个如何判别入射波方向的问题,下面给出判定波入射方向的第一种方法。
1.2 波入射方向的判定———波阵面法根据波阵面理论,在任何时刻,同一波阵面(波前)上各质点的运动状态相同,并且波阵面的外法线方向一定是波的传播方向。
具体地,为了确定在第(n +1)Δt 时刻边界点M 的位移u (n +1)(M )的值,用第n Δt 和(n -1)Δt 时刻所有各点以及第(n +1)Δt 时刻内部节点的位移值表示,搜索满足如下条件的点m 0。
min m ∈Ω(M)|u n (M )-u n (m )|=|u n (M )-u n (m 0)|(7)可以近似地认为m 0点和边界点M 在同一波阵面上,则二者连线垂直方向可以认为是入射于点M 的波传播方向。
实际应用表明,用这种方法确定入射方向,对平面波比较精确,而对于球面波相对差一点。
2 波动方程旋转校正吸收边界条件文献[8]中对各种近似的边界进行了详细的分析,用Pad e 近似得到的边界条件,具有良好的局部逼近性优点和不好的整体逼近性缺点。
鉴于此,我们通过旋转Pad e 近似边界得到旋转校正边界条件。
2.1 旋转校正吸收边界条件的建立在二维空间中旋转具有如下的形式x ′=x cos θ+y sin θy ′=-x sin θ+y cos θ(8)x =x ′cos θ-y ′sin θy =x ′sin θ+y ′cos θ(9)二阶Pad e 近似的右边界条件为u x t +1c u t t -c2u zz=0(10)方程(10)是方程(6)在θ=0的特殊情况,对正入射的地震波吸收较好,而对非正入射的地震波吸收不好。
将方程(10)变换到频率波数域中,有ωc k x -ωc 2+12k 2z=0(11)旋转方程(11),得ωc(k x cos θ+k z sin θ)-ωc2+12(k z cos θ-k x sin θ)2=0(12)对应到时间空间域得到旋转后的右边界条件-u xt -tan θu tz -1c cos θu t t +c2tan θsin θu x x -c sin θu xz +12c cos θu zz =0(13) 同理可以得到左边界、底边界条件。
用这种方法可以对弹性波的近似边界进行旋转得到弹性波旋转校正边界。
方程(13)中也存在一个如何判别入射波方向的问题,下面给出判定波入射方向的另一种方法。
2.2 入射波方向的确定———能流密度法能流密度矢量的方向和波传播方向相同,在数量上等于单位时间内流经垂直于波传播方向的单位截面的流量。
弹性波的能流密度为I =I 1i +I 2j +I 3k(14)式中,I 1=σ11u 1+σ12u 2+σ13u 3I 2=σ21u 1+σ22u 2+σ23u 3I 3=σ31u 1+σ32u 2+σ33u 3u 1,u 2和u 3为位移向量U 的各个分量。
在二维(i ・k )情况下得到的弹性波能流密度具体表达式如下I =(2μ+λ)5u 15x 1 u 1+λ5u 35x 3 u 3+μ5u 15x 3+ 5u 35x 1 u 3i +μ5u 15x 3+5u 35x 1 u 1+(2μ+λ)5u 35x 3 u 3+λ5u 15x 1 u 3k (15)式中,λ和μ为拉梅系数。
声波的能流密度为ω=P U (16)式中,P 为声压。
只要确定出波在某一点能流密度矢量各个分量值的大小,就确定了传播到该点波的传播方向。
2.3 数值计算实现步骤在地震波正演的数值模拟中,不论是采用有限差分法还是有限元法,人工边界处理是关系到正演模拟效果好坏的关键。
在实际处理中,首先是解决如何判别入射波传播方向的问题,即方程(6)和(13)中的入射方向角θ,可以采用方程(7)或方程(15)来求取,方法(7)易于实现,但不十分准确,方法(15)是一种判别入射波方向的较好方法,但需要给出地质模型的拉梅系数λ和μ;采用自适应校正吸收边界方程(6)或旋转校正吸收边界方程(13)进・354・第4期崔兴福等1地震波方程人工边界的两种处理方法行边界处理,前者较为简明,后者烦琐,但效果较好。
下面分别给出方程(6)和方程(13)的差分离散格式。
将波场函数进行差分离散u (x ,z ,t )=u mi ,j ,i =1,2,…,I ;j =1,2,…,J ;m =1,2,…,M ,i ,j 和k 分别表示x ,z 和t 方向的离散节点号,则方程(6)和方程(13)的离散差分格式为u m +1i +1,j -u m +1i ,j -u mi +1,j +u mi ,jΔt Δx+ 1cu m +1i ,j -2u m i ,j +u m -1i ,jΔt 2-c1+cos θu mi ,j +1-2u mi ,j +u mi ,j -1Δz 2=0(17)-u m +1i +1,j -u m +1i ,j -u mi +1,j +u mi ,jΔt Δx -tan θu m +1i ,j +1-u m +1i ,j -u mi ,j +1+u mi ,jΔt Δz-1c cos θu m +1i ,j -2u m i ,j +u m -1i ,j Δt 2+c2tan θsin θu m i +1,j -2u m i ,j +u m i -1,jΔx 2-c sin θu mi +1,j +1-u mi +1,j -u mi ,j +1+u mi ,jΔx Δz +12c cos θu m i ,j +1-2u m i ,j +u m i ,j -1Δz 2=0(18)3 数值计算比较图1是方向波数与波前关系示意图。
