地震波方程人工边界的两种处理方法

第42卷第4期2003年12月

石　油　物　探

GEOPHYSICAL PROSPECTIN G FOR PETROL EUM

Vol.42,No.4

Dec.,2003

文章编号:100021441(2003)0420452204

地震波方程人工边界的两种处理方法

崔兴福1,2,张关泉2

(1.中国石油辽河油田分公司勘探开发研究院计算所,辽宁盘锦124010;2.中国科学院数学与系

统科学研究院计算数学与科学工程计算研究所,北京100080)

摘要:在分析以往人工边界处理优缺点的基础上,提出了利用波动方程的频散关系式,得到可以吸收任何方向入射波的自适应校正吸收边界条件和旋转校正吸收边界条件。同时,采用波阵面法和能流密度法判别入射波方向,克服了Pad

e近似吸收边界只对正入射波具有较好吸收性,而对非正入射的波吸收性不好的缺点。数值试验结果表明了本方法的有效性。

关键词:自适应校正;旋转校正;波阵面;能流密度

中图分类号:P63114 文献标识码:A

Two processing methods for artif icial boundary of seismic w ave equation

Cui Xingfu1,2,Zhang Guanquan2

(https://www.360docs.net/doc/9311503098.html,puting Center of Exploration&Development Research Institute of CNPC Liaohe Oilfield Com pany,Panjin

124010,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]、移动边界等,用Pad

e近似法得到的吸收边界条件[4,5]是目前边界处理效果比较好的一种。我们正是在分析这种边界处理优缺点的基础上,从不同的侧面提出了这种边界的2种校正方法,即自适应校正法和旋转校正法。

1　声波方程自适应校正吸收边界条件

1.1　二阶精度吸收边界的推导

由二维声波方程

u t t=c2(u x x+u zz)(1)进行三维Fourier变换得到频散关系式

ω2=c2(k2

+k2z)=k2c2(2)式中,k2=k2x+k2z,k x=k cosθ,k z=k sinθ,θ是波前k和k z的夹角。引进中间参量

a(θ)=

1-cosθ

1-cos

(3)由图1表示的方向波数与波前关系可以推得

cos

(α-1)k+k x

αk(4)利用k x=k cosθ,cosθ=2cos2

1和2

α2-4α+4=

1+cosθ

及方程(4)得到在频率波数域右行波的边界条件

ωk

ω2+1

1+cosθk

=0(5)

收稿日期:20021113;改回日期:20030323。

作者简介:崔兴福(1965—),男,高级工程师,博士在读,现从事地震资料处理方法研究工作。

基金项目:本文得到国家973重点基础研究项目(G199903280)资助。

由Fourier 逆变换得到时间空间域中右行波的边界条件

u x t +1c u t t -c

1+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 -c

2u zz

(10)方程(10)是方程(6)在θ=0的特殊情况,对正入射的地震波吸收较好,而对非正入射的地震波吸收不好。将方程(10)变换到频率波数域中,有

ωc k x -ωc 2+12k 2

(11)旋转方程(11),得

ωc

(k x cos θ+k z sin θ)-ωc

(k z cos θ-k x sin θ)2=0(12)

对应到时间空间域得到旋转后的右边界条件

-u xt -tan θu tz -

1c cos θu t t +c

tan θsin θu x x -c sin θu xz +

c 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 3

u 1,u 2和u 3为位移向量U 的各个分量。

在二维(i ?k )情况下得到的弹性波能流密度具体表达式如下

I =

(2μ+λ)

5u 15x 1 u 1+λ5u 3

5x 3 u 3+μ5u 1

5x 3

+　5u 35x 1 u 3i +μ5u 15x 3+5u 35x 1 u 1

+(2μ+λ)

5u 35x 3 u 3+λ5u 1

5x 1 u 3

