不稳定地温场的有限差分模拟
寒区隧道温度场综述
寒区隧道温度场综述摘要：近年来，随着隧道工程施工技术发展的越来越成熟，公路与铁路的建设在路线选择上更加多元化，由原来的避山绕路到现在的凿山取直，隧道工程发挥着重要的作用。
但是部分地区冬天寒冷的气温对隧道衬砌结构的稳定性造成极大的威胁，同时为行车安全留下隐患。
因此对季冻区隧道围岩温度场分布情况的掌握，就是我们需要亟待解决重要问题之一。
关键词：寒区隧道；温度场；冻害研究；综述引言：随着对气温寒冷地区隧道的研究逐渐深入，地下工程研究人员逐渐发现隧道工程中产生失稳现象以及冻害现象的原因是因为在工程实际中水分、温度和应力这三种状态相互耦合，相互影响。
但由于我国在季冻区隧道温度场研究领域起步较晚，对这一问题的认识不足，在实际工程中不能切实有效的施以更好的措施防止冻害现象的发生；再者，伴随着一带一路政策的实施，国内寒区隧道建设日益增多，加快对其的研究工作有助于更好地提高寒区工程建设的效率以及安全度。
对于季冻区隧道温度场的研究工作近几年主要围绕温度场的分布规律、热物性等对温度场影响研究及借助软件进行数值模拟的研究上。
对寒区隧道温度场的研究有助于认识隧道冻害机理，更有针对性的采取相应的防治冻害措施，确定防排水设施类型，为此国、内外隧道工作者进行了大量研究工作。
1、国外研究美国一实验研究所[1]从上世纪60年代开始进行一段较长时间的监测专门观察一座多年冻土隧道，通过此次观测在寒区隧道温度场分布规律方面的研究取得了一系列科研成果，也对进一步开展冻土研究有着重要价值。
1973年，国外学者Bonacina[2]和Cominit[3]分别用数值求解的方法解决温度场问题，是地下工程领域首次对寒区隧道温度场多场耦合数值模拟及非线性分析进行研究[4][5]。
2、国内研究近些年，在我国西部与东北部等寒冷地区修建公路、铁路隧道越来越多，我国许多隧道工作者对各个不同的隧道进行了围岩温度的现场监测，通过长期观测隧道围岩温度场，其变化规律已经被掌握。
有限差分法-导热模拟
有限差分法-导热模拟有限差分法（Finite Differential Method ）是基于差分原理的一种数值计算法。
其基本思想：将场域离散为许多小网格，应用差分原理，将求解连续函数ϕ的泊松方程的问题转换为求解网格节点上ϕ的差分方程组的问题。
一、利用有限差分法离散三维傅立叶热传导微分方程：T z T y T xT t T 2222222∇=⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂=∂∂αα （1-1）解：将三维温度场域划分为足够小的正方体网格，网格之间距离为h ，图一显示为节0(i,j,k)及其周围的节点1(i-1,j,k)、2(i+1,j,k)、3(i,j-1,k)、4(i,j+1,k)、5(i,j,k-1)、6(i,j,k+1)。
节点上的电位分别用6543210T T T T T T T ，，，，，，表示由有限差分法得：2220122)1()(2)1(2)(0hk j i T k j i T k j i T h T T T x T x x ，，，，，，++--=+-≈∂∂= （1-2） 同理：2240322)1()(2)1(2)(0hk j i T k j i T k j i T h T T T y T y y ，，，，，，++--=+-≈∂∂=（1-3） 2260522)1()(2)1(2)(0hk j i T k j i T k j i T h T T T z T z z ++--=+-≈∂∂=，，，，，，（1-4） 将时间t 划分为足够小的时间段，时间节点之间的距离为g ，则采用有限差分法的后向差分法得：g T T dt dT n n 1--≈ （1-5） Z YX 1(i-1,j,k)0(i,j,k)2(i+1,j,k)3(i,j-1,k)4(i,j+1,k) 5(i,j,k-1) 6(i,j,k)图1 三维节点图将式（1-2）、（1-3）、（1-4）、（1-5）代入式（1-1）得：()2121)()1()1()1()1()1()1()(61)](6)1()1()1()1()1()1([)()(h gr k j i T k j i rT k j i rT k j i rT k j i rT k j i rT k j i rT k j i T r k j i T k j i T k j i T k j i T k j i T k j i T k j i T hg k j i T k j i T n n n n n n n n n n n n n n n n n αα==+---+---+---+⇒-++-+++-+++-=---其中：，，，，，，，，，，，，，，，，传导差分公式上式整理可推出三维热，，，，，，，，，，，，，，，，，，求解完毕。
温度场分布仿真计算方法
温度场分布仿真计算方法温度场分布仿真计算方法温度场分布仿真计算方法是一种通过数值模拟和计算机仿真来预测和分析温度分布的方法。
它在工程设计、热力学研究和环境保护等领域中得到广泛应用。
本文将介绍温度场分布仿真计算方法的基本原理和常用技术。
温度场分布仿真计算方法的基本原理是建立一套数学模型来描述温度场的变化规律，并通过计算机程序对模型进行求解和模拟。
根据具体问题的需求和实际情况，可以选择不同的数学模型和计算方法。
常见的数学模型包括传热方程、能量守恒方程和流体动力学方程等。
计算方法主要包括有限差分法、有限元法和边界元法等。
有限差分法是最常用的一种计算方法。
它将温度场划分为若干个网格点，并通过计算相邻网格点之间的温度差来近似描述温度场的变化。
有限差分法的优点是计算简单，适用于各种尺度和几何形状的问题。
但是，它需要较密集的网格划分，以获得较精确的结果。
有限元法是一种更精确的计算方法。
它将温度场划分为若干个有限元素，通过求解每个元素上的温度分布来近似描述整个温度场。
有限元法的优点是可以灵活地处理复杂的几何形状和边界条件。
但是，它需要对模型进行离散化处理，计算量较大。
边界元法是一种特殊的计算方法。
它通过求解温度场的边界值来推导出整个温度场的分布。
边界元法的优点是计算量较小，适用于二维和三维问题。
但是，它对边界条件的要求较高，需要较精确的输入数据。
除了上述常用的计算方法外，还有一些其他的技术和方法可以用于温度场分布仿真计算，如Monte Carlo方法、遗传算法和人工神经网络等。
这些方法可以根据具体问题的需求进行选择和组合，以获得更准确和可靠的结果。
综上所述，温度场分布仿真计算方法是一种重要的工程分析工具。
它通过数值模拟和计算机仿真来预测和分析温度场的分布规律，为工程设计和科学研究提供了有力的支持。
随着计算机技术的不断发展和进步，温度场分布仿真计算方法将更加精确和高效，为解决实际问题提供更好的解决方案。
研究报告有限差分格式稳定性的其他方法-报告
研究有限差分格式稳定性的其他方法摘要偏微分方程的求解一直是大家比较关心的一个问题，而有限差分格式则是求解偏微分方程时常用并且有效的一个方法。
因此，研究有限差分格式的性质就显得尤为重要。
在课上我们已经跟着老师学习了运用Fourier方法研究有限差分格式的稳定性，但是在很多研究有限差分格式稳定性的问题中仅仅会用Fourier方法是不够的，所以在本篇论文中，将会介绍其他三种常用的研究有限差分格式稳定性的方法，分别是：Hirt启示型方法、直接方法（或称矩阵方法）和能量不等式方法。
关键字：偏微分方程；有限差分格式；稳定性AbstractThe solution of partial differential equations has been more concerned with a problem, and the finite difference scheme is a mon and effective method for solving partial differential equations. Therefore, it is very important to study the character of the finite difference scheme. We have followed the teacher to learn the use of Fourier method of finite difference scheme stability, but in a lot of research on the stability of finite difference scheme is only used Fourier method is not enough, so in this paper, will introduce the other three kinds of monly used in the study of finite difference scheme stability method, respectively is: Hirt enlightenment method, direct method (or matrix method) and energy inequality method.Key words: partial differential equation; finite difference scheme; stability1 前言微分方程的定解问题就是在满足某些定解条件下求微分方程的解。
一种基于全局优化的交错网格有限差分法
一种基于全局优化的交错网格有限差分法印兴耀;刘博;杨凤英【摘要】在地震波场数值模拟中，交错网格有限差分技术得到了广泛的应用，但是在弹性模量变化较大时，通常会因插值而导致模拟误差增大。
旋转交错网格可以很好地克服这个缺点，因而适合于各向异性介质正演模拟。
但是对于同样大小的网格单元，旋转交错网格需要的步长比常规交错网格要大，这会使梯度和散度算子的误差增大因而更易产生空间数值频散。
针对这些问题，本文提出了旋转交错网格与紧致有限差分相结合的方法，并基于模拟退火算法进行全局优化，压制数值频散，拓宽波数范围。
数值模拟结果表明，此方法可以有效地压制数值频散，且具有较高的模拟精度。
%Staggered grid finite difference techniques have been widely used in numerical simulation of seismic wave field,but the interpolation might enlarge the simulation errors when the elastic moduli change dramatically.Rotated staggered grid can overcome this drawback thus be more applicable for modeling the anisotropic media.However,for the grid cells with the same size,rotated staggered grid needs larger step than conventional staggered grid does,conse-quently,the gradient operator and divergence operator generate larger errors which are more likely to cause spatial numerical dispersion.Aiming at the above problems,a combination of rotated staggered grid and compact finite difference method is proposed,simultaneously,global optimization is preformed based on simulated annealing algorithm to suppress numerical dispersion and to broaden the range of wavenumber.Numerical simulation results show thatthis method can effectively suppress the numerical dispersion and exhibits high simulation accuracy at the same time.【期刊名称】《地震学报》【年(卷),期】2015(000)002【总页数】11页(P278-288)【关键词】全局优化;旋转交错网格;紧致有限差分;数值频散;正演模拟【作者】印兴耀;刘博;杨凤英【作者单位】中国山东青岛 266580 中国石油大学华东地球科学与技术学院;中国山东青岛 266580 中国石油大学华东地球科学与技术学院;中国山东青岛 266580 中国石油大学华东地球科学与技术学院【正文语种】中文【中图分类】P315.3+1引言地震波数值模拟是在地震波传播理论的基础上，通过数值计算来模拟地震波在地下介质中的传播（董良国，2003），是研究地震波传播特性与地球介质参数关系的重要手段.交错网格有限差分技术已经广泛应用于各向异性介质（Igel et al，1995；Collino，Tsogka，2001；裴正林，王尚旭，2005；殷文等，2006；何燕，2008）、孔隙介质（Dai et al，1995；Zeng et al，2001；孙林洁，2012）、黏弹性介质（白晓寅，2008；孙成禹等，2010）的波动方程模拟中，但是常规交错网格在模拟非均匀性较强的介质和各向异性介质时需对介质参数进行平均或内插（陈浩等，2006），低阶插值可能会导致精度降低，高阶插值则会导致计算量增大.Gold等（1997）提出了旋转交错网格有限差分算法，Saenger等进一步研究了旋转交错网格的稳定性与频散的关系（Saenger et al，2000；Saenger，Bohlen，2004）.在旋转交错网格中采用沿网格对角线差分的方式，可以避免弹性模量的平均与插值，因而更准确.随后一些研究人员通过旋转交错网格的完全匹配层（perfectly matched layer，简写为PML）边界条件的实现，利用旋转交错网格技术对各向异性介质（李敏，刘洋，2012）和黏弹介质（严红勇，刘洋，2012）等进行了数值模拟.但由于旋转交错网格中差分是沿对角线进行，步长是同等条件下常规交错网格的■2倍（网格是正方形的情况下），若用常规的差分格式会导致离散的梯度算子和散度算子在计算时误差增大，因而容易出现数值频散.而如果使用高阶有限差分则会因需要更多的网格点而增大计算量与内存.紧致有限差分是一种隐式差分格式，在同等网格数的条件下，具有比常规差分格式更高的精度以及更低的数值频散，故成为目前关注的热点.Lele（1992）分析了紧致差分格式的分辨率特性；王书强等（2002）对弹性波方程的紧致差分格式及其与中心差分的误差进行了研究；Du等（2009）基于交错网格紧致有限差分进行了横向各向同性介质的正演模拟；Boersma（2011）利用6阶紧致有限差分求解了Navier-Stokes方程.实现有限差分时，差分系数的确定会影响数值模拟的精度.有限差分系数的求取一般包括泰勒展开和最优化算法两种方法（Liu，2013）.基于泰勒展开得到的差分系数在波数较小时精度高，但随波数增大其精度会降低；基于最优化算法的差分系数求取方法是在给定精度与差分算子长度的情况下尽可能拓宽波数范围（Holberg，1987）.Shan（2009），Kosloff等（2010），Zhou和Zhang（2011）以及Liu（2013）等对基于最小二乘算法的差分算子优化方法进行了深入研究；Zhang 和Yao（2012）则基于模拟退火算法对中心差分算子进行了优化.本文拟结合旋转交错网格与紧致有限差分技术，基于模拟退火全局优化算法对紧致差分算子进行优化，拓宽数值模拟的波数范围，在此基础上进行频散分析，最后通过数值模拟验证该方法的可行性和有效性.1 旋转交错网格紧致有限差分方法1.1 一阶速度－应力方程在弹性波正演模拟中，除了使用二阶方程外，还常常采用一阶速度－应力弹性波方程，其主要优点是无需对弹性常数进行空间差分（Virieux，1984）.这里根据运动平衡方程和本构关系给出二维情况下体力为零时各向异性介质的一阶速度－应力方程：式中：v，σ分别代表速度和应力，1表示x方向，3表示z方向；cij为介质弹性张量矩阵的元素，当极端各向异性介质时，该矩阵有21个独立参数；当二维时，不考虑y分量，故张量矩阵中与y分量有关参数为0，只有5个独立参数（c11，c13＝c31，c15＝c51，c33，c35＝c53，c55）.1.2 紧致有限差分在地震波正演模拟中，采用常规的中心有限差分往往需要较多的网格点才能比较有效地压制数值频散，不利于边界的处理（王书强等，2002）.紧致有限差分是一种隐式差分格式，只需要较少的网格点即可有效地压制数值频散，且同等阶次的精度要比常规中心有限差分格式高，因此本文采用紧致有限差分格式，弥补旋转交错网格因差分步长大导致梯度和散度算子在计算时误差增大因而更容易产生数值频散的不足.常规的2（M－1）点2 M阶紧致有限差分格式（Kim，Lee，1996）为式中，Δx为网格的空间步长，a和bn分别为差分的权系数，可以通过在第i个网格点的泰勒展开，然后对比对应系数来求解即可得到2（M－1）点2 M阶紧致有限差分格式的差分系数（Liu，Sen，2009），例如当M＝5时，可以得到8点10阶的紧致差分格式系数为a＝0.257 894 74，b1＝0.889 871 16，b2＝0.216 121 16，b3＝－0.004 701 2，b4＝0.000 151 55.1.3 旋转交错网格及完全匹配层吸收边界条件在二维旋转交错网格上，密度和速度的各个分量定义在相同位置，弹性模量和应力的各个分量定义在另外的相同位置上，如图1所示.旋转交错网格在计算时分为两步（陈浩等，2006）：首先，计算沿着对角线方向即˜x，˜z方向进行差分，得到相关物理量的对角线方向的空间一阶导数；然后，通过坐标旋转的换算关系，将得到的两个对角线方向的差分进行线性组合从而获得坐标轴方向即x，z方向的一阶空间导数，其换算关系为（Saenger et al，2000）：图1 旋转交错网格及完全匹配层吸收边界示意图x，z为坐标轴方向，x˜，z˜为对角线方向；灰色区域表示沿着x和z方向均进行衰减，白色区域表示只沿着x方向或者z方向进行衰减Fig.1 Schematic diagram of rotated staggeredgridand perfectly matched layer（PML）absorbing boundary xand zare the coordinate directions，andx˜andz˜ are the diagonal directions.Gray areas indicate that waves are absorbed along both xand z directions and white areas indicate that waves are absorbed only along the xor zdirection式中，∂／∂x˜和∂／∂z˜是沿对角线方向的导数；∂／∂x和∂／∂z是沿坐标轴方向的导数.利用波动方程进行数值模拟，一个关键问题就是边界条件.为了解决边界反射问题，引入弹性波方程的PML边界条件.陈浩等（2006）指出，在旋转交错网格中，PML吸收边界条件的处理方式及吸收效果与常规交错网格中几乎是一样的.为此，引入图1所示的PML吸收层.PML边界的基本做法是在研究区域四周引入PML，波在PML中传播时不会产生反射，并且随传播距离按一定规律衰减.当波传播到PML边界时，波场近似为零，也不会产生反射（王守东，2003）.以式（1）中的vx为例来说明旋转交错网格紧致有限差分算法PML的实现方法：式中表示n＋1／2时刻vx 在网格点（i，j）处沿坐标轴x方向的分量，表示n－1／2时刻的x方向分量，di为衰减因子，表示n 时刻σxx在网格点（i，j）处沿坐标轴x方向的导数.根据式（5）可推导出式中沿着所在对角线˜x与斜对角线˜z方向的导数值，根据求出.其中类似地可以推导出其它分量的旋转交错网格紧致有限差分下的PML控制方程. 1.4 紧致有限差分的全局优化根据平面波理论，令式中，β＝kΔx／2，0≤β≤π／2.为了使式（2）左右两边的误差尽可能小，须满足β－β＊≈0的β取值范围应尽可能地大.本文中，采用模拟退火全局最优算法（Kirkpatrick et al，1983）来优化旋转交错网格紧致有限差分算子.为此，建立以下目标函数：Kim和Lee（1996）指出，当θ接近最大值π／2时，会出现很多难以控制的误差，从而导致优化效果不佳，因此本文选择θ＝rπ／2，其中0＜r＜1.表1列出了优化的8点10阶，10点12阶，12点14阶以及14点16阶的紧致有限差分系数.表1 优化的10—16阶紧致有限差分系数Table 1 Optimized coefficients of 10—16-order compact finite-difference阶次 r a b1 b2 b3 b4／10－2 b5／10－5 b6／10－3 b7／10－5 10 0.85 0.384 0.728 0.369 －0.015 0.152 12 0.87 0.390 0.720 0.376 －0.016 0.171 －6.49 14 0.89 0.434 0.658 0.434 －0.022 0.361 －74 0.138 16 0.85 0.452 0.633 0.458 －0.0252 0.452 －114 0.322 －7.792 频散分析由式（12）可以看出，β接近于β＊的程度表征了优化的紧致有限差分算子的精度，故定义使α≈1的β取值范围越大，数值频散越小.图2给出了常规的基于泰勒展开得到的差分算子与本文提出的全局优化算子的精度对比.可以看出，对于常规的基于泰勒展开得到的紧致有限差分算子，随着空间阶次的增大，使α≈1的波数范围变得越宽，即数值模拟的精度随阶次的增大而增大.显然，本文得到的全局优化的紧致有限差分算子在相同最大误差下使α≈1的波数范围比常规紧致有限差分算子还要大.从图2中还可以看出，优化的8点10阶紧致有限差分算子的频散比常规14点16阶紧致有限差分算子小，即优化的10阶精度要高于常规的16阶精度，从而可以使用较少的网格点来达到较高的精度，节省计算内存.图2 10阶常规有限差分，10—16阶常规紧致有限差分及对应的优化紧致有限差分的频散对比（k代表波数，Δx表示空间网格间距）Fig.2 Dispersion comparison by 10-order conventional finite-difference，10—16-order conventional compact finite-difference and the corresponding optimized compact finitedifference.k denotes the wavenumber，andΔxis spatial grid spacing由图2中的10阶常规有限差分、10阶常规紧致有限差分与10阶优化紧致有限差分的频散曲线可以看出，10阶常规紧致有限差分格式的波数范围较10阶常规有限差分要宽，而经过优化的10阶紧致有限差分格式的波数范围则更进一步拓宽.设最大允许误差为εmax＝｜α－1｜0≤β≤βmax，取εmax＝0.2%，分别计算得到10阶常规有限差分、10阶常规紧致有限差分及10阶优化紧致有限差分的最大波数范围βmax分别为0.847 5，1.037 5，1.458 0，即10阶常规紧致有限差分的波数范围比常规差分格式拓宽了1.22倍，而经过全局优化后，波数范围提升至1.72倍.考虑到正演模拟中最常采用Δx＝Δz的情况，则对于旋转交错网格，Δr比常规交错网格差分步长Δx及Δz增大了2倍，因而可以采用优化系数的紧致有限差分来弥补这一不足，从而有效地压制数值频散.需要指出的是，优化系数的紧致有限差分格式也可以很方便地应用于PML边界条件，而不需要特殊的处理.3 模型试算3.1 模型1为了验证全局优化的旋转交错网格紧致有限差分的精度比常规旋转交错网格有限差分方法有所提升，首先采用一个简单的各向同性模型进行测试（模型大小为3 000m×3 000m，在测试中并未添加任何边界条件），震源采用主频30Hz的雷克子波，z方向集中力源激发，网格大小为Δx＝Δz＝10m，时间采样为Δt＝1ms.纵横波速度分别为4 000m／s和2 600m／s，模型密度为2 300kg／m3.图3给出了3种有限差分条件下400ms时刻x分量波场快照，均采用旋转交错网格系统.图3a采用常规10阶常规有限差分格式；图3b采用8点10阶，10点12阶，12点14阶，14点16阶常规紧致有限差分格式；图3c采用优化的8点10阶紧致有限差分格式.可以看出，10阶常规有限差分与10阶紧致有限差分对比，后者频散要小；而对于图3b中的紧致有限差分格式，随着空间阶次的增大，其精度越高，频散越小；对比图3b与图3c可以看出，优化的8点10阶紧致有限差分的数值频散甚至比常规14点16阶紧致有限差分还要小，这与图2所示是一致的.图4给出了从图3中提取的z＝500m处（即50个网格点处）的波形曲线，以此来对比以上3种差分格式的数值频散.可以看出，10阶常规有限差分算法在波形外出现了数值抖动，说明数值频散比较高，模拟精度相对较低，而10阶紧致有限差分和优化的10阶紧致差分精度较高，从虚线框中可以进一步看出优化的10阶紧致有限差分算法的精度更高一些.图3 10阶常规有限差分（a），10—16阶常规紧致有限差分（b）与优化的10阶紧致有限差分（c）得到的波场快照Fig.3 Snapshots of wave field by using 10-order conventional finite-difference（a），10－16-order conventional compact finite-difference（b）and optimized 10-order compact finite-difference（c）3.2 模型2图4 10阶常规有限差分，10阶常规紧致有限差分及10阶优化紧致有限差分波形对比图Fig.4 Comparison of waveforms by using 10-order conventional finite-difference（FD）（dark line），10-order conventional compact FD （blue line）and optimized 10-order compact FD （red line）模型2的几何构造如图5所示，其中凹陷构造上部及下部均为各向同性介质，最下层为横向各向同性（VTI）介质.模型纵向和横向长为3 000m，具体参数见表2.采用优化的旋转交错网格紧致有限差分方法进行正演模拟，网格大小为Δx＝Δz＝10m，时间采样为Δt＝1ms.震源采用30Hz雷克子波，在模型的（1 150m，1 500m）处以纵波源形式激发.模拟过程中，采用PML边界条件（PML层部分同样采用优化的旋转交错网格紧致有限差分算法），厚度层为30个网格点.为体现吸收效果包含了PML层，图6给出了600ms时刻x及z分量的波场快照.可以看出，震源激发的纵波传播至分界面时，发生了反射与透射，可以见到PP波、PS波等，并可在凹陷处看到断面波，波场复杂，但波前面清晰，数值频散小.从图7所示的单炮记录中也可以看出优化的旋转交错网格紧致有限差分算法具有较高的模拟精度，各类波均可得到清晰体现，并且能够证明所加PML边界条件的吸收效果很好，无明显人为边界反射产生.图5 凹陷模型示意图Fig.5 Schematic diagram of the depression model图6 600ms时刻优化的旋转交错网格紧致有限差分正演x分量（左）和z分量（右）的波场快照（虚线框外围代表PML层）Fig.6 Wavefield snapshots of x-（left）and z-component（right）at 600ms using optimized rotated staggered-grid compact finite-difference method（PML layers are outside the dashed box）表2 凹陷模型参数Table 2 Parameters of the depression model序号介质类型 c11／GPa c13／GPa c33／GPa c55／GPa 密度／（kg·m－3）31.77 42.68 31.77 13.75 2200第二层各向同性 42.33 47.04 42.34 18.82 2400第三层横向各向同性第一层各向同性87.88 21.22 67.60 22.56 2500图7 凹陷模型优化的旋转交错网格紧致有限差分正演x分量（左）和z分量（右）单炮记录Fig.7 Records of x-（left）and z-component（right）of the depression model using optimized rotated staggered-grid compact finite-difference method4 讨论与结论旋转交错网格在地震波正演模拟中由于其可以避免在模拟中因弹性模量插值或平均而产生的误差和其对边界条件处理的便捷性而被广泛应用.但由于其差分方向是沿网格的对角线方向，因此相对于常规的交错网格而言，其差分步长较大，梯度算子、散度算子在计算时更容易产生误差，这意味着旋转交错网格更容易产生数值频散，从而影响模拟的精度.若采用小网格，虽然可以压制数值频散，却又大大增加了运算量及内存.因此本文采用旋转交错网格与紧致有限差分技术相结合的方法来提高模拟精度，并在此基础上，基于模拟退火算法对紧致有限差分进行全局优化，进一步提高计算精度，压制数值频散.数值频散分析结果表明，优化的旋转交错网格紧致有限差分算法相对于普通旋转交错网格有限差分算法具有更宽的波数范围，这意味着优化的旋转交错网格有限差分算法可以在采用更大的空间采样间隔或者更高的震源主频时依然有较高的精度和较低的数值频散；且经过全局优化的10阶紧致有限差分算子比常规的16阶紧致有限差分算子具有更高的精度，即可以采用较少的网格点来达到较高的模拟精度，因此可以节省计算内存.数值模拟实验进一步证实了本文提出的优化的旋转交错网格紧致有限差分算法的正确性与可行性.凹陷模型中，各类波均具有比较清晰的响应，揭示了波场的传播规律.同时波场快照与单炮记录中无明显的人为边界反射产生，也证明了模拟采用的PML吸收边界条件的有效性.但是本文提出的方法亦存在不足之处.在求解紧致差分格式时需要解三角矩阵线性方程，较常规显式有限差分格式增加额外的计算负担，且不利于在GPU等平台直接并行，例如图3a和图3c中的400ms时长的计算时间分别为46.5s和85.6s，后者计算耗时约为前者的1.8倍.在科学计算中，效率与计算精度常常不可兼得，随着计算机科学计算能力的提高，可以采用GPU加速的LU分解算法或者使用LAPACK等优化程序包来克服求解大型矩阵的计算效率问题.在实际应用中需权衡利弊，选择适合的算法.衷心感谢审稿专家提出的宝贵意见和建议.参考文献白晓寅.2008.基于地震波衰减理论的地层吸收参数提取方法研究［D］.东营：中国石油大学（华东）地球科学与技术学院：24-54.Bai X Y.2008.The Study on Method of Layer Absorbing Parameters Extraction Based on the Theory of Seismic Wave Attenuation［D］.Dongying：School of Geosciences，China University of Petroleum：24-54（in Chinese）.陈浩，王秀明，赵海波.2006.旋转交错网格有限差分及其完全匹配层吸收边界条件［J］.科学通报，51（17）：1985-1994.Chen H，Wang X M，Zhao H B.2006.Rotated staggered grid finite-difference and the PML boundary［J］.Chinese Science Bulletin，51（17）：1985-1994（in Chinese）.董良国.2003.地震波数值模拟与反演中几个关键问题研究［D］.上海：同济大学海洋与地球科学学院：1-2.Dong L G.2003.Several Key Problems in Seismic Wave Numerical Simulation and Inversion［D］.Shanghai：School of Ocean and Earth Science，Tongji University：1-2（in Chinese）.何燕.2008.正交各向异性弹性波高阶有限差分正演模拟研究［D］.东营：中国石油大学（华东）地球科学与技术学院：45-80.He Y.2008.High-Order Finite-Difference Forward Modeling of Elastic-Wavein Orthorhombic Anisotropic Media［D］.Dongying：School of Geosciences，China University of Petroleum：45-80（in Chinese）.李敏，刘洋.2012.高阶旋转交错网格有限差分方法模拟TTI介质中横波分裂［J］.物探与化探，36（6）：934-940.Li M，Liu Y.2012.Modeling of the S-wave splitting in TTI media using high-order rotated staggered grid scheme［J］.Geophysical and Geochemical Exploration，36（6）：934-940（in Chinese）.裴正林，王尚旭.2005.任意倾斜各向异性介质中弹性波波场交错网格高阶有限差分法模拟［J］.地震学报，27（4）：441-451.Pei Z L，Wang S X.2005.A staggered-grid high-order finite-difference modeling for elastic wave field in arbitrary tilt anisotropic media［J］.Acta Seismologica Sinica，27（4）：441-451（in Chinese）.孙成禹，肖云飞，印兴耀，彭洪超.2010.黏弹介质波动方程有限差分解的稳定性研究［J］.地震学报，32（2）：147-156.Sun C Y，Xiao Y F，Yin X Y，Peng H C.2010.Study on the stability of finite difference solution of visco-elastic wave equations［J］.Acta Seismologica Sinica，32（2）：147-156（in Chinese）.孙林洁.2012.非均匀孔隙介质有限差分正演模拟方法研究［D］.青岛：中国石油大学（华东）地球科学与技术学院：64-108.Sun L J.2012.Finite-Difference Modeling Research in Heterogeneous Porous Media［D］.Qingdao：School of Geosciences，China University of Petroleum：64-108（in Chinese）.王守东.2003.声波方程完全匹配层吸收边界［J］.石油地球物理勘探，38（1）：31-34.Wang S D.2003.Absorbing boundary condition for acoustic wave equation by perfectly matched layer［J］.Oil Geophysi-cal Prospecting，38（1）：31-34（in Chinese）.王书强，杨顶辉，杨宽德.2002.弹性波方程的紧致差分方法［J］.清华大学学报：自然科学版，42（8）：1128-1131.Wang S Q，Yang D H，Yang K pact finite difference scheme for elastic equations［J］.Journal of Tsinghua University：Science and Technology，42（8）：1128-1131（in Chinese）.严红勇，刘洋.2012.黏弹TTI介质中旋转交错网格高阶有限差分数值模拟［J］.地球物理学报，55（4）：1354-1365.Yan H Y，Liu Y.2012.Rotated staggered grid high-order finite-difference numerical modeling for wave propagation in viscoelastic TTI media ［J］.Chinese Journal of Geophysics，55（4）：1354-1365（in Chinese）. 殷文，印兴耀，吴国忱，梁锴.2006.高精度频率域弹性波方程有限差分方法及波场模拟［J］.地球物理学报，49（2）：561-568.Yin W，Yin X Y，Wu G C，Liang K.2006.The method of finite difference of high precision elastic wave equations in the frequency domain and wave field simulation［J］.Chinese Journal of Geophysics，49（2）：561-568（in Chinese）.Boersma B J.2011.A 6th order staggered compact finite difference method for the incompressible Navier-Stokes and scalar transport equations［J］.J Comput Phys，230（12）：4940-4954.Collino F，Tsogka C.2001.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropicheterogeneous media［J］.Geophysics，66（1）：294-307.Dai N，Vafidis A，Kanasewich E R.1995.Wave propagation in heterogeneous，porous media：A velocity-stress，finitedifference method ［J］.Geophysics，60（2）：327-340.Du Q Z，Li B，Hou B.2009.Numerical modeling of seismic wavefields in transversely isotropic media with a compact staggered-grid finite difference scheme［J］.Appl Geophys，6（1）：42-49.Gold N，Shapiro S A，Burr E.1997.Modelling of high contrasts in elastic media using a modified finite difference scheme［C］∥SEG Technical Program Expanded Abstracts.Dallas，Texas：SEG：1850-1853.Holberg putational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena［J］.Geophys Prospect，35（6）：629-655.Igel H，Mora P，Riollet B.1995.Anisotropic wave propagation through finite-difference grids［J］.Geophysics，60（4）：1203-1216.Kim J W，Lee D J.1996.Optimized compact finite difference schemes with maximum resolution［J］.AIAA J，34（5）：887-893.Kirkpatrick S，Gelatt C D Jr，Vecchi M P.1983.Optimization by simulated annealing［J］.Science，220（4598）：671-680.Kosloff D，Pestana R C，Tal-Ezer H.2010.Acoustic and elastic numerical wave simulations by recursive spatial derivative operators［J］.Geophysics，75（6）：T167-T174.Lele S pact finite difference schemes with spectral-like resolution［J］.J Comput Phys，103（1）：16-42.Liu Y，Sen M K.2009.An implicit staggered-grid finite-difference method for seismic modelling［J］.Geophys J Int，179（1）：459-474.Liu Y.2013.Globally optimal finite-difference schemes based on least squares［J］.Geophysics，78（4）：T113-T132.Saenger E H，Bohlen T.2004.Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid［J］.Geophysics，69（2）：583-591.Saenger E H，Gold N，Shapiro S A.2000.Modeling the propagation of elastic waves using a modified finite-difference grid［J］.Wave Motion，31（1）：77-92.Shan G.2009.Optimized implicit finite-difference and Fourier finite-difference migration for VTI media［J］.Geophysics，74（6）：WCA189-WCA197.Virieux J.1984.SH-wave propagation in heterogeneous media：Velocity-stress finite-difference method［J］.Geophysics，49（11）：1933-1957. Zeng Y Q，He J Q，Liu Q H.2001.The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media ［J］.Geophysics，66（4）：1258-1266.Zhang J H，Yao Z X.2012.Optimized finite-difference operator for broadband seismic wave modeling［J］.Geophysics，78（1）：A13-A18. Zhou H B，Zhang G Q.2011.Prefactored optimized compact finite-difference schemes for second spatial derivatives［J］.Geophysics，76（5）：WB87-WB95.。
三维波动方程有限差分正演方法
三维波动方程有限差分正演方法三维波动方程是描述地震波传播和地震波场的方程，是地震勘探、地震监测等领域中常用的数值模拟方法。
三维波动方程的有限差分正演方法是一种常用的数值解法，通过将连续的偏微分方程离散化，转化为差分方程进行求解，可以得到地震波场的时空分布情况。
∂²P/∂t²=c²(∂²P/∂x²+∂²P/∂y²+∂²P/∂z²)其中，P是地震波场的压力波动，t是时间，c是地震波传播速度，x、y、z是空间坐标。
有限差分正演方法通过离散化空间和时间进行数值求解。
对空间进行离散化，将地震波场的x、y、z坐标分别划分为Nx、Ny、Nz个网格点，其中Nx、Ny、Nz分别表示x、y、z方向的网格数。
对时间进行离散化，将t划分为Nt个时间步长，其中Nt表示时间步数。
将地震波场的压力P表示为P(i,j,k,n)，其中i、j、k表示网格点坐标，n表示时间步长。
在有限差分正演方法中，采用中央差分格式对空间导数进行离散化，采用二阶精度的差分格式对时间导数进行离散化，得到如下差分方程：P(i,j,k,n+1)=2P(i,j,k,n)-P(i,j,k,n-1)+c²Δt²(∂²P/∂x²+∂²P/∂y²+∂²P/∂z²)其中，Δt是时间步长。
上述差分方程可以通过迭代求解，从而得到地震波场的时空分布。
有限差分正演方法的优点是简单、直观，容易编程实现。
但是，由于空间和时间的离散化会引入数值误差，因此需要合理选择网格大小和时间步长，以保证数值解的精度和稳定性。
此外，由于三维方程的求解计算量较大，所以需要进行高效的并行计算。
总结来说，三维波动方程的有限差分正演方法是一种常用的地震波场数值模拟方法，通过离散化地震波场的空间和时间，将偏微分方程转化为差分方程进行求解。
地震波数值模拟与分析
地震波数值模拟与分析地震波是地震活动中最重要的研究对象之一。
而地震波数值模拟和分析则是地震学领域中的重要研究方向之一。
在地震波数值模拟和分析的过程中，人们可以通过计算机模拟地震波的传播过程，并从中获取有关地震特征及其引起的地表破坏和建筑物结构变形等各种信息。
这对于地震灾害的预防、预测和减轻有着重要的意义。
地震波的数值模拟方法主要有有限差分法、有限元法、边界元法和谱元法等。
其中，有限差分法是目前地震波数值模拟中应用最为广泛的一种方法。
有限差分法在解决非线性、多维度和非静态问题方面表现尤为出色。
其基本思想是将地震波场离散成网格，并利用二阶精度差分公式计算各个时刻在网格点处的地震波场值。
有限差分法的优点在于精度高、计算速度快，同时可以对复杂地质构造及其他复杂条件进行模拟分析。
地震波的数值分析方法主要有PTA和TFI等。
其中，PTA是计算地震波传播中频谱组成的一种方法。
PTA方法基于傅里叶变换，将地震波在频域中进行分析，主要考虑波振幅和频率之间的关系。
通过对地震波的频谱进行分析，可以得出波传播路径、应变速率及层间的速度等信息。
而TFI则是通过时间域内的雷克子波分析地震波的能量分布，从而得出地表加速度和地震破坏信息。
当我们研究地震波数值模拟的同时，还要重视地震波分析的意义。
地震波的分析能够帮助我们对地震发生的原因、机制及它们对地表的影响进行研究。
同时，地震波分析也可以帮助我们评估地震对建筑物和基础设施的破坏。
这项工作通常涉及结构动力学模拟、震害评估、震害预测等研究领域。
此外，通过地震波分析，我们也可以了解地震所带来的生态影响和异常现象（如水波、地陷等）。
在地震波数值模拟和分析过程中，实际数据采集十分必要。
地震数据采集主要分为地震观测和近场强动观测两种方法。
地震观测是通过装置地震仪器等方法获得的数据。
而近场强动观测则是通过现场安装观测设备，获取地震波传播的信息。
同时，人工模拟地震波也是一种可行的方法，但其对于地震波的形态和波速等方面需进行较为精确的估计。
二维地震波场有限差分法数值模拟研究的开题报告
二维地震波场有限差分法数值模拟研究的开题报告一、研究背景地震学是研究地震现象及其产生的原因、规律、预测方法和应对措施的一门学科。
在地震学中，地震波场模拟是一个重要的研究方向。
通过数值模拟地震波场，可以模拟地震过程中的地震波传播规律，预测地震波的传播方向和影响范围，为地震预测和预警提供参考。
二、研究目的和意义本研究旨在采用有限差分法对二维地震波场进行数值模拟研究，探究地震波传播规律和影响范围，为地震预测和预警提供重要参考。
三、研究方法和步骤本研究将采用有限差分法对二维地震波场进行数值模拟研究。
具体步骤如下：1. 建立二维地震波场的数学模型。
2. 设计有限差分算法，对地震波场进行数值模拟计算。
3. 根据模拟结果进行分析，并与实际地震数据进行对比。
4. 优化模型和算法，提高模拟精度和计算效率。
四、预期结果通过本研究，我们将得到二维地震波场的传播规律和影响范围，并能够对地震波进行预测。
该研究结果将为地震预测和预警提供重要参考，具有重要的科学研究和实践应用价值。
五、研究难点和关键技术本研究面临的主要难点是算法的设计和优化，需要针对地震波场的特性进行合理的数学建模和算法设计，并进行多次优化和调试，以提高模拟精度和计算效率。
关键技术包括数值计算方法、地震波场建模、数值优化等方面的知识。
六、研究计划本研究的时间安排和具体任务分配如下：第一年1. 确定研究方向和目标。
2. 学习和掌握地震波场建模和有限差分法算法。
3. 设计二维地震波场模型和有限差分算法。
第二年1. 完成地震波场的数值模拟计算。
2. 分析模拟结果，与实际地震数据进行对比。
3. 优化模型和算法，提高模拟精度和计算效率。
第三年1. 进一步优化模型和算法。
2. 撰写论文并准备发表。
七、研究团队本研究由一名硕士研究生和一名导师组成。
硕士研究生主要负责具体的研究和实验操作，导师负责指导和协助研究工作。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
#! ・ ・ ・
# 4 )"6
#
差分模拟举例
根据所建立的差分方程，编制了岩
岩层 #"3
#" 4 )"6
( 4 "L # 8 &" ’ ! 3# 9 :
床侵入煤系时不稳定地温场有限差分 模拟的计算程序， 并进行了运算。 设某一研究区域， 见图 # 。 岩床半宽
! &
岩床 ) 3 图#
#" 4 5)"6
通常， 煤系主要接受地壳深部传导上来的热量， 并通过地面向大气散发； 大都情况下， 煤系 中的热交换已处于动态平衡状态， 可采用稳定温度场来研究。 但在有些情况下， 如： 由于岩浆侵 入、 煤的自燃等， 煤系温度的升高和降低， 就不能用稳定温度场来研究， 这时必须采用不稳定地 温场来研究。对不稳定地温场的研究具有意义， 如岩浆侵入的地温场， 对于分析岩浆侵入的影 响范围和程度、 分析地热历史、 分析有机质的成熟度等都具有意义。 目前， 应用广泛的数值模拟方法包括有限差分法和有限单元法； 本文以岩床为例， 采用有 限差分法对不稳定地温场作一分析。
图% 各结点温度随时间变化图
123456 % 7852892:; :< 68=) ;:/6 96>?6589456 @29) 92>6
图-
.、 0、 !% 、 "" 结点温度随时间变化 ( 时间不取对数 * 图.
96>?6589456’ @29) 92>6
123456 - 7852892:; :< A:B .C 0C !% 8;/ "" ;:/6
收稿日期： -../ 5 .8 5 .9 基金项目： 安徽省高等学校优秀青年教师资助计划项目 1 编号： -...=>/-/ 2 作者简介： 赵志根 1 /86? 5 2 ， 男， 浙江金华人， 博士， 副教授， 从事煤田地质教学、 研究
1/2
!"
安
徽
地
质
#""# 年
!
差分方程的建立
所谓差分方法， 就是把热传导的微分方程近似地改用差分方程 $ 代数方程 % 来表示， 把求解
$&$ $
4 #" ; $ ( ( % $ $
$ #" ( &; $ ( #" ’ &; $ ’ # #"; $ % !#
$!%
图&
差分方程推导示意图
.E AF@ C0EE@/@BJ@ @K?+A0.B
=0>?/@ & G0+>/+3 :F.H0B> AF@ C@/0I+A0.B
"
差分方程的求解
根据 $ ! % 式， 在参数 (%、 只要已知各结点的初始温度， 根据结点 " ’ &、 $ $、 ! 已确定时， "、 " ( & 在 $ 时刻的温度，就可求得结点 " 在 $ $ & $ $ % 时刻的温度；同理，可以求出所有结点在 $ $ ( $ $ % 时刻的温度。
!"#$% &’ #’!"()*+ "+,-+.("#.+ /0+*$ )% /0’0"+ $0//+.+’1+ ,+"2&$
LMNO L)2P36;
( "#$%&’(#)’ *+ ,#-*.&/#- %)0 1)23&*)(#)’%4 1)53)##&3)56 7.%3)%) 8)-’3’.’# *+ !#/9)*4*5:6 7.%3)%) "%"&&! 6 ;93)% * NQ’958=9R S)6 ’94/T :; 96>?6589456 <26D/ :< =:8D >2;6 2’ 4’6<4D 2; >8;T 8’?6=9’, U;’98QD6 96>?6589456 <26D/ 2’ 56’685=)6/ QT <2;296 /2<<656;=6 >69):/ 2; 9)2’ ?8?65C 2;=D4/2;3 9)6 6’98QD2’)>6;9 :< /2<<656;=6 6V4892:; 8;/ 9)6 /6W2=6 :< ?5:358>, N;/ 8D’:C 8; 6X8>?D6 8Q:49 ’2DD 2’ 32W6; 9: 2DD4’95896 9)6 56’685=), ;4>652=8D >:/6D； <2;296 /2<<656;=6 >69):/ ； ’2DD Y6T @:5/R 4;’98QD6 96>?6589456 <26D/ ；
$&$ $
’ #" ; $$
$
$2%
"(& " "’&
— 第 0 结点在 $ $ ( $ $ % 时刻的温度 式中： #"，$ ( $ $—— — 第 " ( &、 #" & & ，$， #"，$， #" ’ &，$—— "、 " ’ & 结点在 $ 时刻的温度 — 时间步长 $ $—— $ 2 % 代入 $ & % 式， 可得结点 " 的差分公式： 将$1%、 #" ;
( 4 &7 ! 8 &" ’ ! 3# 9 :
研究区特征及剖分图
为 23，初始温度 #" 4 5)"6 ，热导率 =0>?/@ # =@+A?/@ +BC D/.E0-@ .E AF@ :A?C0@C +/@+ ’! # ， 下边界为对称边界； 岩层厚度为 初始温度 #" 4 )"6 ， 热导率 ( 4 ( 4 &7 ! 8 &" 3 9: #"3， 岩层外环境为常温 )"6 。 取网格间距 ! 为 & 3， 从下往上编号， 共 #! 个结点。 "7 # 8 &" ’ ! 3# 9 :，
微分方程的问题改换为求解代数方程的问题。 用相隔等间距 ! 而平行于坐标轴的一组平行线织成差分网络，见图 &；温度为一连续函 数，" ’ &、"、" ( & 为连续的 ) 个结点，可将温度 # 在 $ 时刻的结点 " ’ & 和 " ( & 沿 % 方向以 *+,-./ 级数展开： #" ( &; $ 4 #"; $ ( ! $ #" ’ &; $ 4 #"; $ ’ ! $ &# % !# $ & # # % !) $ & ) # % !1 $ & 1 # % " ( " ( " ( " ( …… # ) &% #< & % )< & % 1< & % 1 &# % !# $ & # # % !) $ & ) # % !1 $ & 1 # % " ( " ’ " ( " ’ …… # ) &% #< & % )< & % 1< & % 1 $#% $)%
-..- 年
文章编号： /..4 5 6/47 1 -..- 2 ./ 5 ..48 5 .9
安 徽 地 质 !"#$#%& #’ ()*+,
第 /- 卷 第 / 期
不稳定地温场的有限差分模拟
赵志根
1 淮南工业学院资源与环境工程系， 安徽 淮南 -9-../ 2
摘要： 对矿井地温场的研究具有多方面的意义。本文采用有限差分法对不稳定地温场进行了研究， 根据数学模型建立起差分方程，编制了计算程序，并通过一个实例分析了岩床侵入煤系地层时的温度 变化规律。 关键词： 不稳定地温场: 数学模型: 有限差分法: 岩床 中图分类号： 0;/69 文献标识码： (
E ! F 丁国玺 , 工程热力学及传热学 E G F , 北京： 煤炭工业出版社 , !H0#, I0 J 0& E " F 杨绪灿， 金建三 , 弹性力学 E G F , 北京： 高等教育出版社， !H0I, "&# J ""! E % F K, 邦特巴思 , 地热学导论 E G F , 易志新， 熊亮萍译 , 北京： 地震出版社 , !H00
对于所研究的区域来说，假定 ! 是充分小的，可以不计三次幂及更高次幂的各项； $ # % 、 $ ) % 式相加、 整理， 可得： # #" ( &; $ ( #" ’ &; $ ’ # #"; $ $1% $& # %"4 &% # !# ! 对时间进行差分， 可得： $ #" ; &# %"4 &%
.、 0、 !% 、 "" 结点温度随时间变化 ( 时间取对数 *
由图 %、 岩床的温度 ( . 结点 * 一 -、 . 可知： 直处于降低， 而岩层的温度 ( 0、 !%、 "" 结点 * 经
123456 . 7852892:; :< A:B .C 0C !% 8;/ "" ;:/6 96>?6589456’ @29) 9)6 D:38529)> :< 92>6