k (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 m

i ,j ,

i =1,2,…,I ;j =1,2,…,J ;m =1,2,…,M ,i ,j 和k 分别表示x ,z 和t 方向的离散节点号,则

方程(6)和方程(13)的离散差分格式为

u m +1

i +1,j -u m +1

i ,j -u m

i +1,j +u m

i ,j

Δt Δx

+ 1c

u m +1i ,j -2u m i ,j +u m -1

i ,j

Δt 2

-c

1+cos θ

u m

i ,j +1-2u m

i ,j +u m

i ,j -1

Δz 2

(17)

-u m +1

i +1,j -u m +1

i ,j -u m

i +1,j +u m

i ,j

Δt Δx -tan θu m +1

i ,j +1-u m +1

i ,j -u m

i ,j +1+u m

i ,j

Δt Δz

-1

c cos θu m +1i ,j -2u m i ,j +u m -1

i ,j Δt 2

tan θsin θ

u m i +1,j -2u m i ,j +u m i -1,j

Δx 2

c sin θu m

i +1,j +1-u m

i +1,j -u m

i ,j +1+u m

i ,j

Δx Δz +

2c cos θu m i ,j +1-2u m i ,j +u m i ,j -1Δz 2

=0(18)

3　数值计算比较

图1是方向波数与波前关系示意图。图2a 是

模型示意图,2层介质的速度分别为v 1=3000m/s

和v 2=2300m/s ,密度分别为ρ1=2.3g/cm 3

和ρ2=2.0g/cm 3,

Δx =Δz =5m ,Δt =0.00075s ;图2b ～图2e 为加入不同边界的平面震源时间

剖面,比较可以看出,自适应校正边界具有良好的边界吸收效果。图3是在斜层介质模型上进行点震源模拟得到的仿真时间剖面,模型为2层,介质的速度和密度分别为v 1=3600m/s ,v 2=2500

m/s ,ρ1=2.1g/cm 3,ρ2=2.0g/cm 3。

在图3a 中2边加了Neumann 边界,在图3b 中2边加了旋转

校正吸收边界,比较2图可以看到,旋转校正吸收

边界能吸收大部分的入射波能量。图4是在均匀介质中得到的点震源各个时刻的波场图,速度v =3600m/s ,密度ρ=2.0g/cm 3,Δx =Δz =5m ,Δt =0.00075s 。图4a ～图4e 为2边加Neu 2mann 边界的波场图,图4f ～图4j 为2

边加旋转校

图1　

波数与波前关系示意图

图2　Neumann 边界和自动校正吸收边界平面震源时间剖面

(a )模型;(b )左右都加Neumann 边界;(c )左边加自适应校

正边界,右边加Neumann 边界;(d )左边加Neumann 边界,右边加自适应校正边界;(e )2

边都加自适应校正边界

图3　Neumann 边界(a )和旋转校正吸收边界(b )点源时间剖面

?454?石　油　物　探第42卷

正边界的波场图,从图中可以看到,采用吸收边界对地震波具有较好的吸收作用,消除了人工边界的虚假反射,净化了仿真地震剖面

。

图4　Neumann 边界和旋转校正吸收边界点源各个时刻的波场

(a )40Δt ;(b )120Δt ;(c )160Δt ;(d )180Δt ;(e )250Δt ;(f )80Δt ;(g )120Δt ;(h )170Δt ;(i )200Δt ;(j )250Δt

4　结束语

在地震波正演的数值模拟中,不论是采用有限差分法还是有限元法,人工边界处理都是关系到正演模拟效果好坏的关键,若处理不好将直接影响正演模拟的效果。我们在深入研究了前人成果的基础上,提出了利用频散关系来得到可以吸收任何方向入射波的自适应校正吸收的边界条件和旋转校正吸收的边界条件,从而克服了Pad e 近似吸收边界只对正入射波具有较好吸收性,而对非正入射的波吸收性不好的缺点。从旋转坐标的角度,给出了Pad e 近似边界的另一种修正方案,充分利用了Pad e 近似极好的局部逼近性之优点,同时克服了Pad e 近似不好的整体逼近性之缺点,对入射波具

有很好的吸收性。

参　考　文　献

1　Lysmer J ,Kuhlemeyer R L.Finite dynamic model for in 2finite media[J ].Engineering Mechanics Division of ASEC ,1969,95(7):858～877

2　Cerjan C ,K osloff D ,K osloff R ,et al.A nonreflecting

boundary conditions for discrete acoustic and elastic wave equations[J ].G eophysics ,1980,50(4):705～7083　S ochacki J ,Kubichick R ,G eorge J ,et al.Absorbing

boundary conditions and surface waves [J ].G eophysics ,1987,52(1):60～71

4　Clayton R ,Engquist B.Absorbing boundary conditions for

acoustic and elastic wave equations[J ].Bulletin of Seismic S ociety in America ,1977,66(11):1529～15405　Engquist B ,Majda A.Absorbing boundary conditions for

the numerical simulation of waves[J ].Mathematics Com 2putation ,1977,31(9):629～651

6　Engquist B ,Majda A.Radiation boundary conditions for a 2coustic wave calculations [J ].Communication of Pure and Applied Mathematics ,1979,32(1):313～357

7　Liao Z P ,Wong H L ,Y ang B P ,et al.A transmitting

boundary for transient wave analysis [J ].Scienta Sinica.

1984,27(5):1063～1076

8　Lauence Halpern.Wide 2angle one 2way wave equations

[J ].Acoustics S ociety in America ,1984,35(4):1397～

1404

554?第4期崔兴福等1地震波方程人工边界的两种处理方法

地震波使用说明

地震波使用说明此目录下提供了四类场地土的地震波时程曲线和上海人工波。按照场地土类型（1，2，3或4），选择时程曲线。在定义时程工况时，对于多遇或罕遇地震，按比例调整时程曲线的最大值。中国抗震规范规定，作为抗震计算中底部剪力法和振型分解反应谱法的补充方法，对于特别不规则，特别重要的和较高的结构应采用时程分析法进行多遇地震下的补充计算。可取多条时程曲线的计算结果的平均值与振型分解反应谱法计算结果的较大值。采用时程分析法时，应咱建筑场地类别和设计地震分组选用不少于二组的实际强震记录和一组人工模拟的加速度时程曲线，其平均地震影响系数曲线应与振型分解反应谱法所采用的地震影响系数曲线在统计意义上相符。其加速度时程最大值可按规范中对于多遇和罕遇地震在不同烈度下的值。弹性时程分析时，每条时程曲线计算所得结构底部剪力不应小于振型分解反应谱法计算结果的65％，多条时程曲线计算所得结构底部剪力的平均值不应小于振型分解反应谱法计算结果的80％。可使用弹塑性时程分析法计算罕遇地震下结构的变形。时程分析是一个承受随时间变化的指定荷载结构的逐步动态反应分析，可以是线性或非线性的。此章对时程分析进行一般的描述，特别是线性时程分析。定义时程函数用户可使用“从文件中添加函数”，导入已定义的文本文件，即实测的时程曲线；也可使用程序内置的时程函数。

时程函数定义对话框时程函数定义对话框中的条目解释如下： ?函数名通过在编辑框中直接键入以指定或修改时程函数的名称。 ?函数文件 1.在函数文件域点击浏览按钮以调出一个对话框，在此可找出包含时程函数的文本文件名。注意文件名显示在文件名框中 2.在 "要跳过的标题行" 编辑框中输入一个希望ETABS在文本文件中跳过的行数。 3.在 "每行要跳过的前缀字符" 编辑框中输入一个希望ETABS在文本文件中每行要跳过的字符数。 4.在 "每行的点数" 编辑框中输入一个数告诉ETABS文本文件每行的绘图点数。

地震波的选取方法 (MIDAS内部技术资料)

地震波的选取方法（MIDAS内部技术资料） (GB50011-2001)的 5.1.2条文说明中规定，正确选择输入的地震加速度时程曲线，要满足地震动三要素的要求，即频谱特性、有效峰值和持续时间要符合规定。频谱特性可用地震影响系数曲线表征，依据所处的场地类别和设计地震分组确定。这句话的含义是选择的实际地震波所处场地的设计分组(震中距离、震级大小)和场地类别(场地条件)应与要分析的结构物所处场地的相同，简单的说两者的特征周期Tg值应接近或相同。特征周期Tg 值的计算方法见下面公式(1)、(2)、(3)。加速度有效峰值按建筑抗震设计规范(GB50011-2001)中的表5.1.2-2采用。地震波的加速度有效峰值的计算方法见下面公式(1)及下面说明。持续时间的概念不是指地震波数据中总的时间长度。持时Td的定义可分为两大类，一类是以地震动幅值的绝对值来定义的绝对持时，即指地震地面加速度值大于某值的时间总和，即绝对值｜a(t)｜＞k*g的时间总和，k常取为0.05；另一类为以相对值定义的相对持时，即最先与最后一个k*amax之间的时段长度，k一般取0.3～0.5。不论实际的强震记录还是人工模拟波形，一般持续时间取结构基本周期的5～10倍。说明: 有效峰值加速度EPA＝Sa/2.5（1）有效峰值速度EPV＝Sv/2.5（2）特征周期Tg=2*EPV/EPA（3）

1978年美国ATC－3规范中将阻尼比为5％的加速度反应谱取周期为0.1-0.5秒之间的值平均为Sa，将阻尼比为5%的速度反应谱取周期为0.5-2秒之间的值平均为Sv(或取1s附近的平均速度反应谱)，上面公式中常数2.5为0.05组尼比加速度反应谱的平均放大系数。上述方法使用的是将频段固定的方法来求EPA和EPV，1990年的《中国地震烈度区划图》采用了不固定频段的方法分析各条反应谱确定其相应的平台频段。具体做法是：在对数坐标系中同时做出绝对加速度反应谱和拟速度反应谱，找出加速度反应谱平台段的起始周期T0和结束周期T1，然后在拟速度反应谱上选定平台段，其起始周期为T1(即加速度反应谱平台段的结束周期T1)，结束周期为T2，将加速度反应谱在T0至T1之间的谱值求平均得Sa，拟速度反应谱在T1至T2之间的谱值求平均得Sv，加速度反应谱和拟速度反应谱在平台段的放大系数采用2.5，按公式(1)、（2）、（3）求得EPA、EPV、Tg。在MIDAS程序中提供将地震波转换为绝对加速度反应谱和拟速度反应谱的功能(工具地震波数据生成器，生成后保存为SGS文件)，用户可利用保存的SGS文件(文本格式文件)根据上面所述方法计算Sv、Sa、Tg。通过Tg值可判断该地震波是否适合当地场地和地震设计分组，然后将抗震规范中表5.1.2-2中的EPA值与Sa相比求出调整系数，将其代入到地震波调整系数中。将地震波转换为绝对加速度反应谱和拟速度反应谱时注意周期范围要到6秒(建筑抗震规范规定)。建筑抗震设计规范5.1.2条中规定，采用时程分析方法时，应按照场地类别和设计地震分组选用不少于二组的实际强震记录和一组人工模拟

6.1 电磁场边界积分方程

第六章边界单元法有限元法属于偏微分方程法。对于求解有界电磁场域的场分布，尤其是有复杂边界和多种媒质、线性或非线性、静态或时变场的数值计算都是十分成功的，有的文献认为有限元法是应用最广，最重要的数值分析方法。当然，任何一种数值分析方法都不是万能的，有限元法的不足之处主要表现为： 1. 对于无界求解区域的处理比较困难； 2. 所求得的数值解是位函数值，再通过求导，一般比位值的精度低一个数量级，所以计算精度较低； 3. 对时变电磁场的求解，计算量太大。在以上这几点所反映的问题上，边界单元法解决得比较好，有明显优势。此外，边界单元法还具有能降低所研究问题的维数，离散剖分和数据准备简单等特点，它已成为计算场的重要方法，我们需要进行学习。 6.1 电磁场边界积分方程 6.1.1电磁场边界元方程的基本关系设三维线性泊松方程为所求场的控制方程，D 是具有边界面S 的求解区域。在S 上含有给定的第一和第二类边界条件的边界1S 和2S ，21S S S +=。对于这类恒定场，定解问题可表示为：式中：u 表示位函数，f 是场源密度函数（如ε ρ-）。若已求得近似解u ~ ，带入边值问题，用R 、1R 和2R 分别表示方程余量及边界余量：

f u R -?=~2 u u R S ~-=1 S q q R -=2 取权函数w ，按加权余量法，令误差分配的加权积分为： 021>=<->??<->

反应谱理论与人工模拟地震波技术简介

第３３卷第２６期?１０６?２００７年９月山西建筑ＳＨＡＮＸＩＡＲＣＨＩＴＥＣＴＵＲＥＶｄ３３Ｎｏ．２６Ｓｅｐ．２００７文章编号：１００９—６８２５｛２００７）２６—０１０６—０３反应谱理论与人工模拟地震波技术简介邱玉国王玉富摘要：介绍了反应谱理论的发展历程和国内外研究现状，分析了研究问题的思路，指出了利用反应谱理论来解决实际工程时遇到的问题，并简单介绍了国外对人工模拟地震波技术的应用和研究，为抗震理论提供了参考依据。关键词：反应谱理论，地震波，随机振动，非弹性地震波中图分类号：ＴＵ３５２文献标识码：Ａ１概述反应谱理论是建筑结构抗震设计的重要理论基础之一。从２０世纪５０年代开始，反应谱理论逐渐成为结构抗震设计的重要方法，经过５０多年的发展，目前这种方法已经为世界上大多数国家的设计规范所采用。但是，由于地震产生机理和作用效果的复杂性，采用反应谱理论进行分析和设计与工程实践还存在很多与实际不相符合之处。此外，对于反应地震重要特性的时间问题，反应谱法也无能为力。人工模拟地震波技术是近年来才发展起来的一项新的结构抗震设计的技术手段，目前主要用于计算机模拟和特别重要结构模型的振动台试验。它能够通过模拟地震波的特性来用于对结构进行时程分析，是～种新兴的、具有革命性意义的试验手段。图２数值模拟结果２．３计算结果分析通过数值模拟和试验得到瓦斯管承载力等数值如表２所示。表２数值模拟和试验结果Ｉ研究方法承载力仆但ａ最大应变／％最大剪应力／ＳＰａＩ数值模拟７．１４Ｏ．０８４２１６０室内试验６．６２０．０９６４３结语通过对丁集煤矿瓦斯管材质和整体抗外压的试验研究以及数值模拟分析，可以获得如下重要结论：１）通过对管材材质的试验研究表明：工作管材质采用Ｑ３４５，尺寸为柘３０ｒｆｌｌＴｌ×１４ｉｎｌｎ，能够满足强度和稳定性要求。２）瓦斯管整体抗外压试验结果表明：工作管抗外压承载力为６，６２ＭＰａ；通过大变形有限元数值计算，采用变形稳定性控制其承载力，结果为７．１４ＭＰａ，两者数值十分接近，说明用文中方法模拟大直径瓦斯管的承载力是可行的。参考文献：［１］李正来．瓦斯抽排钻孔定向技术的改进［Ｊ］．安徽科技，２００６（３）：４９—５０．［２］汪东生．瓦斯抽排技术治理本煤层采空区瓦斯涌出的实践［Ｊ］．煤矿安全，２００６（１）：１３—１５．［３］张敦伍，任胜杰．瓦斯抽排钻孔防偏斜实践［Ｊ］．矿业安全与环保，２００５（８）：６７—６８．［４］刘克功，范再良，赵新华．采空区瓦斯抽排法治理综放面瓦斯超限［Ｊ］．煤，１９９８（２）：４８—５０．ＳｔｕｄｙｉｎｇｏｎｒａｄｉａｌｓｔａｂｉｌｉｔｙｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎｏｆｌａｒｇｅｐｉｐｅｉｎｍｉｎｅＴＯＮＧＷｅｎ－ｌｉｎＡｂｓｔｒａｃｔ：ＴｈｅｅｘｐｅｒｉｍｅｎｔａｌａｎｄｖａｌｕｅｓｉｍｕｌａｔｉｏｎｍｅｔｈｏｄｓｈａｖｅｓｔｕｄｉｅｄｔｈｅＤｉｎｇｉｉｃｏａｌｍｉｎｅｌａｒｇｅｄｉａｍｅｔｅｒｇａｓｔｕｂｅｕｎｄｅｒｍｅｃｈａｎｉｃｓｃｈａｒａｃｔｅｒ—ｉｓｔｉｅ．Ｒｅｓｕｌｔｉｎｄｉｃａｔｅｄ：ｔｈｅｌａｒｇｅｄｉａｍｅｔｅｒｇａｓｔｕｂｅｉｓｐｒｅｓｅｎｔｅｄｓｔａｂｉｌｉｔｙｆａｉｌｕｒｅｍｏｄｅｌｉｎｅｎｃｉｒｃｌｅｓｐｒｅｓｓｅｓｓｈａｐｅ，ｉｔｓｓａｆｅｔｙｆａｃｔｏｒｒｅａｃｈｅｓ３．０，ｉｔｉｓｄｅｓｉｇｎｔｈｅｌａｒｇｅｄｉａｍｅｔｅｒｇａｓｔｕｂｅａｎｄｔｈｅｃｏｎｓｔｒｕｃｔｐｍｖｉｄｅｓｔｈｅｒｅｆｅｒｅｎｃｅ．Ｋｅｙｗｏｒｄｓ：ｌａｒｇｅｄｉａｍｅｔｅｒｇａｓｔｕｂｅ，ｅｘｐｅｒｉｍｅｎｔａｌｉｎｌａｂ，ｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎ，ｓｔａｂｉｌｉｔｙｆａｉｌｕｒｅｍｏｄｅｌ收稿日期：２００７．０４．０６作者简介：邱玉国（１９７３。），男，工程师，辽宁工程技术大学软件学院，辽宁阜新１２３０００王玉富（１９７０．），男，工程师，中铁十九局集团第三工程有限公司，辽宁辽阳１１１０００

6.3 边界积分方程的离散化方程

6.3 离散化边界积分方程的建立以二维边界离散化方程的建立为例，重点突出离散化方法的学习。 6.3.1建立Laplace 场的边界离散化方程电磁场边界元法的通用积分方程 (4) 其中： ?????? ?∈∈∈=域外光滑的边界上域内D D c i 0 211 设在Laplace 场中的二维边界上一点i 处，有方程：在二维场的边界线l 上进行离散，将l 划分为许多小段，每段以直线段或曲线段逼近，作为一个单元。设l 点共被分为0N 个单元，其中在第一类边界1l 段上划分了1N 个单元，在第二类边界2l 段上划分了2N 个单元： 210N N N += 作为单元待求量的插值计算方式，可分为几种： ① 恒值单元同一单元中的待求量u 和 n u ??都设为恒定值 (或称零次插值)，实际上是取单元中点的u 值(或 n u ??值)作为单元的u 值(或n u ??值)。这样，取单元中点为节点，所以求解变量数等于节点数。 ② 线性单元

它也是直线单元，其u 值在单元两端点之间按线性变化(即线性插值)。单元两端点为单元的节点。 ③ 曲线单元每单元上的节点数大于2，以多节点拟合的曲线逼近边界单元，以单元节点上的高阶插值函数作为待求位函数近似解。取最简单的单元——恒值单元为例，介绍边界元离散方法。按上面的方程对i 单元的“i ”节点离散化 ∑? ∑? ==??= ??+ o j o j N j l N j l i l n u F l n F u u 1 1 d d 2 1 ∑? ∑?=== ??+ 1 1 d d 2 1N j l N j l j j i j o j l F q l n F u u ，?= j l ij l F G d ，上式表示为：设i 点为i 单元的中点（021N i 、、、 =），有 ()∑∑==== 1 01 21N j N j j ij j ij N i q G u H ,,, 式中：于是上述0N 个方程写为矩阵形式 GQ HU = 由定解问题中的第一类边界1l ，对应有1N 个单元的位值u s 是已知的，2l 是第二类边界，对应有2N 个单元n u q s ??= 位是已知。所以上述矩阵方程中，有2N 个单元的u 值和 1N 个单元的q 值是未知的，即是说矩阵方程有021N N N =+个未知数。设单元排列顺序在1l 边界上为1，2，……，1N ，在2l 边是上为11+N 、21+N 、…、0N ，则上述矩阵方

0为什么能用地震波来探测地球内部的构造

为什么能用地震波来探测地球内部的构造？地震波是地震发生时，地下岩石受到强烈冲击所产生的弹性震动传播波。地震波是弹性波，它能穿过包括地核在内，在整个地球传播。地震波可分为纵波、横波、面波和界面波四种类型。纵波（P波），也称疏密波，通过物体时，物体质点的震动方向与地震波传播的方向一致，传播速度最快，周期短，振幅小，能通过固体、液体和气体传播。地震发生后，纵波最先到达地面，引起地面上下颠簸。横波（S波），通过物体时，物体的质点震动方向与地震波传播方向垂直，在地壳中传播速度比纵波慢，周期较长，振幅较大，只能通过固体介质传播，比纵波到达地面晚，横波能引起地面摇晃。纵波、横波合称体波，体波在地球体内部可以向任何方向传播。面波（L波），也称地面波，是纵波或横波到达地面后，从震中沿地面表层向四周传播的次生波。面波振幅较体波显著，波速比体波小，周期较体波长。利用面波的波散现象，可推算相应地区的地壳和上地幔的结构状况和性质。界面波是在两个弹性层之间的平界面附近传播的地震波。由于不同的地震波，具有不同的性质和传播特点，因此可以利用地震波来探测地

球的内部构造。目前世界上最深的钻井只有10公里多一点，能直接取样观察的最深矿井仅有3公里。目前人们还不能对地球整个内部进行直接观察研究，主要是利用地震波研究地球的内部结构。在地球内部地震波传播曲线图上，从地球大陆的地表面往下到33公里深处，横波速度每秒约4公里，纵波速度每秒约8公里。从33公里往下到2900公里深处，横波速度由每秒4公里多增快到每秒7公里以上，纵波速度由每秒8公里左右增快到每秒13公里以上。从2900公里往下到5000公里深处，横波完全消失，纵波传播速度突然下降到每秒8～10公里左右。从5000公里往下到地心，无横波传播，纵波速度又逐渐增快到每秒约11公里左右。从地震波在地球内传播的情况表明，在大陆33公里深处以下，横波和纵波的速度明显加快，证明是密度很大的可塑性固体层，因此地下33公里深处是地震波传播的一个不连续面，这个不连续面是莫霍洛维奇发现的，所以叫莫霍面。在2900公里深处往下，横波完全消失，纵波速度突然下降，证明到了液态层，这个地震波传播的不连续面，是古登堡最早研究的，所以叫古登堡面。5000公里以下纵波速度又加快，证明是固态层。根据地震波的传播情况，说地球内部构造是不同的物质圈层组成的。据此，人们以莫霍面和古登堡面为分界面，把地球的内部构造划分为地壳、地幔和地核三个圈层，并将地下2900～5000公里深处，推测

Matlab求解线性方程组非线性方程组

求解线性方程组 solve，linsolve 例： A=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1]; %矩阵的行之间用分号隔开，元素之间用逗号或空格 B=[3;1;1;0] X=zeros(4,1);%建立一个4元列向量 X=linsolve(A,B) diff（fun，var，n）：对表达式fun中的变量var求n阶导数。例如：F=sym（'u(x,y)*v(x,y)'）; %sym（）用来定义一个符号表达式 diff(F); %matlab区分大小写 pretty(ans) %pretty（）：用习惯书写方式显示变量；ans是答案表达式非线性方程求解 fsolve(fun,x0,options) 为待解方程或方程组的文件名；fun其中 x0位求解方程的初始向量或矩阵； option为设置命令参数建立文件fun.m： function y=fun(x) y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ... x(2) - 0.5*cos(x(1))+0.3*sin(x(2))]; >>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve')) 注： ...为续行符 m文件必须以function为文件头，调用符为@；文件名必须与定义的函数名相同；fsolve（）主要求解复杂非线性方程和方程组，求解过程是一个逼近过程。Matlab求解线性方程组 AX=B或XA=B 在MATLAB中，求解线性方程组时，主要采用前面章节介绍的除法运算符“/”和“\”。如： X=A\B表示求矩阵方程AX＝B的解；的解。XA=B表示矩阵方程B/A＝X．对方程组X＝A\B，要求A和B用相同的行数，X和B有相同的列数，它的行数等于矩阵A的列数，方程X＝B/A同理。如果矩阵A不是方阵，其维数是m×n，则有： m＝n 恰定方程，求解精确解； m>n 超定方程，寻求最小二乘解； m

固体力学中的边界积分方程及其边界元法综述

计算固体力学读书报告固体力学中的边界积分方程及其边界元法综述 Review of the Boundary Integral Equation and Boundary Element Method in Solid Mechanics 土木工程系 2014年03月17日

评语

目录摘要 (2) A BSTRACT (2) 一、引言 (3) 1)什么是边界元法[1] (3) 2)积分方程和边界元法的发展历史[2] (3) 二、边界元法[5] (4) 1)概述 (4) 2)基本解 (4) 3)拉普拉斯（Laplace）积分方程 (5) 4)拉普拉斯（Laplace）边界积分方程 (6) 5)拉普拉斯（Laplace）积分方程离散化与解法 (6) 6)泊松（Poisson）边界积分方程 (7) 三、结束语 (8) 参考文献 (9)

摘要本文综述了边界元法的历史、现状及发展，并对积分方程和边界元法的原理进行了简单推导。边界元法是在经典的积分方程的基础上，吸收了有限元法的离散技术而发展起来的计算方法，具有计算简单、适应性强、精度高的优点。它以边界积分方程为数学基础，同时采用了与有限元法相似的划分单元离散技术，通过将边界离散为边界元，将边界积分方程离散为代数方程组，再用数值方法求解代数方程组，从而得到原问题边界积分方程的解。用传统的有限单元法求解不可压缩材料会遇到严重困难，但是用边界元法求解这类材料不会有任何问题。近年来随着将快速多级算法引入边界元法，使边界元法的计算效率和解题规模都有了几个数量级的提高。关键词：边界元法积分方程边界离散快速多级算法 Abstract This paper reviews the history, current situation and development of the boundary element method and deduced the integral equation. The boundary element method is based on the integral equation and absorbed the discrete technology of finite element method. It has the advantages of simple calculation, strong adaptability and high accuracy. It is based on the boundary integral equation, though boundary discretization discrete boundary integral equations into algebraic equations, and then by the numerical method solving algebraic equations, thus obtain the original problem solution of boundary integral equations. The solution of nearly or exactly incompressible material problems presents serious difficulties and errors when using the conventional displacement-based finite element method, because the general stress-strain equations of elasticity contain terms that become infinite as Poisson’s ratio reaches 0.5, while the boundary element method accommodates such problems without any difficulty due to the nature of the integral equations used in the analysis. In recent years, the fast multi-pole boundary element method has received much attention because some large-scale engineering design and analysis problems were analyzed faster using boundary element method than with finite element method. This new trend suggests future prospects for boundary element method applications. Keywords：Boundary Element Method; Integral Equation; Boundary Discretization Method; Fast Multipole Algorithm

Midas地震波的选取方法

地震波的选取方法建筑抗震设计规范(GB 50011-2001)的5.1.2条文说明中规定，正确选择输入的地震加速度时程曲线，要满足地震动三要素的要求，即频谱特性、有效峰值和持续时间要符合规定。频谱特性可用地震影响系数曲线表征，依据所处的场地类别和设计地震分组确定。这句话的含义是选择的实际地震波所处场地的设计分组(震中距离、震级大小)和场地类别(场地条件)应与要分析的结构物所处场地的相同，简单的说两者的特征周期Tg值应接近或相同。特征周期Tg值的计算方法见下面公式(1)、(2)、(3)。加速度有效峰值按建筑抗震设计规范(GB 50011-2001)中的表5.1.2-2采用。地震波的加速度有效峰值的计算方法见下面公式(1)及下面说明。持续时间的概念不是指地震波数据中总的时间长度。持时T d的定义可分为两大类，一类是以地震动幅值的绝对值来定义的绝对持时，即指地震地面加速度值大于某值的时间总和，即绝对值｜a(t)｜＞k*g的时间总和，k常取为0.05；另一类为以相对值定义的相对持时，即最先与最后一个k*a max之间的时段长度，k一般取0.3～0.5。不论实际的强震记录还是人工模拟波形，一般持续时间取结构基本周期的5～10倍。说明: 有效峰值加速度 EPA＝Sa/2.5 （1）有效峰值速度 EPV＝Sv/2.5 （2）特征周期 Tg = 2π*EPV/EPA（3） 1978年美国ATC－3规范中将阻尼比为5％的加速度反应谱取周期为0.1-0.5秒之间的值平均为Sa，将阻尼比为5%的速度反应谱取周期为0.5-2秒之间的值平均为Sv(或取1s附近的平均速度反应谱)，上面公式中常数2.5为0.05组尼比加速度反应谱的平均放大系数。上述方法使用的是将频段固定的方法来求EPA和EPV，1990年的《中国地震烈度区划图》采用了不固定频段的方法分析各条反应谱确定其相应的平台频段。具体做法是：在对数坐标系中同时做出绝对加速度反应谱和拟速度反应谱，找出加速度反应谱平台段的起始周期T0和结束周期T1，然后在拟速度反应谱上选定平台段，其起始周期为T1(即加速度反应谱平台段的结束周期T1)，结束周期为T2，将加速度反应谱在T0至T1之间的谱值求平均得Sa，拟速度反应谱在T1至T2之间的谱值求平均得Sv（注：生成谱的时候一定要用对数谱），加速度反应谱和拟速度反应谱在平台段的放大系数采用2.5，按公式(1)、（2）、（3）求得EPA、EPV、Tg。在MIDAS程序中提供将地震波转换为绝对加速度反应谱和拟速度反应谱的功能(工具>地震波数据生成器，生成后保存为SGS文件)，用户可利用保存的SGS文件(文本格式文件)根据上面所述方法计算Sv、Sa、Tg＝Sv/Sa。通过Tg值可判断该地震波是否适合当地场地和地震设计分组，然后将抗震规范中表5.1.2-2中的EPA值与Sa相比求出调整系数（即放大系数），将其代入到地震波调整系数中。将地震波转换为绝对加速度反应谱和拟速度反应谱时注意周期范围要到6秒(建筑抗震规范规定)。建筑抗震设计规范5.1.2条中规定，采用时程分析方法时，应按照场地类别和设计地震分组选用不少于二组的实际强震记录和一组人工模拟的加速度时程曲线，其平均地震影响系数曲线应与振型分解反应谱法所采用的地震影响系数曲线在统计意义上相符。所谓“在统计意义上相符”指的是，其平均影响系数曲线与振型分解反应谱法所用的地震影响系数曲线相比，在各周期点上相差不大于20%。在MIDAS程序中，可选取两组实际强震记录生成两个SGS文件(调整Sa后的)，然后将一组人

人工地震波生成程序简介

姓名：郭勇学号：022******* 人工地震波生成程序简介一、程序设计内容及方法 1、程序内容本程序根据特征周期、水平地震波影响系数最大值和地震波幅值等初始条件生成人工地震波，为结构动力分析的时程分析法提供地震波来源。 2、程序设计方法 (1) 理论依据本程序采用三角级数法生成人工地震波。对于给定的功率谱密度函数，按照下面的公式可以方便的生成以为功率谱密度函数、均值为零的高斯平稳过程。 (1) 式中： (2) 为内均匀分布的随机相角；，分别为正域内的上、下限值，即认为的有效功率在范围内，而范围外的值可视为零。为了反映地面运动的非平稳性，采用包络函数乘以平稳过程， (3) (3)式即为人工地震波模型。可根据下式确定： (4) 式中：为衰减系数，通常取值范围为0.1~1.0，本程序取0.15；，和根据不同实际情况取值，为地震波持时，本程序取，分别为4s，15s，和均为40s。本程序采用《建筑抗震设计规范》(GB50011-2001)中的反应谱作为目标谱，通过Kaul 提出的平稳过程反应谱与功率谱的近似关系 (5) 式中：为规范反应谱；为阻尼比；为地震动持时；为反应不超过反应谱值的概率，本程序取0.85。通过(3)式和(5)式即可生成人工地震波。 (2) 程序实现方法首先建立基于对话框的应用程序框架，添加的主要控件为3个编辑框和4个按钮。3个编辑框分别作为程序中的特征周期(对应成员变量为m_dTg)、水平地震影响系数最大值(对应成员变量为m_dAmax)和地震波幅值(对应成员变量为m_pd)3个数据的交互输入处；4个按钮分别为"生成地震波"、"输出地震波"、"输入地震波"和"退出"。添加的成员函数有：Wavegener()(生成地震波)、Wavedrawing()(绘制地震波加速度时程曲线)、OnSTART()(对应"生成地震波"按钮，实现生成地震波的功能)、OnOutput()(对应"输出地震波"按钮，实现输出数字化的地震波记录的功能)和OnInput(对应"输入地震波"按钮，实现输入数字化的地震波记录并绘制其加速度时程曲线的功能)。几点说明： a 生成随机相角的程序如下： srand((unsigned)time( NULL ));

差分方程的解法分析及MATLAB实现(程序)

差分方程的解法分析及MATLAB 实现（程序）摘自：张登奇,彭仕玉.差分方程的解法分析及其MATLAB 实现[J]. 湖南理工学院学报.2014(03) 引言线性常系数差分方程是描述线性时不变离散时间系统的数学模型,求解差分方程是分析离散时间系统的重要内容.在《信号与系统》课程中介绍的求解方法主要有迭代法、时域经典法、双零法和变换域法[1]. 1 迭代法例1 已知离散系统的差分方程为)1(3 1)()2(81)1(43)(-+=-+--n x n x n y n y n y ,激励信号为)()4 3()(n u n x n =,初始状态为21)2(4)1(=-=-y y ，.求系统响应. 根据激励信号和初始状态,手工依次迭代可算出24 59)1(,25)0(==y y . 利用MATLAB 中的filter 函数实现迭代过程的m 程序如下: clc;clear;format compact; a=[1,-3/4,1/8],b=[1,1/3,0], %输入差分方程系数向量,不足补0对齐 n=0:10;xn=(3/4).^n, %输入激励信号 zx=[0,0],zy=[4,12], %输入初始状态 zi=filtic(b,a,zy,zx),%计算等效初始条件 [yn,zf]=filter(b,a,xn,zi),%迭代计算输出和后段等效初始条件 2 时域经典法用时域经典法求解差分方程:先求齐次解；再将激励信号代入方程右端化简得自由项,根据自由项形式求特解；然后根据边界条件求完全解[3].用时域经典法求解例1的基本步骤如下. （1）求齐次解.特征方程为081432=+-αα,可算出4 1 , 2121==αα.高阶特征根可用MATLAB 的roots 函数计算.齐次解为. 0 , )4 1()21()(21≥+=n C C n y n n h （2）求方程的特解.将)()4 3()(n u n x n =代入差分方程右端得自由项为 ?????≥?==-?+-1,)4 3(9130 ,1)1()43(31)()43(1n n n u n u n n n 当1≥n 时,特解可设为n p D n y )4 3()(=,代入差分方程求得213=D . （3）利用边界条件求完全解.当n =0时迭代求出25)0(=y ,当n ≥1时,完全解的形式为 ,)4 3(213 )41()21()(21n n n C C n y ?++=选择求完全解系数的边界条件可参考文[4]选)1(),0(-y y .根据边界条件求得35,31721=-=C C .注意完全解的表达式只适于特解成立的n 取值范围,其他点要用 )(n δ及其延迟表示,如果其值符合表达式则可合并处理.差分方程的完全解为

反应谱生成人工地震波

反应谱生成人工地震波一、软件SIMQKE_GR使用说明 1.先安装程序 2.使用方法双击，打开程序，可以得到如图1界面。图1 程序开始界面如图1所示，由于程序本身提供的反应谱是适用于欧洲规范的，不适合于我国的规范反应谱，因此不能通过调整参数来获得符合我国规范的反应谱。可以采用导入的方法来输入反应谱。 3.点击菜单栏“file”—“Import spectra data”，出现打开对话框，如图2所示，要求打开一个已经存在的反应谱文件（如 1.srf）。

图2 导入反应谱文件对话框 4.文件格式如下所示（红字部分不能修改，注意反应谱单位为g），下面部分可以替换。 response spectrum time(s) acc(g) 0 0.1215 0.01 0.13635 0.02 0.1512 0.03 0.16605 0.04 0.1809 0.05 0.19575 0.06 0.2106 0.07 0.22545 0.08 0.2403 0.09 0.25515 0.1 0.27 0.15 0.27 0.2 0.27 0.25 0.27 0.3 0.27 0.35 0.27 0.4 0.27 0.45 0.27

0.5 0.243 0.6 0.2025 0.7 0.173571429 0.8 0.151875 0.9 0.135 1 0.1215 1.1 0.110454545 1.2 0.10125 1.3 0.093461538 1.4 0.086785714 1.5 0.081 1.6 0.0759375 1.7 0.071470588 1.8 0.0675 1.9 0.063947368 2 0.06075 2.1 0.057857143 2.2 0.055227273 2.3 0.052826087 2.4 0.050625 2.5 0.0486 2.6 0.046730769 2.7 0.045 2.8 0.043392857 2.9 0.041896552 3 0.0405 3.1 0.039193548 3.2 0.03796875 3.3 0.036818182 3.4 0.035735294 3.5 0.034714286 3.6 0.03375 3.7 0.032837838 3.8 0.031973684 3.9 0.031153846 4 0.030375 4.1 0.029634146 4.2 0.028928571 4.3 0.028255814 4.4 0.027613636 4.5 0.027 4.6 0.026413043 4.7 0.025851064 4.8 0.0253125

MATLAB应用求解非线性方程

第7章求解非线性方程 7.1 多项式运算在MATLAB 中的实现一、多项式的表达 n 次多项式表达为：n a +??++=x a x a x a p(x )1-n 1-n 1n 0，是n+1项之和在MATLAB 中，n 次多项式可以用n 次多项式系数构成的长度为n+1的行向量表示 [a0, a1,……an-1,an] 二、多项式的加减运算设有两个多项式 n a +??++=x a x a x a p1(x )1-n 1-n 1n 0和 m b +??++=x b x b x b p2(x )1-m 1-m 1m 0。它们的加减运算实际上就是它们的对应系数的加减运算。当它们的次数相同时，可以直接对多项式的系数向量进行加减运算。当它们的次数不同时，应该把次数低的多项式无高次项部分用0系数表示。例2 计算()()1635223-+++-x x x x a=[1, -2, 5, 3]; b=[0, 0, 6, -1]; c=a+b 例3 设()6572532345++-+-=x x x x x x f ，()3532-+=x x x g ，求f(x)+g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; g1=[0, 0, 0, g];%为了和f 的次数找齐 f+g1, f-g1 三、多项式的乘法运算 conv(p1,p2) 例4 在上例中，求f(x)*g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; conv(f, g) 四、多项式的除法运算 [Q, r]=deconv(p1, p2) 表示p1除以p2，给出商式Q(x)，余式r(x)。Q,和r 仍为多项式系数向量例4 在上例中，求f(x)/g(x) f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; [Q, r]=deconv(f, g) 五、多项式的导函数 p=polyder(P)：求多项式P 的导函数 p=polyder(P,Q)：求P·Q 的导函数

八层框架的地震响应计算和人工波生成的matlab实现及所需曲线图的自动存储

一、作业概况结构基本参数：层间剪切型结构，采用Rayleigh 阻尼，第一、第二阶阻尼比分别取3%、5%。图1 结构基本形状表1 各层集中质量 ( 105kg) 层号 1 2 3 4 5 6 7 8 质量 3.40 3.40 3.20 3.20 2.80 2.80 2.70 2.60 表2 各层层间刚度 (×108N/m) 层号 1 2 3 4 5 6 7 8 层间刚度 2.00 2.00 1.80 1.80 1.80 1.80 1.60 1.60 m m m m m m m m &&g x

二、频率及振型计算根据层间模型的假定，可以建立结构的质量矩阵以及刚度矩阵如下。 1234567800000000000000000000000000000000000000000000000000000000 3.40 0000000 3.400000000 3.200000000 3.20000 =0000 2.800000000 2.800000000 2.700000000 2.6m m m m m m m m ?? ? ? ? ? ?= ? ? ? ? ? ?? ? ?? ?M 510kg ????????????? 1112131415161718212223242526272831323334353637384142 4344454647485152535455565758616263646566676871727374757677788182838485868788k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k ? =?K 8420000002 3.8 1.8000000 1.8 3.6 1.8000000 1.8 3.6 1.8000 =10/000 1.8 3.6 1.8000000 1.8 3.4 1.6000000 1.6 3.2 1.6000000 1.6 1.6N m ????? ? ? ? ? ? ? ??-?? ?-- ? ?-- ?-- ?? ?-- ?-- ? ?-- ? ?-??

地震波方程人工边界的两种处理方法

推荐-Broyden方法求解非线性方程组的Matlab实现 精品

地震波使用说明

地震波的选取方法 (MIDAS内部技术资料)

6.1 电磁场边界积分方程

反应谱理论与人工模拟地震波技术简介

6.3 边界积分方程的离散化方程

0为什么能用地震波来探测地球内部的构造

Matlab求解线性方程组非线性方程组

固体力学中的边界积分方程及其边界元法综述

Midas地震波的选取方法

人工地震波生成程序简介

差分方程的解法分析及MATLAB实现(程序)

反应谱生成人工地震波

MATLAB应用 求解非线性方程

八层框架的地震响应计算和人工波生成的matlab实现及所需曲线图的自动存储

推荐-Broyden方法求解非线性方程组的Matlab实现精品

MATLAB应用求解非线性方程